remove warnings in scaler and use m_cut_solver_cycle_on_var

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

detect slow propagations

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

fiddle with the stop conditions

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

get rid of constraint->m_predecessors, fix a bug in push/pop with lemmas

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

clean detection of stale lemmas in pop

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

add constraints lazily to cut_solver

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

refactor some of cut_solver classes into include files

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

prepare to index constraint from 0 to to n - 1, where n is the number of constraints

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

prepare for constraint priority

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

use priorities in active_set

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

remove unnecesessary parameters

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

speedup bound propagations

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

restore tactics

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

speedup bound propagation by avoiding some calls to propagate_constraint_only_one_unlim

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

fixes by Nikolaj

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

fix print lp_core_solver

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

work on gomory test, subs terms indices correctly

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

correct const_iterator for lar_term

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

improve static_matrix with iterators

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

make row_strip a struct

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

move row_strip outside of static_matrix

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

add const_iterator to row_strip

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

remove the hierarchy of iterators - use std::iterators

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

adding gcd_test stats and taking care of for iterators

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

restore qflia_tactic.cpp

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

run gcd_test according to settings()

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

experiment with picking a narrow or random branch

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

src/smt/theory_lra.cpp

@@ -307,7 +307,7 @@ class theory_lra::imp {
         m_solver->settings().bound_propagation() = BP_NONE != propagation_mode();
         m_solver->set_track_pivoted_rows(lp.bprop_on_pivoted_rows());
         m_solver->settings().m_int_branch_cut_gomory_threshold = ctx().get_fparams().m_arith_branch_cut_ratio;
-        m_solver->settings().m_int_branch_cut_solver = std::max(4u, ctx().get_fparams().m_arith_branch_cut_ratio);
+        m_solver->settings().m_int_branch_cut_solver = std::max(8u, ctx().get_fparams().m_arith_branch_cut_ratio);
         m_solver->settings().m_run_gcd_test = ctx().get_fparams().m_arith_gcd_test;
         m_solver->settings().set_random_seed(ctx().get_fparams().m_random_seed);
         //m_solver->settings().set_ostream(0);
@@ -672,6 +672,10 @@ class theory_lra::imp {
     }
 
     theory_var internalize_def(app* term, scoped_internalize_state& st) {
+        if (ctx().e_internalized(term)) {
+            IF_VERBOSE(0, verbose_stream() << "repeated term\n";);
+            return mk_var(term, false);
+        }
         linearize_term(term, st);
         if (is_unit_var(st)) {
             return st.vars()[0];
@@ -1252,17 +1256,21 @@ public:
     }
 
     // create a bound atom representing term <= k
-    app_ref mk_bound(lp::lar_term const& term, rational const& k) {
+    app_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound) {
         app_ref t = mk_term(term, k.is_int());
-        app_ref atom(a.mk_le(t, a.mk_numeral(k, k.is_int())), m);
+        app_ref atom(m);
+        if (lower_bound) {
+            atom = a.mk_ge(t, a.mk_numeral(k, k.is_int()));
+        }
+        else {
+            atom = a.mk_le(t, a.mk_numeral(k, k.is_int()));
+        }
         expr_ref atom1(m);
         proof_ref atomp(m);
         ctx().get_rewriter()(atom, atom1, atomp);
         atom = to_app(atom1);
         TRACE("arith", tout << atom << "\n";
-              m_solver->print_term(term, tout << "bound atom: "); tout << " <= " << k << "\n";
+              m_solver->print_term(term, tout << "bound atom: "); tout << " <= " << k << "\n";);
-              display(tout);
-              );
         ctx().internalize(atom, true);
         ctx().mark_as_relevant(atom.get());
         return atom;
@@ -1276,12 +1284,15 @@ public:
         lp::lar_term term;
         lp::mpq k;
         lp::explanation ex; // TBD, this should be streamlined accross different explanations
-        switch(m_lia->check(term, k, ex)) {
+        bool upper;
+        switch(m_lia->check(term, k, ex, upper)) {
         case lp::lia_move::ok:
             return l_true;
         case lp::lia_move::branch: {
-            (void)mk_bound(term, k);
+            app_ref b = mk_bound(term, k, upper);
+            // branch on term >= k + 1
             // branch on term <= k
+            // TBD: ctx().force_phase(ctx().get_literal(b));
             // at this point we have a new unassigned atom that the 
             // SAT core assigns a value to
             return l_false;
@@ -1289,7 +1300,7 @@ public:
         case lp::lia_move::cut: {
             ++m_stats.m_gomory_cuts;
             // m_explanation implies term <= k
-            app_ref b = mk_bound(term, k);
+            app_ref b = mk_bound(term, k, upper);
             m_eqs.reset();
             m_core.reset();
             m_params.reset();
@@ -1298,7 +1309,11 @@ public:
                     set_evidence(ev.second);
                 }
             }
-            assign(literal(ctx().get_bool_var(b), false));
+            literal lit(ctx().get_bool_var(b), false);
+            TRACE("arith", 
+                  ctx().display_lemma_as_smt_problem(tout << "new cut:\n", m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), lit);
+                  display(tout););
+            assign(lit);
             return l_false;
         }
         case lp::lia_move::conflict:
@@ -2725,6 +2740,8 @@ public:
         if (m_solver) {
             m_solver->print_constraints(out);
             m_solver->print_terms(out);
+            // auto pp = lp ::core_solver_pretty_printer<lp::mpq, lp::impq>(m_solver->m_mpq_lar_core_solver.m_r_solver, out);
+            // pp.print();
         }
         unsigned nv = th.get_num_vars();
         for (unsigned v = 0; v < nv; ++v) {
@@ -2793,6 +2810,8 @@ public:
         st.update("cut_solver-true", m_solver->settings().st().m_cut_solver_true);
         st.update("cut_solver-false", m_solver->settings().st().m_cut_solver_false);
         st.update("cut_solver-undef", m_solver->settings().st().m_cut_solver_undef);
+        st.update("gcd_calls", m_solver->settings().st().m_gcd_calls);
+        st.update("gcd_conflict", m_solver->settings().st().m_gcd_conflicts);
     }        
 };
     

src/test/lp/gomory_test.h

@@ -0,0 +1,236 @@
+namespace lp {
+struct gomory_test {
+    gomory_test(
+        std::function<std::string (unsigned)> name_function_p,
+        std::function<mpq (unsigned)> get_value_p,
+        std::function<bool (unsigned)> at_low_p,
+        std::function<bool (unsigned)> at_upper_p,
+        std::function<impq (unsigned) > lower_bound_p,
+        std::function<impq (unsigned) > upper_bound_p,
+        std::function<unsigned (unsigned) > column_lower_bound_constraint_p,
+        std::function<unsigned (unsigned) > column_upper_bound_constraint_p
+                ) :
+        m_name_function(name_function_p),
+        get_value(get_value_p),
+        at_low(at_low_p),
+        at_upper(at_upper_p),
+        lower_bound(lower_bound_p),
+        upper_bound(upper_bound_p),
+        column_lower_bound_constraint(column_lower_bound_constraint_p),
+        column_upper_bound_constraint(column_upper_bound_constraint_p)
+    {}
+
+    std::function<std::string (unsigned)> m_name_function;
+    std::function<mpq (unsigned)> get_value;
+    std::function<bool (unsigned)> at_low;
+    std::function<bool (unsigned)> at_upper;
+    std::function<impq (unsigned) > lower_bound;
+    std::function<impq (unsigned) > upper_bound;
+    std::function<unsigned (unsigned) > column_lower_bound_constraint;
+    std::function<unsigned (unsigned) > column_upper_bound_constraint;
+    
+    bool is_real(unsigned) { return false; } // todo: test real case
+    void print_row(std::ostream& out, vector<std::pair<mpq, unsigned>> & row ) {
+        bool first = true;
+        for (const auto & it : row) {
+            auto val = it.first;
+            if (first) {
+                first = false;
+            } else {
+                if (numeric_traits<mpq>::is_pos(val)) {
+                    out << " + ";
+                } else {
+                    out << " - ";
+                    val = -val;
+                }
+            }
+            if (val == -numeric_traits<mpq>::one())
+                out << " - ";
+            else if (val != numeric_traits<mpq>::one())
+                out << T_to_string(val);
+        
+            out << m_name_function(it.second);
+        }
+    }
+
+    void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, const mpq& f_0, const mpq& one_minus_f_0) {
+        TRACE("gomory_cut_detail_real", tout << "real\n";);
+        mpq new_a;
+        if (at_low(x_j)) {
+            if (a.is_pos()) {
+                new_a  =  a / (1 - f_0);
+            }
+            else {
+                new_a  =  a / f_0;
+                new_a.neg();
+            }
+            k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation than
+            // k += lower_bound(x_j).x * new_a;  
+            expl.push_justification(column_lower_bound_constraint(x_j), new_a);
+        }
+        else {
+            lp_assert(at_upper(x_j));
+            if (a.is_pos()) {
+                new_a =   a / f_0; 
+                new_a.neg(); // the upper terms are inverted.
+            }
+            else {
+                new_a =   a / (mpq(1) - f_0); 
+            }
+            k.addmul(new_a, upper_bound(x_j).x); //  k += upper_bound(x_j).x * new_a; 
+            expl.push_justification(column_upper_bound_constraint(x_j), new_a);
+        }
+        TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << k << "\n";);
+        pol.add_monomial(new_a, x_j);
+    }
+    
+    void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) {
+        lp_assert(is_int(x_j));
+        lp_assert(!a.is_int());
+             lp_assert(f_0 > zero_of_type<mpq>() && f_0 < one_of_type<mpq>());
+        mpq f_j =  int_solver::fractional_part(a);
+        TRACE("gomory_cut_detail", 
+              tout << a << " x_j = " << x_j << ", k = " << k << "\n";
+              tout << "f_j: " << f_j << "\n";
+              tout << "f_0: " << f_0 << "\n";
+              tout << "1 - f_0: " << one_minus_f_0 << "\n";
+              tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl;
+              );
+        lp_assert (!f_j.is_zero());
+        mpq new_a;
+        if (at_low(x_j)) {
+            if (f_j <= one_minus_f_0) {
+                new_a = f_j / one_minus_f_0;
+            }
+            else {
+                new_a = (1 - f_j) / f_0;
+            }
+            k.addmul(new_a, lower_bound(x_j).x);
+            expl.push_justification(column_lower_bound_constraint(x_j), new_a);
+        }
+        else {
+            lp_assert(at_upper(x_j));
+            if (f_j <= f_0) {
+                new_a = f_j / f_0;
+            }
+            else {
+                new_a = (mpq(1) - f_j) / (one_minus_f_0);
+            }
+            new_a.neg(); // the upper terms are inverted
+            k.addmul(new_a, upper_bound(x_j).x);
+            expl.push_justification(column_upper_bound_constraint(x_j), new_a);
+        }
+        TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";);
+        t.add_monomial(new_a, x_j);
+        lcm_den = lcm(lcm_den, denominator(new_a));
+    }
+
+
+    void report_conflict_from_gomory_cut(mpq &k) {
+        lp_assert(false);
+    }
+
+    void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq &lcm_den) {
+        lp_assert(!t.is_empty());
+        auto pol = t.coeffs_as_vector();
+        t.clear();
+        if (pol.size() == 1) {
+            TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
+            unsigned v = pol[0].second;
+            lp_assert(is_int(v));
+            const mpq& a = pol[0].first;
+            k /= a;
+            if (a.is_pos()) { // we have av >= k
+                if (!k.is_int())
+                    k = ceil(k);
+                // switch size
+                t.add_monomial(- mpq(1), v);
+                k.neg();
+            } else {
+                if (!k.is_int())
+                    k = floor(k);
+                t.add_monomial(mpq(1), v);
+            }
+        } else {
+            TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;);
+            lcm_den = lcm(lcm_den, denominator(k));
+            TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n";
+                  for (unsigned i = 0; i < pol.size(); i++) {
+                      tout << pol[i].first << " " << pol[i].second << "\n";
+                  }
+                  tout << "k: " << k << "\n";);
+            lp_assert(lcm_den.is_pos());
+            if (!lcm_den.is_one()) {
+                // normalize coefficients of integer parameters to be integers.
+                for (auto & pi: pol) {
+                    pi.first *= lcm_den;
+                    SASSERT(!is_int(pi.second) || pi.first.is_int());
+                }
+                k *= lcm_den;
+            }
+            TRACE("gomory_cut_detail", tout << "after *lcm\n";
+                  for (unsigned i = 0; i < pol.size(); i++) {
+                      tout << pol[i].first << " * v" << pol[i].second << "\n";
+                  }
+                  tout << "k: " << k << "\n";);
+            
+            // negate everything to return -pol <= -k
+            for (const auto & pi: pol)
+                t.add_monomial(-pi.first, pi.second);
+            k.neg();
+        }
+        TRACE("gomory_cut_detail", tout << "k = " << k << std::endl;);
+        lp_assert(k.is_int());
+    }
+
+    void print_term(lar_term & t, std::ostream & out) {
+        lp_assert(is_zero(t.m_v));
+        vector<std::pair<mpq, unsigned>>  row;
+        for (auto p : t.m_coeffs)
+            row.push_back(std::make_pair(p.second, p.first));
+        print_row(out, row);
+    }
+    
+    void mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, vector<std::pair<mpq, unsigned>> & row) {
+        enable_trace("gomory_cut");
+        enable_trace("gomory_cut_detail");
+
+        TRACE("gomory_cut",
+              tout << "applying cut at:\n"; print_row(tout, row); 
+              tout << std::endl << "inf_col = " << inf_col << std::endl;
+              );
+        
+        // gomory will be   t >= k
+        k = 1;
+        mpq lcm_den(1);
+        unsigned x_j;
+        mpq a;
+        bool some_int_columns = false;
+        mpq f_0  = int_solver::fractional_part(get_value(inf_col));
+        mpq one_min_f_0 = 1 - f_0;
+        for ( auto pp : row) {
+            a = pp.first;
+            x_j = pp.second;
+            if (x_j == inf_col)
+                continue;
+            // make the format compatible with the format used in: Integrating Simplex with DPLL(T)
+            a.neg();  
+            if (is_real(x_j)) 
+                real_case_in_gomory_cut(a, x_j, k, t, expl, f_0, one_min_f_0);
+            else {
+                if (a.is_int()) continue; // f_j will be zero and no monomial will be added
+                some_int_columns = true;
+                int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, f_0, one_min_f_0);
+            }
+        }
+
+        if (t.is_empty())
+            return report_conflict_from_gomory_cut(k);
+        if (some_int_columns)
+            adjust_term_and_k_for_some_ints_case_gomory(t, k, lcm_den);
+
+        TRACE("gomory_cut", tout<<"new cut :"; print_term(t, tout); tout << " >= " << k << std::endl;);
+
+    }
+};
+}

src/test/lp/lp.cpp

@@ -50,6 +50,8 @@ Revision History:
 #include "util/lp/integer_domain.h"
 #include "util/lp/stacked_map.h"
 #include <cstdlib>
+#include "test/lp/gomory_test.h"
+
 namespace lp {
 unsigned seed = 1;
 
@@ -1883,6 +1885,7 @@ void test_replace_column() {
 
 
 void setup_args_parser(argument_parser & parser) {
+    parser.add_option_with_help_string("-gomory", "gomory");
     parser.add_option_with_help_string("-intd", "test integer_domain");
     parser.add_option_with_help_string("-xyz_sample", "run a small interactive scenario");
     parser.add_option_with_after_string_with_help("--density", "the percentage of non-zeroes in the matrix below which it is not dense");
@@ -3255,6 +3258,153 @@ void test_resolve(cut_solver& cs, unsigned constraint_index, unsigned i0)  {
  }
 
 
+void test_gomory_cut_0() {
+        gomory_test g(
+            [](unsigned j) { return "v" + T_to_string(j);} // name_function_p
+            ,
+            [](unsigned j) { //get_value_p
+                if (j == 1)
+                    return mpq(2730, 1727);
+                if (j == 2)
+                    return zero_of_type<mpq>();
+                if (j == 3) return mpq(3);
+                lp_assert(false);
+                return zero_of_type<mpq>();
+            },
+            [](unsigned j) { // at_low_p
+                if (j == 1)
+                    return false;
+                if (j == 2)
+                    return true;
+                if (j == 3)
+                    return true;
+                lp_assert(false);
+                return false;
+            },
+            [](unsigned j) { // at_upper
+                if (j == 1)
+                    return false;
+                if (j == 2)
+                    return true;
+                if (j == 3)
+                    return false;
+                lp_assert(false);
+                return false;
+            },
+            [](unsigned j) { // lower_bound
+                if (j == 1) {
+                    lp_assert(false); //unlimited from below
+                    return 0;
+                }
+                if (j == 2)
+                    return 0;
+                if (j == 3)
+                    return 3;
+                lp_assert(false);
+                return 0;
+            },
+            [](unsigned j) { // upper
+                if (j == 1) {
+                    lp_assert(false); //unlimited from above
+                    return 0;
+                }
+                if (j == 2)
+                    return 0;
+                if (j == 3)
+                    return 10;
+                lp_assert(false);
+                return 0;
+            },
+            [] (unsigned) { return 0; },
+            [] (unsigned) { return 0; }
+                      );
+        lar_term t;
+        mpq k;
+        explanation expl;
+        unsigned inf_col = 1;
+        vector<std::pair<mpq, unsigned>> row;
+        row.push_back(std::make_pair(mpq(1), 1));
+        row.push_back(std::make_pair(mpq(2731, 1727), 2));
+        row.push_back(std::make_pair(mpq(-910, 1727), 3));
+        g.mk_gomory_cut(t,  k, expl, inf_col, row);
+}
+
+void test_gomory_cut_1() {
+        gomory_test g(
+            [](unsigned j) { return "v" + T_to_string(j);} // name_function_p
+            ,
+            [](unsigned j) { //get_value_p
+                if (j == 1)
+                    return mpq(-2);
+                if (j == 2)
+                    return mpq(4363334, 2730001);
+                if (j == 3)
+                    return mpq(1);
+                lp_assert(false);
+                return zero_of_type<mpq>();
+            },
+            [](unsigned j) { // at_low_p
+                if (j == 1)
+                    return false;
+                if (j == 2)
+                    return false;
+                if (j == 3)
+                    return true;
+                lp_assert(false);
+                return false;
+            },
+            [](unsigned j) { // at_upper
+                if (j == 1)
+                    return true;
+                if (j == 2)
+                    return false;
+                if (j == 3)
+                    return true;
+                lp_assert(false);
+                return false;
+            },
+            [](unsigned j) { // lower_bound
+                if (j == 1) {
+                    lp_assert(false); //unlimited from below
+                    return 0;
+                }
+                if (j == 2)
+                    return 1;
+                if (j == 3)
+                    return 1;
+                lp_assert(false);
+                return 0;
+            },
+            [](unsigned j) { // upper
+                if (j == 1) {
+                    return -2;
+                }
+                if (j == 2)
+                    return 3333;
+                if (j == 3)
+                    return 10000;
+                lp_assert(false);
+                return 0;
+            },
+            [] (unsigned) { return 0; },
+            [] (unsigned) { return 0; }
+                      );
+        lar_term t;
+        mpq k;
+        explanation expl;
+        unsigned inf_col = 2;
+        vector<std::pair<mpq, unsigned>> row;
+        row.push_back(std::make_pair(mpq(1726667, 2730001), 1));
+        row.push_back(std::make_pair(mpq(-910000, 2730001), 3));
+        row.push_back(std::make_pair(mpq(1), 2));
+        g.mk_gomory_cut(t,  k, expl, inf_col, row);
+}
+
+void test_gomory_cut() {
+    test_gomory_cut_0();
+    test_gomory_cut_1();
+}
+
 void test_lp_local(int argn, char**argv) {
     
     // initialize_util_module();
@@ -3270,7 +3420,10 @@ void test_lp_local(int argn, char**argv) {
     }
 
     args_parser.print();
-
+    if (args_parser.option_is_used("-gomory")) {
+        test_gomory_cut();
+        return finalize(0);
+    }
     if (args_parser.option_is_used("-intd")) {
         test_integer_domain();
         return finalize(0);

src/util/lp/active_set.h

@@ -0,0 +1,76 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+
+#pragma once
+#include "util/lp/binary_heap_priority_queue.h"
+namespace lp {
+class active_set {
+    std::unordered_set<constraint*, constraint_hash, constraint_equal> m_cs;
+    binary_heap_priority_queue<int>                                    m_q;
+    std::unordered_map<unsigned, constraint *>                         m_id_to_constraint;
+public:
+    std::unordered_set<constraint*, constraint_hash, constraint_equal> cs() const { return m_cs;}
+
+    bool contains(const constraint* c)  const {
+        return m_id_to_constraint.find(c->id()) != m_id_to_constraint.end();
+    }
+    
+    bool is_empty() const { return m_cs.size() == 0; }
+    // low priority will be dequeued first
+    void add_constraint(constraint* c, int priority) {
+        if (contains(c))
+            return;
+        m_cs.insert(c);
+        m_id_to_constraint[c->id()] = c;
+        m_q.enqueue(c->id(), priority);
+    }
+
+    void clear() {
+        m_cs.clear();
+        m_id_to_constraint.clear();
+        m_q.clear();
+    }
+        
+        
+    constraint* remove_constraint() {
+        if (m_cs.size() == 0)
+            return nullptr;
+        unsigned id = m_q.dequeue();
+        auto it = m_id_to_constraint.find(id);
+        lp_assert(it != m_id_to_constraint.end());
+        constraint* c = it->second;
+        m_cs.erase(c);
+        m_id_to_constraint.erase(it);
+        return c;
+    }
+        
+    unsigned size() const {
+        return static_cast<unsigned>(m_cs.size());
+    }
+
+    void remove_constraint(constraint * c) {
+        if (! contains(c)) return;
+        
+        m_cs.erase(c);
+        m_id_to_constraint.erase(c->id());
+        m_q.remove(c->id());
+    }
+};
+}

src/util/lp/binary_heap_priority_queue.h

@@ -45,7 +45,7 @@ public:
     unsigned size() const { return m_heap_size; }
     binary_heap_priority_queue(): m_heap(1), m_heap_size(0) {} // the empty constructror
     // n is the initial queue capacity.
-    // The capacity will be enlarged two times automatically if needed
+    // The capacity will be enlarged each time twice if needed
     binary_heap_priority_queue(unsigned n);
 
     void clear() {

src/util/lp/binary_heap_priority_queue_def.h

@@ -130,8 +130,12 @@ template <typename T> void binary_heap_priority_queue<T>::enqueue_new(unsigned o
 // In this case the priority will be changed and the queue adjusted.
 template <typename T> void binary_heap_priority_queue<T>::enqueue(unsigned o, const T & priority) {
     if (o >= m_priorities.size()) {
-        resize(o << 1); // make the size twice larger
+        if (o == 0)
+            resize(2);
+        else 
+            resize(o << 1); // make the size twice larger
     }
+    
     if (m_heap_inverse[o] == -1)
         enqueue_new(o, priority);
     else

src/util/lp/bound_analyzer_on_row.h

@@ -19,7 +19,6 @@ Revision History:
 --*/
 #pragma once
 #include "util/vector.h"
-#include "util/lp/linear_combination_iterator.h"
 #include "implied_bound.h"
 #include "test_bound_analyzer.h"
 #include "util/lp/bound_propagator.h"
@@ -28,27 +27,85 @@ Revision History:
 // In a loop we drive the partial sum down, denoting the variables of this process by _u.
 // In the same loop trying to pin variables by pushing the partial sum up, denoting the variable related to it by _l
 namespace lp {
-
+template <typename C> // C plays a role of a container
 class bound_analyzer_on_row {
-    
+    struct term_with_basis_col {
-    linear_combination_iterator<mpq> & m_it;
+        const C &  m_row;
+        unsigned   m_bj;
+        struct ival {
+            unsigned m_var;
+            const mpq & m_coeff;
+            ival(unsigned var, const mpq & val) : m_var(var), m_coeff(val) {
+            }
+            unsigned var() const { return m_var;}
+            const mpq & coeff() const { return m_coeff; }
+        };
+        
+        term_with_basis_col(const C& row, unsigned bj) : m_row(row), m_bj(bj) {}
+        struct const_iterator {
+            // fields
+            typename C::const_iterator m_it;
+            unsigned m_bj;
+            
+
+            //typedefs
+            
+            
+            typedef const_iterator self_type;
+            typedef ival value_type;
+            typedef ival reference;
+            typedef int difference_type;
+            typedef std::forward_iterator_tag iterator_category;
+
+            reference operator*() const {
+                if (m_bj == static_cast<unsigned>(-1))
+                    return ival((*m_it).var(), (*m_it).coeff());
+                return ival(m_bj, - 1);
+            }        
+            self_type operator++() {  self_type i = *this; operator++(1); return i;  }
+            
+            self_type operator++(int) {
+                if (m_bj == static_cast<unsigned>(-1))
+                    m_it++;
+                else
+                    m_bj = static_cast<unsigned>(-1);
+                return *this;
+            }
+
+            // constructor
+            const_iterator(const typename C::const_iterator& it, unsigned bj) :
+                m_it(it),
+                m_bj(bj)
+            {}
+            bool operator==(const self_type &other) const {
+                return m_it == other.m_it && m_bj == other.m_bj ;
+            }
+            bool operator!=(const self_type &other) const { return !(*this == other); }
+        };
+        const_iterator begin() const {
+            return const_iterator( m_row.begin(), m_bj);
+        }
+        const_iterator end() const { return const_iterator(m_row.end(), m_bj); }
+    };
+    term_with_basis_col                m_row;
     bound_propagator &                 m_bp;
-    unsigned           m_row_or_term_index;
+    unsigned                           m_row_or_term_index;
-    int                m_column_of_u; // index of an unlimited from above monoid
+    int                                m_column_of_u; // index of an unlimited from above monoid
     // -1 means that such a value is not found, -2 means that at least two of such monoids were found
-    int                m_column_of_l; // index of an unlimited from below monoid
+    int                                m_column_of_l; // index of an unlimited from below monoid
-    impq               m_rs;
+    impq                               m_rs;
 
 public :
     // constructor
     bound_analyzer_on_row(
-        linear_combination_iterator<mpq> &it,
+        const C & it,
+        unsigned  bj, // basis column for the row
         const numeric_pair<mpq>& rs,
         unsigned row_or_term_index,
         bound_propagator & bp
                           )
         :
-        m_it(it),
+        m_row(it, bj),
         m_bp(bp),
         m_row_or_term_index(row_or_term_index),
         m_column_of_u(-1),
@@ -59,11 +116,11 @@ public :
 
     unsigned j;
     void analyze() {
-        
+        for (auto  c : m_row) {
-        mpq a; unsigned j;
+            if ((m_column_of_l == -2) && (m_column_of_u == -2))
-        while (((m_column_of_l != -2) || (m_column_of_u != -2)) && m_it.next(a, j))
+                break;
-            analyze_bound_on_var_on_coeff(j, a);
+            analyze_bound_on_var_on_coeff(c.var(), c.coeff());
-
+        }
         if (m_column_of_u >= 0)
             limit_monoid_u_from_below();
         else if (m_column_of_u == -1)
@@ -168,25 +225,23 @@ public :
         int strict = 0;
         mpq total;
         lp_assert(is_zero(total));
-        m_it.reset();
+        for (auto p : m_row) {
-        mpq a; unsigned j;
-        while (m_it.next(a, j)) {
             bool str;
-            total -= monoid_min(a, j, str);
+            total -= monoid_min(p.coeff(), p.var(), str);
             if (str)
                 strict++;
         }
 
-        m_it.reset();
+        
-        while (m_it.next(a, j)) {
+        for (auto p : m_row) {
             bool str;
-            bool a_is_pos = is_pos(a);
+            bool a_is_pos = is_pos(p.coeff());
-            mpq bound = total / a + monoid_min_no_mult(a_is_pos, j, str);
+            mpq bound = total / p.coeff() + monoid_min_no_mult(a_is_pos, p.var(), str);
             if (a_is_pos) {
-                limit_j(j, bound, true, false, strict - static_cast<int>(str) > 0);
+                limit_j(p.var(), bound, true, false, strict - static_cast<int>(str) > 0);
             }
             else {
-                limit_j(j, bound, false, true, strict - static_cast<int>(str) > 0);
+                limit_j(p.var(), bound, false, true, strict - static_cast<int>(str) > 0);
             }
         }
     }
@@ -195,25 +250,23 @@ public :
         int strict = 0;
         mpq total;
         lp_assert(is_zero(total));
-        m_it.reset();
+        for (auto p : m_row) {
-        mpq a; unsigned j;
-        while (m_it.next(a, j)) {
             bool str;
-            total -= monoid_max(a, j, str);
+            total -= monoid_max(p.coeff(), p.var(), str);
             if (str)
                 strict++;
         }
-        m_it.reset();
+
-        while (m_it.next(a, j)) {
+        for (auto p : m_row) {
             bool str;
-            bool a_is_pos = is_pos(a);
+            bool a_is_pos = is_pos(p.coeff());
-            mpq bound = total / a + monoid_max_no_mult(a_is_pos, j, str);
+            mpq bound = total / p.coeff() + monoid_max_no_mult(a_is_pos, p.var(), str);
             bool astrict = strict - static_cast<int>(str) > 0; 
             if (a_is_pos) {
-                limit_j(j, bound, true, true, astrict);
+                limit_j(p.var(), bound, true, true, astrict);
             }
             else {
-                limit_j(j, bound, false, false, astrict);
+                limit_j(p.var(), bound, false, false, astrict);
             }
         }
     }
@@ -222,18 +275,18 @@ public :
     void limit_monoid_u_from_below() {
         // we are going to limit from below the monoid m_column_of_u,
         // every other monoid is impossible to limit from below
-        mpq u_coeff, a;
+        mpq u_coeff;
         unsigned j;
         mpq bound = -m_rs.x;
-        m_it.reset();
         bool strict = false;
-        while (m_it.next(a, j)) {
+        for (auto p : m_row) {
+            j = p.var();
             if (j == static_cast<unsigned>(m_column_of_u)) {
-                u_coeff = a;
+                u_coeff = p.coeff();
                 continue;
             }
             bool str;
-            bound -= monoid_max(a, j, str);
+            bound -= monoid_max(p.coeff(), j, str);
             if (str)
                 strict = true;
         }
@@ -251,19 +304,19 @@ public :
     void limit_monoid_l_from_above() {
         // we are going to limit from above the monoid m_column_of_l,
         // every other monoid is impossible to limit from above
-        mpq l_coeff, a;
+        mpq l_coeff;
         unsigned j;
         mpq bound = -m_rs.x;
         bool strict = false;
-        m_it.reset();
+        for (auto p : m_row) {
-        while (m_it.next(a, j)) {
+            j = p.var();
             if (j == static_cast<unsigned>(m_column_of_l)) {
-                l_coeff = a;
+                l_coeff = p.coeff();
                 continue;
             }
 
             bool str;
-            bound -= monoid_min(a, j, str);
+            bound -= monoid_min(p.coeff(), j, str);
             if (str)
                 strict = true;
         }
@@ -282,7 +335,7 @@ public :
     //     bool lower_bound = be.m_lower_bound;
     //     if (!coeff_is_pos)
     //         lower_bound = !lower_bound;
-    //     auto it = m_it.clone();
+    //     auto it = m_row.clone();
     //     mpq a; unsigned j;
     //     while (it->next(a, j)) {
     //         if (be.m_j == j) continue;
@@ -336,14 +389,16 @@ public :
         }
     }
 
-    static void analyze_row(linear_combination_iterator<mpq> &it,
+    static void analyze_row(const C & row,
+                            unsigned bj, // basis column for the row
                             const numeric_pair<mpq>& rs,
                             unsigned row_or_term_index,
                             bound_propagator & bp
                             ) {
-        bound_analyzer_on_row a(it, rs, row_or_term_index, bp);
+        bound_analyzer_on_row a(row, bj, rs, row_or_term_index, bp);
         a.analyze();
     }
 
 };
 }
+

src/util/lp/column_namer.h

@@ -19,31 +19,19 @@ Revision History:
 
 --*/
 #include <string>
-#include "util/lp/linear_combination_iterator.h"
+#include "util/lp/static_matrix.h"
 namespace lp {
 class column_namer {
 public:
     virtual std::string get_column_name(unsigned j) const = 0;
     template <typename T>
-    void print_linear_iterator(linear_combination_iterator<T>* it, std::ostream & out) const {
+    void print_row(const row_strip<T> & row, std::ostream & out) const {
         vector<std::pair<T, unsigned>> coeff;
-        T a;
+        for (auto p : row) {
-        unsigned i;
+            coeff.push_back(std::make_pair(p.coeff(), p.var()));
-        while (it->next(a, i)) {
-            coeff.push_back(std::make_pair(a, i));
         }
         print_linear_combination_of_column_indices(coeff, out);
     }
-    template <typename T>
-    void print_linear_iterator_indices_only(linear_combination_iterator<T>* it, std::ostream & out) const {
-        vector<std::pair<T, unsigned>> coeff;
-        T a;
-        unsigned i;
-        while (it->next(a, i)) {
-            coeff.push_back(std::make_pair(a, i));
-        }
-        print_linear_combination_of_column_indices_only(coeff, out);
-    }
     
     template <typename T>
     void print_linear_combination_of_column_indices_only(const vector<std::pair<T, unsigned>> & coeffs, std::ostream & out) const {

src/util/lp/constraint.h

@@ -0,0 +1,99 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+
+#pragma once
+namespace lp {
+class constraint; // forward definition
+struct constraint_hash {
+    size_t operator() (const constraint* c) const;
+};
+
+struct constraint_equal {
+    bool operator() (const constraint * a, const constraint * b) const;
+};
+
+class constraint { // we only have less or equal for the inequality sign, which is enough for integral variables
+    int                                              m_id;  
+    bool                                             m_is_ineq;
+    polynomial                                       m_poly;
+    mpq                                              m_d; // the divider for the case of a divisibility constraint
+    std::unordered_set<constraint_index>             m_assert_origins; // these indices come from the client and get collected during tightening
+public :
+    unsigned id() const { return m_id; }
+    const polynomial & poly() const { return m_poly; }
+    polynomial & poly() { return m_poly; }
+    std::unordered_set<constraint_index> & assert_origins() { return m_assert_origins;}
+    const std::unordered_set<constraint_index> & assert_origins() const { return m_assert_origins;}
+    bool is_lemma() const { return !is_assert(); }
+    bool is_assert() const { return m_assert_origins.size() == 1; }
+    bool is_ineq() const { return m_is_ineq; }
+    const mpq & divider() const { return m_d; }
+public:
+    constraint(
+        unsigned id,
+        constraint_index assert_origin,
+        const polynomial & p,
+        bool is_ineq):
+        m_id(id),
+        m_is_ineq(is_ineq),
+        m_poly(p)
+    { // creates an assert
+        m_assert_origins.insert(assert_origin);
+    }
+    constraint(
+        unsigned id,
+        const std::unordered_set<constraint_index>& origins,
+        const polynomial & p,
+        bool is_ineq):
+        m_id(id),
+        m_is_ineq(is_ineq),
+        m_poly(p),
+        m_assert_origins(origins)
+    {}
+
+
+    
+    constraint(
+        unsigned id,
+        const polynomial & p,
+        bool is_ineq):
+        m_id(id),
+        m_is_ineq(is_ineq),
+        m_poly(p) { // creates a lemma
+    }
+        
+public:
+    constraint() {}
+
+    const mpq & coeff(var_index j) const {
+        return m_poly.coeff(j);
+    }
+    const vector<monomial>& coeffs() const { return m_poly.m_coeffs;}
+
+    bool is_tight(unsigned j) const {
+        const mpq & a = m_poly.coeff(j);
+        return a == 1 || a == -1;
+    }
+    void add_predecessor(const constraint* p) {
+        lp_assert(p != nullptr);
+        for (auto m : p->assert_origins())
+            m_assert_origins.insert(m); }
+};
+}

src/util/lp/core_solver_pretty_printer_def.h

@@ -105,6 +105,8 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_
                 
             string name = m_core_solver.column_name(column);
             for (unsigned row = 0; row < nrows(); row ++) {
+                m_A[row].resize(ncols(), "");
+                m_signs[row].resize(ncols(),"");
                 set_coeff(
                           m_A[row],
                           m_signs[row],

src/util/lp/cut_solver.cpp

@@ -19,6 +19,5 @@ namespace lp {
         c.s.print_constraint(out, c.c);
         return out;
     }
-
 }
 

src/util/lp/eta_matrix_def.h

@@ -74,7 +74,7 @@ void eta_matrix<T, X>::apply_from_right(vector<T> & w) {
         t += w[it.first] * it.second;
     }
     w[m_column_index] = t;
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(vectors_are_equal<T>(clone_w, w, get_number_of_rows()));
     // delete clone_w;
 #endif
@@ -114,7 +114,7 @@ void eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
         }
     }
     
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(w.is_OK());
     // lp_assert(vectors_are_equal<T>(wcopy, w.m_data));
 #endif
@@ -144,7 +144,7 @@ void eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p) {
     for (auto & pair : m_column_vector.m_data) {
         pair.first = p.get_rev(pair.first);
     }
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(deb == *this);
 #endif
 }

src/util/lp/hash_helper.h

@@ -1,54 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-#include <utility>
-#include <functional>
-#include "util/numerics/mpq.h"
-#ifdef __CLANG__
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wmismatched-tags"
-#endif
-namespace std {
-template<>
-struct hash<lp::mpq> {
-    inline size_t operator()(const lp::mpq & v) const {
-        return v.hash();
-    }
-};
-}
-
-template <class T>
-inline void hash_combine(std::size_t & seed, const T & v) {
-    seed ^= std::hash<T>()(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
-}
-
-namespace std {
-template<typename S, typename T> struct hash<pair<S, T>> {
-    inline size_t operator()(const pair<S, T> & v) const {
-        size_t seed = 0;
-        hash_combine(seed, v.first);
-        hash_combine(seed, v.second);
-        return seed;
-    }
-};
-}
-#ifdef __CLANG__
-#pragma clang diagnostic pop
-#endif

src/util/lp/indexed_vector.h

@@ -166,6 +166,55 @@ public:
             }
         }
     }
+
+    struct ival {
+        unsigned m_var;
+        const T & m_coeff;
+        ival(unsigned var, const T & val) : m_var(var), m_coeff(val) {
+        }
+        unsigned var() const { return m_var;}
+        const T & coeff() const { return m_coeff; }
+    };
+    
+    struct const_iterator {
+            // fields
+        const unsigned *m_i;
+        const indexed_vector& m_v;
+        
+        //typedefs
+            
+            
+        typedef const_iterator self_type;
+        typedef ival value_type;
+        typedef const ival reference;
+        //        typedef const column_cell* pointer;
+        typedef int difference_type;
+        typedef std::forward_iterator_tag iterator_category;
+
+        reference operator*() const {
+            return ival(*m_i, m_v[*m_i]);
+        }        
+        self_type operator++() {  self_type i = *this; m_i++; return i;  }
+        self_type operator++(int) { m_i++; return *this; }
+
+        const_iterator(const unsigned* it, const indexed_vector& v) :
+            m_i(it),
+            m_v(v)
+        {}
+        bool operator==(const self_type &other) const {
+            return m_i == other.m_i;
+        }
+        bool operator!=(const self_type &other) const { return !(*this == other); }
+    };
+
+    const_iterator begin() const {
+        return const_iterator(m_index.begin(), *this);
+    }
+        
+    const_iterator end() const {
+        return const_iterator(m_index.end(), *this);
+    }
+
     
 #ifdef Z3DEBUG
     bool is_OK() const;

src/util/lp/indexer_of_constraints.h

@@ -0,0 +1,45 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+
+#pragma once
+#include "util/lp/binary_heap_priority_queue.h"
+namespace lp {
+
+class indexer_of_constraints {
+    binary_heap_priority_queue<unsigned> m_queue_of_released_indices;
+    unsigned m_max;
+public:
+    indexer_of_constraints() :m_max(0) {}
+    unsigned get_new_index() {
+        unsigned ret;
+        if (m_queue_of_released_indices.is_empty()) {
+            ret = m_max++;
+        }
+        else {
+            ret = m_queue_of_released_indices.dequeue();
+        }
+        return ret;
+    };
+    void release_index(unsigned i) {
+        m_queue_of_released_indices.enqueue(i, i);
+    };
+    unsigned max() const { return m_max; }
+};
+}

src/util/lp/int_solver.cpp

@@ -31,17 +31,17 @@ void int_solver::trace_inf_rows() const {
         }
     }
     
-    num = 0;
+          num = 0;
-    for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) {
+          for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) {
-        unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i];
+              unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i];
-        if (column_is_int_inf(j)) {
+              if (column_is_int_inf(j)) {
-            num++;
+                  num++;
-            iterator_on_row<mpq> it(m_lar_solver->A_r().m_rows[i]);
+                  m_lar_solver->print_row(m_lar_solver->A_r().m_rows[i], tout);
-            m_lar_solver->print_linear_iterator(&it, tout);
+                  tout << "\n";
-            tout << "\n";
+              }
-        }
+          }
-    }
+          tout << "num of int infeasible: " << num << "\n";
-    tout << "num of int infeasible: " << num << "\n";
+          );
 }
 
 int_set& int_solver::inf_int_set() {
@@ -106,32 +106,29 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
     
 }
 
-bool int_solver::is_gomory_cut_target(linear_combination_iterator<mpq> &iter) {
+bool int_solver::is_gomory_cut_target(const row_strip<mpq>& row) {
-    unsigned j;
-    lp_assert(iter.is_reset());
     // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed).
-    while (iter.next(j)) {
+    unsigned j;
+    for (auto p : row) {
+        j = p.var();
         if (is_base(j)) continue;
         if (!is_zero(get_value(j).y)) {
             TRACE("gomory_cut", tout << "row is not gomory cut target:\n";
                   display_column(tout, j);
                   tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";);
-            iter.reset();
             return false;
         }
     }
-    iter.reset();
     return true;
 }
 
 
-void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, unsigned gomory_cut_inf_column) {
+void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, const mpq& f_0, const mpq& one_minus_f_0) {
     TRACE("gomory_cut_detail_real", tout << "real\n";);
-    mpq f_0  = fractional_part(get_value(gomory_cut_inf_column));
     mpq new_a;
     if (at_low(x_j)) {
         if (a.is_pos()) {
-            new_a  =  a / (1 - f_0);
+            new_a  =  a / one_minus_f_0;
         }
         else {
             new_a  =  a / f_0;
@@ -148,7 +145,7 @@ void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, l
             new_a.neg(); // the upper terms are inverted.
         }
         else {
-            new_a =   a / (mpq(1) - f_0); 
+            new_a =   a / one_minus_f_0; 
         }
         k.addmul(new_a, upper_bound(x_j).x); //  k += upper_bound(x_j).x * new_a; 
         expl.push_justification(column_upper_bound_constraint(x_j), new_a);
@@ -166,11 +163,9 @@ constraint_index int_solver::column_lower_bound_constraint(unsigned j) const {
 }
 
 
-void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, unsigned inf_column) {
+void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) {
     lp_assert(is_int(x_j));
     lp_assert(!a.is_int());
-    mpq f_0  = fractional_part(get_value(inf_column));
-    lp_assert(f_0 > zero_of_type<mpq>() && f_0 < one_of_type<mpq>());
     mpq f_j =  fractional_part(a);
     TRACE("gomory_cut_detail", 
           tout << a << " x_j" << x_j << " k = " << k << "\n";
@@ -182,7 +177,6 @@ void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, la
     lp_assert (!f_j.is_zero());
     mpq new_a;
     if (at_low(x_j)) {
-        auto one_minus_f_0 = 1 - f_0;
         if (f_j <= one_minus_f_0) {
             new_a = f_j / one_minus_f_0;
         }
@@ -198,7 +192,7 @@ void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, la
             new_a = f_j / f_0;
         }
         else {
-            new_a = (mpq(1) - f_j) / (1 - f_0);
+            new_a = (mpq(1) - f_j) / one_minus_f_0;
         }
         new_a.neg(); // the upper terms are inverted
         k.addmul(new_a, upper_bound(x_j).x);
@@ -319,18 +313,15 @@ void int_solver::adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k
 
 
 
-lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, linear_combination_iterator<mpq>& iter) {
+lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, const row_strip<mpq> & row) {
 
     lp_assert(column_is_int_inf(inf_col));
 
     TRACE("gomory_cut",
-          tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator_indices_only(&iter, tout); tout << std::endl;
+          tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl;
-          iter.reset();
+          for (auto p : row) {
-          unsigned j;
+              m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout);
-          while(iter.next(j)) {
-              m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, tout);
           }
-          iter.reset();
           tout << "inf_col = " << inf_col << std::endl;
           );
         
@@ -340,18 +331,21 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi
     unsigned x_j;
     mpq a;
     bool some_int_columns = false;
-    lp_assert(iter.is_reset());
+    mpq f_0  = int_solver::fractional_part(get_value(inf_col));
-    while (iter.next(a, x_j)) {
+    mpq one_min_f_0 = 1 - f_0;
+    for (auto p : row) {
+        x_j = p.var();
         if (x_j == inf_col)
             continue;
         // make the format compatible with the format used in: Integrating Simplex with DPLL(T)
+        a = p.coeff();
         a.neg();  
         if (is_real(x_j)) 
-            real_case_in_gomory_cut(a, x_j, k, t, expl, inf_col);
+            real_case_in_gomory_cut(a, x_j, k, t, expl, f_0, one_min_f_0);
         else {
             if (a.is_int()) continue; // f_j will be zero and no monomial will be added
             some_int_columns = true;
-            int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, inf_col);
+            int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, f_0, one_min_f_0);
         }
     }
 
@@ -362,6 +356,7 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi
 
     lp_assert(current_solution_is_inf_on_cut(t, k));
     m_lar_solver->subs_term_columns(t);
+    TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(t, tout); tout << ">= " << k << std::endl;);
     return lia_move::cut;
 }
 
@@ -371,28 +366,30 @@ void int_solver::init_check_data() {
     m_old_values_data.resize(n);
 }
 
-int int_solver::find_free_var_in_gomory_row(linear_combination_iterator<mpq>& iter) {
+int int_solver::find_free_var_in_gomory_row(const row_strip<mpq>& row) {
     unsigned j;
-    while(iter.next(j)) {
+    for (auto p : row) {
+        j = p.var();
         if (!is_base(j) && is_free(j))
             return static_cast<int>(j);
     }
-    iter.reset();
     return -1;
 }
 
-lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j) {
+lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, bool & upper) {
     lia_move ret;
-    linear_combination_iterator<mpq>* iter = m_lar_solver->get_iterator_on_row(row_of_basic_column(j));
+    
-    int free_j = find_free_var_in_gomory_row(*iter);
+    const row_strip<mpq>& row = m_lar_solver->get_row(row_of_basic_column(j));
+    int free_j = find_free_var_in_gomory_row(row);
     if (free_j != -1) {
-        ret = create_branch_on_column(j, t, k, true);
+        ret = create_branch_on_column(j, t, k, true, upper);
-    } else if (!is_gomory_cut_target(*iter)) {
+    } else if (!is_gomory_cut_target(row)) {
-        ret = create_branch_on_column(j, t, k, false);
+        bool upper;
+        ret = create_branch_on_column(j, t, k, false, upper);
     } else {
-        ret = mk_gomory_cut(t, k, ex, j, *iter);
+        upper = false;
+        ret = mk_gomory_cut(t, k, ex, j, row);
     }
-    delete iter;
     return ret;
 }
 
@@ -478,7 +475,14 @@ void int_solver::copy_values_from_cut_solver() {
     }
 }
 
-lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
+void int_solver::catch_up_in_adding_constraints_to_cut_solver() {
+    lp_assert(m_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size());
+    for (unsigned j = m_cut_solver.number_of_asserts(); j < m_lar_solver->constraints().size(); j++) {
+        add_constraint_to_cut_solver(j, m_lar_solver->constraints()[j]);
+    }
+}
+
+lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
     init_check_data();
     lp_assert(inf_int_set_is_correct());
     // it is a reimplementation of 
@@ -486,9 +490,16 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
     // from theory_arith_int.h with the addition of cut_solver
     if (!has_inf_int()) 
         return lia_move::ok;
-    if (settings().m_run_gcd_test)
+    if (settings().m_run_gcd_test) {
-        if (!gcd_test(ex))
+        settings().st().m_gcd_calls++;
+        if (!gcd_test(ex)) {
+            TRACE("gcd_test", tout << "conflict";);
+            settings().st().m_gcd_conflicts++;
             return lia_move::conflict;
+        }
+    } else {
+        TRACE("gcd_test", tout << "no test";);
+    }
     pivoted_rows_tracking_control pc(m_lar_solver);
     /* if (m_params.m_arith_euclidean_solver) apply_euclidean_solver();  */
     //m_lar_solver->pivot_fixed_vars_from_basis();
@@ -498,20 +509,20 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
 
     if ((++m_branch_cut_counter) % settings().m_int_branch_cut_solver == 0) {
         TRACE("check_main_int", tout<<"cut_solver";);
+        catch_up_in_adding_constraints_to_cut_solver();
         auto check_res = m_cut_solver.check();
         settings().st().m_cut_solver_calls++;
         switch (check_res) {
-        case lbool::l_false:
+        case cut_solver::lbool::l_false:
             copy_explanations_from_cut_solver(ex); 
             settings().st().m_cut_solver_false++;
             return lia_move::conflict;
-        case lbool::l_true:
+        case cut_solver::lbool::l_true:
             settings().st().m_cut_solver_true++;
             copy_values_from_cut_solver();
             return lia_move::ok;
-        case lbool::l_undef:
+        case cut_solver::lbool::l_undef:
             settings().st().m_cut_solver_undef++;
-            settings().m_int_branch_cut_solver *= (settings().m_int_branch_cut_solver); // take a square
             break;
         default:
             return lia_move::give_up;
@@ -530,10 +541,11 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
         int j = find_inf_int_base_column();
         if (j == -1) return lia_move::ok;
         TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";);
-        return proceed_with_gomory_cut(t, k, ex, j);
+        
+        return proceed_with_gomory_cut(t, k, ex, j, upper);
     }
     TRACE("check_main_int", tout << "branch"; );
-    return create_branch_on_column(find_inf_int_base_column(), t, k, false);
+    return create_branch_on_column(find_inf_int_base_column(), t, k, false, upper);
 }
 
 bool int_solver::move_non_basic_column_to_bounds(unsigned j) {
@@ -668,26 +680,24 @@ void int_solver::patch_int_infeasible_nbasic_columns() {
     lp_assert(is_feasible() && inf_int_set_is_correct());
 }
 
-mpq get_denominators_lcm(iterator_on_row<mpq> &it) {
+mpq get_denominators_lcm(const row_strip<mpq> & row) {
     mpq r(1);
-    mpq a;
+    for  (auto c : row) {
-    unsigned j;
+        r = lcm(r, denominator(c.coeff()));
-    while (it.next(a, j)) {
-        r = lcm(r, denominator(a));
     }
     return r;
 }
     
 bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, unsigned i, explanation & ex) {
-    iterator_on_row<mpq> it(A.m_rows[i]);
+    mpq lcm_den = get_denominators_lcm(A.m_rows[i]);
-    mpq lcm_den = get_denominators_lcm(it);
     mpq consts(0);
     mpq gcds(0);
     mpq least_coeff(0);
     bool    least_coeff_is_bounded = false;
-    mpq a;
     unsigned j;
-    while (it.next(a, j)) {
+    for (auto &c : A.m_rows[i]) {
+        j = c.var();
+        const mpq& a = c.coeff();
         if (m_lar_solver->column_is_fixed(j)) {
             mpq aux = lcm_den * a;
             consts += aux * m_lar_solver->column_lower_bound(j).x;
@@ -725,8 +735,11 @@ bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, uns
         return true;
     }
         
-    if (!(consts / gcds).is_int())
+    if (!(consts / gcds).is_int()) {
-        fill_explanation_from_fixed_columns(it, ex);
+        TRACE("gcd_test", tout << "row failed the GCD test:\n"; display_row_info(tout, i););
+        fill_explanation_from_fixed_columns(A.m_rows[i], ex);
+        return false;
+    }
         
     if (least_coeff.is_one() && !least_coeff_is_bounded) {
         SASSERT(gcds.is_one());
@@ -734,7 +747,7 @@ bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, uns
     }
         
     if (least_coeff_is_bounded) {
-        return ext_gcd_test(it, least_coeff, lcm_den, consts, ex);
+        return ext_gcd_test(A.m_rows[i], least_coeff, lcm_den, consts, ex);
     }
     return true;
 }
@@ -745,13 +758,11 @@ void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j, expla
     ex.m_explanation.push_back(std::make_pair(mpq(1), lc));
     ex.m_explanation.push_back(std::make_pair(mpq(1), uc));
 }
-void int_solver::fill_explanation_from_fixed_columns(iterator_on_row<mpq> & it, explanation & ex) {
+void int_solver::fill_explanation_from_fixed_columns(const row_strip<mpq> & row, explanation & ex) {
-    it.reset();
+    for (const auto & c : row) {
-    unsigned j;
+        if (!m_lar_solver->column_is_fixed(c.var()))
-    while (it.next(j)) {
-        if (!m_lar_solver->column_is_fixed(j))
             continue;
-        add_to_explanation_from_fixed_or_boxed_column(j, ex);
+        add_to_explanation_from_fixed_or_boxed_column(c.var(), ex);
     }
 }
     
@@ -759,14 +770,13 @@ bool int_solver::gcd_test(explanation & ex) {
     auto & A = m_lar_solver->A_r(); // getting the matrix
     for (unsigned i = 0; i < A.row_count(); i++)
         if (!gcd_test_for_row(A, i, ex)) {
-            std::cout << "false from gcd_test\n" ;
             return false;
         }
         
     return true;
 }
 
-bool int_solver::ext_gcd_test(iterator_on_row<mpq> & it,
+bool int_solver::ext_gcd_test(const row_strip<mpq> & row,
                               mpq const & least_coeff, 
                               mpq const & lcm_den,
                               mpq const & consts, explanation& ex) {
@@ -774,10 +784,11 @@ bool int_solver::ext_gcd_test(iterator_on_row<mpq> & it,
     mpq l(consts);
     mpq u(consts);
 
-    it.reset();
     mpq a;
     unsigned j;
-    while (it.next(a, j)) {
+    for (const auto & c : row) {
+        j = c.var();
+        const mpq & a = c.coeff();
         if (m_lar_solver->column_is_fixed(j))
             continue;
         SASSERT(!m_lar_solver->column_is_real(j));
@@ -817,20 +828,20 @@ bool int_solver::ext_gcd_test(iterator_on_row<mpq> & it,
     mpq u1 = floor(u/gcds);
         
     if (u1 < l1) {
-        fill_explanation_from_fixed_columns(it, ex);
+        fill_explanation_from_fixed_columns(row, ex);
         return false;
     }
         
     return true;
 
 }
-
+/*
 linear_combination_iterator<mpq> * int_solver::get_column_iterator(unsigned j) {
     if (m_lar_solver->use_tableau())
         return new iterator_on_column<mpq, impq>(m_lar_solver->A_r().m_columns[j], m_lar_solver->A_r());
     return new iterator_on_indexed_vector<mpq>(m_lar_solver->get_column_in_lu_mode(j));
 }
-
+*/
 
 int_solver::int_solver(lar_solver* lar_slv) :
     m_lar_solver(lar_slv),
@@ -893,8 +904,7 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
         return false;
 
     impq const & xj = get_value(j);
-    linear_combination_iterator<mpq> *it = get_column_iterator(j);
+    
-
     inf_l = true;
     inf_u = true;
     l = u = zero_of_type<impq>();
@@ -909,7 +919,12 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
 
     mpq a; // the coefficient in the column
     unsigned row_index;
-    while (it->next(a, row_index)) {
+    lp_assert(settings().use_tableau());
+    const auto & A = m_lar_solver->A_r();
+    for (auto c : A.column(j)) {
+        row_index = c.var();
+        const mpq & a = c.coeff();
+        
         unsigned i = lcs.m_r_basis[row_index];
         impq const & xi = get_value(i);
         if (is_int(i) && is_int(j) && !a.is_int())
@@ -930,7 +945,6 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
         if (!inf_l && !inf_u && l == u) break;;                
     }
 
-    delete it;
     TRACE("freedom_interval",
           tout << "freedom variable for:\n";
           tout << m_lar_solver->get_column_name(j);
@@ -1062,21 +1076,16 @@ lp_settings& int_solver::settings() {
 
 void int_solver::display_row_info(std::ostream & out, unsigned row_index) const  {
     auto & rslv = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
-    auto it = m_lar_solver->get_iterator_on_row(row_index);
+    for (auto &c: rslv.m_A.m_rows[row_index]) {
-    mpq a;
+        if (numeric_traits<mpq>::is_pos(c.coeff()))
-    unsigned j;
-    while (it->next(a, j)) {
-        if (numeric_traits<mpq>::is_pos(a))
             out << "+";
-        out << a << rslv.column_name(j) << " ";
+        out << c.coeff() << rslv.column_name(c.var()) << " ";
     }
 
-    it->reset();
+    for (auto& c: rslv.m_A.m_rows[row_index]) {
-    while(it->next(j)) {
+        rslv.print_column_bound_info(c.var(), out);
-        rslv.print_column_bound_info(j, out);
     }
     rslv.print_column_bound_info(rslv.m_basis[row_index], out);
-    delete it;
 }
 
 unsigned int_solver::random() {
@@ -1171,11 +1180,18 @@ const impq& int_solver::lower_bound(unsigned j) const {
     return m_lar_solver->column_lower_bound(j);
 }
 
-lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const {
+lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper) {
     lp_assert(t.is_empty());
     lp_assert(j != -1);
     t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
-    k = free_column? mpq(0) : floor(get_value(j));
+    if (free_column) {
+        upper = true;
+        k = mpq(0);
+    } else {
+        upper = left_branch_is_more_narrow_than_right(j);
+        k = upper? floor(get_value(j)) : ceil(get_value(j));        
+    }
+
     TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n";
           display_column(tout, j);
           tout << "k = " << k << std::endl;
@@ -1184,6 +1200,25 @@ lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool fr
 
 }
 
+bool int_solver::left_branch_is_more_narrow_than_right(unsigned j) {
+    return settings().random_next() % 2;
+    switch (m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_column_types[j] ) {
+    case column_type::fixed:
+        return false;
+    case column_type::boxed:
+        {
+            auto k = floor(get_value(j));
+            return k - lower_bound(j).x < upper_bound(j).x - (k + mpq(1));
+        }
+    case column_type::lower_bound: 
+        return true;
+    case column_type::upper_bound:
+        return false;
+    default:
+        return false;
+    }       
+}
+    
 const impq& int_solver::upper_bound(unsigned j) const {
     return m_lar_solver->column_upper_bound(j);
 }
@@ -1205,19 +1240,14 @@ void int_solver::add_constraint_to_cut_solver(unsigned ci, const lar_base_constr
     vector<mono> coeffs;
     mpq rs;
     get_int_coeffs_from_constraint<mpq>(c, coeffs, rs);
-    svector<constraint_index> explanation;
+    m_cut_solver.add_ineq(coeffs, -rs, ci);
-    explanation.push_back(ci);
-    m_cut_solver.add_ineq(coeffs, -rs, explanation);
-}
-
-void int_solver::notify_on_last_added_constraint() {
-    unsigned ci = m_lar_solver->constraints().size() - 1;
-    const lar_base_constraint* c = m_lar_solver->constraints()[ci];
-    add_constraint_to_cut_solver(ci, c);
 }
 
 void int_solver::pop(unsigned k) {
-    m_cut_solver.pop(k);
+    m_cut_solver.pop_trail(k);
+    while (m_cut_solver.number_of_asserts() > m_lar_solver->constraints().size())
+        m_cut_solver.pop_last_assert();
+    m_cut_solver.pop_constraints();
 }
 
 void int_solver::push() {

src/util/lp/int_solver.h

@@ -1,11 +1,25 @@
-/*
+/*++
-  Copyright (c) 2017 Microsoft Corporation
+Copyright (c) 2017 Microsoft Corporation
-  Author: Lev Nachmanson
+
-*/
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
 #pragma once
 #include "util/lp/lp_settings.h"
 #include "util/lp/static_matrix.h"
-#include "util/lp/iterator_on_row.h"
 #include "util/lp/int_set.h"
 #include "util/lp/lar_term.h"
 #include "util/lp/cut_solver.h"
@@ -49,7 +63,7 @@ public:
     const int_set& inf_int_set() const;
     // main function to check that the solution provided by lar_solver is valid for integral values,
     // or provide a way of how it can be adjusted.
-    lia_move check(lar_term& t, mpq& k, explanation& ex);
+    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);
 private:
@@ -76,17 +90,16 @@ private:
     // creates a fresh inequality.
 
     bool branch(const lp_constraint<mpq, mpq> & new_inequality);
-    bool ext_gcd_test(iterator_on_row<mpq> & it,
+    bool ext_gcd_test(const row_strip<mpq>& row,
                       mpq const & least_coeff, 
                       mpq const & lcm_den,
                       mpq const & consts,
                       explanation & ex);
-    void fill_explanation_from_fixed_columns(iterator_on_row<mpq> & it, explanation &);
+    void fill_explanation_from_fixed_columns(const row_strip<mpq> & row, explanation &);
     void add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation &);
     void patch_int_infeasible_non_basic_column(unsigned j);
     void patch_int_infeasible_nbasic_columns();
     bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m);
-    linear_combination_iterator<mpq> * get_column_iterator(unsigned j);
     const impq & lower_bound(unsigned j) const;
     const impq & upper_bound(unsigned j) const;
     bool is_int(unsigned j) const;
@@ -112,14 +125,13 @@ private:
     lp_settings& settings();
     bool move_non_basic_columns_to_bounds();
     void branch_infeasible_int_var(unsigned);
-    lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, linear_combination_iterator<mpq>& iter);
+    lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, const row_strip<mpq>& row);
     lia_move report_conflict_from_gomory_cut(mpq & k);
     void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq& lcm_den);
     void init_check_data();
-    bool constrain_free_vars(linear_combination_iterator<mpq> *  r);
+    lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, bool & upper);
-    lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j);
+    int find_free_var_in_gomory_row(const row_strip<mpq>& );
-    int find_free_var_in_gomory_row(linear_combination_iterator<mpq>& iter);
+    bool is_gomory_cut_target(const row_strip<mpq>&);
-    bool is_gomory_cut_target(linear_combination_iterator<mpq> &iter);
     bool at_bound(unsigned j) const;
     bool at_low(unsigned j) const;
     bool at_upper(unsigned j) const;
@@ -130,13 +142,15 @@ private:
         return is_zero(n.y);  
     }
 
+public:
     inline static
     mpq fractional_part(const impq & n) {
         lp_assert(is_rational(n));
         return n.x - floor(n.x);
     }
-    void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex, unsigned inf_column);
+private:
-    void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den, unsigned inf_column);
+    void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex, const mpq& f_0, const mpq& one_minus_f_0);
+    void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0);
     constraint_index column_upper_bound_constraint(unsigned j) const;
     constraint_index column_lower_bound_constraint(unsigned j) const;
     void display_row_info(std::ostream & out, unsigned row_index) const;
@@ -147,7 +161,8 @@ public:
 private:
     unsigned random();
     bool has_inf_int() const;
-    lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const;
+    lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper);
+    void catch_up_in_adding_constraints_to_cut_solver();
 public:
     void display_inf_or_int_inf_columns(std::ostream & out) const;
     template <typename T>
@@ -155,11 +170,11 @@ public:
     template <typename T>
     void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector<cut_solver::monomial>& coeff, T & rs);
     bool is_term(unsigned j) const;
-    void notify_on_last_added_constraint();
     void add_constraint_to_cut_solver(unsigned,const lar_base_constraint*);
     void copy_explanations_from_cut_solver(explanation &);
     void pop(unsigned);
     void push();
     void copy_values_from_cut_solver();
+    bool left_branch_is_more_narrow_than_right(unsigned);
 };
 }

src/util/lp/iterator_on_column.h

@@ -1,67 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-#include "util/lp/linear_combination_iterator.h"
-#include "util/lp/static_matrix.h"
-#include "util/lp/lar_term.h"
-namespace lp {
-template <typename T, typename X>
-struct iterator_on_column:linear_combination_iterator<T> {
-    const vector<column_cell>& m_column; // the offset in term coeffs
-    const static_matrix<T, X> & m_A;
-    int m_i; // the initial offset in the column
-    unsigned size() const override { return m_column.size(); }
-    iterator_on_column(const vector<column_cell>& column, const static_matrix<T,X> & A) // the offset in term coeffs
-        :
-        m_column(column),
-        m_A(A),
-        m_i(-1) {}
-    
-    bool next(mpq & a, unsigned & i) override {
-        if (++m_i >= static_cast<int>(m_column.size()))
-            return false;
-
-        const column_cell& c = m_column[m_i];
-        a = m_A.get_val(c);
-        i = c.m_i;
-        return true;
-    }
-
-    bool next(unsigned & i) override {
-        if (++m_i >= static_cast<int>(m_column.size()))
-            return false;
-
-        const column_cell& c = m_column[m_i];
-        i = c.m_i;
-        return true;
-    }
-    
-    void reset() override {
-        m_i = -1;
-    }
-
-    bool is_reset() const { return m_i == -1;}
-
-    linear_combination_iterator<mpq> * clone() {
-        iterator_on_column * r = new iterator_on_column(m_column, m_A);
-        return r;
-    }
-};
-}

src/util/lp/iterator_on_indexed_vector.h

@@ -1,56 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-#include "util/lp/linear_combination_iterator.h"
-namespace lp {
-template <typename T>
-struct iterator_on_indexed_vector:linear_combination_iterator<T> {
-    const indexed_vector<T> & m_v;
-    unsigned m_offset;
-    iterator_on_indexed_vector(const indexed_vector<T> & v) :
-        m_v(v),
-        m_offset(0)
-    {}
-    unsigned size() const override { return m_v.m_index.size(); }
-    bool next(T & a, unsigned & i) override {
-        if (m_offset >= m_v.m_index.size())
-            return false;
-        i = m_v.m_index[m_offset++];
-        a = m_v.m_data[i];
-        return true;
-    }
-    
-    bool next(unsigned & i) override {
-        if (m_offset >= m_v.m_index.size())
-            return false;
-        i = m_v.m_index[m_offset++];
-        return true;
-    }
-    void reset() override {
-        m_offset = 0;
-    }
-
-    bool is_reset() const { return m_offset == 0;}
-
-    linear_combination_iterator<T>* clone() {
-        return new iterator_on_indexed_vector(m_v);
-    }
-};
-}

src/util/lp/iterator_on_pivot_row.h

@@ -1,62 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-#include "util/lp/iterator_on_indexed_vector.h"
-namespace lp {
-template <typename T>
-struct iterator_on_pivot_row:linear_combination_iterator<T> {
-    bool m_basis_returned;
-    const indexed_vector<T> & m_v;
-    unsigned m_basis_j;
-    iterator_on_indexed_vector<T> m_it;
-    unsigned size() const override { return m_it.size(); }
-    iterator_on_pivot_row(const indexed_vector<T> & v, unsigned basis_j) :
-        m_basis_returned(false),
-        m_v(v), m_basis_j(basis_j), m_it(v) {}
-    bool next(T & a, unsigned & i) override {
-        if (m_basis_returned == false) {
-            m_basis_returned = true;
-            a = one_of_type<T>();
-            i = m_basis_j;
-            return true;
-        }
-        return m_it.next(a, i);
-    }
-    bool next(unsigned & i) override {
-        if (m_basis_returned == false) {
-            m_basis_returned = true;
-            i = m_basis_j;
-            return true;
-        }
-        return m_it.next(i);
-    }
-    void reset() override {
-        m_basis_returned = false;
-        m_it.reset();
-    }
-    
-    bool is_reset() const { return m_basis_returned == false && m_it.is_reset();}
-
-    linear_combination_iterator<T> * clone() {
-        iterator_on_pivot_row * r = new iterator_on_pivot_row(m_v, m_basis_j);
-        return r;
-    }
-};
-}

src/util/lp/iterator_on_row.h

@@ -1,55 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-#include "util/lp/linear_combination_iterator.h"
-namespace lp {
-template <typename T>
-struct iterator_on_row:linear_combination_iterator<T> {
-    const vector<row_cell<T>> & m_row;
-    unsigned m_i; // offset
-    iterator_on_row(const vector<row_cell<T>> & row) : m_row(row), m_i(0)
-    {}
-    unsigned size() const override { return m_row.size(); }
-    bool next(T & a, unsigned & i) override {
-        if (m_i == m_row.size())
-            return false;
-        auto &c = m_row[m_i++];
-        i = c.m_j;
-        a = c.get_val();
-        return true;
-    }
-    bool next(unsigned & i) override {
-        if (m_i == m_row.size())
-            return false;
-        auto &c = m_row[m_i++];
-        i = c.m_j;
-        return true;
-    }
-    void reset() override {
-        m_i = 0;
-    }
-
-    bool is_reset() const { return m_i == 0;}
-
-    linear_combination_iterator<T>* clone() {
-        return new iterator_on_row(m_row);
-    }
-};
-}

src/util/lp/iterator_on_term_with_basis_var.h

@@ -1,75 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-#include "util/lp/linear_combination_iterator.h"
-#include "util/lp/numeric_pair.h"
-#include "util/lp/lar_term.h"
-namespace lp {
-struct iterator_on_term_with_basis_var:linear_combination_iterator<mpq> {
-    const lar_term & m_term;
-    std::unordered_map<unsigned, mpq>::const_iterator m_i; // the offset in term coeffs
-    bool             m_term_j_returned;
-    unsigned         m_term_j;
-    unsigned size() const override {return static_cast<unsigned>(m_term.m_coeffs.size() + 1);}
-    iterator_on_term_with_basis_var(const lar_term & t, unsigned term_j) :
-        m_term(t),
-        m_i(t.m_coeffs.begin()),
-        m_term_j_returned(false),
-        m_term_j(term_j) {}
-
-    bool next(mpq & a, unsigned & i) override {
-        if (m_term_j_returned == false) {
-            m_term_j_returned = true;
-            a = - one_of_type<mpq>();
-            i = m_term_j;
-            return true;
-        }
-        if (m_i == m_term.m_coeffs.end())
-            return false;
-        i = m_i->first;
-        a = m_i->second;
-        m_i++;
-        return true;
-    }
-    bool next(unsigned & i) override {
-        if (m_term_j_returned == false) {
-            m_term_j_returned = true;
-            i = m_term_j;
-            return true;
-        }
-        if (m_i == m_term.m_coeffs.end())
-            return false;
-        i = m_i->first;
-        m_i++;
-        return true;
-    }
-
-    bool is_reset() const { return m_term_j_returned == false && m_i == m_term.m_coeffs.begin();}
-
-    void reset() {
-        m_term_j_returned = false;
-        m_i = m_term.m_coeffs.begin();
-    }
-    linear_combination_iterator<mpq> * clone() override {
-        iterator_on_term_with_basis_var * r = new iterator_on_term_with_basis_var(m_term, m_term_j);
-        return r;
-    }
-};
-}

src/util/lp/lar_core_solver.h

@@ -30,8 +30,6 @@ Revision History:
 #include "util/lp/lp_primal_core_solver.h"
 #include "util/lp/stacked_vector.h"
 #include "util/lp/lar_solution_signature.h"
-#include "util/lp/iterator_on_column.h"
-#include "util/lp/iterator_on_indexed_vector.h"
 #include "util/lp/stacked_value.h"
 namespace lp {
 
@@ -801,15 +799,6 @@ public:
         m_r_solver.init_column_row_non_zeroes();
     }
 
-    linear_combination_iterator<mpq> * get_column_iterator(unsigned j) {
-        if (settings().use_tableau()) {
-            return new iterator_on_column<mpq, numeric_pair<mpq>>(m_r_solver.m_A.m_columns[j], m_r_solver.m_A);
-        } else {
-            m_r_solver.solve_Bd(j);
-            return new iterator_on_indexed_vector<mpq>(m_r_solver.m_ed);
-        }
-    }
-
     bool column_is_fixed(unsigned j) const {
         return m_column_types()[j] == column_type::fixed ||
             ( m_column_types()[j] == column_type::boxed &&

src/util/lp/lar_solver.cpp

@@ -151,9 +151,10 @@ void lar_solver::analyze_new_bounds_on_row(
     unsigned row_index,
     bound_propagator & bp) {
     lp_assert(!use_tableau());
-    iterator_on_pivot_row<mpq> it(m_mpq_lar_core_solver.get_pivot_row(), m_mpq_lar_core_solver.m_r_basis[row_index]);
+    unsigned j =  m_mpq_lar_core_solver.m_r_basis[row_index]; // basis column for the row
-     
+    bound_analyzer_on_row<indexed_vector<mpq>>
-    bound_analyzer_on_row ra_pos(it,
+        ra_pos(m_mpq_lar_core_solver.get_pivot_row(),
+                                 j,
                                  zero_of_type<numeric_pair<mpq>>(),
                                  row_index,
                                  bp
@@ -163,14 +164,14 @@ void lar_solver::analyze_new_bounds_on_row(
 
 void lar_solver::analyze_new_bounds_on_row_tableau(
     unsigned row_index,
-    bound_propagator & bp
+    bound_propagator & bp ) {
-                                                   ) {
 
     if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation)
         return;
-    iterator_on_row<mpq> it(A_r().m_rows[row_index]);
+
     lp_assert(use_tableau());
-    bound_analyzer_on_row::analyze_row(it,
+    bound_analyzer_on_row<row_strip<mpq>>::analyze_row(A_r().m_rows[row_index],
+                                       static_cast<unsigned>(-1),
                                        zero_of_type<numeric_pair<mpq>>(),
                                        row_index,
                                        bp
@@ -200,13 +201,6 @@ void lar_solver::calculate_implied_bounds_for_row(unsigned i, bound_propagator &
     }
 }
 
-  
-linear_combination_iterator<mpq> * lar_solver::create_new_iter_from_term(unsigned term_index) const {
-    lp_assert(false); // not implemented
-    return nullptr;
-    //        new linear_combination_iterator_on_vector<mpq>(m_terms[adjust_term_index(term_index)]->coeffs_as_vector());
-}
-
 unsigned lar_solver::adjust_column_index_to_term_index(unsigned j) const {
     unsigned ext_var_or_term = m_columns_to_ext_vars_or_term_indices[j];
     return ext_var_or_term < m_terms_start_index ? j : ext_var_or_term;
@@ -407,6 +401,9 @@ void lar_solver::pop(unsigned k) {
     m_constraints.resize(m_constraint_count);
     m_term_count.pop(k);
     for (unsigned i = m_term_count; i < m_terms.size(); i++) {
+#if Z3DEBUG_CHECK_UNIQUE_TERMS
+        m_set_of_terms.erase(m_terms[i]);
+#endif
         delete m_terms[i];
     }
     m_terms.resize(m_term_count);
@@ -595,7 +592,8 @@ void lar_solver::substitute_terms_in_linear_expression(const vector<std::pair<mp
     }
 
     for (auto & p : coeffs)
-        left_side.push_back(std::make_pair(p.second, p.first));
+        if (!is_zero(p.second))
+            left_side.push_back(std::make_pair(p.second, p.first));
 }
 
 
@@ -1023,7 +1021,7 @@ bool lar_solver::the_right_sides_do_not_sum_to_zero(const vector<std::pair<mpq,
 }
 
 bool lar_solver::explanation_is_correct(const vector<std::pair<mpq, unsigned>>& explanation) const {
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     lconstraint_kind kind;
     lp_assert(the_relations_are_of_same_type(explanation, kind));
     lp_assert(the_left_sides_sum_to_zero(explanation));
@@ -1048,7 +1046,7 @@ bool lar_solver::explanation_is_correct(const vector<std::pair<mpq, unsigned>>&
 }
 
 bool lar_solver::inf_explanation_is_correct() const {
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     vector<std::pair<mpq, unsigned>> explanation;
     get_infeasibility_explanation(explanation);
     return explanation_is_correct(explanation);
@@ -1190,6 +1188,7 @@ void lar_solver::print_constraint(constraint_index ci, std::ostream & out) const
 }
 
 void lar_solver::print_constraints(std::ostream& out) const  {
+    out << "number of constraints = " << m_constraints.size() << std::endl;
     for (auto c : m_constraints) {
         print_constraint(c, out);
     }
@@ -1214,7 +1213,26 @@ void lar_solver::print_term(lar_term const& term, std::ostream & out) const {
     if (!numeric_traits<mpq>::is_zero(term.m_v)) {
         out << term.m_v << " + ";
     }
-    print_linear_combination_of_column_indices(term.coeffs_as_vector(), out);
+    bool first = true;
+    for (const auto p : term) {
+        mpq val = p.coeff();
+        if (first) {
+            first = false;
+        } else {
+            if (is_pos(val)) {
+                out << " + ";
+            } else {
+                out << " - ";
+                val = -val;
+            }
+        }
+        if (val == -numeric_traits<mpq>::one())
+            out << " - ";
+        else if (val != numeric_traits<mpq>::one())
+            out << T_to_string(val);
+        out << this->get_column_name(p.var());
+    }
+
 }
 
 void lar_solver::print_term_as_indices(lar_term const& term, std::ostream & out) const {
@@ -1584,17 +1602,59 @@ void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
 
 var_index lar_solver::add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
                                          const mpq &m_v) {
-    m_terms.push_back(new lar_term(coeffs, m_v));
+    push_and_register_term(new lar_term(coeffs, m_v));
     return m_terms_start_index + m_terms.size() - 1;
 }
 
+#if Z3DEBUG_CHECK_UNIQUE_TERMS
+bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs, const mpq& v) {
+    if (coeffs.empty()) {
+        return is_zero(v);
+    }
+
+    for (const auto & p : coeffs) {
+        if (column_is_real(p.second))
+            return true;
+    }
+    
+    mpq g;
+    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) {
+            g_is_set = true;
+            g = p.first;
+        } else {
+            g = gcd(g, p.first);
+        }
+    }
+    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
+void lar_solver::push_and_register_term(lar_term* t) {
+#if Z3DEBUG_CHECK_UNIQUE_TERMS
+    lp_assert(m_set_of_terms.find(t) == m_set_of_terms.end());
+    m_set_of_terms.insert(t);
+#endif
+    m_terms.push_back(t);
+}
+
 // terms
 var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
                                const mpq &m_v) {
     if (strategy_is_undecided())
         return add_term_undecided(coeffs, m_v);
 
-    m_terms.push_back(new lar_term(coeffs, m_v));
+    push_and_register_term(new lar_term(coeffs, m_v));
     unsigned adjusted_term_index = m_terms.size() - 1;
     var_index ret = m_terms_start_index + adjusted_term_index;
     if (use_tableau() && !coeffs.empty()) {
@@ -1607,6 +1667,7 @@ var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
 }
 
 void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
+    TRACE("dump_terms", print_term(*term, tout); tout << std::endl;);
     register_new_ext_var_index(term_ext_index, term_is_int(term));
     // j will be a new variable
     unsigned j = A_r().column_count();
@@ -1614,8 +1675,8 @@ void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned
     m_columns_to_ul_pairs.push_back(ul);
     add_basic_var_to_core_fields();
     if (use_tableau()) {
-        auto it = iterator_on_term_with_basis_var(*term, j);
+        A_r().fill_last_row_with_pivoting(*term,
-        A_r().fill_last_row_with_pivoting(it,
+                                          j,
                                           m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
         m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
     }
@@ -1653,9 +1714,6 @@ constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, c
         lp_assert(bound_is_integer_for_integer_column(j, right_side));
         auto vc = new lar_var_constraint(j, kind, right_side);
         m_constraints.push_back(vc);
-        if (has_int_var()) {
-            m_int_solver->notify_on_last_added_constraint();
-        }
         update_column_type_and_bound(j, kind, right_side, ci);
     }
     else {
@@ -1697,8 +1755,6 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k
         unsigned term_j = it->second.internal_j();
         mpq rs = right_side - m_terms[adjusted_term_index]->m_v;
         m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side));
-        if (has_int_var())
-            m_int_solver->notify_on_last_added_constraint();
         update_column_type_and_bound(term_j, kind, rs, ci);
     }
     else {

src/util/lp/lar_solver.h

@@ -37,15 +37,13 @@ Revision History:
 #include "util/lp/stacked_value.h"
 #include "util/lp/stacked_vector.h"
 #include "util/lp/stacked_unordered_set.h"
-#include "util/lp/iterator_on_pivot_row.h"
 #include "util/lp/implied_bound.h"
 #include "util/lp/bound_analyzer_on_row.h"
-#include "util/lp/iterator_on_term_with_basis_var.h"
-#include "util/lp/iterator_on_row.h"
 #include "util/lp/quick_xplain.h"
 #include "util/lp/conversion_helper.h"
 #include "util/lp/int_solver.h"
 #include "util/lp/nra_solver.h"
+#include "util/lp/bound_propagator.h"
 
 namespace lp {
 
@@ -60,30 +58,71 @@ class lar_solver : public column_namer {
         unsigned internal_j() const { return m_internal_j;}
         bool is_integer() const {return m_is_integer;}
     };
+#if Z3DEBUG_CHECK_UNIQUE_TERMS
+    struct term_hasher {
+        std::size_t operator()(const lar_term *t) const
+        {
+            using std::size_t;
+            using std::hash;
+            using std::string;
+            size_t seed = 0;
+            for (const auto& p : t->m_coeffs) {
+                hash_combine(seed, p);
+            }
+            return seed;
+        }
+    };
+
+    struct term_ls_comparer {
+        bool operator()(const lar_term *a, const lar_term* b) const
+        {
+            // a is contained in b
+            for (auto & p : a->m_coeffs) {
+                auto t = b->m_coeffs.find(p.first);
+                if (t == b->m_coeffs.end())
+                    return false;
+                if (p.second != t->second)
+                    return false;
+            }
+            // zz is contained in b
+            for (auto & p : b->m_coeffs) {
+                auto t = a->m_coeffs.find(p.first);
+                if (t == a->m_coeffs.end())
+                    return false;
+                if (p.second != t->second)
+                    return false;
+            }
+            return true;
+        }
+    };
+    std::unordered_set<lar_term*, term_hasher, term_ls_comparer>    m_set_of_terms;
+#endif
+    
     //////////////////// fields //////////////////////////
-    lp_settings m_settings;
+    lp_settings                                         m_settings;
-    lp_status m_status;
+    lp_status                                           m_status;
-    stacked_value<simplex_strategy_enum> m_simplex_strategy;
+    stacked_value<simplex_strategy_enum>                m_simplex_strategy;
-    std::unordered_map<unsigned, ext_var_info> m_ext_vars_to_columns;
+    std::unordered_map<unsigned, ext_var_info>          m_ext_vars_to_columns;
-    vector<unsigned> m_columns_to_ext_vars_or_term_indices;
+    vector<unsigned>                                    m_columns_to_ext_vars_or_term_indices;
-    stacked_vector<ul_pair> m_columns_to_ul_pairs;
+    stacked_vector<ul_pair>                             m_columns_to_ul_pairs;
-    vector<lar_base_constraint*> m_constraints;
+    vector<lar_base_constraint*>                        m_constraints;
+private:
+    stacked_value<unsigned>                             m_constraint_count;
+    // the set of column indices j such that bounds have changed for j
+    int_set                                             m_columns_with_changed_bound;
+    int_set                                             m_rows_with_changed_bounds;
+    int_set                                             m_basic_columns_with_changed_cost;
+    stacked_value<int>                                  m_infeasible_column_index; // such can be found at the initialization step
+    stacked_value<unsigned>                             m_term_count;
+    vector<lar_term*>                                   m_terms;
+    const var_index                                     m_terms_start_index;
+    indexed_vector<mpq>                                 m_column_buffer;
+
+    
 public :
     const vector<lar_base_constraint*>& constraints() const {
         return m_constraints;
     }
-private:
-    stacked_value<unsigned> m_constraint_count;
-    // the set of column indices j such that bounds have changed for j
-    int_set m_columns_with_changed_bound;
-    int_set m_rows_with_changed_bounds;
-    int_set m_basic_columns_with_changed_cost;
-    stacked_value<int> m_infeasible_column_index; // such can be found at the initialization step
-    stacked_value<unsigned> m_term_count;
-    vector<lar_term*> m_terms;
-    const var_index m_terms_start_index;
-    indexed_vector<mpq> m_column_buffer;
-public:
     lar_core_solver m_mpq_lar_core_solver;
 private:
     std::function<void (unsigned)> m_tracker_of_x_change;
@@ -138,13 +177,17 @@ public:
     void add_new_var_to_core_fields_for_mpq(bool register_in_basis);
 
 
-    var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
-                                 const mpq &m_v);
 
     // terms
     var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
                        const mpq &m_v);
 
+    var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
+                                 const mpq &m_v);
+
+    bool term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs, const mpq& v);
+    void push_and_register_term(lar_term* t);
+
     void add_row_for_term(const lar_term * term, unsigned term_ext_index);
 
     void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index);
@@ -218,9 +261,6 @@ public:
     
     void calculate_implied_bounds_for_row(unsigned i, bound_propagator & bp);
 
-  
-    linear_combination_iterator<mpq> * create_new_iter_from_term(unsigned term_index) const;
-
     unsigned adjust_column_index_to_term_index(unsigned j) const;
     
     void propagate_bounds_on_a_term(const lar_term& t, bound_propagator & bp, unsigned term_offset);
@@ -301,8 +341,6 @@ public:
 
     void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j);
 
-
-    
     void detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j);
 
     bool use_tableau() const;
@@ -479,10 +517,12 @@ public:
     }
     
     bool bound_is_integer_for_integer_column(unsigned j, const mpq & right_side) const;
-    linear_combination_iterator<mpq> * get_iterator_on_row(unsigned i) {
+    
-        return m_mpq_lar_core_solver.m_r_solver.get_iterator_on_row(i);
+    const row_strip<mpq> &  get_row(unsigned i) {
+        return A_r().m_rows[i];
     }
 
+    
     unsigned get_base_column_in_row(unsigned row_index) const {
         return m_mpq_lar_core_solver.m_r_solver.get_base_column_in_row(row_index);
     }
@@ -495,17 +535,20 @@ public:
         return m_columns_to_ul_pairs()[j].lower_bound_witness();
     }
 
-    void subs_terms_for_debugging(lar_term& t) {
+    void subs_term_columns(lar_term& t) {
-        vector<std::pair<mpq, unsigned>> pol;
+        vector<std::pair<unsigned,unsigned>> columns_to_subs;
         for (const auto & m : t.m_coeffs) {
-            pol.push_back(std::make_pair(m.second, adjust_column_index_to_term_index(m.first)));
+            unsigned tj = adjust_column_index_to_term_index(m.first);
+            if (tj == m.first) continue;
+            columns_to_subs.push_back(std::make_pair(m.first, tj));
+        }
+        for (const auto & p : columns_to_subs) {
+            auto it = t.m_coeffs.find(p.first);
+            lp_assert(it != t.m_coeffs.end());
+            mpq v = it->second;
+            t.m_coeffs.erase(it);
+            t.m_coeffs[p.second] = v;
         }
-        mpq v = t.m_v;
-        
-        vector<std::pair<mpq, unsigned>> pol_after_subs;
-        substitute_terms_in_linear_expression(pol, pol_after_subs, v);
-        t.clear();
-        t = lar_term(pol_after_subs, v);
     }
 
     bool inf_int_set_is_correct_for_column(unsigned j) const {

src/util/lp/lar_term.h

@@ -1,22 +1,22 @@
 /*++
-Copyright (c) 2017 Microsoft Corporation
+  Copyright (c) 2017 Microsoft Corporation
 
-Module Name:
+  Module Name:
 
-    <name>
+  <name>
 
-Abstract:
+  Abstract:
 
-    <abstract>
+  <abstract>
 
-Author:
+  Author:
 
-    Lev Nachmanson (levnach)
+  Lev Nachmanson (levnach)
 
-Revision History:
+  Revision History:
 
 
---*/
+  --*/
 #pragma once
 #include "util/lp/indexed_vector.h"
 namespace lp {
@@ -85,11 +85,55 @@ struct lar_term {
             t.second.neg();
     }
 
+    template <typename T>
+    T apply(const vector<T>& x) const {
+        T ret = T(m_v);
+        for (const auto & t : m_coeffs) {
+            ret += t.second * x[t.first];
+        }
+        return ret;
+    }
+   
     void clear() {
         m_coeffs.clear();
         m_v = zero_of_type<mpq>();
     }
 
+    struct ival {
+        unsigned m_var;
+        const mpq & m_coeff;
+        ival(unsigned var, const mpq & val) : m_var(var), m_coeff(val) {
+        }
+        unsigned var() const { return m_var;}
+        const mpq & coeff() const { return m_coeff; }
+    };
     
+    struct const_iterator {
+        //fields
+        std::unordered_map<unsigned, mpq>::const_iterator m_it;
+
+        typedef const_iterator self_type;
+        typedef ival value_type;
+        typedef ival reference;
+        //        typedef std::pair<const unsigned, mpq>* pointer;
+        typedef int difference_type;
+        typedef std::forward_iterator_tag iterator_category;
+
+        reference operator*() const {
+            return ival(m_it->first, m_it->second);
+        }
+        
+        self_type operator++() {  self_type i = *this; m_it++; return i;  }
+        self_type operator++(int) { m_it++; return *this; }
+
+        const_iterator(std::unordered_map<unsigned, mpq>::const_iterator it) : m_it(it) {}
+        bool operator==(const self_type &other) const {
+            return m_it == other.m_it;
+        }
+        bool operator!=(const self_type &other) const { return !(*this == other); }
+    };
+
+    const_iterator begin() const { return m_coeffs.begin();}
+    const_iterator end() const { return m_coeffs.end(); }
 };
 }

src/util/lp/linear_combination_iterator.h

@@ -1,66 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-namespace lp {
-template <typename T>
-struct linear_combination_iterator {
-    virtual bool next(T & a, unsigned & i) = 0;
-    virtual bool next(unsigned & i) = 0;
-    virtual void reset() = 0;
-    virtual linear_combination_iterator * clone() = 0;
-    virtual ~linear_combination_iterator(){}
-    virtual unsigned size() const = 0;
-    virtual bool is_reset() const = 0;
-};
-template <typename T>
-struct linear_combination_iterator_on_vector : linear_combination_iterator<T> {
-    vector<std::pair<T, unsigned>> & m_vector;
-    int m_offset;
-    bool next(T & a, unsigned & i) {
-        if(static_cast<unsigned>(m_offset) >= m_vector.size())
-            return false;
-        auto & p = m_vector[m_offset];
-        a = p.first;
-        i = p.second;
-        m_offset++;
-        return true;
-    }
-
-    bool next(unsigned & i) {
-        if(static_cast<unsigned>(m_offset) >= m_vector.size())
-            return false;
-        auto & p = m_vector[m_offset];
-        i = p.second;
-        m_offset++;
-        return true;
-    }
-    bool is_reset() const { return m_offset == 0;}
-    void reset() {m_offset = 0;}
-    linear_combination_iterator<T> * clone() {
-        return new linear_combination_iterator_on_vector(m_vector);
-    }
-    linear_combination_iterator_on_vector(vector<std::pair<T, unsigned>> & vec):
-        m_vector(vec),
-        m_offset(0)
-    {}
-    unsigned size() const { return m_vector.size(); }
-};
-
-}

src/util/lp/lp_bound_propagator.cpp

@@ -1,67 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#include "util/lp/lar_solver.h"
-namespace lp {
-    lp_bound_propagator::lp_bound_propagator(lar_solver & ls):
-    m_lar_solver(ls) {}
-column_type lp_bound_propagator::get_column_type(unsigned j) const {
-    return m_lar_solver.m_mpq_lar_core_solver.m_column_types()[j];
-}
-const impq & lp_bound_propagator::get_low_bound(unsigned j) const {
-    return m_lar_solver.m_mpq_lar_core_solver.m_r_low_bounds()[j];
-}
-const impq & lp_bound_propagator::get_upper_bound(unsigned j) const {
-    return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j];
-}
-void lp_bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
-    unsigned term_j = m_lar_solver.adjust_column_index_to_term_index(j);
-    mpq w = v;
-    if (term_j != j) {
-        j = term_j;
-        w += m_lar_solver.get_term(term_j).m_v; // when terms are turned into the columns they "lose" the right side, at this moment they aquire it back
-    }
-    lconstraint_kind kind = is_low? GE : LE;
-    if (strict)
-        kind = static_cast<lconstraint_kind>(kind / 2);
-
-    if (!bound_is_interesting(j, kind, w))
-        return;
-     unsigned k; // index to ibounds
-     if (is_low) {
-         if (try_get_val(m_improved_low_bounds, j, k)) {
-             auto & found_bound = m_ibounds[k];
-             if (w > found_bound.m_bound || (w == found_bound.m_bound && found_bound.m_strict == false && strict))
-                 found_bound = implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
-         } else {
-             m_improved_low_bounds[j] = m_ibounds.size();
-             m_ibounds.push_back(implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
-         }
-     } else { // the upper bound case
-         if (try_get_val(m_improved_upper_bounds, j, k)) {
-             auto & found_bound = m_ibounds[k];
-             if (w < found_bound.m_bound || (w == found_bound.m_bound && found_bound.m_strict == false && strict))
-                 found_bound = implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
-         } else {
-             m_improved_upper_bounds[j] = m_ibounds.size();
-             m_ibounds.push_back(implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
-        }
-     }
-}
-}

src/util/lp/lp_bound_propagator.h

@@ -1,42 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-#include "util/lp/lp_settings.h"
-namespace lp {
-class lar_solver;
-class bound_propagator {
-    std::unordered_map<unsigned, unsigned> m_improved_lower_bounds; // these maps map a column index to the corresponding index in ibounds
-    std::unordered_map<unsigned, unsigned> m_improved_upper_bounds;
-    lar_solver & m_lar_solver;
-public:
-    vector<implied_bound> m_ibounds;
-public:
-    lp_bound_propagator(lar_solver & ls);
-    column_type get_column_type(unsigned) const;
-    const impq & get_lower_bound(unsigned) const;
-    const impq & get_upper_bound(unsigned) const;
-    void try_add_bound(mpq  v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict);
-    virtual bool bound_is_interesting(unsigned vi,
-                                      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"; }
-};
-}

src/util/lp/lp_core_solver_base.h

@@ -28,8 +28,6 @@ Revision History:
 #include "util/lp/lu.h"
 #include "util/lp/permutation_matrix.h"
 #include "util/lp/column_namer.h"
-#include "util/lp/iterator_on_row.h"
-#include "util/lp/iterator_on_pivot_row.h"
 
 namespace lp {
 
@@ -751,13 +749,6 @@ public:
         return m_iters_with_no_cost_growing;
     }
 
-    linear_combination_iterator<T> * get_iterator_on_row(unsigned i) {
-        if (m_settings.use_tableau())
-            return new iterator_on_row<T>(m_A.m_rows[i]);
-        calculate_pivot_row(i);
-        return new iterator_on_pivot_row<T>(m_pivot_row, m_basis[i]);
-    }
-
     void calculate_pivot_row(unsigned i);
     unsigned get_base_column_in_row(unsigned row_index) const {
         return m_basis[row_index];

src/util/lp/lp_core_solver_base_def.h

@@ -159,7 +159,7 @@ solve_Bd(unsigned entering) {
         m_columns_nz[entering] = m_ed.m_index.size();
     lp_assert(m_ed.is_OK());
     lp_assert(m_w.is_OK());
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // auto B = get_B(*m_factorization, m_basis);
     // vector<T>  a(m_m());
     // m_A.copy_column_to_vector(entering, a);
@@ -978,8 +978,8 @@ template <typename T, typename X>  void lp_core_solver_base<T, X>::pivot_fixed_v
         if (get_column_type(basic_j) != column_type::fixed) continue;
         T a;
         unsigned j;
-        auto * it = get_iterator_on_row(i);
+        for (auto &c : m_A.m_rows[i]) {
-        while (it->next(a, j)) {
+            j = c.var();
             if (j == basic_j)
                 continue;
             if (get_column_type(j) != column_type::fixed) {
@@ -987,7 +987,6 @@ template <typename T, typename X>  void lp_core_solver_base<T, X>::pivot_fixed_v
                     break;
             }
         }
-        delete it;
     }
 }
 

src/util/lp/lp_primal_core_solver.h

@@ -37,7 +37,6 @@ Revision History:
 #include "util/lp/breakpoint.h"
 #include "util/lp/binary_heap_priority_queue.h"
 #include "util/lp/int_set.h"
-#include "util/lp/iterator_on_row.h"
 namespace lp {
 
 // This core solver solves (Ax=b, lower_bound_values \leq x \leq upper_bound_values, maximize costs*x )

src/util/lp/lp_primal_simplex.h

@@ -26,7 +26,6 @@ Revision History:
 #include "util/lp/column_info.h"
 #include "util/lp/lp_primal_core_solver.h"
 #include "util/lp/lp_solver.h"
-#include "util/lp/iterator_on_row.h"
 namespace lp {
 template <typename T, typename X>
 class lp_primal_simplex: public lp_solver<T, X> {

src/util/lp/lp_settings.h

@@ -106,6 +106,8 @@ struct stats {
     unsigned m_cut_solver_true;
     unsigned m_cut_solver_false;
     unsigned m_cut_solver_undef;
+    unsigned m_gcd_calls;
+    unsigned m_gcd_conflicts;
     stats() { reset(); }
     void reset() { memset(this, 0, sizeof(*this)); }
 };
@@ -226,9 +228,8 @@ public:
                     backup_costs(true),
                     column_number_threshold_for_using_lu_in_lar_solver(4000),
                     m_int_branch_cut_gomory_threshold(4),
-                    m_int_branch_cut_solver(4),
+                    m_int_branch_cut_solver(8),
                     m_run_gcd_test(true),
-                    m_cut_solver_bound_propagation_factor(5),
                     m_cut_solver_cycle_on_var(10)
     {}
 
@@ -339,7 +340,6 @@ public:
     unsigned m_int_branch_cut_gomory_threshold;
     unsigned m_int_branch_cut_solver;
     bool m_run_gcd_test;
-    unsigned m_cut_solver_bound_propagation_factor;
     unsigned m_cut_solver_cycle_on_var;
 }; // end of lp_settings class
 

src/util/lp/lp_solver.h

@@ -28,7 +28,6 @@ Revision History:
 #include "util/lp/static_matrix.h"
 #include "util/lp/lp_core_solver_base.h"
 #include "util/lp/scaler.h"
-#include "util/lp/linear_combination_iterator.h"
 #include "util/lp/bound_analyzer_on_row.h"
 namespace lp {
 enum lp_relation  {

src/util/lp/lp_utils.h

@@ -38,7 +38,7 @@ bool contains(const std::unordered_map<A, B> & map, const A& key) {
 #ifdef lp_for_z3
 
 #ifdef Z3DEBUG
-#define LEAN_DEBUG 1
+#define Z3DEBUG 1
 #endif
 
 namespace lp {

src/util/lp/lu.cpp

@@ -39,7 +39,7 @@ template lp::mpq lp::dot_product<lp::mpq, lp::mpq>(vector<lp::mpq > const&, vect
 template void lp::init_factorization<double, double>(lp::lu<double, double>*&, lp::static_matrix<double, double>&, vector<unsigned int>&, lp::lp_settings&);
 template void lp::init_factorization<lp::mpq, lp::mpq>(lp::lu<lp::mpq, lp::mpq>*&, lp::static_matrix<lp::mpq, lp::mpq>&, vector<unsigned int>&, lp::lp_settings&);
 template void lp::init_factorization<lp::mpq, lp::numeric_pair<lp::mpq> >(lp::lu<lp::mpq, lp::numeric_pair<lp::mpq> >*&, lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&, vector<unsigned int>&, lp::lp_settings&);
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
 template void lp::print_matrix<double, double>(lp::sparse_matrix<double, double>&, std::ostream & out);
 template void lp::print_matrix<lp::mpq, lp::mpq>(lp::static_matrix<lp::mpq, lp::mpq>&, std::ostream&);
 template void lp::print_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >(lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&, std::ostream&);

src/util/lp/lu.h

@@ -34,7 +34,7 @@ Revision History:
 #include "util/lp/square_dense_submatrix.h"
 #include "util/lp/dense_matrix.h"
 namespace lp {
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
 template <typename T, typename X> // print the nr x nc submatrix at the top left corner
 void print_submatrix(sparse_matrix<T, X> & m, unsigned mr, unsigned nc);
 
@@ -113,7 +113,7 @@ public:
 #endif
         m_i = p.apply_reverse(m_i);
 
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
         // lp_assert(*this == deb);
 #endif
     }

src/util/lp/lu_def.h

@@ -25,7 +25,7 @@ Revision History:
 #include "util/debug.h"
 #include "util/lp/lu.h"
 namespace lp {
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
 template <typename T, typename X> // print the nr x nc submatrix at the top left corner
 void print_submatrix(sparse_matrix<T, X> & m, unsigned mr, unsigned nc, std::ostream & out) {
     vector<vector<std::string>> A;
@@ -138,12 +138,12 @@ lu<T, X>::lu(static_matrix<T, X> const & A,
     m_row_eta_work_vector(A.row_count()),
     m_refactor_counter(0) {
     lp_assert(!(numeric_traits<T>::precise() && settings.use_tableau()));
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     debug_test_of_basis(A, basis);
 #endif
     ++m_settings.st().m_num_factorizations;
     create_initial_factorization();
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(check_correctness());
 #endif
 }

src/util/lp/monomial.h

@@ -0,0 +1,33 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+namespace lp {
+struct monomial {
+    mpq           m_coeff; // the coefficient of the monomial
+    var_index     m_var; // the variable index
+public:
+    monomial(const mpq& coeff, var_index var) : m_coeff(coeff), m_var(var) {}
+    monomial(var_index var) : monomial(one_of_type<mpq>(), var) {}
+    const mpq & coeff() const { return m_coeff; }
+    mpq & coeff() { return m_coeff; }
+    var_index var() const { return m_var; }
+    std::pair<mpq, var_index> to_pair() const { return std::make_pair(coeff(), var());}
+};
+}

src/util/lp/permutation_matrix.h

@@ -28,7 +28,7 @@ Revision History:
 #include "util/lp/matrix.h"
 #include "util/lp/tail_matrix.h"
 namespace lp {
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     inline bool is_even(int k) {  return (k/2)*2 == k; }
 #endif
 

src/util/lp/permutation_matrix_def.h

@@ -74,7 +74,7 @@ void permutation_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings & ) {
     while (i-- > 0) {
         w[i] = m_X_buffer[i];
     }
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(vectors_are_equal<L>(deb_w, w, row_count()));
     // delete [] deb_w;
 #endif
@@ -109,7 +109,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::apply_from_righ
     for (unsigned i = 0; i < size(); i++) {
         w[i] = m_T_buffer[i];
     }
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(vectors_are_equal<T>(deb_w, w, row_count()));
     // delete [] deb_w;
 #endif
@@ -133,7 +133,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::apply_from_righ
         w.set_value(buffer[i], pj);
     }
     lp_assert(w.is_OK());
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     lp_assert(vectors_are_equal(wcopy, w.m_data));
 #endif
 }
@@ -180,7 +180,7 @@ void permutation_matrix<T, X>::apply_reverse_from_left(indexed_vector<L> & w) {
         w[j] = t[i];
         w.m_index[i] = j;
     }
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(vectors_are_equal<L>(deb_w, w.m_data, row_count()));
     // delete [] deb_w;
 #endif

src/util/lp/polynomial.h

@@ -0,0 +1,106 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+namespace lp {
+struct polynomial {
+    static mpq m_local_zero;
+    // the polynomial evaluates to m_coeffs + m_a
+    vector<monomial>        m_coeffs;
+    mpq                     m_a; // the free coefficient
+    polynomial(const vector<monomial>& p, const mpq & a) : m_coeffs(p), m_a(a) {}
+    polynomial(const vector<monomial>& p) : polynomial(p, zero_of_type<mpq>()) {}
+    polynomial(): m_a(zero_of_type<mpq>()) {}
+    polynomial(const polynomial & p) : m_coeffs(p.m_coeffs), m_a(p.m_a) {} 
+            
+    const mpq & coeff(var_index j) const {
+        for (const auto & t : m_coeffs) {
+            if (j == t.var()) {
+                return t.coeff();
+            }
+        }
+        return m_local_zero;
+    }
+
+    polynomial &  operator+=(const polynomial & p) {
+        m_a += p.m_a;
+        for (const auto & t: p.m_coeffs)
+            *this += monomial(t.coeff(), t.var());
+        return *this;
+    }
+
+    void add(const mpq & c, const polynomial &p) {
+        m_a += p.m_a * c;
+            
+        for (const auto & t: p.m_coeffs)
+            *this += monomial(c * t.coeff(), t.var());
+    }
+        
+    void clear() {
+        m_coeffs.clear();
+        m_a = zero_of_type<mpq>();
+    }
+        
+    bool is_empty() const { return m_coeffs.size() == 0 && numeric_traits<mpq>::is_zero(m_a); }
+
+    unsigned number_of_monomials() const { return m_coeffs.size();}
+
+    void add(const monomial &m ){
+        if (is_zero(m.coeff())) return;
+        for (unsigned k = 0; k < m_coeffs.size(); k++) {
+            auto & l = m_coeffs[k];
+            if (m.var() == l.var()) {
+                l.coeff() += m.coeff();
+                if (l.coeff() == 0)
+                    m_coeffs.erase(m_coeffs.begin() + k);
+                return;
+            }
+        }
+        m_coeffs.push_back(m);
+        lp_assert(is_correct());
+    }
+        
+    polynomial & operator+=(const monomial &m ){
+        add(m);
+        return *this;
+    }
+
+    polynomial & operator+=(const mpq &c ){
+        m_a += c;
+        return *this;
+    }
+
+        
+    bool is_correct() const {
+        std::unordered_set<var_index> s;
+        for (auto & l : m_coeffs) {
+            if (l.coeff() == 0)
+                return false;
+            s.insert(l.var());
+        }
+        return m_coeffs.size() == s.size();
+    }
+
+    bool var_coeff_is_unit(unsigned j) const {
+        const mpq & a = coeff(j);
+        return a == 1 || a == -1;
+    }
+    const vector<monomial> & coeffs() const { return m_coeffs; }
+};
+}

src/util/lp/random_updater.h

@@ -25,7 +25,6 @@ Revision History:
 #include <string>
 #include <algorithm>
 #include "util/lp/lp_settings.h"
-#include "util/lp/linear_combination_iterator.h"
 // see http://research.microsoft.com/projects/z3/smt07.pdf
 // The class searches for a feasible solution with as many different values of variables as it can find
 namespace lp {

src/util/lp/random_updater_def.h

@@ -83,7 +83,7 @@ void random_updater::add_column_to_sets(unsigned j) {
         add_value(m_lar_solver.get_core_solver().m_r_x[j]);
     } else {
         unsigned row = m_lar_solver.get_core_solver().m_r_heading[j];
-        for (auto row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row]) {
+        for (auto & row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row]) {
             unsigned cj = row_c.m_j;
             if (m_lar_solver.get_core_solver().m_r_heading[cj] < 0) {
                 m_var_set.insert(cj);

src/util/lp/row_eta_matrix_def.h

@@ -33,7 +33,7 @@ void row_eta_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings &) {
         w_at_row += w[it.first] * it.second;
     }
     // w[m_row] = w_at_row;
-    // #ifdef LEAN_DEBUG
+    // #ifdef Z3DEBUG
     //         lp_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
     //         delete [] clone_w;
     // #endif
@@ -95,7 +95,7 @@ void row_eta_matrix<T, X>::apply_from_right(vector<T> & w) {
     for (auto & it : m_row_vector.m_data) {
         w[it.first] += w_row * it.second;
     }
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
     // delete clone_w;
 #endif
@@ -144,7 +144,7 @@ void row_eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
             }
         }
     }
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(vectors_are_equal(wcopy, w.m_data));
 
 #endif
@@ -165,7 +165,7 @@ void row_eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p
         columns.push_back(it.first);
     for (unsigned i = static_cast<unsigned>(columns.size()); i-- > 0;)
         m_row_vector.m_data[i].first = p.get_rev(columns[i]);
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     // lp_assert(deb == *this);
 #endif
 }

src/util/lp/scaler_def.h

@@ -220,18 +220,18 @@ template <typename T, typename X>    void scaler<T, X>::scale_row(unsigned i) {
     if (numeric_traits<T>::is_zero(row_max)) {
         return;
     }
-    if (numeric_traits<T>::get_double(row_max) < m_scaling_minimum) {
+    if (row_max < m_scaling_minimum) {
         do {
             alpha *= 2;
             row_max *= 2;
-        } while (numeric_traits<T>::get_double(row_max) < m_scaling_minimum);
+        } while (row_max < m_scaling_minimum);
         m_A.multiply_row(i, alpha);
         m_b[i] *= alpha;
-    } else if (numeric_traits<T>::get_double(row_max) > m_scaling_maximum) {
+    } else if (row_max > m_scaling_maximum) {
         do {
             alpha /= 2;
             row_max /= 2;
-        } while (numeric_traits<T>::get_double(row_max) > m_scaling_maximum);
+        } while (row_max > m_scaling_maximum);
         m_A.multiply_row(i, alpha);
         m_b[i] *= alpha;
     }
@@ -245,16 +245,16 @@ template <typename T, typename X>    void scaler<T, X>::scale_column(unsigned i)
         return; // the column has zeros only
     }
     std::cout << "f";
-    if (numeric_traits<T>::get_double(column_max) < m_scaling_minimum) {
+    if (column_max < m_scaling_minimum) {
         do {
             alpha *= 2;
             column_max *= 2;
-        } while (numeric_traits<T>::get_double(column_max) < m_scaling_minimum);
+        } while (column_max < m_scaling_minimum);
-    } else if (numeric_traits<T>::get_double(column_max) > m_scaling_maximum) {
+    } else if (column_max > m_scaling_maximum) {
         do {
             alpha /= 2;
             column_max /= 2;
-        } while (numeric_traits<T>::get_double(column_max) > m_scaling_maximum);
+        } while (column_max > m_scaling_maximum);
     } else {
         return;
     }

src/util/lp/sparse_matrix.cpp

@@ -82,7 +82,7 @@ template void sparse_matrix<mpq, numeric_pair<mpq>>::double_solve_U_y<mpq>(index
 template void sparse_matrix<mpq, numeric_pair<mpq> >::double_solve_U_y<numeric_pair<mpq> >(indexed_vector<numeric_pair<mpq>>&, const lp_settings&);
 template void lp::sparse_matrix<double, double>::solve_U_y_indexed_only<double>(lp::indexed_vector<double>&, const lp_settings&, vector<unsigned> &);
 template void lp::sparse_matrix<lp::mpq, lp::mpq>::solve_U_y_indexed_only<lp::mpq>(lp::indexed_vector<lp::mpq>&, const lp_settings &, vector<unsigned> &);
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
 template bool sparse_matrix<double, double>::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;
 template bool sparse_matrix<mpq, mpq>::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;
 template bool sparse_matrix<mpq, numeric_pair<mpq> >::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;

src/util/lp/sparse_matrix_def.h

@@ -490,7 +490,7 @@ void sparse_matrix<T, X>::solve_y_U_indexed(indexed_vector<T> & y, const lp_sett
     }
 
     lp_assert(y.is_OK());
-#if 0 && LEAN_DEBUG
+#if 0 && Z3DEBUG
     if (numeric_traits<T>::precise() == false)
         lp_assert(vectors_are_equal(ycopy, y.m_data));
 #endif
@@ -700,7 +700,7 @@ void sparse_matrix<T, X>::solve_U_y_indexed_only(indexed_vector<L> & y, const lp
     }
 
     lp_assert(y.is_OK());
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
      // dense_matrix<T,X> deb(this);
      // vector<T> clone_y(y.m_data);
      // deb.apply_from_left(clone_y);

src/util/lp/square_dense_submatrix_def.h

@@ -300,7 +300,7 @@ void square_dense_submatrix<T, X>::apply_from_left_to_vector(vector<L> & w) {
 }
 
 template <typename T, typename X>    bool square_dense_submatrix<T, X>::is_L_matrix() const {
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     lp_assert(m_row_permutation.is_identity());
     for (unsigned i = 0; i < m_parent->dimension(); i++) {
         if (i < m_index_start) {
@@ -341,7 +341,7 @@ template <typename T, typename X> void square_dense_submatrix<T, X>::apply_from_
         t[adjust_column_inverse(j)] = column_by_vector_product(j, w);
     }
     w = t;
-#ifdef LEAN_DEBUG
+#ifdef Z3DEBUG
     //  lp_assert(vector_are_equal<T>(deb_w, w));
 #endif
 }

src/util/lp/static_matrix.h

@@ -26,7 +26,6 @@ Revision History:
 #include "util/lp/sparse_vector.h"
 #include "util/lp/indexed_vector.h"
 #include "util/lp/permutation_matrix.h"
-#include "util/lp/linear_combination_iterator.h"
 #include <stack>
 namespace lp {
 
@@ -35,20 +34,26 @@ struct column_cell {
     unsigned m_offset;  // the offset of the element in the matrix row
     column_cell(unsigned i, unsigned offset) : m_i(i), m_offset(offset) {
     }
+    
 };
 
 template <typename T>
 struct row_cell {
     unsigned m_j; // points to the column
     unsigned m_offset; // offset in column
+    T        m_value;
     row_cell(unsigned j, unsigned offset, T const & val) : m_j(j), m_offset(offset), m_value(val) {
     }
     const T & get_val() const { return m_value;}
     T & get_val() { return m_value;}
     void set_val(const T& v) { m_value = v;  }
-    T m_value;
+    unsigned var() const { return m_j;}
+    const T & coeff() const { return m_value;}
 };
 
+template <typename T>
+using row_strip = vector<row_cell<T>>; 
+
 // each assignment for this matrix should be issued only once!!!
 template <typename T, typename X>
 class static_matrix
@@ -63,11 +68,10 @@ class static_matrix
     };
     std::stack<dim> m_stack;
 public:
-    typedef vector<row_cell<T>> row_strip;
     typedef vector<column_cell> column_strip;
     vector<int> m_vector_of_row_offsets;
     indexed_vector<T> m_work_vector;
-    vector<row_strip> m_rows;
+    vector<row_strip<T>> m_rows;
     vector<column_strip> m_columns;
     // starting inner classes
     class ref {
@@ -97,10 +101,6 @@ public:
         return m_rows[c.m_i][c.m_offset].get_val();
     }
 
-    row_cell<T> & get_row_cell(const column_cell & c) {
-        return m_rows[c.m_i][c.m_offset];
-    }
-
     column_cell & get_column_cell(const row_cell<T> &rc) {
         return m_columns[rc.m_j][rc.m_offset];
     }
@@ -132,7 +132,7 @@ public:
     void add_columns_at_the_end(unsigned delta);
     void add_new_element(unsigned i, unsigned j, const T & v);
 
-    void add_row() {m_rows.push_back(row_strip());}
+    void add_row() {m_rows.push_back(row_strip<T>());}
     void add_column() {
         m_columns.push_back(column_strip());
         m_vector_of_row_offsets.push_back(-1);
@@ -195,7 +195,7 @@ public:
     //
     void fix_row_indices_in_each_column_for_crossed_row(unsigned k);
 
-    void cross_out_row_from_columns(unsigned k, row_strip & row);
+    void cross_out_row_from_columns(unsigned k, row_strip<T> & row);
 
     void cross_out_row_from_column(unsigned col, unsigned k);
 
@@ -294,7 +294,7 @@ public:
 
     // pivot row i to row ii
     bool pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned);
-    void scan_row_ii_to_offset_vector(const row_strip & rvals);
+    void scan_row_ii_to_offset_vector(const row_strip<T> & rvals);
 
     void transpose_rows(unsigned i, unsigned ii) {
         auto t = m_rows[i];
@@ -313,46 +313,54 @@ public:
         }
     
     }
+    void fill_last_row_with_pivoting_loop_block(unsigned j, const vector<int> & basis_heading) {
+        int row_index = basis_heading[j];
+        if (row_index < 0)
+            return;
+        T & alpha = m_work_vector[j]; // the pivot alpha
+        if (is_zero(alpha))
+            return;
+        
+        for (const auto & c : m_rows[row_index]) {
+            if (c.m_j == j) {
+                continue;
+            }
+            T & wv = m_work_vector.m_data[c.m_j];
+            bool was_zero = is_zero(wv);
+            wv -= alpha * c.m_value;
+            if (was_zero)
+                m_work_vector.m_index.push_back(c.m_j);
+            else {
+                if (is_zero(wv)) {
+                    m_work_vector.erase_from_index(c.m_j);
+                }
+            }
+        }
+        alpha = zero_of_type<T>();
+        m_work_vector.erase_from_index(j);
+    }
 
-    void fill_last_row_with_pivoting(linear_combination_iterator<T> & it, const vector<int> & basis_heading) {
+    
+    
+    template <typename term>
+    void fill_last_row_with_pivoting(const term& row,
+                                     unsigned bj, // the index of the basis column
+                                     const vector<int> & basis_heading) {
         lp_assert(numeric_traits<T>::precise());
         lp_assert(row_count() > 0);
         m_work_vector.resize(column_count());
         T a;
-        unsigned j;
+         // we use the form -it + 1 = 0
-        while (it.next(a, j)) {
+        m_work_vector.set_value(one_of_type<T>(), bj);
-            m_work_vector.set_value(-a, j); // we use the form -it + 1 = 0
+        for (auto p : row) {
+            m_work_vector.set_value(-p.coeff(), p.var());
             // but take care of the basis 1 later
         }
     
-        it.reset();
+        // now iterate with pivoting
-        // not iterate with pivoting
+        fill_last_row_with_pivoting_loop_block(bj, basis_heading);
-        while (it.next(j)) {
+        for (auto p : row) {
-            int row_index = basis_heading[j];
+            fill_last_row_with_pivoting_loop_block(p.var(), basis_heading);
-            if (row_index < 0)
-                continue;
-
-            T & alpha = m_work_vector[j]; // the pivot alpha
-            if (is_zero(alpha))
-                continue;
-        
-            for (const auto & c : m_rows[row_index]) {
-                if (c.m_j == j) {
-                    continue;
-                }
-                T & wv = m_work_vector.m_data[c.m_j];
-                bool was_zero = is_zero(wv);
-                wv -= alpha * c.m_value;
-                if (was_zero)
-                    m_work_vector.m_index.push_back(c.m_j);
-                else {
-                    if (is_zero(wv)) {
-                        m_work_vector.erase_from_index(c.m_j);
-                    }
-                }
-            }
-            alpha = zero_of_type<T>();
-            m_work_vector.erase_from_index(j);
         }
         lp_assert(m_work_vector.is_OK());
         unsigned last_row = row_count() - 1;
@@ -382,5 +390,66 @@ public:
         }
         return ret;
     }
+
+    struct column_cell_plus {
+        const column_cell & m_c;
+        const static_matrix& m_A;
+        // constructor
+        column_cell_plus(const column_cell & c, const static_matrix& A) :
+            m_c(c), m_A(A) {}
+        unsigned var() const { return m_c.m_i; }
+        const T & coeff() const { return m_A.m_rows[var()][m_c.m_offset].get_val(); }
+        
+    };
+    
+    struct column_container {
+        unsigned              m_j; // the column index
+        const static_matrix & m_A;
+        column_container(unsigned j, const static_matrix& A) : m_j(j), m_A(A) {
+        }
+        struct const_iterator {
+            // fields
+            const column_cell *m_c;
+            const static_matrix& m_A;
+
+            //typedefs
+            
+            
+            typedef const_iterator self_type;
+            typedef column_cell_plus value_type;
+            typedef const column_cell_plus reference;
+            //            typedef const column_cell* pointer;
+            typedef int difference_type;
+            typedef std::forward_iterator_tag iterator_category;
+
+            reference operator*() const {
+                return column_cell_plus(*m_c, m_A);
+            }        
+            self_type operator++() {  self_type i = *this; m_c++; return i;  }
+            self_type operator++(int) { m_c++; return *this; }
+
+            const_iterator(const column_cell* it, const static_matrix& A) :
+                m_c(it),
+                m_A(A)
+            {}
+            bool operator==(const self_type &other) const {
+                return m_c == other.m_c;
+            }
+            bool operator!=(const self_type &other) const { return !(*this == other); }
+        };
+
+        const_iterator begin() const {
+            return const_iterator(m_A.m_columns[m_j].begin(), m_A);
+        }
+        
+        const_iterator end() const {
+            return const_iterator(m_A.m_columns[m_j].end(), m_A);
+        }
+    };
+
+    column_container column(unsigned j) const {
+        return column_container(j, *this);
+    }
+    
 };
 }

src/util/lp/static_matrix_def.h

@@ -27,7 +27,7 @@ template <typename T, typename X>
 void  static_matrix<T, X>::init_row_columns(unsigned m, unsigned n) {
     lp_assert(m_rows.size() == 0 && m_columns.size() == 0);
     for (unsigned i = 0; i < m; i++){
-        m_rows.push_back(row_strip());
+        m_rows.push_back(row_strip<T>());
     }
     for (unsigned j = 0; j < n; j++){
         m_columns.push_back(column_strip());
@@ -35,7 +35,7 @@ void  static_matrix<T, X>::init_row_columns(unsigned m, unsigned n) {
 }
 
 
-template <typename T, typename X> void static_matrix<T, X>::scan_row_ii_to_offset_vector(const row_strip & rvals) {
+template <typename T, typename X> void static_matrix<T, X>::scan_row_ii_to_offset_vector(const row_strip<T> & rvals) {
     for (unsigned j = 0; j < rvals.size(); j++)
         m_vector_of_row_offsets[rvals[j].m_j] = j;
 }
@@ -166,7 +166,7 @@ template <typename T, typename X>
 std::set<std::pair<unsigned, unsigned>>  static_matrix<T, X>::get_domain() {
     std::set<std::pair<unsigned, unsigned>> ret;
     for (unsigned i = 0; i < m_rows.size(); i++) {
-        for (auto it : m_rows[i]) {
+        for (auto &it : m_rows[i]) {
             ret.insert(std::make_pair(i, it.m_j));
         }
     }
@@ -296,7 +296,7 @@ template <typename T, typename X>    void static_matrix<T, X>::fix_row_indices_i
     }
 }
 
-template <typename T, typename X>    void static_matrix<T, X>::cross_out_row_from_columns(unsigned k, row_strip & row) {
+template <typename T, typename X>    void static_matrix<T, X>::cross_out_row_from_columns(unsigned k, row_strip<T> & row) {
     for (auto & t : row) {
         cross_out_row_from_column(t.m_j, k);
     }

src/util/lp/test_bound_analyzer.h

@@ -17,9 +17,9 @@ Revision History:
 
 
 --*/
+#if 0
 #pragma once
 #include "util/vector.h"
-#include "util/lp/linear_combination_iterator.h"
 #include "util/lp/implied_bound.h"
 #include "util/lp/lp_settings.h"
 #include <functional>
@@ -273,3 +273,4 @@ public :
 };
 
 }
+#endif