move to interval arithmetic

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

src/math/polysat/fixplex.h

@@ -29,6 +29,45 @@ namespace polysat {
 
     typedef unsigned var_t;
 
+    template<typename Numeral>
+    struct pp {
+        Numeral n;
+        pp(Numeral const& n):n(n) {}
+    };
+
+    template<typename Numeral>
+    inline std::ostream& operator<<(std::ostream& out, pp<Numeral> const& p) {
+        if ((0 - p.n) < p.n)
+            return out << "-" << (0 - p.n);
+        return out << p.n;
+    }
+
+    /**
+     * Modular interval arithmetic
+     */
+    template<typename Numeral>
+    struct interval {
+        Numeral lo { 0 };
+        Numeral hi { 0 };
+        interval() {}
+        interval(Numeral const& l, Numeral const& h): lo(l), hi(h) {}
+        static interval free() { return interval(0, 0); }
+        bool is_free() const { return lo == hi; }
+        interval operator+(interval const& other) const;
+        interval operator-(interval const& other) const;
+        interval operator-() const;
+        interval operator*(Numeral const& n) const;
+        interval operator+(Numeral const& n) const { return interval(lo + n, hi + n); }
+        interval operator-(Numeral const& n) const { return interval(lo - n, hi - n); }
+        interval& operator+=(interval const& other) { *this = *this + other; return *this; }
+        std::ostream& display(std::ostream& out) const { return out << "[" << pp(lo) << ", " << pp(hi) << "["; }
+    };
+
+    template<typename Numeral>
+    inline std::ostream& operator<<(std::ostream& out, interval<Numeral> const& i) {
+        return i.display(out);
+    }
+
     template<typename Ext>
     class fixplex {
     public:
@@ -62,11 +101,9 @@ namespace polysat {
             S_DEFAULT
         };
 
-        struct var_info {
+        struct var_info : public interval<numeral> {
             unsigned    m_base2row:29;
             unsigned    m_is_base:1;
-            numeral     m_lo { 0 };
-            numeral     m_hi { 0 };
             numeral     m_value { 0 };
             var_info():
                 m_base2row(0),
@@ -124,11 +161,11 @@ namespace polysat {
 
 
         void set_bounds(var_t v, numeral const& lo, numeral const& hi);
-        void unset_bounds(var_t v) { m_vars[v].m_lo = m_vars[v].m_hi; }
+        void unset_bounds(var_t v) { m_vars[v].lo = m_vars[v].hi; }
 
         var_t get_base_var(row const& r) const { return m_rows[r.id()].m_base; }
-        numeral const& lo(var_t var) const { return m_vars[var].m_lo; }
+        numeral const& lo(var_t var) const { return m_vars[var].lo; }
-        numeral const& hi(var_t var) const { return m_vars[var].m_hi; }
+        numeral const& hi(var_t var) const { return m_vars[var].hi; }
         numeral const& value(var_t var) const { return m_vars[var].m_value; }
         void set_max_iterations(unsigned n) { m_max_iterations = n; }
         unsigned get_num_vars() const { return m_vars.size(); }
@@ -215,9 +252,11 @@ namespace polysat {
 
     };
 
+
     struct uint64_ext {
         typedef uint64_t numeral;
         static const uint64_t max_numeral = 0; // std::limits<uint64_t>::max();
+
         struct manager {
             typedef uint64_t numeral;
             struct hash {
@@ -258,16 +297,18 @@ namespace polysat {
                 return x <= r;
             }
             std::ostream& display(std::ostream& out, numeral const& x) const { 
-                return out << x; 
+                return out << pp(x); 
             }
         };
         typedef _scoped_numeral<manager> scoped_numeral;
     };
 
+
     template<typename Ext>
     inline std::ostream& operator<<(std::ostream& out, fixplex<Ext> const& fp) {
         return fp.display(out);
     }
 
+
 };
 

src/math/polysat/fixplex_def.h

@@ -23,6 +23,54 @@ Author:
 
 namespace polysat {
 
+    template<typename Numeral>
+    interval<Numeral> interval<Numeral>::operator+(interval<Numeral> const& other) const {
+        if (is_free())
+            return *this;
+        if (other.is_free())
+            return interval::free();
+        Numeral sz = (hi - lo) + (other.hi - other.lo);
+        if (sz < (hi - lo))
+            return interval::free();
+        return interval(lo + other.lo, hi + other.hi);
+    }
+
+    template<typename Numeral>
+    interval<Numeral> interval<Numeral>::operator-(interval<Numeral> const& other) const {
+        return *this + (-other);
+    }
+
+    template<typename Numeral>
+    interval<Numeral> interval<Numeral>::operator-() const {
+        if (is_free())
+            return *this;
+        return interval(1 - hi, 1 - lo);
+    }
+
+    template<typename Numeral>
+    interval<Numeral> interval<Numeral>::operator*(Numeral const& n) const {
+        if (n == 0)
+            return interval(0, 1);
+        if (n == 1)
+            return *this;
+        if (is_free())
+            return *this;
+        Numeral sz = hi - lo;
+        if (0 - n < n) {
+            Numeral mn = 0 - n;
+            Numeral mz = mn * sz;            
+            if (mz / mn != sz)
+                return interval::free();
+            return interval((hi - 1) * n, n * lo + 1);
+        }
+        else {
+            Numeral mz = n * sz;
+            if (mz / n != sz)
+                return interval::free();
+            return interval(n * lo, n * (hi - 1) + 1);
+        }
+    }
+
     template<typename Ext>
     fixplex<Ext>::~fixplex() {
         reset();
@@ -309,7 +357,7 @@ namespace polysat {
 
     /**
      * Compute delta to add to the value, such that value + delta is either lo(v), or hi(v) - 1
-     * A pre-condition is that value is not in the range [lo(v),hi(v)[, 
+     * A pre-condition is that value is not in the interval [lo(v),hi(v)[, 
      * and therefore as a consequence lo(v) != hi(v).
      */
     template<typename Ext>
@@ -333,8 +381,8 @@ namespace polysat {
      */
     template<typename Ext>
     void fixplex<Ext>::set_bounds(var_t v, numeral const& lo, numeral const& hi) {         
-        m_vars[v].m_lo = lo; 
+        m_vars[v].lo = lo; 
-        m_vars[v].m_hi = hi; 
+        m_vars[v].hi = hi; 
         if (in_bounds(v))
             return;
         if (is_base(v)) 
@@ -383,21 +431,15 @@ namespace polysat {
     bool fixplex<Ext>::is_infeasible_row(var_t x) {
         SASSERT(is_base(x));
         auto r = base2row(x);
-        numeral lo_sum = 0, hi_sum = 0, diff = 0;
+        interval<numeral> range(0, 1);
         for (auto const& e : M.row_entries(row(r))) {
             var_t v = e.var();
             numeral const& c = e.coeff();
-            if (lo(v) == hi(v))
+            range += m_vars[v] * c;
-                return false;
+            if (range.is_free())
-            lo_sum += lo(v) * c;
-            hi_sum += (hi(v) - 1) * c;
-            numeral range = hi(v) - lo(v);
-            if (!m.signed_mul(range, c, range))
-                return false;
-            if (!m.signed_add(diff, diff, range))
                 return false;
         }        
-        return 0 < lo_sum && lo_sum <= hi_sum;
+        return 0 < range.lo && range.lo < range.hi;
     }
 
     /**
@@ -744,7 +786,8 @@ namespace polysat {
 
     template<typename Ext>
     void fixplex<Ext>::propagate_bounds(row const& r) {
-        numeral lo_sum = 0, hi_sum = 0, diff = 0, free_c = 0;
+        interval<numeral> range(0, 1);
+        numeral free_c = 0;
         var_t free_v = null_var;
         for (auto const& e : M.row_entries(r)) {
             var_t v = e.var();
@@ -752,34 +795,22 @@ namespace polysat {
             if (is_free(v)) {
                 if (free_v != null_var)
                     return;
-                if (c != 1 && c + 1 != 0)
-                    return;
                 free_v = v;
                 free_c = c;
                 continue;
             }
-            lo_sum += lo(v) * c;
+            range += m_vars[v] * c;
-            hi_sum += (hi(v) - 1) * c;
+            if (range.is_free())
-            numeral range = hi(v) - lo(v);
-            if (!m.signed_mul(range, c, range))
-                return;
-            if (!m.signed_add(diff, diff, range))
                 return;
         }        
 
+        std::cout << "bounds " << free_v << " range " << range << "\n";
+
         if (free_v != null_var) {
-
+            range = (-range) * free_c;
-            //
+            if (range.is_free())
-            // c = 1:
+                return;
-            //   free_v in [0 - hi_sum, 1 - lo_sum[
+            new_bound(r, free_v, range.lo, range.hi);
-            // c = -1:
-            //   free_v in [lo_sum, hi_sum + 1[
-            //
-
-            if (free_c == 1) 
-                new_bound(r, free_v, 0 - hi_sum, 1 - lo_sum);
-            else 
-                new_bound(r, free_v, lo_sum, hi_sum + 1);
             SASSERT(in_bounds(free_v));
             return;
         }
@@ -787,23 +818,19 @@ namespace polysat {
             var_t v = e.var();
             SASSERT(!is_free(v));
             numeral const& c = e.coeff();
-            if (c != 1 && c + 1 != 0)
+            auto range1 = range - m_vars[v] * c;
-                continue;
+            new_bound(r, v, range1.lo, range1.hi);
-            numeral lo_sum1 = lo_sum - lo(v) * c;
+            // SASSERT(in_bounds(v));
-            numeral hi_sum1 = hi_sum - (hi(v) - 1) * c;
-            if (c == 1)
-                new_bound(r, v, 0 - hi_sum1, 1 - lo_sum1);
-            else
-                new_bound(r, v, lo_sum1, hi_sum1 + 1);
-            SASSERT(in_bounds(v));
         }
     }
 
     template<typename Ext>
     void fixplex<Ext>::new_bound(row const& r, var_t x, numeral const& l, numeral const& h) {
+        // TBD: 
+        std::cout << is_free(x) << " " << l << " " << lo(x) << " " << hi(x) << " " << h << "\n";
         if (!is_free(x) && l <= lo(x) && hi(x) <= h)
             return;            
-        IF_VERBOSE(0, verbose_stream() << "new-bound v" << x << " [" << l << "," << h << "[\n");
+        IF_VERBOSE(0, verbose_stream() << "new-bound v" << x << " [" << pp(l) << "," << pp(h) << "[\n");
     }
 
     template<typename Ext>    
@@ -811,7 +838,7 @@ namespace polysat {
         M.display(out);
         for (unsigned i = 0; i < m_vars.size(); ++i) {
             var_info const& vi = m_vars[i];
-            out << "v" << i << " " << value(i) << " [" << lo(i) << ", " << hi(i) << "[ ";
+            out << "v" << i << " " << pp(value(i)) << " [" << pp(lo(i)) << ", " << pp(hi(i)) << "[ ";
             if (vi.m_is_base) out << "b:" << vi.m_base2row << " ";
             out << "\n";
         }
@@ -823,10 +850,10 @@ namespace polysat {
         for (auto const& e : M.row_entries(r)) {
             var_t v = e.var();
             if (e.coeff() != 1)
-                out << e.coeff() << " * ";
+                out << pp(e.coeff()) << " * ";
             out << "v" << v << " ";
             if (values) 
-                out << value(v) << " [" << lo(v) << ", " << hi(v) << "[ ";            
+                out << pp(value(v)) << " [" << pp(lo(v)) << ", " << pp(hi(v)) << "[ ";            
         }
         return out << "\n";
     }

src/test/fixplex.cpp

@@ -71,6 +71,25 @@ namespace polysat {
         fp.set_bounds(y, 3, 6);
         fp.run();
         fp.propagate_bounds();
+        fp.reset();
+        fp.add_row(x, 3, ys, coeffs);
+        fp.set_bounds(x, 3, 4);
+        fp.set_bounds(y, 3, 6);
+        fp.set_bounds(z, 1, 8);
+        fp.run();
+        fp.propagate_bounds();
+        fp.reset();
+    }
+
+    static void test_interval1() {
+        interval<uint64_t> i1(1, 2);
+        interval<uint64_t> i2(3, 6);
+        std::cout << i1 << " " << i2 << "\n";
+        std::cout << i1 << " * 4 := " << (i1 * 4) << "\n";
+        std::cout << i2 << " * 3 := " << (i2 * 3) << "\n";
+        std::cout << i1 << " * -4 := " << (i1 * (0 - 4)) << "\n";
+        std::cout << i2 << " * -3 := " << (i2 * (0 - 3)) << "\n";
+        std::cout << "-" << i2 << " := " << (-i2) << "\n";
     }
 
 }
@@ -80,4 +99,5 @@ void tst_fixplex() {
     polysat::test2();
     polysat::test3();
     polysat::test4();
+    polysat::test_interval1();
 }