avoid going creating hnf_cuts if all involved vars have integral values

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

add explanations to hnf cuts

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

nits and virtual methods (#68)

* local

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

* virtual method in bound propagator

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

cleanup from std::cout

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

handle the case when the number of terms is greater than the number of variables in hnf

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

method name's fix

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

restore hnf_cutter to work with m_row_count <= m_column_count

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

tune addition of rational numbers (#70)

* log quantifiers only if present

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

* merge and fix some warnings

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

* set new arith as default for LIA

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

* local

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

* virtual method in bound propagator

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

* prepare for mixed integer-real

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

* fix default tactic usage

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

give shorter explanations, call hnf only when have a not integral var

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

overhaul of mpq (#71)

* log quantifiers only if present

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

* merge and fix some warnings

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

* set new arith as default for LIA

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

* local

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

* virtual method in bound propagator

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

* prepare for mixed integer-real

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

* fix default tactic usage

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

* overhaul of mpz/mpq

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

* disabled temporary setting

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

* remove prints

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

fix for 32 bit build (#72)

* log quantifiers only if present

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

* merge and fix some warnings

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

* set new arith as default for LIA

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

* local

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

* virtual method in bound propagator

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

* prepare for mixed integer-real

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

* fix default tactic usage

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

* overhaul of mpz/mpq

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

* disabled temporary setting

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

* remove prints

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

* customize for 64 bit

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

yes (#74)

* log quantifiers only if present

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

* merge and fix some warnings

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

* set new arith as default for LIA

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

* local

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

* virtual method in bound propagator

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

* prepare for mixed integer-real

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

* fix default tactic usage

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

* overhaul of mpz/mpq

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

* disabled temporary setting

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

* remove prints

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

* customize for 64 bit

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

* customize for 64 bit

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

* more refactor

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

fix the merge

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

fixes in maximize_term untested

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

fix compilation (#75)

* log quantifiers only if present

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

* merge and fix some warnings

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

* set new arith as default for LIA

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

* local

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

* virtual method in bound propagator

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

* prepare for mixed integer-real

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

* fix default tactic usage

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

* overhaul of mpz/mpq

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

* disabled temporary setting

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

* remove prints

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

* customize for 64 bit

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

* customize for 64 bit

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

* more refactor

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

* merge

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

* relax check

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

* change for gcc

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

* merge

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

src/smt/params/smt_params_helper.pyg

@@ -40,7 +40,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 4 - utvpi, 5 - infinitary lra, 6 - lra solver'),
+                          ('arith.solver', UINT, 6, '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/smt_setup.cpp

@@ -818,10 +818,7 @@ namespace smt {
                 m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params));
             break;
         default:
-            if (m_params.m_arith_int_only && int_only)
+            setup_i_arith();
-                setup_i_arith();
-            else
-                m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params));
             break;
         }
     }

src/smt/theory_lra.cpp

@@ -751,12 +751,14 @@ public:
         init_solver();
     }
         
-    bool internalize_atom(app * atom, bool gate_ctx) {
+    void internalize_is_int(app * n) {
-        return internalize_atom_strict(atom, gate_ctx);
+        SASSERT(a.is_is_int(n));
-            
+        (void) mk_enode(n);
+        if (!ctx().relevancy())
+            mk_is_int_axiom(n);        
     }
 
-    bool internalize_atom_strict(app * atom, bool gate_ctx) {
+    bool internalize_atom(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());
@@ -772,6 +774,10 @@ public:
             v = internalize_def(to_app(n1));
             k = lp_api::lower_t;
         }    
+        else if (a.is_is_int(atom)) {
+            internalize_is_int(atom);
+            return true;
+        }
         else {
             TRACE("arith", tout << "Could not internalize " << mk_pp(atom, m) << "\n";);
             found_not_handled(atom);
@@ -1527,7 +1533,7 @@ public:
             return m_imp.bound_is_interesting(j, kind, v);
         }
 
-        virtual void consume(rational const& v, unsigned j) {
+        void consume(rational const& v, lp::constraint_index j) override {
             m_imp.set_evidence(j);
             m_imp.m_explanation.push_back(std::make_pair(v, j));
         }
@@ -2619,12 +2625,17 @@ public:
             coeff = rational::zero();
         }
         lp::impq term_max;
-        if (m_solver->maximize_term(coeffs, term_max)) {
+        lp::lp_status st = m_solver->maximize_term(coeffs, term_max);
+        if (st == lp::lp_status::OPTIMAL) {
             blocker = mk_gt(v);
             inf_rational val(term_max.x + coeff, term_max.y);
             return inf_eps(rational::zero(), val);
+        } if (st == lp::lp_status::FEASIBLE) {
+            lp_assert(false); // not implemented
+            // todo: what to return?
         }
         else {
+            SASSERT(st == lp::lp_status::UNBOUNDED);
             TRACE("arith", tout << "Unbounded " << v << "\n";);
             has_shared = false;
             blocker = m.mk_false();

src/tactic/smtlogics/qflia_tactic.cpp

@@ -228,18 +228,6 @@ tactic * mk_qflia_tactic(ast_manager & m, params_ref const & p) {
 #endif
                                main_p);
 
-    //
-    
-    
-    // tactic * st = using_params(and_then(preamble_st,
-    //                                     or_else(mk_ilp_model_finder_tactic(m),
-    //                                             mk_pb_tactic(m),
-    //                                             and_then(fail_if_not(mk_is_quasi_pb_probe()), 
-    //                                                      using_params(mk_lia2sat_tactic(m), quasi_pb_p),
-    //                                                      mk_fail_if_undecided_tactic()),
-    //                                             mk_bounded_tactic(m),
-    //                                             mk_smt_tactic())),
-    //                            main_p);
     
     st->updt_params(p);
     return st;

src/test/lp/lp.cpp

@@ -44,10 +44,8 @@
 #include "util/lp/numeric_pair.h"
 #include "util/lp/binary_heap_upair_queue.h"
 #include "util/lp/stacked_value.h"
-#include "util/lp/stacked_unordered_set.h"
 #include "util/lp/int_set.h"
 #include "util/stopwatch.h"
-#include "util/lp/stacked_map.h"
 #include <cstdlib>
 #include "test/lp/gomory_test.h"
 #include "util/lp/matrix.h"
@@ -55,9 +53,16 @@
 #include "util/lp/square_sparse_matrix_def.h"
 #include "util/lp/lu_def.h"
 #include "util/lp/general_matrix.h"
+#include "util/lp/bound_propagator.h"
 namespace lp {
 unsigned seed = 1;
 
+    class my_bound_propagator : public bound_propagator {
+    public:
+        my_bound_propagator(lar_solver & ls): bound_propagator(ls) {}
+        void consume(mpq const& v, lp::constraint_index j) override {}
+    };
+
 random_gen g_rand;
 static unsigned my_random() {
     return g_rand();
@@ -1942,44 +1947,6 @@ void setup_args_parser(argument_parser & parser) {
 
 struct fff { int a; int b;};
 
-void test_stacked_map_itself() {
-    vector<int> v(3,0);
-    for(auto u : v)
-        std::cout << u << std::endl;    
-    
-    std::unordered_map<int, fff> foo;
-    fff l;
-    l.a = 0;
-    l.b =1;
-    foo[1] = l;
-    int r = 1;
-    int k = foo[r].a;
-    std::cout << k << std::endl;
-    
-    stacked_map<int, double> m;
-    m[0] = 3;
-    m[1] = 4;
-    m.push();
-    m[1] = 5;
-    m[2] = 2;
-    m.pop();
-    m.erase(2);
-    m[2] = 3;
-    m.erase(1);
-    m.push();
-    m[3] = 100;
-    m[4] = 200;
-    m.erase(1);
-    m.push();
-    m[5] = 300;
-    m[6] = 400;
-    m[5] = 301;
-    m.erase(5);
-    m[3] = 122;
-
-    m.pop(2);
-    m.pop();
-}
 
 void test_stacked_unsigned() {
     std::cout << "test stacked unsigned" << std::endl;
@@ -2038,36 +2005,10 @@ void test_stacked_vector() {
         
 }
 
-void test_stacked_set() {
-#ifdef Z3DEBUG
-    std::cout << "test_stacked_set" << std::endl;
-    stacked_unordered_set<int> s;
-    s.insert(1);
-    s.insert(2);
-    s.insert(3);
-    std::unordered_set<int> scopy = s();
-    s.push();
-    s.insert(4);
-    s.pop();
-    lp_assert(s() == scopy);
-    s.push();
-    s.push();
-    s.insert(4);
-    s.insert(5);
-    s.push();
-    s.insert(4);
-    s.pop(3);
-    lp_assert(s() == scopy);
-#endif
-}
 
 void test_stacked() {
-    std::cout << "test_stacked_map()" << std::endl;
-    test_stacked_map_itself();
     test_stacked_value();
     test_stacked_vector();
-    test_stacked_set();
-    
 }
 
 char * find_home_dir() {
@@ -2815,7 +2756,7 @@ void test_bound_propagation_one_small_sample1() {
     vector<implied_bound> ev;
     ls.add_var_bound(a, LE, mpq(1));
     ls.solve();
-    bound_propagator bp(ls);
+    my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
     std::cout << " bound ev from test_bound_propagation_one_small_sample1" << std::endl;
     for (auto & be : bp.m_ibounds)  {
@@ -2868,7 +2809,7 @@ void test_bound_propagation_one_row() {
     vector<implied_bound> ev;
     ls.add_var_bound(x0, LE, mpq(1));
     ls.solve();
-    bound_propagator bp(ls);
+    my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
 } 
 void test_bound_propagation_one_row_with_bounded_vars() {
@@ -2884,7 +2825,7 @@ void test_bound_propagation_one_row_with_bounded_vars() {
     ls.add_var_bound(x0, LE, mpq(3));
     ls.add_var_bound(x0, LE, mpq(1));
     ls.solve();
-    bound_propagator bp(ls);
+    my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
 }
 void test_bound_propagation_one_row_mixed() {
@@ -2898,7 +2839,7 @@ void test_bound_propagation_one_row_mixed() {
     vector<implied_bound> ev;
     ls.add_var_bound(x1, LE, mpq(1));
     ls.solve();
-    bound_propagator bp(ls);
+    my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
 } 
 
@@ -2921,7 +2862,7 @@ void test_bound_propagation_two_rows() {
     vector<implied_bound> ev;
     ls.add_var_bound(y, LE, mpq(1));
     ls.solve();
-    bound_propagator bp(ls);
+    my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
 } 
 
@@ -2941,7 +2882,7 @@ void test_total_case_u() {
     vector<implied_bound> ev;
     ls.add_var_bound(z, GE, zero_of_type<mpq>());
     ls.solve();
-    bound_propagator bp(ls);
+    my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
 }
 bool contains_j_kind(unsigned j, lconstraint_kind kind, const mpq & rs, const vector<implied_bound> & ev) {
@@ -2968,7 +2909,7 @@ void test_total_case_l(){
     vector<implied_bound> ev;
     ls.add_var_bound(z, LE, zero_of_type<mpq>());
     ls.solve();
-    bound_propagator bp(ls);
+    my_bound_propagator bp(ls);
     ls.propagate_bounds_for_touched_rows(bp);
     lp_assert(ev.size() == 4);
     lp_assert(contains_j_kind(x, GE, - one_of_type<mpq>(), ev));

src/test/simplex.cpp

@@ -4,7 +4,6 @@ Copyright (c) 2015 Microsoft Corporation
 
 --*/
 
-#include "util/lp/sparse_matrix.h"
 #include "math/simplex/sparse_matrix_def.h"
 #include "math/simplex/simplex.h"
 #include "math/simplex/simplex_def.h"

src/util/lp/bound_propagator.h

@@ -22,6 +22,6 @@ public:
                                       lp::lconstraint_kind kind,
                                       const rational & bval) {return true;}
     unsigned number_of_found_bounds() const { return m_ibounds.size(); }
-    virtual void consume(mpq const& v, unsigned j) { std::cout << "doh\n"; }
+    virtual void consume(mpq const& v, lp::constraint_index j) = 0;
 };
 }

src/util/lp/explanation.h

@@ -20,6 +20,7 @@ Revision History:
 #pragma once
 namespace lp {
 struct explanation {
+    void clear() { m_explanation.clear(); }
     vector<std::pair<mpq, constraint_index>> m_explanation;
     void push_justification(constraint_index j, const mpq& v) {
         m_explanation.push_back(std::make_pair(v, j));

src/util/lp/hnf.h

@@ -32,7 +32,7 @@ Revision History:
 namespace lp {
 namespace hnf_calc {
 
-    // d = u * a + b * b and the sum of abs(u) + abs(v) is minimal, d is positive
+    // d = u * a + v * b and the sum of abs(u) + abs(v) is minimal, d is positive
 inline
 void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq & v) {
     if (is_zero(a)) {
@@ -47,7 +47,11 @@ void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq
         d = a;
         return;
     }
+#if 1
+    d = gcd(a, b, u, v);
+#else
     extended_gcd(a, b, d, u, v);
+#endif
     if (is_neg(d)) {
         d = -d;
         u = -u;
@@ -103,7 +107,8 @@ void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq
 
 
 
-template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
+template <typename M> 
+bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
     for (unsigned i = r; i < m.row_count(); i++) {
         for (unsigned j = r; j < m.column_count(); j++) {
             if (!is_zero(m[i][j])) {
@@ -120,13 +125,13 @@ template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
     return false;
 }
 
-template <typename M> void pivot_column_non_fractional(M &m, unsigned r) {
+template <typename M> 
+void pivot_column_non_fractional(M &m, unsigned r) {
     lp_assert(!is_zero(m[r][r]));
-    lp_assert(m.row_count() <= m.column_count());
     for (unsigned j = r + 1; j < m.column_count(); j++) {
         for (unsigned i = r + 1; i  < m.row_count(); i++) {
             m[i][j] =
-                (r > 0 )?
+                (r > 0) ?
                 (m[r][r]*m[i][j] - m[i][r]*m[r][j]) / m[r-1][r-1] :
                 (m[r][r]*m[i][j] - m[i][r]*m[r][j]);
             lp_assert(is_int(m[i][j]));
@@ -140,8 +145,8 @@ template <typename M> void pivot_column_non_fractional(M &m, unsigned r) {
 }
 
 // returns the rank of the matrix
-template <typename M> unsigned to_lower_triangle_non_fractional(M &m) {
+template <typename M> 
-    lp_assert(m.row_count() <= m.column_count());
+unsigned to_lower_triangle_non_fractional(M &m) {
     unsigned i = 0;
     for (; i < m.row_count(); i++) {
         if (!prepare_pivot_for_lower_triangle(m, i)) {
@@ -152,6 +157,8 @@ template <typename M> unsigned to_lower_triangle_non_fractional(M &m) {
     lp_assert(i == m.row_count());
     return i; 
 }
+
+// returns gcd of values below diagonal i,i
 template <typename M>
 mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) {
     mpq g = zero_of_type<mpq>();
@@ -170,28 +177,28 @@ mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) {
     return g;
 }
 
-
+// It fills "r" - the basic rows of m.
-// it fills "r" - the basic rows of m
+// The plan is to transform m to the lower triangular form by using non-fractional Gaussian Elimination by columns.
-template <typename M> mpq determinant_of_rectangular_matrix(const M& m, svector<unsigned> & basis_rows) {
+// Then the trailing after the diagonal elements of the following elements of the last non-zero row of the matrix,
-    if (m.column_count() < m.row_count())
+// namely, m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r.
-        throw "not implemented"; // consider working with the transposed m, or create a "transposed" version if needed
+// The gcd of these minors is the return value.
-    // The plan is to transform m to the lower triangular form by using non-fractional Gaussian Elimination by columns.
+ 
-    // Then the elements of the following elements of the last non-zero row of the matrix
+template <typename M> 
-    // m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r.
+mpq determinant_of_rectangular_matrix(const M& m, svector<unsigned> & basis_rows) {
-    // The gcd of these minors is the return value
+    auto m_copy = m;
-    auto mc = m;
+    unsigned rank = to_lower_triangle_non_fractional(m_copy);
-    unsigned rank = to_lower_triangle_non_fractional(mc);
     if (rank == 0)
         return one_of_type<mpq>();
 
     for (unsigned i = 0; i < rank; i++) {
-        basis_rows.push_back(mc.adjust_row(i));
+        basis_rows.push_back(m_copy.adjust_row(i));
     }
-    TRACE("hnf_calc", tout << "basis_rows = "; print_vector(basis_rows, tout); mc.print(tout, "mc = "););
+    TRACE("hnf_calc", tout << "basis_rows = "; print_vector(basis_rows, tout); m_copy.print(tout, "m_copy = "););
-    return gcd_of_row_starting_from_diagonal(mc, rank - 1);
+    return gcd_of_row_starting_from_diagonal(m_copy, rank - 1);
 }
 
-template <typename M> mpq determinant(const M& m) {
+template <typename M> 
+mpq determinant(const M& m) {
     lp_assert(m.row_count() == m.column_count());
     auto mc = m;
     svector<unsigned> basis_rows;
@@ -229,10 +236,8 @@ class hnf {
 
     mpq mod_R(const mpq & a) const {
        mpq t = a % m_R;
-       t = is_neg(t) ? t + m_R : t;
+       t = is_neg(t) ? t + m_R : t;       
-       if(is_neg(t)){
+       CTRACE("hnf", is_neg(t), tout << "a=" << a << ", m_R= " << m_R << std::endl;);
-           std::cout << "a=" << a << ", m_R= " << m_R << std::endl;
-       }
        return t;
            
     }
@@ -319,8 +324,8 @@ class hnf {
         }
     }
 
- void handle_column_ij_in_row_i(unsigned i, unsigned j) {
+    void handle_column_ij_in_row_i(unsigned i, unsigned j) {
-     lp_assert(is_correct_modulo());
+        lp_assert(is_correct_modulo());
         const mpq& aii = m_H[i][i];
         const mpq& aij = m_H[i][j];
         mpq p,q,r;
@@ -470,15 +475,14 @@ class hnf {
 
     bool is_correct_final() const {
         if (!is_correct()) {
-            std::cout << "m_H =            "; m_H.print(std::cout, 17);
+            TRACE("hnf_calc",
-            
+                  tout << "m_H =            "; m_H.print(tout, 17);
-            std::cout << "\nm_A_orig * m_U = ";  (m_A_orig * m_U).print(std::cout, 17);
+                  tout << "\nm_A_orig * m_U = ";  (m_A_orig * m_U).print(tout, 17);
-            
+                  tout << "is_correct() does not hold" << std::endl;);
-            std::cout << "is_correct() does not hold" << std::endl;
             return false;
         }
         if (!is_correct_form()) {
-            std::cout << "is_correct_form() does not hold" << std::endl;
+            TRACE("hnf_calc", tout << "is_correct_form() does not hold" << std::endl;);
             return false;
         }
         return true;
@@ -538,10 +542,6 @@ private:
         copy_buffer_to_col_i_W_modulo();
     }
 
-    void endl() const {
-        std::cout << std::endl;
-    }
-    
     void fix_row_under_diagonal_W_modulo() {
         mpq d, u, v;
         if (is_zero(m_W[m_i][m_i])) {
@@ -616,12 +616,11 @@ public:
 #endif
         calculate_by_modulo();
 #ifdef Z3DEBUG
-        if (m_H != m_W) {
+        CTRACE("hnf_calc", m_H != m_W,
-            std::cout << "A = "; m_A_orig.print(std::cout, 4); endl();
+               tout << "A = "; m_A_orig.print(tout, 4); tout << std::endl;
-            std::cout << "H = "; m_H.print(std::cout, 4); endl();
+               tout << "H = "; m_H.print(tout, 4);  tout << std::endl;
-            std::cout << "W = "; m_W.print(std::cout, 4); endl();
+               tout << "W = "; m_W.print(tout, 4);  tout << std::endl;);
-            lp_assert(false);
+        lp_assert (m_H == m_W);
-        }
 #endif
     }
 

src/util/lp/hnf_cutter.h

@@ -29,10 +29,12 @@ class hnf_cutter {
     var_register               m_var_register;
     general_matrix             m_A;
     vector<const lar_term*>    m_terms;
+    svector<constraint_index>  m_constraints_for_explanation;
     vector<mpq>                m_right_sides;
     unsigned                   m_row_count;
     unsigned                   m_column_count;
     lp_settings &              m_settings;
+    
 public:
     hnf_cutter(lp_settings & settings) : m_row_count(0), m_column_count(0), m_settings(settings) {}
 
@@ -41,25 +43,29 @@ public:
     }
 
     const vector<const lar_term*>& terms() const { return m_terms; }
+    const svector<unsigned>& constraints_for_explanation() const {
+        return m_constraints_for_explanation;
+    }
     const vector<mpq> & right_sides() const { return m_right_sides; }
     void clear() {
         // m_A will be filled from scratch in init_matrix_A
         m_var_register.clear();
         m_terms.clear();
-        m_right_sides.clear();
+        m_constraints_for_explanation.clear();
+        m_right_sides.clear();        
         m_row_count = m_column_count = 0;
     }
-    void add_term(const lar_term* t, const mpq &rs) {
+    void add_term(const lar_term* t, const mpq &rs, constraint_index ci) {
         m_terms.push_back(t);
+        m_right_sides.push_back(rs);
+        m_constraints_for_explanation.push_back(ci);
         for (const auto &p : *t) {
             m_var_register.add_var(p.var());
         }
-        m_right_sides.push_back(rs);
         if (m_terms.size() <= m_var_register.size()) { // capture the maximal acceptable sizes
             m_row_count = m_terms.size();
             m_column_count = m_var_register.size();
         }
-        
     }
     void print(std::ostream & out) {
         out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl;
@@ -93,8 +99,9 @@ public:
         if (basis_rows.size() == m_right_sides.size())
             return m_right_sides;
         vector<mpq> b;
-        for (unsigned i : basis_rows)
+        for (unsigned i : basis_rows) {
             b.push_back(m_right_sides[i]);
+        }
         return b;
     }
 
@@ -161,19 +168,30 @@ public:
         return ret;
     }
 #endif
+    void shrink_explanation(const svector<unsigned>& basis_rows) {
+        svector<unsigned> new_expl;
+        for (unsigned i : basis_rows) {
+            new_expl.push_back(m_constraints_for_explanation[i]);
+        }
+        m_constraints_for_explanation = new_expl;
+    }
+
     lia_move create_cut(lar_term& t, mpq& k, explanation& ex, bool & upper
                         #ifdef Z3DEBUG
                         ,
                         const vector<mpq> & x0
                         #endif
                         ) {
+        // we suppose that x0 has at least one non integer element 
         init_matrix_A();
         svector<unsigned> basis_rows;
         mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows);
         if (m_settings.get_cancel_flag())
             return lia_move::undef;
-        if (basis_rows.size() < m_A.row_count())
+        if (basis_rows.size() < m_A.row_count()) {
             m_A.shrink_to_rank(basis_rows);
+            shrink_explanation(basis_rows);
+        }
         
         hnf<general_matrix> h(m_A, d);
         //  general_matrix A_orig = m_A;
@@ -184,22 +202,10 @@ public:
         find_h_minus_1_b(h.W(), b);
         // lp_assert(bcopy == h.W().take_first_n_columns(b.size()) * b);
         int cut_row = find_cut_row_index(b);
-        if (cut_row == -1) {
+        if (cut_row == -1)
             return lia_move::undef;
-        }
+        // the matrix is not square - we can get
-        // test region
+        // all integers in b's projection
-        /*
-        general_matrix U(m_A.column_count());
-        vector<mpq> rt(m_A.column_count());
-        for (unsigned i = 0; i < U.row_count(); i++) {
-            get_ei_H_minus_1(i, h.W(), rt);
-            vector<mpq> ft = rt * A_orig;
-            for (unsigned j = 0; j < ft.size(); j++)
-                U[i][j] = ft[j];
-        }
-        std::cout << "U reverse = "; U.print(std::cout, 12); std::cout << std::endl;
-        */
-        // end test region
         
         vector<mpq> row(m_A.column_count());
         get_ei_H_minus_1(cut_row, h.W(), row);
@@ -209,5 +215,7 @@ public:
         upper = true;
         return lia_move::cut;
     }
+
+    svector<unsigned> vars() const { return m_var_register.vars(); }
 };
 }

src/util/lp/int_solver.cpp

@@ -409,6 +409,8 @@ struct pivoted_rows_tracking_control {
     }
 };
 
+
+
 impq int_solver::get_cube_delta_for_term(const lar_term& t) const {
     if (t.size() == 2) {
         bool seen_minus = false;
@@ -528,34 +530,37 @@ lia_move int_solver::gomory_cut() {
 }
 
 
-bool int_solver::try_add_term_to_A_for_hnf(unsigned i) {
+void int_solver::try_add_term_to_A_for_hnf(unsigned i) {
     mpq rs;
     const lar_term* t = m_lar_solver->terms()[i];
+#if Z3DEBUG
     for (const auto & p : *t) {
         if (!is_int(p.var())) {
             lp_assert(false); 
-            return false; // todo : the mix case!
         }
     }
-    bool has_bounds;
+#endif
-    if (m_lar_solver->get_equality_and_right_side_for_term_on_corrent_x(i, rs, has_bounds)) {
+    constraint_index ci;
-        m_hnf_cutter.add_term(t, rs);
+    if (m_lar_solver->get_equality_and_right_side_for_term_on_current_x(i, rs, ci)) {
-        return true;
+        m_hnf_cutter.add_term(t, rs, ci);
     }
-    return !has_bounds;
 }
 
-bool int_solver::hnf_matrix_is_empty() const { return true; }
+bool int_solver::hnf_has_non_integral_var() const {
+    for (unsigned j : m_hnf_cutter.vars())
+        if (get_value(j).is_int() == false)
+            return true;
+    return false;
+}
 
 bool int_solver::init_terms_for_hnf_cut() {
     m_hnf_cutter.clear();
-    for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) {
+    for (unsigned i = 0; i < m_lar_solver->terms().size() && m_hnf_cutter.row_count() < settings().limit_on_rows_for_hnf_cutter; i++) {
         try_add_term_to_A_for_hnf(i);
     }
-    return m_hnf_cutter.row_count() < settings().limit_on_rows_for_hnf_cutter;
+    return hnf_has_non_integral_var();
 }
 
-
 lia_move int_solver::make_hnf_cut() {
     if (!init_terms_for_hnf_cut()) {
         return lia_move::undef;
@@ -574,6 +579,10 @@ lia_move int_solver::make_hnf_cut() {
     if (r == lia_move::cut) {        
         lp_assert(current_solution_is_inf_on_cut());
         settings().st().m_hnf_cuts++;
+        m_ex->clear();        
+        for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
+             m_ex->push_justification(i);
+        }
     }
     return r;
 }
@@ -599,7 +608,6 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
 
     r = patch_nbasic_columns();
     if (r != lia_move::undef) return r;
-
     ++m_branch_cut_counter;
 
     r = find_cube();

src/util/lp/int_solver.h

@@ -51,7 +51,9 @@ public:
     // or provide a way of how it can be adjusted.
     lia_move check(lar_term& t, mpq& k, explanation& ex, bool & upper);
     bool move_non_basic_column_to_bounds(unsigned j);
-    lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex);
+    lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex);    
+    bool is_base(unsigned j) const;
+    
 private:
 
     // how to tighten bounds for integer variables.
@@ -75,14 +77,12 @@ private:
                       mpq const & consts);
     void fill_explanation_from_fixed_columns(const row_strip<mpq> & row);
     void add_to_explanation_from_fixed_or_boxed_column(unsigned j);
-    void patch_nbasic_column(unsigned j);
     lia_move patch_nbasic_columns();
     bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m);
     const impq & lower_bound(unsigned j) const;
     const impq & upper_bound(unsigned j) const;
     bool is_int(unsigned j) const;
     bool is_real(unsigned j) const;
-    bool is_base(unsigned j) const;
     bool is_boxed(unsigned j) const;
     bool is_fixed(unsigned j) const;
     bool is_free(unsigned j) const;
@@ -156,6 +156,8 @@ public:
     lia_move make_hnf_cut();
     bool init_terms_for_hnf_cut();
     bool hnf_matrix_is_empty() const;
-    bool try_add_term_to_A_for_hnf(unsigned term_index);
+    void try_add_term_to_A_for_hnf(unsigned term_index);
+    bool hnf_has_non_integral_var() const;
+    void patch_nbasic_column(unsigned j);
 };
 }

src/util/lp/lar_core_solver.h

@@ -26,7 +26,6 @@ Revision History:
 #include "util/lp/indexed_vector.h"
 #include "util/lp/binary_heap_priority_queue.h"
 #include "util/lp/breakpoint.h"
-#include "util/lp/stacked_unordered_set.h"
 #include "util/lp/lp_primal_core_solver.h"
 #include "util/lp/stacked_vector.h"
 #include "util/lp/lar_solution_signature.h"
@@ -611,7 +610,6 @@ public:
         }
 
         if (no_r_lu()) { // it is the case where m_d_solver gives a degenerated basis, we need to roll back
-            //            std::cout << "no_r_lu" << std::endl;
             catch_up_in_lu_in_reverse(changes_of_basis, m_r_solver);
             m_r_solver.find_feasible_solution();
             m_d_basis = m_r_basis;

src/util/lp/lar_solver.cpp

@@ -370,7 +370,6 @@ void lar_solver::pop(unsigned k) {
     unsigned m = A_r().row_count();
     clean_popped_elements(m, m_rows_with_changed_bounds);
     clean_inf_set_of_r_solver_after_pop();
-    m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN;
     lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
               (!use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
         
@@ -508,14 +507,45 @@ bool lar_solver::maximize_term_on_corrected_r_solver(const vector<std::pair<mpq,
         lp_unreachable(); // wrong mode
     }
     return false;
-}    
+}
+
+
+bool lar_solver::remove_from_basis(unsigned j) {
+    return m_mpq_lar_core_solver.m_r_solver.remove_from_basis(j);
+}
+
+
 // starting from a given feasible state look for the maximum of the term
 // return true if found and false if unbounded
-bool lar_solver::maximize_term(const vector<std::pair<mpq, var_index>> & term,
+lp_status lar_solver::maximize_term(const vector<std::pair<mpq, var_index>> & term,
                                impq &term_max) {
     lp_assert(m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible());
     m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false;
-    return maximize_term_on_corrected_r_solver(term, term_max);
+    if (!maximize_term_on_corrected_r_solver(term, term_max))
+        return lp_status::UNBOUNDED;
+
+    bool change = false;
+    for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_x.size(); j++) {
+        if (!column_is_int(j))
+            continue;
+        if (column_value_is_integer(j))
+            continue;
+        if (m_int_solver->is_base(j)) {
+            if (!remove_from_basis(j)) // consider a special version of remove_from_basis that would not remove inf_int columns
+                return lp_status::FEASIBLE; // it should not happen
+        }
+        m_int_solver->patch_nbasic_column(j);
+        if (!column_value_is_integer(j))
+            return lp_status::FEASIBLE;
+        change = true;
+    }
+    if (change) {
+        term_max = zero_of_type<impq>();
+        for (const auto& p : term) {
+            term_max += p.first * m_mpq_lar_core_solver.m_r_x[p.second];
+        }
+    }
+    return lp_status::OPTIMAL;
 }
     
 
@@ -792,18 +822,11 @@ void lar_solver::add_last_rows_to_lu(lp_primal_core_solver<K,L> & s) {
     
 bool lar_solver::x_is_correct() const {
     if (m_mpq_lar_core_solver.m_r_x.size() != A_r().column_count()) {
-        //            std::cout << "the size is off " << m_r_solver.m_x.size() << ", " << A().column_count() << std::endl;
         return false;
     }
     for (unsigned i = 0; i < A_r().row_count(); i++) {
         numeric_pair<mpq> delta =  A_r().dot_product_with_row(i, m_mpq_lar_core_solver.m_r_x);
         if (!delta.is_zero()) {
-            // std::cout << "x is off (";
-            // std::cout << "m_b[" << i  << "] = " << m_b[i] << " ";
-            // std::cout << "left side = " << A().dot_product_with_row(i, m_r_solver.m_x) << ' ';
-            // std::cout << "delta = " << delta << ' ';
-            // std::cout << "iters = " << total_iterations() << ")" << std::endl;
-            // std::cout << "row " << i << " is off" << std::endl;
             return false;
         }
     }
@@ -906,7 +929,6 @@ bool lar_solver::all_constraints_hold() const {
     
     for (unsigned i = 0; i < m_constraints.size(); i++) {
         if (!constraint_holds(*m_constraints[i], var_map)) {
-            print_constraint(i, std::cout);
             return false;
         }
     }
@@ -973,13 +995,6 @@ bool lar_solver::the_left_sides_sum_to_zero(const vector<std::pair<mpq, unsigned
     }
 
     if (!coeff_map.empty()) {
-        std::cout << "left side = ";
-        vector<std::pair<mpq, var_index>> t;
-        for (auto & it : coeff_map) {
-            t.push_back(std::make_pair(it.second, it.first));
-        }
-        print_linear_combination_of_column_indices(t, std::cout);
-        std::cout << std::endl;
         return false;
     }
 
@@ -1043,7 +1058,7 @@ mpq lar_solver::sum_of_right_sides_of_explanation(const vector<std::pair<mpq, un
     return ret;
 }
 
-bool lar_solver::has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) {
+bool lar_solver::has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const {
 
     if (var >= m_columns_to_ul_pairs.size()) {
         // TBD: bounds on terms could also be used, caller may have to track these.
@@ -1062,7 +1077,7 @@ bool lar_solver::has_lower_bound(var_index var, constraint_index& ci, mpq& value
     }
 }
     
-bool lar_solver::has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) {
+bool lar_solver::has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const {
 
     if (var >= m_columns_to_ul_pairs.size()) {
         // TBD: bounds on terms could also be used, caller may have to track these.
@@ -1143,9 +1158,11 @@ void lar_solver::get_model(std::unordered_map<var_index, mpq> & variable_values)
 }
 
 void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map<var_index, mpq> & variable_values) const {
+    mpq delta = mpq(1);
+    delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(delta);
     for (unsigned i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++ ) {
         const impq & rp = m_mpq_lar_core_solver.m_r_x[i];
-        variable_values[i] = rp.x + rp.y;
+        variable_values[i] = rp.x + delta * rp.y;
     }
 }
 
@@ -1436,7 +1453,7 @@ bool lar_solver::model_is_int_feasible() const {
 
 bool lar_solver::term_is_int(const lar_term * t) const {
     for (auto const & p :  t->m_coeffs)
-        if (!column_is_int(p.first) || p.second.is_int())
+        if (! (column_is_int(p.first)  && p.second.is_int()))
             return false;
     return t->m_v.is_int();
 }
@@ -1481,9 +1498,6 @@ var_index lar_solver::add_var(unsigned ext_j, bool is_int) {
     TRACE("add_var", tout << "adding var " << ext_j << (is_int? " int" : " nonint") << std::endl;);
     var_index i;
     lp_assert(ext_j < m_terms_start_index);
-
-    if (ext_j >= m_terms_start_index)
-        throw 0; // todo : what is the right way to exit?
     auto it = m_ext_vars_to_columns.find(ext_j);
     if (it != m_ext_vars_to_columns.end()) {
         return it->second.internal_j();
@@ -1580,7 +1594,6 @@ bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & co
     bool g_is_set = false;
     for (const auto & p : coeffs) {
         if (!p.first.is_int()) {
-            std::cout << "p.first = " << p.first << " is not an int\n";
             return false;
         }
         if (!g_is_set) {
@@ -1593,9 +1606,6 @@ bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & co
     if (g == one_of_type<mpq>())
         return true;
 
-    std::cout << "term is not ok: g = " << g << std::endl;
-    this->print_linear_combination_of_column_indices_only(coeffs, std::cout);
-    std::cout << " rs = " << v << std::endl;
     return false;
 }
 #endif
@@ -2207,36 +2217,34 @@ void lar_solver::round_to_integer_solution() {
     }
 }
 
-bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & rs, bool & has_bounds) const {
+bool lar_solver::get_equality_and_right_side_for_term_on_current_x(unsigned term_index, mpq & rs, constraint_index& ci) const {
     unsigned tj = term_index + m_terms_start_index;
     auto it = m_ext_vars_to_columns.find(tj);
-    has_bounds = false;
     if (it == m_ext_vars_to_columns.end())
         return false;
     unsigned j = it->second.internal_j();
-    auto & slv = m_mpq_lar_core_solver.m_r_solver;
     impq term_val;
     bool term_val_ready = false;
-    if (slv.column_has_upper_bound(j)) {
+    mpq b;
-        has_bounds = true;
+    bool is_strict;
-        const impq & b = slv.m_upper_bounds[j];
+    if (has_upper_bound(j, ci, b, is_strict)) {
-        lp_assert(is_zero(b.y) && is_int(b.x));
+        lp_assert(!is_strict);
+        lp_assert(b.is_int());
         term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x);
         term_val_ready = true;
-        if (term_val == b) {
+        if (term_val.x == b) {
-            rs = b.x;
+            rs = b;
             return true;
         }
     }
-    if (slv.column_has_lower_bound(j)) {
+    if (has_lower_bound(j, ci, b, is_strict)) {
-        has_bounds = true;
+        lp_assert(!is_strict);
         if (!term_val_ready)
             term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x);
-        const impq & b = slv.m_lower_bounds[j];
+        lp_assert(b.is_int());
-        lp_assert(is_zero(b.y) && is_int(b.x));
         
-        if (term_val == b) {
+        if (term_val.x == b) {
-            rs = b.x;
+            rs = b;
             return true;
         }
     }

src/util/lp/lar_solver.h

@@ -36,7 +36,6 @@ Revision History:
 #include <stack>
 #include "util/lp/stacked_value.h"
 #include "util/lp/stacked_vector.h"
-#include "util/lp/stacked_unordered_set.h"
 #include "util/lp/implied_bound.h"
 #include "util/lp/bound_analyzer_on_row.h"
 #include "util/lp/conversion_helper.h"
@@ -326,8 +325,8 @@ public:
                                              impq &term_max);    
     // starting from a given feasible state look for the maximum of the term
     // return true if found and false if unbounded
-    bool maximize_term(const vector<std::pair<mpq, var_index>> & term,
+    lp_status maximize_term(const vector<std::pair<mpq, var_index>> & term,
-                       impq &term_max);
+                            impq &term_max);
     
 
     
@@ -435,9 +434,9 @@ public:
 
     mpq sum_of_right_sides_of_explanation(const vector<std::pair<mpq, unsigned>> & explanation) const;
 
-    bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict);
+    bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const;
     
-    bool has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict);
+    bool has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const;
 
 
     void get_infeasibility_explanation(vector<std::pair<mpq, constraint_index>> & explanation) const;
@@ -587,6 +586,7 @@ public:
     unsigned column_count() const { return A_r().column_count(); }
     const vector<unsigned> & r_basis() const { return m_mpq_lar_core_solver.r_basis(); }
     const vector<unsigned> & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); }
-    bool get_equality_and_right_side_for_term_on_corrent_x(unsigned i, mpq &rs, bool & has_bounds) const;
+    bool get_equality_and_right_side_for_term_on_current_x(unsigned i, mpq &rs, constraint_index& ci) const;
+    bool remove_from_basis(unsigned);
 };
 }

src/util/lp/lp_core_solver_base.cpp

@@ -145,3 +145,4 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::infe
 template bool lp::lp_core_solver_base<lp::mpq, lp::mpq >::infeasibility_costs_are_correct() const;
 template bool lp::lp_core_solver_base<double, double >::infeasibility_costs_are_correct() const;
 template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::calculate_pivot_row(unsigned int);
+template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::remove_from_basis(unsigned int);

src/util/lp/lp_core_solver_base.h

@@ -225,7 +225,6 @@ public:
         lp_assert(m_A.is_correct());
         if (m_using_infeas_costs) {
             if (infeasibility_costs_are_correct() == false) {
-                std::cout << "infeasibility_costs_are_correct() does not hold" << std::endl;
                 return false;
             }
         }
@@ -234,9 +233,6 @@ public:
         for (unsigned j = 0; j < n; j++) {
             if (m_basis_heading[j] >= 0) {
                 if (!is_zero(m_d[j])) {
-                    
-                    std::cout << "case a\n";
-                    print_column_info(j, std::cout);
                     return false;
                 }
             } else {
@@ -245,8 +241,6 @@ public:
                     d -= this->m_costs[this->m_basis[cc.m_i]] * this->m_A.get_val(cc);
                 }
                 if (m_d[j] != d) {
-                    std::cout << "case b\n";
-                    print_column_info(j, std::cout);
                     return false;
                 }
             }
@@ -453,6 +447,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 remove_from_basis(unsigned j);
     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);
@@ -552,7 +547,6 @@ public:
     bool non_basic_columns_are_set_correctly() const {
         for (unsigned j : this->m_nbasis)
             if (!column_is_feasible(j)) {
-                print_column_info(j, std::cout);
                 return false;
             }
         return true;
@@ -601,10 +595,6 @@ public:
             out << " base\n";
         else
             out << " nbas\n";
-                
-        /*
-        std::cout << "cost = " << m_costs[j] << std::endl;
-        std:: cout << "m_d = " << m_d[j] << std::endl;*/
     }
 
     bool column_is_free(unsigned j) const { return this->m_column_type[j] == free; }

src/util/lp/lp_core_solver_base_def.h

@@ -228,11 +228,6 @@ A_mult_x_is_off() const {
 		for (unsigned i = 0; i < m_m(); i++) {
             X delta = m_b[i] - m_A.dot_product_with_row(i, m_x);
             if (delta != numeric_traits<X>::zero()) {
-                std::cout << "precise x is off (";
-                std::cout << "m_b[" << i << "] = " << m_b[i] << " ";
-                std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
-                std::cout << "delta = " << delta << ' ';
-                std::cout << "iters = " << total_iterations() << ")" << std::endl;
                 return true;
             }
         }
@@ -265,11 +260,6 @@ A_mult_x_is_off_on_index(const vector<unsigned> & index) const {
     for (unsigned i : index) {
         X delta = m_b[i] - m_A.dot_product_with_row(i, m_x);
         if (delta != numeric_traits<X>::zero()) {
-            // std::cout << "x is off (";
-            // std::cout << "m_b[" << i  << "] = " << m_b[i] << " ";
-            // std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
-            // std::cout << "delta = " << delta << ' ';
-            // std::cout << "iters = " << total_iterations() << ")" << std::endl;
             return true;
         }
     }
@@ -416,13 +406,10 @@ column_is_dual_feasible(unsigned j) const {
     case column_type::lower_bound:
         return x_is_at_lower_bound(j) && d_is_not_negative(j);
     case column_type::upper_bound:
-        LP_OUT(m_settings,  "upper_bound type should be switched to lower_bound" << std::endl);
         lp_assert(false); // impossible case
     case column_type::free_column:
         return numeric_traits<X>::is_zero(m_d[j]);
     default:
-        LP_OUT(m_settings,  "column = " << j << std::endl);
-        LP_OUT(m_settings,  "unexpected column type = " << column_type_to_string(m_column_types[j]) << std::endl);
         lp_unreachable();
     }
     lp_unreachable();
@@ -530,9 +517,6 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::inf_set_is_cor
         bool is_feas = column_is_feasible(j);
     
         if (is_feas == belongs_to_set) {
-            print_column_info(j, std::cout);
-            std::cout << "belongs_to_set = " << belongs_to_set << std::endl;
-            std::cout <<( is_feas? "feas":"inf") << std::endl;
             return false;
         }
     }
@@ -714,22 +698,18 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::
     lp_assert(m_basis.size() == m_A.row_count());
     lp_assert(m_nbasis.size() <= m_A.column_count() - m_A.row_count()); // for the dual the size of non basis can be smaller
     if (!basis_has_no_doubles()) {
-        //        std::cout << "basis_has_no_doubles" << std::endl;
         return false;
     }
 
     if (!non_basis_has_no_doubles()) {
-        // std::cout << "non_basis_has_no_doubles" << std::endl;
         return false;
     }
 
     if (!basis_is_correctly_represented_in_heading()) {
-        // std::cout << "basis_is_correctly_represented_in_heading" << std::endl;
         return false;
     }
 
     if (!non_basis_is_correctly_represented_in_heading()) {
-        // std::cout << "non_basis_is_correctly_represented_in_heading" << std::endl;
         return false;
     }
 
@@ -989,6 +969,17 @@ template <typename T, typename X>  void lp_core_solver_base<T, X>::pivot_fixed_v
     }
 }
 
+template <typename T, typename X> bool lp_core_solver_base<T, X>::remove_from_basis(unsigned basic_j) {
+    indexed_vector<T> w(m_basis.size()); // the buffer
+    unsigned i = m_basis_heading[basic_j];
+    for (auto &c : m_A.m_rows[i]) {
+        if (pivot_column_general(c.var(), basic_j, w))
+            return true;
+    }
+    return false;
+}
+
+
 template <typename T, typename X> bool 
 lp_core_solver_base<T, X>::infeasibility_costs_are_correct() const {
     if (! this->m_using_infeas_costs)
@@ -996,13 +987,9 @@ lp_core_solver_base<T, X>::infeasibility_costs_are_correct() const {
     lp_assert(costs_on_nbasis_are_zeros());
     for (unsigned j :this->m_basis) {
         if (!infeasibility_cost_is_correct_for_column(j)) {
-            std::cout << "infeasibility_cost_is_correct_for_column does not hold\n";
-            print_column_info(j, std::cout);
             return false;
         }
         if (!is_zero(m_d[j])) {
-            std::cout << "m_d is not zero\n";
-            print_column_info(j, std::cout);
             return false;
         }
     }
@@ -1052,9 +1039,9 @@ void lp_core_solver_base<T, X>::calculate_pivot_row(unsigned i) {
     m_pivot_row.clear();
     m_pivot_row.resize(m_n());
     if (m_settings.use_tableau()) {
-        unsigned basis_j = m_basis[i];
+        unsigned basic_j = m_basis[i];
         for (auto & c : m_A.m_rows[i]) {
-            if (c.m_j != basis_j)
+            if (c.m_j != basic_j)
                 m_pivot_row.set_value(c.get_val(), c.m_j);
         }
         return;

src/util/lp/lp_dual_core_solver_def.h

@@ -535,12 +535,6 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::snap_runaway_n
 template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dual_feasible() const {
     for (unsigned j : this->non_basis()){
         if (!this->column_is_dual_feasible(j)) {
-            // std::cout << "column " << j << " is not dual feasible" << std::endl;
-            // std::cout << "m_d[" << j << "] = " << this->m_d[j] << std::endl;
-            // std::cout << "x[" << j << "] = " << this->m_x[j] << std::endl;
-            // std::cout << "type = " << column_type_to_string(this->m_column_type[j]) << std::endl;
-            // std::cout << "bounds = " << this->m_lower_bounds[j] << "," << this->m_upper_bounds[j] << std::endl;
-            // std::cout << "total_iterations = " << this->total_iterations() << std::endl;
             return false;
         }
     }

src/util/lp/lp_primal_core_solver.h

@@ -504,7 +504,8 @@ public:
         }
         X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
         advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta );
-        lp_assert(this->m_x[leaving] == new_val_for_leaving);
+        X xleaving = this->m_x[leaving];
+        lp_assert(xleaving == new_val_for_leaving);
         if (this->current_x_is_feasible())
             this->set_status(lp_status::OPTIMAL);
     }

src/util/lp/lp_primal_core_solver_tableau_def.h

@@ -322,7 +322,6 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
     return m_leaving_candidates[k];
 }
 template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tableau() {
-        //        print_matrix(&(this->m_A), std::cout);
     CASSERT("A_off", this->A_mult_x_is_off() == false);
         lp_assert(basis_columns_are_set_correctly());
         this->m_basis_sort_counter = 0; // to initiate the sort of the basis

src/util/lp/lp_primal_simplex_def.h

@@ -278,7 +278,6 @@ template <typename T, typename X> bool lp_primal_simplex<T, X>::bounds_hold(std:
         }
 
         if (!it.second->bounds_hold(sol_it->second)) {
-            //            std::cout << "bounds do not hold for " << it.second->get_name() << std::endl;
             it.second->bounds_hold(sol_it->second);
             return false;
         }

src/util/lp/lp_settings.h

@@ -214,7 +214,7 @@ public:
     
     lp_settings() : m_default_resource_limit(*this),
                     m_resource_limit(&m_default_resource_limit),
-                    m_debug_out( &std::cout),
+                    m_debug_out(&std::cout),
                     m_message_out(&std::cout),
                     reps_in_scaler(20),
                     pivot_epsilon(0.00000001),
@@ -231,7 +231,6 @@ public:
                     drop_tolerance ( 1e-14),
                     tolerance_for_artificials ( 1e-4),
                     can_be_taken_to_basis_tolerance ( 0.00001),
-                  
                     percent_of_entering_to_check ( 5),// we try to find a profitable column in a percentage of the columns
                     use_scaling ( true),
                     scaling_maximum ( 1),
@@ -265,7 +264,7 @@ public:
                     m_chase_cut_solver_cycle_on_var(10),
                     m_int_pivot_fixed_vars_from_basis(false),
                     m_int_patch_only_integer_values(true),
-                    limit_on_rows_for_hnf_cutter(100)
+                    limit_on_rows_for_hnf_cutter(75)
     {}
 
     void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; }

src/util/lp/lp_settings_def.h

@@ -74,14 +74,12 @@ bool vectors_are_equal(T * a, vector<T>  &b, unsigned n) {
     if (numeric_traits<T>::precise()) {
         for (unsigned i = 0; i < n; i ++){
             if (!numeric_traits<T>::is_zero(a[i] - b[i])) {
-                // std::cout << "a[" << i <<"]" << a[i] << ", " << "b[" << i <<"]" << b[i] << std::endl;
                 return false;
             }
         }
     } else {
         for (unsigned i = 0; i < n; i ++){
             if (std::abs(numeric_traits<T>::get_double(a[i] - b[i])) > 0.000001) {
-                // std::cout << "a[" << i <<"]" << a[i] << ", " << "b[" << i <<"]" << b[i] << std::endl;
                 return false;
             }
         }
@@ -97,7 +95,6 @@ bool vectors_are_equal(const vector<T> & a, const vector<T>  &b) {
     if (numeric_traits<T>::precise()) {
         for (unsigned i = 0; i < n; i ++){
             if (!numeric_traits<T>::is_zero(a[i] - b[i])) {
-                // std::cout << "a[" << i <<"]" << a[i] << ", " << "b[" << i <<"]" << b[i] << std::endl;
                 return false;
             }
         }
@@ -112,7 +109,6 @@ bool vectors_are_equal(const vector<T> & a, const vector<T>  &b) {
             }
                 
             if (fabs(da - db) > 0.000001) {
-                // std::cout << "a[" << i <<"] = " << a[i] << ", but " << "b[" << i <<"] = " << b[i] << std::endl;
                 return false;
             }
         }

src/util/lp/lu_def.h

@@ -144,8 +144,6 @@ lu<M>::lu(const M& A,
     m_row_eta_work_vector(A.row_count()),
     m_refactor_counter(0) {
     lp_assert(A.row_count() == A.column_count());
-    std::cout << "m_U = \n";
-    print_matrix(&m_U, std::cout);
     create_initial_factorization();
 #ifdef Z3DEBUG
     lp_assert(is_correct());
@@ -659,13 +657,10 @@ template <typename M>
 bool lu<M>::is_correct() {
 #ifdef Z3DEBUG
     if (get_status() != LU_status::OK) {
-        std::cout << " status" << std::endl;
         return false;
     }
     dense_matrix<T, X> left_side = get_left_side();
-    std::cout << "ls = "; print_matrix(left_side, std::cout);
     dense_matrix<T, X> right_side = get_right_side();
-    std::cout << "rs = "; print_matrix(right_side, std::cout);
     return left_side == right_side;
 #else
     return true;

src/util/lp/mps_reader.h

@@ -303,7 +303,6 @@ class mps_reader {
         do {
             read_line();
             if (m_line.find("RHS") == 0) {
-                //  cout << "found RHS" << std::endl;
                 break;
             }
             if (m_line.size() < 22) {

src/util/lp/permutation_matrix_def.h

@@ -64,7 +64,6 @@ void permutation_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings & ) {
     // L * deb_w = clone_vector<L>(w, row_count());
     // deb.apply_from_left(deb_w);
 #endif
-    // std::cout << " apply_from_left " << std::endl;
     lp_assert(m_X_buffer.size() == w.size());
     unsigned i = size();
     while (i-- > 0) {

src/util/lp/scaler_def.h

@@ -214,7 +214,6 @@ template <typename T, typename X>    void scaler<T, X>::scale_rows() {
 }
 
 template <typename T, typename X>    void scaler<T, X>::scale_row(unsigned i) {
-    std::cout << "t" << "\n";
     T row_max = std::max(m_A.get_max_abs_in_row(i), abs(convert_struct<T, X>::convert(m_b[i])));
     T alpha = numeric_traits<T>::one();
     if (numeric_traits<T>::is_zero(row_max)) {
@@ -244,7 +243,6 @@ template <typename T, typename X>    void scaler<T, X>::scale_column(unsigned i)
     if (numeric_traits<T>::is_zero(column_max)){
         return; // the column has zeros only
     }
-    std::cout << "f";
     if (column_max < m_scaling_minimum) {
         do {
             alpha *= 2;

src/util/lp/square_sparse_matrix.h

@@ -254,11 +254,7 @@ public:
     void delete_column(int i);
 
     void swap_columns(unsigned a, unsigned b) {
-        // cout << "swaapoiiin" << std::endl;
-        // dense_matrix<T, X> d(*this);
         m_column_permutation.transpose_from_left(a, b);
-        // d.swap_columns(a, b);
-        // lp_assert(*this == d);
     }
 
     void swap_rows(unsigned a, unsigned b) {

src/util/lp/square_sparse_matrix_def.h

@@ -1214,7 +1214,6 @@ bool square_sparse_matrix<T, X>::is_upper_triangular_and_maximums_are_set_correc
                 return false;
             if (aj == ai) {
                 if (iv.m_value != 1) {
-                    // std::cout << "value at diagonal = " << iv.m_value << std::endl;
                     return false;
                 }
             }
@@ -1245,18 +1244,14 @@ void square_sparse_matrix<T, X>::check_column_vs_rows(unsigned col) {
     for (indexed_value<T> & column_iv : mc) {
         indexed_value<T> & row_iv = column_iv_other(column_iv);
         if (row_iv.m_index != col) {
-            //            std::cout << "m_other in row does not belong to column " << col << ", but to column  " << row_iv.m_index << std::endl;
             lp_assert(false);
         }
 
         if (& row_iv_other(row_iv) != &column_iv) {
-            // std::cout << "row and col do not point to each other" << std::endl;
             lp_assert(false);
         }
 
         if (row_iv.m_value != column_iv.m_value) {
-            // std::cout << "the data from col " << col << " for row " << column_iv.m_index << " is different in the column " << std::endl;
-            // std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl;
             lp_assert(false);
         }
     }
@@ -1269,18 +1264,14 @@ void square_sparse_matrix<T, X>::check_row_vs_columns(unsigned row) {
         indexed_value<T> & column_iv = row_iv_other(row_iv);
 
         if (column_iv.m_index != row) {
-            // std::cout << "col_iv does not point to correct row " << row << " but to " << column_iv.m_index << std::endl;
             lp_assert(false);
         }
 
         if (& row_iv != & column_iv_other(column_iv)) {
-            // std::cout << "row and col do not point to each other" << std::endl;
             lp_assert(false);
         }
 
         if (row_iv.m_value != column_iv.m_value) {
-            // std::cout << "the data from col " << column_iv.m_index << " for row " << row << " is different in the column " << std::endl;
-            // std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl;
             lp_assert(false);
         }
     }

src/util/lp/stacked_map.h

@@ -1,255 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-
-#pragma once
-// this class implements a map with some stack functionality
-#include <unordered_map>
-#include <unordered_set>
-#include <set>
-#include <stack>
-#include <map>
-namespace lp {
-template<class A,
-		class B,
-		class _Pr = std::less<A>,
-		class _Alloc = std::allocator<std::pair<const A, B> > >
-class stacked_map {
-    struct delta {
-        std::unordered_set<A> m_new;
-        std::unordered_map<A, B> m_original_changed;
-        //        std::unordered_map<A,B, Hash, KeyEqual, Allocator > m_deb_copy;
-    };
-    std::map<A,B,_Pr, _Alloc> m_map;
-    std::stack<delta> m_stack;
-public:
-    class ref {
-        stacked_map<A,B,_Pr, _Alloc> & m_map;
-        const A & m_key;
-    public:
-        ref(stacked_map<A,B,_Pr, _Alloc> & m, const A & key) :m_map(m), m_key(key) {}
-        ref & operator=(const B & b) {
-            m_map.emplace_replace(m_key, b);
-            return *this;
-        }
-        ref & operator=(const ref & b) { lp_assert(false); return *this; }
-        operator const B&() const {
-            auto it = m_map.m_map.find(m_key);
-            lp_assert(it != m_map.m_map.end());
-            return it->second;
-        }
-    };
-	typename std::map < A, B, _Pr, _Alloc>::iterator upper_bound(const A& k) {
-		return m_map.upper_bound(k);
-	}
-	typename std::map < A, B, _Pr, _Alloc>::const_iterator upper_bound(const A& k) const {
-		return m_map.upper_bound(k);
-	}
-	typename std::map < A, B, _Pr, _Alloc>::iterator lower_bound(const A& k) {
-		return m_map.lower_bound(k);
-	}
-	typename std::map < A, B, _Pr, _Alloc>::const_iterator lower_bound(const A& k) const {
-		return m_map.lower_bound(k);
-	}
-	typename std::map < A, B, _Pr, _Alloc>::iterator end() {
-		return m_map.end();
-	}
-	typename std::map < A, B, _Pr, _Alloc>::const_iterator end() const {
-		return m_map.end();
-	}
-	typename std::map < A, B, _Pr, _Alloc>::reverse_iterator rend() {
-		return m_map.rend();
-	}
-	typename std::map < A, B, _Pr, _Alloc>::const_reverse_iterator rend() const {
-		return m_map.rend();
-	}
-	typename std::map < A, B, _Pr, _Alloc>::iterator begin() {
-		return m_map.begin();
-	}
-	typename std::map < A, B, _Pr, _Alloc>::const_iterator begin() const {
-		return m_map.begin();
-	}
-	typename std::map < A, B, _Pr, _Alloc>::reverse_iterator rbegin() {
-		return m_map.rbegin();
-	}
-	typename std::map < A, B, _Pr, _Alloc>::const_reverse_iterator rbegin() const {
-		return m_map.rbegin();
-	}
-private:
-	void emplace_replace(const A & a, const B & b)  {
-        if (!m_stack.empty()) {
-            delta & d = m_stack.top();
-            auto it = m_map.find(a);
-            if (it == m_map.end()) {
-                d.m_new.insert(a);
-                m_map.emplace(a, b);
-            } else if (it->second != b) {
-                auto nit = d.m_new.find(a);
-                if (nit == d.m_new.end()) { // we do not have the old key
-                    auto & orig_changed= d.m_original_changed;
-                    auto itt = orig_changed.find(a);
-                    if (itt == orig_changed.end()) {
-                        orig_changed.emplace(a, it->second);
-                    } else if (itt->second == b) {
-                        orig_changed.erase(itt);
-                    }
-                }
-                it->second = b;
-            }            
-        } else { // there is no stack: just emplace or replace
-            m_map[a] = b;
-        }
-    }
-public:    
-    ref operator[] (const A & a) {
-        return ref(*this, a);
-    }
-
-    const B & operator[]( const A & a) const {
-        auto it = m_map.find(a);
-        if (it == m_map.end()) {
-            lp_assert(false);
-        }
-
-        return it->second;
-    }
-
-    bool try_get_value(const A& key, B& val) const  {
-        auto it = m_map.find(key);
-        if (it == m_map.end())
-            return false;
-        
-        val = it->second;
-        return true;
-    }
-    bool try_get_value(const A&& key, B& val) const  {
-        auto it = m_map.find(std::move(key));
-        if (it == m_map.end())
-            return false;
-        
-        val = it->second;
-        return true;
-    }
-    
-    unsigned size() const {
-        return static_cast<unsigned>(m_map.size());
-    }
-
-    bool contains(const A & key) const {
-        return m_map.find(key) != m_map.end();
-    }
-
-    bool contains(const A && key) const {
-        return m_map.find(std::move(key)) != m_map.end();
-    }
-    
-    void push() {
-        delta d;
-        //        d.m_deb_copy = m_map;
-        m_stack.push(d);
-    }
-
-    void revert() {
-        if (m_stack.empty()) return;
-
-        delta & d = m_stack.top();
-        for (auto & t : d.m_new) {
-            m_map.erase(t);
-        }
-        for (auto & t: d.m_original_changed) {
-            m_map[t.first] = t.second;
-        }
-        d.clear();
-    }
-    
-    void pop() {
-        pop(1);
-    }
-    void pop(unsigned k) {
-        while (k-- > 0) {
-            if (m_stack.empty())
-                return;
-            delta & d = m_stack.top();
-            for (auto & t : d.m_new) {
-                m_map.erase(t);
-            }
-            for (auto & t: d.m_original_changed) {
-                m_map[t.first] = t.second;
-            }
-            //            lp_assert(d.m_deb_copy == m_map);
-            m_stack.pop();
-        }
-    }
-    void erase(const typename std::map < A, B, _Pr, _Alloc>::iterator & it) {
-        erase(it->first);
-    }
-    void erase(const typename std::map < A, B, _Pr, _Alloc>::const_iterator & it) {
-        erase(it->first);
-    }
-    void erase(const typename std::map < A, B, _Pr, _Alloc>::reverse_iterator & it) {
-        erase(it->first);
-    }
-    void erase(const A & key) {
-        if (m_stack.empty()) {
-            m_map.erase(key);
-            return;
-        }
-        
-        delta & d = m_stack.top();
-        auto it = m_map.find(key);
-        if (it == m_map.end()) {
-            lp_assert(d.m_new.find(key) == d.m_new.end());
-            return;
-        }
-        auto &orig_changed = d.m_original_changed;
-        auto nit = d.m_new.find(key);
-        if (nit == d.m_new.end()) { // key is old
-            if (orig_changed.find(key) == orig_changed.end())
-                orig_changed.emplace(it->first, it->second); // need to restore
-        } else { // k is new
-            lp_assert(orig_changed.find(key) == orig_changed.end());
-            d.m_new.erase(nit);
-        }
-
-        m_map.erase(it);
-    }
-    
-    void clear() {
-        if (m_stack.empty()) {
-            m_map.clear();
-            return;
-        }
-
-        delta & d = m_stack.top();
-        auto & oc = d.m_original_changed;
-        for (auto & p : m_map) {
-            const auto & it = oc.find(p.first);
-            if (it == oc.end() && d.m_new.find(p.first) == d.m_new.end())
-                oc.emplace(p.first, p.second);
-        }
-        m_map.clear();
-    }
-
-    unsigned stack_size() const { return m_stack.size(); }
-    
-    bool empty() const { return m_map.empty(); }
-
-    const std::map<A, B, _Pr, _Alloc>& operator()() const { return m_map;}
-};
-}

src/util/lp/stacked_unordered_set.h

@@ -1,107 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-
-#pragma once
-// this class implements an unordered_set with some stack functionality
-#include <unordered_set>
-#include <set>
-#include <stack>
-namespace lp {
-
-template <typename A, 
-          typename Hash = std::hash<A>,
-          typename KeyEqual = std::equal_to<A>,
-          typename Allocator = std::allocator<A>
-          > class stacked_unordered_set {
-    struct delta {
-        std::unordered_set<A, Hash, KeyEqual> m_inserted;
-        std::unordered_set<A, Hash, KeyEqual> m_erased;
-        std::unordered_set<A, Hash, KeyEqual> m_deb_copy;
-    };
-    std::unordered_set<A, Hash, KeyEqual, Allocator> m_set;
-    std::stack<delta> m_stack;
-public:
-    void insert(const A & a) {
-        if (m_stack.empty()) {
-            m_set.insert(a);
-        } else if (m_set.find(a) == m_set.end()) {
-            m_set.insert(a);
-            size_t in_erased = m_stack.top().m_erased.erase(a);
-            if (in_erased == 0) {
-                m_stack.top().m_inserted.insert(a);
-            }
-        }
-    }
-
-    void erase(const A &a) {
-        if (m_stack.empty()) {
-            m_set.erase(a);
-            return;
-        }
-        auto erased = m_set.erase(a);
-        if (erased == 1) {
-            auto was_new = m_stack.top().m_inserted.erase(a);
-            if (was_new == 0) {
-                m_stack.top().m_erased.insert(a);
-            }
-        }
-    }
-    
-    unsigned size() const {
-        return m_set.size();
-    }
-
-    bool contains(A & key) const {
-        return m_set.find(key) != m_set.end();
-    }
-
-    bool contains(A && key) const {
-        return m_set.find(std::move(key)) != m_set.end();
-    }
-    
-    void push() {
-        delta d;
-        d.m_deb_copy = m_set;
-        m_stack.push(d);
-    }
-    
-    void pop() {
-        pop(1);
-    }
-    void pop(unsigned k) {
-        while (k-- > 0) {
-            if (m_stack.empty())
-                return;
-            delta & d = m_stack.top();
-            for (auto & t : d.m_inserted) {
-                m_set.erase(t);
-            }
-            for (auto & t : d.m_erased) {
-                m_set.insert(t);
-            }
-            lp_assert(d.m_deb_copy == m_set);
-            m_stack.pop();
-        }
-    }
-
-    const std::unordered_set<A, Hash, KeyEqual, Allocator>& operator()() const { return m_set;}
-};
-}
-

src/util/lp/static_matrix_def.h

@@ -260,10 +260,6 @@ template <typename T, typename X>    void static_matrix<T, X>::check_consistency
 
     for (auto & t : by_rows) {
         auto ic = by_cols.find(t.first);
-        if (ic == by_cols.end()){
-            //std::cout << "rows have pair (" << t.first.first <<"," << t.first.second
-            //         << "), but columns don't " << std::endl;
-        }
         lp_assert(ic != by_cols.end());
         lp_assert(t.second == ic->second);
     }
@@ -347,7 +343,6 @@ template <typename T, typename X> bool static_matrix<T, X>::is_correct() const {
         std::unordered_set<unsigned> s;
         for (auto & rc : r) {
             if (s.find(rc.m_j) != s.end()) {
-                std::cout << "found column " << rc.m_j << " twice in a row " << i << std::endl;
                 return false;
             }
             s.insert(rc.m_j);
@@ -358,7 +353,6 @@ template <typename T, typename X> bool static_matrix<T, X>::is_correct() const {
             if (rc.get_val() != get_val(m_columns[rc.m_j][rc.m_offset]))
                 return false;
             if (is_zero(rc.get_val())) {
-                std::cout << "found zero column " << rc.m_j << " in row " << i << std::endl;
                 return false;
             }
             
@@ -370,7 +364,6 @@ template <typename T, typename X> bool static_matrix<T, X>::is_correct() const {
         std::unordered_set<unsigned> s;
         for (auto & cc : c) {
             if (s.find(cc.m_i) != s.end()) {
-                std::cout << "found row " << cc.m_i << " twice in a column " << j << std::endl;
                 return false;
             }
             s.insert(cc.m_i);

src/util/lp/test_bound_analyzer.h

@@ -115,7 +115,6 @@ public :
                 }
             }
         default:
-            // std::cout << " got an upper bound with " << T_to_string(l) << "\n";
             m_implied_bounds.push_back(implied_bound(l, j, false, is_pos(m_coeffs[i]), m_row_or_term_index, strict));
         }
     }
@@ -207,7 +206,6 @@ public :
                 }
             }
         default:
-            // std::cout << " got a lower bound with " << T_to_string(l) << "\n";
             m_implied_bounds.push_back(implied_bound(l, j, true, is_pos(m_coeffs[i]), m_row_or_term_index, strict));
         }
     }

src/util/lp/var_register.h

@@ -34,6 +34,8 @@ public:
         return ret;
     }
 
+    const svector<unsigned> & vars() const { return m_local_vars_to_external; }
+    
     unsigned local_var_to_user_var(unsigned local_var) const {
         return m_local_vars_to_external[local_var];
     }

src/util/mpff.cpp

@@ -1070,7 +1070,7 @@ bool mpff_manager::is_power_of_two(mpff const & a) const {
 
 template<bool SYNCH>
 void mpff_manager::significand_core(mpff const & n, mpz_manager<SYNCH> & m, mpz & t) {
-    m.set(t, m_precision, sig(n));
+    m.set_digits(t, m_precision, sig(n));
 }
 
 void mpff_manager::significand(mpff const & n, unsynch_mpz_manager & m, mpz & t) {
@@ -1090,10 +1090,10 @@ void mpff_manager::to_mpz_core(mpff const & n, mpz_manager<SYNCH> & m, mpz & t)
         to_buffer(0, n);
         unsigned * b = m_buffers[0].c_ptr();
         shr(m_precision, b, -exp, m_precision, b);
-        m.set(t, m_precision, b);
+        m.set_digits(t, m_precision, b);
     }
     else {
-        m.set(t, m_precision, sig(n));
+        m.set_digits(t, m_precision, sig(n));
         if (exp > 0) {
             _scoped_numeral<mpz_manager<SYNCH> > p(m);
             m.set(p, 2);

src/util/mpfx.cpp

@@ -705,7 +705,7 @@ template<bool SYNCH>
 void mpfx_manager::to_mpz_core(mpfx const & n, mpz_manager<SYNCH> & m, mpz & t) {
     SASSERT(is_int(n));
     unsigned * w = words(n);
-    m.set(t, m_int_part_sz, w+m_frac_part_sz);
+    m.set_digits(t, m_int_part_sz, w+m_frac_part_sz);
     if (is_neg(n))
         m.neg(t);
 }

src/util/mpq.cpp

@@ -315,6 +315,112 @@ unsigned mpq_manager<SYNCH>::prev_power_of_two(mpq const & a) {
     return prev_power_of_two(_tmp);
 }
 
+
+template<bool SYNCH>
+template<bool SUB>
+void mpq_manager<SYNCH>::lin_arith_op(mpq const& a, mpq const& b, mpq& c, mpz& g, mpz& tmp1, mpz& tmp2, mpz& tmp3) {
+    gcd(a.m_den, b.m_den, g);                   
+    if (is_one(g)) {                            
+       mul(a.m_num, b.m_den, tmp1);             
+       mul(b.m_num, a.m_den, tmp2); 
+       if (SUB) sub(tmp1, tmp2, c.m_num); else add(tmp1, tmp2, c.m_num);
+       mul(a.m_den, b.m_den, c.m_den);          
+    }                                           
+    else {                                      
+        div(a.m_den, g, tmp3);                  
+        mul(tmp3, b.m_den, c.m_den);            
+        mul(tmp3, b.m_num, tmp2);               
+        div(b.m_den, g, tmp3);                  
+        mul(tmp3, a.m_num, tmp1);               
+        if (SUB) sub(tmp1, tmp2, tmp3); else add(tmp1, tmp2, tmp3);
+        gcd(tmp3, g, tmp1);                     
+        if (is_one(tmp1)) {                     
+            set(c.m_num, tmp3);                 
+        }                                       
+        else {                                  
+            div(tmp3, tmp1, c.m_num);           
+            div(c.m_den, tmp1, c.m_den);        
+        }                                       
+    }                                           
+}
+
+template<bool SYNCH>
+void mpq_manager<SYNCH>::rat_mul(mpq const & a, mpq const & b, mpq & c, mpz& g1, mpz& g2, mpz& tmp1, mpz& tmp2) {
+    gcd(a.m_den, b.m_num, g1);
+    gcd(a.m_num, b.m_den, g2);
+    div(a.m_num, g2, tmp1);
+    div(b.m_num, g1, tmp2);
+    mul(tmp1, tmp2, c.m_num);
+    div(b.m_den, g2, tmp1);
+    div(a.m_den, g1, tmp2);
+    mul(tmp1, tmp2, c.m_den);
+}
+
+template<bool SYNCH>
+void mpq_manager<SYNCH>::rat_mul(mpz const & a, mpq const & b, mpq & c) {
+    STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); 
+    mul(a, b.m_num, c.m_num);
+    set(c.m_den, b.m_den);
+    normalize(c);
+    STRACE("rat_mpq", tout << to_string(c) << "\n";);
+}
+
+
+template<bool SYNCH>
+void mpq_manager<SYNCH>::rat_mul(mpq const & a, mpq const & b, mpq & c) {
+    STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); 
+    if (SYNCH) {
+        mpz g1, g2, tmp1, tmp2;
+        rat_mul(a, b, c, g1, g2, tmp1, tmp2);
+        del(g1);
+        del(g2);
+        del(tmp1);
+        del(tmp2);
+    }
+    else {
+        mpz& g1 = m_n_tmp, &g2 = m_addmul_tmp.m_num, &tmp1 = m_add_tmp1, &tmp2 = m_add_tmp2;
+        rat_mul(a, b, c, g1, g2, tmp1, tmp2);
+    }
+    STRACE("rat_mpq", tout << to_string(c) << "\n";);
+}
+
+template<bool SYNCH>
+void mpq_manager<SYNCH>::rat_add(mpq const & a, mpq const & b, mpq & c) {
+    STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " + " << to_string(b) << " == ";); 
+    if (SYNCH) {
+        mpz_stack tmp1, tmp2, tmp3, g;
+        lin_arith_op<false>(a, b, c, g, tmp1, tmp2, tmp3);
+        del(tmp1);
+        del(tmp2);
+        del(tmp3);
+        del(g);
+    }
+    else {
+        mpz& g = m_n_tmp, &tmp1 = m_add_tmp1, &tmp2 = m_add_tmp2, &tmp3 = m_addmul_tmp.m_num;
+        lin_arith_op<false>(a, b, c, g, tmp1, tmp2, tmp3);
+    }
+    STRACE("rat_mpq", tout << to_string(c) << "\n";);
+}
+
+template<bool SYNCH>
+void mpq_manager<SYNCH>::rat_sub(mpq const & a, mpq const & b, mpq & c) {
+    STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " - " << to_string(b) << " == ";); 
+    if (SYNCH) {
+        mpz tmp1, tmp2, tmp3, g;
+        lin_arith_op<true>(a, b, c, g, tmp1, tmp2, tmp3);
+        del(tmp1);
+        del(tmp2);
+        del(tmp3);
+        del(g);
+    }
+    else {
+        mpz& g = m_n_tmp, &tmp1 = m_add_tmp1, &tmp2 = m_add_tmp2, &tmp3 = m_addmul_tmp.m_num;
+        lin_arith_op<true>(a, b, c, g, tmp1, tmp2, tmp3);
+    }
+    STRACE("rat_mpq", tout << to_string(c) << "\n";);
+}
+
+
 template class mpq_manager<true>;
 template class mpq_manager<false>;
 

src/util/mpq.h

@@ -74,27 +74,7 @@ class mpq_manager : public mpz_manager<SYNCH> {
         }
     }
 
-    void rat_add(mpq const & a, mpq const & b, mpq & c) {
+    void rat_add(mpq const & a, mpq const & b, mpq & c);
-        STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " + " << to_string(b) << " == ";); 
-        if (SYNCH) {
-            mpz tmp1, tmp2;
-            mul(a.m_num, b.m_den, tmp1);
-            mul(b.m_num, a.m_den, tmp2);
-            mul(a.m_den, b.m_den, c.m_den);
-            add(tmp1, tmp2, c.m_num);
-            normalize(c);
-            del(tmp1);
-            del(tmp2);
-        }
-        else {
-            mul(a.m_num, b.m_den, m_add_tmp1);
-            mul(b.m_num, a.m_den, m_add_tmp2);
-            mul(a.m_den, b.m_den, c.m_den);
-            add(m_add_tmp1, m_add_tmp2, c.m_num);
-            normalize(c);
-        }
-        STRACE("rat_mpq", tout << to_string(c) << "\n";);
-    }
 
     void rat_add(mpq const & a, mpz const & b, mpq & c) {
         STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " + " << to_string(b) << " == ";); 
@@ -115,46 +95,19 @@ class mpq_manager : public mpz_manager<SYNCH> {
         STRACE("rat_mpq", tout << to_string(c) << "\n";);
     }
 
-    void rat_sub(mpq const & a, mpq const & b, mpq & c) {
+    void rat_sub(mpq const & a, mpq const & b, mpq & c);
-        STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " - " << to_string(b) << " == ";); 
-        if (SYNCH) {
-            mpz tmp1, tmp2;
-            mul(a.m_num, b.m_den, tmp1);
-            mul(b.m_num, a.m_den, tmp2);
-            mul(a.m_den, b.m_den, c.m_den);
-            sub(tmp1, tmp2, c.m_num);
-            normalize(c);
-            del(tmp1);
-            del(tmp2);
-        }
-        else {
-            mul(a.m_num, b.m_den, m_add_tmp1);
-            mul(b.m_num, a.m_den, m_add_tmp2);
-            mul(a.m_den, b.m_den, c.m_den);
-            sub(m_add_tmp1, m_add_tmp2, c.m_num);
-            normalize(c);
-        }
-        STRACE("rat_mpq", tout << to_string(c) << "\n";);
-    }
 
-    void rat_mul(mpq const & a, mpq const & b, mpq & c) {
+    void rat_mul(mpq const & a, mpq const & b, mpq & c);
-        STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); 
-        mul(a.m_num, b.m_num, c.m_num);
-        mul(a.m_den, b.m_den, c.m_den);
-        normalize(c);
-        STRACE("rat_mpq", tout << to_string(c) << "\n";);
-    }
 
-    void rat_mul(mpz const & a, mpq const & b, mpq & c) {
+    void rat_mul(mpz const & a, mpq const & b, mpq & c);
-        STRACE("rat_mpq", tout << "[mpq] " << to_string(a) << " * " << to_string(b) << " == ";); 
-        mul(a, b.m_num, c.m_num);
-        set(c.m_den, b.m_den);
-        normalize(c);
-        STRACE("rat_mpq", tout << to_string(c) << "\n";);
-    }
 
     bool rat_lt(mpq const & a, mpq const & b);
 
+    template<bool SUB>
+    void lin_arith_op(mpq const& a, mpq const& b, mpq& c, mpz& g, mpz& tmp1, mpz& tmp2, mpz& tmp3);
+
+    void rat_mul(mpq const & a, mpq const & b, mpq & c, mpz& g1, mpz& g2, mpz& tmp1, mpz& tmp2);
+
 public:
     typedef mpq numeral;
     typedef mpq rational;
@@ -746,10 +699,12 @@ public:
         reset_denominator(a);
     }
 
-    void set(mpz & a, unsigned sz, digit_t const * digits) { mpz_manager<SYNCH>::set(a, sz, digits); }
+    void set(mpz & a, unsigned sz, digit_t const * digits) { 
+        mpz_manager<SYNCH>::set_digits(a, sz, digits); 
+    }
 
     void set(mpq & a, unsigned sz, digit_t const * digits) { 
-        mpz_manager<SYNCH>::set(a.m_num, sz, digits); 
+        mpz_manager<SYNCH>::set_digits(a.m_num, sz, digits); 
         reset_denominator(a); 
     }
 

src/util/mpz.h

@@ -64,6 +64,7 @@ class mpz_cell {
     digit_t   m_digits[0];
     friend class mpz_manager<true>;
     friend class mpz_manager<false>;
+    friend class mpz_stack;
 };
 #else
 #include<gmp.h>
@@ -72,17 +73,25 @@ class mpz_cell {
 /**
    \brief Multi-precision integer.
    
-   If m_ptr == 0, the it is a small number and the value is stored at m_val.
+   If m_kind == mpz_small, it is a small number and the value is stored in m_val.
-   Otherwise, m_val contains the sign (-1 negative, 1 positive), and m_ptr points to a mpz_cell that
+   If m_kind == mpz_ptr,   the value is stored in m_ptr and m_ptr != nullptr.
-   store the value. <<< This last statement is true only in Windows.
+                           m_val contains the sign (-1 negative, 1 positive)   
+                           under winodws, m_ptr points to a mpz_cell that store the value. 
 */
+
+enum mpz_kind { mpz_small = 0, mpz_ptr = 1};
+enum mpz_owner { mpz_self = 0, mpz_ext = 1};
+
 class mpz {
-    int        m_val; 
 #ifndef _MP_GMP
-    mpz_cell * m_ptr;
+    typedef mpz_cell mpz_type;
 #else
-    mpz_t * m_ptr;
+    typedef mpz_t mpz_type;
 #endif
+    int        m_val; 
+    unsigned   m_kind:1;
+    unsigned   m_owner:1;
+    mpz_type * m_ptr;
     friend class mpz_manager<true>;
     friend class mpz_manager<false>;
     friend class mpq_manager<true>;
@@ -90,19 +99,38 @@ class mpz {
     friend class mpq;
     friend class mpbq;
     friend class mpbq_manager;
+    friend class mpz_stack;
     mpz & operator=(mpz const & other) { UNREACHABLE(); return *this; }
 public:
-    mpz(int v):m_val(v), m_ptr(nullptr) {}
+    mpz(int v):m_val(v), m_kind(mpz_small), m_owner(mpz_self), m_ptr(nullptr) {}
-    mpz():m_val(0), m_ptr(nullptr) {}
+    mpz():m_val(0), m_kind(mpz_small), m_owner(mpz_self), m_ptr(nullptr) {}
-    mpz(mpz && other) : m_val(other.m_val), m_ptr(nullptr) {
+    mpz(mpz_type* ptr): m_val(0), m_kind(mpz_small), m_owner(mpz_ext), m_ptr(ptr) { SASSERT(ptr);}
+    mpz(mpz && other) : m_val(other.m_val), m_kind(mpz_small), m_owner(mpz_self), m_ptr(nullptr) {
         std::swap(m_ptr, other.m_ptr);
+        unsigned o = m_owner; m_owner = other.m_owner; other.m_owner = o;
+        unsigned k = m_kind; m_kind = other.m_kind; other.m_kind = k;
     }
     void swap(mpz & other) { 
         std::swap(m_val, other.m_val);
         std::swap(m_ptr, other.m_ptr);
+        unsigned o = m_owner; m_owner = other.m_owner; other.m_owner = o;
+        unsigned k = m_kind; m_kind = other.m_kind; other.m_kind = k;
     }
 };
 
+#ifndef _MP_GMP
+class mpz_stack : public mpz {
+    static const unsigned capacity = 8;
+    unsigned char m_bytes[sizeof(mpz_cell) + sizeof(digit_t) * capacity];
+public:
+    mpz_stack():mpz(reinterpret_cast<mpz_cell*>(m_bytes)) {
+        m_ptr->m_capacity = capacity;
+    }
+};
+#else
+class mpz_stack : public mpz {};
+#endif
+
 inline void swap(mpz & m1, mpz & m2) { m1.swap(m2); }
 
 template<bool SYNCH = true>
@@ -115,57 +143,56 @@ class mpz_manager {
 
 #ifndef _MP_GMP
     unsigned                m_init_cell_capacity;
-    mpz_cell *              m_tmp[2];
-    mpz_cell *              m_arg[2];
     mpz                     m_int_min;
     
-    static unsigned cell_size(unsigned capacity) { return sizeof(mpz_cell) + sizeof(digit_t) * capacity; }
+    static unsigned cell_size(unsigned capacity) { 
+        return sizeof(mpz_cell) + sizeof(digit_t) * capacity; 
+    }
 
     mpz_cell * allocate(unsigned capacity) {
         SASSERT(capacity >= m_init_cell_capacity);
+        MPZ_BEGIN_CRITICAL();
         mpz_cell * cell  = reinterpret_cast<mpz_cell *>(m_allocator.allocate(cell_size(capacity)));
+        MPZ_END_CRITICAL();
         cell->m_capacity = capacity;
         return cell;
     }
+
+    void deallocate(mpz& n) {
+        if (n.m_ptr) {
+            deallocate(n.m_owner == mpz_self, n.m_ptr);
+            n.m_ptr = nullptr;
+            n.m_kind = mpz_small;
+        }
+    }
     
     // make sure that n is a big number and has capacity equal to at least c.
     void allocate_if_needed(mpz & n, unsigned c) {
-        if (c < m_init_cell_capacity)
+        c = std::max(c, m_init_cell_capacity);
-            c = m_init_cell_capacity;
+        if (n.m_ptr == nullptr || capacity(n) < c) {
-        if (is_small(n)) {
+            deallocate(n);
             n.m_val             = 1;
+            n.m_kind            = mpz_ptr;
+            n.m_owner           = mpz_self;
             n.m_ptr             = allocate(c);
-            n.m_ptr->m_capacity = c;
-        }
-        else if (capacity(n) < c) {
-            deallocate(n.m_ptr);
-            n.m_val             = 1;
-            n.m_ptr             = allocate(c);
-            n.m_ptr->m_capacity = c;
         }
     }
 
-    void deallocate(mpz_cell * ptr) { 
+    void deallocate(bool is_heap, mpz_cell * ptr) { 
-        m_allocator.deallocate(cell_size(ptr->m_capacity), ptr); 
+        if (is_heap) {
-    }
+            MPZ_BEGIN_CRITICAL();
-
+            m_allocator.deallocate(cell_size(ptr->m_capacity), ptr); 
-    /**
+            MPZ_END_CRITICAL();
-      \brief Make sure that m_tmp[IDX] can hold the given number of digits
+        }
-    */
-    template<int IDX>
-    void ensure_tmp_capacity(unsigned capacity) {
-        if (m_tmp[IDX]->m_capacity >= capacity)
-            return;
-        deallocate(m_tmp[IDX]);
-        unsigned new_capacity = (3 * capacity + 1) >> 1;
-        m_tmp[IDX] = allocate(new_capacity);
-        SASSERT(m_tmp[IDX]->m_capacity >= capacity);
     }
     
     // Expand capacity of a while preserving its content.
     void ensure_capacity(mpz & a, unsigned sz);
 
     void normalize(mpz & a);
+
+    void clear(mpz& n) { }
+
 #else
     // GMP code
     mpz_t     m_tmp, m_tmp2;
@@ -175,22 +202,34 @@ class mpz_manager {
     mpz_t     m_int64_max;
     mpz_t     m_int64_min;
 
-    mpz_t * allocate() {
+    mpz_t * allocate() {        
+        MPZ_BEGIN_CRITICAL();
         mpz_t * cell = reinterpret_cast<mpz_t*>(m_allocator.allocate(sizeof(mpz_t)));
+        MPZ_END_CRITICAL();
         mpz_init(*cell);
         return cell;
     }
 
-    void deallocate(mpz_t * ptr) { mpz_clear(*ptr); m_allocator.deallocate(sizeof(mpz_t), ptr); }
+    void deallocate(bool is_heap, mpz_t * ptr) { 
+        mpz_clear(*ptr); 
+        if (is_heap) {
+            MPZ_BEGIN_CRITICAL();
+            m_allocator.deallocate(sizeof(mpz_t), ptr); 
+            MPZ_END_CRITICAL();
+        }
+    }
+
+    void clear(mpz& n) { if (n.m_ptr) mpz_clear(n.m_ptr); }
 #endif
+
+
     mpz                     m_two64;
 
     /**
-       \brief Set \c a with the value stored at m_tmp[IDX], and the given sign.
+       \brief Set \c a with the value stored at src, and the given sign.
-       \c sz is an overapproximation of the size of the number stored at \c tmp.
+       \c sz is an overapproximation of the size of the number stored at \c src.
     */
-    template<int IDX>
+    void set(mpz_cell& src, mpz & a, int sign, unsigned sz); 
-    void set(mpz & a, int sign, unsigned sz);
 
     static int64_t i64(mpz const & a) { return static_cast<int64_t>(a.m_val); }
 
@@ -198,19 +237,19 @@ class mpz_manager {
 
     void set_i64(mpz & c, int64_t v) {
         if (v >= INT_MIN && v <= INT_MAX) {
-            del(c);
             c.m_val = static_cast<int>(v); 
+            c.m_kind = mpz_small;
         }
         else {
-            MPZ_BEGIN_CRITICAL();
             set_big_i64(c, v);
-            MPZ_END_CRITICAL();
         }
     }
 
     void set_big_ui64(mpz & c, uint64_t v);
 
+
 #ifndef _MP_GMP
+
     static unsigned capacity(mpz const & c) { return c.m_ptr->m_capacity; }
 
     static unsigned size(mpz const & c) { return c.m_ptr->m_size; }
@@ -241,16 +280,28 @@ class mpz_manager {
             return ((static_cast<uint64_t>(digits(a)[1]) << 32) | (static_cast<uint64_t>(digits(a)[0])));
     }
 
-    template<int IDX>
+    class sign_cell {
-    void get_sign_cell(mpz const & a, int & sign, mpz_cell * & cell) {
+        static const unsigned capacity = 2;
+        unsigned char m_bytes[sizeof(mpz_cell) + sizeof(digit_t) * capacity];
+        mpz m_local;
+        mpz const& m_a;
+        int m_sign;
+        mpz_cell* m_cell;
+    public:
+        sign_cell(mpz_manager& m, mpz const& a);
+        int             sign() { return m_sign; }
+        mpz_cell const* cell() { return m_cell; }
+    };
+
+    void get_sign_cell(mpz const & a, int & sign, mpz_cell * & cell, mpz_cell* reserve) {
         if (is_small(a)) {
             if (a.m_val == INT_MIN) {
                 sign = -1;
                 cell = m_int_min.m_ptr;
             }
             else {
-                cell = m_arg[IDX];
+                cell = reserve;
-                SASSERT(cell->m_size == 1);
+                cell->m_size = 1;
                 if (a.m_val < 0) {
                     sign = -1;
                     cell->m_digits[0] = -a.m_val;
@@ -266,26 +317,29 @@ class mpz_manager {
             cell = a.m_ptr;
         }
     }
+
 #else
     // GMP code
 
-    template<int IDX>
+    class ensure_mpz_t {
-    void get_arg(mpz const & a, mpz_t * & result) {
+        mpz_t m_local;
-        if (is_small(a)) {
+        mpz_t* m_result;
-            result = m_arg[IDX];
+    public:
-            mpz_set_si(*result, a.m_val);
+        ensure_mpz_t(mpz const& a);
-        }
+        ~ensure_mpz_t();
-        else {
+        mpz_t& operator()() { return *m_result; }
-            result = a.m_ptr;
+    };
-        }
-    }
     
     void mk_big(mpz & a) {
-        if (a.m_ptr == 0) {
+        if (a.m_ptr == null) {
             a.m_val = 0;
             a.m_ptr = allocate();
+            a.m_owner = mpz_self;
         }
+        a.m_kind = mpz_ptr;
     }
+
+
 #endif 
 
 #ifndef _MP_GMP
@@ -333,245 +387,49 @@ public:
 
     ~mpz_manager();
 
-    static bool is_small(mpz const & a) { return a.m_ptr == nullptr; }
+    static bool is_small(mpz const & a) { return a.m_kind == mpz_small; }
 
     static mpz mk_z(int val) { return mpz(val); }
     
-    void del(mpz & a) { 
+    void del(mpz & a);
-        if (a.m_ptr != nullptr) {
-            MPZ_BEGIN_CRITICAL();
-            deallocate(a.m_ptr); 
-            MPZ_END_CRITICAL();
-            a.m_ptr = nullptr;
-        } 
-    }
     
-    void add(mpz const & a, mpz const & b, mpz & c) {
+    void add(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz] " << to_string(a) << " + " << to_string(b) << " == ";); 
-        if (is_small(a) && is_small(b)) {
-            set_i64(c, i64(a) + i64(b));
-        }
-        else {
-            MPZ_BEGIN_CRITICAL();
-            big_add(a, b, c);
-            MPZ_END_CRITICAL();
-        }
-        STRACE("mpz", tout << to_string(c) << "\n";);
-    }
 
-    void sub(mpz const & a, mpz const & b, mpz & c) {
+    void sub(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz] " << to_string(a) << " - " << to_string(b) << " == ";); 
-        if (is_small(a) && is_small(b)) {
-            set_i64(c, i64(a) - i64(b));
-        }
-        else {
-            MPZ_BEGIN_CRITICAL();
-            big_sub(a, b, c);
-            MPZ_END_CRITICAL();
-        }
-        STRACE("mpz", tout << to_string(c) << "\n";);
-    }
     
     void inc(mpz & a) { add(a, mpz(1), a); }
 
     void dec(mpz & a) { add(a, mpz(-1), a); }
 
-    void mul(mpz const & a, mpz const & b, mpz & c) {
+    void mul(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz] " << to_string(a) << " * " << to_string(b) << " == ";); 
-        if (is_small(a) && is_small(b)) {
-            set_i64(c, i64(a) * i64(b));
-        }
-        else {
-            MPZ_BEGIN_CRITICAL();
-            big_mul(a, b, c);
-            MPZ_END_CRITICAL();
-        }
-        STRACE("mpz", tout << to_string(c) << "\n";);
-    }
 
     // d <- a + b*c
-    void addmul(mpz const & a, mpz const & b, mpz const & c, mpz & d) {
+    void addmul(mpz const & a, mpz const & b, mpz const & c, mpz & d);
-        if (is_one(b)) {
-            add(a, c, d);
-        }
-        else if (is_minus_one(b)) {
-            sub(a, c, d);
-        }
-        else {
-            mpz tmp;
-            mul(b,c,tmp);
-            add(a,tmp,d);
-            del(tmp);
-        }
-    }
-
 
     // d <- a - b*c
-    void submul(mpz const & a, mpz const & b, mpz const & c, mpz & d) {
+    void submul(mpz const & a, mpz const & b, mpz const & c, mpz & d);
-        if (is_one(b)) {
-            sub(a, c, d);
-        }
-        else if (is_minus_one(b)) {
-            add(a, c, d);
-        }
-        else {
-            mpz tmp;
-            mul(b,c,tmp);
-            sub(a,tmp,d);
-            del(tmp);
-        }
-    }
 
+    void machine_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r);
 
-    void machine_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r) {
+    void machine_div(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz-ext] divrem(" << to_string(a) << ",  " << to_string(b) << ") == ";); 
-        if (is_small(a) && is_small(b)) {
-            int64_t _a = i64(a);
-            int64_t _b = i64(b);
-            set_i64(q, _a / _b);
-            set_i64(r, _a % _b);
-        }
-        else {
-            MPZ_BEGIN_CRITICAL();
-            big_div_rem(a, b, q, r);
-            MPZ_END_CRITICAL();
-        }
-        STRACE("mpz", tout << "(" << to_string(q) << ", " << to_string(r) << ")\n";);
-    }
 
-    void machine_div(mpz const & a, mpz const & b, mpz & c) {
+    void rem(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz-ext] machine-div(" << to_string(a) << ",  " << to_string(b) << ") == ";); 
-        if (is_small(a) && is_small(b)) {
-            set_i64(c, i64(a) / i64(b));
-        }
-        else {
-            MPZ_BEGIN_CRITICAL();
-            big_div(a, b, c);
-            MPZ_END_CRITICAL();
-        }
-        STRACE("mpz", tout << to_string(c) << "\n";);
-    }
 
-    void rem(mpz const & a, mpz const & b, mpz & c) {
+    void div_gcd(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz-ext] rem(" << to_string(a) << ",  " << to_string(b) << ") == ";); 
-        if (is_small(a) && is_small(b)) {
-            set_i64(c, i64(a) % i64(b));
-        }
-        else {
-            MPZ_BEGIN_CRITICAL();
-            big_rem(a, b, c);
-            MPZ_END_CRITICAL();
-        }
-        STRACE("mpz", tout << to_string(c) << "\n";);
-    }
 
-    void div(mpz const & a, mpz const & b, mpz & c) {
+    void div(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz-ext] div(" << to_string(a) << ",  " << to_string(b) << ") == ";); 
-        if (is_neg(a)) {
-            mpz tmp;
-            machine_div_rem(a, b, c, tmp);
-            if (!is_zero(tmp)) {
-                if (is_neg(b))
-                    add(c, mk_z(1), c);
-                else
-                    sub(c, mk_z(1), c);
-            }
-            del(tmp);
-        }
-        else {
-            machine_div(a, b, c);
-        }
-        STRACE("mpz", tout << to_string(c) << "\n";);
-    }
 
-    void mod(mpz const & a, mpz const & b, mpz & c) {
+    void mod(mpz const & a, mpz const & b, mpz & c);
-        STRACE("mpz", tout << "[mpz-ext] mod(" << to_string(a) << ",  " << to_string(b) << ") == ";); 
-        rem(a, b, c);
-        if (is_neg(c)) {
-            if (is_pos(b))
-                add(c, b, c);
-            else
-                sub(c, b, c);
-        }
-        STRACE("mpz", tout << to_string(c) << "\n";);
-    }
 
-    void neg(mpz & a) {
+    void neg(mpz & a);
-        STRACE("mpz", tout << "[mpz] 0 - " << to_string(a) << " == ";); 
-        if (is_small(a) && a.m_val == INT_MIN) {
-            // neg(INT_MIN) is not a small int
-            MPZ_BEGIN_CRITICAL();
-            set_big_i64(a, - static_cast<long long>(INT_MIN)); 
-            MPZ_END_CRITICAL();
-            return;
-        }
-#ifndef _MP_GMP
-        a.m_val = -a.m_val;
-#else
-        if (is_small(a)) {
-            a.m_val = -a.m_val;
-        }
-        else {
-            mpz_neg(*a.m_ptr, *a.m_ptr);
-        }
-#endif
-        STRACE("mpz", tout << to_string(a) << "\n";); 
-    }
 
-    void abs(mpz & a) {
+    void abs(mpz & a);
-        if (is_small(a)) {
-            if (a.m_val < 0) {
-                if (a.m_val == INT_MIN) {
-                    // abs(INT_MIN) is not a small int
-                    MPZ_BEGIN_CRITICAL();
-                    set_big_i64(a, - static_cast<long long>(INT_MIN)); 
-                    MPZ_END_CRITICAL();
-                }
-                else
-                    a.m_val = -a.m_val;
-            }
-        }
-        else {
-#ifndef _MP_GMP
-            a.m_val = 1;
-#else
-            mpz_abs(*a.m_ptr, *a.m_ptr);
-#endif
-        }
-    }
 
-    static bool is_pos(mpz const & a) { 
+    static bool is_pos(mpz const & a) { return sign(a) > 0; }
-#ifndef _MP_GMP
-        return a.m_val > 0; 
-#else
-        if (is_small(a))
-            return a.m_val > 0;
-        else
-            return mpz_sgn(*a.m_ptr) > 0;
-#endif
-    }
 
-    static bool is_neg(mpz const & a) { 
+    static bool is_neg(mpz const & a) { return sign(a) < 0; }
-#ifndef _MP_GMP
-        return a.m_val < 0; 
-#else
-        if (is_small(a))
-            return a.m_val < 0;
-        else
-            return mpz_sgn(*a.m_ptr) < 0;
-#endif
-    }
     
-    static bool is_zero(mpz const & a) { 
+    static bool is_zero(mpz const & a) { return sign(a) == 0; }
-#ifndef _MP_GMP
-        return a.m_val == 0; 
-#else
-        if (is_small(a))
-            return a.m_val == 0;
-        else
-            return mpz_sgn(*a.m_ptr) == 0;
-#endif
-    }
 
     static int sign(mpz const & a) {
 #ifndef _MP_GMP
@@ -593,10 +451,7 @@ public:
             return a.m_val == b.m_val;
         }
         else {
-            MPZ_BEGIN_CRITICAL();
+            return big_compare(a, b) == 0;
-            bool res = big_compare(a, b) == 0;
-            MPZ_END_CRITICAL();
-            return res;
         }
     }
 
@@ -614,10 +469,7 @@ public:
             return a.m_val < b.m_val;
         }
         else {
-            MPZ_BEGIN_CRITICAL();
+            return big_compare(a, b) < 0;
-            bool res = big_compare(a, b) < 0;
-            MPZ_END_CRITICAL();
-            return res;
         }
     }
 
@@ -661,25 +513,24 @@ public:
 
     void set(mpz & target, mpz const & source) {
         if (is_small(source)) {
-            del(target);
             target.m_val = source.m_val;
+            target.m_kind = mpz_small;
         }
         else {
-            MPZ_BEGIN_CRITICAL();
             big_set(target, source);
-            MPZ_END_CRITICAL();
         }
     }
 
     void set(mpz & target, mpz && source) {
-        del(target);
         target.m_val = source.m_val;
         std::swap(target.m_ptr, source.m_ptr);
+        auto o = target.m_owner; target.m_owner = source.m_owner; source.m_owner = o;
+        auto k = target.m_kind; target.m_kind = source.m_kind; source.m_kind = k;
     }
 
     void set(mpz & a, int val) {
-        del(a);
         a.m_val = val;
+        a.m_kind = mpz_small;
     }
 
     void set(mpz & a, unsigned val) {
@@ -697,17 +548,15 @@ public:
 
     void set(mpz & a, uint64_t val) {
         if (val < INT_MAX) {
-            del(a);
             a.m_val = static_cast<int>(val);
+            a.m_kind = mpz_small;
         }
         else {
-            MPZ_BEGIN_CRITICAL();
             set_big_ui64(a, val);
-            MPZ_END_CRITICAL();
         }
     }
 
-    void set(mpz & target, unsigned sz, digit_t const * digits);
+    void set_digits(mpz & target, unsigned sz, digit_t const * digits);
 
     mpz dup(const mpz & source) {
         mpz temp;
@@ -716,13 +565,15 @@ public:
     }
 
     void reset(mpz & a) {
-        del(a);
         a.m_val = 0;
+        a.m_kind = mpz_small;
     }
 
     void swap(mpz & a, mpz & b) {
         std::swap(a.m_val, b.m_val);
         std::swap(a.m_ptr, b.m_ptr);
+        auto o = a.m_owner; a.m_owner = b.m_owner; b.m_owner = o;
+        auto k = a.m_kind; a.m_kind = b.m_kind; b.m_kind = k;
     }
 
     bool is_uint64(mpz const & a) const;