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.
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)
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}});
}
}
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;
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)
{
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)
{
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)
{
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));
}
}
}
{
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;
{
if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
{
- l.emplace_back(p);
+ l.add(p, true);
}
}
d_assignment.unset(d_variableOrdering[cur_variable]);
{
if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
{
- u.emplace_back(p);
+ u.add(p, true);
}
}
d_assignment.unset(d_variableOrdering[cur_variable]);
* 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.
* 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)
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))
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
{
else
{
// Push to d_downPolys
- interval.d_downPolys.emplace_back(subst.second);
+ interval.d_downPolys.add(subst.second);
}
}
Assert(replaced)
}
}
}
- 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
#include <vector>
#include "expr/node.h"
+#include "theory/arith/nl/cad/projections.h"
namespace CVC4 {
namespace theory {
/** 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;
};
#ifdef CVC4_POLY_IMP
+#include "base/check.h"
+
namespace CVC4 {
namespace theory {
namespace arith {
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;
}
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
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
--- /dev/null
+; 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)
--- /dev/null
+; 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)
projects / cvc5.git / commitdiff