From 497a685f14ff12eb05e4aa6ad7b05682609bf7a9 Mon Sep 17 00:00:00 2001 From: Gereon Kremer Date: Thu, 7 Jan 2021 22:55:31 +0100 Subject: [PATCH] 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. --- src/theory/arith/nl/cad/cdcac.cpp | 75 +++++++------------ src/theory/arith/nl/cad/cdcac.h | 13 ++-- src/theory/arith/nl/cad/cdcac_utils.cpp | 22 +++--- src/theory/arith/nl/cad/cdcac_utils.h | 9 ++- src/theory/arith/nl/cad/projections.cpp | 71 ++++++++++-------- src/theory/arith/nl/cad/projections.h | 48 ++++++++---- test/regress/CMakeLists.txt | 2 + .../regress0/nl/issue5726-downpolys.smt2 | 9 +++ .../regress0/nl/issue5726-sqfactor.smt2 | 8 ++ 9 files changed, 141 insertions(+), 116 deletions(-) create mode 100644 test/regress/regress0/nl/issue5726-downpolys.smt2 create mode 100644 test/regress/regress0/nl/issue5726-sqfactor.smt2 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 -void removeDuplicates(std::vector& 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& ordering) @@ -125,10 +115,11 @@ std::vector CDCAC::getUnsatIntervals( for (const auto& i : intervals) { Trace("cdcac") << "-> " << i << std::endl; - std::vector 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& infeasible, return sampleOutside(infeasible, sample); } -std::vector CDCAC::requiredCoefficients( - const poly::Polynomial& p) const +PolyVector CDCAC::requiredCoefficients(const poly::Polynomial& p) const { - std::vector 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 CDCAC::requiredCoefficients( return res; } -std::vector CDCAC::constructCharacterization( - std::vector& intervals) +PolyVector CDCAC::constructCharacterization(std::vector& intervals) { Assert(!intervals.empty()) << "A covering can not be empty"; Trace("cdcac") << "Constructing characterization now" << std::endl; - std::vector res; - + PolyVector res; for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i) { @@ -198,20 +186,20 @@ std::vector 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 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 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 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& characterization, + const PolyVector& characterization, std::size_t cur_variable, const poly::Value& sample) { - std::vector l; - std::vector u; - std::vector m; - std::vector 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 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 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 constructCharacterization( - std::vector& intervals); + PolyVector constructCharacterization(std::vector& intervals); /** * Constructs an infeasible interval from a characterization. * Implements Algorithm 5. */ - CACInterval intervalFromCharacterization( - const std::vector& 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& polys, +void replace_polynomial(PolyVector& polys, std::size_t id, std::pair subst, CACInterval& interval) @@ -301,7 +301,7 @@ void replace_polynomial(std::vector& 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& 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& 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 #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 d_lowerPolys; + PolyVector d_lowerPolys; /** The polynomials characterizing the upper bound. */ - std::vector d_upperPolys; + PolyVector d_upperPolys; /** The characterizing polynomials in the main variable. */ - std::vector d_mainPolys; + PolyVector d_mainPolys; /** The characterizing polynomials in lower variables. */ - std::vector d_downPolys; + PolyVector d_downPolys; /** The constraints used to derive this interval. */ std::vector 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& 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& 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::emplace_back(p); } } -void addPolynomials(std::vector& polys, - const std::vector& 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& 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 projection_mccallum( - const std::vector& polys) +PolyVector projection_mccallum(const std::vector& polys) { - std::vector 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& 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 that ensures that all + * polynomials are properly factorized and pruned when added to the list. */ -void addPolynomial(std::vector& polys, - const poly::Polynomial& poly); - -/** Adds a list of polynomials using add_polynomial(). */ -void addPolynomials(std::vector& polys, - const std::vector& p); +class PolyVector : public std::vector +{ + 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& polys); + public: + PolyVector() {} + /** Construct from a set of polynomials */ + PolyVector(std::initializer_list 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 projectionMcCallum( - const std::vector& polys); +PolyVector projectionMcCallum(const std::vector& polys); } // namespace cad } // namespace nl diff --git a/test/regress/CMakeLists.txt b/test/regress/CMakeLists.txt index 509427d45..72757ef32 100644 --- a/test/regress/CMakeLists.txt +++ b/test/regress/CMakeLists.txt @@ -684,6 +684,8 @@ set(regress_0_tests regress0/nl/issue3991.smt2 regress0/nl/issue4007-rint-uf.smt2 regress0/nl/issue5534-no-assertions.smt2 + regress0/nl/issue5726-downpolys.smt2 + regress0/nl/issue5726-sqfactor.smt2 regress0/nl/magnitude-wrong-1020-m.smt2 regress0/nl/mult-po.smt2 regress0/nl/nia-wrong-tl.smt2 diff --git a/test/regress/regress0/nl/issue5726-downpolys.smt2 b/test/regress/regress0/nl/issue5726-downpolys.smt2 new file mode 100644 index 000000000..75e6d8cc9 --- /dev/null +++ b/test/regress/regress0/nl/issue5726-downpolys.smt2 @@ -0,0 +1,9 @@ +; COMMAND-LINE: --no-nl-ext --nl-cad +; REQUIRES: poly +; EXPECT: unsat +(set-logic QF_NRA) +(declare-fun x () Real) +(declare-fun y () Real) +(declare-fun z () Real) +(assert (and (> x 0.0) (> 1.0 (+ 1.0 z)) (= y (+ y (* z x))))) +(check-sat) diff --git a/test/regress/regress0/nl/issue5726-sqfactor.smt2 b/test/regress/regress0/nl/issue5726-sqfactor.smt2 new file mode 100644 index 000000000..4608746f6 --- /dev/null +++ b/test/regress/regress0/nl/issue5726-sqfactor.smt2 @@ -0,0 +1,8 @@ +; COMMAND-LINE: --no-nl-ext --nl-cad +; REQUIRES: poly +; EXPECT: sat +(set-logic QF_NRA) +(declare-fun x () Real) +(declare-fun y () Real) +(assert (and (> (* y y y y y y) 0) (> (- x (+ x (* x x (+ (* x x) (* y (- 1.0)))))) 0))) +(check-sat) -- 2.30.2