wip

src/math/polysat/conflict.cpp

@@ -72,8 +72,19 @@ namespace polysat {
     void conflict::set(signed_constraint c) {
         LOG("Conflict: " << c);
         SASSERT(empty());
+        if (c.bvalue(s) == l_false) {
+            auto* cl = s.m_bvars.reason(c.blit().var());
+            if (cl)
+                set(*cl);
+            else
+                insert(c);
+        }
+        else {
+            SASSERT(c.is_currently_false(s));
+            SASSERT(c.bvalue(s) == l_true);
             c->set_var_dependent();
             insert(c);
+        }
         SASSERT(!empty());
     }
 
@@ -102,11 +113,8 @@ namespace polysat {
     void conflict::set(clause const& cl) {
         LOG("Conflict: " << cl);
         SASSERT(empty());
-        for (auto lit : cl) {
+        for (auto lit : cl)          
-            auto c = s.lit2cnstr(lit);           
+            insert(~s.lit2cnstr(lit));        
-            // no c->set_var_dependent();
-            insert(~c);
-        }
         SASSERT(!empty());
     }
 
@@ -117,11 +125,13 @@ namespace polysat {
      *   should appear, otherwise the lemma would be a tautology
      */
     void conflict::insert(signed_constraint c) {
+
         if (c.is_always_true())
             return;
         if (c->is_marked())
             return;
         LOG("inserting: " << c);
+        SASSERT(!c->vars().empty());
         set_mark(c);
         if (c->has_bvar())
             insert_literal(c.blit());
@@ -230,8 +240,8 @@ namespace polysat {
                 continue;
             auto diseq = ~s.eq(s.var(v), s.get_value(v));
             cm().ensure_bvar(diseq.get());
-            //if (diseq.bvalue(s) == l_undef) 
+            if (diseq.bvalue(s) == l_undef) 
-            //    s.assign_bool(s.get_level(v), ~diseq.blit(), nullptr, nullptr);            
+                s.assign_bool(s.get_level(v), ~diseq.blit(), nullptr, nullptr);            
             lemma.push(diseq);
         }        
 

src/math/polysat/constraint.cpp

@@ -85,6 +85,9 @@ namespace polysat {
         for (unsigned i = 0; i < cl.size(); ++i) {
             if (m_bvars.is_false(cl[i]))
                 continue;
+            signed_constraint sc = s.lit2cnstr(cl[i]);
+            if (sc.is_currently_false(s)) 
+                continue;            
             m_bvars.watch(cl[i]).push_back(&cl);
             std::swap(cl[!first], cl[i]);
             if (!first)

src/math/polysat/forbidden_intervals.cpp

@@ -28,7 +28,7 @@ namespace polysat {
 
     struct fi_record {
         eval_interval interval;
-        signed_constraint neg_cond;  // could be multiple constraints later
+        vector<signed_constraint> side_cond;  
         signed_constraint src;
     };
 
@@ -75,10 +75,10 @@ namespace polysat {
     * We assume that neg_cond is a consequence of src that  
     * does not mention the variable v to be eliminated.
     */
-    void forbidden_intervals::full_interval_conflict(signed_constraint src, signed_constraint neg_cond, conflict& core) {
+    void forbidden_intervals::full_interval_conflict(signed_constraint src, vector<signed_constraint> const& side_cond, conflict& core) {
-        SASSERT(neg_cond);
         core.reset();
-        core.insert(~neg_cond);
+        for (auto c : side_cond)
+            core.insert(c);
         core.insert(src);
         core.set_bailout();
     }
@@ -92,16 +92,16 @@ namespace polysat {
         for (signed_constraint c : just) {
             LOG_H3("Computing forbidden interval for: " << c);
             eval_interval interval = eval_interval::full();
-            signed_constraint neg_cond;
+            vector<signed_constraint> side_cond;
-            if (get_interval(c, v, interval, neg_cond)) {
+            if (get_interval(c, v, interval, side_cond)) {
                 LOG("interval: " << interval);
-                LOG("neg_cond: " << neg_cond);
+                LOG("neg_cond: " << side_cond);
                 if (interval.is_currently_empty())
                     continue;
                 if (interval.is_full()) {
                     // We have a single interval covering the whole domain
                     // => the side conditions of that interval are enough to produce a conflict
-                    full_interval_conflict(c, neg_cond, core);
+                    full_interval_conflict(c, side_cond, core);
                     revert_core(core);
                     return true;
                 }
@@ -112,7 +112,7 @@ namespace polysat {
                         longest_i = records.size();
                     }
                 }
-                records.push_back({ std::move(interval), std::move(neg_cond), c });
+                records.push_back({ std::move(interval), std::move(side_cond), c });
             }
         }
 
@@ -164,9 +164,8 @@ namespace polysat {
             core.insert(c);
             // Side conditions
             // TODO: check whether the condition is subsumed by c?  maybe at the end do a "lemma reduction" step, to try and reduce branching?
-            signed_constraint& neg_cond = records[i].neg_cond;
+            for (auto sc : records[i].side_cond)
-            if (neg_cond)
+                core.insert(sc);
-                core.insert(~neg_cond);
             
             core.insert(records[i].src);
         }
@@ -183,7 +182,7 @@ namespace polysat {
      * \returns True iff a forbidden interval exists and the output parameters were set.
      */
 
-    bool forbidden_intervals::get_interval(signed_constraint const& c, pvar v, eval_interval& out_interval, signed_constraint& out_neg_cond)
+    bool forbidden_intervals::get_interval(signed_constraint const& c, pvar v, eval_interval& out_interval, vector<signed_constraint>& out_side_cond)
     {
         if (!c->is_ule())
             return false;
@@ -232,10 +231,19 @@ namespace polysat {
         // test_monot_bounds(8)
         // 
         // Currently only works if coefficient is a power of two
-        if (!p1.is_val())
+        if (!p1.is_val()) {
+            pdd q1 = p1.subst_val(s.assignment());
+            if (!q1.is_val())
                 return false;
-        if (!p2.is_val())
+            out_side_cond.push_back(s.eq(q1, p1));
+        }
+        if (!p2.is_val()) {
+            pdd q2 = p2.subst_val(s.assignment());
+            if (!q2.is_val())
                 return false;
+            out_side_cond.push_back(s.eq(q2, p2));
+            p2 = q2;
+        }
         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.
@@ -320,12 +328,11 @@ 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());
-            out_neg_cond = nullptr;
         }
         else if (is_trivial)
-            out_neg_cond = ~s.m_constraints.eq(condition_body);
+            out_side_cond.push_back(s.m_constraints.eq(condition_body));
         else
-            out_neg_cond = s.m_constraints.eq(condition_body);
+            out_side_cond.push_back(~s.m_constraints.eq(condition_body));
 
         if (is_trivial) {
             if (c.is_positive())

src/math/polysat/forbidden_intervals.h

@@ -22,8 +22,8 @@ namespace polysat {
     class forbidden_intervals {
         solver& s;
         void revert_core(conflict& core);
-        void full_interval_conflict(signed_constraint c, signed_constraint neg_cond, conflict& core);
+        void full_interval_conflict(signed_constraint c, vector<signed_constraint> const & side_cond, conflict& core);
-        bool get_interval(signed_constraint const& c, pvar v, eval_interval& out_interval, signed_constraint& out_neg_cond);
+        bool get_interval(signed_constraint const& c, pvar v, eval_interval& out_interval, vector<signed_constraint>& out_side_cond);
     public:
         forbidden_intervals(solver& s) :s(s) {}
         bool perform(pvar v, vector<signed_constraint> const& just, conflict& core);

src/math/polysat/saturation.cpp

@@ -115,7 +115,7 @@ namespace polysat {
         if ((x_lo + 1) * y_lo <= bound) {
             x_hi = x_max;
             while (x_lo < x_hi) {
-                rational x_mid = div(x_hi + x_lo, two);
+                rational x_mid = div(x_hi + x_lo + 1, two);
                 if (x_mid * y_lo > bound)
                     x_hi = x_mid - 1;
                 else
@@ -125,7 +125,7 @@ namespace polysat {
         else if (x_lo * (y_lo + 1) <= bound) {
             y_hi = y_max;
             while (y_lo < y_hi) {
-                rational y_mid = div(y_hi + y_lo, two);
+                rational y_mid = div(y_hi + y_lo + 1, two);
                 if (y_mid * x_lo > bound)
                     y_hi = y_mid - 1;
                 else
@@ -144,6 +144,7 @@ namespace polysat {
         auto c2 = s.ule(y, pddm.mk_val(y_lo));
         new_constraints.insert(c1);
         new_constraints.insert(c2);
+        LOG("bounded " << bound << " : " << c1 << " " << c2);
     }
 
     // determine worst case upper bounds for x, y
@@ -390,8 +391,14 @@ namespace polysat {
 
     // [x] p(x) <= q(x) where value(p) > value(q)
     //     ==> q <= value(q) => p <= value(q)
+    // 
+    // for strict?
+    //     p(x) < q(x) where value(p) >= value(q)
+    //     ==> value(p) <= p => value(p) < q
 
     bool inf_saturate::try_tangent(pvar v, conflict& core, inequality const& c) {        
+        if (c.is_strict)
+            return false;
         if (!c.src->contains_var(v))
             return false;
         if (c.lhs.is_val() || c.rhs.is_val())
@@ -407,8 +414,14 @@ namespace polysat {
         SASSERT(!c.is_strict || l_val >= r_val);
         vector<signed_constraint> new_constraints;
         new_constraints.push_back(c.as_signed_constraint());
+        if (c.is_strict) {
+            new_constraints.push_back(s.ule(l_val, c.lhs));
+            return propagate(core, c, c, s.ult(r_val, c.rhs), new_constraints);
+        }
+        else {
             new_constraints.push_back(s.ule(c.rhs, r_val));
-        return propagate(core, c, c, c.is_strict ? s.ult(c.lhs, r_val) : s.ule(c.lhs, r_val), new_constraints);
+            return propagate(core, c, c, s.ule(c.lhs, r_val), new_constraints);
+        }
     }
 
 }

src/math/polysat/solver.cpp

@@ -122,13 +122,16 @@ namespace polysat {
         m_constraints.ensure_bvar(c.get());
         sat::literal lit = c.blit();
         LOG("New constraint: " << c);
-        if (m_bvars.is_false(lit) || c.is_currently_false(*this)) {
+        if (m_bvars.is_false(lit)) {
             set_conflict(c /*, dep*/);
             return;
         }
         m_bvars.assign(lit, m_level, nullptr, nullptr, dep);
         m_trail.push_back(trail_instr_t::assign_bool_i);
         m_search.push_boolean(lit);
+        if (c.is_currently_false(*this)) 
+            set_conflict(c /*, dep */);
+        
 
 #if ENABLE_LINEAR_SOLVER
         m_linear_solver.new_constraint(*c.get());
@@ -197,13 +200,13 @@ namespace polysat {
         if (m_bvars.is_true(cl[idx]))
             return false;
         unsigned i = 2;
-        for (; i < cl.size() && m_bvars.is_false(cl[i]); ++i);
+        for (; i < cl.size() && (m_bvars.is_false(cl[i]) || lit2cnstr(cl[i]).is_currently_false(*this)); ++i);
         if (i < cl.size()) {
             m_bvars.watch(cl[i]).push_back(&cl);
             std::swap(cl[1 - idx], cl[i]);
             return true;
         }
-        if (m_bvars.is_false(cl[idx]))
+        if (m_bvars.is_false(cl[idx]) || lit2cnstr(cl[idx]).is_currently_false(*this))
             set_conflict(cl);
         else
             assign_bool(level(cl), cl[idx], &cl, nullptr);

src/math/polysat/solver.h

@@ -268,12 +268,15 @@ namespace polysat {
         signed_constraint eq(pdd const& p, rational const& q) { return eq(p - q); }
         signed_constraint diseq(pdd const& p, rational const& q) { return diseq(p - q); }
         signed_constraint ule(pdd const& p, pdd const& q) { return m_constraints.ule(p, q); }
-        signed_constraint ule(pdd const& p, rational const& q) { return m_constraints.ule(p, p.manager().mk_val(q)); }
+        signed_constraint ule(pdd const& p, rational const& q) { return ule(p, p.manager().mk_val(q)); }
+        signed_constraint ule(rational const& p, pdd const& q) { return ule(q.manager().mk_val(p), q); }
         signed_constraint ult(pdd const& p, pdd const& q) { return m_constraints.ult(p, q); }
-        signed_constraint ult(pdd const& p, rational const& q) { return m_constraints.ult(p, p.manager().mk_val(q)); }
+        signed_constraint ult(pdd const& p, rational const& q) { return ult(p, p.manager().mk_val(q)); }
+        signed_constraint ult(rational const& p, pdd const& q) { return ult(q.manager().mk_val(p), q); }
         signed_constraint sle(pdd const& p, pdd const& q) { return m_constraints.sle(p, q); }
         signed_constraint slt(pdd const& p, pdd const& q) { return m_constraints.slt(p, q); }
 
+
         /** Create and activate polynomial constraints. */
         void add_eq(pdd const& p, unsigned dep = null_dependency)                { assign_eh(eq(p), dep); }
         void add_diseq(pdd const& p, unsigned dep = null_dependency)             { assign_eh(diseq(p), dep); }

src/test/polysat.cpp

@@ -781,32 +781,7 @@ namespace polysat {
 
     // xy < xz and !Omega(x*y) => y < z
     static void test_ineq_axiom1(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw-1);
+        auto const bound = rational::power_of_two(bw/2);
-
-        {
-            scoped_solver s(__func__);
-            auto x = s.var(s.add_var(bw));
-            auto y = s.var(s.add_var(bw));
-            auto z = s.var(s.add_var(bw));
-            s.add_ult(x * y, x * z);
-            s.add_ule(z, y);
-            s.add_ult(x, bound);
-            s.add_ult(y, bound);
-            s.check();
-            s.expect_unsat();
-        }
-        {
-            scoped_solver s(__func__);
-            auto y = s.var(s.add_var(bw));
-            auto x = s.var(s.add_var(bw));
-            auto z = s.var(s.add_var(bw));
-            s.add_ult(x * y, x * z);
-            s.add_ule(z, y);
-            s.add_ult(x, bound);
-            s.add_ult(y, bound);
-            s.check();
-            s.expect_unsat();
-        }
 
         for (unsigned i = 0; i < 6; ++i) {
             scoped_solver s(__func__);
@@ -823,9 +798,27 @@ namespace polysat {
         }
     }
 
+    static void test_ineq_non_axiom1(unsigned bw = 32) {
+        auto const bound = rational::power_of_two(bw - 1);
+
+        for (unsigned i = 0; i < 6; ++i) {
+            scoped_solver s(__func__);
+            auto x = s.var(s.add_var(bw));
+            auto y = s.var(s.add_var(bw));
+            auto z = s.var(s.add_var(bw));
+            permute_args(i, x, y, z);
+            s.add_ult(x * y, x * z);
+            s.add_ule(z, y);
+            //s.add_ult(x, bound);
+            //s.add_ult(y, bound);
+            s.check();
+            s.expect_sat();
+        }
+    }
+
     // xy <= xz & !Omega(x*y) => y <= z or x = 0
     static void test_ineq_axiom2(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw - 1);
+        auto const bound = rational::power_of_two(bw/2);
         for (unsigned i = 0; i < 6; ++i) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
@@ -844,7 +837,7 @@ namespace polysat {
 
     // xy < b & a <= x  & !Omega(x*y) => a*y < b
     static void test_ineq_axiom3(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw - 1);
+        auto const bound = rational::power_of_two(bw/2);
         for (unsigned i = 0; i < 24; ++i) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
@@ -862,9 +855,9 @@ namespace polysat {
         }
     }
 
-    // xy <= b & a <= x & !Omega(x*y) => a*y <= b
+    // x*y <= b & a <= x & !Omega(x*y) => a*y <= b
     static void test_ineq_axiom4(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw - 1);
+        auto const bound = rational::power_of_two(bw/2);
         for (unsigned i = 0; i < 24; ++i) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
@@ -882,9 +875,29 @@ namespace polysat {
         }
     }
 
+    // x*y <= b & a <= x & !Omega(x*y) => a*y <= b
+    static void test_ineq_non_axiom4(unsigned bw = 32) {
+        auto const bound = rational::power_of_two(bw - 1);
+        for (unsigned i = 0; i < 24; ++i) {
+            scoped_solver s(__func__);
+            auto x = s.var(s.add_var(bw));
+            auto y = s.var(s.add_var(bw));
+            auto a = s.var(s.add_var(bw));
+            auto b = s.var(s.add_var(bw));
+            permute_args(i, x, y, a, b);
+            s.add_ule(x * y, b);
+            s.add_ule(a, x);
+            s.add_ult(x, bound);
+            s.add_ult(y, bound);
+            s.add_ult(b, a * y);
+            s.check();
+            s.expect_sat();
+        }
+    }
+
     // a < xy & x <= b  & !Omega(x*y) => a < b*y
     static void test_ineq_axiom5(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw - 1);
+        auto const bound = rational::power_of_two(bw/2);
         for (unsigned i = 0; i < 24; ++i) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
@@ -904,7 +917,7 @@ namespace polysat {
 
     // a <= xy & x <= b & !Omega(x*y) => a <= b*y
     static void test_ineq_axiom6(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw - 1);
+        auto const bound = rational::power_of_two(bw/2);
         for (unsigned i = 0; i < 24; ++i) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
@@ -1047,15 +1060,14 @@ namespace polysat {
 
 void tst_polysat() {
 
-    polysat::test_cjust();
+//    polysat::test_ineq_axiom1();
-
+//    polysat::test_ineq_axiom2();
-    
+//    polysat::test_ineq_axiom3();
-
+    polysat::test_ineq_non_axiom1();
-
+    polysat::test_ineq_non_axiom4();
-
+    polysat::test_ineq_axiom4();
-    
+    polysat::test_ineq_axiom5();
-
+    polysat::test_ineq_axiom6();
-
 
     // not working
     // polysat::test_fixed_point_arith_div_mul_inverse();