(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.

diff --git a/src/theory/arith/nl/cad/cdcac.cpp b/src/theory/arith/nl/cad/cdcac.cpp
index 0dcf7f7a7bf3f6dccc622d1b54b630a6a528e334..f1ae77e2e21cc41f7df60cf4f92fd7f1c62e37c5 100644 (file)
--- 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 a9b4c6128769681c614840b25c1cdb5c63c2d476..23eaff0335065a5e944229c3b4e400bfd25984ad 100644 (file)
--- 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 ed1c0d1c2fed2f2e5fda3e1c8033dc04e9b21dc6..c0f80017001719e51a3680385422b7221bf7a077 100644 (file)
--- 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 488a1f1c62a57f244cccd5a70c1d048dc5820865..276494afd806a67a4bd78f6d3d081184cfa96a8d 100644 (file)
--- 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 c4a34ee562f93f565b5541112bf10852a355c1b9..afed5b1e9a06ee2b6a71cfd66edfdc7ae246bd24 100644 (file)
--- a/src/theory/arith/nl/cad/projections.h
+++ b/src/theory/arith/nl/cad/projections.h
@@ -52,13 +52,6 @@ void addPolynomials(std::vector<poly::Polynomial>& polys,
 /** Make a set of polynomials a finest square-free basis. */
 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.
  */