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>

contrib/cmake/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/CMakeLists.txt

@@ -35,10 +35,10 @@ endforeach()
 # raised if you try to declare a component is dependent on another component
 # that has not yet been declared.
 add_subdirectory(util)
-add_subdirectory(util/lp)
 add_subdirectory(math/polynomial)
 add_subdirectory(sat)
 add_subdirectory(nlsat)
+add_subdirectory(util/lp)
 add_subdirectory(math/hilbert)
 add_subdirectory(math/simplex)
 add_subdirectory(math/automata)

src/ast/ast.cpp

@@ -1727,7 +1727,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

@@ -738,9 +738,54 @@ br_status arith_rewriter::mk_idiv_core(expr * arg1, expr * arg2, expr_ref & resu
         result = m_util.mk_idiv0(arg1);
         return BR_REWRITE1;
     }
+    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;
 }
 
+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;
+    }
+    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;
+                }            
+            }
+        }
+    }
+    return false;
+}
+
 br_status arith_rewriter::mk_mod_core(expr * arg1, expr * arg2, expr_ref & result) {
     set_curr_sort(m().get_sort(arg1));
     numeral v1, v2;

src/ast/rewriter/arith_rewriter.h

@@ -90,6 +90,7 @@ class arith_rewriter : public poly_rewriter<arith_rewriter_core> {
     bool is_pi_integer_offset(expr * t, expr * & m);
     expr * mk_sin_value(rational const & k);
     app * mk_sqrt(rational const & k);
+    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.cpp

@@ -267,10 +267,55 @@ void arith_simplifier_plugin::mk_idiv(expr * arg1, expr * arg2, expr_ref & resul
     bool is_int;
     if (m_util.is_numeral(arg1, v1, is_int) && m_util.is_numeral(arg2, v2, is_int) && !v2.is_zero())
         result = m_util.mk_numeral(div(v1, v2), is_int);
+    else if (divides(arg2, arg1, result)) {
+        result = m_util.mk_mul(result, m_util.mk_idiv(arg2, arg2));
+    }
     else 
         result = m_util.mk_idiv(arg1, arg2);
 }
 
+bool arith_simplifier_plugin::divides(expr* d, expr* n, expr_ref& quot) {
+    ast_manager& m = m_manager;
+    if (d == n) {
+        quot = m_util.mk_numeral(rational(1), m_util.is_int(d));
+        return true;
+    }
+    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;
+                }            
+            }
+        }
+    }
+    return false;
+}
+
+
 void arith_simplifier_plugin::prop_mod_const(expr * e, unsigned depth, numeral const& k, expr_ref& result) {
     SASSERT(m_util.is_int(e));
     SASSERT(k.is_int() && k.is_pos());

src/ast/simplifier/arith_simplifier_plugin.h

@@ -48,6 +48,7 @@ protected:
 
     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);

src/math/polynomial/polynomial.h

@@ -19,16 +19,16 @@ Notes:
 #ifndef POLYNOMIAL_H_
 #define POLYNOMIAL_H_
 
-#include"mpz.h"
-#include"rational.h"
-#include"obj_ref.h"
-#include"ref_vector.h"
-#include"z3_exception.h"
-#include"scoped_numeral.h"
-#include"scoped_numeral_vector.h"
-#include"params.h"
-#include"mpbqi.h"
-#include"rlimit.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.cpp

@@ -623,6 +623,7 @@ namespace nlsat {
             unsigned sz = cs.size();
             for (unsigned i = 0; i < sz; i++)
                 del_clause(cs[i]);
+			cs.reset();
         }
 
         void del_clauses() {

src/nlsat/nlsat_solver.h

@@ -21,10 +21,10 @@ Revision History:
 #ifndef NLSAT_SOLVER_H_
 #define NLSAT_SOLVER_H_
 
-#include"nlsat_types.h"
-#include"params.h"
-#include"statistics.h"
-#include"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"polynomial.h"
-#include"buffer.h"
-#include"sat_types.h"
-#include"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,12 +19,12 @@ Revision History:
 #ifndef SAT_TYPES_H_
 #define SAT_TYPES_H_
 
-#include"debug.h"
-#include"approx_set.h"
-#include"lbool.h"
-#include"z3_exception.h"
-#include"common_msgs.h"
-#include"vector.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

@@ -80,8 +80,7 @@ void run_solver(lp_params & params, char const * mps_file_name) {
     solver->settings().set_message_ostream(&std::cout);
     solver->settings().report_frequency = params.rep_freq();
     solver->settings().print_statistics = params.print_stats();
-    solver->settings().simplex_strategy() = lean:: simplex_strategy_enum::lu;
-    
+    solver->settings().simplex_strategy() = lean::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/CMakeLists.txt

@@ -70,6 +70,7 @@ z3_add_component(smt
     euclid
     fpa
     grobner
+    nlsat
     lp
     macros
     normal_forms

src/smt/params/smt_params_helper.pyg

@@ -36,7 +36,7 @@ def_module_params(module_name='smt',
                           ('bv.reflect', BOOL, True, 'create enode for every bit-vector term'),
                           ('bv.enable_int2bv', BOOL, True, 'enable support for int2bv and bv2int operators'),
                           ('arith.random_initial_value', BOOL, False, 'use random initial values in the simplex-based procedure for linear arithmetic'),
-                          ('arith.solver', UINT, 2, 'arithmetic solver: 0 - no solver, 1 - bellman-ford based solver (diff. logic only), 2 - simplex based solver, 3 - floyd-warshall based solver (diff. logic only) and no theory combination'),
+                          ('arith.solver', UINT, 2, 'arithmetic solver: 0 - no solver, 1 - bellman-ford based solver (diff. logic only), 2 - simplex based solver, 3 - floyd-warshall based solver (diff. logic only) and no theory combination 4 - utvpi, 5 - infinitary lra, 6 - lra solver'),
                           ('arith.nl', BOOL, True, '(incomplete) nonlinear arithmetic support based on Groebner basis and interval propagation'),
                           ('arith.nl.gb', BOOL, True, 'groebner Basis computation, this option is ignored when arith.nl=false'),
                           ('arith.nl.branching', BOOL, True, 'branching on integer variables in non linear clusters'),

src/smt/params/theory_arith_params.h

@@ -23,12 +23,13 @@ Revision History:
 #include"params.h"
 
 enum arith_solver_id {
-    AS_NO_ARITH,
-    AS_DIFF_LOGIC,
-    AS_ARITH,
-    AS_DENSE_DIFF_LOGIC,
-    AS_UTVPI,
-    AS_OPTINF
+    AS_NO_ARITH,              // 0
+    AS_DIFF_LOGIC,            // 1
+    AS_ARITH,                 // 2
+    AS_DENSE_DIFF_LOGIC,      // 3
+    AS_UTVPI,                 // 4
+    AS_OPTINF,                // 5
+    AS_LRA                    // 6
 };
 
 enum bound_prop_mode {

src/smt/smt_model_generator.cpp

@@ -388,6 +388,7 @@ namespace smt {
             enode * n = *it3;
             if (is_uninterp_const(n->get_owner()) && m_context->is_relevant(n)) {
                 func_decl * d = n->get_owner()->get_decl();
+                TRACE("mg_top_sort", tout << d->get_name() << " " << (m_hidden_ufs.contains(d)?"hidden":"visible") << "\n";);
                 if (m_hidden_ufs.contains(d)) continue;
                 expr * val    = get_value(n);
                 m_model->register_decl(d, val);

src/smt/smt_setup.cpp

@@ -724,8 +724,6 @@ namespace smt {
     }
 
     void setup::setup_r_arith() {
-        // to disable theory lra
-        // m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params));        
         m_context.register_plugin(alloc(smt::theory_lra, m_manager, m_params));
     }
 
@@ -789,6 +787,9 @@ namespace smt {
         case AS_OPTINF:
             m_context.register_plugin(alloc(smt::theory_inf_arith, m_manager, m_params));            
             break;
+        case AS_LRA:
+            setup_r_arith();
+            break;
         default:
             if (m_params.m_arith_int_only && int_only)
                 m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params));

src/smt/theory_lra.cpp

@@ -35,7 +35,10 @@ Revision History:
 #include "smt/smt_model_generator.h"
 #include "smt/arith_eq_adapter.h"
 #include "util/nat_set.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 lp {
     enum bound_kind { lower_t, upper_t };
@@ -87,16 +90,12 @@ namespace 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));
@@ -118,9 +117,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;
@@ -144,10 +140,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;                     
@@ -198,14 +194,6 @@ namespace smt {
         };
        
         typedef vector<std::pair<rational, lean::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) {}
-        };
 
         svector<lean::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  
@@ -224,11 +212,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;
         expr*                  m_not_handled;
         ptr_vector<app>        m_underspecified;
         unsigned_vector        m_var_trail;
@@ -248,16 +232,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<lean::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;
 
@@ -289,14 +293,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) {
@@ -366,6 +381,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;
@@ -415,6 +438,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);
@@ -459,6 +520,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);                
@@ -487,7 +556,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);
@@ -549,14 +618,6 @@ 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, lean::EQ, -d.m_coeff), d.m_var);
-        }
         
         void internalize_eq(theory_var v1, theory_var v2) {                  
             enode* n1 = get_enode(v1);
@@ -650,15 +711,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(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(0),
-            m_resource_limit(*this) {
+            m_resource_limit(*this),
+            m_use_nra_model(false) {
         }
         
         ~imp() {
@@ -671,12 +731,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) {
@@ -710,54 +766,11 @@ namespace smt {
             //add_use_lists(b);
             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, *n2;
-            rational r;
-            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 = 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 = lp::lower_t;
-            }    
-            else {
-                TRACE("arith", tout << "Could not internalize " << mk_pp(atom, m) << "\n";);
-                found_not_handled(atom);
-                return true;
-            }
-            lp::bound* b = alloc(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);
             }
@@ -783,13 +796,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 {
@@ -808,13 +816,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) {
@@ -835,18 +841,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_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";);
         }
 
@@ -1080,7 +1084,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););
             
@@ -1130,54 +1136,69 @@ namespace smt {
                 enode* n2 = get_enode(v2);
                 m_assume_eq_head++;
                 CTRACE("arith", 
-                       get_ivalue(v1) == get_ivalue(v2) && n1->get_root() != n2->get_root(),
+                       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() != lean::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 != 0) {                    
-                    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;
@@ -1190,6 +1211,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
@@ -1259,14 +1344,13 @@ namespace smt {
 #else
                 propagate_bound(bv, is_true, b);
 #endif
-                if (!m_delay_constraints) {
-                    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;
             }
@@ -2133,18 +2217,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););
-            lean::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 lean::lp_status::INFEASIBLE:
                 return l_false;
@@ -2197,11 +2271,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;
@@ -2250,9 +2328,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) {
@@ -2278,6 +2390,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();
@@ -2290,6 +2403,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);
@@ -2301,6 +2415,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())));
@@ -2496,7 +2611,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);
@@ -2504,9 +2619,9 @@ 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);
         }        
     };
     

src/test/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_subdirectory(fuzzing)
+add_subdirectory(lp)
 ################################################################################
 # z3-test executable
 ################################################################################
@@ -117,7 +118,7 @@ add_executable(test-z3
   upolynomial.cpp
   var_subst.cpp
   vector.cpp
-  lp.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

@@ -2695,8 +2695,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));
@@ -2709,9 +2709,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) {
@@ -2723,8 +2730,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;
@@ -2758,9 +2765,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));
@@ -2823,8 +2830,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));
@@ -2837,8 +2844,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));
@@ -2853,8 +2860,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));
@@ -2868,9 +2875,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));
@@ -2892,9 +2899,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));
@@ -2918,9 +2925,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

@@ -376,7 +376,7 @@ namespace lean {
         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,8 +20,9 @@ z3_add_component(lp
     lp_solver_instances.cpp
     lu_instances.cpp
     matrix_instances.cpp
+    nra_solver.cpp
     permutation_matrix_instances.cpp
-	quick_xplain.cpp
+    quick_xplain.cpp
     row_eta_matrix_instances.cpp
     scaler_instances.cpp
     sparse_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/column_namer.h

@@ -15,7 +15,7 @@ public:
         T a;
         unsigned i;
         while (it->next(a, i)) {
-            coeff.emplace_back(a, i);
+            coeff.push_back(std::make_pair(a, i));
         }
         print_linear_combination_of_column_indices(coeff, out);
     }

src/util/lp/indexed_vector.h

@@ -75,16 +75,7 @@ public:
     }
 
     void set_value(const T& value, unsigned index);
-    void set_value_as_in_dictionary(unsigned index) {
-        lean_assert(index < m_data.size());
-        T & loc = m_data[index];
-        if (is_zero(loc)) {
-            m_index.push_back(index);
-            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,576 +0,0 @@
-/*
-  Copyright (c) 2017 Microsoft Corporation
-  Author: Lev Nachmanson
-*/
-
-// here we are inside lean::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;
-    lean_assert (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;
-    }
-    lean_assert(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);
-    lean_assert(sizes_are_correct());
-    return i;
-}
-
-void register_new_ext_var_index(unsigned ext_v) {
-    lean_assert(!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;
-    lean_assert(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();
-    lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
-    //        lean_assert(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);
-    lean_assert(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();
-    lean_assert(m_mpq_lar_core_solver.m_r_x.size() == j);
-    //        lean_assert(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);
-    lean_assert(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);
-    } 
-    lean_assert(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) {
-    lean_assert(sizes_are_correct());
-    add_row_from_term_no_constraint(term, term_ext_index);
-    lean_assert(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();
-    lean_assert(!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);
-    }
-    lean_assert(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:
-        lean_assert(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) {
-    lean_assert(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));
-    lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
-}
-
-void decide_on_strategy_and_adjust_initial_state() {
-	lean_assert(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:
-        lean_assert(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();
-    lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
-    //        lean_assert(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);
-    lean_assert(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) {
-    lean_assert(A.row_count() > 0);
-    lean_assert(A.column_count() > 0);
-    unsigned last_row = A.row_count() - 1;
-    lean_assert(A.m_rows[last_row].empty());
-
-    for (auto & t : ls->m_coeffs) {
-        lean_assert(!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;
-        lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
-        lean_assert(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;
-        lean_assert(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:
-        lean_unreachable();
-                
-    }
-    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) {
-    lean_assert(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:
-        lean_unreachable();
-                
-    }
-}
-    
-void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    lean_assert(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;
-                lean_assert(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:
-        lean_unreachable();
-                
-    }
-}
-void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    lean_assert(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:
-        lean_unreachable();
-                
-    }
-}
-
-void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
-    lean_assert(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]));
-    lean_assert(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:
-        lean_unreachable();
-                
-    }
-}
-    

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

@@ -796,6 +796,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_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

@@ -429,6 +429,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() {
@@ -568,8 +569,8 @@ public:
         default:
             lean_assert(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

@@ -923,7 +923,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
@@ -943,22 +963,9 @@ 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;
-            }
-        }
-        lean_assert(m_factorization->get_status()== LU_status::OK);
+        } 
     }
 }
 

src/util/lp/lp_settings.h

@@ -16,13 +16,16 @@ namespace lean {
 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,
@@ -75,11 +78,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)); }
 };
@@ -198,7 +204,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; }
@@ -278,13 +285,13 @@ public:
         return m_simplex_strategy;
     }
 
-	bool use_lu() const {
-		return m_simplex_strategy == simplex_strategy_enum::lu;
-	}
+    bool use_lu() const {
+        return m_simplex_strategy == simplex_strategy_enum::lu;
+    }
 
     bool use_tableau() const {
-		return m_simplex_strategy == simplex_strategy_enum::tableau_rows ||
-			m_simplex_strategy == simplex_strategy_enum::tableau_costs;
+        return m_simplex_strategy == simplex_strategy_enum::tableau_rows ||
+            m_simplex_strategy == simplex_strategy_enum::tableau_costs;
     }
 
     bool use_tableau_rows() const {
@@ -305,6 +312,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/lp_settings_instances.cpp

@@ -2,8 +2,8 @@
   Copyright (c) 2017 Microsoft Corporation
   Author: Lev Nachmanson
 */
-#include "util/vector.h"
 #include <memory>
+#include "util/vector.h"
 #include "util/lp/lp_settings.hpp"
 template bool lean::vectors_are_equal<double>(vector<double> const&, vector<double> const&);
 template bool lean::vectors_are_equal<lean::mpq>(vector<lean::mpq > const&, vector<lean::mpq> const&);

src/util/lp/mps_reader.h

@@ -810,7 +810,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);
@@ -828,20 +828,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);
@@ -850,7 +850,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

@@ -199,6 +199,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();
+    }
+    
 };
 
 
@@ -324,4 +329,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

@@ -20,7 +20,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));
@@ -94,7 +94,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

@@ -32,7 +32,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

@@ -2,8 +2,8 @@
   Copyright (c) 2017 Microsoft Corporation
   Author: Lev Nachmanson
 */
-#include "util/vector.h"
 #include <memory>
+#include "util/vector.h"
 #include <set>
 #include <utility>
 #include "util/lp/static_matrix.hpp"