From 3f0fd456553223a59001a6bcdc714c3b56192787 Mon Sep 17 00:00:00 2001 From: Gereon Kremer Date: Thu, 24 Jun 2021 13:27:30 +0200 Subject: [PATCH] Add CoCoA implementation (#6733) This PR adds the actual implementation for the Lazard evaluation based on CoCoALib. It is only used if CoCoALib is available and falls back to a default libpoly-based implementation otherwise. --- src/theory/arith/nl/cad/lazard_evaluation.cpp | 928 ++++++++++++++++++ 1 file changed, 928 insertions(+) diff --git a/src/theory/arith/nl/cad/lazard_evaluation.cpp b/src/theory/arith/nl/cad/lazard_evaluation.cpp index 2fee21cbf..aec0d46e3 100644 --- a/src/theory/arith/nl/cad/lazard_evaluation.cpp +++ b/src/theory/arith/nl/cad/lazard_evaluation.cpp @@ -4,6 +4,933 @@ #include "base/check.h" #include "base/output.h" +#include "smt/smt_statistics_registry.h" +#include "util/statistics_stats.h" + +#ifdef CVC5_USE_COCOA + +#include + +#include + +namespace cvc5::theory::arith::nl::cad { + +struct LazardEvaluationStats +{ + IntStat d_directAssignments = + smtStatisticsRegistry().registerInt("theory::arith::cad::lazard-direct"); + IntStat d_ranAssignments = + smtStatisticsRegistry().registerInt("theory::arith::cad::lazard-rans"); + IntStat d_evaluations = + smtStatisticsRegistry().registerInt("theory::arith::cad::lazard-evals"); + IntStat d_reductions = + smtStatisticsRegistry().registerInt("theory::arith::cad::lazard-reduce"); +}; + +struct LazardEvaluationState; +std::ostream& operator<<(std::ostream& os, const LazardEvaluationState& state); + +/** + * This class holds and implements all the technicalities required to map + * polynomials from libpoly into CoCoALib, perform these computations properly + * within CoCoALib and map the result back to libpoly. + * + * We need to be careful to perform all computations in the proper polynomial + * rings, both to be correct and because CoCoALib explicitly requires it. As we + * change the ring we are computing it all the time, we also need appropriate + * ring homomorphisms to map polynomials from one into the other. We first give + * a short overview of our approach, then describe the various polynomial rings + * that are used, and then discuss which rings are used where. + * + * Inputs: + * - (real) variables x_0, ..., x_n + * - real algebraic numbers a_0, ..., a_{n-1} with + * - defining polynomials p_0, ..., p_{n-1}; p_i from Q[x_i] + * - a polynomial q over all variables x_0, ..., x_n + * + * We first iteratively build the field extensions Q(a_0), Q(a_0, a_2) ... + * Instead of the extension field Q(a_0), we use the isomorphic quotient ring + * Q[x_0]/ and recursively extend it with a_1, etc, in the same way. Doing + * this recursive construction naively fails: (Q[x_0]/)[x_1]/ is not + * necessarily a proper field as p_1 (though a minimal polynomial in Q[x_1]) may + * factor over Q[x_0]/. Consider p_0 = x_0*x_0-2 and p_1 = + * x_1*x_1*x_1*x_1-2 as an example, where p_1 factors into + * (x_1*x_1-x_0)*(x_1*x_1+x_0) over Q[x_0]/. We overcome this by explicitly + * computing this factorization and using the factor that vanishes over {x_0 -> + * a_0, x_1 -> a_1 } as the minimal polynomial of a_1 over Q[x_0]/. + * + * After we have built the field extensions in that way, we iteratively push q + * through the field extensions, each one extended to a polynomial ring over all + * x_0, ..., x_n. When in the k'th field extension, we check whether the k'th + * minimal polynomial divides q. If so, q would vanish in the next step and we + * instead set q = q/p_{k}. Only then we map q into K_{k+1}. + * + * Eventually, we end up with q in Q(a_0, ..., a_{n-1})[x_n]. This polynomial is + * univariate conceptually, and we want to compute its roots. However, it is not + * technically univariate and we need to make it so. We can do this by computing + * the Gröbner basis of the q and all minimal polynomials p_i with an + * elimination order with x_n at the bottom over Q[x_0, ..., x_n]. + * We then collect the polynomials + * that are univariate in x_n from the Gröbner basis. We can show that the roots + * of these polynomials are a superset of the roots we are looking for. + * + * + * To implement all that, we construct the following polynomial rings: + * - K_i: K_0 = Q, K_{i+1} = K_{i}[x_i]/ (with p_i reduced w.r.t. K_i) + * - R_i = K_i[x_i] + * - J_i = K_i[x_i, ..., x_n] = R_i[x_{i+1}, ..., x_n] + * + * While p_i conceptually live in Q[x_i], we immediately convert them from + * libpoly into R_i. We then factor it there, obtaining the actual minimal + * polynomial p_i that we use to construct K_{i+1}. We do this to construct all + * K_i and R_i. We then reduce q, initially in Q[x_0, ..., x_n] = J_0. We check + * in J_i whether p_i divides q (and if so divide q by p_i). To do + * this, we need to embed p_i into J_i. We then + * map q from J_i to J_{i+1}. While obvious in theory, this is somewhat tricky + * in practice as J_i and J_{i+1} have no direct relationship. + * Finally, we need to push all p_i and the final q back into J_0 = Q[x_0, ..., + * x_n] to compute the Gröbner basis. + * + * We thus furthermore store the following ring homomorphisms: + * - phom_i: R_i -> J_i (canonical embedding) + * - qhom_i: J_i -> J_{i+1} (hand-crafted homomorphism) + * + * We can sometimes avoid this construction for individual variables, i.e., if + * the assignment for x_i already lives (algebraically) in K_i. This can be the + * case if a_i is rational; in general, we check whether the vanishing factor + * of p_i is linear. If so, it has the form x_i-r where is some term in lower + * variables. We store r as the "direct assignment" in d_direct[i] and use it + * to directly replace x_i when appropriate. Also, we have K_i = K_{i-1}. + * + */ +struct LazardEvaluationState +{ + CoCoA::GlobalManager d_gm; + static std::unique_ptr 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. + */ + std::map, poly::Variable> d_varCoCoA; + + /** + * The minimal polynomials p_i used for constructing d_K. + * If a variable x_i has a rational assignment, p_i holds no value (i.e. + * d_p[i] == CoCoA::RingElem()). + */ + std::vector d_p; + + /** + * The sequence of extension fields. + * K_0 = Q, K_{i+1} = K_i[x_i]/ + * Every K_i is a field. + */ + std::vector d_K = {CoCoA::RingQQ()}; + /** + * R_i = K_i[x_i] + * Every R_i is a univariate polynomial ring over the field K_i. + */ + std::vector d_R; + /** + * J_i = K_i[x_i, ..., x_n] + * All J_i are constructed with CoCoA::lex ordering, just to make sure that + * the Gröbner basis of J_0 is computed as necessary. + */ + std::vector d_J; + + /** + * Custom homomorphism from R_i to J_i. PolyAlgebraHom with + * Indets(R_i) = (x_i) --> (x_i) + */ + std::vector d_phom; + /** + * Custom homomorphism from J_i to J_{i+1} + * If assignment of x_i is rational a PolyAlgebraHom with + * Indets(J_i) = (x_i,...,x_n) --> (a_i,x_{i+1},...,x_n) + * Otherwise a PolyRingHom with: + * - CoeffHom: K_{i-1} --> R_{i-1} --> K_i + * - (x_i,...,x_n) --> (x_i,x_{i+1},...,x_n), x_i = Indet(R_{i-1}) + */ + std::vector d_qhom; + + /** + * The base ideal for the Gröbner basis we compute in the end. Contains all + * p_i pushed into J_0. + */ + std::vector d_GBBaseIdeal; + + /** + * The current assignment, used to identify the vanishing factor to construct + * K_i. + */ + poly::Assignment d_assignment; + /** + * The libpoly variables in proper order. Directly correspond to x_0,...,x_n. + */ + std::vector d_variables; + /** + * Direct assignments for variables x_i as polynomials in lower variables. + * If the assignment for x_i is no direct assignment, d_direct[i] holds no + * value. + */ + std::vector> d_direct; + + LazardEvaluationState() + { + if (!d_stats) + { + d_stats = std::make_unique(); + } + } + + /** + * 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 + * highest variable x_i known yet. It thus directly constructs p_i in R_i. + */ + 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; + } + + /** + * Checks whether the given CoCoA polynomial evaluates to zero over the + * current libpoly assignment. The polynomial should live over the current + * R_i. + */ + bool evaluatesToZero(const CoCoA::RingElem& cp) const + { + Assert(CoCoA::owner(cp) == d_R.back()); + poly::Polynomial pp = convert(cp); + return poly::evaluate_constraint(pp, d_assignment, poly::SignCondition::EQ); + } + + /** + * Maps p from J_i to J_{i-1}. There can be no suitable homomorphism, and we + * thus manually decompose p into its terms and reconstruct them in J_{i-1}. + * If a_{i-1} is rational, we know that the coefficient rings of J_i and + * J_{i-1} are identical (K_{i-1} and K_{i-2}, respectively). We can thus + * immediately use coefficients from J_i as coefficients in J_{i-1}. + * Otherwise, we map coefficients from K_{i-1} to their canonical + * representation in R_{i-1} and then use d_phom[i-1] to map those into + * J_{i-1}. Afterwards, we iterate over the power product of the term + * reconstruct it in J_{i-1}. + */ + CoCoA::RingElem pushDownJ(const CoCoA::RingElem& p, size_t i) const + { + Trace("cad::lazard") << "Push " << p << " from " << d_J[i] << " to " + << d_J[i - 1] << std::endl; + Assert(CoCoA::owner(p) == d_J[i]); + CoCoA::RingElem res(d_J[i - 1]); + for (CoCoA::SparsePolyIter it = CoCoA::BeginIter(p); !CoCoA::IsEnded(it); + ++it) + { + CoCoA::RingElem coeff = CoCoA::coeff(it); + Assert(CoCoA::owner(coeff) == d_K[i]); + if (d_direct[i - 1]) + { + Assert(CoCoA::CoeffRing(d_J[i]) == CoCoA::CoeffRing(d_J[i - 1])); + coeff = CoCoA::CoeffEmbeddingHom(d_J[i - 1])(coeff); + } + else + { + coeff = CoCoA::CanonicalRepr(coeff); + Assert(CoCoA::owner(coeff) == d_R[i - 1]); + coeff = d_phom[i - 1](coeff); + } + Assert(CoCoA::owner(coeff) == d_J[i - 1]); + auto pp = CoCoA::PP(it); + std::vector indets = CoCoA::IndetsIn(pp); + for (size_t k = 0; k < indets.size(); ++k) + { + long exp = CoCoA::exponent(pp, indets[k]); + auto ind = CoCoA::indet(d_J[i - 1], indets[k] + 1); + coeff *= CoCoA::power(ind, exp); + } + res += coeff; + } + return res; + } + + /** + * Uses pushDownJ repeatedly to map p from J_{i+1} to J_0. + * Is used to map the minimal polynomials p_i and the reduced polynomial q + * into J_0 to eventually compute the Gröbner basis. + */ + CoCoA::RingElem pushDownJ0(const CoCoA::RingElem& p, size_t i) const + { + CoCoA::RingElem res = p; + for (; i > 0; --i) + { + Trace("cad::lazard") << "Pushing " << p << " from J" << i << " to J" + << i - 1 << std::endl; + res = pushDownJ(res, i); + } + return res; + } + + /** + * Add the next R_i: + * - add variable x_i to d_variables + * - extract the variable name + * - construct R_i = K_i[x_i] + * - add new variable to d_varCoCoA + */ + void addR(const poly::Variable& var) + { + d_variables.emplace_back(var); + if (Trace.isOn("cad::lazard")) + { + std::string vname = lp_variable_db_get_name( + poly::Context::get_context().get_variable_db(), var.get_internal()); + d_R.emplace_back(CoCoA::NewPolyRing(d_K.back(), {CoCoA::symbol(vname)})); + } + else + { + d_R.emplace_back(CoCoA::NewPolyRing(d_K.back(), {CoCoA::NewSymbol()})); + } + Trace("cad::lazard") << "R" << d_R.size() - 1 << " = " << d_R.back() + << std::endl; + d_varCoCoA.emplace(std::make_pair(CoCoA::RingID(d_R.back()), 0), var); + } + + /** + * Add the next K_{i+1} from a minimal polynomial: + * - store dummy value in d_direct + * - store the minimal polynomial p_i in d_p + * - construct K_{i+1} = R_i/ + */ + void addK(const poly::Variable& var, const CoCoA::RingElem& p) + { + d_direct.emplace_back(); + d_p.emplace_back(p); + Trace("cad::lazard") << "p" << d_p.size() - 1 << " = " << d_p.back() + << std::endl; + d_K.emplace_back(CoCoA::NewQuotientRing(d_R.back(), CoCoA::ideal(p))); + Trace("cad::lazard") << "K" << d_K.size() - 1 << " = " << d_K.back() + << std::endl; + } + + /** + * Add the next K_{i+1} from a rational assignment: + * - store assignment a_i in d_direct + * - store a dummy minimal polynomial in d_p + * - construct K_{i+1} as copy of K_i + */ + void addKRational(const poly::Variable& var, const CoCoA::RingElem& r) + { + d_direct.emplace_back(r); + d_p.emplace_back(); + Trace("cad::lazard") << "x" << d_p.size() - 1 << " = " << r << std::endl; + d_K.emplace_back(d_K.back()); + Trace("cad::lazard") << "K" << d_K.size() - 1 << " = " << d_K.back() + << std::endl; + } + + /** + * Finish the whole construction by adding the free variable: + * - add R_n by calling addR(var) + * - 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) + * - fill d_GBBaseIdeal with p_i mapped to J_0 + */ + void addFreeVariable(const poly::Variable& var) + { + Trace("cad::lazard") << "Add free variable " << var << std::endl; + addR(var); + std::vector symbols; + for (size_t i = 0; i < d_R.size(); ++i) + { + symbols.emplace_back(CoCoA::symbols(d_R[i]).back()); + } + for (size_t i = 0; i < d_R.size(); ++i) + { + d_J.emplace_back(CoCoA::NewPolyRing(d_K[i], symbols, CoCoA::lex)); + Trace("cad::lazard") << "J" << d_J.size() - 1 << " = " << d_J.back() + << std::endl; + symbols.erase(symbols.begin()); + // R_i --> J_i + d_phom.emplace_back( + CoCoA::PolyAlgebraHom(d_R[i], d_J[i], {CoCoA::indet(d_J[i], 0)})); + Trace("cad::lazard") << "R" << i << " --> J" << i << ": " << d_phom.back() + << std::endl; + if (i > 0) + { + Trace("cad::lazard") + << "Constructing J" << i - 1 << " --> J" << i << ": " << std::endl; + Trace("cad::lazard") << "Constructing " << d_J[i - 1] << " --> " + << d_J[i] << ": " << std::endl; + if (d_direct[i - 1]) + { + Trace("cad::lazard") << "Using " << d_variables[i - 1] << " for " + << CoCoA::indet(d_J[i - 1], 0) << std::endl; + Assert(CoCoA::CoeffRing(d_J[i]) == CoCoA::owner(*d_direct[i - 1])); + std::vector indets = { + CoCoA::RingElem(d_J[i], *d_direct[i - 1])}; + for (size_t j = 0; j < d_R.size() - i; ++j) + { + indets.push_back(CoCoA::indet(d_J[i], j)); + } + d_qhom.emplace_back( + CoCoA::PolyAlgebraHom(d_J[i - 1], d_J[i], indets)); + } + else + { + // K_{i-1} --> R_{i-1} + auto K2R = CoCoA::CoeffEmbeddingHom(d_R[i - 1]); + Assert(CoCoA::domain(K2R) == d_K[i - 1]); + Assert(CoCoA::codomain(K2R) == d_R[i - 1]); + // R_{i-1} --> K_i + auto R2K = CoCoA::QuotientingHom(d_K[i]); + Assert(CoCoA::domain(R2K) == d_R[i - 1]); + Assert(CoCoA::codomain(R2K) == d_K[i]); + // K_i --> J_i + auto K2J = CoCoA::CoeffEmbeddingHom(d_J[i]); + Assert(CoCoA::domain(K2J) == d_K[i]); + Assert(CoCoA::codomain(K2J) == d_J[i]); + // J_{i-1} --> J_i, consisting of + // - a homomorphism for the coefficients + // - a mapping for the indets + // Constructs [phom_i(x_i), x_i+1, ..., x_n] + std::vector indets = { + K2J(R2K(CoCoA::indet(d_R[i - 1], 0)))}; + for (size_t j = 0; j < d_R.size() - i; ++j) + { + indets.push_back(CoCoA::indet(d_J[i], j)); + } + d_qhom.emplace_back( + CoCoA::PolyRingHom(d_J[i - 1], d_J[i], R2K(K2R), indets)); + } + Trace("cad::lazard") << "J" << i - 1 << " --> J" << i << ": " + << d_qhom.back() << std::endl; + } + } + 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_variables[i]); + } + + d_GBBaseIdeal.clear(); + for (size_t i = 0; i < d_p.size(); ++i) + { + if (d_direct[i]) continue; + Trace("cad::lazard") << "Apply " << d_phom[i] << " to " << d_p[i] + << " from " << CoCoA::owner(d_p[i]) << std::endl; + d_GBBaseIdeal.emplace_back(pushDownJ0(d_phom[i](d_p[i]), i)); + } + + Trace("cad::lazard") << "Finished construction" << std::endl + << *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("cad::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("cad::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); + } + + /** + * Now reduce the polynomial qpoly: + * - convert qpoly into J_0 and factor it + * - for every factor q: + * - for every x_i: + * - if a_i is rational: + * - while q[x_i -> a_i] == 0 + * - q = q / (x_i - a_i) + * - set q = q[x_i -> a_i] + * - otherwise + * - obtain tmp = phom_i(p_i) + * - while tmp divides q + * - q = q / tmp + * - embed q = qhom_i(q) + * - compute (reduced) GBasis(p_0, ..., p_{n-i}, q) + * - collect and convert basis elements univariate in the free variable + */ + std::vector reduce(const poly::Polynomial& qpoly) const + { + d_stats->d_evaluations++; + std::vector res; + Trace("cad::lazard") << "Reducing " << qpoly << std::endl; + auto input = convertQ(qpoly); + Assert(CoCoA::owner(input) == d_J[0]); + auto factorization = CoCoA::factor(input); + for (const auto& f : factorization.myFactors()) + { + Trace("cad::lazard") << "-> factor " << f << std::endl; + CoCoA::RingElem q = f; + for (size_t i = 0; i < d_J.size() - 1; ++i) + { + Trace("cad::lazard") << "i = " << i << std::endl; + if (d_direct[i]) + { + Trace("cad::lazard") + << "Substitute " << d_variables[i] << " = " << *d_direct[i] + << " into " << q << " from " << CoCoA::owner(q) << std::endl; + auto indets = CoCoA::indets(d_J[i]); + auto var = indets[0]; + Assert(CoCoA::CoeffRing(d_J[i]) == CoCoA::owner(*d_direct[i])); + indets[0] = CoCoA::RingElem(d_J[i], *d_direct[i]); + auto hom = CoCoA::PolyAlgebraHom(d_J[i], d_J[i], indets); + while (CoCoA::IsZero(hom(q))) + { + q = q / (var - indets[0]); + d_stats->d_reductions++; + } + // substitute x_i -> a_i + q = hom(q); + Trace("cad::lazard") + << "-> " << q << " from " << CoCoA::owner(q) << std::endl; + } + else + { + auto tmp = d_phom[i](d_p[i]); + while (CoCoA::IsDivisible(q, tmp)) + { + q /= tmp; + d_stats->d_reductions++; + } + } + q = d_qhom[i](q); + } + Trace("cad::lazard") << "-> reduced to " << q << std::endl; + Assert(CoCoA::owner(q) == d_J.back()); + std::vector ideal = d_GBBaseIdeal; + ideal.emplace_back(pushDownJ0(q, d_J.size() - 1)); + Trace("cad::lazard") << "-> ideal " << ideal << std::endl; + auto basis = CoCoA::ReducedGBasis(CoCoA::ideal(ideal)); + Trace("cad::lazard") << "-> basis " << basis << std::endl; + for (const auto& belem : basis) + { + Trace("cad::lazard") << "-> retrieved " << belem << std::endl; + auto pres = convert(belem); + Trace("cad::lazard") << "-> converted " << pres << std::endl; + // These checks are orthogonal! + if (poly::is_univariate(pres) + && poly::is_univariate_over_assignment(pres, d_assignment)) + { + res.emplace_back(pres); + } + } + } + return res; + } +}; + +std::ostream& operator<<(std::ostream& os, const LazardEvaluationState& state) +{ + for (size_t i = 0; i < state.d_K.size(); ++i) + { + os << "K" << i << " = " << state.d_K[i] << std::endl; + os << "R" << i << " = " << state.d_R[i] << std::endl; + os << "J" << i << " = " << state.d_J[i] << std::endl; + + os << "R" << i << " --> J" << i << ": " << state.d_phom[i] << std::endl; + if (i > 0) + { + os << "J" << (i - 1) << " --> J" << i << ": " << state.d_qhom[i - 1] + << std::endl; + } + } + os << "GBBaseIdeal: " << state.d_GBBaseIdeal << std::endl; + os << "Done" << std::endl; + return os; +} +std::unique_ptr LazardEvaluationState::d_stats; + +LazardEvaluation::LazardEvaluation() + : d_state(std::make_unique()) +{ +} + +LazardEvaluation::~LazardEvaluation() {} + +/** + * Add a new variable with real algebraic number: + * - add var = ran to the assignment + * - add the next R_i by calling addR(var) + * - if ran is actually rational: + * - obtain the rational and call addKRational() + * - otherwise: + * - convert the minimal polynomial and identify vanishing factor + * - add the next K_i with the vanishing factor by valling addK() + */ +void LazardEvaluation::add(const poly::Variable& var, const poly::Value& val) +{ + Trace("cad::lazard") << "Adding " << var << " -> " << val << std::endl; + try + { + d_state->d_assignment.set(var, val); + d_state->addR(var); + + std::optional rational; + poly::UPolynomial polymipo; + if (poly::is_algebraic_number(val)) + { + const poly::AlgebraicNumber& ran = poly::as_algebraic_number(val); + const poly::DyadicInterval& di = poly::get_isolating_interval(ran); + if (poly::is_point(di)) + { + rational = d_state->convert(poly::get_point(di)); + } + else + { + Trace("cad::lazard") << "\tis proper ran" << std::endl; + polymipo = poly::get_defining_polynomial(ran); + } + } + else + { + Assert(poly::is_dyadic_rational(val) || poly::is_integer(val) + || poly::is_rational(val)); + if (poly::is_dyadic_rational(val)) + { + rational = d_state->convert(poly::as_dyadic_rational(val)); + } + else if (poly::is_integer(val)) + { + rational = CoCoA::BigRat(d_state->convert(poly::as_integer(val)), 1); + } + else if (poly::is_rational(val)) + { + rational = d_state->convert(poly::as_rational(val)); + } + } + + if (rational) + { + d_state->addKRational(var, + CoCoA::RingElem(d_state->d_K.back(), *rational)); + d_state->d_stats->d_directAssignments++; + return; + } + Trace("cad::lazard") << "Got mipo " << polymipo << std::endl; + auto mipo = d_state->convertMiPo(polymipo, var); + Trace("cad::lazard") << "Factoring " << mipo << " from " + << CoCoA::owner(mipo) << std::endl; + auto factorization = CoCoA::factor(mipo); + Trace("cad::lazard") << "-> " << factorization << std::endl; + bool used_factor = false; + for (const auto& f : factorization.myFactors()) + { + if (d_state->evaluatesToZero(f)) + { + Assert(CoCoA::deg(f) > 0 && CoCoA::NumTerms(f) <= 2); + if (CoCoA::deg(f) == 1) + { + auto rat = -CoCoA::ConstantCoeff(f) / CoCoA::LC(f); + Trace("cad::lazard") << "Using linear factor " << f << " -> " << var + << " = " << rat << std::endl; + d_state->addKRational(var, rat); + d_state->d_stats->d_directAssignments++; + } + else + { + Trace("cad::lazard") << "Using nonlinear factor " << f << std::endl; + d_state->addK(var, f); + d_state->d_stats->d_ranAssignments++; + } + used_factor = true; + break; + } + else + { + Trace("cad::lazard") << "Skipping " << f << std::endl; + } + } + Assert(used_factor); + } + catch (CoCoA::ErrorInfo& e) + { + e.myOutputSelf(std::cerr); + throw; + } +} + +void LazardEvaluation::addFreeVariable(const poly::Variable& var) +{ + try + { + d_state->addFreeVariable(var); + } + catch (CoCoA::ErrorInfo& e) + { + e.myOutputSelf(std::cerr); + throw; + } +} + +std::vector LazardEvaluation::reducePolynomial( + const poly::Polynomial& p) const +{ + try + { + return d_state->reduce(p); + } + catch (CoCoA::ErrorInfo& e) + { + e.myOutputSelf(std::cerr); + throw; + } + return {p}; +} + +/** + * Compute the infeasible regions of the given polynomial according to a sign + * condition. We first reduce the polynomial and isolate the real roots of every + * resulting polynomial. We store all roots (except for -infty, +infty and none) + * in a set. Then, we transform the set of roots into a list of infeasible + * regions by generating intervals between -infty and the first root, in between + * every two consecutive roots and between the last root and +infty. While doing + * this, we only keep those intervals that are actually infeasible for the + * original polynomial q over the partial assignment. Finally, we go over the + * intervals and aggregate consecutive intervals that connect. + */ +std::vector LazardEvaluation::infeasibleRegions( + const poly::Polynomial& q, poly::SignCondition sc) const +{ + poly::Assignment a; + std::set roots; + // reduce q to a set of reduced polynomials p + for (const auto& p : reducePolynomial(q)) + { + // collect all real roots except for -infty, none, +infty + Trace("cad::lazard") << "Isolating roots of " << p << std::endl; + Assert(poly::is_univariate(p) && poly::is_univariate_over_assignment(p, a)); + std::vector proots = poly::isolate_real_roots(p, a); + for (const auto& r : proots) + { + if (poly::is_minus_infinity(r)) continue; + if (poly::is_none(r)) continue; + if (poly::is_plus_infinity(r)) continue; + roots.insert(r); + } + } + + // generate all intervals + // (-infty,root_0), [root_0], (root_0,root_1), ..., [root_m], (root_m,+infty) + // if q is true over d_assignment x interval (represented by a sample) + std::vector res; + poly::Value last = poly::Value::minus_infty(); + for (const auto& r : roots) + { + poly::Value sample = poly::value_between(last, true, r, true); + d_state->d_assignment.set(d_state->d_variables.back(), sample); + if (!poly::evaluate_constraint(q, d_state->d_assignment, sc)) + { + res.emplace_back(last, true, r, true); + } + d_state->d_assignment.set(d_state->d_variables.back(), r); + if (!poly::evaluate_constraint(q, d_state->d_assignment, sc)) + { + res.emplace_back(r); + } + last = r; + } + poly::Value sample = + poly::value_between(last, true, poly::Value::plus_infty(), true); + d_state->d_assignment.set(d_state->d_variables.back(), sample); + if (!poly::evaluate_constraint(q, d_state->d_assignment, sc)) + { + res.emplace_back(last, true, poly::Value::plus_infty(), true); + } + // clean up assignment + d_state->d_assignment.unset(d_state->d_variables.back()); + + Trace("cad::lazard") << "Shrinking:" << std::endl; + for (const auto& i : res) + { + Trace("cad::lazard") << "-> " << i << std::endl; + } + std::vector combined; + if (res.empty()) + { + // nothing to do if there are no intervals to start with + // return combined to simplify return value optimization + return combined; + } + for (size_t i = 0; i < res.size() - 1; ++i) + { + // Invariant: the intervals do not overlap. Check for our own sanity. + Assert(poly::get_upper(res[i]) <= poly::get_lower(res[i + 1])); + if (poly::get_upper_open(res[i]) && poly::get_lower_open(res[i + 1])) + { + // does not connect, both are open + combined.emplace_back(res[i]); + continue; + } + if (poly::get_upper(res[i]) != poly::get_lower(res[i + 1])) + { + // does not connect, there is some space in between + combined.emplace_back(res[i]); + continue; + } + // combine them into res[i+1], do not copy res[i] over to combined + Trace("cad::lazard") << "Combine " << res[i] << " and " << res[i + 1] + << std::endl; + Assert(poly::get_lower(res[i]) <= poly::get_lower(res[i + 1])); + res[i + 1].set_lower(poly::get_lower(res[i]), poly::get_lower_open(res[i])); + } + + // always use the last one, it is never dropped + combined.emplace_back(res.back()); + Trace("cad::lazard") << "To:" << std::endl; + for (const auto& i : combined) + { + Trace("cad::lazard") << "-> " << i << std::endl; + } + return combined; +} + +} // namespace cvc5::theory::arith::nl::cad + +#else namespace cvc5::theory::arith::nl::cad { @@ -48,3 +975,4 @@ std::vector LazardEvaluation::infeasibleRegions( } // namespace cvc5::theory::arith::nl::cad #endif +#endif -- 2.30.2