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.

1 files changed, 928 insertions, 0 deletions

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 <CoCoA/library.H>
+
+#include <optional>
+
+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]/<p_0> and recursively extend it with a_1, etc, in the same way. Doing
+ * this recursive construction naively fails: (Q[x_0]/<p_0>)[x_1]/<p_1> is not
+ * necessarily a proper field as p_1 (though a minimal polynomial in Q[x_1]) may
+ * factor over Q[x_0]/<p_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]/<p_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]/<p_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]/<p_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<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<poly::Variable, CoCoA::RingElem> 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<std::pair<long, size_t>, 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<CoCoA::RingElem> d_p;
+
+  /**
+   * The sequence of extension fields.
+   * K_0 = Q, K_{i+1} = K_i[x_i]/<p_i>
+   * Every K_i is a field.
+   */
+  std::vector<CoCoA::ring> 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<CoCoA::ring> 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<CoCoA::ring> d_J;
+
+  /**
+   * Custom homomorphism from R_i to J_i. PolyAlgebraHom with
+   * Indets(R_i) = (x_i) --> (x_i)
+   */
+  std::vector<CoCoA::RingHom> 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<CoCoA::RingHom> 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<CoCoA::RingElem> 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<poly::Variable> 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<std::optional<CoCoA::RingElem>> d_direct;
+
+  LazardEvaluationState()
+  {
+    if (!d_stats)
+    {
+      d_stats = std::make_unique<LazardEvaluationStats>();
+    }
+  }
+
+  /**
+   * 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<poly::Integer> 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<long> 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/<p_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<CoCoA::symbol> 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<CoCoA::RingElem> 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<CoCoA::RingElem> 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<CoCoAPolyConstructor*>(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<long> 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<poly::Polynomial> reduce(const poly::Polynomial& qpoly) const
+  {
+    d_stats->d_evaluations++;
+    std::vector<poly::Polynomial> 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<CoCoA::RingElem> 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<LazardEvaluationStats> LazardEvaluationState::d_stats;
+
+LazardEvaluation::LazardEvaluation()
+    : d_state(std::make_unique<LazardEvaluationState>())
+{
+}
+
+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<CoCoA::BigRat> 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<poly::Polynomial> 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<poly::Interval> LazardEvaluation::infeasibleRegions(
+    const poly::Polynomial& q, poly::SignCondition sc) const
+{
+  poly::Assignment a;
+  std::set<poly::Value> 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<poly::Value> 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<poly::Interval> 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<poly::Interval> 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<poly::Interval> LazardEvaluation::infeasibleRegions(
 }  // namespace cvc5::theory::arith::nl::cad
 
 #endif
+#endif