(cad-solver) Fix square-free-basis computation (#5085)

This PR fixes a subtle issue when making the polynomials of two subsequent (overlapping) intervals relative square-free. We now make sure that the resulting polynomials are ignored (if constant) or pushed to the lower dimension (if lower main variable). Also, we now appropriately update the main polynomials of the corresponding intervals.

5 files changed, 115 insertions, 31 deletions

diff --git a/src/theory/arith/nl/cad/cdcac.cpp b/src/theory/arith/nl/cad/cdcac.cpp
index 0dcf7f7a7..f1ae77e2e 100644
--- a/src/theory/arith/nl/cad/cdcac.cpp
+++ b/src/theory/arith/nl/cad/cdcac.cpp
@@ -180,6 +180,12 @@ std::vector<poly::Polynomial> CDCAC::constructCharacterization(
   Trace("cdcac") << "Constructing characterization now" << std::endl;
   std::vector<poly::Polynomial> res;
 
+
+  for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
+  {
+    cad::makeFinestSquareFreeBasis(intervals[i], intervals[i + 1]);
+  }
+
   for (const auto& i : intervals)
   {
     Trace("cdcac") << "Considering " << i.d_interval << std::endl;
@@ -229,8 +235,6 @@ std::vector<poly::Polynomial> CDCAC::constructCharacterization(
   for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
   {
     // Add resultants of consecutive intervals.
-    cad::makeFinestSquareFreeBasis(intervals[i].d_upperPolys,
-                                   intervals[i + 1].d_lowerPolys);
     for (const auto& p : intervals[i].d_upperPolys)
     {
       for (const auto& q : intervals[i + 1].d_lowerPolys)
diff --git a/src/theory/arith/nl/cad/cdcac_utils.cpp b/src/theory/arith/nl/cad/cdcac_utils.cpp
index a9b4c6128..23eaff033 100644
--- a/src/theory/arith/nl/cad/cdcac_utils.cpp
+++ b/src/theory/arith/nl/cad/cdcac_utils.cpp
@@ -18,6 +18,8 @@
 
 #ifdef CVC4_POLY_IMP
 
+#include "theory/arith/nl/cad/projections.h"
+
 namespace CVC4 {
 namespace theory {
 namespace arith {
@@ -264,6 +266,106 @@ bool sampleOutside(const std::vector<CACInterval>& infeasible, Value& sample)
   return false;
 }
 
+namespace {
+/**
+ * Replace a polynomial at polys[id] with the given pair of polynomials.
+ * Also update d_mainPolys in the given interval accordingly.
+ * If one of the factors (from the pair) is from a lower level (has a different
+ * main variable), push this factor to the d_downPolys.
+ * The first factor needs to be a proper polynomial (!is_constant(subst.first)),
+ * but the second factor may be anything.
+ */
+void replace_polynomial(std::vector<poly::Polynomial>& polys,
+                        std::size_t id,
+                        std::pair<poly::Polynomial, poly::Polynomial> subst,
+                        CACInterval& interval)
+{
+  Assert(polys[id] == subst.first * subst.second);
+  Assert(!is_constant(subst.first));
+  // Whether polys[id] has already been replaced
+  bool replaced = false;
+  poly::Variable var = main_variable(polys[id]);
+  // The position within interval.d_mainPolys to be replaced.
+  auto mit = std::find(
+      interval.d_mainPolys.begin(), interval.d_mainPolys.end(), polys[id]);
+  if (main_variable(subst.first) == var)
+  {
+    // Replace in polys[id] and *mit
+    polys[id] = subst.first;
+    if (mit != interval.d_mainPolys.end())
+    {
+      *mit = subst.first;
+    }
+    replaced = true;
+  }
+  else
+  {
+    // Push to d_downPolys
+    interval.d_downPolys.emplace_back(subst.first);
+  }
+  // Skip constant poly
+  if (!is_constant(subst.second))
+  {
+    if (main_variable(subst.second) == var)
+    {
+      if (replaced)
+      {
+        // Append to polys and d_mainPolys
+        polys.emplace_back(subst.second);
+        interval.d_mainPolys.emplace_back(subst.second);
+      }
+      else
+      {
+        // Replace in polys[id] and *mit
+        polys[id] = subst.second;
+        if (mit != interval.d_mainPolys.end())
+        {
+          *mit = subst.second;
+        }
+        replaced = true;
+      }
+    }
+    else
+    {
+      // Push to d_downPolys
+      interval.d_downPolys.emplace_back(subst.second);
+    }
+  }
+  Assert(replaced)
+      << "At least one of the factors should have the main variable";
+}
+}  // namespace
+
+void makeFinestSquareFreeBasis(CACInterval& lhs, CACInterval& rhs)
+{
+  auto& l = lhs.d_upperPolys;
+  auto& r = rhs.d_lowerPolys;
+  if (l.empty()) return;
+  for (std::size_t i = 0, ln = l.size(); i < ln; ++i)
+  {
+    for (std::size_t j = 0, rn = r.size(); j < rn; ++j)
+    {
+      // All main variables should be the same
+      Assert(main_variable(l[i]) == main_variable(r[j]));
+      if (l[i] == r[j]) continue;
+      Polynomial g = gcd(l[i], r[j]);
+      if (!is_constant(g))
+      {
+        auto newl = div(l[i], g);
+        auto newr = div(r[j], g);
+        replace_polynomial(l, i, std::make_pair(g, newl), lhs);
+        replace_polynomial(r, j, std::make_pair(g, newr), rhs);
+      }
+    }
+  }
+  reduceProjectionPolynomials(l);
+  reduceProjectionPolynomials(r);
+  reduceProjectionPolynomials(lhs.d_mainPolys);
+  reduceProjectionPolynomials(rhs.d_mainPolys);
+  reduceProjectionPolynomials(lhs.d_downPolys);
+  reduceProjectionPolynomials(rhs.d_downPolys);
+}
+
 }  // namespace cad
 }  // namespace nl
 }  // namespace arith
diff --git a/src/theory/arith/nl/cad/cdcac_utils.h b/src/theory/arith/nl/cad/cdcac_utils.h
index ed1c0d1c2..c0f800170 100644
--- a/src/theory/arith/nl/cad/cdcac_utils.h
+++ b/src/theory/arith/nl/cad/cdcac_utils.h
@@ -95,6 +95,13 @@ std::vector<Node> collectConstraints(const std::vector<CACInterval>& intervals);
 bool sampleOutside(const std::vector<CACInterval>& infeasible,
                    poly::Value& sample);
 
+/**
+ * Compute the finest square of the upper polynomials of lhs and the lower
+ * polynomials of rhs. Also pushes reduced polynomials to lower level if
+ * necessary.
+ */
+void makeFinestSquareFreeBasis(CACInterval& lhs, CACInterval& rhs);
+
 }  // namespace cad
 }  // namespace nl
 }  // namespace arith
diff --git a/src/theory/arith/nl/cad/projections.cpp b/src/theory/arith/nl/cad/projections.cpp
index 488a1f1c6..276494afd 100644
--- a/src/theory/arith/nl/cad/projections.cpp
+++ b/src/theory/arith/nl/cad/projections.cpp
@@ -70,28 +70,6 @@ void makeFinestSquareFreeBasis(std::vector<Polynomial>& polys)
   reduceProjectionPolynomials(polys);
 }
 
-void makeFinestSquareFreeBasis(std::vector<poly::Polynomial>& lhs,
-                               std::vector<poly::Polynomial>& rhs)
-{
-  for (std::size_t i = 0, ln = lhs.size(); i < ln; ++i)
-  {
-    for (std::size_t j = 0, rn = rhs.size(); j < rn; ++j)
-    {
-      if (lhs[i] == rhs[j]) continue;
-      Polynomial g = gcd(lhs[i], rhs[j]);
-      if (!is_constant(g))
-      {
-        lhs[i] = div(lhs[i], g);
-        rhs[j] = div(rhs[j], g);
-        lhs.emplace_back(g);
-        rhs.emplace_back(g);
-      }
-    }
-  }
-  reduceProjectionPolynomials(lhs);
-  reduceProjectionPolynomials(rhs);
-}
-
 std::vector<Polynomial> projection_mccallum(
     const std::vector<Polynomial>& polys)
 {
diff --git a/src/theory/arith/nl/cad/projections.h b/src/theory/arith/nl/cad/projections.h
index c4a34ee56..afed5b1e9 100644
--- a/src/theory/arith/nl/cad/projections.h
+++ b/src/theory/arith/nl/cad/projections.h
@@ -53,13 +53,6 @@ void addPolynomials(std::vector<poly::Polynomial>& polys,
 void makeFinestSquareFreeBasis(std::vector<poly::Polynomial>& polys);
 
 /**
- * Ensure that two sets of polynomials are finest square-free basis relative to
- * each other.
- */
-void makeFinestSquareFreeBasis(std::vector<poly::Polynomial>& lhs,
-                               std::vector<poly::Polynomial>& rhs);
-
-/**
  * Computes McCallum's projection operator.
  */
 std::vector<poly::Polynomial> projectionMcCallum(