Make sure polynomials are properly factorized in nl-cad (#5733)

CAD theory (used in nl-cad) requires that polynomials are properly factorized (a finest square-free basis). This PR replaces usage of raw std::vector by a new wrapper PolyVector that ensures proper factorization whenever a polynomial is added. This fixes one piece of code that omitted factorization, leading to soundness issues as in #5726.
Fixes #5726.

6 files changed, 122 insertions, 116 deletions

diff --git a/src/theory/arith/nl/cad/cdcac.cpp b/src/theory/arith/nl/cad/cdcac.cpp
index 57a8b51df..e17e946cd 100644
--- a/src/theory/arith/nl/cad/cdcac.cpp
+++ b/src/theory/arith/nl/cad/cdcac.cpp
@@ -39,16 +39,6 @@ namespace arith {
 namespace nl {
 namespace cad {
 
-namespace {
-/** Removed duplicates from a vector. */
-template <typename T>
-void removeDuplicates(std::vector<T>& v)
-{
-  std::sort(v.begin(), v.end());
-  v.erase(std::unique(v.begin(), v.end()), v.end());
-}
-}  // namespace
-
 CDCAC::CDCAC() {}
 
 CDCAC::CDCAC(const std::vector<poly::Variable>& ordering)
@@ -125,10 +115,11 @@ std::vector<CACInterval> CDCAC::getUnsatIntervals(
     for (const auto& i : intervals)
     {
       Trace("cdcac") << "-> " << i << std::endl;
-      std::vector<poly::Polynomial> l, u, m, d;
-      if (!is_minus_infinity(get_lower(i))) l.emplace_back(p);
-      if (!is_plus_infinity(get_upper(i))) u.emplace_back(p);
-      m.emplace_back(p);
+      PolyVector l, u, m, d;
+      m.add(p);
+      m.pushDownPolys(d, d_variableOrdering[cur_variable]);
+      if (!is_minus_infinity(get_lower(i))) l = m;
+      if (!is_plus_infinity(get_upper(i))) u = m;
       res.emplace_back(CACInterval{i, l, u, m, d, {n}});
     }
   }
@@ -158,15 +149,14 @@ bool CDCAC::sampleOutsideWithInitial(const std::vector<CACInterval>& infeasible,
   return sampleOutside(infeasible, sample);
 }
 
-std::vector<poly::Polynomial> CDCAC::requiredCoefficients(
-    const poly::Polynomial& p) const
+PolyVector CDCAC::requiredCoefficients(const poly::Polynomial& p) const
 {
-  std::vector<poly::Polynomial> res;
+  PolyVector res;
   for (long deg = degree(p); deg >= 0; --deg)
   {
     auto coeff = coefficient(p, deg);
     if (lp_polynomial_is_constant(coeff.get_internal())) break;
-    res.emplace_back(coeff);
+    res.add(coeff);
     if (evaluate_constraint(coeff, d_assignment, poly::SignCondition::NE))
     {
       break;
@@ -175,13 +165,11 @@ std::vector<poly::Polynomial> CDCAC::requiredCoefficients(
   return res;
 }
 
-std::vector<poly::Polynomial> CDCAC::constructCharacterization(
-    std::vector<CACInterval>& intervals)
+PolyVector CDCAC::constructCharacterization(std::vector<CACInterval>& intervals)
 {
   Assert(!intervals.empty()) << "A covering can not be empty";
   Trace("cdcac") << "Constructing characterization now" << std::endl;
-  std::vector<poly::Polynomial> res;
-
+  PolyVector res;
 
   for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
   {
@@ -198,20 +186,20 @@ std::vector<poly::Polynomial> CDCAC::constructCharacterization(
     for (const auto& p : i.d_downPolys)
     {
       // Add all polynomial from lower levels.
-      addPolynomial(res, p);
+      res.add(p);
     }
     for (const auto& p : i.d_mainPolys)
     {
       Trace("cdcac") << "Discriminant of " << p << " -> " << discriminant(p)
                      << std::endl;
       // Add all discriminants
-      addPolynomial(res, discriminant(p));
+      res.add(discriminant(p));
 
       for (const auto& q : requiredCoefficients(p))
       {
         // Add all required coefficients
         Trace("cdcac") << "Coeff of " << p << " -> " << q << std::endl;
-        addPolynomial(res, q);
+        res.add(q);
       }
       for (const auto& q : i.d_lowerPolys)
       {
@@ -220,7 +208,7 @@ std::vector<poly::Polynomial> CDCAC::constructCharacterization(
         if (!hasRootBelow(q, get_lower(i.d_interval))) continue;
         Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
                        << resultant(p, q) << std::endl;
-        addPolynomial(res, resultant(p, q));
+        res.add(resultant(p, q));
       }
       for (const auto& q : i.d_upperPolys)
       {
@@ -229,7 +217,7 @@ std::vector<poly::Polynomial> CDCAC::constructCharacterization(
         if (!hasRootAbove(q, get_upper(i.d_interval))) continue;
         Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
                        << resultant(p, q) << std::endl;
-        addPolynomial(res, resultant(p, q));
+        res.add(resultant(p, q));
       }
     }
   }
@@ -243,39 +231,34 @@ std::vector<poly::Polynomial> CDCAC::constructCharacterization(
       {
         Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
                        << resultant(p, q) << std::endl;
-        addPolynomial(res, resultant(p, q));
+        res.add(resultant(p, q));
       }
     }
   }
 
-  removeDuplicates(res);
-  makeFinestSquareFreeBasis(res);
+  res.reduce();
+  res.makeFinestSquareFreeBasis();
 
   return res;
 }
 
 CACInterval CDCAC::intervalFromCharacterization(
-    const std::vector<poly::Polynomial>& characterization,
+    const PolyVector& characterization,
     std::size_t cur_variable,
     const poly::Value& sample)
 {
-  std::vector<poly::Polynomial> l;
-  std::vector<poly::Polynomial> u;
-  std::vector<poly::Polynomial> m;
-  std::vector<poly::Polynomial> d;
+  PolyVector l;
+  PolyVector u;
+  PolyVector m;
+  PolyVector d;
 
   for (const auto& p : characterization)
   {
-    // Add polynomials to either main or down
-    if (main_variable(p) == d_variableOrdering[cur_variable])
-    {
-      m.emplace_back(p);
-    }
-    else
-    {
-      d.emplace_back(p);
-    }
+    // Add polynomials to main
+    m.add(p);
   }
+  // Push lower-dimensional polys to down
+  m.pushDownPolys(d, d_variableOrdering[cur_variable]);
 
   // Collect -oo, all roots, oo
   std::vector<poly::Value> roots;
@@ -316,7 +299,7 @@ CACInterval CDCAC::intervalFromCharacterization(
     {
       if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
       {
-        l.emplace_back(p);
+        l.add(p, true);
       }
     }
     d_assignment.unset(d_variableOrdering[cur_variable]);
@@ -329,7 +312,7 @@ CACInterval CDCAC::intervalFromCharacterization(
     {
       if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
       {
-        u.emplace_back(p);
+        u.add(p, true);
       }
     }
     d_assignment.unset(d_variableOrdering[cur_variable]);
diff --git a/src/theory/arith/nl/cad/cdcac.h b/src/theory/arith/nl/cad/cdcac.h
index a2e7ae682..1b01d0bf6 100644
--- a/src/theory/arith/nl/cad/cdcac.h
+++ b/src/theory/arith/nl/cad/cdcac.h
@@ -103,25 +103,22 @@ class CDCAC
    * Collects the coefficients required for projection from the given
    * polynomial. Implements Algorithm 6.
    */
-  std::vector<poly::Polynomial> requiredCoefficients(
-      const poly::Polynomial& p) const;
+  PolyVector requiredCoefficients(const poly::Polynomial& p) const;
 
   /**
    * Constructs a characterization of the given covering.
    * A characterization contains polynomials whose roots bound the region around
    * the current assignment. Implements Algorithm 4.
    */
-  std::vector<poly::Polynomial> constructCharacterization(
-      std::vector<CACInterval>& intervals);
+  PolyVector constructCharacterization(std::vector<CACInterval>& intervals);
 
   /**
    * Constructs an infeasible interval from a characterization.
    * Implements Algorithm 5.
    */
-  CACInterval intervalFromCharacterization(
-      const std::vector<poly::Polynomial>& characterization,
-      std::size_t cur_variable,
-      const poly::Value& sample);
+  CACInterval intervalFromCharacterization(const PolyVector& characterization,
+                                           std::size_t cur_variable,
+                                           const poly::Value& sample);
 
   /**
    * Main method that checks for the satisfiability of the constraints.
diff --git a/src/theory/arith/nl/cad/cdcac_utils.cpp b/src/theory/arith/nl/cad/cdcac_utils.cpp
index f36ec775f..e58ba3730 100644
--- a/src/theory/arith/nl/cad/cdcac_utils.cpp
+++ b/src/theory/arith/nl/cad/cdcac_utils.cpp
@@ -275,7 +275,7 @@ namespace {
  * 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,
+void replace_polynomial(PolyVector& polys,
                         std::size_t id,
                         std::pair<poly::Polynomial, poly::Polynomial> subst,
                         CACInterval& interval)
@@ -301,7 +301,7 @@ void replace_polynomial(std::vector<poly::Polynomial>& polys,
   else
   {
     // Push to d_downPolys
-    interval.d_downPolys.emplace_back(subst.first);
+    interval.d_downPolys.add(subst.first);
   }
   // Skip constant poly
   if (!is_constant(subst.second))
@@ -311,8 +311,8 @@ void replace_polynomial(std::vector<poly::Polynomial>& polys,
       if (replaced)
       {
         // Append to polys and d_mainPolys
-        polys.emplace_back(subst.second);
-        interval.d_mainPolys.emplace_back(subst.second);
+        polys.add(subst.second);
+        interval.d_mainPolys.add(subst.second);
       }
       else
       {
@@ -328,7 +328,7 @@ void replace_polynomial(std::vector<poly::Polynomial>& polys,
     else
     {
       // Push to d_downPolys
-      interval.d_downPolys.emplace_back(subst.second);
+      interval.d_downPolys.add(subst.second);
     }
   }
   Assert(replaced)
@@ -358,12 +358,12 @@ void makeFinestSquareFreeBasis(CACInterval& lhs, CACInterval& rhs)
       }
     }
   }
-  reduceProjectionPolynomials(l);
-  reduceProjectionPolynomials(r);
-  reduceProjectionPolynomials(lhs.d_mainPolys);
-  reduceProjectionPolynomials(rhs.d_mainPolys);
-  reduceProjectionPolynomials(lhs.d_downPolys);
-  reduceProjectionPolynomials(rhs.d_downPolys);
+  l.reduce();
+  r.reduce();
+  lhs.d_mainPolys.reduce();
+  rhs.d_mainPolys.reduce();
+  lhs.d_downPolys.reduce();
+  rhs.d_downPolys.reduce();
 }
 
 }  // namespace cad
diff --git a/src/theory/arith/nl/cad/cdcac_utils.h b/src/theory/arith/nl/cad/cdcac_utils.h
index 43bef32aa..335735574 100644
--- a/src/theory/arith/nl/cad/cdcac_utils.h
+++ b/src/theory/arith/nl/cad/cdcac_utils.h
@@ -28,6 +28,7 @@
 #include <vector>
 
 #include "expr/node.h"
+#include "theory/arith/nl/cad/projections.h"
 
 namespace CVC4 {
 namespace theory {
@@ -51,13 +52,13 @@ struct CACInterval
   /** The actual interval. */
   poly::Interval d_interval;
   /** The polynomials characterizing the lower bound. */
-  std::vector<poly::Polynomial> d_lowerPolys;
+  PolyVector d_lowerPolys;
   /** The polynomials characterizing the upper bound. */
-  std::vector<poly::Polynomial> d_upperPolys;
+  PolyVector d_upperPolys;
   /** The characterizing polynomials in the main variable. */
-  std::vector<poly::Polynomial> d_mainPolys;
+  PolyVector d_mainPolys;
   /** The characterizing polynomials in lower variables. */
-  std::vector<poly::Polynomial> d_downPolys;
+  PolyVector d_downPolys;
   /** The constraints used to derive this interval. */
   std::vector<Node> d_origins;
 };
diff --git a/src/theory/arith/nl/cad/projections.cpp b/src/theory/arith/nl/cad/projections.cpp
index 162e9f7be..8d33bf794 100644
--- a/src/theory/arith/nl/cad/projections.cpp
+++ b/src/theory/arith/nl/cad/projections.cpp
@@ -18,6 +18,8 @@
 
 #ifdef CVC4_POLY_IMP
 
+#include "base/check.h"
+
 namespace CVC4 {
 namespace theory {
 namespace arith {
@@ -26,72 +28,77 @@ namespace cad {
 
 using namespace poly;
 
-void reduceProjectionPolynomials(std::vector<Polynomial>& polys)
+void PolyVector::add(const poly::Polynomial& poly, bool assertMain)
 {
-  std::sort(polys.begin(), polys.end());
-  auto it = std::unique(polys.begin(), polys.end());
-  polys.erase(it, polys.end());
-}
-
-void addPolynomial(std::vector<Polynomial>& polys, const Polynomial& poly)
-{
-  for (const auto& p : square_free_factors(poly))
+  for (const auto& p : poly::square_free_factors(poly))
   {
-    if (is_constant(p)) continue;
-    polys.emplace_back(p);
+    if (poly::is_constant(p)) continue;
+    if (assertMain)
+    {
+      Assert(main_variable(poly) == main_variable(p));
+    }
+    std::vector<poly::Polynomial>::emplace_back(p);
   }
 }
 
-void addPolynomials(std::vector<Polynomial>& polys,
-                    const std::vector<Polynomial>& p)
+void PolyVector::reduce()
 {
-  for (const auto& q : p) addPolynomial(polys, q);
+  std::sort(begin(), end());
+  erase(std::unique(begin(), end()), end());
 }
 
-void makeFinestSquareFreeBasis(std::vector<Polynomial>& polys)
+void PolyVector::makeFinestSquareFreeBasis()
 {
-  for (std::size_t i = 0, n = polys.size(); i < n; ++i)
+  for (std::size_t i = 0, n = size(); i < n; ++i)
   {
     for (std::size_t j = i + 1; j < n; ++j)
     {
-      Polynomial g = gcd(polys[i], polys[j]);
+      Polynomial g = gcd((*this)[i], (*this)[j]);
       if (!is_constant(g))
       {
-        polys[i] = div(polys[i], g);
-        polys[j] = div(polys[j], g);
-        polys.emplace_back(g);
+        (*this)[i] = div((*this)[i], g);
+        (*this)[j] = div((*this)[j], g);
+        add(g);
       }
     }
   }
-  auto it = std::remove_if(polys.begin(), polys.end(), [](const Polynomial& p) {
-    return is_constant(p);
-  });
-  polys.erase(it, polys.end());
-  reduceProjectionPolynomials(polys);
+  auto it = std::remove_if(
+      begin(), end(), [](const Polynomial& p) { return is_constant(p); });
+  erase(it, end());
+  reduce();
+}
+void PolyVector::pushDownPolys(PolyVector& down, poly::Variable var)
+{
+  auto it =
+      std::remove_if(begin(), end(), [&down, &var](const poly::Polynomial& p) {
+        if (main_variable(p) == var) return false;
+        down.add(p);
+        return true;
+      });
+  erase(it, end());
 }
 
-std::vector<Polynomial> projection_mccallum(
-    const std::vector<Polynomial>& polys)
+PolyVector projection_mccallum(const std::vector<Polynomial>& polys)
 {
-  std::vector<Polynomial> res;
+  PolyVector res;
 
   for (const auto& p : polys)
   {
     for (const auto& coeff : coefficients(p))
     {
-      addPolynomial(res, coeff);
+      res.add(coeff);
     }
-    addPolynomial(res, discriminant(p));
+    res.add(discriminant(p));
   }
   for (std::size_t i = 0, n = polys.size(); i < n; ++i)
   {
     for (std::size_t j = i + 1; j < n; ++j)
     {
-      addPolynomial(res, resultant(polys[i], polys[j]));
+      res.add(resultant(polys[i], polys[j]));
     }
   }
 
-  reduceProjectionPolynomials(res);
+  res.reduce();
   return res;
 }
 
diff --git a/src/theory/arith/nl/cad/projections.h b/src/theory/arith/nl/cad/projections.h
index 71c2d5e7f..dd485b372 100644
--- a/src/theory/arith/nl/cad/projections.h
+++ b/src/theory/arith/nl/cad/projections.h
@@ -35,28 +35,46 @@ namespace arith {
 namespace nl {
 namespace cad {
 
-/** Sort and remove duplicates from the list of polynomials. */
-void reduceProjectionPolynomials(std::vector<poly::Polynomial>& polys);
-
 /**
- * Adds a polynomial to the list of projection polynomials.
- * Before adding, it factorizes the polynomials and removed constant factors.
+ * A simple wrapper around std::vector<poly::Polynomial> that ensures that all
+ * polynomials are properly factorized and pruned when added to the list.
  */
-void addPolynomial(std::vector<poly::Polynomial>& polys,
-                   const poly::Polynomial& poly);
-
-/** Adds a list of polynomials using add_polynomial(). */
-void addPolynomials(std::vector<poly::Polynomial>& polys,
-                    const std::vector<poly::Polynomial>& p);
+class PolyVector : public std::vector<poly::Polynomial>
+{
+ private:
+  /** Disable all emplace() */
+  void emplace() {}
+  /** Disable all emplace_back() */
+  void emplace_back() {}
+  /** Disable all insert() */
+  void insert() {}
+  /** Disable all push_back() */
+  void push_back() {}
 
-/** Make a set of polynomials a finest square-free basis. */
-void makeFinestSquareFreeBasis(std::vector<poly::Polynomial>& polys);
+ public:
+  PolyVector() {}
+  /** Construct from a set of polynomials */
+  PolyVector(std::initializer_list<poly::Polynomial> i)
+  {
+    for (const auto& p : i) add(p);
+  }
+  /**
+   * Adds a polynomial to the list of projection polynomials.
+   * Before adding, it factorizes the polynomials and removed constant factors.
+   */
+  void add(const poly::Polynomial& poly, bool assertMain = false);
+  /** Sort and remove duplicates from the list of polynomials. */
+  void reduce();
+  /** Make this list of polynomials a finest square-free basis. */
+  void makeFinestSquareFreeBasis();
+  /** Push polynomials with a lower main variable to another PolyVector. */
+  void pushDownPolys(PolyVector& down, poly::Variable var);
+};
 
 /**
  * Computes McCallum's projection operator.
  */
-std::vector<poly::Polynomial> projectionMcCallum(
-    const std::vector<poly::Polynomial>& polys);
+PolyVector projectionMcCallum(const std::vector<poly::Polynomial>& polys);
 
 }  // namespace cad
 }  // namespace nl