Polysat: check test results, forbidden intervals for coefficient -1 (#5241)

* Use scoped_ptr for condition

* Check solver result in unit tests

* Add test for unusual cjust

* Add solver::get_value

* Broken assertion

* Support forbidden interval for coefficient -1

src/math/polysat/constraint.h

@@ -68,8 +68,13 @@ namespace polysat {
         bool is_negative() const { return m_status == l_false; }
         bool is_undef() const { return m_status == l_undef; }
 
-        /** Precondition: all variables other than v are assigned. */
+        /** Precondition: all variables other than v are assigned.
-        virtual bool forbidden_interval(solver& s, pvar v, eval_interval& i, constraint*& neg_condition) { return false; }
+         *
+         * \param[out] out_interval     The forbidden interval for this constraint
+         * \param[out] out_neg_cond     Negation of the side condition (the side condition is true when the forbidden interval is trivial). May be NULL if the condition is constant.
+         * \returns True iff a forbidden interval exists and the output parameters were set.
+         */
+        virtual bool forbidden_interval(solver& s, pvar v, eval_interval& out_interval, scoped_ptr<constraint>& out_neg_cond) { return false; }
     };
 
     inline std::ostream& operator<<(std::ostream& out, constraint const& c) { return c.display(out); }

src/math/polysat/eq_constraint.cpp

@@ -113,11 +113,17 @@ namespace polysat {
      */
 
     constraint* eq_constraint::eq_resolve(solver& s, pvar v) {
-        SASSERT(is_currently_true(s));
+        LOG("Resolve " << *this << " upon v" << v);
         if (s.m_conflict.size() != 1)
             return nullptr;
         constraint* c = s.m_conflict[0];
         SASSERT(c->is_currently_false(s));
+        // 'c == this' can happen if propagation was from decide() with only one value left
+        // (e.g., if there's an unsatisfiable clause and we try all values).
+        // Resolution would give us '0 == 0' in this case, which is useless.
+        if (c == this)
+            return nullptr;
+        SASSERT(is_currently_true(s));  // TODO: might not always hold (due to similar case as in comment above?)
         if (c->is_eq()) {
             pdd a = c->to_eq().p();
             pdd b = p();
@@ -146,7 +152,7 @@ namespace polysat {
 
 
     /// Compute forbidden interval for equality constraint by considering it as p <=u 0 (or p >u 0 for disequality)
-    bool eq_constraint::forbidden_interval(solver& s, pvar v, eval_interval& i, constraint*& neg_condition)
+    bool eq_constraint::forbidden_interval(solver& s, pvar v, eval_interval& out_interval, scoped_ptr<constraint>& out_neg_cond)
     {
         SASSERT(!is_undef());
 
@@ -201,8 +207,8 @@ namespace polysat {
             swap(lo, hi);
             lo_val.swap(hi_val);
         }
-        i = eval_interval::proper(lo, lo_val, hi, hi_val);
+        out_interval = eval_interval::proper(lo, lo_val, hi, hi_val);
-        neg_condition = nullptr;
+        out_neg_cond = nullptr;
         return true;
     }
 

src/math/polysat/eq_constraint.h

@@ -32,7 +32,7 @@ namespace polysat {
         bool is_currently_false(solver& s) override;
         bool is_currently_true(solver& s) override;
         void narrow(solver& s) override;
-        bool forbidden_interval(solver& s, pvar v, eval_interval& i, constraint*& neg_condition) override;
+        bool forbidden_interval(solver& s, pvar v, eval_interval& out_interval, scoped_ptr<constraint>& out_neg_cond) override;
 
     private:
         constraint* eq_resolve(solver& s, pvar v);

src/math/polysat/forbidden_intervals.cpp

@@ -70,16 +70,14 @@ namespace polysat {
         rational longest_len;
         unsigned longest_i = UINT_MAX;
         for (constraint* c : conflict) {
-            LOG("constraint: " << *c);
+            LOG_H3("Computing forbidden interval for: " << *c);
             eval_interval interval = eval_interval::full();
-            constraint* neg_cond = nullptr;  // TODO: change to scoped_ptr
+            scoped_ptr<constraint> neg_cond;
             if (c->forbidden_interval(s, v, interval, neg_cond)) {
-                LOG("~> interval: " << interval);
+                LOG("interval: " << interval);
                 LOG("neg_cond: " << show_deref(neg_cond));
-                if (interval.is_currently_empty()) {
+                if (interval.is_currently_empty())
-                    dealloc(neg_cond);
                     continue;
-                }
                 if (interval.is_full())
                     has_full = true;
                 else {
@@ -89,7 +87,7 @@ namespace polysat {
                         longest_i = records.size();
                     }
                 }
-                records.push_back({std::move(interval), neg_cond, c});
+                records.push_back({std::move(interval), std::move(neg_cond), c});
                 if (has_full)
                     break;
             }
@@ -119,6 +117,7 @@ namespace polysat {
             return false;
         }
         LOG("seq: " << seq);
+        SASSERT(seq.size() >= 2);  // otherwise has_full should have been true
 
         p_dependency* d = nullptr;
         unsigned lemma_lvl = 0;

src/math/polysat/solver.h

@@ -272,6 +272,11 @@ namespace polysat {
          */
         pdd var(pvar v) { return m_vars[v]; }
 
+        /**
+         * Return value of v in the current model (only meaningful if check_sat() returned l_true).
+         */
+        rational get_value(pvar v) const { SASSERT(!m_justification[v].is_unassigned()); return m_value[v]; }
+
         /**
          * Create polynomial constraints (but do not activate them).
          * Each constraint is tracked by a dependency.

src/math/polysat/ule_constraint.cpp

@@ -118,7 +118,7 @@ namespace polysat {
             return p.is_val() && q.is_val() && p.val() > q.val();
     }
 
-    bool ule_constraint::forbidden_interval(solver& s, pvar v, eval_interval& i, constraint*& neg_condition)
+    bool ule_constraint::forbidden_interval(solver& s, pvar v, eval_interval& out_interval, scoped_ptr<constraint>& out_neg_cond)
     {
         SASSERT(!is_undef());
 
@@ -136,6 +136,8 @@ namespace polysat {
 
         unsigned const sz = s.size(v);
         dd::pdd_manager& m = s.sz2pdd(sz);
+        rational const pow2 = rational::power_of_two(sz);
+        rational const minus_one = pow2 - 1;
 
         pdd p1 = m.zero();
         pdd e1 = m.zero();
@@ -163,10 +165,10 @@ namespace polysat {
         rational a1 = p1.val();
         rational a2 = p2.val();
         // TODO: to express the interval for coefficient 2^i symbolically, we need right-shift/upper-bits-extract in the language.
-        // So currently we can only do it if the coefficient is 1.
+        // So currently we can only do it if the coefficient is 1 or -1.
-        if (!a1.is_zero() && !a1.is_one())
+        if (!a1.is_zero() && !a1.is_one() && a1 != minus_one)
             return false;
-        if (!a2.is_zero() && !a2.is_one())
+        if (!a2.is_zero() && !a2.is_one() && a2 != minus_one)
             return false;
         /*
         unsigned j1 = 0;
@@ -177,6 +179,9 @@ namespace polysat {
             return false;
         */
 
+        rational const y_coeff = a1.is_zero() ? a2 : a1;
+        SASSERT(!y_coeff.is_zero());
+
         // Concrete values of evaluable terms
         auto e1s = e1.subst_val(s.m_search);
         auto e2s = e2.subst_val(s.m_search);
@@ -219,6 +224,7 @@ namespace polysat {
         else {
             SASSERT(!a1.is_zero());
             SASSERT(!a2.is_zero());
+            SASSERT_EQ(a1, a2);
             // e1 + t <= e2 + t, with t = 2^j1*y = 2^j2*y
             // condition for empty/full: e1 == e2
             is_trivial = e1s.val() == e2s.val();
@@ -234,22 +240,44 @@ namespace polysat {
         if (condition_body.is_val()) {
             // Condition is trivial; no need to create a constraint for that.
             SASSERT(is_trivial == condition_body.is_zero());
-            neg_condition = nullptr;
+            out_neg_cond = nullptr;
         }
         else
-            neg_condition = constraint::eq(level(), s.m_next_bvar++, is_trivial ? neg_t : pos_t, condition_body, m_dep);
+            out_neg_cond = constraint::eq(level(), s.m_next_bvar++, is_trivial ? neg_t : pos_t, condition_body, m_dep);
 
         if (is_trivial) {
             if (is_positive())
-                i = eval_interval::empty(m);
+                // TODO: we cannot use empty intervals for interpolation. So we
+                // can remove the empty case (make it represent 'full' instead),
+                // and return 'false' here. Then we do not need the proper/full
+                // tag on intervals.
+                out_interval = eval_interval::empty(m);
             else
-                i = eval_interval::full();
+                out_interval = eval_interval::full();
         } else {
+            if (y_coeff == minus_one) {
+                // Transform according to:  y \in [l;u[  <=>  -y \in [1-u;1-l[
+                //      -y \in [1-u;1-l[
+                //      <=>  -y - (1 - u) < (1 - l) - (1 - u)    { by: y \in [l;u[  <=>  y - l < u - l }
+                //      <=>  u - y - 1 < u - l                   { simplified }
+                //      <=>  (u-l) - (u-y-1) - 1 < u-l           { by: a < b  <=>  b - a - 1 < b }
+                //      <=>  y - l < u - l                       { simplified }
+                //      <=>  y \in [l;u[.
+                lo = 1 - lo;
+                hi = 1 - hi;
+                swap(lo, hi);
+                lo_val = mod(1 - lo_val, pow2);
+                hi_val = mod(1 - hi_val, pow2);
+                lo_val.swap(hi_val);
+            }
+            else
+                SASSERT(y_coeff.is_one());
+
             if (is_negative()) {
                 swap(lo, hi);
                 lo_val.swap(hi_val);
             }
-            i = eval_interval::proper(lo, lo_val, hi, hi_val);
+            out_interval = eval_interval::proper(lo, lo_val, hi, hi_val);
         }
 
         return true;

src/math/polysat/ule_constraint.h

@@ -38,8 +38,7 @@ namespace polysat {
         bool is_currently_false(solver& s) override;
         bool is_currently_true(solver& s) override;
         void narrow(solver& s) override;
-
+        bool forbidden_interval(solver& s, pvar v, eval_interval& out_interval, scoped_ptr<constraint>& out_neg_cond) override;
-        bool forbidden_interval(solver& s, pvar v, eval_interval& i, constraint*& neg_condition) override;
     };
 
 }

src/test/polysat.cpp

@@ -1,5 +1,7 @@
+#include "math/polysat/log.h"
 #include "math/polysat/solver.h"
 #include "ast/ast.h"
+#include <vector>
 
 namespace polysat {
     // test resolve, factoring routines
@@ -10,7 +12,10 @@ namespace polysat {
     };
 
     struct scoped_solver : public solver_scope, public solver {
-        scoped_solver(): solver(lim) {}
+        scoped_solver(std::string name): solver(lim), m_name(name) {}
+
+        std::string m_name;
+        lbool m_last_result = l_undef;
 
         lbool check_rec() {
             lbool result = check_sat();
@@ -38,12 +43,39 @@ namespace polysat {
         }
 
         void check() {
-            std::cout << check_rec() << "\n";
+            m_last_result = check_rec();
+            std::cout << m_name << ": " << m_last_result << "\n";
             statistics st;
             collect_statistics(st);
             std::cout << st << "\n";
             std::cout << *this << "\n";
         }
+
+        void expect_unsat() {
+            if (m_last_result != l_false) {
+                LOG_H1("FAIL: " << m_name << ": expected UNSAT, got " << m_last_result << "!");
+                VERIFY(false);
+            }
+        }
+
+        void expect_sat(std::vector<std::pair<dd::pdd, unsigned>> const& expected_assignment = {}) {
+            if (m_last_result == l_true) {
+                for (auto const& p : expected_assignment) {
+                    auto const& v_pdd = p.first;
+                    auto const expected_value = p.second;
+                    SASSERT(v_pdd.is_monomial() && !v_pdd.is_val());
+                    auto const v = v_pdd.var();
+                    if (get_value(v) != expected_value) {
+                        LOG_H1("FAIL: " << m_name << ": expected assignment v" << v << " := " << expected_value << ", got value " << get_value(v) << "!");
+                        VERIFY(false);
+                    }
+                }
+            }
+            else {
+                LOG_H1("FAIL: " << m_name << ": expected SAT, got " << m_last_result << "!");
+                VERIFY(false);
+            }
+        }
     };
 
 
@@ -53,7 +85,7 @@ namespace polysat {
 
     /// Creates two separate conflicts (from narrowing) before solving loop is started.
     static void test_add_conflicts() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(3));
         auto b = s.var(s.add_var(3));
         s.add_eq(a + 1);
@@ -61,11 +93,12 @@ namespace polysat {
         s.add_eq(b + 1);
         s.add_eq(b + 2);
         s.check();
+        s.expect_unsat();
     }
 
     /// Has constraints which must be inserted into other watchlist to discover UNSAT
     static void test_wlist() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(3));
         auto b = s.var(s.add_var(3));
         auto c = s.var(s.add_var(3));
@@ -76,8 +109,24 @@ namespace polysat {
         s.add_eq(d + c);
         s.add_eq(d);
         s.check();
+        s.expect_unsat();
     }
 
+    /// Has a constraint in cjust[a] where a does not occur.
+    static void test_cjust() {
+        scoped_solver s(__func__);
+        auto a = s.var(s.add_var(3));
+        auto b = s.var(s.add_var(3));
+        auto c = s.var(s.add_var(3));
+        // 1. Decide a = 0.
+        s.add_eq(a*a + b + 7);                  // 2. Propagate b = 1
+        s.add_eq(b*b + c*c*c*(b+7) + c + 5);    // 3. Propagate c = 2
+        s.add_eq(b*b + c*c);                    // 4. Conflict
+        // Resolution fails because second constraint has c*c*c
+        // => cjust[a] += b*b + c*c
+        s.check();
+        s.expect_unsat();
+    }
 
     /**
      * most basic linear equation solving.
@@ -88,50 +137,49 @@ namespace polysat {
      */
 
     static void test_l1() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(2));
         s.add_eq(a + 1);
         s.check();
-        // Expected result: SAT with a = 3
+        s.expect_sat({{a, 3}});
     }
     
     static void test_l2() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(2));
         auto b = s.var(s.add_var(2));
         s.add_eq(2*a + b + 1);
         s.add_eq(2*b + a);
         s.check();
-        // Expected result: SAT with a = 2, b = 3
+        s.expect_sat({{a, 2}, {b, 3}});
     }
 
     static void test_l3() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(2));
         auto b = s.var(s.add_var(2));
         s.add_eq(3*b + a + 2);
         s.check();
-        // Expected result: SAT
+        s.expect_sat();
     }
 
     static void test_l4() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(3));
-        // auto b = s.var(s.add_var(3));
         s.add_eq(4*a + 2);
         s.check();
-        // Expected result: UNSAT
+        s.expect_unsat();
     }
 
-    // Goal: test propagate_eq in case of 2*a*x + 2*b == 0
     static void test_l5() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(3));
         auto b = s.var(s.add_var(3));
+        s.add_diseq(b);
         s.add_eq(a + 2*b + 4);
         s.add_eq(a + 4*b + 4);
         s.check();
-        // Expected result: UNSAT
+        s.expect_sat({{a, 4}, {b, 4}});
     }
 
     
@@ -140,22 +188,24 @@ namespace polysat {
      * is 0 for all values of a.
      */
     static void test_p1() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(2));
         auto p = a*a*(a*a - 1) + 1;
         s.add_eq(p);
         s.check();
+        s.expect_unsat();
     }
 
     /**
-     * has solution a = 3
+     * has solutions a = 2 and a = 3
      */
     static void test_p2() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto a = s.var(s.add_var(2));
         auto p = a*(a-1) + 2;
         s.add_eq(p);
         s.check();
+        s.expect_sat();
     }
 
     /**
@@ -164,7 +214,7 @@ namespace polysat {
      * - this forces x == 5, which means the first constraint is unsatisfiable by parity.
      */
     static void test_p3() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto x = s.var(s.add_var(4));
         auto y = s.var(s.add_var(4));
         auto z = s.var(s.add_var(4));
@@ -172,32 +222,35 @@ namespace polysat {
         s.add_eq(2*y + z + 8);
         s.add_eq(3*x + 4*y*z + 2*z*z + 1);
         s.check();
+        s.expect_unsat();
     }
 
 
     // Unique solution: u = 5
     static void test_ineq_basic1() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto u = s.var(s.add_var(4));
         auto zero = u - u;
         s.add_ule(u, zero + 5);
         s.add_ule(zero + 5, u);
         s.check();
+        s.expect_sat({{u, 5}});
     }
 
     // Unsatisfiable
     static void test_ineq_basic2() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto u = s.var(s.add_var(4));
         auto zero = u - u;
         s.add_ult(u, zero + 5);
         s.add_ule(zero + 5, u);
         s.check();
+        s.expect_unsat();
     }
 
     // Solutions with u = v = w
     static void test_ineq_basic3() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto u = s.var(s.add_var(4));
         auto v = s.var(s.add_var(4));
         auto w = s.var(s.add_var(4));
@@ -205,11 +258,14 @@ namespace polysat {
         s.add_ule(v, w);
         s.add_ule(w, u);
         s.check();
+        s.expect_sat();
+        SASSERT_EQ(s.get_value(u.var()), s.get_value(v.var()));
+        SASSERT_EQ(s.get_value(u.var()), s.get_value(w.var()));
     }
 
     // Unsatisfiable
     static void test_ineq_basic4() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto u = s.var(s.add_var(4));
         auto v = s.var(s.add_var(4));
         auto w = s.var(s.add_var(4));
@@ -217,29 +273,32 @@ namespace polysat {
         s.add_ult(v, w);
         s.add_ule(w, u);
         s.check();
+        s.expect_unsat();
     }
 
     // Satisfiable
     // Without forbidden intervals, we just try values for u until it works
     static void test_ineq_basic5() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto u = s.var(s.add_var(4));
         auto v = s.var(s.add_var(4));
         auto zero = u - u;
         s.add_ule(zero + 12, u + v);
         s.add_ule(v, zero + 2);
         s.check();
+        s.expect_sat();  // e.g., u = 12, v = 0
     }
 
-    // Like 5 but the other forbidden interval will be the longest
+    // Like test_ineq_basic5 but the other forbidden interval will be the longest
     static void test_ineq_basic6() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto u = s.var(s.add_var(4));
         auto v = s.var(s.add_var(4));
         auto zero = u - u;
         s.add_ule(zero + 14, u + v);
         s.add_ule(v, zero + 2);
         s.check();
+        s.expect_sat();
     }
 
 
@@ -250,7 +309,7 @@ namespace polysat {
      * v*q > u
      */
     static void test_ineq1() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto u = s.var(s.add_var(5));
         auto v = s.var(s.add_var(5));
         auto q = s.var(s.add_var(5));
@@ -259,6 +318,7 @@ namespace polysat {
         s.add_ult(r, u);
         s.add_ult(u, v*q);
         s.check();
+        s.expect_unsat();
     }
 
     /**
@@ -268,7 +328,7 @@ namespace polysat {
      * n > r2 > 0
      */
     static void test_ineq2() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto n = s.var(s.add_var(5));
         auto q1 = s.var(s.add_var(5));
         auto a = s.var(s.add_var(5));
@@ -281,6 +341,7 @@ namespace polysat {
         s.add_ult(r2, n);
         s.add_diseq(n);
         s.check();
+        s.expect_unsat();
     }
 
 
@@ -289,7 +350,7 @@ namespace polysat {
      *  expected: unsat
      */
     static void test_monot1() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto bw = 5;
 
         auto tb1 = s.var(s.add_var(bw));
@@ -331,6 +392,7 @@ namespace polysat {
 
         s.add_ult(tb1, tb2);
         s.check();
+        s.expect_unsat();
     }
 
     /**
@@ -340,7 +402,7 @@ namespace polysat {
      *  only difference to (1) is the inequality right before the check
      */
     static void test_monot2() {
-        scoped_solver s;
+        scoped_solver s(__func__);
         auto bw = 5;
 
         auto tb1 = s.var(s.add_var(bw));
@@ -382,13 +444,14 @@ namespace polysat {
 
         s.add_ult(tb1 + err, tb2);
         s.check();
+        s.expect_unsat();
     }
 
 
     // Goal: we probably mix up polysat variables and PDD variables at several points; try to uncover such cases
     // NOTE: actually, add_var seems to keep them in sync, so this is not an issue at the moment (but we should still test it later)
     // static void test_mixed_vars() {
-    //     scoped_solver s;
+    //     scoped_solver s(__func__);
     //     auto a = s.var(s.add_var(2));
     //     auto b = s.var(s.add_var(4));
     //     auto c = s.var(s.add_var(2));
@@ -409,6 +472,7 @@ namespace polysat {
 void tst_polysat() {
     polysat::test_add_conflicts();
     polysat::test_wlist();
+    polysat::test_cjust();
     polysat::test_l1();
     polysat::test_l2();
     polysat::test_l3();