From 0ac85165d7550f5206593b0ea6440979050fa736 Mon Sep 17 00:00:00 2001 From: Gereon Kremer Date: Wed, 2 Mar 2022 18:21:32 +0100 Subject: [PATCH] Move libpoly <-> CoCoA conversion to new utility (#8199) MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit This PR factors the conversion between libpoly and CoCoA out of the LazardEvaluation class into a new utility. The goal is to make it reusable for other applications of CoCoA within cvc5, for example to replace the (rather naive) equality substitution utility by a new simplification technique based on CoCoAs Gröbner bases. This PR merely moves code around and should not actually change anything. --- src/CMakeLists.txt | 2 + .../arith/nl/coverings/cocoa_converter.cpp | 126 ++++++++++++ .../arith/nl/coverings/cocoa_converter.h | 141 ++++++++++++++ .../arith/nl/coverings/lazard_evaluation.cpp | 179 +++--------------- 4 files changed, 291 insertions(+), 157 deletions(-) create mode 100644 src/theory/arith/nl/coverings/cocoa_converter.cpp create mode 100644 src/theory/arith/nl/coverings/cocoa_converter.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e29b234ba..94c077ced 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -416,6 +416,8 @@ libcvc5_add_sources( theory/arith/nl/coverings/cdcac.h theory/arith/nl/coverings/cdcac_utils.cpp theory/arith/nl/coverings/cdcac_utils.h + theory/arith/nl/coverings/cocoa_converter.cpp + theory/arith/nl/coverings/cocoa_converter.h theory/arith/nl/coverings/constraints.cpp theory/arith/nl/coverings/constraints.h theory/arith/nl/coverings/lazard_evaluation.cpp diff --git a/src/theory/arith/nl/coverings/cocoa_converter.cpp b/src/theory/arith/nl/coverings/cocoa_converter.cpp new file mode 100644 index 000000000..4165c1815 --- /dev/null +++ b/src/theory/arith/nl/coverings/cocoa_converter.cpp @@ -0,0 +1,126 @@ +/****************************************************************************** + * Top contributors (to current version): + * Gereon Kremer + * + * This file is part of the cvc5 project. + * + * Copyright (c) 2009-2021 by the authors listed in the file AUTHORS + * in the top-level source directory and their institutional affiliations. + * All rights reserved. See the file COPYING in the top-level source + * directory for licensing information. + * **************************************************************************** + * + * Generic utility to convert between libpoly polynomials and CoCoALib + * polynomials. + */ + +#include "theory/arith/nl/coverings/cocoa_converter.h" + +#ifdef CVC5_POLY_IMP +#ifdef CVC5_USE_COCOA + +namespace cvc5 { +namespace theory { +namespace arith { +namespace nl { +namespace coverings { +CoCoA::RingElem CoCoAConverter::operator()(const poly::UPolynomial& p, + const poly::Variable& var, + const CoCoA::ring& ring) const +{ + std::vector coeffs = poly::coefficients(p); + CoCoA::RingElem res(ring); + CoCoA::RingElem v = CoCoA::indet(ring, 0); + CoCoA::RingElem mult(ring, 1); + for (const auto& c : coeffs) + { + if (!poly::is_zero(c)) + { + res += (*this)(c)*mult; + } + mult *= v; + } + Trace("nl-cov::cocoa") << "Converted " << p << " to " << res << std::endl; + return res; +} + +CoCoA::RingElem CoCoAConverter::operator()(const poly::Polynomial& q, + const CoCoA::ring& ring) const +{ + CoCoAPolyConstructor cmd{*this, ring}; + // Do the actual conversion + cmd.d_result = CoCoA::RingElem(ring); + lp_polynomial_traverse_f f = + [](const lp_polynomial_context_t* ctx, lp_monomial_t* m, void* data) { + CoCoAPolyConstructor* d = static_cast(data); + CoCoA::BigInt coeff = (d->d_state)(*poly::detail::cast_from(&m->a)); + CoCoA::RingElem re(d->d_ring, coeff); + for (size_t i = 0; i < m->n; ++i) + { + // variable exponent pair + CoCoA::RingElem var = d->d_state.d_varPC.at(m->p[i].x); + re *= CoCoA::power(var, m->p[i].d); + } + d->d_result += re; + }; + lp_polynomial_traverse(q.get_internal(), f, &cmd); + Trace("nl-cov::cocoa") << "Converted " << q << " to " << cmd.d_result + << std::endl; + return cmd.d_result; +} + +poly::Polynomial CoCoAConverter::convertImpl(const CoCoA::RingElem& p, + poly::Integer& denominator) const +{ + denominator = poly::Integer(1); + poly::Polynomial res; + for (CoCoA::SparsePolyIter i = CoCoA::BeginIter(p); !CoCoA::IsEnded(i); ++i) + { + poly::Polynomial coeff; + poly::Integer denom(1); + CoCoA::BigRat numcoeff; + if (CoCoA::IsRational(numcoeff, CoCoA::coeff(i))) + { + poly::Rational rat(mpq_class(CoCoA::mpqref(numcoeff))); + denom = poly::denominator(rat); + coeff = poly::numerator(rat); + } + else + { + coeff = convertImpl(CoCoA::CanonicalRepr(CoCoA::coeff(i)), denom); + } + if (!CoCoA::IsOne(CoCoA::PP(i))) + { + std::vector exponents; + CoCoA::exponents(exponents, CoCoA::PP(i)); + for (size_t vid = 0; vid < exponents.size(); ++vid) + { + if (exponents[vid] == 0) continue; + const auto& ring = CoCoA::owner(p); + poly::Variable v = d_varCP.at(std::make_pair(CoCoA::RingID(ring), vid)); + coeff *= poly::Polynomial(poly::Integer(1), v, exponents[vid]); + } + } + if (denom != denominator) + { + poly::Integer g = gcd(denom, denominator); + res = res * (denom / g) + coeff * (denominator / g); + denominator *= (denom / g); + } + else + { + res += coeff; + } + } + Trace("nl-cov::cocoa") << "Converted " << p << " to " << res << std::endl; + return res; +} + +} // namespace coverings +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace cvc5 + +#endif +#endif diff --git a/src/theory/arith/nl/coverings/cocoa_converter.h b/src/theory/arith/nl/coverings/cocoa_converter.h new file mode 100644 index 000000000..432464fc9 --- /dev/null +++ b/src/theory/arith/nl/coverings/cocoa_converter.h @@ -0,0 +1,141 @@ +/****************************************************************************** + * Top contributors (to current version): + * Gereon Kremer + * + * This file is part of the cvc5 project. + * + * Copyright (c) 2009-2021 by the authors listed in the file AUTHORS + * in the top-level source directory and their institutional affiliations. + * All rights reserved. See the file COPYING in the top-level source + * directory for licensing information. + * **************************************************************************** + * + * Generic utility to convert between libpoly polynomials and CoCoALib + * polynomials. + */ + +#include "cvc5_private.h" + +#ifndef CVC5__THEORY__ARITH__NL__COVERINGS__COCOA_CONVERTER_H +#define CVC5__THEORY__ARITH__NL__COVERINGS__COCOA_CONVERTER_H + +#ifdef CVC5_POLY_IMP +#ifdef CVC5_USE_COCOA + +#include +#include + +#include + +#include "base/output.h" + +namespace cvc5 { +namespace theory { +namespace arith { +namespace nl { +namespace coverings { + +/** + * This class implements a generic converter utility between libpoly polynomials + * and CoCoA polynomials. + */ +class CoCoAConverter +{ + public: + /** Add mapping from libpoly variable to a CoCoA variable */ + void addVar(const poly::Variable& pv, const CoCoA::RingElem& cv) + { + d_varPC.emplace(pv, cv); + } + /** + * Add mapping from a CoCoA variable, identified by its ring id and the + * indet identifier, to a libpoly polynomial. + */ + void addVar(const std::pair& cv, const poly::Variable& pv) + { + d_varCP.emplace(cv, pv); + } + + /** + * Converts a libpoly integer to a CoCoA::BigInt. + */ + CoCoA::BigInt operator()(const poly::Integer& i) const + { + return CoCoA::BigIntFromMPZ(poly::detail::cast_to_gmp(&i)->get_mpz_t()); + } + /** + * Converts a libpoly dyadic rational to a CoCoA::BigRat. + */ + CoCoA::BigRat operator()(const poly::DyadicRational& dr) const + { + return CoCoA::BigRat((*this)(poly::numerator(dr)), + (*this)(poly::denominator(dr))); + } + /** + * Converts a libpoly rational to a CoCoA::BigRat. + */ + CoCoA::BigRat operator()(const poly::Rational& r) const + { + return CoCoA::BigRatFromMPQ(poly::detail::cast_to_gmp(&r)->get_mpq_t()); + } + /** + * Converts a univariate libpoly polynomial p in variable var to a CoCoA + * polynomial in the given CoCoA ring. + */ + CoCoA::RingElem operator()(const poly::UPolynomial& p, + const poly::Variable& var, + const CoCoA::ring& ring) const; + + /** + * Converts a libpoly polynomial p in variable var to a CoCoA polynomial in + * the given CoCoA ring. + */ + CoCoA::RingElem operator()(const poly::Polynomial& q, + const CoCoA::ring& ring) const; + + /** + * Converts a CoCoA polynomial to a libpoly polynomial. + */ + poly::Polynomial operator()(const CoCoA::RingElem& p) const + { + poly::Integer denom; + return convertImpl(p, denom); + } + + private: + + /** Helper class for the conversion of a libpoly polynomial to CoCoA. */ + struct CoCoAPolyConstructor + { + const CoCoAConverter& d_state; + const CoCoA::ring& d_ring; + CoCoA::RingElem d_result; + }; + /** Helper method for the conversion of a CoCoA polynomial to libpoly */ + poly::Polynomial convertImpl(const CoCoA::RingElem& p, + poly::Integer& denominator) const; + + /** Some global state that CoCoA needs to be around whenever it is used */ + CoCoA::GlobalManager d_gm; + + /** + * Maps libpoly variables to indets in CoCoA. + */ + std::map d_varPC; + + /** + * Maps CoCoA variables (identified by ring id and indet id) to libpoly + * variables. + */ + std::map, poly::Variable> d_varCP; +}; + +} // namespace coverings +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace cvc5 + +#endif +#endif +#endif diff --git a/src/theory/arith/nl/coverings/lazard_evaluation.cpp b/src/theory/arith/nl/coverings/lazard_evaluation.cpp index 2332a8fa3..53d0898dd 100644 --- a/src/theory/arith/nl/coverings/lazard_evaluation.cpp +++ b/src/theory/arith/nl/coverings/lazard_evaluation.cpp @@ -13,6 +13,8 @@ #include +#include "theory/arith/nl/coverings/cocoa_converter.h" + namespace cvc5::theory::arith::nl::coverings { struct LazardEvaluationStats @@ -105,8 +107,6 @@ std::ostream& operator<<(std::ostream& os, const LazardEvaluationState& state); */ struct LazardEvaluationState { - CoCoA::GlobalManager d_gm; - /** * Statistics about the lazard evaluation. * Although this class is short-lived, there is no need to make the statistics @@ -117,19 +117,12 @@ struct LazardEvaluationState mutable LazardEvaluationStats d_stats; /** - * Maps libpoly variables to indets in J0. Used when constructing the input - * polynomial q in the first polynomial ring J0. - */ - std::map d_varQ; - /** - * Maps CoCoA indets back to to libpoly variables. - * Use when converting CoCoA RingElems to libpoly polynomials, either when - * checking whether a factor vanishes or when returning the univariate - * elements of the final Gröbner basis. The CoCoA indets are identified by the - * pair of the ring id and the indet identifier. Hence, we can put all of them - * in one map, no matter which ring they belong to. + * Converter between libpoly and CoCoA. + * + * d_converter.d_varPC * Maps libpoly variables to indets in J0. Used when + * constructing the input polynomial q in the first polynomial ring J0. */ - std::map, poly::Variable> d_varCoCoA; + CoCoAConverter d_converter; /** * The minimal polynomials p_i used for constructing d_K. @@ -193,28 +186,6 @@ struct LazardEvaluationState */ std::vector> d_direct; - /** - * Converts a libpoly integer to a CoCoA::BigInt. - */ - CoCoA::BigInt convert(const poly::Integer& i) const - { - return CoCoA::BigIntFromMPZ(poly::detail::cast_to_gmp(&i)->get_mpz_t()); - } - /** - * Converts a libpoly dyadic rational to a CoCoA::BigRat. - */ - CoCoA::BigRat convert(const poly::DyadicRational& dr) const - { - return CoCoA::BigRat(convert(poly::numerator(dr)), - convert(poly::denominator(dr))); - } - /** - * Converts a libpoly rational to a CoCoA::BigRat. - */ - CoCoA::BigRat convert(const poly::Rational& r) const - { - return CoCoA::BigRatFromMPQ(poly::detail::cast_to_gmp(&r)->get_mpq_t()); - } /** * Converts a univariate libpoly polynomial p in variable var to CoCoA. It * assumes that p is a minimal polynomial p_i over variable x_i for the @@ -223,19 +194,7 @@ struct LazardEvaluationState CoCoA::RingElem convertMiPo(const poly::UPolynomial& p, const poly::Variable& var) const { - std::vector coeffs = poly::coefficients(p); - CoCoA::RingElem res(d_R.back()); - CoCoA::RingElem v = CoCoA::indet(d_R.back(), 0); - CoCoA::RingElem mult(d_R.back(), 1); - for (const auto& c : coeffs) - { - if (!poly::is_zero(c)) - { - res += convert(c) * mult; - } - mult *= v; - } - return res; + return d_converter(p, var, d_R.back()); } /** @@ -246,7 +205,7 @@ struct LazardEvaluationState bool evaluatesToZero(const CoCoA::RingElem& cp) const { Assert(CoCoA::owner(cp) == d_R.back()); - poly::Polynomial pp = convert(cp); + poly::Polynomial pp = d_converter(cp); return poly::evaluate_constraint(pp, d_assignment, poly::SignCondition::EQ); } @@ -319,7 +278,7 @@ struct LazardEvaluationState * - add variable x_i to d_variables * - extract the variable name * - construct R_i = K_i[x_i] - * - add new variable to d_varCoCoA + * - add new variable mapping to d_converter */ void addR(const poly::Variable& var) { @@ -336,7 +295,7 @@ struct LazardEvaluationState } Trace("nl-cov::lazard") << "R" << d_R.size() - 1 << " = " << d_R.back() << std::endl; - d_varCoCoA.emplace(std::make_pair(CoCoA::RingID(d_R.back()), 0), var); + d_converter.addVar(std::make_pair(CoCoA::RingID(d_R.back()), 0), var); } /** @@ -378,8 +337,8 @@ struct LazardEvaluationState * - construct all J_i * - construct all p homomorphisms (R_i --> J_i) * - construct all q homomorphisms (J_i --> J_{i+1}) - * - fill the mapping d_varQ (libpoly -> J_0) - * - fill the mapping d_varCoCoA (J_n -> libpoly) + * - add the variable mapping d_converter (libpoly -> J_0) + * - add the variable mapping d_converter (J_n -> libpoly) * - fill d_GBBaseIdeal with p_i mapped to J_0 */ void addFreeVariable(const poly::Variable& var) @@ -455,11 +414,8 @@ struct LazardEvaluationState } for (size_t i = 0; i < d_variables.size(); ++i) { - d_varQ.emplace(d_variables[i], CoCoA::indet(d_J[0], i)); - } - for (size_t i = 0; i < d_variables.size(); ++i) - { - d_varCoCoA.emplace(std::make_pair(CoCoA::RingID(d_J[0]), i), + d_converter.addVar(d_variables[i], CoCoA::indet(d_J[0], i)); + d_converter.addVar(std::make_pair(CoCoA::RingID(d_J[0]), i), d_variables[i]); } @@ -476,104 +432,12 @@ struct LazardEvaluationState << *this << std::endl; } - /** - * Helper class for conversion from libpoly to CoCoA polynomials. - * The lambda can not capture anything, as it needs to be of type - * lp_polynomial_traverse_f. - */ - struct CoCoAPolyConstructor - { - const LazardEvaluationState& d_state; - CoCoA::RingElem d_result; - }; - /** * Convert the polynomial q to CoCoA into J_0. */ CoCoA::RingElem convertQ(const poly::Polynomial& q) const { - CoCoAPolyConstructor cmd{*this}; - // Do the actual conversion - cmd.d_result = CoCoA::RingElem(d_J[0]); - lp_polynomial_traverse_f f = [](const lp_polynomial_context_t* ctx, - lp_monomial_t* m, - void* data) { - CoCoAPolyConstructor* d = static_cast(data); - CoCoA::BigInt coeff = d->d_state.convert(*poly::detail::cast_from(&m->a)); - CoCoA::RingElem re(d->d_state.d_J[0], coeff); - for (size_t i = 0; i < m->n; ++i) - { - // variable exponent pair - CoCoA::RingElem var = d->d_state.d_varQ.at(m->p[i].x); - re *= CoCoA::power(var, m->p[i].d); - } - d->d_result += re; - }; - lp_polynomial_traverse(q.get_internal(), f, &cmd); - return cmd.d_result; - } - /** - * Actual (recursive) implementation of converting a CoCoA polynomial to a - * libpoly polynomial. As libpoly polynomials only have integer coefficients, - * we need to maintain an integer denominator to normalize all terms to the - * same denominator. - */ - poly::Polynomial convertImpl(const CoCoA::RingElem& p, - poly::Integer& denominator) const - { - Trace("nl-cov::lazard") << "Converting " << p << std::endl; - denominator = poly::Integer(1); - poly::Polynomial res; - for (CoCoA::SparsePolyIter i = CoCoA::BeginIter(p); !CoCoA::IsEnded(i); ++i) - { - poly::Polynomial coeff; - poly::Integer denom(1); - CoCoA::BigRat numcoeff; - if (CoCoA::IsRational(numcoeff, CoCoA::coeff(i))) - { - poly::Rational rat(mpq_class(CoCoA::mpqref(numcoeff))); - denom = poly::denominator(rat); - coeff = poly::numerator(rat); - } - else - { - coeff = convertImpl(CoCoA::CanonicalRepr(CoCoA::coeff(i)), denom); - } - if (!CoCoA::IsOne(CoCoA::PP(i))) - { - std::vector exponents; - CoCoA::exponents(exponents, CoCoA::PP(i)); - for (size_t vid = 0; vid < exponents.size(); ++vid) - { - if (exponents[vid] == 0) continue; - const auto& ring = CoCoA::owner(p); - poly::Variable v = - d_varCoCoA.at(std::make_pair(CoCoA::RingID(ring), vid)); - coeff *= poly::Polynomial(poly::Integer(1), v, exponents[vid]); - } - } - if (denom != denominator) - { - poly::Integer g = gcd(denom, denominator); - res = res * (denom / g) + coeff * (denominator / g); - denominator *= (denom / g); - } - else - { - res += coeff; - } - } - Trace("nl-cov::lazard") << "-> " << res << std::endl; - return res; - } - /** - * Actually convert a CoCoA RingElem to a libpoly polynomial. - * Requires d_varCoCoA to be filled appropriately. - */ - poly::Polynomial convert(const CoCoA::RingElem& p) const - { - poly::Integer denom; - return convertImpl(p, denom); + return d_converter(q, d_J[0]); } /** @@ -649,7 +513,7 @@ struct LazardEvaluationState for (const auto& belem : basis) { Trace("nl-cov::lazard") << "-> retrieved " << belem << std::endl; - auto pres = convert(belem); + auto pres = d_converter(belem); Trace("nl-cov::lazard") << "-> converted " << pres << std::endl; // These checks are orthogonal! if (poly::is_univariate(pres) @@ -716,7 +580,7 @@ void LazardEvaluation::add(const poly::Variable& var, const poly::Value& val) const poly::DyadicInterval& di = poly::get_isolating_interval(ran); if (poly::is_point(di)) { - rational = d_state->convert(poly::get_point(di)); + rational = d_state->d_converter(poly::get_point(di)); } else { @@ -730,15 +594,16 @@ void LazardEvaluation::add(const poly::Variable& var, const poly::Value& val) || poly::is_rational(val)); if (poly::is_dyadic_rational(val)) { - rational = d_state->convert(poly::as_dyadic_rational(val)); + rational = d_state->d_converter(poly::as_dyadic_rational(val)); } else if (poly::is_integer(val)) { - rational = CoCoA::BigRat(d_state->convert(poly::as_integer(val)), 1); + rational = + CoCoA::BigRat(d_state->d_converter(poly::as_integer(val)), 1); } else if (poly::is_rational(val)) { - rational = d_state->convert(poly::as_rational(val)); + rational = d_state->d_converter(poly::as_rational(val)); } } -- 2.30.2