3 files changed, 157 insertions, 0 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 1edafb977..bec3a8345 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -391,6 +391,8 @@ libcvc5_add_sources(
   theory/arith/nl/cad/cdcac_utils.h
   theory/arith/nl/cad/constraints.cpp
   theory/arith/nl/cad/constraints.h
+  theory/arith/nl/cad/lazard_evaluation.cpp
+  theory/arith/nl/cad/lazard_evaluation.h
   theory/arith/nl/cad/projections.cpp
   theory/arith/nl/cad/projections.h
   theory/arith/nl/cad/proof_checker.cpp
diff --git a/src/theory/arith/nl/cad/lazard_evaluation.cpp b/src/theory/arith/nl/cad/lazard_evaluation.cpp
new file mode 100644
index 000000000..82b127ed0
--- /dev/null
+++ b/src/theory/arith/nl/cad/lazard_evaluation.cpp
@@ -0,0 +1,46 @@
+#include "theory/arith/nl/cad/lazard_evaluation.h"
+
+#include "base/check.h"
+#include "base/output.h"
+
+namespace cvc5::theory::arith::nl::cad {
+
+/**
+ * Do a very simple wrapper around the regular poly::infeasible_regions.
+ * Warn the user about doing this.
+ * This allows for a graceful fallback (albeit with a warning) if CoCoA is not
+ * available.
+ */
+struct LazardEvaluationState
+{
+  poly::Assignment d_assignment;
+};
+LazardEvaluation::LazardEvaluation()
+    : d_state(std::make_unique<LazardEvaluationState>())
+{
+}
+LazardEvaluation::~LazardEvaluation() {}
+
+void LazardEvaluation::add(const poly::Variable& var, const poly::Value& val)
+{
+  d_state->d_assignment.set(var, val);
+}
+
+void LazardEvaluation::addFreeVariable(const poly::Variable& var) {}
+
+std::vector<poly::Polynomial> LazardEvaluation::reducePolynomial(
+    const poly::Polynomial& p) const
+{
+  return {p};
+}
+std::vector<poly::Interval> LazardEvaluation::infeasibleRegions(
+    const poly::Polynomial& q, poly::SignCondition sc) const
+{
+  WarningOnce()
+      << "CAD::LazardEvaluation is disabled because CoCoA is not available. "
+         "Falling back to regular calculation of infeasible regions."
+      << std::endl;
+  return poly::infeasible_regions(q, d_state->d_assignment, sc);
+}
+
+}  // namespace cvc5::theory::arith::nl::cad
diff --git a/src/theory/arith/nl/cad/lazard_evaluation.h b/src/theory/arith/nl/cad/lazard_evaluation.h
new file mode 100644
index 000000000..0edb12010
--- /dev/null
+++ b/src/theory/arith/nl/cad/lazard_evaluation.h
@@ -0,0 +1,109 @@
+/******************************************************************************
+ * 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.
+ * ****************************************************************************
+ *
+ * Implements the CDCAC approach as described in
+ * https://arxiv.org/pdf/2003.05633.pdf.
+ */
+
+#include "cvc5_private.h"
+
+#ifndef CVC5__THEORY__ARITH__NL__CAD__LAZARD_EVALUATION_H
+#define CVC5__THEORY__ARITH__NL__CAD__LAZARD_EVALUATION_H
+
+#include <poly/polyxx.h>
+
+#include <memory>
+
+namespace cvc5::theory::arith::nl::cad {
+
+struct LazardEvaluationState;
+/**
+ * This class does the heavy lifting for the modified lifting procedure that is
+ * required for Lazard's projection.
+ *
+ * Let p \in Q[x_0, ..., x_n] a multivariate polynomial whose roots we want to
+ * isolate over the partial sample point A = [x_0 = a_0, ... x_{n-1} = a_{n-1}]
+ * where a_0, ... a_{n-1} are real algebraic numbers and x_n is the last free
+ * variable.
+ *
+ * The modified lifting procedure conceptually works as follows:
+ *
+ * for (x = a) in A:
+ *    while p[x // a] = 0:
+ *       p = p / (x - a)
+ *    p = p[x // a]
+ * return roots(p)
+ *
+ * As the assignment contains real algebraic numbers, though, we can not do any
+ * of the computations directly, as our polynomials only support coefficients
+ * from Z or Q, but not extensions (in the algebraic sense) thereof.
+ *
+ * Our approach is as follows:
+ * Let pk be the minimal polynomial for a_k.
+ * Instead of substituting p[x_k // a_k] we (canonically) embed p into the
+ * quotient ring Q[x_k]/<p_k> and recursively build a tower of such quotient
+ * rings that is isomorphic to nesting the corresponding field extensions
+ * Q(a_1)(a_2)... When we have done that, we obtain p that is reduced with
+ * respect to all minimal polynomials, but may still contain x_0,... x_{n-1}. To
+ * get rid of these, we compute a Gröbner basis of p and the minimal polynomials
+ * (using a suitable elimination order) and extract the polynomial in x_n. This
+ * polynomial has all roots (and possibly some more) that we are looking for.
+ * Instead of a Gröbner basis, we can also compute the iterated resultant as
+ * follows: Res(Res(p, p_{n-1}, x_{n-1}), p_{n-2}, x_{n-2})...
+ *
+ * Consider
+ * http://sunsite.informatik.rwth-aachen.de/Publications/AIB/2020/2020-04.pdf
+ * Section 2.5.1 for a full discussion.
+ *
+ * !!! WARNING !!!
+ * If CoCoALib is not available, this class will simply fall back to
+ * poly::infeasible_regions and issue a warning about this.
+ */
+class LazardEvaluation
+{
+ public:
+  LazardEvaluation();
+  ~LazardEvaluation();
+
+  /**
+   * Add the next assigned variable x_k = a_k to this construction.
+   */
+  void add(const poly::Variable& var, const poly::Value& val);
+  /**
+   * Finish by adding the free variable x_n.
+   */
+  void addFreeVariable(const poly::Variable& var);
+  /**
+   * Reduce the polynomial q. Compared to the above description, there may
+   * actually be multiple polynomials in the Gröbner basis and instead of
+   * loosing this knowledge and returning their product, we return them as a
+   * vector.
+   */
+  std::vector<poly::Polynomial> reducePolynomial(
+      const poly::Polynomial& q) const;
+
+  /**
+   * Compute the infeasible regions of q under the given sign condition.
+   * Uses reducePolynomial and then performs real root isolation on the
+   * resulting polynomials to obtain the intervals. Mimics
+   * poly::infeasible_regions, but uses Lazard's evaluation.
+   */
+  std::vector<poly::Interval> infeasibleRegions(const poly::Polynomial& q,
+                                                poly::SignCondition sc) const;
+
+ private:
+  std::unique_ptr<LazardEvaluationState> d_state;
+};
+
+}  // namespace cvc5::theory::arith::nl::cad
+
+#endif
+\ No newline at end of file