move arithmetical mbp functionality to model_based_opt

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

src/api/z3_api.h

@@ -3068,7 +3068,8 @@ extern "C" {
        \brief Create a numeral of a given sort.
 
        \param c logical context.
-       \param numeral A string representing the numeral value in decimal notation. If the given sort is a real, then the numeral can be a rational, that is, a string of the form \ccode{[num]* / [num]*}.
+       \param numeral A string representing the numeral value in decimal notation. The string may be of the form \code{[num]*[.[num]*][E[+|-][num]+]}.
+                      If the given sort is a real, then the numeral can be a rational, that is, a string of the form \ccode{[num]* / [num]*}.                      
        \param ty The sort of the numeral. In the current implementation, the given sort can be an int, real, finite-domain, or bit-vectors of arbitrary size.
 
        \sa Z3_mk_int

src/ast/arith_decl_plugin.h

@@ -422,5 +422,103 @@ public:
     expr_ref mk_add_simplify(unsigned sz, expr* const* args);
 };
 
+
+inline app_ref mk_numeral(rational const& r, app_ref const& x) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_numeral(r, r.is_int() && a.is_int(x)), x.get_manager());
+}
+
+inline app_ref operator+(app_ref const& x, app_ref const& y) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_add(x, y), x.get_manager());
+}
+
+inline app_ref operator+(app_ref const& x, rational const& y) {
+    return x + mk_numeral(y, x);
+}
+
+inline app_ref operator+(app_ref const& x, int y) {
+    return x + rational(y);
+}
+
+inline app_ref operator+(rational const& x, app_ref const& y) {
+    return mk_numeral(x, y) + y;
+}
+
+inline app_ref operator+(int x, app_ref const& y) {
+    return rational(x) + y;
+}
+
+inline app_ref operator-(app_ref const& x, app_ref const& y) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_sub(x, y), x.get_manager());
+}
+
+inline app_ref operator-(app_ref const& x, rational const& y) {
+    return x - mk_numeral(y, x);
+}
+
+inline app_ref operator-(app_ref const& x, int y) {
+    return x - rational(y);
+}
+
+inline app_ref operator-(rational const& x, app_ref const& y) {
+    return mk_numeral(x, y) - y;
+}
+
+inline app_ref operator-(int x, app_ref const& y) {
+    return rational(x) - y;
+}
+
+
+inline app_ref operator*(app_ref const& x, app_ref const& y) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_mul(x, y), x.get_manager());
+}
+
+inline app_ref operator*(app_ref const& x, rational const& y) {
+    return x * mk_numeral(y, x);
+}
+
+inline app_ref operator*(rational const& x, app_ref const& y) {
+    return mk_numeral(x, y) * y;
+}
+
+inline app_ref operator*(app_ref const& x, int y) {
+    return x * rational(y);
+}
+
+inline app_ref operator*(int x, app_ref const& y) {
+    return rational(x) * y;
+}
+
+inline app_ref operator<=(app_ref const& x, app_ref const& y) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_le(x, y), x.get_manager());
+}
+
+inline app_ref operator<=(app_ref const& x, rational const& y) {
+    return x <= mk_numeral(y, x);
+}
+
+inline app_ref operator<=(app_ref const& x, int y) {
+    return x <= rational(y);
+}
+
+inline app_ref operator>=(app_ref const& x, app_ref const& y) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_ge(x, y), x.get_manager());
+}
+
+inline app_ref operator<(app_ref const& x, app_ref const& y) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_lt(x, y), x.get_manager());
+}
+
+inline app_ref operator>(app_ref const& x, app_ref const& y) {
+    arith_util a(x.get_manager());
+    return app_ref(a.mk_gt(x, y), x.get_manager());
+}
+
 #endif /* ARITH_DECL_PLUGIN_H_ */
 

src/math/simplex/model_based_opt.cpp

@@ -7,7 +7,7 @@ Module Name:
 
 Abstract:
 
-    Model-based optimization for linear real arithmetic.
+    Model-based optimization and projection for linear real, integer arithmetic.
 
 Author:
 
@@ -26,6 +26,7 @@ std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) {
     case opt::t_eq: return out << " = ";
     case opt::t_lt: return out << " < ";
     case opt::t_le: return out << " <= ";        
+    case opt::t_mod: return out << " mod ";
     }
     return out;
 }
@@ -55,8 +56,8 @@ namespace opt {
         vector<var> const& vars = r.m_vars;
         for (unsigned i = 0; i < vars.size(); ++i) {
             // variables in each row are sorted and have non-zero coefficients
-            SASSERT(i + 1 == vars.size() || vars[i].m_id < vars[i+1].m_id);
-            SASSERT(!vars[i].m_coeff.is_zero());
+            PASSERT(i + 1 == vars.size() || vars[i].m_id < vars[i+1].m_id);
+            PASSERT(!vars[i].m_coeff.is_zero());
         }
         
         PASSERT(r.m_value == get_row_value(r));
@@ -64,6 +65,7 @@ namespace opt {
         // values satisfy constraints
         PASSERT(index == 0 || r.m_type != t_lt ||  r.m_value.is_neg());
         PASSERT(index == 0 || r.m_type != t_le || !r.m_value.is_pos());        
+        PASSERT(index == 0 || r.m_type != t_mod || (mod(r.m_value, r.m_mod).is_zero()));
         return true;
     }
         
@@ -117,7 +119,7 @@ namespace opt {
                 // objective + t2*coeff/a2 <= ub
 
                 mul_add(false, m_objective_id, - coeff/bound_coeff, bound_row_index);
-                m_rows[bound_row_index].m_alive = false;
+                retire_row(bound_row_index);
                 bound_trail.push_back(bound_row_index);
                 bound_vars.push_back(x);
             }
@@ -281,6 +283,12 @@ namespace opt {
         }
         return bound_row_index != UINT_MAX;
     }
+
+    void model_based_opt::retire_row(unsigned row_id) {
+        m_rows[row_id].m_alive = false;
+        m_retired_rows.push_back(row_id);
+    }
+
        
     rational model_based_opt::get_row_value(row const& r) const {
         vector<var> const& vars = r.m_vars;
@@ -308,7 +316,7 @@ namespace opt {
             }
             if (id < var_id) {
                 lo = mid + 1;
-            }			
+            }
             else {
                 hi = mid;
             }
@@ -355,10 +363,161 @@ namespace opt {
                 
         if (m_rows[row_dst].m_alive) {
             rational a2 = get_coefficient(row_dst, x);
-            mul_add(row_dst != m_objective_id && a1.is_pos() == a2.is_pos(), row_dst, -a2/a1, row_src);            
+            if (is_int(x)) {
+                TRACE("opt", 
+                      tout << a1 << " " << a2 << ": ";
+                      display(tout, m_rows[row_dst]);
+                      display(tout, m_rows[row_src]););
+                if (a1.is_pos() != a2.is_pos()) {  
+                    mul_add(x, a1, row_src, a2, row_dst);
+                }
+                else {
+                    mul(row_dst, abs(a1));
+                    mul_add(false, row_dst, -abs(a2), row_src);
+                }
+                TRACE("opt", display(tout, m_rows[row_dst]););
+                normalize(row_dst);
+            }
+            else {
+                mul_add(row_dst != m_objective_id && a1.is_pos() == a2.is_pos(), row_dst, -a2/a1, row_src);            
+            }
         }
     }
-    
+
+    // resolution for integer rows.
+    void model_based_opt::mul_add(
+        unsigned x, rational const& src_c, unsigned row_src, rational const& dst_c, unsigned row_dst) {
+        row& dst = m_rows[row_dst];
+        row const& src = m_rows[row_src];
+        SASSERT(is_int(x));
+        SASSERT(t_le == dst.m_type && t_le == src.m_type);
+        SASSERT(src_c.is_int());
+        SASSERT(dst_c.is_int());
+
+        rational abs_src_c = abs(src_c);
+        rational abs_dst_c = abs(dst_c);            
+        rational x_val = m_var2value[x];
+        rational slack = (abs_src_c - rational::one()) * (abs_dst_c - rational::one());
+        rational dst_val = dst.m_value - x_val*dst_c;
+        rational src_val = src.m_value - x_val*src_c;
+        bool use_case1 = 
+            (src_c * dst_val + dst_c * src_val + slack).is_nonpos() 
+            || abs_src_c.is_one() 
+            || abs_dst_c.is_one();
+
+        if (use_case1) {
+            // dst <- abs_src_c*dst + abs_dst_c*src - slack
+            mul(row_dst, abs_src_c);
+            sub(row_dst, slack);
+            mul_add(false, row_dst, abs_dst_c, row_src);
+            return;
+        }
+
+        //
+        // create finite disjunction for |b|.                                
+        //    exists x, z in [0 .. |b|-2] . b*x + s + z = 0 && ax + t <= 0 && bx + s <= 0
+        // <=> 
+        //    exists x, z in [0 .. |b|-2] . b*x = -z - s && ax + t <= 0 && bx + s <= 0
+        // <=>
+        //    exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && bx + s <= 0
+        // <=>
+        //    exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z - s + s <= 0
+        // <=>
+        //    exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z <= 0
+        // <=>
+        //    exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0
+        // <=>
+        //    exists x, z in [0 .. |b|-2] . b*x = -z - s && a*n_sign(b)(s + z) + |b|t <= 0
+        // <=>
+        //    exists z in [0 .. |b|-2] . |b| | (z + s) && a*n_sign(b)(s + z) + |b|t <= 0
+        //
+
+        vector<var> coeffs;
+        if (abs_dst_c <= abs_src_c) {
+            rational z = mod(dst_val, abs_dst_c);
+            if (!z.is_zero()) z = abs_dst_c - z;
+            mk_coeffs_without(coeffs, dst.m_vars, x); 
+            add_divides(coeffs, dst.m_coeff + z, abs_dst_c);
+            add(row_dst, z);
+            mul(row_dst, src_c * n_sign(dst_c));
+            mul_add(false, row_dst, abs_dst_c, row_src);            
+        }
+        else {
+            // z := b - (s + bx) mod b 
+            //   := b - s mod b
+            // b | s + z <=> b | s + b - s mod b <=> b | s - s mod b
+            rational z = mod(src_val, abs_src_c);
+            if (!z.is_zero()) z = abs_src_c - z;      
+            mk_coeffs_without(coeffs, src.m_vars, x);
+            add_divides(coeffs, src.m_coeff + z, abs_src_c);
+            mul(row_dst, abs_src_c);
+            add(row_dst, z * dst_c * n_sign(src_c));            
+            mul_add(false, row_dst, dst_c * n_sign(src_c), row_src);
+        }
+    }
+
+    void model_based_opt::mk_coeffs_without(vector<var>& dst, vector<var> const src, unsigned x) {
+        for (unsigned i = 0; i < src.size(); ++i) {
+            if (src[i].m_id != x) dst.push_back(src[i]);
+        }
+    }
+
+
+    rational model_based_opt::n_sign(rational const& b) const {
+        return rational(b.is_pos()?-1:1);
+    }
+   
+    void model_based_opt::mul(unsigned dst, rational const& c) {
+        if (c.is_one()) return;
+        row& r = m_rows[dst];
+        for (unsigned i = 0; i < r.m_vars.size(); ++i) {
+            r.m_vars[i].m_coeff *= c;
+        }
+        r.m_coeff *= c;
+        r.m_value *= c;
+    }
+
+    void model_based_opt::add(unsigned dst, rational const& c) {
+        row& r = m_rows[dst];
+        r.m_coeff += c;
+        r.m_value += c;
+    }
+
+    void model_based_opt::sub(unsigned dst, rational const& c) {
+        row& r = m_rows[dst];
+        r.m_coeff -= c;
+        r.m_value -= c;
+    }
+
+    void model_based_opt::normalize(unsigned row_id) {
+        row& r = m_rows[row_id];
+        if (r.m_vars.empty()) return;
+        if (r.m_type == t_mod) return;
+        rational g(abs(r.m_vars[0].m_coeff));
+        bool all_int = g.is_int();
+        for (unsigned i = 1; all_int && !g.is_one() && i < r.m_vars.size(); ++i) {
+            rational const& coeff = r.m_vars[i].m_coeff;
+            if (coeff.is_int()) {
+                g = gcd(g, abs(coeff));
+            }
+            else {
+                all_int = false;
+            }
+        }
+        if (all_int && !r.m_coeff.is_zero()) {
+            if (r.m_coeff.is_int()) {
+                g = gcd(g, abs(r.m_coeff));
+            }
+            else {
+                all_int = false;
+            }
+        }
+        if (all_int && !g.is_one()) {
+            SASSERT(!g.is_zero());
+            mul(row_id, rational::one()/g);
+        }
+    }
+ 
     //
     // set row1 <- row1 + c*row2
     //
@@ -452,12 +611,18 @@ namespace opt {
         else if (r.m_coeff.is_neg()) {
             out << r.m_coeff << " ";
         }        
-        out << r.m_type << " 0; value: " << r.m_value  << "\n";
+        if (r.m_type == opt::t_mod) {
+            out << r.m_type << " " << r.m_mod << " = 0; value: " << r.m_value  << "\n";
+        }
+        else {
+            out << r.m_type << " 0; value: " << r.m_value  << "\n";
+        }
     }
 
-    unsigned model_based_opt::add_var(rational const& value) {
+    unsigned model_based_opt::add_var(rational const& value, bool is_int) {
         unsigned v = m_var2value.size();
         m_var2value.push_back(value);
+        m_var2is_int.push_back(is_int);
         m_var2row_ids.push_back(unsigned_vector());
         return v;
     }
@@ -466,34 +631,65 @@ namespace opt {
         return m_var2value[var];
     }
 
-    void model_based_opt::set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, ineq_type rel) {
+    void model_based_opt::set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel) {
         row& r = m_rows[row_id];
         rational val(c);
         SASSERT(r.m_vars.empty());
         r.m_vars.append(coeffs.size(), coeffs.c_ptr());
+        bool is_int_row = true;
         std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
         for (unsigned i = 0; i < coeffs.size(); ++i) {
             val += m_var2value[coeffs[i].m_id] * coeffs[i].m_coeff;
+            SASSERT(!is_int(coeffs[i].m_id) || coeffs[i].m_coeff.is_int());
+            is_int_row &= is_int(coeffs[i].m_id);
         }
         r.m_alive = true;
         r.m_coeff = c;
         r.m_value = val;
         r.m_type = rel;
+        r.m_mod = m;
+        if (is_int_row && rel == t_lt) {
+            r.m_type = t_le;
+            r.m_coeff  += rational::one();
+            r.m_value  += rational::one();
+        }
         SASSERT(invariant(row_id, r));
     }
 
+    unsigned model_based_opt::new_row() {
+        unsigned row_id = 0;
+        if (m_retired_rows.empty()) {
+            row_id = m_rows.size();
+            m_rows.push_back(row());
+        }
+        else {
+            row_id = m_retired_rows.back();
+            m_retired_rows.pop_back();
+            m_rows[row_id].reset();
+            m_rows[row_id].m_alive = true;
+        }
+        return row_id;
+    }
+
     void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, ineq_type rel) {
-        rational val(c);        
-        unsigned row_id = m_rows.size();
-        m_rows.push_back(row());
-        set_row(row_id, coeffs, c, rel);
+        add_constraint(coeffs, c, rational::zero(), rel);
+    }
+
+    void model_based_opt::add_divides(vector<var> const& coeffs, rational const& c, rational const& m) {
+        add_constraint(coeffs, c, m, t_mod);
+    }
+
+    void model_based_opt::add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel) {
+        unsigned row_id = new_row();
+        set_row(row_id, coeffs, c, m, rel);
         for (unsigned i = 0; i < coeffs.size(); ++i) {
             m_var2row_ids[coeffs[i].m_id].push_back(row_id); 
         }
+        SASSERT(invariant(row_id, m_rows[row_id]));
     }
 
     void model_based_opt::set_objective(vector<var> const& coeffs, rational const& c) {
-        set_row(m_objective_id, coeffs, c, t_le);
+        set_row(m_objective_id, coeffs, c, rational::zero(), t_le);
     }
 
     void model_based_opt::get_live_rows(vector<row>& rows) {
@@ -524,7 +720,8 @@ namespace opt {
     // 
     void model_based_opt::project(unsigned x) {
         unsigned_vector& lub_rows = m_lub;
-        unsigned_vector& glb_rows = m_glb;        
+        unsigned_vector& glb_rows = m_glb;
+        unsigned_vector& mod_rows = m_mod;
         unsigned lub_index = UINT_MAX, glb_index = UINT_MAX;
         bool     lub_strict = false, glb_strict = false;
         rational lub_val, glb_val;
@@ -533,6 +730,8 @@ namespace opt {
         uint_set visited;
         lub_rows.reset();
         glb_rows.reset();
+        mod_rows.reset();
+        bool lub_is_unit = false, glb_is_unit = false;
         // select the lub and glb.
         for (unsigned i = 0; i < row_ids.size(); ++i) {
             unsigned row_id = row_ids[i];
@@ -552,7 +751,10 @@ namespace opt {
                 solve_for(row_id, x);
                 return;
             }
-            if (a.is_pos()) {
+            if (r.m_type == t_mod) {
+                mod_rows.push_back(row_id);
+            }
+            else if (a.is_pos()) {
                 rational lub_value = x_val - (r.m_value/a);
                 if (lub_rows.empty() || 
                     lub_value < lub_val ||
@@ -562,6 +764,7 @@ namespace opt {
                     lub_strict = r.m_type == t_lt;                    
                 }
                 lub_rows.push_back(row_id);
+                lub_is_unit &= a.is_one();
             }
             else {
                 SASSERT(a.is_neg());
@@ -574,36 +777,187 @@ namespace opt {
                     glb_strict = r.m_type == t_lt;                    
                 }
                 glb_rows.push_back(row_id);
+                glb_is_unit &= a.is_minus_one();
             }
         }
-        unsigned row_index = (lub_rows.size() <= glb_rows.size())? lub_index : glb_index;
-        glb_rows.append(lub_rows);            
+
+        if (!mod_rows.empty()) {
+            solve_mod(x, mod_rows);
+            return;
+        }
+        
+        unsigned lub_size = lub_rows.size();
+        unsigned glb_size = glb_rows.size();
+        unsigned row_index = (lub_size <= glb_size) ? lub_index : glb_index;
+        glb_rows.append(lub_rows);  
+        
+        // There are only upper or only lower bounds.
         if (row_index == UINT_MAX) {
             for (unsigned i = 0; i < glb_rows.size(); ++i) {
                 unsigned row_id = glb_rows[i];
                 SASSERT(m_rows[row_id].m_alive);
                 SASSERT(!get_coefficient(row_id, x).is_zero());
-                m_rows[row_id].m_alive = false;
+                retire_row(row_id);
             }
+            return;
         }
-        else {
-            rational coeff = get_coefficient(row_index, x);
-            for (unsigned i = 0; i < glb_rows.size(); ++i) {
-                unsigned row_id = glb_rows[i];
-                if (row_id != row_index) {
-                    resolve(row_index, coeff, row_id, x);
+
+        // The number of matching lower and upper bounds is small.
+        if ((lub_size <= 2 || glb_size <= 2) &&
+            (lub_size <= 3 && glb_size <= 3) && 
+            (!is_int(x) || lub_is_unit || glb_is_unit)) {
+            for (unsigned i = 0; i < lub_size; ++i) {
+                unsigned row_id1 = lub_rows[i];
+                bool last = i + 1 == lub_rows.size();
+                rational coeff = get_coefficient(row_id1, x);
+                for (unsigned j = 0; j < glb_size; ++j) {
+                    unsigned row_id2 = glb_rows[j];
+                    if (last) {
+                        resolve(row_id1, coeff, row_id2, x);
+                    }
+                    else {
+                        row r(m_rows[row_id2]);
+                        m_rows.push_back(r);
+                        resolve(row_id1, coeff, m_rows.size()-1, x);
+                    }
                 }
             }
-            m_rows[row_index].m_alive = false;
+            for (unsigned i = 0; i < lub_size; ++i) {
+                retire_row(lub_rows[i]);
+            }
+            return;
         }
+
+        // General case.
+        rational coeff = get_coefficient(row_index, x);
+        for (unsigned i = 0; i < glb_rows.size(); ++i) {
+            unsigned row_id = glb_rows[i];
+            if (row_id != row_index) {
+                resolve(row_index, coeff, row_id, x);
+            }
+        }
+        retire_row(row_index);
     }
 
+    // 
+    // compute D and u.
+    //
+    // D = lcm(d1, d2)
+    // u = eval(x) mod D
+    // 
+    //   d1 | (a1x + t1) & d2 | (a2x + t2)
+    // = 
+    //   d1 | (a1(D*x' + u) + t1) & d2 | (a2(D*x' + u) + t2)
+    // =
+    //   d1 | (a1*u + t1) & d2 | (a2*u + t2)
+    // 
+    // x := D*x' + u
+    // 
+
+    void model_based_opt::solve_mod(unsigned x, unsigned_vector const& mod_rows) {
+        SASSERT(!mod_rows.empty());
+        rational D(1);
+        for (unsigned i = 0; i < mod_rows.size(); ++i) {
+            D = lcm(D, m_rows[mod_rows[i]].m_mod);            
+        }
+        TRACE("opt", display(tout << "lcm: " << D << " tableau\n"););
+        rational val_x = m_var2value[x];
+        rational u = mod(val_x, D);
+        SASSERT(u.is_nonneg() && u < D);
+        for (unsigned i = 0; i < mod_rows.size(); ++i) {
+            replace_var(mod_rows[i], x, u);            
+            SASSERT(invariant(mod_rows[i], m_rows[mod_rows[i]]));
+        }
+        //
+        // update inequalities such that u is added to t and
+        // D is multiplied to coefficient of x.
+        // the interpretation of the new version of x is (x-u)/D
+        //
+        // a*x + t <= 0
+        // a*(D*x' + u) + t <= 0
+        // a*D*x' + a*u + t <= 0
+        //
+        rational new_val = (val_x - u) / D;
+        SASSERT(new_val.is_int());
+        unsigned y = add_var(new_val, true);
+        unsigned_vector const& row_ids = m_var2row_ids[x];
+        uint_set visited;
+        for (unsigned i = 0; i < row_ids.size(); ++i) {
+            unsigned row_id = row_ids[i];
+            if (!visited.contains(row_id)) {
+                // x |-> D*y + u
+                replace_var(row_id, x, D, y, u);
+                visited.insert(row_id);
+            }
+        }
+        project(y);
+    }
+
+    // update row with: x |-> C
+    void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& C) {
+        row& r = m_rows[row_id];
+        SASSERT(!get_coefficient(row_id, x).is_zero());
+        unsigned sz = r.m_vars.size();
+        unsigned i = 0, j = 0;
+        rational coeff(0);
+        for (; i < sz; ++i) {
+            if (r.m_vars[i].m_id == x) {
+                coeff = r.m_vars[i].m_coeff;
+            }
+            else {
+                if (i != j) {
+                    r.m_vars[j] = r.m_vars[i];
+                }
+                ++j;
+            }
+        }
+        if (j != sz) {
+            r.m_vars.shrink(j);
+        }
+        r.m_coeff += coeff*C;
+        r.m_value += coeff*(C - m_var2value[x]);
+    }
+
+    // update row with: x |-> A*y + B
+    void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B) {
+        row& r = m_rows[row_id];
+        rational coeff = get_coefficient(row_id, x);
+        if (coeff.is_zero()) return;
+        if (!r.m_alive) return;
+        replace_var(row_id, x, B);        
+        r.m_vars.push_back(var(y, coeff*A));
+        r.m_value += coeff*A*m_var2value[y];
+        if (!r.m_vars.empty() && r.m_vars.back().m_id > y) {
+            std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare());
+        }
+        m_var2row_ids[y].push_back(row_id);
+        SASSERT(invariant(row_id, r));
+    }
+
+    // 3x + t = 0 & 7 | (c*x + s) & ax <= u 
+    // 3 | -t  & 21 | (-ct + 3s) & a-t <= 3u
+
+#if 0
+    void solve_for_int(unsigned row_id1, unsigned x) {
+            
+        for (unsigned i = 0; i < num_divs(); ++i) {
+            add_lit(model, lits, mk_divides(c*div_divisor(i), 
+                                            mk_sub(mk_mul(c, div_term(i)), mk_mul(div_coeff(i), t))));
+        }
+    }
+#endif
+
     void model_based_opt::solve_for(unsigned row_id1, unsigned x) {
-        rational a = get_coefficient(row_id1, x);
-        row& r1 = m_rows[row_id1];
+        rational a = get_coefficient(row_id1, x), b;
         SASSERT(!a.is_zero());
-        SASSERT(r1.m_type == t_eq);
-        SASSERT(r1.m_alive);
+        SASSERT(m_rows[row_id1].m_type == t_eq);
+        SASSERT(m_rows[row_id1].m_alive);
+        if (m_var2is_int[x] && !abs(a).is_one()) {
+            row& r1 = m_rows[row_id1];
+            vector<var> coeffs;
+            mk_coeffs_without(coeffs, r1.m_vars, x);
+            add_divides(coeffs, r1.m_coeff, abs(a));
+        }
         unsigned_vector const& row_ids = m_var2row_ids[x];
         uint_set visited;
         visited.insert(row_id1);
@@ -611,15 +965,20 @@ namespace opt {
             unsigned row_id2 = row_ids[i];
             if (!visited.contains(row_id2)) {                
                 visited.insert(row_id2);
-                resolve(row_id1, a, row_id2, x);            
+                row const& r2 = m_rows[row_id2];
+                b = get_coefficient(row_id2, x);
+                if (!b.is_zero()) {
+                    resolve(row_id1, a, row_id2, x);            
+                }
             }
         }
-        r1.m_alive = false;
+        retire_row(row_id1);
     }
 
     void model_based_opt::project(unsigned num_vars, unsigned const* vars) {
         for (unsigned i = 0; i < num_vars; ++i) {
             project(vars[i]);
+            TRACE("opt", display(tout << "After projecting: v" << vars[i] << "\n"););
         }
     }
 

src/math/simplex/model_based_opt.h

@@ -30,7 +30,8 @@ namespace opt {
     enum ineq_type {
         t_eq,
         t_lt,
-        t_le
+        t_le,
+        t_mod
     };
 
 
@@ -52,9 +53,11 @@ namespace opt {
             row(): m_type(t_le), m_value(0), m_alive(false) {}
             vector<var> m_vars;         // variables with coefficients
             rational    m_coeff;        // constant in inequality
+            rational    m_mod;          // value the term divide
             ineq_type   m_type;         // inequality type
             rational    m_value;        // value of m_vars + m_coeff under interpretation of m_var2value.
             bool        m_alive;        // rows can be marked dead if they have been processed.
+            void reset() { m_vars.reset(); m_coeff.reset(); m_value.reset(); }
         };
 
     private:
@@ -63,9 +66,11 @@ namespace opt {
         static const unsigned   m_objective_id = 0;
         vector<unsigned_vector> m_var2row_ids;
         vector<rational>        m_var2value;
+        svector<bool>           m_var2is_int;
         vector<var>             m_new_vars;
-        unsigned_vector         m_lub, m_glb;
+        unsigned_vector         m_lub, m_glb, m_mod;
         unsigned_vector         m_above, m_below;
+        unsigned_vector         m_retired_rows;
         
         bool invariant();
         bool invariant(unsigned index, row const& r);
@@ -82,7 +87,29 @@ namespace opt {
 
         void mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2);
 
-        void set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, ineq_type rel);       
+        void mul_add(unsigned x, rational const& a1, unsigned row_src, rational const& a2, unsigned row_dst);
+
+        void mul(unsigned dst, rational const& c);
+
+        void add(unsigned dst, rational const& c);
+
+        void sub(unsigned dst, rational const& c);
+
+        void set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, rational const& m, ineq_type rel);       
+
+        void add_constraint(vector<var> const& coeffs, rational const& c, rational const& m, ineq_type r);
+
+        void replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B);
+
+        void replace_var(unsigned row_id, unsigned x, rational const& C);
+
+        void normalize(unsigned row_id);
+
+        void mk_coeffs_without(vector<var>& dst, vector<var> const src, unsigned x);
+
+        unsigned new_row();
+
+        rational n_sign(rational const& b) const;
 
         void update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail);
 
@@ -92,12 +119,18 @@ namespace opt {
 
         void solve_for(unsigned row_id, unsigned x);
 
+        void solve_mod(unsigned x, unsigned_vector const& mod_rows);
+
+        bool is_int(unsigned x) const { return m_var2is_int[x]; }
+
+        void retire_row(unsigned row_id);
+
     public:
 
         model_based_opt();
 
         // add a fresh variable with value 'value'.
-        unsigned add_var(rational const& value);
+        unsigned add_var(rational const& value, bool is_int = false);
 
         // retrieve updated value of variable.
         rational get_value(unsigned var_id);
@@ -106,6 +139,9 @@ namespace opt {
         // satisfied under the values provided to the variables.
         void add_constraint(vector<var> const& coeffs, rational const& c, ineq_type r);
 
+        // add a divisibility constraint. The row should divide m.
+        void add_divides(vector<var> const& coeffs, rational const& c, rational const& m);
+
         // Set the objective function (linear).
         void set_objective(vector<var> const& coeffs, rational const& c);
 

src/qe/qe_mbp.cpp

@@ -270,6 +270,8 @@ class mbp::impl {
                 sub(fmls[i].get(), val);
                 m_rw(val);
                 if (!m.is_true(val)) {
+                    TRACE("qe", tout << mk_pp(fmls[i].get(), m) << " -> " << val << "\n";);
+                    fmls[i] = val;
                     if (j != i) {
                         fmls[j] = fmls[i].get();
                     }

src/test/model_based_opt.cpp

@@ -287,9 +287,76 @@ static void test7() {
     mbo.display(std::cout);
 }
 
+static void test8() {
+    opt::model_based_opt mbo;
+    unsigned x0 = mbo.add_var(rational(2));
+    unsigned x = mbo.add_var(rational(1));
+    unsigned y = mbo.add_var(rational(3));
+    unsigned z = mbo.add_var(rational(4));
+    unsigned u = mbo.add_var(rational(5));
+    unsigned v = mbo.add_var(rational(6));
+    unsigned w = mbo.add_var(rational(6));
+
+    add_ineq(mbo, x0, 1, y, -1, 0, opt::t_le);
+    add_ineq(mbo, x, 1, y, -1, 0, opt::t_lt);
+    add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
+    add_ineq(mbo, y, 1, z, -1, 1, opt::t_le);
+    add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
+
+    mbo.display(std::cout);
+    mbo.project(1, &y);
+    mbo.display(std::cout);
+}
+
+
+static void test9() {
+    opt::model_based_opt mbo;
+    unsigned x0 = mbo.add_var(rational(2), true);
+    unsigned x = mbo.add_var(rational(1), true);
+    unsigned y = mbo.add_var(rational(3), true);
+    unsigned z = mbo.add_var(rational(4), true);
+    unsigned u = mbo.add_var(rational(5), true);
+    unsigned v = mbo.add_var(rational(6), true);
+
+    add_ineq(mbo, x0, 1, y, -1, 0, opt::t_le);
+    add_ineq(mbo, x, 1, y, -1, 0, opt::t_lt);
+    add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
+    add_ineq(mbo, y, 1, z, -1, 1, opt::t_le);
+    add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
+
+    mbo.display(std::cout);
+    mbo.project(1, &y);
+    mbo.display(std::cout);
+}
+
+
+static void test10() {
+    opt::model_based_opt mbo;
+    unsigned x0 = mbo.add_var(rational(2), true);
+    unsigned x = mbo.add_var(rational(1), true);
+    unsigned y = mbo.add_var(rational(3), true);
+    unsigned z = mbo.add_var(rational(4), true);
+    unsigned u = mbo.add_var(rational(5), true);
+    unsigned v = mbo.add_var(rational(6), true);
+
+    add_ineq(mbo, x0, 1, y, -2, 0, opt::t_le);
+    add_ineq(mbo, x, 1, y, -2, 0, opt::t_lt);
+    add_ineq(mbo, y, 3, u, -4, 0, opt::t_le);
+    add_ineq(mbo, y, 3, z, -5, 1, opt::t_le);
+    add_ineq(mbo, y, 3, v, -6, 1, opt::t_le);
+
+    mbo.display(std::cout);
+    mbo.project(1, &y);
+    mbo.display(std::cout);
+    mbo.project(1, &x0);
+    mbo.display(std::cout);
+}
+
 // test with mix of upper and lower bounds
 
 void tst_model_based_opt() {
+    test10();
+    return;
     check_random_ineqs();
     test1();
     test2();
@@ -298,5 +365,6 @@ void tst_model_based_opt() {
     test5();
     test6();
     test7();
-
+    test8();
+    test9();
 }

src/test/qe_arith.cpp

@@ -438,10 +438,158 @@ static void check_random_ineqs() {
     }
 }
 
+static void test_project() {
+    ast_manager m;
+    reg_decl_plugins(m);    
+    qe::arith_project_plugin plugin(m);    
+    arith_util a(m);
+    app_ref_vector vars(m);
+    expr_ref_vector lits(m), ds(m);
+    model mdl(m);    
+    app_ref x(m), y(m), z(m), u(m);
+    x = m.mk_const(symbol("x"), a.mk_int());
+    y = m.mk_const(symbol("y"), a.mk_int());
+    z = m.mk_const(symbol("z"), a.mk_int());
+    u = m.mk_const(symbol("u"), a.mk_int());
+    func_decl_ref f(m);
+    sort* int_sort = a.mk_int();
+    f = m.mk_func_decl(symbol("f"), 1, &int_sort, int_sort);
 
+    // test non-projection
+    mdl.register_decl(x->get_decl(), a.mk_int(0));
+    mdl.register_decl(y->get_decl(), a.mk_int(1));
+    mdl.register_decl(z->get_decl(), a.mk_int(2));
+    mdl.register_decl(u->get_decl(), a.mk_int(3));
+    func_interp* fi = alloc(func_interp, m, 1);
+    expr_ref_vector nums(m);
+    nums.push_back(a.mk_int(0));
+    nums.push_back(a.mk_int(1));
+    nums.push_back(a.mk_int(2));
+    fi->insert_new_entry(nums.c_ptr(),   a.mk_int(1));
+    fi->insert_new_entry(nums.c_ptr()+1, a.mk_int(2));
+    fi->insert_new_entry(nums.c_ptr()+2, a.mk_int(3));
+    fi->set_else(a.mk_int(10));
+    mdl.register_decl(f, fi);
+    vars.reset();
+    lits.reset();
+    vars.push_back(x);
+    lits.push_back(x <= app_ref(m.mk_app(f, (expr*)x), m));
+    lits.push_back(x < y);
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+    // test not-equals
+    vars.reset();
+    lits.reset();
+    vars.push_back(x);
+    lits.push_back(m.mk_not(m.mk_eq(x, y)));
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+    // test negation of distinct using bound variables
+    mdl.register_decl(x->get_decl(), a.mk_int(0));
+    mdl.register_decl(y->get_decl(), a.mk_int(1));
+    mdl.register_decl(z->get_decl(), a.mk_int(0));
+    mdl.register_decl(u->get_decl(), a.mk_int(6));
+    vars.reset();
+    lits.reset();
+    ds.reset();
+    vars.push_back(x);
+    vars.push_back(y);
+    ds.push_back(x);
+    ds.push_back(y);
+    ds.push_back(z + 2);
+    ds.push_back(u);
+    ds.push_back(z);
+    lits.push_back(m.mk_not(m.mk_distinct(ds.size(), ds.c_ptr())));
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+    // test negation of distinct, not using bound variables
+    mdl.register_decl(x->get_decl(), a.mk_int(0));
+    mdl.register_decl(y->get_decl(), a.mk_int(1));
+    mdl.register_decl(z->get_decl(), a.mk_int(0));
+    mdl.register_decl(u->get_decl(), a.mk_int(6));
+    vars.reset();
+    lits.reset();
+    ds.reset();
+    vars.push_back(x);
+    vars.push_back(y);
+    ds.push_back(x);
+    ds.push_back(y);
+    ds.push_back(z + 2);
+    ds.push_back(u);
+    ds.push_back(z + 10);
+    ds.push_back(u + 4);
+    lits.push_back(m.mk_not(m.mk_distinct(ds.size(), ds.c_ptr())));
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+
+    // test distinct
+    mdl.register_decl(x->get_decl(), a.mk_int(0));
+    mdl.register_decl(y->get_decl(), a.mk_int(1));
+    mdl.register_decl(z->get_decl(), a.mk_int(0));
+    mdl.register_decl(u->get_decl(), a.mk_int(6));
+    vars.reset();
+    lits.reset();
+    ds.reset();
+    vars.push_back(x);
+    vars.push_back(y);
+    ds.push_back(x);
+    ds.push_back(y);
+    ds.push_back(z + 2);
+    ds.push_back(u);
+    lits.push_back(m.mk_distinct(ds.size(), ds.c_ptr()));
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+    // equality over modulus
+    mdl.register_decl(y->get_decl(), a.mk_int(4));
+    mdl.register_decl(z->get_decl(), a.mk_int(8));
+    lits.reset();
+    vars.reset();
+    vars.push_back(y);
+    lits.push_back(m.mk_eq(a.mk_mod(y, a.mk_int(3)), a.mk_int(1)));
+    lits.push_back(m.mk_eq(2*y, z));
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+    // inequality test
+    mdl.register_decl(x->get_decl(), a.mk_int(0));
+    mdl.register_decl(y->get_decl(), a.mk_int(1));
+    mdl.register_decl(z->get_decl(), a.mk_int(0));
+    mdl.register_decl(u->get_decl(), a.mk_int(6));
+    vars.reset();
+    lits.reset();
+    vars.push_back(x);
+    vars.push_back(y);
+    lits.push_back(z <= (x + (2*y)));
+    lits.push_back(2*x < u + 3);
+    lits.push_back(2*y <= u);
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+    // non-unit equalities
+    mdl.register_decl(x->get_decl(), a.mk_int(1));
+    mdl.register_decl(z->get_decl(), a.mk_int(2));
+    mdl.register_decl(u->get_decl(), a.mk_int(3));
+    mdl.register_decl(y->get_decl(), a.mk_int(4));
+    lits.reset();
+    vars.reset();
+    vars.push_back(x);
+    lits.push_back(m.mk_eq(2*x, z));
+    lits.push_back(m.mk_eq(3*x, u));
+    plugin(mdl, vars, lits);
+    std::cout << lits << "\n";
+
+
+}
 
 
 void tst_qe_arith() {
+    test_project();
+    return;
     check_random_ineqs();
     return;
 //    enable_trace("qe");