Test forbidden intervals, disequal case

src/math/polysat/forbidden_intervals.cpp

@@ -244,6 +244,7 @@ namespace polysat {
             fi.coeff = -1;
             fi.interval = to_interval(c, false, fi.coeff, lo_val, lo, hi_val, hi);
             add_non_unit_side_conds(fi, b1, e1, b2, e2);
+            SASSERT(!fi.interval.is_currently_empty());
             return true;
         }
         return false;

src/math/polysat/interval.h

@@ -70,6 +70,8 @@ namespace polysat {
             return {interval::full(), rational::zero(), rational::zero()};
         }
         static eval_interval proper(pdd const &lo, rational const &lo_val, pdd const &hi, rational const &hi_val) {
+            SASSERT(0 <= lo_val && lo_val <= lo.manager().max_value());
+            SASSERT(0 <= hi_val && hi_val <= hi.manager().max_value());
             return {interval::proper(lo, hi), lo_val, hi_val};
         }
 

src/math/polysat/solver.h

@@ -77,6 +77,7 @@ namespace polysat {
         friend class constraint_manager;
         friend class scoped_solverv;
         friend class test_polysat;
+        friend class test_fi;
 
         reslimit&                m_lim;
         params_ref               m_params;

src/math/polysat/viable.cpp

@@ -206,6 +206,10 @@ namespace polysat {
             return floor((e->interval.hi_val() - coeff_val - 1) / e->coeff);
         };
 
+// TODO: looks like there's an infinite loop for:
+// [match_linear3]        a1 = 3; b1 = 0; e1 = 0
+// [match_linear3]        a2 = 3; b2 = -2; e2 = -2
+
         // naive widening. TODO: can we accelerate this?
         // condition e->interval.hi_val() < coeff_val is
         // to ensure that widening is performed on the same interval
@@ -255,7 +259,8 @@ namespace polysat {
                     increase_hi(hi);
                 }
                 LOG("forbidden interval " << e->coeff << " * " << e->interval << " [" << lo << ", " << hi << "[");
-                SASSERT(hi <= max_value);
+                SASSERT(hi <= mod_value);
+                if (hi == mod_value) hi = 0;
                 pdd lop = s.var2pdd(v).mk_val(lo);
                 pdd hip = s.var2pdd(v).mk_val(hi);
                 entry* ne = alloc_entry();
@@ -291,10 +296,11 @@ namespace polysat {
             rational const& b1 = e->interval.lo().val();
             rational const& a2 = e->interval.hi_val();
             rational const& b2 = e->interval.hi().val();
-            rational lhs = a1 * val + b1;
-            rational rhs = a2 * val + b2;
             SASSERT(a1 != a2 && a1 != 0 && a2 != 0);
 
+            rational lhs = mod(a1 * val + b1, mod_value);
+            rational rhs = mod(a2 * val + b2, mod_value);
+
             auto delta_l = [&](rational const& val) {
                 rational m1 = ceil((rhs + 1) / a2);
                 int corr = e->src.is_negative() ? 1 : 0;
@@ -304,7 +310,8 @@ namespace polysat {
                 else
                     m3 = ceil(m3);
 
-                return std::min(m1, m3) - 1;
+                // return std::min(m1, m3) - 1;
+                return std::min(val, std::min(m1, m3) - 1);
             };
             auto delta_u = [&](rational const& val) {
                 rational m1 = ceil((mod_value - lhs) / a1);
@@ -326,8 +333,10 @@ namespace polysat {
                 // TODO: increase interval
                 LOG("refine-disequal-lin: " << " [" << lo << ", " << hi << "[");
 
-                SASSERT(0 <= lo);
-                SASSERT(hi <= max_value);
+                SASSERT(0 <= lo && lo <= val);
+                // SASSERT(val <= hi && hi <= max_value);
+                SASSERT(val <= hi && hi <= mod_value);
+                if (hi == mod_value) hi = 0;
                 pdd lop = s.var2pdd(v).mk_val(lo);
                 pdd hip = s.var2pdd(v).mk_val(hi);
                 entry* ne = alloc_entry();

src/math/polysat/viable.h

@@ -29,6 +29,8 @@ namespace polysat {
     class solver;
 
     class viable {
+        friend class test_fi;
+
         solver& s;
         
         struct entry : public dll_base<entry>, public fi_record { 

src/test/polysat.cpp

@@ -2,6 +2,7 @@
 #include "math/polysat/solver.h"
 #include "ast/ast.h"
 #include "parsers/smt2/smt2parser.h"
+#include "util/util.h"
 #include <vector>
 
 namespace {
@@ -793,21 +794,24 @@ public:
     }
 
     // xy < xz and !Omega(x*y) => y < z
-    static void test_ineq_axiom1(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw/2);
-
-        for (unsigned i = 0; i < 6; ++i) {
-            scoped_solver s(__func__);
+    static void test_ineq_axiom1(unsigned bw = 32, std::optional<unsigned> perm = std::nullopt) {
+        if (perm) {
+            scoped_solver s(std::string(__func__) + " perm=" + std::to_string(*perm));
+            auto const bound = rational::power_of_two(bw/2);
             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);
+            permute_args(*perm, 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_unsat();
+        } else {
+            for (unsigned i = 0; i < 6; ++i) {
+                test_ineq_axiom1(bw, i);
+            }
         }
     }
 
@@ -1062,8 +1066,41 @@ public:
             s.check();
             s.expect_unsat();
         }
- 
+    }
 
+    static void test_fi_linear4(unsigned bw = 4) {
+        {
+            scoped_solver s(__func__);
+            auto y = s.var(s.add_var(bw));
+            s.add_ule(3*y + 1, 10*y);
+            s.check();
+            s.expect_sat();
+        }
+        {
+            scoped_solver s(__func__);
+            auto y = s.var(s.add_var(bw));
+            s.add_ult(3*y + 1, 10*y);
+            s.check();
+            s.expect_sat();
+        }
+        {
+            scoped_solver s(__func__);
+            auto y = s.var(s.add_var(bw));
+            s.add_ule(y-y+8, y);
+            s.add_ule(y, y-y+12);
+            s.add_ult(9*y + 3, 7*y + 1);
+            s.check();
+            s.expect_sat();
+        }
+        {
+            scoped_solver s(__func__);
+            auto y = s.var(s.add_var(bw));
+            s.add_ule(y-y+8, y);
+            s.add_ule(y, y-y+12);
+            s.add_ule(9*y + 3, 7*y + 1);
+            s.check();
+            s.expect_sat();
+        }
     }
 
     // Goal: we probably mix up polysat variables and PDD variables at several points; try to uncover such cases
@@ -1082,6 +1119,107 @@ public:
 };  // class test_polysat
 
 
+/// Here we deal with linear constraints of the form
+///
+///     a1*x + b1 <= a2*x + b2   (mod m = 2^bw)
+///
+/// and their negation.
+class test_fi {
+
+    static bool is_violated(rational const& a1, rational const& b1, rational const& a2, rational const& b2, rational const& val, bool negated, rational const& m) {
+        rational const lhs = (a1*val + b1) % m;
+        rational const rhs = (a2*val + b2) % m;
+        if (negated)
+            return lhs <= rhs;
+        else
+            return lhs > rhs;
+    }
+
+    // Returns true if the input is valid and the test did useful work
+    static bool check_one(rational const& a1, rational const& b1, rational const& a2, rational const& b2, rational const& val, bool negated, unsigned bw) {
+        rational const m = rational::power_of_two(bw);
+        if (a1.is_zero() && a2.is_zero())
+            return false;
+        if (!is_violated(a1, b1, a2, b2, val, negated, m))
+            return false;
+
+        scoped_solver s(__func__);
+        auto x = s.var(s.add_var(bw));
+        signed_constraint c = s.ule(a1*x + b1, a2*x + b2);
+        if (negated)
+            c.negate();
+        viable& v = s.m_viable;
+        v.intersect(x.var(), c);
+        // Trigger forbidden interval refinement
+        v.is_viable(x.var(), val);
+        auto* e = v.m_units[x.var()];
+        if (!e) {
+            std::cout << "test_fi: no interval for a1=" << a1 << " b1=" << b1 << " a2=" << a2 << " b2=" << b2 << " val=" << val << " neg=" << negated << std::endl;
+            // VERIFY(false);
+            return false;
+        }
+        VERIFY(e);
+        auto* first = e;
+        SASSERT(e->next() == e);  // the result is expected to be a single interval (although for this check it doesn't really matter if there's more...)
+        do {
+            rational const& lo = e->interval.lo_val();
+            rational const& hi = e->interval.hi_val();
+            for (rational x = lo; x != hi; x = (x + 1) % m) {
+                // LOG("lo=" << lo << " hi=" << hi << " x=" << x);
+                if (!is_violated(a1, b1, a2, b2, val, negated, m)) {
+                    std::cout << "test_fi: unsound for a1=" << a1 << " b1=" << b1 << " a2=" << a2 << " b2=" << b2 << " val=" << val << " neg=" << negated << std::endl;
+                    VERIFY(false);
+                }
+            }
+            e = e->next();
+        }
+        while (e != first);
+        return true;
+    }
+
+public:
+    static void exhaustive(unsigned bw = 3) {
+        rational const m = rational::power_of_two(bw);
+        for (rational a1(1); a1 < m; ++a1) {
+            for (rational a2(1); a2 < m; ++a2) {
+                // TODO: remove this to test other cases
+                if (a1 == a2)
+                    continue;
+                for (rational b1(0); b1 < m; ++b1)
+                    for (rational b2(0); b2 < m; ++b2)
+                        for (rational val(0); val < m; ++val)
+                            for (bool negated : {true, false})
+                                check_one(a1, b1, a2, b2, val, negated, bw);
+            }
+        }
+    }
+
+    static void randomized(unsigned num_rounds = 10'000, unsigned bw = 16) {
+        rational const m = rational::power_of_two(bw);
+        VERIFY(bw <= 32 && "random_gen generates 32-bit numbers");
+        random_gen rng;
+        unsigned round = num_rounds;
+        while (round) {
+            // rational a1 = (rational(rng()) % (m - 1)) + 1;
+            // rational a2 = (rational(rng()) % (m - 1)) + 1;
+            rational a1 = rational(rng()) % m;
+            rational a2 = rational(rng()) % m;
+            if (a1.is_zero() || a2.is_zero() || a1 == a2)
+                continue;
+            rational b1 = rational(rng()) % m;
+            rational b2 = rational(rng()) % m;
+            rational val = rational(rng()) % m;
+            bool useful =
+                check_one(a1, b1, a2, b2, val, true, bw)
+                || check_one(a1, b1, a2, b2, val, false, bw);
+            if (useful)
+                round--;
+        }
+    }
+
+};  // class test_fi
+
+
     // convert assertions into internal solver state
     // support small grammar of formulas.
     pdd to_pdd(ast_manager& m, solver& s, obj_map<expr, pdd*>& expr2pdd, expr* e) {
@@ -1203,6 +1341,20 @@ public:
 void tst_polysat() {
     using namespace polysat;
 
+    test_fi::exhaustive();
+    test_fi::randomized();
+    return;
+
+    test_polysat::test_fi_linear4();
+    return;
+
+    test_polysat::test_ineq_axiom1(32, 2);  // crashes
+    return;
+
+    // looks like a fishy conflict lemma?
+    test_polysat::test_monot_bounds();
+    return;
+
     test_polysat::test_quot_rem_incomplete();
     test_polysat::test_quot_rem_fixed();
     return;