Dev (#56)

* introduce int_solver.h

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* add int_solver class

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* track which var is an integer

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* add queries for integrality of vars

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* resurrect lp_tst in its own director lp

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* add file

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* add_constraint has got a body

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* fix add_constraint and substitute_terms_in_linear_expression

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* after merge with Z3Prover

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* adding stub check_int_feasibility()

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* Dev (#50)

* initial skeletons for nra solver

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

* initial skeletons for nra solver

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

* small fix in lar_solver.cpp

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* adding some content to the new check_int_feasibility()

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* Dev (#51)

* initial skeletons for nra solver

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

* initial skeletons for nra solver

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

* adding more nlsat

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

* nlsat integration

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

* adding constraints

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

* adding nra solver

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

* add missing initialization

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

* adding nra solver

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

* test

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* Dev (#53)

* change in a comment

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* Disabled debug output

* removing FOCI2 interface from interp

* remove foci reference from cmakelist.txt

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

* initial skeletons for nra solver

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

* initial skeletons for nra solver

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

* adding more nlsat

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

* nlsat integration

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

* adding constraints

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

* adding nra solver

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

* add missing initialization

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

* adding nra solver

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

* adding nra

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

* debugging nra

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

* updates to nra_solver integration to call it directly from theory_lra instead of over lar_solver

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

* n/a

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

* integrate nlsat

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

* tidy

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

* preserve is_int flag

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* remove a debug printout

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* Dev (#54)

* change in a comment

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* Disabled debug output

* removing FOCI2 interface from interp

* remove foci reference from cmakelist.txt

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

* initial skeletons for nra solver

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

* initial skeletons for nra solver

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

* adding more nlsat

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

* nlsat integration

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

* adding constraints

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

* adding nra solver

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

* add missing initialization

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

* adding nra solver

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

* adding nra

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

* debugging nra

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

* updates to nra_solver integration to call it directly from theory_lra instead of over lar_solver

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

* n/a

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

* integrate nlsat

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

* tidy

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

* use integer test from lra solver, updated it to work on term variables

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

* fix equality check in assume-eq

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

* fix model_is_int_feasible

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* untested gcd_test()

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* call fill_explanation_from_fixed_columns()

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* add the call to pivot_fixed_vars_from_basis() to int_solver.cpp::check()

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* port more of theory_arith_int.h

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* use statistics of lar_solver by theory_lra.cpp

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* port more code to int_solver.cpp

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* add an assert

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* more int porting

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* fix a bug in pivot_fixed_vars_from_basis

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* small change

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

* implement find_inf_int_base_column()

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* catch unregistered vars in add_var_bound

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* add a file

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* compile for vs2012

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* fix asserts in add_var_bound

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* fix the lp_solver init when workig on an mps file

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* towards int_solver::check()

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* change in int_solver::check() signature

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

* add handlers for lia moves

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

* spacing

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

src/ast/ast.cpp

@@ -1714,7 +1714,6 @@ ast * ast_manager::register_node_core(ast * n) {
 
 
     n->m_id   = is_decl(n) ? m_decl_id_gen.mk() : m_expr_id_gen.mk();
-
     TRACE("ast", tout << "Object " << n->m_id << " was created.\n";);
     TRACE("mk_var_bug", tout << "mk_ast: " << n->m_id << "\n";);
     // increment reference counters

src/ast/rewriter/arith_rewriter.cpp

@@ -800,105 +800,56 @@ br_status arith_rewriter::mk_idiv_core(expr * arg1, expr * arg2, expr_ref & resu
         result = m_util.mk_numeral(div(v1, v2), is_int);
         return BR_DONE;
     }
-    if (m_util.is_numeral(arg2, v2, is_int) && v2.is_one()) {
-        result = arg1;
-        return BR_DONE;
-    }
     if (m_util.is_numeral(arg2, v2, is_int) && v2.is_zero()) {
-        return BR_FAILED;
+        result = m_util.mk_idiv0(arg1);
+        return BR_REWRITE1;
     }
-    if (arg1 == arg2) {
-        expr_ref zero(m_util.mk_int(0), m());
-        result = m().mk_ite(m().mk_eq(arg1, zero), m_util.mk_idiv(zero, zero), m_util.mk_int(1));
-        return BR_REWRITE3;
-    }
-    if (divides(arg1, arg2, result)) {
-        return BR_REWRITE_FULL;
+    expr_ref quot(m());
+    if (divides(arg1, arg2, quot)) {
+        result = m_util.mk_mul(quot, m_util.mk_idiv(arg1, arg1));
+        return BR_REWRITE2;
     }
     return BR_FAILED;
 }
 
-
-// 
-// implement div ab ac = floor( ab / ac) = floor (b / c) = div b c
-//
-bool arith_rewriter::divides(expr* num, expr* den, expr_ref& result) {
-    expr_fast_mark1 mark;
-    rational num_r(1), den_r(1);
-    expr* num_e = nullptr, *den_e = nullptr;
-    ptr_buffer<expr> args1, args2;
-    flat_mul(num, args1);
-    flat_mul(den, args2);
-    for (expr * arg : args1) {
-        mark.mark(arg);
-        if (m_util.is_numeral(arg, num_r)) num_e = arg;
-    }
-    for (expr* arg : args2) {
-        if (mark.is_marked(arg)) {
-            result = remove_divisor(arg, num, den);
-            return true;
-        }
-        if (m_util.is_numeral(arg, den_r)) den_e = arg;
-    }
-    rational g = gcd(num_r, den_r);
-    if (!g.is_one()) {
-        SASSERT(g.is_pos());
-        // replace num_e, den_e by their gcd reduction.
-        for (unsigned i = 0; i < args1.size(); ++i) {
-            if (args1[i] == num_e) {
-                args1[i] = m_util.mk_numeral(num_r / g, true);
-                break;
-            }
-        }
-        for (unsigned i = 0; i < args2.size(); ++i) {
-            if (args2[i] == den_e) {
-                args2[i] = m_util.mk_numeral(den_r / g, true);
-                break;
-            }
-        }
-        
-        num = m_util.mk_mul(args1.size(), args1.c_ptr());
-        den = m_util.mk_mul(args2.size(), args2.c_ptr());
-        result = m_util.mk_idiv(num, den);
+bool arith_rewriter::divides(expr* d, expr* n, expr_ref& quot) {
+    if (d == n) {
+        quot = m_util.mk_numeral(rational(1), m_util.is_int(d));
         return true;
     }
-    return false;
-}
-
-expr_ref arith_rewriter::remove_divisor(expr* arg, expr* num, expr* den) {
-    ptr_buffer<expr> args1, args2;
-    flat_mul(num, args1);
-    flat_mul(den, args2);
-    remove_divisor(arg, args1);
-    remove_divisor(arg, args2);
-    expr_ref zero(m_util.mk_int(0), m());
-    num = args1.empty() ? m_util.mk_int(1) : m_util.mk_mul(args1.size(), args1.c_ptr());
-    den = args2.empty() ? m_util.mk_int(1) : m_util.mk_mul(args2.size(), args2.c_ptr());
-    return expr_ref(m().mk_ite(m().mk_eq(zero, arg), m_util.mk_idiv(zero, zero), m_util.mk_idiv(num, den)), m());
-}
-
-void arith_rewriter::flat_mul(expr* e, ptr_buffer<expr>& args) {
-    args.push_back(e);
-    for (unsigned i = 0; i < args.size(); ++i) {
-        e = args[i];
-        if (m_util.is_mul(e)) {
-            args.append(to_app(e)->get_num_args(), to_app(e)->get_args());
-            args[i] = args.back();
-            args.shrink(args.size()-1);
-            --i;
-        }                
-    }
-}
-
-void arith_rewriter::remove_divisor(expr* d, ptr_buffer<expr>& args) {
-    for (unsigned i = 0; i < args.size(); ++i) {
-        if (args[i] == d) {
-            args[i] = args.back();
-            args.shrink(args.size()-1);
-            return;
+    if (m_util.is_mul(n)) {
+        expr_ref_vector muls(m());
+        muls.push_back(n);
+        expr* n1, *n2;
+        rational r1, r2;
+        for (unsigned i = 0; i < muls.size(); ++i) {
+            if (m_util.is_mul(muls[i].get(), n1, n2)) {
+                muls[i] = n1;
+                muls.push_back(n2);
+                --i;
+            }
+        }
+        if (m_util.is_numeral(d, r1) && !r1.is_zero()) {
+            for (unsigned i = 0; i < muls.size(); ++i) {
+                if (m_util.is_numeral(muls[i].get(), r2) && (r2 / r1).is_int()) {
+                    muls[i] = m_util.mk_numeral(r2 / r1, m_util.is_int(d));
+                    quot = m_util.mk_mul(muls.size(), muls.c_ptr());
+                    return true;
+                }            
+            }
+        }
+        else {
+            for (unsigned i = 0; i < muls.size(); ++i) {
+                if (d == muls[i].get()) {
+                    muls[i] = muls.back();
+                    muls.pop_back();
+                    quot = m_util.mk_mul(muls.size(), muls.c_ptr());
+                    return true;
+                }            
+            }
         }
     }
-    UNREACHABLE();
+    return false;
 }
 
 br_status arith_rewriter::mk_mod_core(expr * arg1, expr * arg2, expr_ref & result) {

src/ast/rewriter/arith_rewriter.h

@@ -95,10 +95,7 @@ class arith_rewriter : public poly_rewriter<arith_rewriter_core> {
     expr_ref neg_monomial(expr * e) const;
     expr * mk_sin_value(rational const & k);
     app * mk_sqrt(rational const & k);
-    bool divides(expr* d, expr* n, expr_ref& result);
-    expr_ref remove_divisor(expr* arg, expr* num, expr* den);
-    void flat_mul(expr* e, ptr_buffer<expr>& args);
-    void remove_divisor(expr* d, ptr_buffer<expr>& args);
+    bool divides(expr* d, expr* n, expr_ref& quot);
 
 public:
     arith_rewriter(ast_manager & m, params_ref const & p = params_ref()):

src/ast/simplifier/arith_simplifier_plugin.h

@@ -0,0 +1,97 @@
+/*++
+Copyright (c) 2007 Microsoft Corporation
+
+Module Name:
+
+    arith_simplifier_plugin.h
+
+Abstract:
+
+    Simplifier for the arithmetic family.
+
+Author:
+
+    Leonardo (leonardo) 2008-01-08
+    
+--*/
+#ifndef ARITH_SIMPLIFIER_PLUGIN_H_
+#define ARITH_SIMPLIFIER_PLUGIN_H_
+
+#include"basic_simplifier_plugin.h"
+#include"poly_simplifier_plugin.h"
+#include"arith_decl_plugin.h"
+#include"arith_simplifier_params.h"
+
+/**
+   \brief Simplifier for the arith family.
+*/
+class arith_simplifier_plugin : public poly_simplifier_plugin {
+public:
+    enum op_kind {
+        LE, GE, EQ
+    };
+protected:
+    arith_simplifier_params & m_params;
+    arith_util                m_util;
+    basic_simplifier_plugin & m_bsimp;
+    expr_ref                  m_int_zero;
+    expr_ref                  m_real_zero;
+
+    bool is_neg_poly(expr * t) const;
+
+    template<op_kind k>
+    void mk_le_ge_eq_core(expr * arg1, expr * arg2, expr_ref & result);
+
+    void prop_mod_const(expr * e, unsigned depth, numeral const& k, expr_ref& result);
+
+    void gcd_reduce_monomial(expr_ref_vector& monomials, numeral& k);
+
+    void div_monomial(expr_ref_vector& monomials, numeral const& g);
+    void get_monomial_gcd(expr_ref_vector& monomials, numeral& g);
+    bool divides(expr* d, expr* n, expr_ref& quot);
+
+public:
+    arith_simplifier_plugin(ast_manager & m, basic_simplifier_plugin & b, arith_simplifier_params & p);
+    ~arith_simplifier_plugin();
+    arith_util & get_arith_util() { return m_util; }
+    virtual numeral norm(const numeral & n) { return n; }
+    virtual bool is_numeral(expr * n, rational & val) const { bool f; return m_util.is_numeral(n, val, f); }
+    bool is_numeral(expr * n) const { return m_util.is_numeral(n); }
+    virtual bool is_minus_one(expr * n) const { numeral tmp; return is_numeral(n, tmp) && tmp.is_minus_one(); }
+    virtual expr * get_zero(sort * s) const { return m_util.is_int(s) ? m_int_zero.get() : m_real_zero.get(); }
+
+    virtual app * mk_numeral(numeral const & n) { return m_util.mk_numeral(n, m_curr_sort->get_decl_kind() == INT_SORT); }
+    app * mk_numeral(numeral const & n, bool is_int) { return m_util.mk_numeral(n, is_int); }
+    bool is_int_sort(sort const * s) const { return m_util.is_int(s); }
+    bool is_real_sort(sort const * s) const { return m_util.is_real(s); }
+    bool is_arith_sort(sort const * s) const { return is_int_sort(s) || is_real_sort(s); }
+    bool is_int(expr const * n) const { return m_util.is_int(n); }
+    bool is_le(expr const * n) const { return m_util.is_le(n); }
+    bool is_ge(expr const * n) const { return m_util.is_ge(n); }
+
+    virtual bool is_le_ge(expr * n) const { return is_le(n) || is_ge(n); }
+    
+    void mk_le(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_ge(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_lt(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_gt(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_arith_eq(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_div(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_idiv(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_mod(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_rem(expr * arg1, expr * arg2, expr_ref & result);
+    void mk_to_real(expr * arg, expr_ref & result);
+    void mk_to_int(expr * arg, expr_ref & result);
+    void mk_is_int(expr * arg, expr_ref & result);
+    void mk_abs(expr * arg, expr_ref & result);
+
+    virtual bool reduce(func_decl * f, unsigned num_args, expr * const * args, expr_ref & result);
+    virtual bool reduce_eq(expr * lhs, expr * rhs, expr_ref & result);
+
+    bool is_arith_term(expr * n) const;
+
+    void gcd_normalize(numeral & coeff, expr_ref& term);
+
+};
+
+#endif /* ARITH_SIMPLIFIER_PLUGIN_H_ */

src/math/polynomial/polynomial.h

@@ -19,17 +19,16 @@ Notes:
 #ifndef POLYNOMIAL_H_
 #define POLYNOMIAL_H_
 
-#include "util/mpz.h"
-#include "util/rational.h"
-#include "util/obj_ref.h"
-#include "util/ref_vector.h"
-#include "util/z3_exception.h"
-#include "util/scoped_numeral.h"
-#include "util/scoped_numeral_vector.h"
-#include "util/params.h"
-#include "util/mpbqi.h"
-#include "util/rlimit.h"
-#include "util/lbool.h"
+#include"util/mpz.h"
+#include"util/rational.h"
+#include"util/obj_ref.h"
+#include"util/ref_vector.h"
+#include"util/z3_exception.h"
+#include"util/scoped_numeral.h"
+#include"util/scoped_numeral_vector.h"
+#include"util/params.h"
+#include"util/mpbqi.h"
+#include"util/rlimit.h"
 
 class small_object_allocator;
 

src/nlsat/nlsat_solver.h

@@ -21,10 +21,10 @@ Revision History:
 #ifndef NLSAT_SOLVER_H_
 #define NLSAT_SOLVER_H_
 
-#include "nlsat/nlsat_types.h"
-#include "util/params.h"
-#include "util/statistics.h"
-#include "util/rlimit.h"
+#include"nlsat/nlsat_types.h"
+#include"util/params.h"
+#include"util/statistics.h"
+#include"util/rlimit.h"
 
 namespace nlsat {
 

src/nlsat/nlsat_types.h

@@ -19,10 +19,10 @@ Revision History:
 #ifndef NLSAT_TYPES_H_
 #define NLSAT_TYPES_H_
 
-#include "math/polynomial/polynomial.h"
-#include "util/buffer.h"
-#include "sat/sat_types.h"
-#include "util/z3_exception.h"
+#include"math/polynomial/polynomial.h"
+#include"util/buffer.h"
+#include"sat/sat_types.h"
+#include"util/z3_exception.h"
 
 namespace algebraic_numbers {
     class anum;

src/sat/sat_types.h

@@ -19,13 +19,12 @@ Revision History:
 #ifndef SAT_TYPES_H_
 #define SAT_TYPES_H_
 
-#include "util/debug.h"
-#include "util/approx_set.h"
-#include "util/lbool.h"
-#include "util/z3_exception.h"
-#include "util/common_msgs.h"
-#include "util/vector.h"
-#include "util/uint_set.h"
+#include"util/debug.h"
+#include"util/approx_set.h"
+#include"util/lbool.h"
+#include"util/z3_exception.h"
+#include"util/common_msgs.h"
+#include"util/vector.h"
 #include<iomanip>
 
 namespace sat {

src/shell/lp_frontend.cpp

@@ -81,7 +81,6 @@ void run_solver(lp_params & params, char const * mps_file_name) {
     solver->settings().report_frequency = params.rep_freq();
     solver->settings().print_statistics = params.print_stats();
     solver->settings().simplex_strategy() = lp:: simplex_strategy_enum::lu;
-
     solver->find_maximal_solution();
 
     *(solver->settings().get_message_ostream()) << "status is " << lp_status_to_string(solver->get_status()) << std::endl;

src/smt/theory_lra.cpp

@@ -36,7 +36,10 @@ Revision History:
 #include "smt/smt_model_generator.h"
 #include "smt/arith_eq_adapter.h"
 #include "util/nat_set.h"
-#include "tactic/generic_model_converter.h"
+#include "util/lp/nra_solver.h"
+#include "tactic/filter_model_converter.h"
+#include "math/polynomial/algebraic_numbers.h"
+#include "math/polynomial/polynomial.h"
 
 namespace lra_lp {
     enum bound_kind { lower_t, upper_t };
@@ -88,16 +91,12 @@ namespace lra_lp {
         unsigned m_bounds_propagations;
         unsigned m_num_iterations;
         unsigned m_num_iterations_with_no_progress;
-        unsigned m_num_factorizations;
         unsigned m_need_to_solve_inf;
         unsigned m_fixed_eqs;
         unsigned m_conflicts;
         unsigned m_bound_propagations1;
         unsigned m_bound_propagations2;
         unsigned m_assert_diseq;
-        unsigned m_make_feasible;
-        unsigned m_max_cols;
-        unsigned m_max_rows;
         stats() { reset(); }
         void reset() {
             memset(this, 0, sizeof(*this));
@@ -119,9 +118,6 @@ namespace smt {
             unsigned m_bounds_lim;
             unsigned m_asserted_qhead;
             unsigned m_asserted_atoms_lim;
-            unsigned m_delayed_terms_lim;
-            unsigned m_delayed_equalities_lim;
-            unsigned m_delayed_defs_lim;
             unsigned m_underspecified_lim;
             unsigned m_var_trail_lim;
             expr*    m_not_handled;
@@ -145,10 +141,10 @@ namespace smt {
         ast_manager&         m;
         theory_arith_params& m_arith_params;
         arith_util           a;
-
         arith_eq_adapter     m_arith_eq_adapter;
+        vector<rational>     m_columns;
+      
 
-        vector<rational>    m_columns;
         // temporary values kept during internalization
         struct internalize_state {
             expr_ref_vector     m_terms;
@@ -197,16 +193,8 @@ namespace smt {
                 coeffs().pop_back();
             }
         };
-
-        typedef vector<std::pair<rational, lp::var_index>> var_coeffs;
-        struct delayed_def {
-            vector<rational>    m_coeffs;
-            svector<theory_var> m_vars;
-            rational            m_coeff;
-            theory_var          m_var;
-            delayed_def(svector<theory_var> const& vars, vector<rational> const& coeffs, rational const& r, theory_var v):
-                m_coeffs(coeffs), m_vars(vars), m_coeff(r), m_var(v) {}
-        };
+       
+        typedef vector<std::pair<rational, lean::var_index>> var_coeffs;
 
         svector<lp::var_index> m_theory_var2var_index;   // translate from theory variables to lar vars
         svector<theory_var>      m_var_index2theory_var;   // reverse map from lp_solver variables to theory variables
@@ -225,11 +213,7 @@ namespace smt {
         svector<enode_pair>                           m_equalities;      // asserted rows corresponding to equalities.
         svector<theory_var>                           m_definitions;     // asserted rows corresponding to definitions
 
-        bool                   m_delay_constraints;    // configuration
-        svector<delayed_atom>  m_asserted_atoms;
-        app_ref_vector         m_delayed_terms;
-        svector<std::pair<theory_var, theory_var>> m_delayed_equalities;
-        vector<delayed_def>    m_delayed_defs;
+        svector<delayed_atom>  m_asserted_atoms;        
         expr*                  m_not_handled;
         ptr_vector<app>        m_underspecified;
         unsigned_vector        m_var_trail;
@@ -249,16 +233,36 @@ namespace smt {
 
         unsigned               m_num_conflicts;
 
+        // non-linear arithmetic
+        scoped_ptr<nra::solver> m_nra;
+        bool                    m_use_nra_model;
+        scoped_ptr<scoped_anum> m_a1, m_a2;
+
+        // integer arithmetic
+        scoped_ptr<lean::int_solver> m_lia;
+
 
         struct var_value_eq {
             imp & m_th;
             var_value_eq(imp & th):m_th(th) {}
-            bool operator()(theory_var v1, theory_var v2) const { return m_th.get_ivalue(v1) == m_th.get_ivalue(v2) && m_th.is_int(v1) == m_th.is_int(v2); }
+            bool operator()(theory_var v1, theory_var v2) const { 
+                if (m_th.is_int(v1) != m_th.is_int(v2)) {
+                    return false;
+                }
+                return m_th.is_eq(v1, v2);
+            }
         };
         struct var_value_hash {
             imp & m_th;
             var_value_hash(imp & th):m_th(th) {}
-            unsigned operator()(theory_var v) const { return (unsigned)std::hash<lp::impq>()(m_th.get_ivalue(v)); }
+            unsigned operator()(theory_var v) const { 
+                if (m_th.m_use_nra_model) {
+                    return m_th.is_int(v);
+                }
+                else {
+                    return (unsigned)std::hash<lean::impq>()(m_th.get_ivalue(v)); 
+                }
+            }
         };
         int_hashtable<var_value_hash, var_value_eq>   m_model_eqs;
 
@@ -290,14 +294,25 @@ namespace smt {
             m_solver->settings().bound_propagation() = BP_NONE != propagation_mode();
             m_solver->set_propagate_bounds_on_pivoted_rows_mode(lp.bprop_on_pivoted_rows());
             //m_solver->settings().set_ostream(0);
+            m_lia = alloc(lean::int_solver, m_solver.get());
         }
 
+        void ensure_nra() {
+            if (!m_nra) {
+                m_nra = alloc(nra::solver, *m_solver.get(), m.limit(), ctx().get_params());
+                for (unsigned i = 0; i < m_scopes.size(); ++i) {
+                    m_nra->push();
+                }
+            }
+        }
+
+
         void found_not_handled(expr* n) {
             m_not_handled = n;
             if (is_app(n) && is_underspecified(to_app(n))) {
+                TRACE("arith", tout << "Unhandled: " << mk_pp(n, m) << "\n";);
                 m_underspecified.push_back(to_app(n));
             }
-            TRACE("arith", tout << "Unhandled: " << mk_pp(n, m) << "\n";);
         }
 
         bool is_numeral(expr* term, rational& r) {
@@ -367,6 +382,14 @@ namespace smt {
                     terms[index] = n1;
                     st.terms_to_internalize().push_back(n2);
                 }
+                else if (a.is_mul(n)) {
+                    theory_var v;
+                    internalize_mul(to_app(n), v, r);
+                    coeffs[index] *= r;
+                    coeffs[vars.size()] = coeffs[index];
+                    vars.push_back(v);
+                    ++index;
+                }
                 else if (a.is_numeral(n, r)) {
                     coeff += coeffs[index]*r;
                     ++index;
@@ -421,6 +444,44 @@ namespace smt {
             }
         }
 
+        void internalize_mul(app* t, theory_var& v, rational& r) {
+            SASSERT(a.is_mul(t));
+            bool _has_var = has_var(t);
+            if (!_has_var) {
+                internalize_args(t);
+                mk_enode(t);
+            }
+            r = rational::one();
+            rational r1;
+            v = mk_var(t);
+            svector<lean::var_index> vars;
+            ptr_vector<expr> todo;
+            todo.push_back(t);
+            while (!todo.empty()) {
+                expr* n = todo.back();
+                todo.pop_back();
+                expr* n1, *n2;
+                if (a.is_mul(n, n1, n2)) {
+                    todo.push_back(n1);
+                    todo.push_back(n2);
+                }
+                else if (a.is_numeral(n, r1)) {
+                    r *= r1;
+                }
+                else {
+                    if (!ctx().e_internalized(n)) {
+                        ctx().internalize(n, false);
+                    }
+                    vars.push_back(get_var_index(mk_var(n)));
+                }
+            }
+            TRACE("arith", tout << mk_pp(t, m) << "\n";);
+            if (!_has_var) {
+                ensure_nra();
+                m_nra->add_monomial(get_var_index(v), vars.size(), vars.c_ptr());
+            }
+        }
+
         enode * mk_enode(app * n) {
             if (ctx().e_internalized(n)) {
                 return get_enode(n);
@@ -465,6 +526,14 @@ namespace smt {
             return m_arith_params.m_arith_reflect || is_underspecified(n);
         }
 
+        bool has_var(expr* n) {
+            if (!ctx().e_internalized(n)) {
+                return false;
+            }
+            enode* e = get_enode(n);
+            return th.is_attached_to_var(e);
+        }
+
         theory_var mk_var(expr* n, bool internalize = true) {
             if (!ctx().e_internalized(n)) {
                 ctx().internalize(n, false);
@@ -493,7 +562,7 @@ namespace smt {
                 result = m_theory_var2var_index[v];
             }
             if (result == UINT_MAX) {
-                result = m_solver->add_var(v); // TBD: is_int(v);
+                result = m_solver->add_var(v, is_int(v));
                 m_theory_var2var_index.setx(v, result, UINT_MAX);
                 m_var_index2theory_var.setx(result, v, UINT_MAX);
                 m_var_trail.push_back(v);
@@ -555,16 +624,8 @@ namespace smt {
             m_definitions.setx(index, v, null_theory_var);
             ++m_stats.m_add_rows;
         }
-
-        void internalize_eq(delayed_def const& d) {
-            scoped_internalize_state st(*this);
-            st.vars().append(d.m_vars);
-            st.coeffs().append(d.m_coeffs);
-            init_left_side(st);
-            add_def_constraint(m_solver->add_constraint(m_left_side, lp::EQ, -d.m_coeff), d.m_var);
-        }
-
-        void internalize_eq(theory_var v1, theory_var v2) {
+        
+        void internalize_eq(theory_var v1, theory_var v2) {                  
             enode* n1 = get_enode(v1);
             enode* n2 = get_enode(v2);
             scoped_internalize_state st(*this);
@@ -656,15 +717,14 @@ namespace smt {
             a(m),
             m_arith_eq_adapter(th, ap, a),
             m_internalize_head(0),
-            m_delay_constraints(false),
-            m_delayed_terms(m),
-            m_not_handled(nullptr),
-            m_asserted_qhead(0),
+            m_not_handled(0),
+            m_asserted_qhead(0), 
             m_assume_eq_head(0),
             m_num_conflicts(0),
             m_model_eqs(DEFAULT_HASHTABLE_INITIAL_CAPACITY, var_value_hash(*this), var_value_eq(*this)),
-            m_solver(nullptr),
-            m_resource_limit(*this) {
+            m_solver(0),
+            m_resource_limit(*this),
+            m_use_nra_model(false) {
         }
 
         ~imp() {
@@ -677,12 +737,8 @@ namespace smt {
         }
 
         bool internalize_atom(app * atom, bool gate_ctx) {
-            if (m_delay_constraints) {
-                return internalize_atom_lazy(atom, gate_ctx);
-            }
-            else {
-                return internalize_atom_strict(atom, gate_ctx);
-            }
+            return internalize_atom_strict(atom, gate_ctx);
+            
         }
 
         bool internalize_atom_strict(app * atom, bool gate_ctx) {
@@ -717,53 +773,10 @@ namespace smt {
             return true;
         }
 
-        bool internalize_atom_lazy(app * atom, bool gate_ctx) {
-            SASSERT(!ctx().b_internalized(atom));
-            bool_var bv = ctx().mk_bool_var(atom);
-            ctx().set_var_theory(bv, get_id());
-            expr* n1 = nullptr, *n2 = nullptr;
-            rational r;
-            lra_lp::bound_kind k;
-            theory_var v = null_theory_var;
-            scoped_internalize_state st(*this);
-            if (a.is_le(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) {
-                v = internalize_def(to_app(n1), st);
-                k = lra_lp::upper_t;
-            }
-            else if (a.is_ge(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) {
-                v = internalize_def(to_app(n1), st);
-                k = lra_lp::lower_t;
-            }
-            else {
-                TRACE("arith", tout << "Could not internalize " << mk_pp(atom, m) << "\n";);
-                found_not_handled(atom);
-                return true;
-            }
-            lra_lp::bound* b = alloc(lra_lp::bound, bv, v, r, k);
-            m_bounds[v].push_back(b);
-            updt_unassigned_bounds(v, +1);
-            m_bounds_trail.push_back(v);
-            m_bool_var2bound.insert(bv, b);
-            TRACE("arith", tout << "Internalized " << mk_pp(atom, m) << "\n";);
-            if (!is_unit_var(st) && m_bounds[v].size() == 1) {
-                m_delayed_defs.push_back(delayed_def(st.vars(), st.coeffs(), st.coeff(), v));
-            }
-            return true;
-        }
-
         bool internalize_term(app * term) {
             if (ctx().e_internalized(term) && th.is_attached_to_var(ctx().get_enode(term))) {
                 // skip
             }
-            else if (m_delay_constraints) {
-                scoped_internalize_state st(*this);
-                linearize_term(term, st);  // ensure that a theory_var was created.
-                SASSERT(ctx().e_internalized(term));
-                if(!th.is_attached_to_var(ctx().get_enode(term))) {
-                    mk_var(term);
-                }
-                m_delayed_terms.push_back(term);
-            }
             else {
                 internalize_def(term);
             }
@@ -789,13 +802,8 @@ namespace smt {
         }
 
         void new_eq_eh(theory_var v1, theory_var v2) {
-            if (m_delay_constraints) {
-                m_delayed_equalities.push_back(std::make_pair(v1, v2));
-            }
-            else {
-                // or internalize_eq(v1, v2);
-                m_arith_eq_adapter.new_eq_eh(v1, v2);
-            }
+            // or internalize_eq(v1, v2);
+            m_arith_eq_adapter.new_eq_eh(v1, v2);            
         }
 
         bool use_diseqs() const {
@@ -814,13 +822,11 @@ namespace smt {
             s.m_bounds_lim = m_bounds_trail.size();
             s.m_asserted_qhead = m_asserted_qhead;
             s.m_asserted_atoms_lim = m_asserted_atoms.size();
-            s.m_delayed_terms_lim = m_delayed_terms.size();
-            s.m_delayed_equalities_lim = m_delayed_equalities.size();
-            s.m_delayed_defs_lim = m_delayed_defs.size();
             s.m_not_handled = m_not_handled;
             s.m_underspecified_lim = m_underspecified.size();
             s.m_var_trail_lim = m_var_trail.size();
-            if (!m_delay_constraints) m_solver->push();
+            m_solver->push();
+            if (m_nra) m_nra->push();
         }
 
         void pop_scope_eh(unsigned num_scopes) {
@@ -841,18 +847,16 @@ namespace smt {
                 m_theory_var2var_index[m_var_trail[i]] = UINT_MAX;
             }
             m_asserted_atoms.shrink(m_scopes[old_size].m_asserted_atoms_lim);
-            m_delayed_terms.shrink(m_scopes[old_size].m_delayed_terms_lim);
-            m_delayed_defs.shrink(m_scopes[old_size].m_delayed_defs_lim);
-            m_delayed_equalities.shrink(m_scopes[old_size].m_delayed_equalities_lim);
             m_asserted_qhead = m_scopes[old_size].m_asserted_qhead;
             m_underspecified.shrink(m_scopes[old_size].m_underspecified_lim);
             m_var_trail.shrink(m_scopes[old_size].m_var_trail_lim);
             m_not_handled = m_scopes[old_size].m_not_handled;
-            m_scopes.resize(old_size);
-            if (!m_delay_constraints) m_solver->pop(num_scopes);
+            m_scopes.resize(old_size);            
+            m_solver->pop(num_scopes);
             // VERIFY(l_false != make_feasible());
             m_new_bounds.reset();
             m_to_check.reset();
+            if (m_nra) m_nra->pop(num_scopes);
             TRACE("arith", tout << "num scopes: " << num_scopes << " new scope level: " << m_scopes.size() << "\n";);
         }
 
@@ -1086,7 +1090,9 @@ namespace smt {
                   }
                   tout << "\n";
                   );
-            m_solver->random_update(vars.size(), vars.c_ptr());
+            if (!m_use_nra_model) {
+                m_solver->random_update(vars.size(), vars.c_ptr());
+            }
             m_model_eqs.reset();
             TRACE("arith", display(tout););
 
@@ -1135,55 +1141,70 @@ namespace smt {
                 enode* n1 = get_enode(v1);
                 enode* n2 = get_enode(v2);
                 m_assume_eq_head++;
-                CTRACE("arith",
-                       get_ivalue(v1) == get_ivalue(v2) && n1->get_root() != n2->get_root(),
+                CTRACE("arith", 
+                       is_eq(v1, v2) && n1->get_root() != n2->get_root(),
                        tout << "assuming eq: v" << v1 << " = v" << v2 << "\n";);
-                if (get_ivalue(v1) == get_ivalue(v2) &&  n1->get_root() != n2->get_root() && th.assume_eq(n1, n2)) {
+                if (is_eq(v1, v2) &&  n1->get_root() != n2->get_root() && th.assume_eq(n1, n2)) {
                     return true;
                 }
             }
             return false;
         }
 
+        bool is_eq(theory_var v1, theory_var v2) {
+            if (m_use_nra_model) {
+                return m_nra->am().eq(nl_value(v1, *m_a1), nl_value(v2, *m_a2));
+            }
+            else {
+                return get_ivalue(v1) == get_ivalue(v2); 
+            }
+        }
+
         bool has_delayed_constraints() const {
-            return !(m_asserted_atoms.empty() && m_delayed_terms.empty() && m_delayed_equalities.empty());
+            return !m_asserted_atoms.empty();
         }
 
         final_check_status final_check_eh() {
+            m_use_nra_model = false;
             lbool is_sat = l_true;
-            if (m_delay_constraints) {
-                init_solver();
-                for (unsigned i = 0; i < m_asserted_atoms.size(); ++i) {
-                    bool_var bv = m_asserted_atoms[i].m_bv;
-                    assert_bound(bv, m_asserted_atoms[i].m_is_true, *m_bool_var2bound.find(bv));
-                }
-                for (unsigned i = 0; i < m_delayed_terms.size(); ++i) {
-                    internalize_def(m_delayed_terms[i].get());
-                }
-                for (unsigned i = 0; i < m_delayed_defs.size(); ++i) {
-                    internalize_eq(m_delayed_defs[i]);
-                }
-                for (unsigned i = 0; i < m_delayed_equalities.size(); ++i) {
-                    std::pair<theory_var, theory_var> const& eq = m_delayed_equalities[i];
-                    internalize_eq(eq.first, eq.second);
-                }
-                is_sat = make_feasible();
-            }
-            else if (m_solver->get_status() != lp::lp_status::OPTIMAL) {
+            if (m_solver->get_status() != lean::lp_status::OPTIMAL) {
                 is_sat = make_feasible();
             }
+            final_check_status st = FC_DONE;
             switch (is_sat) {
             case l_true:
+                
                 if (delayed_assume_eqs()) {
                     return FC_CONTINUE;
                 }
                 if (assume_eqs()) {
                     return FC_CONTINUE;
                 }
-                if (m_not_handled != nullptr) {
-                    return FC_GIVEUP;
+
+                switch (check_lia()) {
+                case l_true:
+                    break;
+                case l_false:
+                    return FC_CONTINUE;
+                case l_undef:
+                    st = FC_GIVEUP;
+                    break;
                 }
-                return FC_DONE;
+                
+                switch (check_nra()) {
+                case l_true:
+                    break;
+                case l_false:
+                    return FC_CONTINUE;
+                case l_undef:
+                    st = FC_GIVEUP;
+                    break;
+                }
+                if (m_not_handled != 0) {                    
+                    st = FC_GIVEUP;
+                }
+
+                return st;
             case l_false:
                 set_conflict();
                 return FC_CONTINUE;
@@ -1196,6 +1217,70 @@ namespace smt {
             return FC_GIVEUP;
         }
 
+        // create a bound atom representing term >= k
+        lp::bound* mk_bound(lean::lar_term const& term, rational const& k) {
+            NOT_IMPLEMENTED_YET();
+            lp::bound_kind bkind = lp::bound_kind::lower_t;
+            bool_var bv = null_bool_var; 
+            theory_var v = null_theory_var;
+            lp::bound* result = alloc(lp::bound, bv, v, k, bkind);
+            return result;
+        }
+
+        lbool check_lia() {
+            std::cout << "called check_lia()\n";
+            lean::lar_term term;
+            lean::mpq k;
+            lean::explanation ex; // TBD, this should be streamlined accross different explanations
+            switch(m_lia->check(term, k, ex)) {
+            case lean::lia_move::ok:
+                return l_true;
+            case lean::lia_move::branch:
+                // branch on term <= k
+                NOT_IMPLEMENTED_YET();
+                return l_false;
+            case lean::lia_move::cut:
+                // m_explanation implies term <= k
+                m_explanation = ex.m_explanation;
+                NOT_IMPLEMENTED_YET();
+                return l_false;
+            case lean::lia_move::conflict:
+                // ex contains unsat core
+                m_explanation = ex.m_explanation;
+                set_conflict1();
+                return l_false;
+            case lean::lia_move::give_up:
+                return l_undef;
+            default:
+                UNREACHABLE();
+            }
+            return l_undef;
+        }
+
+        lbool check_nra() {
+            m_use_nra_model = false;
+            if (m.canceled()) return l_undef;
+            if (!m_nra) return l_true;
+            if (!m_nra->need_check()) return l_true;
+            m_a1 = 0; m_a2 = 0;
+            lbool r = m_nra->check(m_explanation);
+            m_a1 = alloc(scoped_anum, m_nra->am());
+            m_a2 = alloc(scoped_anum, m_nra->am());
+            switch (r) {
+            case l_false:
+                set_conflict1();
+                break;
+            case l_true:
+                m_use_nra_model = true;
+                if (assume_eqs()) {
+                    return l_false;
+                }
+                break;
+            default:
+                break;
+            }
+            return r;
+        }
 
         /**
            \brief We must redefine this method, because theory of arithmetic contains
@@ -1265,14 +1350,12 @@ namespace smt {
 #else
                 propagate_bound(bv, is_true, b);
 #endif
-                if (!m_delay_constraints) {
-                    lra_lp::bound& b = *m_bool_var2bound.find(bv);
-                    assert_bound(bv, is_true, b);
-                }
-
+                lp::bound& b = *m_bool_var2bound.find(bv);
+                assert_bound(bv, is_true, b);
+                
                 ++m_asserted_qhead;
             }
-            if (m_delay_constraints || ctx().inconsistent()) {
+            if (ctx().inconsistent()) {
                 m_to_check.reset();
                 return;
             }
@@ -2141,18 +2224,8 @@ namespace smt {
         }
 
         lbool make_feasible() {
-            reset_variable_values();
-            ++m_stats.m_make_feasible;
-            if (m_solver->A_r().column_count() > m_stats.m_max_cols)
-                m_stats.m_max_cols = m_solver->A_r().column_count();
-            if (m_solver->A_r().row_count() > m_stats.m_max_rows)
-                m_stats.m_max_rows = m_solver->A_r().row_count();
+            auto status = m_solver->find_feasible_solution();
             TRACE("arith_verbose", display(tout););
-            lp::lp_status status = m_solver->find_feasible_solution();
-            m_stats.m_num_iterations = m_solver->settings().st().m_total_iterations;
-            m_stats.m_num_factorizations = m_solver->settings().st().m_num_factorizations;
-            m_stats.m_need_to_solve_inf = m_solver->settings().st().m_need_to_solve_inf;
-
             switch (status) {
             case lp::lp_status::INFEASIBLE:
                 return l_false;
@@ -2205,11 +2278,15 @@ namespace smt {
         }
 
         void set_conflict() {
+            m_explanation.clear();
+            m_solver->get_infeasibility_explanation(m_explanation);
+            set_conflict1();
+        }
+
+        void set_conflict1() {
             m_eqs.reset();
             m_core.reset();
             m_params.reset();
-            m_explanation.clear();
-            m_solver->get_infeasibility_explanation(m_explanation);
             // m_solver->shrink_explanation_to_minimum(m_explanation); // todo, enable when perf is fixed
             /*
             static unsigned cn = 0;
@@ -2258,9 +2335,43 @@ namespace smt {
             TRACE("arith", display(tout););
         }
 
+        nlsat::anum const& nl_value(theory_var v, scoped_anum& r) {
+            SASSERT(m_nra);
+            SASSERT(m_use_nra_model);
+            lean::var_index vi = m_theory_var2var_index[v];
+            if (m_solver->is_term(vi)) {
+                lean::lar_term const& term = m_solver->get_term(vi);
+                scoped_anum r1(m_nra->am());
+                m_nra->am().set(r, term.m_v.to_mpq());
+                
+                for (auto const coeff : term.m_coeffs) {
+                    lean::var_index wi = coeff.first;
+                    m_nra->am().set(r1, coeff.second.to_mpq());
+                    m_nra->am().mul(m_nra->value(wi), r1, r1);
+                    m_nra->am().add(r1, r, r);
+                }
+                return r;
+            }
+            else {
+                return m_nra->value(vi);
+            }
+        }
+
         model_value_proc * mk_value(enode * n, model_generator & mg) {
             theory_var v = n->get_th_var(get_id());
-            return alloc(expr_wrapper_proc, m_factory->mk_value(get_value(v), m.get_sort(n->get_owner())));
+            expr* o = n->get_owner();
+            if (m_use_nra_model) {
+                anum const& an = nl_value(v, *m_a1);
+                if (a.is_int(o) && !m_nra->am().is_int(an)) {
+                    return alloc(expr_wrapper_proc, a.mk_numeral(rational::zero(), a.is_int(o)));
+                }
+                return alloc(expr_wrapper_proc, a.mk_numeral(nl_value(v, *m_a1), a.is_int(o)));
+            }
+            else {
+                rational r = get_value(v);
+                if (a.is_int(o) && !r.is_int()) r = floor(r);
+                return alloc(expr_wrapper_proc, m_factory->mk_value(r,  m.get_sort(o)));
+            }
         }
 
         bool get_value(enode* n, expr_ref& r) {
@@ -2286,6 +2397,7 @@ namespace smt {
             if (dump_lemmas()) {
                 ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), false_literal);
             }
+            if (m_arith_params.m_arith_mode == AS_LRA) return true;
             context nctx(m, ctx().get_fparams(), ctx().get_params());
             add_background(nctx);
             bool result = l_true != nctx.check();
@@ -2298,6 +2410,7 @@ namespace smt {
             if (dump_lemmas()) {
                 ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), lit);
             }
+            if (m_arith_params.m_arith_mode == AS_LRA) return true;
             context nctx(m, ctx().get_fparams(), ctx().get_params());
             m_core.push_back(~lit);
             add_background(nctx);
@@ -2309,6 +2422,7 @@ namespace smt {
         }
 
         bool validate_eq(enode* x, enode* y) {
+            if (m_arith_params.m_arith_mode == AS_LRA) return true;
             context nctx(m, ctx().get_fparams(), ctx().get_params());
             add_background(nctx);
             nctx.assert_expr(m.mk_not(m.mk_eq(x->get_owner(), y->get_owner())));
@@ -2509,7 +2623,7 @@ namespace smt {
             st.update("arith-rows", m_stats.m_add_rows);
             st.update("arith-propagations", m_stats.m_bounds_propagations);
             st.update("arith-iterations", m_stats.m_num_iterations);
-            st.update("arith-factorizations", m_stats.m_num_factorizations);
+            st.update("arith-factorizations", m_solver->settings().st().m_num_factorizations);
             st.update("arith-pivots", m_stats.m_need_to_solve_inf);
             st.update("arith-plateau-iterations", m_stats.m_num_iterations_with_no_progress);
             st.update("arith-fixed-eqs", m_stats.m_fixed_eqs);
@@ -2517,10 +2631,10 @@ namespace smt {
             st.update("arith-bound-propagations-lp", m_stats.m_bound_propagations1);
             st.update("arith-bound-propagations-cheap", m_stats.m_bound_propagations2);
             st.update("arith-diseq", m_stats.m_assert_diseq);
-            st.update("arith-make-feasible", m_stats.m_make_feasible);
-            st.update("arith-max-columns", m_stats.m_max_cols);
-            st.update("arith-max-rows", m_stats.m_max_rows);
-        }
+            st.update("arith-make-feasible", m_solver->settings().st().m_make_feasible);
+            st.update("arith-max-columns", m_solver->settings().st().m_max_cols);
+            st.update("arith-max-rows", m_solver->settings().st().m_max_rows);
+        }        
     };
 
     theory_lra::theory_lra(ast_manager& m, theory_arith_params& ap):

src/test/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_subdirectory(fuzzing)
+add_subdirectory(lp)
 ################################################################################
 # z3-test executable
 ################################################################################
@@ -120,6 +121,7 @@ add_executable(test-z3
   upolynomial.cpp
   var_subst.cpp
   vector.cpp
+  lp/lp.cpp
   ${z3_test_extra_object_files}
 )
 z3_add_install_tactic_rule(${z3_test_deps})

src/test/lp/CMakeLists.txt

@@ -0,0 +1,6 @@
+add_executable(lp_tst lp_main.cpp lp.cpp $<TARGET_OBJECTS:util> $<TARGET_OBJECTS:polynomial> $<TARGET_OBJECTS:nlsat>  $<TARGET_OBJECTS:lp> )
+target_compile_definitions(lp_tst PRIVATE ${Z3_COMPONENT_CXX_DEFINES})
+target_compile_options(lp_tst PRIVATE ${Z3_COMPONENT_CXX_FLAGS})
+target_include_directories(lp_tst PRIVATE ${Z3_COMPONENT_EXTRA_INCLUDE_DIRS})
+target_link_libraries(lp_tst PRIVATE ${Z3_DEPENDENT_LIBS})
+z3_append_linker_flag_list_to_target(lp_tst ${Z3_DEPENDENT_EXTRA_CXX_LINK_FLAGS})

src/test/lp/lp.cpp

@@ -2712,8 +2712,8 @@ void test_term() {
     lar_solver solver;
     unsigned _x = 0;
     unsigned _y = 1;
-    var_index x = solver.add_var(_x);
-    var_index y = solver.add_var(_y);
+    var_index x = solver.add_var(_x, false);
+    var_index y = solver.add_var(_y, false);
 
     vector<std::pair<mpq, var_index>> term_ls;
     term_ls.push_back(std::pair<mpq, var_index>((int)1, x));
@@ -2726,9 +2726,16 @@ void test_term() {
     ls.push_back(std::pair<mpq, var_index>((int)1, z));
     
     solver.add_constraint(ls, lconstraint_kind::EQ, mpq(0));
+    ls.clear();
+    ls.push_back(std::pair<mpq, var_index>((int)1, x));
+    solver.add_constraint(ls, lconstraint_kind::LT, mpq(0));
+    ls.push_back(std::pair<mpq, var_index>((int)2, y));
+    solver.add_constraint(ls, lconstraint_kind::GT, mpq(0));
     auto status = solver.solve();
     std::cout << lp_status_to_string(status) << std::endl;
     std::unordered_map<var_index, mpq> model;
+    if (status != OPTIMAL)
+        return;
     solver.get_model(model);
     
     for (auto & t : model) {
@@ -2740,8 +2747,8 @@ void test_term() {
 
 void test_evidence_for_total_inf_simple(argument_parser & args_parser) {
     lar_solver solver;
-    var_index x = solver.add_var(0);
-    var_index y = solver.add_var(1);
+    var_index x = solver.add_var(0, false);
+    var_index y = solver.add_var(1, false);
     solver.add_var_bound(x, LE, -mpq(1));
     solver.add_var_bound(y, GE, mpq(0));
     vector<std::pair<mpq, var_index>> ls;
@@ -2775,9 +2782,9 @@ If b becomes basic variable, then it is likely the old solver ends up with a row
             return true; 
         };   
     lar_solver ls;
-    unsigned a = ls.add_var(0);
-    unsigned b = ls.add_var(1);
-    unsigned c = ls.add_var(2);
+    unsigned a = ls.add_var(0, false);
+    unsigned b = ls.add_var(1, false);
+    unsigned c = ls.add_var(2, false);
     vector<std::pair<mpq, var_index>> coeffs;
     coeffs.push_back(std::pair<mpq, var_index>(1, a));
     coeffs.push_back(std::pair<mpq, var_index>(-1, c));
@@ -2840,8 +2847,8 @@ If x9 becomes basic variable, then it is likely the old solver ends up with a ro
 }
 void test_bound_propagation_one_row() {
     lar_solver ls;
-    unsigned x0 = ls.add_var(0);
-    unsigned x1 = ls.add_var(1);
+    unsigned x0 = ls.add_var(0, false);
+    unsigned x1 = ls.add_var(1, false);
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(1, x0));
     c.push_back(std::pair<mpq, var_index>(-1, x1));
@@ -2854,8 +2861,8 @@ void test_bound_propagation_one_row() {
 } 
 void test_bound_propagation_one_row_with_bounded_vars() {
     lar_solver ls;
-    unsigned x0 = ls.add_var(0);
-    unsigned x1 = ls.add_var(1);
+    unsigned x0 = ls.add_var(0, false);
+    unsigned x1 = ls.add_var(1, false);
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(1, x0));
     c.push_back(std::pair<mpq, var_index>(-1, x1));
@@ -2870,8 +2877,8 @@ void test_bound_propagation_one_row_with_bounded_vars() {
 }
 void test_bound_propagation_one_row_mixed() {
     lar_solver ls;
-    unsigned x0 = ls.add_var(0);
-    unsigned x1 = ls.add_var(1);
+    unsigned x0 = ls.add_var(0, false);
+    unsigned x1 = ls.add_var(1, false);
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(1, x0));
     c.push_back(std::pair<mpq, var_index>(-1, x1));
@@ -2885,9 +2892,9 @@ void test_bound_propagation_one_row_mixed() {
 
 void test_bound_propagation_two_rows() {
     lar_solver ls;
-    unsigned x = ls.add_var(0);
-    unsigned y = ls.add_var(1);
-    unsigned z = ls.add_var(2);
+    unsigned x = ls.add_var(0, false);
+    unsigned y = ls.add_var(1, false);
+    unsigned z = ls.add_var(2, false);
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(1, x));
     c.push_back(std::pair<mpq, var_index>(2, y));
@@ -2909,9 +2916,9 @@ void test_bound_propagation_two_rows() {
 void test_total_case_u() {
     std::cout << "test_total_case_u\n";
     lar_solver ls;
-    unsigned x = ls.add_var(0);
-    unsigned y = ls.add_var(1);
-    unsigned z = ls.add_var(2);
+    unsigned x = ls.add_var(0, false);
+    unsigned y = ls.add_var(1, false);
+    unsigned z = ls.add_var(2, false);
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(1, x));
     c.push_back(std::pair<mpq, var_index>(2, y));
@@ -2935,9 +2942,9 @@ bool contains_j_kind(unsigned j, lconstraint_kind kind, const mpq & rs, const ve
 void test_total_case_l(){
     std::cout << "test_total_case_l\n";
     lar_solver ls;
-    unsigned x = ls.add_var(0);
-    unsigned y = ls.add_var(1);
-    unsigned z = ls.add_var(2);
+    unsigned x = ls.add_var(0, false);
+    unsigned y = ls.add_var(1, false);
+    unsigned z = ls.add_var(2, false);
     vector<std::pair<mpq, var_index>> c;
     c.push_back(std::pair<mpq, var_index>(1, x));
     c.push_back(std::pair<mpq, var_index>(2, y));

src/test/lp/lp_main.cpp

@@ -0,0 +1,14 @@
+void        gparams_register_modules(){}
+void mem_initialize() {}
+void mem_finalize() {}
+#include "util/rational.h"
+namespace lean {
+void test_lp_local(int argc, char**argv);
+}
+int main(int argn, char**argv){
+    rational::initialize();
+    lean::test_lp_local(argn, argv);
+    rational::finalize();
+    return 0;
+}
+

src/test/lp/smt_reader.h

@@ -389,7 +389,7 @@ namespace lp {
         void add_constraint_to_solver(lar_solver * solver, formula_constraint & fc) {
             vector<std::pair<mpq, var_index>> ls;
             for (auto & it : fc.m_coeffs) {
-                ls.push_back(std::make_pair(it.first, solver->add_var(register_name(it.second))));
+                ls.push_back(std::make_pair(it.first, solver->add_var(register_name(it.second), false)));
             }
             solver->add_constraint(ls, fc.m_kind, fc.m_right_side);
         }

src/util/lp/CMakeLists.txt

@@ -8,6 +8,8 @@ z3_add_component(lp
     dense_matrix_instances.cpp
     eta_matrix_instances.cpp
     indexed_vector_instances.cpp
+    int_solver.cpp
+    lar_solver_instances.cpp
     lar_core_solver_instances.cpp
     lp_core_solver_base_instances.cpp
     lp_dual_core_solver_instances.cpp
@@ -18,6 +20,7 @@ z3_add_component(lp
     lp_solver_instances.cpp
     lu_instances.cpp
     matrix_instances.cpp
+    nra_solver.cpp
     permutation_matrix_instances.cpp
     quick_xplain.cpp
     row_eta_matrix_instances.cpp
@@ -28,6 +31,8 @@ z3_add_component(lp
     random_updater_instances.cpp      
   COMPONENT_DEPENDENCIES
     util
+    polynomial
+    nlsat
   PYG_FILES
     lp_params.pyg
 )

src/util/lp/indexed_vector.h

@@ -98,8 +98,8 @@ public:
             loc = one_of_type<T>(); // use as a characteristic function
         } 
     }
-    
 
+    
     void clear();
     void clear_all();
     const T& operator[] (unsigned i) const {

src/util/lp/init_lar_solver.h

@@ -1,591 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-
-// here we are inside lp::lar_solver class
-
-bool strategy_is_undecided() const {
-    return m_settings.simplex_strategy() == simplex_strategy_enum::undecided;
-}
-
-var_index add_var(unsigned ext_j) {
-    var_index i;
-    SASSERT (ext_j < m_terms_start_index); 
-
-    if (ext_j >= m_terms_start_index)
-        throw 0; // todo : what is the right way to exit?
-            
-    if (try_get_val(m_ext_vars_to_columns, ext_j, i)) {
-        return i;
-    }
-    SASSERT(m_vars_to_ul_pairs.size() == A_r().column_count());
-    i = A_r().column_count();
-    m_vars_to_ul_pairs.push_back (ul_pair(static_cast<unsigned>(-1)));
-    add_non_basic_var_to_core_fields(ext_j);
-    SASSERT(sizes_are_correct());
-    return i;
-}
-
-void register_new_ext_var_index(unsigned ext_v) {
-    SASSERT(!contains(m_ext_vars_to_columns, ext_v));
-    unsigned j = static_cast<unsigned>(m_ext_vars_to_columns.size());
-    m_ext_vars_to_columns[ext_v] = j;
-    SASSERT(m_columns_to_ext_vars_or_term_indices.size() == j);
-    m_columns_to_ext_vars_or_term_indices.push_back(ext_v);
-}
-
-void add_non_basic_var_to_core_fields(unsigned ext_j) {
-    register_new_ext_var_index(ext_j);
-    m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
-    m_columns_with_changed_bound.increase_size_by_one();
-    add_new_var_to_core_fields_for_mpq(false);
-    if (use_lu())
-        add_new_var_to_core_fields_for_doubles(false);
-}
-
-void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
-    unsigned j = A_d().column_count();
-    A_d().add_column();
-    SASSERT(m_mpq_lar_core_solver.m_d_x.size() == j);
-    //        SASSERT(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j);  // restore later
-    m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); 
-    m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
-    m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
-    SASSERT(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
-    if (register_in_basis) {
-        A_d().add_row();
-        m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
-        m_mpq_lar_core_solver.m_d_basis.push_back(j);
-    }else {
-        m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
-        m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
-    }
-}
-
-void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
-    unsigned j = A_r().column_count();
-    A_r().add_column();
-    SASSERT(m_mpq_lar_core_solver.m_r_x.size() == j);
-    //        SASSERT(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j);  // restore later
-    m_mpq_lar_core_solver.m_r_x.resize(j + 1);
-    m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one();
-    m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
-    m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one();
-    m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
-    m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
-    SASSERT(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method
-    if (register_in_basis) {
-        A_r().add_row();
-        m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size());
-        m_mpq_lar_core_solver.m_r_basis.push_back(j);
-        if (m_settings.bound_propagation())
-            m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
-    } else {
-        m_mpq_lar_core_solver.m_r_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1);
-        m_mpq_lar_core_solver.m_r_nbasis.push_back(j);
-    }
-}
-
-
-var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
-                             const mpq &m_v) {
-    m_terms.push_back(new lar_term(coeffs, m_v));
-    m_orig_terms.push_back(new lar_term(coeffs, m_v));
-    return m_terms_start_index +  m_terms.size() - 1;
-}
-
-// terms
-var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
-                   const mpq &m_v) {
-    if (strategy_is_undecided())
-        return add_term_undecided(coeffs, m_v);
-
-    m_terms.push_back(new lar_term(coeffs, m_v));
-    m_orig_terms.push_back(new lar_term(coeffs, m_v));
-    unsigned adjusted_term_index = m_terms.size() - 1;
-    var_index ret = m_terms_start_index + adjusted_term_index;
-    if (use_tableau() && !coeffs.empty()) {
-        add_row_for_term(m_orig_terms.back(), ret);
-        if (m_settings.bound_propagation())
-            m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
-    } 
-    SASSERT(m_ext_vars_to_columns.size() == A_r().column_count());
-    return ret;
-}
-
-void add_row_for_term(const lar_term * term, unsigned term_ext_index) {
-    SASSERT(sizes_are_correct());
-    add_row_from_term_no_constraint(term, term_ext_index);
-    SASSERT(sizes_are_correct());
-}
-
-void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
-    register_new_ext_var_index(term_ext_index);
-    // j will be a new variable
-    unsigned j = A_r().column_count();
-    ul_pair ul(j);
-    m_vars_to_ul_pairs.push_back(ul);
-    add_basic_var_to_core_fields();
-    if (use_tableau()) {
-        auto it = iterator_on_term_with_basis_var(*term, j);
-        A_r().fill_last_row_with_pivoting(it,
-                                          m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
-        m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
-    } else {
-        fill_last_row_of_A_r(A_r(), term);
-    }
-    m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
-    if (use_lu())
-        fill_last_row_of_A_d(A_d(), term);
-}
-
-void add_basic_var_to_core_fields() {
-    bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver();
-    SASSERT(!use_lu || A_r().column_count() == A_d().column_count());
-    m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
-    m_columns_with_changed_bound.increase_size_by_one();
-    m_rows_with_changed_bounds.increase_size_by_one();
-    add_new_var_to_core_fields_for_mpq(true);
-    if (use_lu)
-        add_new_var_to_core_fields_for_doubles(true);
-}
-
-constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side)  {
-    constraint_index ci = m_constraints.size();
-    if (!is_term(j)) { // j is a var
-        auto vc = new lar_var_constraint(j, kind, right_side);
-        m_constraints.push_back(vc);
-        update_column_type_and_bound(j, kind, right_side, ci);
-    } else {
-        add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
-    }
-    SASSERT(sizes_are_correct());
-    return ci;
-}
-
-void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
-    switch(m_mpq_lar_core_solver.m_column_types[j]) {
-    case column_type::free_column:
-        update_free_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::boxed:
-        update_boxed_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::low_bound:
-        update_low_bound_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::upper_bound:
-        update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    case column_type::fixed:
-        update_fixed_column_type_and_bound(j, kind, right_side, constr_index);
-        break;
-    default:
-        SASSERT(false); // cannot be here
-    }
-}
-
-void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    SASSERT(is_term(j));
-    unsigned adjusted_term_index = adjust_term_index(j);
-    unsigned term_j;
-    if (try_get_val(m_ext_vars_to_columns, j, term_j)) {
-        mpq rs = right_side - m_orig_terms[adjusted_term_index]->m_v;
-        m_constraints.push_back(new lar_term_constraint(m_orig_terms[adjusted_term_index], kind, right_side));
-        update_column_type_and_bound(term_j, kind, rs, ci);
-    }
-    else {
-        add_constraint_from_term_and_create_new_column_row(j, m_orig_terms[adjusted_term_index], kind, right_side);
-    }
-}
-
-
-void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
-                                                        lconstraint_kind kind, const mpq & right_side) {
-
-    add_row_from_term_no_constraint(term, term_j);
-    unsigned j = A_r().column_count() - 1;
-    update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
-    m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
-    SASSERT(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
-}
-
-void decide_on_strategy_and_adjust_initial_state() {
-    SASSERT(strategy_is_undecided());
-    if (m_vars_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
-        m_settings.simplex_strategy() = simplex_strategy_enum::lu;
-    } else {
-        m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs?
-    }
-    adjust_initial_state();
-}
-
-void adjust_initial_state() {
-    switch (m_settings.simplex_strategy()) {
-    case simplex_strategy_enum::lu:
-        adjust_initial_state_for_lu();
-        break;
-    case simplex_strategy_enum::tableau_rows:
-        adjust_initial_state_for_tableau_rows();
-        break;
-    case simplex_strategy_enum::tableau_costs:
-        SASSERT(false); // not implemented
-    case simplex_strategy_enum::undecided:
-        adjust_initial_state_for_tableau_rows();
-        break;
-    }
-}
-
-void adjust_initial_state_for_lu() {
-    copy_from_mpq_matrix(A_d());
-    unsigned n = A_d().column_count();
-    m_mpq_lar_core_solver.m_d_x.resize(n);
-    m_mpq_lar_core_solver.m_d_low_bounds.resize(n);
-    m_mpq_lar_core_solver.m_d_upper_bounds.resize(n);
-    m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading;
-    m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis;
-
-    /*
-    unsigned j = A_d().column_count();
-    A_d().add_column();
-    SASSERT(m_mpq_lar_core_solver.m_d_x.size() == j);
-    //        SASSERT(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j);  // restore later
-    m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); 
-    m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
-    m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
-    SASSERT(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
-    if (register_in_basis) {
-        A_d().add_row();
-        m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
-        m_mpq_lar_core_solver.m_d_basis.push_back(j);
-    }else {
-        m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
-        m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
-        }*/
-}
-
-void adjust_initial_state_for_tableau_rows() {
-    for (unsigned j = 0; j < m_terms.size(); j++) {
-        if (contains(m_ext_vars_to_columns, j + m_terms_start_index))
-            continue;
-        add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index);
-    }
-}
-
-// this fills the last row of A_d and sets the basis column: -1 in the last column of the row
-void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
-    SASSERT(A.row_count() > 0);
-    SASSERT(A.column_count() > 0);
-    unsigned last_row = A.row_count() - 1;
-    SASSERT(A.m_rows[last_row].empty());
-
-    for (auto & t : ls->m_coeffs) {
-        SASSERT(!is_zero(t.second));
-        var_index j = t.first;
-        A.set(last_row, j, - t.second.get_double());
-    }
-
-    unsigned basis_j = A.column_count() - 1;
-    A.set(last_row, basis_j, - 1 );
-}
-
-void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
-    mpq y_of_bound(0);
-    switch (kind) {
-    case LT:
-        y_of_bound = -1;
-    case LE:
-        m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound;
-        SASSERT(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
-        SASSERT(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
-        {
-            auto up = numeric_pair<mpq>(right_side, y_of_bound);
-            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-        }
-        set_upper_bound_witness(j, constr_ind);
-        break;
-    case GT:
-        y_of_bound = 1;
-    case GE:
-        m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound;
-        SASSERT(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
-        {
-            auto low = numeric_pair<mpq>(right_side, y_of_bound);
-            m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-        }
-        set_low_bound_witness(j, constr_ind);
-        break;
-    case EQ:
-        m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-        m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-        set_upper_bound_witness(j, constr_ind);
-        set_low_bound_witness(j, constr_ind);
-        break;
-
-    default:
-        SASSERT(false);
-                
-    }
-    m_columns_with_changed_bound.insert(j);
-}
-
-void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    SASSERT(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
-    mpq y_of_bound(0);
-    switch (kind) {
-    case LT:
-        y_of_bound = -1;
-    case LE:
-        {
-            auto up = numeric_pair<mpq>(right_side, y_of_bound);
-            if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) {
-                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-                set_upper_bound_witness(j, ci);
-                m_columns_with_changed_bound.insert(j);
-            }
-        }
-        break;
-    case GT:
-        y_of_bound = 1;
-    case GE:            
-        m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed;
-        {
-            auto low = numeric_pair<mpq>(right_side, y_of_bound);
-            m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-            set_low_bound_witness(j, ci);
-            m_columns_with_changed_bound.insert(j);
-            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-            } else {
-                m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
-            }                     
-        }
-        break;
-    case EQ:
-        {
-            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_status = INFEASIBLE;
-                set_low_bound_witness(j, ci);
-                m_infeasible_column_index = j;
-            } else {
-                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
-                m_columns_with_changed_bound.insert(j);
-                set_low_bound_witness(j, ci);
-                set_upper_bound_witness(j, ci);
-                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-            }
-            break;
-        }
-        break;
-
-    default:
-        SASSERT(false);
-                
-    }
-}
-    
-void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    SASSERT(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
-    mpq y_of_bound(0);
-    switch (kind) {
-    case LT:
-        y_of_bound = -1;
-    case LE:
-        {
-            auto up = numeric_pair<mpq>(right_side, y_of_bound);
-            if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-                set_upper_bound_witness(j, ci);
-                m_columns_with_changed_bound.insert(j);
-            }
-
-            if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_status = INFEASIBLE;
-                SASSERT(false);
-                m_infeasible_column_index = j;
-            } else {
-                if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])
-                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-            }                    
-        }
-        break;
-    case GT:
-        y_of_bound = 1;
-    case GE:            
-        {
-            auto low = numeric_pair<mpq>(right_side, y_of_bound);
-            if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-                m_columns_with_changed_bound.insert(j);
-                set_low_bound_witness(j, ci);
-            }
-            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-            } else if ( low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-            }
-        }
-        break;
-    case EQ:
-        {
-            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_upper_bound_witness(j, ci);                    
-            } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_low_bound_witness(j, ci);                    
-            } else {
-                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
-                set_low_bound_witness(j, ci);
-                set_upper_bound_witness(j, ci);
-                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-                m_columns_with_changed_bound.insert(j);
-            }
-                
-            break;
-        }
-
-    default:
-        SASSERT(false);
-                
-    }
-}
-void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    SASSERT(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound);
-    mpq y_of_bound(0);
-    switch (kind) {
-    case LT:
-        y_of_bound = -1;
-    case LE:
-        {
-            auto up = numeric_pair<mpq>(right_side, y_of_bound);
-            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-            set_upper_bound_witness(j, ci);
-            m_columns_with_changed_bound.insert(j);
-
-            if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-            } else {
-                m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
-            }                    
-        }
-        break;
-    case GT:
-        y_of_bound = 1;
-    case GE:            
-        {
-            auto low = numeric_pair<mpq>(right_side, y_of_bound);
-            if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
-                m_columns_with_changed_bound.insert(j);
-                set_low_bound_witness(j, ci);
-            }
-        }
-        break;
-    case EQ:
-        {
-            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_upper_bound_witness(j, ci);                    
-            } else {
-                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
-                set_low_bound_witness(j, ci);
-                set_upper_bound_witness(j, ci);
-                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-            }
-            m_columns_with_changed_bound.insert(j);
-            break;
-        }
-
-    default:
-        SASSERT(false);
-                
-    }
-}
-
-void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    SASSERT(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
-    SASSERT(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
-    auto v = numeric_pair<mpq>(right_side, mpq(0));
-        
-    mpq y_of_bound(0);
-    switch (kind) {
-    case LT:
-        if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-            m_status = INFEASIBLE;
-            m_infeasible_column_index = j;
-            set_upper_bound_witness(j, ci);
-        }                   
-        break;
-    case LE:
-        {
-            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_upper_bound_witness(j, ci);
-            }                   
-        }
-        break;
-    case GT:
-        {
-            if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index =j;
-                set_low_bound_witness(j, ci);
-            }
-        }
-        break;
-    case GE:            
-        {
-            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_low_bound_witness(j, ci);
-            }
-        }
-        break;
-    case EQ:
-        {
-            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_upper_bound_witness(j, ci);                    
-            } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                m_status = INFEASIBLE;
-                m_infeasible_column_index = j;
-                set_low_bound_witness(j, ci);                    
-            } 
-            break;
-        }
-
-    default:
-        SASSERT(false);
-                
-    }
-}
-    

src/util/lp/int_solver.cpp

@@ -0,0 +1,606 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#include "util/lp/int_solver.h"
+#include "util/lp/lar_solver.h"
+namespace lean {
+
+void int_solver::fix_non_base_columns() {
+    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
+    for (unsigned j : lcs.m_r_nbasis) {
+        if (column_is_int_inf(j)) {
+            set_value(j, floor(lcs.m_r_x[j].x));
+        }
+    }
+    if (m_lar_solver->find_feasible_solution() == INFEASIBLE)
+        failed();
+}
+
+void int_solver::failed() {
+    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
+    
+    for (unsigned j : m_old_values_set.m_index) {
+        lcs.m_r_x[j] = m_old_values_data[j];
+        lean_assert(lcs.m_r_solver.column_is_feasible(j));
+        lcs.m_r_solver.remove_column_from_inf_set(j);
+    }
+    lean_assert(lcs.m_r_solver.calc_current_x_is_feasible_include_non_basis());
+    lean_assert(lcs.m_r_solver.current_x_is_feasible());
+    m_old_values_set.clear();
+}
+
+void int_solver::trace_inf_rows() const {
+    unsigned num = m_lar_solver->A_r().column_count();
+    for (unsigned v = 0; v < num; v++) {
+        if (is_int(v) && !get_value(v).is_int()) {
+            display_column(tout, v);
+        }
+    }
+    
+    num = 0;
+    for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) {
+        unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i];
+        if (column_is_int_inf(j)) {
+            num++;
+            iterator_on_row<mpq> it(m_lar_solver->A_r().m_rows[i]);
+            m_lar_solver->print_linear_iterator(&it, tout);
+            tout << "\n";
+        }
+    }
+    tout << "num of int infeasible: " << num << "\n";
+}
+
+int int_solver::find_inf_int_base_column() {
+    if (m_inf_int_set.is_empty())
+        return -1;
+    int j = find_inf_int_boxed_base_column_with_smallest_range();
+    if (j != -1)
+        return j;
+    unsigned k = settings().random_next() % m_inf_int_set.m_index.size();
+    return m_inf_int_set.m_index[k];
+}
+
+int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
+    int result = -1;
+    mpq range;
+    mpq new_range;
+    mpq small_range_thresold(1024);
+    unsigned n = 0;
+    lar_core_solver & lcs = m_lar_solver->m_mpq_lar_core_solver;
+
+    for (int j : m_inf_int_set.m_index) {
+        lean_assert(is_base(j) && column_is_int_inf(j));
+        if (!is_boxed(j))
+            continue;
+        new_range  = lcs.m_r_upper_bounds()[j].x - lcs.m_r_low_bounds()[j].x;
+        if (new_range > small_range_thresold) 
+            continue;
+        if (result == -1) {
+            result = j;
+            range  = new_range;
+            n      = 1;
+            continue;
+        }
+        if (new_range < range) {
+            n      = 1;
+            result = j;
+            range  = new_range;
+            continue;
+        }
+        if (new_range == range) {
+            n++;
+            if (settings().random_next() % n == 0) {
+                result = j;
+                continue;
+            }
+        }
+    }
+    return result;
+    
+}
+
+lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
+    lean_assert(is_feasible());
+    init_inf_int_set();
+    lean_assert(inf_int_set_is_correct());
+    // currently it is a reimplementation of
+    // final_check_status theory_arith<Ext>::check_int_feasibility()
+    // from theory_arith_int.h
+	if (m_lar_solver->model_is_int_feasible())
+		return lia_move::ok;
+    if (!gcd_test(ex))
+        return lia_move::conflict;
+    /*
+      if (m_params.m_arith_euclidean_solver)
+            apply_euclidean_solver();
+        
+    */
+    m_lar_solver->pivot_fixed_vars_from_basis();
+    patch_int_infeasible_columns();
+    fix_non_base_columns();
+    lean_assert(is_feasible());
+    TRACE("arith_int_rows", trace_inf_rows(););
+    
+    if (find_inf_int_base_column() == -1)
+        return lia_move::ok;
+        
+
+    if ((++m_branch_cut_counter) % settings().m_int_branch_cut_threshold == 0) {
+        move_non_base_vars_to_bounds();
+        /*
+        if (!make_feasible()) {
+            TRACE("arith_int", tout << "failed to move variables to bounds.\n";);
+            failed();
+            return FC_CONTINUE;
+        }
+        int int_var = find_inf_int_base_var();
+        if (int_var != null_int) {
+            TRACE("arith_int", tout << "v" << int_var << " does not have an integer assignment: " << get_value(int_var) << "\n";);
+            SASSERT(is_base(int_var));
+            row const & r = m_rows[get_var_row(int_var)];
+            if (!mk_gomory_cut(r)) {
+                // silent failure
+            }
+            return FC_CONTINUE;
+            }*/
+    }
+    else {
+        int j = find_inf_int_base_column();
+        /*
+        if (j != -1) {
+            TRACE("arith_int", tout << "v" << j << " does not have an integer assignment: " << get_value(j) << "\n";);
+            // apply branching 
+            branch_infeasible_int_var(int_var);
+            return false;
+            }*/
+    }
+    //    return true;
+    return lia_move::give_up;
+}
+
+void int_solver::move_non_base_vars_to_bounds() {
+    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
+    for (unsigned j : lcs.m_r_nbasis) {
+        auto & val = lcs.m_r_x[j];
+        switch (lcs.m_column_types()[j]) {
+        case column_type::boxed:
+            if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j])
+                set_value(j, lcs.m_r_low_bounds()[j]);
+            break;
+        case column_type::low_bound:
+            if (val != lcs.m_r_low_bounds()[j])
+                set_value(j, lcs.m_r_low_bounds()[j]);
+            break;
+        case column_type::upper_bound:
+            if (val != lcs.m_r_upper_bounds()[j])
+                set_value(j, lcs.m_r_upper_bounds()[j]);
+            break;
+        default:
+            if (is_int(j) && !val.is_int()) {
+                set_value(j, impq(floor(val)));
+            }
+        }
+    }
+}
+
+
+
+void int_solver::set_value(unsigned j, const impq & new_val) {
+    auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j];
+    if (!m_old_values_set.contains(j)) {
+        m_old_values_set.insert(j);
+        m_old_values_data[j] = x;
+    }
+    auto delta = new_val - x;
+    x = new_val;
+    m_lar_solver->change_basic_x_by_delta_on_column(j, delta);
+    update_column_in_inf_set_set(j);
+}
+
+void int_solver::patch_int_infeasible_columns() {
+    bool inf_l, inf_u;
+    impq l, u;
+    mpq m;
+    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
+    for (unsigned j : lcs.m_r_nbasis) {
+        if (!is_int(j))
+            continue;
+        get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m);
+        impq & val = lcs.m_r_x[j];
+        bool val_is_int = val.is_int();
+        bool m_is_one = m.is_one();
+        if (m.is_one() && val_is_int)
+            continue;
+        // check whether value of j is already a multiple of m.
+        if (val_is_int && (val.x / m).is_int())
+            continue;
+        TRACE("patch_int",
+              tout << "TARGET j" << j << " -> [";
+              if (inf_l) tout << "-oo"; else tout << l;
+              tout << ", ";
+              if (inf_u) tout << "oo"; else tout << u;
+              tout << "]";
+              tout << ", m: " << m << ", val: " << val << ", is_int: " << m_lar_solver->column_is_int(j) << "\n";);
+        if (!inf_l) {
+            l = m_is_one? ceil(l) : m * ceil(l / m);
+            if (inf_u || l <= u) {
+                TRACE("patch_int",
+                      tout << "patching with l: " << l << '\n';);
+                
+                set_value(j, l);
+            } else {
+                TRACE("patch_int", 
+                      tout << "not patching " << l << "\n";);
+            }
+        } else if (!inf_u) {
+            u = m_is_one? floor(u) : m * floor(u / m);
+            set_value(j, u);
+            TRACE("patch_int",
+                  tout << "patching with u: " << u << '\n';);
+        } else {
+            set_value(j, impq(0));
+            TRACE("patch_int",
+                  tout << "patching with 0\n";);
+        }
+    }
+}
+
+mpq get_denominators_lcm(iterator_on_row<mpq> &it) {
+    mpq r(1);
+    mpq a;
+    unsigned j;
+    while (it.next(a, j)) {
+        r = lcm(r, denominator(a));
+    }
+    return r;
+}
+    
+bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, unsigned i, explanation & ex) {
+    iterator_on_row<mpq> it(A.m_rows[i]);
+    std::cout << "gcd_test_for_row(" << i << ")\n";
+    mpq lcm_den = get_denominators_lcm(it);
+    mpq consts(0);
+    mpq gcds(0);
+    mpq least_coeff(0);
+    bool    least_coeff_is_bounded = false;
+    mpq a;
+    unsigned j;
+    while (it.next(a, j)) {
+            if (m_lar_solver->column_is_fixed(j)) {
+                mpq aux = lcm_den * a;
+                consts += aux * m_lar_solver->column_low_bound(j).x;
+            }
+            else if (m_lar_solver->column_is_real(j)) {
+                return true;
+            }
+            else if (gcds.is_zero()) {
+                gcds = abs(lcm_den * a);
+                least_coeff = gcds;
+                least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
+            }
+            else {
+                mpq aux = abs(lcm_den * a);
+                gcds = gcd(gcds, aux);
+                if (aux < least_coeff) {
+                    least_coeff = aux;
+                    least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
+                }
+                else if (least_coeff_is_bounded && aux == least_coeff) {
+                    least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
+                }
+            }
+            SASSERT(gcds.is_int());
+            SASSERT(least_coeff.is_int());
+            TRACE("gcd_test_bug", tout << "coeff: " << a << ", gcds: " << gcds 
+                  << " least_coeff: " << least_coeff << " consts: " << consts << "\n";);
+        
+        }
+    
+        if (gcds.is_zero()) {
+            // All variables are fixed.
+            // This theory guarantees that the assignment satisfies each row, and
+            // fixed integer variables are assigned to integer values.
+            return true;
+        }
+        
+        if (!(consts / gcds).is_int())
+            fill_explanation_from_fixed_columns(it, ex);
+        
+        if (least_coeff.is_one() && !least_coeff_is_bounded) {
+            SASSERT(gcds.is_one());
+            return true;
+        }
+        
+        if (least_coeff_is_bounded) {
+            return ext_gcd_test(it, least_coeff, lcm_den, consts, ex);
+        }
+        return true;
+}
+
+void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation & ex) {
+    constraint_index lc, uc;
+    m_lar_solver->get_bound_constraint_witnesses_for_column(j, lc, uc);
+    ex.m_explanation.push_back(std::make_pair(mpq(1), lc));
+    ex.m_explanation.push_back(std::make_pair(mpq(1), uc));
+}
+void int_solver::fill_explanation_from_fixed_columns(iterator_on_row<mpq> & it, explanation & ex) {
+    it.reset();
+    unsigned j;
+    while (it.next(j)) {
+        if (!m_lar_solver->column_is_fixed(j))
+            continue;
+        add_to_explanation_from_fixed_or_boxed_column(j, ex);
+    }
+}
+    
+bool int_solver::gcd_test(explanation & ex) {
+	auto & A = m_lar_solver->A_r(); // getting the matrix
+	for (unsigned i = 0; i < A.row_count(); i++)
+        if (!gcd_test_for_row(A, i, ex)) {
+            std::cout << "false from gcd_test\n" ;
+            return false;
+        }
+        
+	return true;
+}
+
+bool int_solver::ext_gcd_test(iterator_on_row<mpq> & it,
+                              mpq const & least_coeff, 
+                              mpq const & lcm_den,
+                              mpq const & consts, explanation& ex) {
+
+    std::cout << "calling ext_gcd_test" << std::endl;
+    mpq gcds(0);
+    mpq l(consts);
+    mpq u(consts);
+
+    it.reset();
+    mpq a;
+    unsigned j;
+    while (it.next(a, j)) {
+        if (m_lar_solver->column_is_fixed(j))
+            continue;
+        SASSERT(!m_lar_solver->column_is_real(j));
+        mpq ncoeff = lcm_den * a;
+        SASSERT(ncoeff.is_int());
+        mpq abs_ncoeff = abs(ncoeff);
+        if (abs_ncoeff == least_coeff) {
+            SASSERT(m_lar_solver->column_is_bounded(j));
+            if (ncoeff.is_pos()) {
+                // l += ncoeff * m_lar_solver->column_low_bound(j).x;
+                l.addmul(ncoeff, m_lar_solver->column_low_bound(j).x);
+                // u += ncoeff * m_lar_solver->column_upper_bound(j).x;
+                u.addmul(ncoeff, m_lar_solver->column_upper_bound(j).x);
+            }
+            else {
+                // l += ncoeff * upper_bound(j).get_rational();
+                l.addmul(ncoeff, m_lar_solver->column_upper_bound(j).x);
+                // u += ncoeff * lower_bound(j).get_rational();
+                u.addmul(ncoeff, m_lar_solver->column_low_bound(j).x);
+            }
+            add_to_explanation_from_fixed_or_boxed_column(j, ex);
+        }
+        else if (gcds.is_zero()) {
+            gcds = abs_ncoeff; 
+        }
+        else {
+            gcds = gcd(gcds, abs_ncoeff);
+        }
+        SASSERT(gcds.is_int());
+    }
+        
+    if (gcds.is_zero()) {
+        return true;
+    }
+        
+    mpq l1 = ceil(l/gcds);
+    mpq u1 = floor(u/gcds);
+        
+    if (u1 < l1) {
+        fill_explanation_from_fixed_columns(it, ex);
+        return false;
+    }
+        
+    return true;
+
+}
+
+linear_combination_iterator<mpq> * int_solver::get_column_iterator(unsigned j) {
+    if (m_lar_solver->use_tableau())
+        return new iterator_on_column<mpq, impq>(m_lar_solver->A_r().m_columns[j], m_lar_solver->A_r());
+   return new iterator_on_indexed_vector<mpq>(m_lar_solver->get_column_in_lu_mode(j));
+}
+
+
+int_solver::int_solver(lar_solver* lar_slv) :
+    m_lar_solver(lar_slv),
+    m_branch_cut_counter(0) {
+    lean_assert(m_old_values_set.size() == 0);
+    m_old_values_set.resize(lar_slv->A_r().column_count());
+    m_old_values_data.resize(lar_slv->A_r().column_count(), zero_of_type<impq>());    
+}
+
+bool int_solver::lower(unsigned j) const {
+    switch (m_lar_solver->m_mpq_lar_core_solver.m_column_types()[j]) {
+    case column_type::fixed:
+    case column_type::boxed:
+    case column_type::low_bound:
+        return true;
+    default:
+        return false;
+    }
+}
+
+bool int_solver::upper(unsigned j) const {
+    switch (m_lar_solver->m_mpq_lar_core_solver.m_column_types()[j]) {
+    case column_type::fixed:
+    case column_type::boxed:
+    case column_type::upper_bound:
+        return true;
+    default:
+        return false;
+    }
+}
+
+const impq& int_solver::lower_bound(unsigned j) const {
+    return m_lar_solver->m_mpq_lar_core_solver.m_r_low_bounds()[j];
+}
+
+const impq& int_solver::upper_bound(unsigned j) const {
+    return m_lar_solver->m_mpq_lar_core_solver.m_r_upper_bounds()[j];
+}
+
+
+void set_lower(impq & l,
+               bool & inf_l,
+               impq const & v ) {
+    if (inf_l || v > l) {
+        l = v;
+        inf_l = false;
+    }
+}
+
+void set_upper(impq & u,
+               bool & inf_u,
+               impq const & v) {
+    if (inf_u || v < u) {
+        u = v;
+        inf_u = false;
+    }
+}
+
+bool int_solver::get_freedom_interval_for_column(unsigned x_j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m) {
+    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
+    if (lcs.m_r_heading[x_j] >= 0) // the basic var 
+        return false;
+
+    impq const & x_j_val = lcs.m_r_x[x_j];
+    linear_combination_iterator<mpq> *it = get_column_iterator(x_j);
+
+    inf_l = true;
+    inf_u = true;
+    l = u = zero_of_type<impq>();
+    m = mpq(1);
+
+    if (lower(x_j)) {
+        set_lower(l, inf_l, lower_bound(x_j));
+    }
+    if (upper(x_j)) {
+        set_upper(u, inf_u, upper_bound(x_j));
+    }
+
+    mpq a_ij; unsigned i;
+    while (it->next(a_ij, i)) {
+        unsigned x_i = lcs.m_r_basis[i];
+        impq const & x_i_val = lcs.m_r_x[x_i];
+        if (is_int(x_i) && is_int(x_j) && !a_ij.is_int())
+            m = lcm(m, denominator(a_ij));
+        bool x_i_lower = lower(x_i);
+        bool x_i_upper = upper(x_i);
+        if (a_ij.is_neg()) {
+            if (x_i_lower) {
+                impq new_l = x_j_val + ((x_i_val - lcs.m_r_low_bounds()[x_i]) / a_ij);
+                set_lower(l, inf_l, new_l);
+                if (!inf_l && !inf_u && l == u) break;;                
+            }
+            if (x_i_upper) {
+                impq new_u = x_j_val + ((x_i_val - lcs.m_r_upper_bounds()[x_i]) / a_ij);
+                set_upper(u, inf_u, new_u);
+                if (!inf_l && !inf_u && l == u) break;;                
+            }
+        }
+        else {
+            if (x_i_upper) {
+                impq new_l = x_j_val + ((x_i_val - lcs.m_r_upper_bounds()[x_i]) / a_ij);
+                set_lower(l, inf_u, new_l);
+                if (!inf_l && !inf_u && l == u) break;;                
+            }
+            if (x_i_lower) {
+                impq new_u = x_j_val + ((x_i_val - lcs.m_r_low_bounds()[x_i]) / a_ij);
+                set_upper(u, inf_u, new_u);
+                if (!inf_l && !inf_u && l == u) break;;                
+            }
+        }
+    }
+
+    delete it;
+    TRACE("freedom_interval",
+          tout << "freedom variable for:\n";
+          tout << m_lar_solver->get_column_name(x_j);
+          tout << "[";
+          if (inf_l) tout << "-oo"; else tout << l;
+          tout << "; ";
+          if (inf_u) tout << "oo";  else tout << u;
+          tout << "]\n";);
+    return true;
+
+}
+
+bool int_solver::is_int(unsigned j) const {
+    return m_lar_solver->column_is_int(j);
+}
+
+bool int_solver::value_is_int(unsigned j) const {
+    return m_lar_solver->m_mpq_lar_core_solver.m_r_x[j].is_int();
+}
+
+    
+
+bool int_solver::is_feasible() const {
+    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
+    lean_assert(
+                lcs.m_r_solver.calc_current_x_is_feasible_include_non_basis() ==
+                lcs.m_r_solver.current_x_is_feasible());
+    return lcs.m_r_solver.current_x_is_feasible();
+}
+const impq & int_solver::get_value(unsigned j) const {
+    return m_lar_solver->m_mpq_lar_core_solver.m_r_x[j];
+}
+
+void int_solver::display_column(std::ostream & out, unsigned j) const {
+    m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, out);
+}
+
+bool int_solver::inf_int_set_is_correct() const {
+    for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) {
+        if (m_inf_int_set.contains(j) != is_int(j) && (!value_is_int(j)))
+            return false;
+    }
+    return true;
+}
+
+bool int_solver::column_is_int_inf(unsigned j) const {
+    return is_int(j) && (!value_is_int(j));
+}
+    
+void int_solver::init_inf_int_set() {
+    m_inf_int_set.clear();
+    m_inf_int_set.resize(m_lar_solver->A_r().column_count());
+    for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_basis) {
+        if (column_is_int_inf(j))
+            m_inf_int_set.insert(j);
+    }
+}
+
+void int_solver::update_column_in_inf_set_set(unsigned j) {
+    if (is_int(j) && (!value_is_int(j)))
+        m_inf_int_set.insert(j);
+    else
+        m_inf_int_set.erase(j);
+}
+
+bool int_solver::is_base(unsigned j) const {
+    return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j] >= 0;
+}
+
+bool int_solver::is_boxed(unsigned j) const {
+    return m_lar_solver->m_mpq_lar_core_solver.m_column_types[j] == column_type::boxed;
+}
+
+lp_settings& int_solver::settings() {
+    return m_lar_solver->settings();
+}
+
+}

src/util/lp/int_solver.h

@@ -0,0 +1,100 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+#pragma once
+#include "util/lp/lp_settings.h"
+#include "util/lp/static_matrix.h"
+#include "util/lp/iterator_on_row.h"
+#include "util/lp/int_set.h"
+#include "util/lp/lar_term.h"
+
+namespace lean {
+class lar_solver;
+template <typename T, typename X>
+struct lp_constraint;
+enum class lia_move {
+        ok,
+        branch,
+        cut,
+        conflict,
+        give_up
+};
+
+struct explanation {
+    vector<std::pair<mpq, constraint_index>> m_explanation;
+};
+
+class int_solver {
+public:
+    // fields
+    lar_solver *m_lar_solver;
+    int_set m_old_values_set;
+    vector<impq> m_old_values_data;
+    int_set m_inf_int_set;
+    unsigned m_branch_cut_counter;
+    // methods
+    int_solver(lar_solver* lp);
+    // main function to check that solution provided by lar_solver is valid for integral values,
+    // or provide a way of how it can be adjusted.
+    lia_move check(lar_term& t, mpq& k, explanation& ex);
+private:
+
+    // how to tighten bounds for integer variables.
+
+    bool gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, unsigned i, explanation &); 
+    
+    // gcd test
+    // 5*x + 3*y + 6*z = 5
+    // suppose x is fixed at 2.
+    // so we have 10 + 3(y + 2z) = 5
+    //             5 = -3(y + 2z)
+    // this is unsolvable because 5/3 is not an integer.
+    // so we create a lemma that rules out this condition.
+    // 
+    bool gcd_test(explanation & ); // returns false in case of failure. Creates a theory lemma in case of failure.
+
+    // create goromy cuts
+    // either creates a conflict or a bound.
+
+    // branch and bound: 
+    // decide what to branch and bound on
+    // creates a fresh inequality.
+
+    bool branch(const lp_constraint<mpq, mpq> & new_inequality);
+    bool ext_gcd_test(iterator_on_row<mpq> & it,
+                      mpq const & least_coeff, 
+                      mpq const & lcm_den,
+                      mpq const & consts,
+                      explanation & ex);
+    void fill_explanation_from_fixed_columns(iterator_on_row<mpq> & it, explanation &);
+    void add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation &);
+    void remove_fixed_vars_from_base();
+    void patch_int_infeasible_columns();
+    bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m);
+    linear_combination_iterator<mpq> * get_column_iterator(unsigned j);
+    bool lower(unsigned j) const;
+    bool upper(unsigned j) const;
+    const impq & lower_bound(unsigned j) const;
+    const impq & upper_bound(unsigned j) const;
+    bool is_int(unsigned j) const;
+    bool is_base(unsigned j) const;
+    bool is_boxed(unsigned j) const;
+    bool value_is_int(unsigned j) const;
+    void set_value(unsigned j, const impq & new_val);
+    void fix_non_base_columns();
+    void failed();
+    bool is_feasible() const;
+    const impq & get_value(unsigned j) const;
+    void display_column(std::ostream & out, unsigned j) const;
+    bool inf_int_set_is_correct() const;
+    void init_inf_int_set();
+    void update_column_in_inf_set_set(unsigned j);
+    bool column_is_int_inf(unsigned j) const;
+    void trace_inf_rows() const;
+    int find_inf_int_base_column();
+    int find_inf_int_boxed_base_column_with_smallest_range();
+    lp_settings& settings();
+    void move_non_base_vars_to_bounds();
+};
+}

src/util/lp/lar_core_solver.h

@@ -811,6 +811,37 @@ public:
             return new iterator_on_indexed_vector<mpq>(m_r_solver.m_ed);
         }
     }
+
+    bool column_is_fixed(unsigned j) const {
+        return m_column_types()[j] == column_type::fixed ||
+            ( m_column_types()[j] == column_type::boxed &&
+              m_r_solver.m_low_bounds[j] == m_r_solver.m_upper_bounds[j]);
+    }
+
+    const impq & low_bound(unsigned j) const {
+        lean_assert(m_column_types()[j] == column_type::fixed ||
+                    m_column_types()[j] == column_type::boxed ||
+                    m_column_types()[j] == column_type::low_bound);
+        return m_r_low_bounds[j];
+    }
+
+    const impq & upper_bound(unsigned j) const {
+        lean_assert(m_column_types()[j] == column_type::fixed ||
+                    m_column_types()[j] == column_type::boxed ||
+                    m_column_types()[j] == column_type::upper_bound);
+        return m_r_upper_bounds[j];
+    }
+
     
+    const bool column_is_bounded(unsigned j) const {
+        switch(m_column_types()[j]) {
+        case column_type::fixed:
+        case column_type::boxed:
+            return true;
+        default:
+            return false;
+        }
+    }
+
 };
 }

src/util/lp/lar_solver.h

@@ -45,19 +45,31 @@ Revision History:
 #include "util/lp/iterator_on_row.h"
 #include "util/lp/quick_xplain.h"
 #include "util/lp/conversion_helper.h"
+
 namespace lp {
 
+
 class lar_solver : public column_namer {
+    class ext_var_info {
+        unsigned m_ext_j; // the external index
+        bool m_is_integer;
+    public:
+        ext_var_info(unsigned j): ext_var_info(j, false) {}
+        ext_var_info(unsigned j , bool is_int) : m_ext_j(j), m_is_integer(is_int) {}
+        unsigned ext_j() const { return m_ext_j;}
+        bool is_integer() const {return m_is_integer;}
+    };
     //////////////////// fields //////////////////////////
-    lp_settings                             m_settings;
-    stacked_value<lp_status>                m_status;
-    stacked_value<simplex_strategy_enum>    m_simplex_strategy;
-    std::unordered_map<unsigned, var_index> m_ext_vars_to_columns;
-    vector<unsigned>                        m_columns_to_ext_vars_or_term_indices;
-    stacked_vector<ul_pair>                 m_vars_to_ul_pairs;
-    vector<lar_base_constraint*>            m_constraints;
-    stacked_value<unsigned>                 m_constraint_count;
+    lp_settings m_settings;
+    stacked_value<lp_status> m_status;
+    stacked_value<simplex_strategy_enum> m_simplex_strategy;
+    std::unordered_map<unsigned, ext_var_info> m_ext_vars_to_columns;
+    vector<unsigned> m_columns_to_ext_vars_or_term_indices;
+    stacked_vector<ul_pair> m_vars_to_ul_pairs;
+    vector<lar_base_constraint*> m_constraints;
+    stacked_value<unsigned> m_constraint_count;
     // the set of column indices j such that bounds have changed for j
+
     int_set                                 m_columns_with_changed_bound;
     int_set                                 m_rows_with_changed_bounds;
     int_set                                 m_basic_columns_with_changed_cost;
@@ -69,12 +81,8 @@ class lar_solver : public column_namer {
     indexed_vector<mpq>                     m_column_buffer;
 public:
     lar_core_solver m_mpq_lar_core_solver;
-    unsigned constraint_count() const {
-        return m_constraints.size();
-    }
-    const lar_base_constraint& get_constraint(unsigned ci) const {
-        return *(m_constraints[ci]);
-    }
+    unsigned constraint_count() const;
+    const lar_base_constraint& get_constraint(unsigned ci) const;
 
     ////////////////// methods ////////////////////////////////
     static_matrix<mpq, numeric_pair<mpq>> & A_r() { return m_mpq_lar_core_solver.m_r_A;}
@@ -84,15 +92,17 @@ public:
 
     static bool valid_index(unsigned j){ return static_cast<int>(j) >= 0;}
 
-
+    bool column_is_int(unsigned j) const;
+    bool column_is_fixed(unsigned j) const;
 public:
-    lp_settings & settings() { return m_settings;}
 
-    lp_settings const & settings() const { return m_settings;}
+    // init region
+    bool strategy_is_undecided() const;
 
     void clear() {SASSERT(false); // not implemented
     }
 
+    void register_new_ext_var_index(unsigned ext_v, bool is_int);
 
     lar_solver() : m_status(OPTIMAL),
                    m_infeasible_column_index(-1),
@@ -104,28 +114,20 @@ public:
         m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows = v? (& m_rows_with_changed_bounds) : nullptr;
     }
 
-    virtual ~lar_solver(){
-        for (auto c : m_constraints)
-            delete c;
-        for (auto t : m_terms)
-            delete t;
-        for (auto t : m_orig_terms)
-            delete t;
-    }
+    bool var_is_int(var_index v) const;
 
 #include "util/lp/init_lar_solver.h"
 
     numeric_pair<mpq> const& get_value(var_index vi) const { return m_mpq_lar_core_solver.m_r_x[vi]; }
 
-    bool is_term(var_index j) const {
-        return j >= m_terms_start_index && j - m_terms_start_index < m_terms.size();
-    }
+    void add_new_var_to_core_fields_for_doubles(bool register_in_basis);
 
     unsigned adjust_term_index(unsigned j) const {
         SASSERT(is_term(j));
         return j - m_terms_start_index;
     }
 
+    void decide_on_strategy_and_adjust_initial_state();
 
     bool use_lu() const { return m_settings.simplex_strategy() == simplex_strategy_enum::lu; }
 
@@ -244,8 +246,6 @@ public:
                                            row_index,
                                            bp
                                            );
-    }
-
 
     void substitute_basis_var_in_terms_for_row(unsigned i) {
         // todo : create a map from term basic vars to the rows where they are used
@@ -428,11 +428,9 @@ public:
 
     lp_status get_status() const { return m_status;}
 
-    void set_status(lp_status s) {m_status = s;}
+    lp_status get_status() const;
 
-    lp_status find_feasible_solution() {
-        if (strategy_is_undecided())
-            decide_on_strategy_and_adjust_initial_state();
+    void set_status(lp_status s);
 
         m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true;
         return solve();
@@ -454,50 +452,19 @@ public:
 
     void fill_explanation_from_infeasible_column(vector<std::pair<mpq, constraint_index>> & evidence) const{
 
-        // this is the case when the lower bound is in conflict with the upper one
-        const ul_pair & ul =  m_vars_to_ul_pairs[m_infeasible_column_index];
-        evidence.push_back(std::make_pair(numeric_traits<mpq>::one(), ul.upper_bound_witness()));
-        evidence.push_back(std::make_pair(-numeric_traits<mpq>::one(), ul.low_bound_witness()));
-    }
-
+    void fill_explanation_from_infeasible_column(explanation_t & evidence) const;
 
     unsigned get_total_iterations() const { return m_mpq_lar_core_solver.m_r_solver.total_iterations(); }
+    
     // see http://research.microsoft.com/projects/z3/smt07.pdf
     // This method searches for a feasible solution with as many different values of variables, reverenced in vars, as it can find
     // Attention, after a call to this method the non-basic variables don't necesserarly stick to their bounds anymore
-    vector<unsigned> get_list_of_all_var_indices() const {
-        vector<unsigned> ret;
-        for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_heading.size(); j++)
-            ret.push_back(j);
-        return ret;
-    }
-    void push() {
-        m_simplex_strategy = m_settings.simplex_strategy();
-        m_simplex_strategy.push();
-        m_status.push();
-        m_vars_to_ul_pairs.push();
-        m_infeasible_column_index.push();
-        m_mpq_lar_core_solver.push();
-        m_term_count = m_terms.size();
-        m_term_count.push();
-        m_constraint_count = m_constraints.size();
-        m_constraint_count.push();
-    }
+    vector<unsigned> get_list_of_all_var_indices() const;
+    void push();
 
-    static void clean_large_elements_after_pop(unsigned n, int_set& set) {
-        vector<int> to_remove;
-        for (unsigned j: set.m_index)
-            if (j >= n)
-                to_remove.push_back(j);
-        for (unsigned j : to_remove)
-            set.erase(j);
-    }
-
-    static void shrink_inf_set_after_pop(unsigned n, int_set & set) {
-        clean_large_elements_after_pop(n, set);
-        set.resize(n);
-    }
+    static void clean_large_elements_after_pop(unsigned n, int_set& set);
 
+    static void shrink_inf_set_after_pop(unsigned n, int_set & set);
 
     void pop(unsigned k) {
         int n_was = static_cast<int>(m_ext_vars_to_columns.size());
@@ -669,27 +636,14 @@ public:
         return *m_terms[j - m_terms_start_index];
     }
 
-    void pop_core_solver_params() {
-        pop_core_solver_params(1);
-    }
+    void pop_core_solver_params();
 
-    void pop_core_solver_params(unsigned k) {
-        A_r().pop(k);
-        A_d().pop(k);
-    }
+    void pop_core_solver_params(unsigned k);
 
 
-    void set_upper_bound_witness(var_index j, constraint_index ci) {
-        ul_pair ul = m_vars_to_ul_pairs[j];
-        ul.upper_bound_witness() = ci;
-        m_vars_to_ul_pairs[j] = ul;
-    }
+    void set_upper_bound_witness(var_index j, constraint_index ci);
 
-    void set_low_bound_witness(var_index j, constraint_index ci) {
-        ul_pair ul = m_vars_to_ul_pairs[j];
-        ul.low_bound_witness() = ci;
-        m_vars_to_ul_pairs[j] = ul;
-    }
+    void set_low_bound_witness(var_index j, constraint_index ci);
 
 
     void substitute_terms(const mpq & mult,
@@ -727,7 +681,7 @@ public:
             m_rows_with_changed_bounds.insert(rc.m_i);
     }
 
-    bool use_tableau() const { return m_settings.use_tableau(); }
+    bool use_tableau() const;
 
     bool use_tableau_costs() const {
         return m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs;
@@ -765,9 +719,7 @@ public:
         return true;
     }
 
-    bool tableau_with_costs() const {
-        return m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs;
-    }
+    bool tableau_with_costs() const;
 
     bool costs_are_used() const {
         return m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows;
@@ -795,23 +747,7 @@ public:
         }
     }
 
-    void update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) {
-        if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) {
-            if (costs_are_used()) {
-                bool was_infeas = m_mpq_lar_core_solver.m_r_solver.m_inf_set.contains(j);
-                m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(j);
-                if (was_infeas != m_mpq_lar_core_solver.m_r_solver.m_inf_set.contains(j))
-                    m_basic_columns_with_changed_cost.insert(j);
-            } else {
-                m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(j);
-            }
-        } else {
-            numeric_pair<mpq> delta;
-            if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta))
-                change_basic_x_by_delta_on_column(j, delta);
-        }
-    }
-
+    void update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j);
 
     void detect_rows_with_changed_bounds_for_column(unsigned j) {
         if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) {
@@ -941,19 +877,9 @@ public:
 
     }
 
-    bool var_is_registered(var_index vj) const {
-        if (vj >= m_terms_start_index) {
-            if (vj - m_terms_start_index >= m_terms.size())
-                return false;
-        } else if ( vj >= A_r().column_count()) {
-            return false;
-        }
-        return true;
-    }
+    bool var_is_registered(var_index vj) const;
 
-    unsigned constraint_stack_size() const {
-        return m_constraint_count.stack_size();
-    }
+    unsigned constraint_stack_size() const;
 
     void fill_last_row_of_A_r(static_matrix<mpq, numeric_pair<mpq>> & A, const lar_term * ls) {
         SASSERT(A.row_count() > 0);
@@ -1080,32 +1006,24 @@ public:
         return false; // it is unreachable
     }
 
+    bool try_to_set_fixed(column_info<mpq> & ci);
+
+    column_type get_column_type(const column_info<mpq> & ci);
+
+    std::string get_column_name(unsigned j) const;
+
+    bool all_constrained_variables_are_registered(const vector<std::pair<mpq, var_index>>& left_side);
+
+    constraint_index add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm);
+    bool all_constraints_hold() const;
+    bool constraint_holds(const lar_base_constraint & constr, std::unordered_map<var_index, mpq> & var_map) const;
+    bool the_relations_are_of_same_type(const vector<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum) const;
+
+    static void register_in_map(std::unordered_map<var_index, mpq> & coeffs, const lar_base_constraint & cn, const mpq & a);
+    static void register_one_coeff_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j);
 
 
-
-
-
-
-    bool the_relations_are_of_same_type(const vector<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum) const {
-        unsigned n_of_G = 0, n_of_L = 0;
-        bool strict = false;
-        for (auto & it : evidence) {
-            mpq coeff = it.first;
-            constraint_index con_ind = it.second;
-            lconstraint_kind kind = coeff.is_pos() ?
-                m_constraints[con_ind]->m_kind :
-                flip_kind(m_constraints[con_ind]->m_kind);
-            if (kind == GT || kind == LT)
-                strict = true;
-            if (kind == GE || kind == GT) n_of_G++;
-            else if (kind == LE || kind == LT) n_of_L++;
-        }
-        the_kind_of_sum = n_of_G ? GE : (n_of_L ? LE : EQ);
-        if (strict)
-            the_kind_of_sum = static_cast<lconstraint_kind>((static_cast<int>(the_kind_of_sum) / 2));
-
-        return n_of_G == 0 || n_of_L == 0;
-    }
+    bool the_left_sides_sum_to_zero(const vector<std::pair<mpq, unsigned>> & evidence) const;
 
     static void register_in_map(std::unordered_map<var_index, mpq> & coeffs, const lar_base_constraint & cn, const mpq & a) {
         for (auto & it : cn.get_left_side_coefficients()) {
@@ -1258,11 +1176,8 @@ public:
     void get_infeasibility_explanation_for_inf_sign(
                                                     vector<std::pair<mpq, constraint_index>> & explanation,
                                                     const vector<std::pair<mpq, unsigned>> & inf_row,
-                                                    int inf_sign) const {
+                                                    int inf_sign) const;
 
-        for (auto & it : inf_row) {
-            mpq coeff = it.first;
-            unsigned j = it.second;
 
             int adj_sign = coeff.is_pos() ? inf_sign : -inf_sign;
             const ul_pair & ul = m_vars_to_ul_pairs[j];
@@ -1273,6 +1188,7 @@ public:
         }
     }
 
+    void get_model(std::unordered_map<var_index, mpq> & variable_values) const;
 
 
     void get_model(std::unordered_map<var_index, mpq> & variable_values) const {
@@ -1306,28 +1222,9 @@ public:
     }
 
     // ********** print region start
-    void print_constraint(constraint_index ci, std::ostream & out) const {
-        if (ci >= m_constraints.size()) {
-            out << "constraint " << T_to_string(ci) << " is not found";
-            out << std::endl;
-            return;
-        }
+    void print_constraint(constraint_index ci, std::ostream & out) const;
 
-        print_constraint(m_constraints[ci], out);
-    }
-
-    void print_constraints(std::ostream& out) const  {
-        for (auto c : m_constraints) {
-            print_constraint(c, out);
-        }
-    }
-
-    void print_terms(std::ostream& out) const  {
-        for (auto it : m_terms) {
-            print_term(*it, out);
-            out << "\n";
-        }
-    }
+    void print_constraints(std::ostream& out) const ;
 
     void print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const {
         print_linear_combination_of_column_indices(c->get_left_side_coefficients(), out);
@@ -1337,12 +1234,7 @@ public:
 
     }
 
-    void print_term(lar_term const& term, std::ostream & out) const {
-        if (!numeric_traits<mpq>::is_zero(term.m_v)) {
-            out << term.m_v << " + ";
-        }
-        print_linear_combination_of_column_indices(term.coeffs_as_vector(), out);
-    }
+    void print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const;
 
     mpq get_left_side_val(const lar_base_constraint &  cns, const std::unordered_map<var_index, mpq> & var_map) const {
         mpq ret = cns.get_free_coeff_of_left_side();
@@ -1355,10 +1247,7 @@ public:
         return ret;
     }
 
-    void print_constraint(const lar_base_constraint * c, std::ostream & out) const {
-        print_left_side_of_constraint(c, out);
-        out << " " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side << std::endl;
-    }
+    mpq get_left_side_val(const lar_base_constraint &  cns, const std::unordered_map<var_index, mpq> & var_map) const;
 
     void fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list) {
         for (unsigned i = 0; i < sz; i++) {
@@ -1373,22 +1262,21 @@ public:
         }
     }
 
-    void random_update(unsigned sz, var_index const * vars) {
-        vector<unsigned> column_list;
-        fill_var_set_for_random_update(sz, vars, column_list);
-        random_updater ru(m_mpq_lar_core_solver, column_list);
-        ru.update();
-    }
+   void fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list);
 
+    void random_update(unsigned sz, var_index const * vars);
+    void pivot_fixed_vars_from_basis();
+    void pop();
+    bool column_represents_row_in_tableau(unsigned j);
+    void make_sure_that_the_bottom_right_elem_not_zero_in_tableau(unsigned i, unsigned j);
+    void remove_last_row_and_column_from_tableau(unsigned j);
+    void remove_last_column_from_tableau(unsigned j);
 
-    void try_pivot_fixed_vars_from_basis() {
-        m_mpq_lar_core_solver.m_r_solver.pivot_fixed_vars_from_basis();
-    }
-
-    void pop() {
-        pop(1);
-    }
-
+    void remove_last_column_from_basis_tableau(unsigned j);
+    void remove_column_from_tableau(unsigned j);
+    void pop_tableau();
+    void clean_inf_set_of_r_solver_after_pop();
+    void shrink_explanation_to_minimum(vector<std::pair<mpq, constraint_index>> & explanation) const;
 
     bool column_represents_row_in_tableau(unsigned j) {
         return m_vars_to_ul_pairs()[j].m_i != static_cast<row_index>(-1);

src/util/lp/lar_solver_instances.cpp

@@ -0,0 +1,13 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#include "util/lp/lar_solver.cpp"
+
+template void  lean::lar_solver::copy_from_mpq_matrix<double,double>(class lean::static_matrix<double,double> &);
+
+
+
+
+

src/util/lp/lp_core_solver_base.h

@@ -444,6 +444,7 @@ public:
     void init_lu();
     int pivots_in_column_and_row_are_different(int entering, int leaving) const;
     void pivot_fixed_vars_from_basis();
+    bool pivot_column_general(unsigned j, unsigned j_basic, indexed_vector<T> & w);
     bool pivot_for_tableau_on_basis();
     bool pivot_row_for_tableau_on_basis(unsigned row);
     void init_basic_part_of_basis_heading() {
@@ -583,8 +584,8 @@ public:
         default:
             SASSERT(false);
         }
-        std::cout << "basis heading = " << m_basis_heading[j] << std::endl;
-        std::cout << "x = " << m_x[j] << std::endl;
+        out << "basis heading = " << m_basis_heading[j] << std::endl;
+        out << "x = " << m_x[j] << std::endl;
         /*
         std::cout << "cost = " << m_costs[j] << std::endl;
         std:: cout << "m_d = " << m_d[j] << std::endl;*/

src/util/lp/lp_core_solver_base.hpp

@@ -938,7 +938,27 @@ template <typename T, typename X>  void lp_core_solver_base<T, X>::transpose_row
     transpose_basis(i, j);
     m_A.transpose_rows(i, j);
 }
-
+// j is the new basic column, j_basic - the leaving column
+template <typename T, typename X> bool lp_core_solver_base<T, X>::pivot_column_general(unsigned j, unsigned j_basic, indexed_vector<T> & w)  {
+    unsigned row_index = m_basis_heading[j_basic];
+    change_basis(j, j_basic);
+    if (m_settings.m_simplex_strategy == simplex_strategy_enum::lu) {
+        if (m_factorization->need_to_refactor()) {
+            init_lu();
+        } else {
+            m_factorization->prepare_entering(j, w); // to init vector w
+            m_factorization->replace_column(zero_of_type<T>(), w, row_index);
+        }
+        if (m_factorization->get_status() != LU_status::OK) {
+            change_basis(j_basic, j);
+            init_lu();
+            return false;
+        }
+    } else { // the tableau case
+        pivot_column_tableau(j, row_index);
+    }
+    return true;
+}
 template <typename T, typename X>  void lp_core_solver_base<T, X>::pivot_fixed_vars_from_basis() {
     // run over basis and non-basis at the same time
     indexed_vector<T> w(m_basis.size()); // the buffer
@@ -958,23 +978,11 @@ template <typename T, typename X>  void lp_core_solver_base<T, X>::pivot_fixed_v
             if (j >= m_nbasis.size())
                 break;
             j++;
-            if (m_factorization->need_to_refactor()) {
-                change_basis(jj, ii);
-                init_lu();
-            } else {
-                m_factorization->prepare_entering(jj, w); // to init vector w
-                m_factorization->replace_column(zero_of_type<T>(), w, m_basis_heading[ii]);
-                change_basis(jj, ii);
-            }
-            if (m_factorization->get_status() != LU_status::OK) {
-                change_basis(ii, jj);
-                init_lu();
-            } else {
+            if (!pivot_column_general(jj, ii, w))
                 break;
             }
         }
         SASSERT(m_factorization->get_status()== LU_status::OK);
-    }
 }
 
 template <typename T, typename X> bool 

src/util/lp/lp_settings.h

@@ -31,13 +31,16 @@ namespace lp {
 typedef unsigned var_index;
 typedef unsigned constraint_index;
 typedef unsigned row_index;
+
+typedef vector<std::pair<mpq, constraint_index>> explanation_t;
+
 enum class column_type  {
     free_column = 0,
-        low_bound = 1,
-        upper_bound = 2,
-        boxed = 3,
-        fixed = 4
-        };
+    low_bound = 1,
+    upper_bound = 2,
+    boxed = 3,
+    fixed = 4
+};
 
 enum class simplex_strategy_enum {
     undecided = 3,
@@ -91,11 +94,14 @@ public:
 };
 
 struct stats {
+    unsigned m_make_feasible;
     unsigned m_total_iterations;
     unsigned m_iters_with_no_cost_growing;
     unsigned m_num_factorizations;
     unsigned m_num_of_implied_bounds;
     unsigned m_need_to_solve_inf;
+    unsigned m_max_cols;
+    unsigned m_max_rows;
     stats() { reset(); }
     void reset() { memset(this, 0, sizeof(*this)); }
 };
@@ -214,7 +220,8 @@ public:
                     use_breakpoints_in_feasibility_search(false),
                     max_row_length_for_bound_propagation(300),
                     backup_costs(true),
-                    column_number_threshold_for_using_lu_in_lar_solver(4000)
+                    column_number_threshold_for_using_lu_in_lar_solver(4000),
+                    m_int_branch_cut_threshold(10000000)
     {}
 
     void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; }
@@ -321,6 +328,7 @@ public:
     unsigned max_row_length_for_bound_propagation;
     bool backup_costs;
     unsigned column_number_threshold_for_using_lu_in_lar_solver;
+    unsigned m_int_branch_cut_threshold;
 }; // end of lp_settings class
 
 

src/util/lp/mps_reader.h

@@ -825,7 +825,7 @@ public:
             auto kind = get_lar_relation_from_row(row->m_type);
             vector<std::pair<mpq, var_index>> ls;
             for (auto s : row->m_row_columns) {
-                var_index i = solver->add_var(get_var_index(s.first));
+                var_index i = solver->add_var(get_var_index(s.first), false);
                 ls.push_back(std::make_pair(s.second, i));
             }
             solver->add_constraint(ls, kind, row->m_right_side);
@@ -843,20 +843,20 @@ public:
 
     void create_low_constraint_for_var(column* col, bound * b, lar_solver *solver) {
         vector<std::pair<mpq, var_index>> ls;
-        var_index i = solver->add_var(col->m_index);
+        var_index i = solver->add_var(col->m_index, false);
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
         solver->add_constraint(ls, GE, b->m_low);
     }
 
     void create_upper_constraint_for_var(column* col, bound * b, lar_solver *solver) {
-        var_index i = solver->add_var(col->m_index);
+        var_index i = solver->add_var(col->m_index, false);
         vector<std::pair<mpq, var_index>> ls;
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
         solver->add_constraint(ls, LE, b->m_upper);
     }
 
     void create_equality_contraint_for_var(column* col, bound * b, lar_solver *solver) {
-        var_index i = solver->add_var(col->m_index);
+        var_index i = solver->add_var(col->m_index, false);
         vector<std::pair<mpq, var_index>> ls;
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
         solver->add_constraint(ls, EQ, b->m_fixed_value);
@@ -865,7 +865,7 @@ public:
     void fill_lar_solver_on_columns(lar_solver * solver) {
         for (auto s : m_columns) {
             mps_reader::column * col = s.second;
-            solver->add_var(col->m_index);
+            solver->add_var(col->m_index, false);
             auto b = col->m_bound;
             if (b == nullptr) return;
 

src/util/lp/nra_solver.cpp

@@ -0,0 +1,264 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#include "util/lp/lar_solver.h"
+#include "util/lp/nra_solver.h"
+#include "nlsat/nlsat_solver.h"
+#include "math/polynomial/polynomial.h"
+#include "math/polynomial/algebraic_numbers.h"
+#include "util/map.h"
+
+
+namespace nra {
+
+    struct mon_eq {
+        mon_eq(lean::var_index v, unsigned sz, lean::var_index const* vs):
+            m_v(v), m_vs(sz, vs) {}
+        lean::var_index          m_v;
+        svector<lean::var_index> m_vs;
+    };
+
+    struct solver::imp {
+        lean::lar_solver&      s;
+        reslimit&              m_limit;  
+        params_ref             m_params; 
+        u_map<polynomial::var> m_lp2nl;  // map from lar_solver variables to nlsat::solver variables
+        scoped_ptr<nlsat::solver>  m_nlsat;
+        vector<mon_eq>         m_monomials;
+        unsigned_vector        m_monomials_lim;
+        mutable std::unordered_map<lean::var_index, rational> m_variable_values; // current model
+
+        imp(lean::lar_solver& s, reslimit& lim, params_ref const& p): 
+            s(s), 
+            m_limit(lim),
+            m_params(p) {
+        }
+
+        bool need_check() {
+            return !m_monomials.empty() && !check_assignments();
+        }
+
+        void add(lean::var_index v, unsigned sz, lean::var_index const* vs) {
+            m_monomials.push_back(mon_eq(v, sz, vs));
+        }
+
+        void push() {
+            m_monomials_lim.push_back(m_monomials.size());
+        }
+
+        void pop(unsigned n) {
+            if (n == 0) return;
+            m_monomials.shrink(m_monomials_lim[m_monomials_lim.size() - n]);
+            m_monomials_lim.shrink(m_monomials_lim.size() - n);       
+        }
+
+        /*
+          \brief Check if polynomials are well defined.
+          multiply values for vs and check if they are equal to value for v.
+          epsilon has been computed.
+        */
+        bool check_assignment(mon_eq const& m) const {
+            rational r1 = m_variable_values[m.m_v];
+            rational r2(1);
+            for (auto w : m.m_vs) {
+                r2 *= m_variable_values[w];
+            }
+            return r1 == r2;
+        }
+
+        bool check_assignments() const {
+            s.get_model(m_variable_values);
+            for (auto const& m : m_monomials) {
+                if (!check_assignment(m)) return false;
+            }
+            return true;
+        }
+
+        /**
+           \brief one-shot nlsat check.
+           A one shot checker is the least functionality that can 
+           enable non-linear reasoning. 
+           In addition to checking satisfiability we would also need
+           to identify equalities in the model that should be assumed
+           with the remaining solver.
+           
+           TBD: use partial model from lra_solver to prime the state of nlsat_solver.
+           TBD: explore more incremental ways of applying nlsat (using assumptions)
+        */
+        lbool check(lean::explanation_t& ex) {
+            SASSERT(need_check());
+            m_nlsat = alloc(nlsat::solver, m_limit, m_params);
+            m_lp2nl.reset();
+            vector<nlsat::assumption, false> core;
+
+            // add linear inequalities from lra_solver
+            for (unsigned i = 0; i < s.constraint_count(); ++i) {
+                add_constraint(i);
+            }
+
+            // add polynomial definitions.
+            for (auto const& m : m_monomials) {
+                add_monomial_eq(m);
+            }
+            // TBD: add variable bounds?
+
+            lbool r = m_nlsat->check(); 
+            TRACE("arith", m_nlsat->display(tout << r << "\n"););
+            switch (r) {
+            case l_true: 
+                break;
+            case l_false: 
+                ex.reset();
+                m_nlsat->get_core(core);
+                for (auto c : core) {
+                    unsigned idx = static_cast<unsigned>(static_cast<imp*>(c) - this);
+                    ex.push_back(std::pair<rational, unsigned>(rational(1), idx));
+                    TRACE("arith", tout << "ex: " << idx << "\n";);
+                }
+                break;
+
+            case l_undef:
+                break;
+            }            
+            return r;
+        }                
+
+        void add_monomial_eq(mon_eq const& m) {
+            polynomial::manager& pm = m_nlsat->pm();
+            svector<polynomial::var> vars;
+            for (auto v : m.m_vs) {
+                vars.push_back(lp2nl(v));
+            }
+            polynomial::monomial_ref m1(pm.mk_monomial(vars.size(), vars.c_ptr()), pm);
+            polynomial::monomial_ref m2(pm.mk_monomial(lp2nl(m.m_v), 1), pm);
+            polynomial::monomial* mls[2] = { m1, m2 };
+            polynomial::scoped_numeral_vector coeffs(pm.m());
+            coeffs.push_back(mpz(1));
+            coeffs.push_back(mpz(-1));
+            polynomial::polynomial_ref p(pm.mk_polynomial(2, coeffs.c_ptr(), mls),  pm);
+            polynomial::polynomial* ps[1] = { p };
+            bool even[1] = { false };
+            nlsat::literal lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, even);
+            m_nlsat->mk_clause(1, &lit, 0);
+        }
+
+        void add_constraint(unsigned idx) {
+            auto& c = s.get_constraint(idx);
+            auto& pm = m_nlsat->pm();
+            auto k = c.m_kind;
+            auto rhs = c.m_right_side;
+            auto lhs = c.get_left_side_coefficients();
+            auto sz = lhs.size();
+            svector<polynomial::var> vars;
+            rational den = denominator(rhs);
+            for (auto kv : lhs) {
+                vars.push_back(lp2nl(kv.second));
+                den = lcm(den, denominator(kv.first));
+            }
+            vector<rational> coeffs;
+            for (auto kv : lhs) {
+                coeffs.push_back(den * kv.first);
+            }
+            rhs *= den;
+            polynomial::polynomial_ref p(pm.mk_linear(sz, coeffs.c_ptr(), vars.c_ptr(), -rhs), pm);
+            polynomial::polynomial* ps[1] = { p };
+            bool is_even[1] = { false };
+            nlsat::literal lit;
+            nlsat::assumption a = this + idx;
+            switch (k) {
+            case lean::lconstraint_kind::LE:
+                lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::GE:
+                lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::LT:
+                lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::GT:
+                lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::EQ:
+                lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even);                
+                break;
+            }
+            m_nlsat->mk_clause(1, &lit, a);
+        }               
+
+        bool is_int(lean::var_index v) {
+            return s.var_is_int(v);
+        }
+
+
+        polynomial::var lp2nl(lean::var_index v) {
+            polynomial::var r;
+            if (!m_lp2nl.find(v, r)) {
+                r = m_nlsat->mk_var(is_int(v));
+                m_lp2nl.insert(v, r);
+            }
+            return r;
+        }
+
+        nlsat::anum const& value(lean::var_index v) const {
+           return m_nlsat->value(m_lp2nl.find(v));
+        }
+
+        nlsat::anum_manager& am() {
+            return m_nlsat->am();
+        }
+
+        std::ostream& display(std::ostream& out) const {
+            for (auto m : m_monomials) {
+                out << "v" << m.m_v << " = ";
+                for (auto v : m.m_vs) {
+                    out << "v" << v << " ";
+                }
+                out << "\n";
+            }
+            return out;
+        }
+    };
+
+    solver::solver(lean::lar_solver& s, reslimit& lim, params_ref const& p) {
+        m_imp = alloc(imp, s, lim, p);
+    }
+
+    solver::~solver() {
+        dealloc(m_imp);
+    }
+
+    void solver::add_monomial(lean::var_index v, unsigned sz, lean::var_index const* vs) {
+        m_imp->add(v, sz, vs);
+    }
+
+    lbool solver::check(lean::explanation_t& ex) {
+        return m_imp->check(ex);
+    }
+
+    bool solver::need_check() {
+        return m_imp->need_check();
+    }
+
+    void solver::push() {
+        m_imp->push();
+    }
+
+    void solver::pop(unsigned n) {
+        m_imp->pop(n);
+    }
+
+    std::ostream& solver::display(std::ostream& out) const {
+        return m_imp->display(out);
+    }
+
+    nlsat::anum const& solver::value(lean::var_index v) const {
+        return m_imp->value(v);
+    }
+
+    nlsat::anum_manager& solver::am() {
+        return m_imp->am();
+    }
+
+}

src/util/lp/nra_solver.h

@@ -0,0 +1,70 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#pragma once
+#include "util/vector.h"
+#include "util/lp/lp_settings.h"
+#include "util/rlimit.h"
+#include "util/params.h"
+#include "nlsat/nlsat_solver.h"
+
+namespace lean {
+    class lar_solver;
+}
+
+
+namespace nra {
+
+
+
+    class solver {
+        struct imp;
+        imp* m_imp;
+
+    public:
+
+        solver(lean::lar_solver& s, reslimit& lim, params_ref const& p = params_ref());
+        
+        ~solver();
+
+        /*
+          \brief Add a definition v = vs[0]*vs[1]*...*vs[sz-1]
+          The variable v is equal to the product of variables vs.
+        */
+        void add_monomial(lean::var_index v, unsigned sz, lean::var_index const* vs);
+
+        /*
+          \brief Check feasiblity of linear constraints augmented by polynomial definitions
+          that are added.
+         */
+        lbool check(lean::explanation_t& ex);
+
+        /*
+          \brief determine whether nra check is needed.
+        */
+        bool need_check();
+
+        /*
+          \brief Access model.
+        */
+        nlsat::anum const& value(lean::var_index v) const;
+
+        nlsat::anum_manager& am();        
+
+        /*
+          \brief push and pop scope. 
+          Monomial definitions are retraced when popping scope.
+        */
+        void push();
+        
+        void pop(unsigned n);
+
+        /*
+          \brief display state
+         */
+        std::ostream& display(std::ostream& out) const;
+
+    };
+}

src/util/lp/numeric_pair.h

@@ -203,6 +203,11 @@ struct numeric_pair {
     std::string to_string() const {
         return std::string("(") + T_to_string(x) + ", "  + T_to_string(y) + ")";
     }
+
+    bool is_int() const {
+        return x.is_int() && y.is_zero();
+    }
+    
 };
 
 
@@ -328,4 +333,26 @@ struct convert_struct<double, double> {
 template <typename X> bool is_epsilon_small(const X & v, const double &eps) { return convert_struct<X, double>::is_epsilon_small(v, eps);}
 template <typename X> bool below_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::below_bound_numeric(x, bound, eps);}
 template <typename X> bool above_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::above_bound_numeric(x, bound, eps);}
+template  <typename T>  T floor(const numeric_pair<T> & r) {
+    if (r.x.is_int()) {
+        if (r.y.is_nonneg()) {
+            return r.x;
+        }
+        return r.x - mpq::one();
+    }
+    
+    return floor(r.x);
+}
+
+template <typename T> T ceil(const numeric_pair<T> & r) {
+    if (r.x.is_int()) {
+        if (r.y.is_nonpos()) {
+            return r.x;
+        }
+        return r.x + mpq::one();
+    }
+    
+    return ceil(r.x);
+}
+
 }

src/util/lp/quick_xplain.cpp

@@ -35,7 +35,7 @@ void quick_xplain::copy_constraint_and_add_constraint_vars(const lar_constraint&
     vector < std::pair<mpq, unsigned>> ls;
     for (auto & p : lar_c.get_left_side_coefficients()) {
         unsigned j = p.second;
-        unsigned lj = m_qsol.add_var(j);
+        unsigned lj = m_qsol.add_var(j, false);
         ls.push_back(std::make_pair(p.first, lj));
     }
     m_constraints_in_local_vars.push_back(lar_constraint(ls, lar_c.m_kind, lar_c.m_right_side));
@@ -109,7 +109,7 @@ bool quick_xplain::is_feasible(const vector<unsigned> & x, unsigned k) const {
         vector < std::pair<mpq, unsigned>> ls;
         const lar_constraint & c = m_constraints_in_local_vars[i];
         for (auto & p : c.get_left_side_coefficients()) {
-            unsigned lj = l.add_var(p.second);
+            unsigned lj = l.add_var(p.second, false);
             ls.push_back(std::make_pair(p.first, lj));
         }
         l.add_constraint(ls, c.m_kind, c.m_right_side);

src/util/lp/stacked_vector.h

@@ -48,7 +48,10 @@ public:
         operator const B&() const {
             return m_vec.m_vector[m_i];
         }
-
+        
+        bool operator==(B const& other) const {
+            return m_vec.m_vector[m_i] == other;
+        }
     };
 
     class ref_const {

src/util/lp/static_matrix_instances.cpp

@@ -17,7 +17,13 @@ Revision History:
 
 
 --*/
+=======
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
 #include <memory>
+#include "util/vector.h"
 #include <set>
 #include <utility>
 #include "util/vector.h"