First try to generalize variable elimination

src/math/polysat/constraint_manager.cpp

@@ -277,10 +277,12 @@ namespace polysat {
         auto& m = a.manager();
         unsigned sz = m.power_of_2();
         if (a.is_val() && b.is_val()) {
-            // TODO: just evaluate?
+            rational r;
+            rational q = machine_div_rem(a.val(), b.val(), r);
+            return { m.mk_val(r), m.mk_val(q) };
         }
 
-        constraint_dedup::quot_rem_args args({a, b});
+        constraint_dedup::quot_rem_args args({ a, b });
         auto it = m_dedup.quot_rem_expr.find_iterator(args);
         if (it != m_dedup.quot_rem_expr.end())
             return { m.mk_var(it->m_value.first), m.mk_var(it->m_value.second) };
@@ -297,13 +299,13 @@ namespace polysat {
         //      b = 0  ==>  q = -1
         s.add_eq(a, b * q + r);
         s.add_umul_noovfl(b, q);
-        s.add_ule(r, b*q+r);
+        s.add_ule(r, b * q + r);
 
         auto c_eq = eq(b);
         s.add_clause(c_eq, ult(r, b), false);
         s.add_clause(~c_eq, eq(q + 1), false);
 
-        return {q, r};
+        return { q, r };
     }
 
     pdd constraint_manager::bnot(pdd const& p) {

src/math/polysat/constraint_manager.h

@@ -81,7 +81,7 @@ namespace polysat {
 
         void ensure_bvar(constraint* c);
         void erase_bvar(constraint* c);
-
+        
         signed_constraint mk_op_constraint(op_constraint::code op, pdd const& p, pdd const& q, pdd const& r);
         pdd mk_op_term(op_constraint::code op, pdd const& p, pdd const& q);
 

src/math/polysat/log.cpp

@@ -32,6 +32,10 @@ void set_log_enabled(bool log_enabled) {
   g_log_enabled = log_enabled;
 }
 
+bool get_log_enabled() {
+  return g_log_enabled;
+}
+
 static LogLevel get_max_log_level(std::string const& fn, std::string const& pretty_fn) {
   (void)fn;
   (void)pretty_fn;

src/math/polysat/log.h

@@ -27,6 +27,7 @@ char const* color_reset();
 #if POLYSAT_LOGGING_ENABLED
 
 void set_log_enabled(bool log_enabled);
+bool get_log_enabled();
 
 class polysat_log_indent
 {

src/math/polysat/op_constraint.h

@@ -48,7 +48,7 @@ namespace polysat {
 
         void narrow_shl(solver& s);
         static lbool eval_shl(pdd const& p, pdd const& q, pdd const& r);
-
+        
         void narrow_and(solver& s);
         static lbool eval_and(pdd const& p, pdd const& q, pdd const& r);
 

src/math/polysat/solver.cpp

@@ -1042,6 +1042,11 @@ namespace polysat {
         }
     }
 
+    void solver::add_clause(signed_constraint c, bool is_redundant) {
+        signed_constraint cs[1] = { c };
+        add_clause(1, cs, is_redundant);
+    }
+
     void solver::add_clause(signed_constraint c1, signed_constraint c2, bool is_redundant) {
         signed_constraint cs[2] = { c1, c2 };
         add_clause(2, cs, is_redundant);

src/math/polysat/solver.h

@@ -224,6 +224,7 @@ namespace polysat {
         void learn_lemma(clause& lemma);
         void backjump(unsigned new_level);
         void add_clause(clause& clause);
+        void add_clause(signed_constraint c, bool is_redundant);
         void add_clause(signed_constraint c1, signed_constraint c2, bool is_redundant);
         void add_clause(signed_constraint c1, signed_constraint c2, signed_constraint c3, bool is_redundant);
         void add_clause(signed_constraint c1, signed_constraint c2, signed_constraint c3, signed_constraint c4, bool is_redundant);
@@ -309,7 +310,7 @@ namespace polysat {
 
         /** Create expression for the logical left shift of p by q. */
         pdd shl(pdd const& p, pdd const& q) { return m_constraints.shl(p, q); }
-
+        
         /** Create expression for the bit-wise negation of p. */
         pdd bnot(pdd const& p) { return m_constraints.bnot(p); }
 

src/math/polysat/variable_elimination.cpp

@@ -18,6 +18,114 @@ Author:
 #include <algorithm>
 
 namespace polysat {
+    
+    pdd free_variable_elimination::get_odd(pdd p) {
+        SASSERT(p.is_val() || p.is_var()); // For now
+        
+        if (p.is_val()) {
+            const rational& v = p.val();
+            unsigned d = v.trailing_zeros(); 
+            if (!d)
+                return s.sz2pdd(p.power_of_2()).mk_val(v);
+            return s.sz2pdd(p.power_of_2()).mk_val(div(v, rational::power_of_two(d))); // TODO: Is there no shift?
+        }
+        pvar v = p.var();
+        if (m_rest_constants.size() > v && m_rest_constants[v] != -1)
+            return s.var(m_rest_constants[v]);
+        
+        get_dyadic_valuation(p);
+        
+        pvar rest = s.add_var(p.power_of_2());
+        m_rest_constants.setx(v, rest, -1);
+        s.add_clause(s.eq(s.var(m_pv_power_constants[v]) * s.var(rest), p), false);
+        return s.var(rest);
+    }
+    
+    optional<pdd> free_variable_elimination::get_inverse(pdd p) {
+        SASSERT(p.is_val() || p.is_var()); // For now
+        
+        if (p.is_val()) {
+            pdd i = s.sz2pdd(p.power_of_2()).zero();
+            if (!inv(p, i))
+                return {};
+            return optional<pdd>(i);
+        }
+        pvar v = p.var();
+        if (m_inverse_constants.size() > v && m_inverse_constants[v] != -1)
+            return optional<pdd>(s.var(m_inverse_constants[v]));
+        
+        pvar inv = s.add_var(s.var(v).power_of_2());
+        m_inverse_constants.setx(v, inv, -1);
+        s.add_clause(s.eq(s.var(inv) * s.var(v), s.sz2pdd(s.var(v).power_of_2()).one()), false);
+        return optional<pdd>(s.var(inv));
+    }
+    
+    // symbolic version of "max_pow2_divisor" for checking if it is exactly "k"
+    void free_variable_elimination::add_dyadic_valuation(pvar v, unsigned k) {
+        pvar pv;
+        pvar pv2;
+        if (m_pv_constants.size() <= v || m_pv_constants[v] == -1) {
+            pv = s.add_var(s.var(v).power_of_2()); // TODO: What's a good value? Unfortunately we cannot use a integer
+            pv2 = s.add_var(s.var(v).power_of_2());
+            m_pv_constants.setx(v, pv, -1);
+            m_pv_power_constants.setx(v, pv2, -1);
+        }
+        else {
+            pv = m_pv_constants[v];
+            pv2 = m_pv_power_constants[v];
+        }
+        
+        // (pv = k && pv2 = 2^k) <==> ((v & (2^(k + 1) - 1)) = 2^k) 
+        
+        rational mask = rational::power_of_two(k + 1) - 1;
+        pdd masked = s.band(s.var(v), s.sz2pdd(s.var(v).power_of_2()).mk_val(mask));
+        std::pair<pdd, pdd> odd_part = s.quot_rem(s.var(v), s.var(pv2));
+        
+        signed_constraint c1 = s.eq(s.var(pv), k);
+        signed_constraint c2 = s.eq(s.var(pv2), rational::power_of_two(k));
+        signed_constraint c3 = s.eq(masked, rational::power_of_two(k));
+        bool e = get_log_enabled();
+        set_log_enabled(false);
+        s.add_clause(c1, ~c3, false);
+        s.add_clause(c2, ~c3, false);
+        s.add_clause(~c1, ~c2, c3, false);
+        
+        s.add_clause(s.eq(odd_part.second, 0), false); // The division has to be exact
+        set_log_enabled(e);
+    }
+    
+    std::pair<pdd, pdd> free_variable_elimination::get_dyadic_valuation(pdd p, unsigned short lower, unsigned short upper) {
+        SASSERT(p.is_val() || p.is_var()); // For now
+        
+        if (p.is_val()) {
+            rational pv(p.val().trailing_zeros());
+            rational pv2 = rational::power_of_two(p.val().trailing_zeros());
+            return { s.sz2pdd(p.power_of_2()).mk_val(pv), s.sz2pdd(p.power_of_2()).mk_val(pv2) };
+        }
+        
+        pvar v = p.var();
+        unsigned short prev_lower = 0, prev_upper = 0;
+        if (m_has_validation_of_range.size() > v) {
+            unsigned range = m_has_validation_of_range[v];
+            prev_lower = range & 0xFFFF;
+            prev_upper = range >> 16;
+            if (lower >= prev_lower && upper <= prev_upper)
+                return { s.var(m_pv_constants[v]), s.var(m_pv_power_constants[v]) }; // exists already in the required range
+        }
+        LOG("Adding valuation function for variable " << v  << " in [" << lower << "; " << upper << ")");
+        m_has_validation_of_range.setx(v, lower | (unsigned)upper << 16, 0);
+        for (unsigned i = lower; i < prev_lower; i++) {
+            add_dyadic_valuation(v, i);
+        }
+        for (unsigned i = prev_upper; i < upper; i++) {
+            add_dyadic_valuation(v, i);
+        }
+        return { s.var(m_pv_constants[v]), s.var(m_pv_power_constants[v]) };
+    }
+    
+    std::pair<pdd, pdd> free_variable_elimination::get_dyadic_valuation(pdd p) {
+        return get_dyadic_valuation(p, 0, p.power_of_2()); 
+    }
 
     void free_variable_elimination::find_lemma(conflict& core) {
         LOG_H1("Free Variable Elimination");
@@ -42,6 +150,27 @@ namespace polysat {
     }
 
     void free_variable_elimination::find_lemma(pvar v, signed_constraint c, conflict& core) {
+        vector<signed_constraint> to_check;
+        // Find another constraint where we want to substitute v
+        for (signed_constraint c_target : core) {
+            if (c == c_target)
+                continue;
+            if (c_target.vars().size() <= 1)
+                continue;
+            if (!c_target.contains_var(v))
+                continue;
+            // TODO: helper method constraint::subst(pvar v, pdd const& p)
+            //       (or rather, add it on constraint_manager since we need to allocate/dedup the new constraint)
+            //  For now, just restrict to ule_constraint.
+            if (!c_target->is_ule()) // TODO: Remove?
+                continue;
+            // TODO: Eliminate without inversion? 2x = y && 2x <= z
+            
+            to_check.push_back(c_target);
+        }
+        if (to_check.empty())
+            return;
+        
         LOG_H3("Free Variable Elimination for v" << v << " using equation " << c);
         pdd const& p = c.eq();
         SASSERT_EQ(p.degree(v), 1);
@@ -60,39 +189,69 @@ namespace polysat {
         LOG("lc: " << lc);
         LOG("rest: " << rest);
 
-        assignment_t a;
+        //pdd rs = rest;
+        
+        if (!lc.is_val() && lc.is_var())
+            // TODO: We could introduce a new variable "new_var = lc" and add the valuation for this new variable
+            return;
+        
+        pdd coeff_odd = get_odd(lc); // a'
+        optional<pdd> lci = get_inverse(coeff_odd); // a'^-1
+        if (!lci)
+            return; // For sure not invertible
+        
+        LOG("lci: " << lci);
+        
+        /*assignment_t a;
         pdd const lcs = eval(lc, core, a);
         LOG("lcs: " << lcs);
         pdd lci = m.zero();
         if (!inv(lcs, lci))
             return;
+            
+        pdd lci = m.zero();
 
         pdd const rs = s.subst(a, rest);
-        pdd const vs = -rs * lci;  // this is the polynomial that computes v
+        */
+        
+        /*pdd const vs = -rs * lci;  // this is the polynomial that computes v
         LOG("vs: " << vs);
-        SASSERT(!vs.free_vars().contains(v));
+        SASSERT(!vs.free_vars().contains(v));*/
 
         // Find another constraint where we want to substitute v
-        for (signed_constraint c_target : core) {
-            if (c == c_target)
-                continue;
-            if (c_target.vars().size() <= 1)
-                continue;
-            if (!c_target.contains_var(v))
-                continue;
-
-            // TODO: helper method constraint::subst(pvar v, pdd const& p)
-            //       (or rather, add it on constraint_manager since we need to allocate/dedup the new constraint)
-            //  For now, just restrict to ule_constraint.
-            if (!c_target->is_ule())
-                continue;
-            // TODO: maybe apply assignment a here as well
-            pdd const lhs = c_target->to_ule().lhs().subst_pdd(v, vs);
-            pdd const rhs = c_target->to_ule().rhs().subst_pdd(v, vs);
-            signed_constraint c_new = s.ule(lhs, rhs);
+        for (signed_constraint c_target : to_check) {
+            
+            // Move the variable to eliminate to one side
+            pdd fac = m.zero(); 
+            pdd fac2 = m.zero(); 
+            pdd rest2 = m.zero(); 
+            pdd rest3 = m.zero(); 
+            c_target->to_ule().lhs().factor(v, 1, fac, rest2);
+            c_target->to_ule().rhs().factor(v, 1, fac2, rest3);
+            
+            rest2 = rest3 - rest2;
+            fac = fac - fac2;
+            
+            LOG("c_target: " << lit_pp(s, c_target));
+            LOG("fac: " << fac);
+            LOG("rest2: " << rest2);
+            
+            if (!fac.is_val() && !fac.is_var())
+                return; // TODO: Again, we might bind this to a variable
+            
+            pdd pv2 = get_dyadic_valuation(fac).first;
+            pdd pv1 = get_dyadic_valuation(lc).first;  
+            pdd odd_fac = get_odd(fac);
+            pdd power_diff = s.shl(m.one(), pv2 - pv1);
+            
+            LOG("pv2: " << pv2);
+            LOG("pv1: " << pv1);
+            LOG("odd_fac: " << odd_fac);
+            LOG("power_diff: " << power_diff);
+            
+            signed_constraint c_new = s.ule(-rest * lci * power_diff * odd_fac, rest2);
             if (c_target.is_negative())
                 c_new.negate();
-            LOG("c_target: " << lit_pp(s, c_target));
             LOG("c_new:    " << lit_pp(s, c_new));
 
             // New constraint is already true (maybe we already derived it previously?)
@@ -102,12 +261,15 @@ namespace polysat {
                 continue;
 
             clause_builder cb(s);
-            for (auto [w, wv] : a)
-                cb.insert(~s.eq(s.var(w), wv));
+            /*for (auto [w, wv] : a)
+                cb.push(~s.eq(s.var(w), wv));*/
             cb.insert(~c);
             cb.insert(~c_target);
+            cb.insert(~s.ule(get_dyadic_valuation(lc).first, get_dyadic_valuation(fac).first));
             cb.insert(c_new);
-            core.add_lemma(cb.build());
+            ref<clause> c = cb.build();
+            if (c) // Can we get tautologies this way?
+                core.add_lemma(c);
         }
     }
 

src/math/polysat/variable_elimination.h

@@ -21,7 +21,19 @@ namespace polysat {
     class conflict;
 
     class free_variable_elimination {
+        
         solver& s;
+        unsigned_vector m_has_validation_of_range; // TODO: Find a better name
+        unsigned_vector m_pv_constants;
+        unsigned_vector m_pv_power_constants;
+        unsigned_vector m_inverse_constants;
+        unsigned_vector m_rest_constants;
+        
+        pdd get_odd(pdd p); // add lemma "2^pv(v) * v' = v"
+        optional<pdd> get_inverse(pdd v); // add lemma "v' * v'^-1 = 1 (where v' is the odd part of v)"
+        void add_dyadic_valuation(pvar v, unsigned k); // add lemma "pv(v) = k" ==> "pv_v = k"
+        std::pair<pdd, pdd> get_dyadic_valuation(pdd p, unsigned short lower, unsigned short upper); // lower bound inclusive; upper exclusive
+        std::pair<pdd, pdd> get_dyadic_valuation(pdd p);
         void find_lemma(pvar v, conflict& core);
         void find_lemma(pvar v, signed_constraint c, conflict& core);
         pdd eval(pdd const& p, conflict& core, assignment_t& out_assignment);