optimizations, fixes, TODO items

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

src/math/dd/dd_pdd.cpp

@@ -34,6 +34,7 @@ namespace dd {
             s = mod2_e;
         m_semantics = s;
         m_mod2N = rational::power_of_two(power_of_2);
+        m_max_value = m_mod2N - 1;
         m_power_of_2 = power_of_2;
         unsigned_vector l2v;
         for (unsigned i = 0; i < num_vars; ++i) l2v.push_back(i);
@@ -1650,6 +1651,20 @@ namespace dd {
         return *this;
     }
 
+    pdd& pdd::operator=(rational const& k) {
+        m.dec_ref(root);
+        root = m.mk_val(k).root;
+        m.inc_ref(root);
+        return *this;
+    }
+
+    rational const& pdd::leading_coefficient() const {
+        pdd p = *this;
+        while (!p.is_val())
+            p = p.hi();
+        return p.val();
+    }
+
     std::ostream& operator<<(std::ostream& out, pdd const& b) { return b.display(out); }
 
     void pdd_iterator::next() {

src/math/dd/dd_pdd.h

@@ -205,8 +205,9 @@ namespace dd {
         unsigned_vector            m_free_values;
         rational                   m_freeze_value;
         rational                   m_mod2N;
-        unsigned                   m_power_of_2 { 0 };
+        unsigned                   m_power_of_2 = 0;
-        unsigned                   m_gc_generation { 0 };  ///< will be incremented on each GC
+        rational                   m_max_value;
+        unsigned                   m_gc_generation = 0;  ///< will be incremented on each GC
 
         void reset_op_cache();
         void init_nodes(unsigned_vector const& l2v);
@@ -364,6 +365,8 @@ namespace dd {
         unsigned num_vars() const { return m_var2pdd.size(); }
 
         unsigned power_of_2() const { return m_power_of_2; }
+        rational const& max_value() const { return m_max_value; }
+        rational const& two_to_N() const { return m_mod2N; }
 
         unsigned_vector const& free_vars(pdd const& p);
 
@@ -387,12 +390,14 @@ namespace dd {
         pdd(pdd && other) noexcept : root(0), m(other.m) { std::swap(root, other.root); }
         pdd& operator=(pdd const& other);
         pdd& operator=(unsigned k);
+        pdd& operator=(rational const& k);
         ~pdd() { m.dec_ref(root); }
         pdd lo() const { return pdd(m.lo(root), m); }
         pdd hi() const { return pdd(m.hi(root), m); }
         unsigned index() const { return root; }
         unsigned var() const { return m.var(root); }
         rational const& val() const { SASSERT(is_val()); return m.val(root); }
+        rational const& leading_coefficient() const;
         bool is_val() const { return m.is_val(root); }
         bool is_one() const { return m.is_one(root); }
         bool is_zero() const { return m.is_zero(root); }
@@ -472,10 +477,13 @@ namespace dd {
 
     inline pdd& operator&=(pdd & p, pdd const& q) { p = p & q; return p; }
     inline pdd& operator^=(pdd & p, pdd const& q) { p = p ^ q; return p; }
-    inline pdd& operator*=(pdd & p, pdd const& q) { p = p * q; return p; }
     inline pdd& operator|=(pdd & p, pdd const& q) { p = p | q; return p; }
+    inline pdd& operator*=(pdd & p, pdd const& q) { p = p * q; return p; }
     inline pdd& operator-=(pdd & p, pdd const& q) { p = p - q; return p; }
     inline pdd& operator+=(pdd & p, pdd const& q) { p = p + q; return p; }
+    inline pdd& operator*=(pdd & p, rational const& q) { p = p * q; return p; }
+    inline pdd& operator-=(pdd & p, rational const& q) { p = p - q; return p; }
+    inline pdd& operator+=(pdd & p, rational const& q) { p = p + q; return p; }
 
     inline void swap(pdd& p, pdd& q) { p.swap(q); }
 

src/math/polysat/conflict.cpp

@@ -16,9 +16,18 @@ Notes:
 
  TODO: try a final core reduction step or other core minimization
 
-TODO: If we have e.g. 4x+y=2 and y=0, then we have a conflict no matter the value of x, so we should drop x=? from the core.
+ TODO: If we have e.g. 4x+y=2 and y=0, then we have a conflict no matter the value of x, so we should drop x=? from the core.
        (works currently if x is unassigned; for other cases we would need extra info from constraint::is_currently_false)
 
+ TODO: keep is buggy. The assert 
+                   SASSERT(premise.is_currently_true(s()) || premise.bvalue(s()) == l_true);
+       does not necessarily hold. A saturation premise could be inserted that is a resolvent that evaluates to false
+       and therefore not a current Boolean literal on the search stack.
+
+ TODO: revert(pvar v) is too weak. 
+       It should apply saturation rules currently only available for for propagated values.
+
+ TODO: dependency tracking for constraints evaluating to false should be minimized.
 --*/
 
 #include "math/polysat/conflict.h"

src/math/polysat/saturation.cpp

@@ -94,48 +94,47 @@ namespace polysat {
     void inf_saturate::push_omega_bisect(vector<signed_constraint>& new_constraints, pdd const& x, rational x_max, pdd const& y, rational y_max) {
         rational x_val, y_val;
         auto& pddm = x.manager();
-        unsigned bit_size = pddm.power_of_2();
+        rational bound = pddm.max_value();
-        rational bound = rational::power_of_two(bit_size);
         VERIFY(s().try_eval(x, x_val));
         VERIFY(s().try_eval(y, y_val));
-        SASSERT(x_val * y_val < bound);
+        SASSERT(x_val * y_val <= bound);
 
         rational x_lo = x_val, x_hi = x_max, y_lo = y_val, y_hi = y_max;
         rational two(2);
         while (x_lo < x_hi || y_lo < y_hi) {
             rational x_mid = div(x_hi + x_lo, two);
             rational y_mid = div(y_hi + y_lo, two);
-            if (x_mid * y_mid >= bound)
+            if (x_mid * y_mid > bound)
                 x_hi = x_mid - 1, y_hi = y_mid - 1;
             else
                 x_lo = x_mid, y_lo = y_mid;
         }
         SASSERT(x_hi == x_lo && y_hi == y_lo);
-        SASSERT(x_lo * y_lo < bound);
+        SASSERT(x_lo * y_lo <= bound);
-        SASSERT((x_lo + 1) * (y_lo + 1) >= bound);
+        SASSERT((x_lo + 1) * (y_lo + 1) > bound);
-        if ((x_lo + 1) * y_lo < bound) {
+        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);
-                if (x_mid * y_lo >= bound)
+                if (x_mid * y_lo > bound)
                     x_hi = x_mid - 1;
                 else
                     x_lo = x_mid;
             }
         }
-        else if (x_lo * (y_lo + 1) < bound) {
+        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);
-                if (y_mid * x_lo >= bound)
+                if (y_mid * x_lo > bound)
                     y_hi = y_mid - 1;
                 else
                     y_lo = y_mid;
             }
         }
-        SASSERT(x_lo * y_lo < bound);
+        SASSERT(x_lo * y_lo <= bound);
-        SASSERT((x_lo + 1) * y_lo >= bound);
+        SASSERT((x_lo + 1) * y_lo > bound);
-        SASSERT(x_lo * (y_lo + 1) >= bound);
+        SASSERT(x_lo * (y_lo + 1) > bound);
 
         // inequalities are justified by current assignments to x, y
         // conflict resolution should be able to pick up this as a valid justification.
@@ -151,17 +150,15 @@ namespace polysat {
     // then extract premises for a non-worst-case bound.
     void inf_saturate::push_omega(vector<signed_constraint>& new_constraints, pdd const& x, pdd const& y) {     
         auto& pddm = x.manager();
-        unsigned bit_size = pddm.power_of_2();
+        rational x_max = pddm.max_value();
-        rational bound = rational::power_of_two(bit_size);
+        rational y_max = pddm.max_value();
-        rational x_max = bound - 1;
-        rational y_max = bound - 1;
         
         if (x.is_var()) 
             x_max = s().m_viable.max_viable(x.var());
         if (y.is_var()) 
             y_max = s().m_viable.max_viable(y.var());        
 
-        if (x_max * y_max  >= bound)            
+        if (x_max * y_max  > pddm.max_value())            
             push_omega_bisect(new_constraints, x, x_max, y, y_max);
         else {
             for (auto c : s().m_cjust[y.var()])

src/math/polysat/ule_constraint.cpp

@@ -12,15 +12,16 @@ Author:
 
 Notes:
 
-   TODO: add rewrite rules to simplify expressions
+   Rewrite rules to simplify expressions
   
    - k1 <= k2   ==> 0 <= 0 if k1 <= k2
    - k1 <= k2   ==> 1 <= 0 if k1 >  k2
    - 0 <= p     ==> 0 <= 0
    - p <= -1    ==> 0 <= 0
-   - k*p <= 0   ==> p <= 0 if k is odd
+   - k*2^n*p <= 0   ==> 2^n*p <= 0 if k is odd, leading coeffient is always a power of 2.
+   - k <= p     ==> p - k <= - k - 1
 
-   - k <= p     ==> p - k <= k - 1
+  TODO:
    - p <= p + q ==> p <= -q - 1
    - p + k <= p ==> p + k <= k - 1   for k > 0
 
@@ -45,11 +46,35 @@ namespace polysat {
     }
 
     void ule_constraint::simplify() {
+        if (m_lhs.is_zero()) {
+            m_rhs = 0;
+            return;
+        }
+        if (m_rhs.is_val() && m_rhs.val() == m_rhs.manager().max_value()) {
+            m_lhs = 0, m_rhs = 0;
+            return;
+        }        
         if (m_lhs.is_val() && m_rhs.is_val()) {
             if (m_lhs.val() <= m_rhs.val())
                 m_lhs = m_rhs = 0;
             else
                 m_lhs = 1, m_rhs = 0;
+            return;
+        }
+        // k <= p => p - k <= - k - 1
+        if (m_lhs.is_val()) {
+            pdd k = m_lhs;
+            m_lhs = m_rhs - k;
+            m_rhs = - k - 1;
+        }
+        // normalize leading coefficient to be a power of 2
+        if (m_rhs.is_zero() && !m_lhs.leading_coefficient().is_power_of_two()) {
+            rational lc = m_lhs.leading_coefficient();
+            rational x, y;
+            gcd(lc, m_lhs.manager().two_to_N(), x, y);
+            if (x.is_neg())
+                x = mod(x, m_lhs.manager().two_to_N());
+            m_lhs *= x;
         }
     }
 

src/util/rational.h

@@ -339,7 +339,12 @@ public:
 
     static rational power_of_two(unsigned k);
 
-    bool is_power_of_two(unsigned & shift) {
+    bool is_power_of_two(unsigned & shift) const {
+        return m().is_power_of_two(m_val, shift);
+    }
+    
+    bool is_power_of_two() const {
+        unsigned shift = 0;
         return m().is_power_of_two(m_val, shift);
     }