move constraint handler functionality to self-contained / separate classes.

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/math/polysat/constraint.cpp

@@ -12,16 +12,126 @@ Author:
 --*/
 
 #include "math/polysat/constraint.h"
+#include "math/polysat/solver.h"
+#include "math/polysat/log.h"
 
 namespace polysat {
 
-
+    eq_constraint& constraint::to_eq() { 
-    constraint* constraint::eq(unsigned lvl, pdd const& p, p_dependency_ref& d) {
+        return *dynamic_cast<eq_constraint*>(this); 
-        return alloc(eq_constraint, lvl, p, dep);
     }
 
-    std::ostream& constraint::display(std::ostream& out) const {
+    eq_constraint const& constraint::to_eq() const { 
+        return *dynamic_cast<eq_constraint const*>(this); 
+    }
+
+    constraint* constraint::eq(unsigned lvl, pdd const& p, p_dependency_ref& d) {
+        return alloc(eq_constraint, lvl, p, d);
+    }
+
+    std::ostream& eq_constraint::display(std::ostream& out) const {
         return out << p() << " == 0";
     }
 
+    bool eq_constraint::propagate(solver& s, pvar v) {
+        LOG_H3("Propagate " << s.m_vars[v] << " in " << *this);
+        SASSERT(!vars().empty());
+        auto var = s.m_vars[v].var();
+        unsigned idx = 0;
+        if (vars()[idx] != v)
+            idx = 1;
+        SASSERT(v == vars()[idx]);
+        // find other watch variable.
+        for (unsigned i = vars().size(); i-- > 2; ) {
+            if (!s.is_assigned(vars()[i])) {
+                std::swap(vars()[idx], vars()[i]);
+                return true;
+            }
+        }        
+        
+
+        LOG("Assignments: " << s.m_search);
+        auto q = p().subst_val(s.m_search);
+        LOG("Substituted: " << p() << " := " << q);
+        TRACE("polysat", tout << p() << " := " << q << "\n";);
+        if (q.is_zero()) 
+            return false;
+        if (q.is_never_zero()) {
+            LOG("Conflict (never zero under current assignment)");
+            // we could tag constraint to allow early substitution before 
+            // swapping watch variable in case we can detect conflict earlier.
+            s.set_conflict(*this);
+            return false;
+        }
+
+        // at most one variable remains unassigned.
+        auto other_var = vars()[1 - idx];
+
+        // Detect and apply unit propagation.
+            
+        if (!q.is_linear())
+            return false;
+
+        // a*x + b == 0
+        rational a = q.hi().val();
+        rational b = q.lo().val();
+        rational inv_a;
+        if (a.is_odd()) {
+            // v1 = -b * inverse(a)
+            unsigned sz = q.power_of_2();
+            VERIFY(a.mult_inverse(sz, inv_a)); 
+            rational val = mod(inv_a * -b, rational::power_of_two(sz));
+            s.m_cjust[other_var].push_back(this);
+            s.propagate(other_var, val, *this);
+            return false;
+        }
+
+        SASSERT(!b.is_odd());  // otherwise p.is_never_zero() would have been true above
+
+        // TBD
+        // constrain viable using condition on x
+        // 2*x + 2 == 0 mod 4 => x is odd
+        //
+        // We have:
+        // 2^j*a'*x + 2^j*b' == 0 mod m, where a' is odd (but not necessarily b')
+        // <=> 2^j*(a'*x + b') == 0 mod m
+        // <=> a'*x + b' == 0 mod (m-j)
+        // <=> x == -b' * inverse_{m-j}(a') mod (m-j)
+        // ( <=> 2^j*x == 2^j * -b' * inverse_{m-j}(a') mod m )
+        //
+        // x == c mod (m-j)
+        // Which x in 2^m satisfy this?
+        // => x \in { c + k * 2^(m-j) | k = 0, ..., 2^j - 1 }
+        unsigned rank_a = a.trailing_zeros();  // j
+        SASSERT(b == 0 || rank_a <= b.trailing_zeros());
+        rational aa = a / rational::power_of_two(rank_a);  // a'
+        rational bb = b / rational::power_of_two(rank_a);  // b'
+        rational inv_aa;
+        unsigned small_sz = q.power_of_2() - rank_a;  // m - j
+        VERIFY(aa.mult_inverse(small_sz, inv_aa));
+        rational cc = mod(inv_aa * -bb, rational::power_of_two(small_sz));
+        LOG(m_vars[other_var] << " = " << cc << " + k * 2^" << small_sz);
+        // TODO: better way to update the BDD, e.g. construct new one (only if rank_a is small?)
+        vector<rational> viable;
+        for (rational k = rational::zero(); k < rational::power_of_two(rank_a); k += 1) {
+            rational val = cc + k * rational::power_of_two(small_sz);
+            viable.push_back(val);
+        }
+        LOG_V("still viable: " << viable);
+        unsigned i = 0;
+        for (rational r = rational::zero(); r < rational::power_of_two(q.power_of_2()); r += 1) {
+            while (i < viable.size() && viable[i] < r)
+                ++i;
+            if (i < viable.size() && viable[i] == r)
+                continue;
+            if (s.is_viable(other_var, r)) {
+                s.add_non_viable(other_var, r);
+            }
+        }
+
+        LOG("TODO");
+        
+        return false;
+    }
+
 }

src/math/polysat/constraint.h

@@ -26,6 +26,8 @@ namespace polysat {
     class ule_constraint;
 
     class constraint {
+        friend class eq_constraint;
+        friend class ule_constraint;
         unsigned         m_level;
         ckind_t          m_kind;
         p_dependency_ref m_dep;
@@ -39,8 +41,9 @@ namespace polysat {
         bool is_sle() const { return m_kind == ckind_t::sle_t; }
         ckind_t kind() const { return m_kind; }
         virtual std::ostream& display(std::ostream& out) const = 0;
-        eq_constraint& to_eq() { return *dynamic_cast<eq_constraint*>(this); }
+        virtual bool propagate(solver& s, pvar v) = 0;
-        eq_constraint const& to_eq() { return *dynamic_cast<eq_constraint const*>(this); }
+        eq_constraint& to_eq();
+        eq_constraint const& to_eq() const;
         p_dependency* dep() const { return m_dep; }
         unsigned_vector& vars() { return m_vars; }
         unsigned level() const { return m_level; }
@@ -49,12 +52,14 @@ namespace polysat {
     class eq_constraint : public constraint {
         pdd m_poly;
     public:
-        eq_constraint(unsigned lvl, pdd const& p, pdd const& q, p_dependency_ref& dep,):
+        eq_constraint(unsigned lvl, pdd const& p, p_dependency_ref& dep):
             constraint(lvl, dep, ckind_t::eq_t), m_poly(p) {
             m_vars.append(p.free_vars());
         }
         pdd const & p() const { return m_poly; }
         std::ostream& display(std::ostream& out) const override;
+
+        bool propagate(solver& s, pvar v) override;
     };
 
     inline std::ostream& operator<<(std::ostream& out, constraint const& c) { return c.display(out); }

src/math/polysat/log.cpp

@@ -66,6 +66,7 @@ std::pair<std::ostream&, bool>
 polysat_log(LogLevel msg_level, std::string fn, std::string /* pretty_fn */)
 {
   std::ostream& os = std::cerr;
+#if 0
   int const fd = fileno(stderr);
 
   size_t width = 20;
@@ -82,6 +83,10 @@ polysat_log(LogLevel msg_level, std::string fn, std::string /* pretty_fn */)
   os << "[" << fn << "] " << std::string(padding, ' ');
   os << std::string(polysat_log_indent_level, ' ');
   return {os, (bool)color};
+#else
+  return {os, false};
+#endif
+
 }
 
 polysat_log_indent::polysat_log_indent(int amount)

src/math/polysat/solver.cpp

@@ -215,7 +215,7 @@ namespace polysat {
     bool solver::propagate(pvar v, constraint& c) {
         switch (c.kind()) {
         case ckind_t::eq_t:
-            return propagate_eq(v, c);
+            return c.to_eq().propagate(*this, v);
         case ckind_t::ule_t:
         case ckind_t::sle_t:
             NOT_IMPLEMENTED_YET();
@@ -225,110 +225,6 @@ namespace polysat {
         return false;
     }
 
-    bool solver::propagate_eq(pvar v, constraint& c) {
-        LOG_H3("Propagate " << m_vars[v] << " in " << c);
-        SASSERT(c.kind() == ckind_t::eq_t);
-        SASSERT(!c.vars().empty());
-        auto var = m_vars[v].var();
-        auto& vars = c.vars();
-        unsigned idx = 0;
-        if (vars[idx] != v)
-            idx = 1;
-        SASSERT(v == vars[idx]);
-        // find other watch variable.
-        for (unsigned i = vars.size(); i-- > 2; ) {
-            if (!is_assigned(vars[i])) {
-                std::swap(vars[idx], vars[i]);
-                return true;
-            }
-        }        
-        
-
-        LOG("Assignments: " << m_search);
-        auto p = c.p().subst_val(m_search);
-        LOG("Substituted: " << c.p() << " := " << p);
-        TRACE("polysat", tout << c.p() << " := " << p << "\n";);
-        if (p.is_zero()) 
-            return false;
-        if (p.is_never_zero()) {
-            LOG("Conflict (never zero under current assignment)");
-            // we could tag constraint to allow early substitution before 
-            // swapping watch variable in case we can detect conflict earlier.
-            set_conflict(c);
-            return false;
-        }
-
-        // at most one variable remains unassigned.
-        auto other_var = vars[1 - idx];
-        // SASSERT(!is_assigned(other_var));
-
-        // Detect and apply unit propagation.
-            
-        if (!p.is_linear())
-            return false;
-
-        // a*x + b == 0
-        rational a = p.hi().val();
-        rational b = p.lo().val();
-        rational inv_a;
-        if (a.is_odd()) {
-            // v1 = -b * inverse(a)
-            unsigned sz = p.power_of_2();
-            VERIFY(a.mult_inverse(sz, inv_a)); 
-            rational val = mod(inv_a * -b, rational::power_of_two(sz));
-            m_cjust[other_var].push_back(&c);
-            propagate(other_var, val, c);
-            return false;
-        }
-
-        SASSERT(!b.is_odd());  // otherwise p.is_never_zero() would have been true above
-
-        // TBD
-        // constrain viable using condition on x
-        // 2*x + 2 == 0 mod 4 => x is odd
-        //
-        // We have:
-        // 2^j*a'*x + 2^j*b' == 0 mod m, where a' is odd (but not necessarily b')
-        // <=> 2^j*(a'*x + b') == 0 mod m
-        // <=> a'*x + b' == 0 mod (m-j)
-        // <=> x == -b' * inverse_{m-j}(a') mod (m-j)
-        // ( <=> 2^j*x == 2^j * -b' * inverse_{m-j}(a') mod m )
-        //
-        // x == c mod (m-j)
-        // Which x in 2^m satisfy this?
-        // => x \in { c + k * 2^(m-j) | k = 0, ..., 2^j - 1 }
-        unsigned rank_a = a.trailing_zeros();  // j
-        SASSERT(b == 0 || rank_a <= b.trailing_zeros());
-        rational aa = a / rational::power_of_two(rank_a);  // a'
-        rational bb = b / rational::power_of_two(rank_a);  // b'
-        rational inv_aa;
-        unsigned small_sz = p.power_of_2() - rank_a;  // m - j
-        VERIFY(aa.mult_inverse(small_sz, inv_aa));
-        rational cc = mod(inv_aa * -bb, rational::power_of_two(small_sz));
-        LOG(m_vars[other_var] << " = " << cc << " + k * 2^" << small_sz);
-        // TODO: better way to update the BDD, e.g. construct new one (only if rank_a is small?)
-        vector<rational> viable;
-        for (rational k = rational::zero(); k < rational::power_of_two(rank_a); k += 1) {
-            rational val = cc + k * rational::power_of_two(small_sz);
-            viable.push_back(val);
-        }
-        LOG_V("still viable: " << viable);
-        unsigned i = 0;
-        for (rational r = rational::zero(); r < rational::power_of_two(p.power_of_2()); r += 1) {
-            while (i < viable.size() && viable[i] < r)
-                ++i;
-            if (i < viable.size() && viable[i] == r)
-                continue;
-            if (is_viable(other_var, r)) {
-                add_non_viable(other_var, r);
-            }
-        }
-
-        LOG("TODO");
-        
-        return false;
-    }
-
     void solver::propagate(pvar v, rational const& val, constraint& c) {
         if (is_viable(v, val)) {
             m_free_vars.del_var_eh(v);
@@ -645,8 +541,8 @@ namespace polysat {
         constraint* c = m_conflict[0];
         constraint* d = m_cjust[v].back();
         if (c->is_eq() && d->is_eq()) {
-            pdd p = c->p();
+            pdd p = c->to_eq().p();
-            pdd q = d->p();
+            pdd q = d->to_eq().p();
             pdd r = p;
             if (!p.resolve(v, q, r)) 
                 return nullptr;
@@ -668,13 +564,13 @@ namespace polysat {
 
     bool solver::is_always_false(constraint& c) {
         if (c.is_eq()) 
-            return c.p().is_never_zero();
+            return c.to_eq().p().is_never_zero();
         return false;
     }
 
     bool solver::eval_to_false(constraint& c) {
         if (c.is_eq()) {
-            pdd r = c.p().subst_val(m_search);
+            pdd r = c.to_eq().p().subst_val(m_search);
             return r.is_never_zero();
         }
         return false;

src/math/polysat/solver.h

@@ -32,6 +32,8 @@ namespace polysat {
             stats() { reset(); }
         };
 
+        friend class eq_constraint;
+
         typedef ptr_vector<constraint> constraints;
 
         trail_stack&             m_trail;
@@ -120,7 +122,6 @@ namespace polysat {
 
         void propagate(pvar v);
         bool propagate(pvar v, constraint& c);
-        bool propagate_eq(pvar v, constraint& c);
         void propagate(pvar v, rational const& val, constraint& c);
         void erase_watch(pvar v, constraint& c);
         void erase_watch(constraint& c);