diff --git a/src/math/lp/.clang-format b/src/math/lp/.clang-format
new file mode 100644
index 000000000..f5c8a41b2
--- /dev/null
+++ b/src/math/lp/.clang-format
@@ -0,0 +1,4 @@
+BasedOnStyle: Google
+IndentWidth: 4
+ColumnLimit: 0
+NamespaceIndentation: All
\ No newline at end of file
diff --git a/src/math/lp/explanation.h b/src/math/lp/explanation.h
index d2e7edc33..b7f721a30 100644
--- a/src/math/lp/explanation.h
+++ b/src/math/lp/explanation.h
@@ -48,6 +48,15 @@ public:
         SASSERT(m_vector.empty());
         m_set.insert(j);
     }
+
+    void remove(constraint_index j) {
+        m_set.remove(j);
+        unsigned i = 0;
+        for (auto& p : m_vector) 
+            if (p.first != j)
+                m_vector[i++] = p;
+        m_vector.shrink(i);
+    }
     
     void add_expl(const explanation& e) {
         if (e.m_vector.empty()) {
diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp
index 6c34ce16f..af4488162 100644
--- a/src/math/lp/int_solver.cpp
+++ b/src/math/lp/int_solver.cpp
@@ -2,8 +2,7 @@
   Copyright (c) 2017 Microsoft Corporation
   Author: Lev Nachmanson
 */
-
-#include <utility>
+// clang-format off
 #include "math/lp/int_solver.h"
 #include "math/lp/lar_solver.h"
 #include "math/lp/lp_utils.h"
@@ -14,56 +13,216 @@
 
 namespace lp {
 
-int_solver::patcher::patcher(int_solver& lia):
-    lia(lia),
-    lra(lia.lra),
-    lrac(lia.lrac),
-    m_num_nbasic_patches(0),
-    m_patch_cost(0),
-    m_next_patch(0),
-    m_delay(0)
-{}
+    int_solver::patcher::patcher(int_solver& lia):
+        lia(lia),
+        lra(lia.lra),
+        lrac(lia.lrac)
+    {}
+    
+    void int_solver::patcher::remove_fixed_vars_from_base() {
+        unsigned num = lra.A_r().column_count();
+        for (unsigned v = 0; v < num; v++) {
+            if (!lia.is_base(v) || !lia.is_fixed(v))
+                continue;
+            auto const & r = lra.basic2row(v);
+            for (auto const& c : r) {
+                if (c.var() != v && !lia.is_fixed(c.var())) {
+                    lra.pivot(c.var(), v); 
+                    break;
+                }
+            }        
+        }
+    }
 
-bool int_solver::patcher::should_apply() {
-#if 1
-    return true;
-#else
-    if (m_delay == 0) {
+
+    unsigned int_solver::patcher::count_non_int() {
+        unsigned non_int = 0;
+        for (auto j : lra.r_basis()) 
+            if (lra.column_is_int(j) && !lra.column_value_is_int(j))
+                ++non_int;
+        return non_int;
+    }
+
+    lia_move int_solver::patcher::patch_basic_columns() {
+        remove_fixed_vars_from_base();
+        lia.settings().stats().m_patches++;
+        lp_assert(lia.is_feasible());
+        
+       // unsigned non_int_before, non_int_after;
+
+       // non_int_before = count_non_int();
+
+
+        // unsigned num = lra.A_r().column_count();
+        for (unsigned j : lra.r_basis()) 
+            if (!lra.get_value(j).is_int())
+                patch_basic_column(j);
+        // non_int_after = count_non_int();
+        // verbose_stream() << non_int_before << " -> " << non_int_after << "\n";
+
+        if (!lia.has_inf_int()) {
+            lia.settings().stats().m_patches_success++;
+            return lia_move::sat;
+        }
+        return lia_move::undef;
+    }
+
+    // clang-format on
+    /**
+     * \brief find integral and minimal, in the absolute values, deltas such that x - alpha*delta is integral too.
+    */
+    bool get_patching_deltas(const rational& x, const rational& alpha,
+                             rational& delta_0, rational& delta_1) {
+        auto a1 = numerator(alpha);
+        auto a2 = denominator(alpha);
+        auto x1 = numerator(x);
+        auto x2 = denominator(x);
+        if (!divides(x2, a2))
+            return false;
+
+        // delta has to be integral.
+        // We need to find delta such that x1/x2 + (a1/a2)*delta is integral (we are going to flip the delta sign later).
+        // Then a2*x1/x2 + a1*delta is integral, but x2 and x1 are coprime:
+        // that means that t = a2/x2 is
+        // integral. We established that a2 = x2*t Then x1 + a1*delta*(x2/a2)  = x1
+        // + a1*(delta/t)  is integral. Taking into account that t and a1 are
+        // coprime we have delta = t*k, where k is an integer.
+        rational t = a2 / x2;
+        // std::cout << "t = " << t << std::endl;
+        // Now we have x1/x2 + (a1/x2)*k is integral, or (x1 + a1*k)/x2 is integral.
+        // It is equivalent to x1 + a1*k = x2*m, where m is an integer
+        // We know that a2 and a1 are coprime, and x2 divides a2, so x2 and a1 are
+        // coprime. We can find u and v such that u*a1 + v*x2 = 1.
+        rational u, v;
+        gcd(a1, x2, u, v);
+        lp_assert(gcd(a1, x2, u, v).is_one());
+        // std::cout << "u = " << u << ", v = " << v << std::endl;
+        // std::cout << "x= " << (x1 / x2) << std::endl;
+        // std::cout << "x + (a1 / a2) * (-u * t) * x1 = "
+        //           << x + (a1 / a2) * (-u * t) * x1 << std::endl;
+        lp_assert((x + (a1 / a2) * (-u * t) * x1).is_int());
+        // 1 = (u- l*x2 ) * a1 + (v + l*a1)*x2, for every integer l.
+        rational l_low, l_high;
+        auto sign = u.is_pos() ? 1 : -1;
+        auto p = sign * floor(abs(u / x2));
+        auto p_ = p + sign;
+        lp_assert(p * x2 <= u && u <= p_ * x2 || p * x2 >= u && u >= p_ * x2);
+        // std::cout << "u = " << u << ", v = " << v << std::endl;
+        // std::cout << "p = " << p << ", p_ = " << p_ << std::endl;
+        // std::cout << "u - p*x2 = " << u - p * x2 << ", u - p_*x2 = " << u - p_ * x2
+        //           << std::endl;
+        mpq d_0 = (u - p * x2) * t * x1;
+        mpq d_1 = (u - p_ * x2) * t * x1;
+        if (abs(d_0) < abs(d_1)) {
+            delta_0 = d_0;
+            delta_1 = d_1;
+        } else {
+            delta_0 = d_1;
+            delta_1 = d_0;
+        }
+        return true;
+        // std::cout << "delta_0 = " << delta_0 << std::endl;
+        // std::cout << "delta_1 = " << delta_1 << std::endl;
+    }
+    // clang-format off
+    /**
+     * \brief try to patch the basic column v
+     */
+    bool int_solver::patcher::patch_basic_column_on_row_cell(unsigned v, row_cell<mpq> const& c) {
+        if (v == c.var())
+            return false;
+        if (!lra.column_is_int(c.var())) // could use real to patch integer
+            return false;
+        if (c.coeff().is_int())
+            return false;
+        mpq a = fractional_part(c.coeff());
+        mpq r = fractional_part(lra.get_value(v));
+        lp_assert(0 < r && r < 1);
+        lp_assert(0 < a && a < 1);
+        mpq delta_0, delta_1;
+        if (!get_patching_deltas(r, a, delta_0, delta_1))
+            return false;
+        
+        if (try_patch_column(v, c.var(), delta_0))
+            return true;
+
+        if (try_patch_column(v, c.var(), delta_1)) 
+            return true;
+       
+        return false;
+    }
+
+    bool int_solver::patcher::try_patch_column(unsigned v, unsigned j, mpq const& delta) {
+        const auto & A = lra.A_r();
+        if (delta < 0) {
+            if (lia.has_lower(j) && lia.get_value(j) + impq(delta) < lra.get_lower_bound(j))
+                return false;
+        }
+        else {
+            if (lia.has_upper(j) && lia.get_value(j) + impq(delta) > lra.get_upper_bound(j))
+                return false;
+        }
+        for (auto const& c : A.column(j)) {
+            unsigned row_index = c.var();
+            unsigned i = lrac.m_r_basis[row_index];
+            auto old_val = lia.get_value(i);
+            auto new_val = old_val - impq(c.coeff()*delta);
+            if (lia.has_lower(i) && new_val < lra.get_lower_bound(i))
+                return false;
+            if (lia.has_upper(i) && new_val > lra.get_upper_bound(i))
+                return false;
+            if (old_val.is_int() && !new_val.is_int()){
+                return false; // do not waste resources on this case
+            }
+            lp_assert(i != v || new_val.is_int())
+        }
+        
+        lra.set_value_for_nbasic_column(j, lia.get_value(j) + impq(delta));
+         
         return true;
     }
-    --m_delay;
-    return false;
-#endif
-}
+    
+    void int_solver::patcher::patch_basic_column(unsigned v) {
+        SASSERT(!lia.is_fixed(v));
+        for (auto const& c : lra.basic2row(v))
+            if (patch_basic_column_on_row_cell(v, c))
+                return;                                       
+    }
 
-lia_move int_solver::patcher::operator()() {
-    return patch_nbasic_columns();
-}
+    lia_move int_solver::patcher::patch_nbasic_columns() {
+        remove_fixed_vars_from_base();
+        lia.settings().stats().m_patches++;
+        lp_assert(lia.is_feasible());
+        m_patch_success = 0;
+        m_patch_fail = 0;
+        m_num_ones = 0;
+        m_num_divides = 0;
+        //unsigned non_int_before = count_non_int();
 
-lia_move int_solver::patcher::patch_nbasic_columns() {
-    lia.settings().stats().m_patches++;
-    lp_assert(lia.is_feasible());
-    m_num_nbasic_patches = 0;
-    m_patch_cost = 0;
-    for (unsigned j : lia.lrac.m_r_nbasis) {
-        patch_nbasic_column(j);
+        unsigned num = lra.A_r().column_count();
+        for (unsigned v = 0; v < num; v++) {
+            if (lia.is_base(v))
+                continue;
+            patch_nbasic_column(v);
+        }
+        unsigned num_fixed = 0;
+        for (unsigned v = 0; v < num; v++) 
+            if (lia.is_fixed(v))
+                ++num_fixed;
+            
+        lp_assert(lia.is_feasible());
+        //verbose_stream() << "patch " << m_patch_success << " fails " << m_patch_fail << " ones " << m_num_ones << " divides " << m_num_divides << " num fixed " << num_fixed << "\n";
+        //lra.display(verbose_stream());
+        //exit(0);
+        //unsigned non_int_after = count_non_int();
+
+        // verbose_stream() << non_int_before << " -> " << non_int_after << "\n";
+        if (!lia.has_inf_int()) {
+            lia.settings().stats().m_patches_success++;
+            return lia_move::sat;
+        }
+        return lia_move::undef;
     }
-    lp_assert(lia.is_feasible());
-    if (!lia.has_inf_int()) {
-        lia.settings().stats().m_patches_success++;
-        m_delay = 0;
-        m_next_patch = 0;
-        return lia_move::sat;
-    }
-    if (m_patch_cost > 0 && m_num_nbasic_patches * 10 < m_patch_cost) {
-        m_delay = std::min(20u, m_next_patch++);
-    }
-    else {
-        m_delay = 0;
-        m_next_patch = 0;
-    }
-    return lia_move::undef;
-}
 
 void int_solver::patcher::patch_nbasic_column(unsigned j) {
     impq & val = lrac.m_r_x[j];
@@ -71,17 +230,48 @@ void int_solver::patcher::patch_nbasic_column(unsigned j) {
     impq l, u;
     mpq m;
     bool has_free = lia.get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m);
-    m_patch_cost += lra.A_r().number_of_non_zeroes_in_column(j);
-    if (!has_free) {
+    if (!has_free) 
         return;
-    }
     bool m_is_one = m.is_one();
     bool val_is_int = lia.value_is_int(j);
 
+#if 0
+const auto & A = lra.A_r();
+#endif
     // check whether value of j is already a multiple of m.
     if (val_is_int && (m_is_one || (val.x / m).is_int())) {
+        if (m_is_one)
+            ++m_num_ones;
+        else
+            ++m_num_divides;
+#if 0
+        for (auto c : A.column(j)) {
+            unsigned row_index = c.var();
+            unsigned i = lrac.m_r_basis[row_index];
+            if (!lia.get_value(i).is_int() ||
+                (lia.has_lower(i) && lia.get_value(i) < lra.get_lower_bound(i)) ||
+                (lia.has_upper(i) && lia.get_value(i) > lra.get_upper_bound(i))) {
+                verbose_stream() << "skip " << j << " " << m << ": ";
+                lia.display_row(verbose_stream(), A.m_rows[row_index]);
+                verbose_stream() << "\n";                    
+            }
+        }
+#endif
         return;
     }
+
+#if 0
+    if (!m_is_one) {
+        // lia.display_column(verbose_stream(), j);
+        for (auto c : A.column(j)) {
+            continue;
+            unsigned row_index = c.var();
+            lia.display_row(verbose_stream(), A.m_rows[row_index]);
+            verbose_stream() << "\n";
+        }
+    }
+#endif
+
     TRACE("patch_int",
           tout << "TARGET j" << j << " -> [";
           if (inf_l) tout << "-oo"; else tout << l;
@@ -89,9 +279,33 @@ void int_solver::patcher::patch_nbasic_column(unsigned j) {
           if (inf_u) tout << "oo"; else tout << u;
           tout << "]";
           tout << ", m: " << m << ", val: " << val << ", is_int: " << lra.column_is_int(j) << "\n";);
+
+#if 0
+    verbose_stream() << "path " << m << " ";
+    if (!inf_l) verbose_stream() << "infl " << l.x << " ";
+    if (!inf_u) verbose_stream() << "infu " << u.x << " ";
+    verbose_stream() << "\n";
     if (m.is_big() || (!inf_l && l.x.is_big()) || (!inf_u && u.x.is_big())) {
         return;
     }
+#endif
+
+#if 0
+    verbose_stream() << "TARGET v" << j << " -> [";
+    if (inf_l) verbose_stream() << "-oo"; else verbose_stream() << ceil(l.x) << " " << l << "\n";
+    verbose_stream() << ", ";
+    if (inf_u) verbose_stream() << "oo"; else verbose_stream() << floor(u.x) << " " << u << "\n";
+    verbose_stream() << "]";
+    verbose_stream() << ", m: " << m << ", val: " << val << ", is_int: " << lra.column_is_int(j) << "\n";
+#endif
+
+#if 0
+    if (!inf_l)
+        l = impq(ceil(l));
+    if (!inf_u)
+        u = impq(floor(u));
+#endif
+    
     if (!inf_l) {
         l = impq(m_is_one ? ceil(l) : m * ceil(l / m));
         if (inf_u || l <= u) {
@@ -99,8 +313,23 @@ void int_solver::patcher::patch_nbasic_column(unsigned j) {
             lra.set_value_for_nbasic_column(j, l);
         }
         else {
-            --m_num_nbasic_patches;
+            //verbose_stream() << "fail: " << j << " " << m << "\n";
+            ++m_patch_fail;
             TRACE("patch_int", tout << "not patching " << l << "\n";);
+#if 0
+            verbose_stream() << "not patched\n";
+            for (auto c : A.column(j)) {
+                unsigned row_index = c.var();
+                unsigned i = lrac.m_r_basis[row_index];
+                if (!lia.get_value(i).is_int() ||
+                    (lia.has_lower(i) && lia.get_value(i) < lra.get_lower_bound(i)) ||
+                    (lia.has_upper(i) && lia.get_value(i) > lra.get_upper_bound(i))) {
+                    lia.display_row(verbose_stream(), A.m_rows[row_index]);
+                    verbose_stream() << "\n";                    
+                }
+            }
+#endif
+            return;
         }
     }
     else if (!inf_u) {
@@ -112,7 +341,21 @@ void int_solver::patcher::patch_nbasic_column(unsigned j) {
         lra.set_value_for_nbasic_column(j, impq(0));
         TRACE("patch_int", tout << "patching with 0\n";);
     }
-    ++m_num_nbasic_patches;
+    ++m_patch_success;
+#if 0
+    verbose_stream() << "patched " << j << "\n";
+    for (auto c : A.column(j)) {
+        unsigned row_index = c.var();
+        unsigned i = lrac.m_r_basis[row_index];
+        if (!lia.get_value(i).is_int() ||
+            (lia.has_lower(i) && lia.get_value(i) < lra.get_lower_bound(i)) ||
+            (lia.has_upper(i) && lia.get_value(i) > lra.get_upper_bound(i))) {
+            lia.display_row(verbose_stream(), A.m_rows[row_index]);
+            verbose_stream() << "\n";                    
+        }
+    }
+#endif
+    
 }
 
 int_solver::int_solver(lar_solver& lar_slv) :
@@ -326,7 +569,7 @@ static void set_upper(impq & u, bool & inf_u, impq const & v) {
 
 // this function assumes that all basic columns dependend on j are feasible
 bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m) {
-    if (lrac.m_r_heading[j] >= 0) // the basic var
+    if (lrac.m_r_heading[j] >= 0 || is_fixed(j)) // basic or fixed var
         return false;
 
     TRACE("random_update", display_column(tout, j) << ", is_int = " << column_is_int(j) << "\n";);
@@ -361,10 +604,9 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
         unsigned i = lrac.m_r_basis[row_index];
         impq const & xi = get_value(i);
         lp_assert(lrac.m_r_solver.column_is_feasible(i));
-        if (column_is_int(i) && !a.is_int())
+        if (column_is_int(i) && !a.is_int() && xi.is_int()) 
             m = lcm(m, denominator(a));
-
-
+        
         if (!inf_l && !inf_u) {
             if (l == u) 
                 continue;            
@@ -372,15 +614,15 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
 
         if (a.is_neg()) {
             if (has_lower(i)) 
-                set_lower(l, inf_l, delta(a, xi, lrac.m_r_lower_bounds()[i]));
+                set_lower(l, inf_l, delta(a, xi, lra.get_lower_bound(i)));
             if (has_upper(i)) 
-                set_upper(u, inf_u, delta(a, xi, lrac.m_r_upper_bounds()[i]));
+                set_upper(u, inf_u, delta(a, xi, lra.get_upper_bound(i)));
         }
         else {
             if (has_upper(i)) 
-                set_lower(l, inf_l, delta(a, xi, lrac.m_r_upper_bounds()[i]));
+                set_lower(l, inf_l, delta(a, xi, lra.get_upper_bound(i)));
             if (has_lower(i)) 
-                set_upper(u, inf_u, delta(a, xi, lrac.m_r_lower_bounds()[i]));
+                set_upper(u, inf_u, delta(a, xi, lra.get_lower_bound(i)));
         }
     }
 
@@ -479,14 +721,11 @@ bool int_solver::at_upper(unsigned j) const {
 }
 
 std::ostream & int_solver::display_row(std::ostream & out, lp::row_strip<rational> const & row) const {
-bool first = true;
+    bool first = true;
     auto & rslv = lrac.m_r_solver;
-for (const auto &c : row)
-    {
-        if (is_fixed(c.var()))
-        {
-            if (!get_value(c.var()).is_zero())
-            {
+    for (const auto &c : row) {
+        if (is_fixed(c.var())) {
+            if (!get_value(c.var()).is_zero()) {
                 impq val = get_value(c.var()) * c.coeff();
                 if (!first && val.is_pos())
                     out << "+";
@@ -505,17 +744,11 @@ for (const auto &c : row)
         }
         else if (c.coeff().is_minus_one())
             out << "-";
-        else
-        {
-            if (c.coeff().is_pos())
-            {
-                if (!first)
-                    out << "+";
-            }
+        else {
+            if (c.coeff().is_pos() && !first)
+                out << "+";
             if (c.coeff().is_big())
-            {
                 out << " b*";
-            }
             else
                 out << c.coeff();
         }
@@ -523,8 +756,7 @@ for (const auto &c : row)
         first = false;
     }
     out << "\n";
-    for (const auto &c : row)
-    {
+    for (const auto &c : row) {
         if (is_fixed(c.var()))
             continue;
         rslv.print_column_info(c.var(), out);
@@ -533,14 +765,13 @@ for (const auto &c : row)
     }
     return out;
 }
+    
 std::ostream& int_solver::display_row_info(std::ostream & out, unsigned row_index) const  {    
     auto & rslv = lrac.m_r_solver;
     auto const& row = rslv.m_A.m_rows[row_index];
     return display_row(out, row);
 }
 
-
-
 bool int_solver::shift_var(unsigned j, unsigned range) {
     if (is_fixed(j) || is_base(j))
         return false;
@@ -549,11 +780,13 @@ bool int_solver::shift_var(unsigned j, unsigned range) {
     bool inf_l = false, inf_u = false;
     impq l, u;
     mpq m;
-    VERIFY(get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m) || settings().get_cancel_flag());
+    if (!get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m))
+        return false;
     if (settings().get_cancel_flag())
         return false;
     const impq & x = get_value(j);
     // x, the value of j column, might be shifted on a multiple of m
+
     if (inf_l && inf_u) {
         impq new_val = m * impq(random() % (range + 1)) + x;
         lra.set_value_for_nbasic_column(j, new_val);
@@ -570,6 +803,7 @@ bool int_solver::shift_var(unsigned j, unsigned range) {
     if (!inf_l && !inf_u && l >= u)
         return false;
 
+
     if (inf_u) {
         SASSERT(!inf_l);
         impq new_val = x + m * impq(random() % (range + 1));
@@ -640,9 +874,6 @@ int int_solver::select_int_infeasible_var() {
     unsigned n = 0;
     lar_core_solver & lcs = lra.m_mpq_lar_core_solver;
     unsigned prev_usage = 0; // to quiet down the compile
-    unsigned k = 0;
-    unsigned usage;
-    unsigned j;
 
     enum state { small_box, is_small_value,  any_value, not_found };
     state st = not_found;
@@ -650,11 +881,10 @@ int int_solver::select_int_infeasible_var() {
     // 1. small box
     // 2. small value
     // 3. any value
-    for (; k < lra.r_basis().size(); k++) {
-        j = lra.r_basis()[k];
+    for (unsigned j : lra.r_basis()) {
         if (!column_is_int_inf(j))
             continue;
-        usage = lra.usage_in_terms(j);
+        unsigned usage = lra.usage_in_terms(j);
         if (is_boxed(j) && (new_range = lcs.m_r_upper_bounds()[j].x - lcs.m_r_lower_bounds()[j].x - rational(2*usage)) <= small_value) {
             SASSERT(!is_fixed(j));
             if (st != small_box) {
@@ -688,12 +918,12 @@ int int_solver::select_int_infeasible_var() {
             continue;
         SASSERT(st == not_found || st == any_value);
         st = any_value;
-        if (n == 0 /*|| usage > prev_usage*/) {
+        if (n == 0 || usage > prev_usage) {
             result = j;
             prev_usage = usage;
             n = 1;
         } 
-        else if (usage > 0 && /*usage == prev_usage && */ (random() % (++n) == 0)) 
+        else if (usage > 0 && usage == prev_usage && (random() % (++n) == 0)) 
             result = j;
     }
     
diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h
index 822e1cf1e..fad040965 100644
--- a/src/math/lp/int_solver.h
+++ b/src/math/lp/int_solver.h
@@ -17,6 +17,7 @@ Revision History:
 
 
 --*/
+// clang-format off
 #pragma once
 #include "math/lp/lp_settings.h"
 #include "math/lp/static_matrix.h"
@@ -44,17 +45,23 @@ class int_solver {
         int_solver&         lia;
         lar_solver&         lra;
         lar_core_solver&    lrac;
-        unsigned            m_num_nbasic_patches;
-        unsigned            m_patch_cost;
-        unsigned            m_next_patch;
-        unsigned            m_delay;
+        unsigned            m_patch_success = 0;
+        unsigned            m_patch_fail = 0;
+        unsigned            m_num_ones = 0;
+        unsigned            m_num_divides = 0;
     public:
         patcher(int_solver& lia);
-        bool should_apply();
-        lia_move operator()();
+        bool should_apply() const { return true; }
+        lia_move operator()() { return patch_basic_columns(); }
         void patch_nbasic_column(unsigned j);
+        bool patch_basic_column_on_row_cell(unsigned v, row_cell<mpq> const& c);
+        void patch_basic_column(unsigned j);
+        bool try_patch_column(unsigned v, unsigned j, mpq const& delta);
+        unsigned count_non_int();
     private:
+        void remove_fixed_vars_from_base();
         lia_move patch_nbasic_columns();
+        lia_move patch_basic_columns();
     };
 
     lar_solver&         lra;
diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h
index 06ef4d50b..10cd53141 100644
--- a/src/math/lp/lar_core_solver.h
+++ b/src/math/lp/lar_core_solver.h
@@ -93,6 +93,8 @@ public:
 
     void solve();
 
+    void pivot(int entering, int leaving) { m_r_solver.pivot(entering, leaving); }
+    
     bool lower_bounds_are_set() const { return true; }
 
     const indexed_vector<mpq> & get_pivot_row() const {
diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index 5f544dff1..9fa9bc540 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -2,7 +2,7 @@
   Copyright (c) 2017 Microsoft Corporation
   Author: Nikolaj Bjorner, Lev Nachmanson
 */
-
+// clang-format off
 #include "math/lp/lar_solver.h"
 #include "smt/params/smt_params_helper.hpp"
 
@@ -41,7 +41,6 @@ namespace lp {
         for (auto t : m_terms)
             delete t;
     }
-
     
     bool lar_solver::sizes_are_correct() const {
         lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_column_types.size());
@@ -50,7 +49,6 @@ namespace lp {
         return true;
     }
 
-
     std::ostream& lar_solver::print_implied_bound(const implied_bound& be, std::ostream& out) const {
         out << "implied bound\n";
         unsigned v = be.m_j;
@@ -1354,8 +1352,8 @@ namespace lp {
     }
 
     bool lar_solver::term_is_int(const vector<std::pair<mpq, unsigned int>>& coeffs) const {
-        for (auto const& p : coeffs)
-            if (!(column_is_int(p.second) && p.first.is_int()))
+        for (auto const& [coeff, v] : coeffs)
+            if (!(column_is_int(v) && coeff.is_int()))
                 return false;
         return true;
     }
@@ -1386,7 +1384,7 @@ namespace lp {
     // the lower/upper bound is not strict.
     // the LP obtained by making the bound strict is infeasible
     // -> the column has to be fixed
-    bool lar_solver::is_fixed_at_bound(column_index const& j) {
+    bool lar_solver::is_fixed_at_bound(column_index const& j, vector<std::tuple<explanation, column_index, bool, mpq>>& bounds) {
         if (column_is_fixed(j))
             return false;
         mpq val;
@@ -1395,50 +1393,56 @@ namespace lp {
         lp::lconstraint_kind k;
         if (column_has_upper_bound(j) &&
             get_upper_bound(j).x == val) {
-            verbose_stream() << "check upper " << j << "\n";
-            push();
-            if (column_is_int(j))
-                k = LE, val -= 1;
-            else
-                k = LT;
-            auto ci = mk_var_bound(j, k, val);
-            update_column_type_and_bound(j, k, val, ci);
+            push();   
+            k = column_is_int(j) ? LE : LT;
+            auto ci = mk_var_bound(j, k, column_is_int(j) ? val - 1 : val);
+            update_column_type_and_bound(j, k, column_is_int(j) ? val - 1 : val, ci);
             auto st = find_feasible_solution();
+            bool infeasible = st == lp_status::INFEASIBLE;
+            if (infeasible) {
+                explanation exp;
+                get_infeasibility_explanation(exp);
+                unsigned_vector cis;
+                exp.remove(ci);
+                verbose_stream() << "tight upper bound " << j << " " << val << "\n";
+                bounds.push_back({exp, j, true, val});
+            }
             pop(1);
-            return st == lp_status::INFEASIBLE;            
+            return infeasible;            
         }
         if (column_has_lower_bound(j) &&
             get_lower_bound(j).x == val) {
-            verbose_stream() << "check lower " << j << "\n";
             push();
-            if (column_is_int(j))
-                k = GE, val += 1;
-            else
-                k = GT;
-            auto ci = mk_var_bound(j, k, val);
-            update_column_type_and_bound(j, k, val, ci);
+            k = column_is_int(j) ? GE : GT;
+            auto ci = mk_var_bound(j, k, column_is_int(j) ? val + 1 : val);
+            update_column_type_and_bound(j, k, column_is_int(j) ? val + 1 : val, ci);
             auto st = find_feasible_solution();
+            bool infeasible = st == lp_status::INFEASIBLE;
+            if (infeasible) {
+                explanation exp;
+                get_infeasibility_explanation(exp);
+                exp.remove(ci);
+                verbose_stream() << "tight lower bound " << j << " " << val << "\n";
+                bounds.push_back({exp, j, false, val});
+            }
             pop(1);
-            return st == lp_status::INFEASIBLE;                        
+            return infeasible;
         }
         
         return false;
     }
 
-    bool lar_solver::has_fixed_at_bound() {
+    bool lar_solver::has_fixed_at_bound(vector<std::tuple<explanation, column_index, bool, mpq>>& bounds) {
         verbose_stream() << "has-fixed-at-bound\n";
-        unsigned num_fixed = 0;
         for (unsigned j = 0; j < A_r().m_columns.size(); ++j) {
             auto ci = column_index(j);
-            if (is_fixed_at_bound(ci)) {
-                ++num_fixed;
+            if (is_fixed_at_bound(ci, bounds)) 
                 verbose_stream() << "fixed " << j << "\n";
-            }
         }
-        verbose_stream() << "num fixed " << num_fixed << "\n";
-        if (num_fixed > 0)
+        verbose_stream() << "num fixed " << bounds.size() << "\n";
+        if (!bounds.empty())
             find_feasible_solution();
-        return num_fixed > 0;
+        return !bounds.empty();
     }
 
 
@@ -1617,6 +1621,18 @@ namespace lp {
         return ret;
     }
 
+    /**
+    * \brief ensure there is a column index corresponding to vi
+    * If vi is already a column, just return vi
+    * If vi is for a term, then create a row that uses the term.
+    */
+    var_index lar_solver::ensure_column(var_index vi) {
+        if (lp::tv::is_term(vi))
+            return to_column(vi);
+        else
+            return vi;
+    }
+
 
     void lar_solver::add_row_from_term_no_constraint(const lar_term* term, unsigned term_ext_index) {
         TRACE("dump_terms", print_term(*term, tout) << std::endl;);
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index f2017acd7..9c45fdb80 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -18,29 +18,30 @@
 
   --*/
 #pragma once
-#include "util/vector.h"
-#include <utility>
-#include "util/debug.h"
-#include "util/buffer.h"
+#include <algorithm>
+#include <functional>
+#include <stack>
+#include <string>
 #include <unordered_map>
 #include <unordered_set>
-#include <string>
-#include <algorithm>
-#include <stack>
-#include <functional>
+#include <utility>
+
+#include "math/lp/bound_analyzer_on_row.h"
+#include "math/lp/implied_bound.h"
+#include "math/lp/int_solver.h"
 #include "math/lp/lar_constraints.h"
 #include "math/lp/lar_core_solver.h"
-#include "math/lp/numeric_pair.h"
-#include "math/lp/lp_primal_core_solver.h"
-#include "math/lp/random_updater.h"
-#include "util/stacked_value.h"
-#include "math/lp/stacked_vector.h"
-#include "math/lp/implied_bound.h"
-#include "math/lp/bound_analyzer_on_row.h"
-#include "math/lp/int_solver.h"
-#include "math/lp/nra_solver.h"
-#include "math/lp/lp_types.h"
 #include "math/lp/lp_bound_propagator.h"
+#include "math/lp/lp_primal_core_solver.h"
+#include "math/lp/lp_types.h"
+#include "math/lp/nra_solver.h"
+#include "math/lp/numeric_pair.h"
+#include "math/lp/random_updater.h"
+#include "math/lp/stacked_vector.h"
+#include "util/buffer.h"
+#include "util/debug.h"
+#include "util/stacked_value.h"
+#include "util/vector.h"
 
 namespace lp {
 
@@ -48,10 +49,9 @@ class int_branch;
 class int_solver;
 class lar_solver : public column_namer {
     struct term_hasher {
-        std::size_t operator()(const lar_term &t) const
-        {            
-            using std::size_t;
+        std::size_t operator()(const lar_term& t) const {
             using std::hash;
+            using std::size_t;
             using std::string;
             size_t seed = 0;
             int i = 0;
@@ -66,110 +66,107 @@ class lar_solver : public column_namer {
     };
 
     struct term_comparer {
-        bool operator()(const lar_term &a, const lar_term& b) const
-        {
-            return a == b;            
+        bool operator()(const lar_term& a, const lar_term& b) const {
+            return a == b;
         }
     };
-    
+
     //////////////////// fields //////////////////////////
-    lp_settings                                         m_settings;
-    lp_status                                           m_status = lp_status::UNKNOWN;
-    stacked_value<simplex_strategy_enum>                m_simplex_strategy;
+    lp_settings m_settings;
+    lp_status m_status = lp_status::UNKNOWN;
+    stacked_value<simplex_strategy_enum> m_simplex_strategy;
     // such can be found at the initialization step: u < l
-    stacked_value<int>                                  m_crossed_bounds_column; 
-    lar_core_solver                                     m_mpq_lar_core_solver;
-    int_solver *                                        m_int_solver = nullptr;
-    bool                                                m_need_register_terms = false;
-    var_register                                        m_var_register;
-    var_register                                        m_term_register;
-    stacked_vector<ul_pair>                             m_columns_to_ul_pairs;
-    constraint_set                                      m_constraints;
+    stacked_value<int> m_crossed_bounds_column;
+    lar_core_solver m_mpq_lar_core_solver;
+    int_solver* m_int_solver = nullptr;
+    bool m_need_register_terms = false;
+    var_register m_var_register;
+    var_register m_term_register;
+    stacked_vector<ul_pair> m_columns_to_ul_pairs;
+    constraint_set m_constraints;
     // the set of column indices j such that bounds have changed for j
-    u_set                                               m_columns_with_changed_bounds;
-    u_set                                               m_rows_with_changed_bounds;
-    unsigned_vector                                     m_row_bounds_to_replay;
-    
-    u_set                                               m_basic_columns_with_changed_cost;
+    u_set m_columns_with_changed_bounds;
+    u_set m_rows_with_changed_bounds;
+    unsigned_vector m_row_bounds_to_replay;
+
+    u_set m_basic_columns_with_changed_cost;
     // these are basic columns with the value changed, so the corresponding row in the tableau
     // does not sum to zero anymore
-    u_set                                               m_incorrect_columns;
+    u_set m_incorrect_columns;
     // copy of m_r_solver.inf_heap()
-    unsigned_vector                                     m_inf_index_copy;
-    stacked_value<unsigned>                             m_term_count;
-    vector<lar_term*>                                   m_terms;
-    indexed_vector<mpq>                                 m_column_buffer;
+    unsigned_vector m_inf_index_copy;
+    stacked_value<unsigned> m_term_count;
+    vector<lar_term*> m_terms;
+    indexed_vector<mpq> m_column_buffer;
     std::unordered_map<lar_term, std::pair<mpq, unsigned>, term_hasher, term_comparer>
-    m_normalized_terms_to_columns;
-    vector<impq>                                        m_backup_x;
-    stacked_vector<unsigned>                            m_usage_in_terms;
+        m_normalized_terms_to_columns;
+    vector<impq> m_backup_x;
+    stacked_vector<unsigned> m_usage_in_terms;
     // ((x[j], is_int(j))->j)  for fixed j, used in equalities propagation
     // maps values to integral fixed vars
-    map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>  m_fixed_var_table_int;
+    map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>> m_fixed_var_table_int;
     // maps values to non-integral fixed vars
-    map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>  m_fixed_var_table_real;
+    map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>> m_fixed_var_table_real;
     // end of fields
 
     ////////////////// methods ////////////////////////////////
-    
-    static bool valid_index(unsigned j) { return static_cast<int>(j) >= 0;}
-    const lar_term & get_term(unsigned j) const;
+
+    static bool valid_index(unsigned j) { return static_cast<int>(j) >= 0; }
+    const lar_term& get_term(unsigned j) const;
     bool row_has_a_big_num(unsigned i) const;
     // init region
     bool strategy_is_undecided() const;
     void register_new_ext_var_index(unsigned ext_v, bool is_int);
-    bool term_is_int(const lar_term * t) const;
-    bool term_is_int(const vector<std::pair<mpq, unsigned int>> & coeffs) const;
+    bool term_is_int(const lar_term* t) const;
+    bool term_is_int(const vector<std::pair<mpq, unsigned int>>& coeffs) const;
     void add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int);
     void add_new_var_to_core_fields_for_mpq(bool register_in_basis);
     mpq adjust_bound_for_int(lpvar j, lconstraint_kind&, const mpq&);
 
     // terms
-    bool all_vars_are_registered(const vector<std::pair<mpq, var_index>> & coeffs);
-    var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs);
-    bool term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs);
+    bool all_vars_are_registered(const vector<std::pair<mpq, var_index>>& coeffs);
+    var_index add_term_undecided(const vector<std::pair<mpq, var_index>>& coeffs);
+    bool term_coeffs_are_ok(const vector<std::pair<mpq, var_index>>& coeffs);
     void push_term(lar_term* t);
-    void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index);
+    void add_row_from_term_no_constraint(const lar_term* term, unsigned term_ext_index);
     void add_basic_var_to_core_fields();
-    bool compare_values(impq const& lhs, lconstraint_kind k, const mpq & rhs);
+    bool compare_values(impq const& lhs, lconstraint_kind k, const mpq& rhs);
 
     inline void clear_columns_with_changed_bounds() { m_columns_with_changed_bounds.clear(); }
     inline void increase_by_one_columns_with_changed_bounds() { m_columns_with_changed_bounds.increase_size_by_one(); }
     inline void insert_to_columns_with_changed_bounds(unsigned j) { m_columns_with_changed_bounds.insert(j); }
-    
-    void update_column_type_and_bound_check_on_equal(unsigned j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, unsigned&);
-    void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
-    void update_column_type_and_bound_with_ub(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
-    void update_column_type_and_bound_with_no_ub(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
-    void update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
-    void update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
-    void update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
-    void update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
+
+    void update_column_type_and_bound_check_on_equal(unsigned j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index, unsigned&);
+    void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index);
+    void update_column_type_and_bound_with_ub(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index);
+    void update_column_type_and_bound_with_no_ub(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index);
+    void update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index);
+    void update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index);
+    void update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index);
+    void update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index);
     void register_in_fixed_var_table(unsigned, unsigned&);
     void remove_non_fixed_from_fixed_var_table();
-    constraint_index add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side);
+    constraint_index add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq& right_side);
     inline void set_infeasible_column(unsigned j) {
         set_status(lp_status::INFEASIBLE);
         m_crossed_bounds_column = j;
     }
     constraint_index add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
-                                                                        lconstraint_kind kind, const mpq & right_side);
+                                                                        lconstraint_kind kind, const mpq& right_side);
     unsigned row_of_basic_column(unsigned) const;
     void decide_on_strategy_and_adjust_initial_state();
     void adjust_initial_state();
     void adjust_initial_state_for_tableau_rows();
     bool sizes_are_correct() const;
-    bool implied_bound_is_correctly_explained(implied_bound const & be, const vector<std::pair<mpq, unsigned>> & explanation) const;
-    
+    bool implied_bound_is_correctly_explained(implied_bound const& be, const vector<std::pair<mpq, unsigned>>& explanation) const;
 
     void substitute_basis_var_in_terms_for_row(unsigned i);
-    
+
     template <typename T>
-    unsigned calculate_implied_bounds_for_row(unsigned row_index, lp_bound_propagator<T> & bp) {
-        
-        if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation  || row_has_a_big_num(row_index))
+    unsigned calculate_implied_bounds_for_row(unsigned row_index, lp_bound_propagator<T>& bp) {
+        if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation || row_has_a_big_num(row_index))
             return 0;
-        
+
         return bound_analyzer_on_row<row_strip<mpq>, lp_bound_propagator<T>>::analyze_row(
             A_r().m_rows[row_index],
             null_ci,
@@ -177,29 +174,28 @@ class lar_solver : public column_namer {
             row_index,
             bp);
     }
-    
+
     static void clean_popped_elements_for_heap(unsigned n, lpvar_heap& set);
     static void clean_popped_elements(unsigned n, u_set& set);
-    
-    bool maximize_term_on_tableau(const lar_term & term,
-                                  impq &term_max);
+    bool maximize_term_on_tableau(const lar_term& term,
+                                  impq& term_max);
     bool costs_are_zeros_for_r_solver() const;
     bool reduced_costs_are_zeroes_for_r_solver() const;
-    void set_costs_to_zero(const lar_term & term);
-    void prepare_costs_for_r_solver(const lar_term & term);
-    bool maximize_term_on_corrected_r_solver(lar_term & term, impq &term_max);    
+    void set_costs_to_zero(const lar_term& term);
+    void prepare_costs_for_r_solver(const lar_term& term);
+    bool maximize_term_on_corrected_r_solver(lar_term& term, impq& term_max);
     void pop_core_solver_params();
     void pop_core_solver_params(unsigned k);
     void set_upper_bound_witness(var_index j, constraint_index ci);
     void set_lower_bound_witness(var_index j, constraint_index ci);
-    void substitute_terms_in_linear_expression( const vector<std::pair<mpq, var_index>>& left_side_with_terms,
-                                                vector<std::pair<mpq, var_index>> &left_side) const;
-    
+    void substitute_terms_in_linear_expression(const vector<std::pair<mpq, var_index>>& left_side_with_terms,
+                                               vector<std::pair<mpq, var_index>>& left_side) const;
+
     void detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j);
     bool use_tableau_costs() const;
     bool tableau_with_costs() const;
     bool costs_are_used() const;
-    void change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq> & delta);
+    void change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq>& delta);
     void update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j);
     unsigned num_changed_bounds() const { return m_rows_with_changed_bounds.size(); }
     void insert_row_with_changed_bounds(unsigned rid);
@@ -211,19 +207,19 @@ class lar_solver : public column_namer {
     numeric_pair<mpq> get_basic_var_value_from_row(unsigned i);
     bool all_constrained_variables_are_registered(const vector<std::pair<mpq, var_index>>& left_side);
     bool all_constraints_hold() const;
-    bool constraint_holds(const lar_base_constraint & constr, std::unordered_map<var_index, mpq> & var_map) const;
-    static void register_in_map(std::unordered_map<var_index, mpq> & coeffs, const lar_base_constraint & cn, const mpq & a);
-    static void register_monoid_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j);
-    bool the_left_sides_sum_to_zero(const vector<std::pair<mpq, unsigned>> & evidence) const;
+    bool constraint_holds(const lar_base_constraint& constr, std::unordered_map<var_index, mpq>& var_map) const;
+    static void register_in_map(std::unordered_map<var_index, mpq>& coeffs, const lar_base_constraint& cn, const mpq& a);
+    static void register_monoid_in_map(std::unordered_map<var_index, mpq>& coeffs, const mpq& a, unsigned j);
+    bool the_left_sides_sum_to_zero(const vector<std::pair<mpq, unsigned>>& evidence) const;
     bool explanation_is_correct(explanation&) const;
     bool inf_explanation_is_correct() const;
-    mpq sum_of_right_sides_of_explanation(explanation &) const;
+    mpq sum_of_right_sides_of_explanation(explanation&) const;
     void get_infeasibility_explanation_for_inf_sign(
-        explanation & exp,
-        const vector<std::pair<mpq, unsigned>> & inf_row,
+        explanation& exp,
+        const vector<std::pair<mpq, unsigned>>& inf_row,
         int inf_sign) const;
-    mpq get_left_side_val(const lar_base_constraint &  cns, const std::unordered_map<var_index, mpq> & var_map) const;
-    void fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list);
+    mpq get_left_side_val(const lar_base_constraint& cns, const std::unordered_map<var_index, mpq>& var_map) const;
+    void fill_var_set_for_random_update(unsigned sz, var_index const* vars, vector<unsigned>& column_list);
     bool column_represents_row_in_tableau(unsigned j);
     void make_sure_that_the_bottom_right_elem_not_zero_in_tableau(unsigned i, unsigned j);
     void remove_last_row_and_column_from_tableau(unsigned j);
@@ -235,22 +231,22 @@ class lar_solver : public column_namer {
     void clean_inf_heap_of_r_solver_after_pop();
     inline bool column_value_is_integer(unsigned j) const { return get_column_value(j).is_int(); }
     bool model_is_int_feasible() const;
-    
-    bool bound_is_integer_for_integer_column(unsigned j, const mpq & right_side) const;
-    inline lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; }
+
+    bool bound_is_integer_for_integer_column(unsigned j, const mpq& right_side) const;
+    inline lar_core_solver& get_core_solver() { return m_mpq_lar_core_solver; }
     var_index to_column(unsigned ext_j) const;
     void fix_terms_with_rounded_columns();
     bool remove_from_basis(unsigned);
     lar_term get_term_to_maximize(unsigned ext_j) const;
-    bool sum_first_coords(const lar_term& t, mpq & val) const;
+    bool sum_first_coords(const lar_term& t, mpq& val) const;
     void register_normalized_term(const lar_term&, lpvar);
     void deregister_normalized_term(const lar_term&);
 
     mutable std::unordered_set<impq> m_set_of_different_pairs;
-    mutable std::unordered_set<mpq>  m_set_of_different_singles;
+    mutable std::unordered_set<mpq> m_set_of_different_singles;
     mutable mpq m_delta;
 
-public:
+   public:
     // this function just looks at the status
     bool is_feasible() const;
 
@@ -258,7 +254,6 @@ public:
         return m_fixed_var_table_int;
     }
 
-    
     const map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>& fixed_var_table_real() const {
         return m_fixed_var_table_real;
     }
@@ -267,11 +262,12 @@ public:
         return m_fixed_var_table_real;
     }
 
-    bool find_in_fixed_tables(const rational& mpq, bool is_int, unsigned & j) const {
-        return is_int? fixed_var_table_int().find(mpq, j) : fixed_var_table_real().find(mpq, j);
+    bool find_in_fixed_tables(const rational& mpq, bool is_int, unsigned& j) const {
+        return is_int ? fixed_var_table_int().find(mpq, j) : fixed_var_table_real().find(mpq, j);
     }
-    
-    template <typename T> void remove_non_fixed_from_table(T&);
+
+    template <typename T>
+    void remove_non_fixed_from_table(T&);
 
     unsigned external_to_column_index(unsigned) const;
 
@@ -284,38 +280,38 @@ public:
     inline void set_column_value_test(unsigned j, const impq& v) {
         set_column_value(j, v);
     }
-    
+
     var_index add_named_var(unsigned ext_j, bool is_integer, const std::string&);
 
-    lp_status maximize_term(unsigned j_or_term, impq &term_max);
+    lp_status maximize_term(unsigned j_or_term, impq& term_max);
 
-    inline core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) const { 
-        return core_solver_pretty_printer<lp::mpq, lp::impq>(m_mpq_lar_core_solver.m_r_solver, out); 
+    inline core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) const {
+        return core_solver_pretty_printer<lp::mpq, lp::impq>(m_mpq_lar_core_solver.m_r_solver, out);
     }
 
-    void get_infeasibility_explanation(explanation &) const;
+    void get_infeasibility_explanation(explanation&) const;
 
     inline void backup_x() { m_backup_x = m_mpq_lar_core_solver.m_r_x; }
 
     inline void restore_x() { m_mpq_lar_core_solver.m_r_x = m_backup_x; }
 
     template <typename T>
-    void explain_implied_bound(const implied_bound & ib, lp_bound_propagator<T> & bp) {
+    void explain_implied_bound(const implied_bound& ib, lp_bound_propagator<T>& bp) {
         unsigned i = ib.m_row_or_term_index;
         int bound_sign = (ib.m_is_lower_bound ? 1 : -1);
         int j_sign = (ib.m_coeff_before_j_is_pos ? 1 : -1) * bound_sign;
         unsigned bound_j = ib.m_j;
-        if (tv::is_term(bound_j)) 
+        if (tv::is_term(bound_j))
             bound_j = m_var_register.external_to_local(bound_j);
 
         for (auto const& r : get_row(i)) {
             unsigned j = r.var();
-            if (j == bound_j) 
+            if (j == bound_j)
                 continue;
             mpq const& a = r.coeff();
             int a_sign = is_pos(a) ? 1 : -1;
             int sign = j_sign * a_sign;
-            const ul_pair & ul =  m_columns_to_ul_pairs[j];
+            const ul_pair& ul = m_columns_to_ul_pairs[j];
             auto witness = sign > 0 ? ul.upper_bound_witness() : ul.lower_bound_witness();
             lp_assert(is_valid(witness));
             bp.consume(a, witness);
@@ -329,13 +325,13 @@ public:
     }
 
     // lp_assert(implied_bound_is_correctly_explained(ib, explanation)); }
-    constraint_index mk_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side);
+    constraint_index mk_var_bound(var_index j, lconstraint_kind kind, const mpq& right_side);
     void activate_check_on_equal(constraint_index, var_index&);
     void activate(constraint_index);
-    void random_update(unsigned sz, var_index const * vars);
+    void random_update(unsigned sz, var_index const* vars);
     void mark_rows_for_bound_prop(lpvar j);
     template <typename T>
-    void propagate_bounds_for_touched_rows(lp_bound_propagator<T> & bp) {
+    void propagate_bounds_for_touched_rows(lp_bound_propagator<T>& bp) {
         unsigned num_prop = 0;
         for (unsigned i : m_rows_with_changed_bounds) {
             num_prop += calculate_implied_bounds_for_row(i, bp);
@@ -349,7 +345,7 @@ public:
             bp.clear_for_eq();
             for (unsigned i : m_rows_with_changed_bounds) {
                 unsigned offset_eqs = stats().m_offset_eqs;
-                bp.cheap_eq_tree(i);                
+                bp.cheap_eq_tree(i);
                 if (settings().get_cancel_flag())
                     return;
                 if (stats().m_offset_eqs > offset_eqs)
@@ -360,75 +356,81 @@ public:
     }
 
     template <typename T>
-    void check_missed_propagations(lp_bound_propagator<T> & bp) {
-        for (unsigned i = 0; i < A_r().row_count(); i++) 
-            if (!m_rows_with_changed_bounds.contains(i))                
+    void check_missed_propagations(lp_bound_propagator<T>& bp) {
+        for (unsigned i = 0; i < A_r().row_count(); i++)
+            if (!m_rows_with_changed_bounds.contains(i))
                 if (0 < calculate_implied_bounds_for_row(i, bp)) {
                     verbose_stream() << i << ": " << get_row(i) << "\n";
                 }
     }
 
-    bool is_fixed_at_bound(column_index const& j);
-    bool has_fixed_at_bound();
-    
-    bool is_fixed(column_index const& j) const { return column_is_fixed(j); }    
+    bool is_fixed_at_bound(column_index const& j, vector<std::tuple<explanation, column_index, bool, mpq>>& bounds);
+    bool has_fixed_at_bound(vector<std::tuple<explanation, column_index, bool, mpq>>& bounds);
+
+    bool is_fixed(column_index const& j) const { return column_is_fixed(j); }
     inline column_index to_column_index(unsigned v) const { return column_index(external_to_column_index(v)); }
     bool external_is_used(unsigned) const;
     void pop(unsigned k);
-    bool compare_values(var_index j, lconstraint_kind kind, const mpq & right_side);
-    var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs, unsigned ext_i);
+    bool compare_values(var_index j, lconstraint_kind kind, const mpq& right_side);
+    var_index add_term(const vector<std::pair<mpq, var_index>>& coeffs, unsigned ext_i);
     void register_existing_terms();
-    constraint_index add_var_bound(var_index, lconstraint_kind, const mpq &);
-    constraint_index add_var_bound_check_on_equal(var_index, lconstraint_kind, const mpq &, var_index&);
-    
+    var_index ensure_column(var_index vi);
+    constraint_index add_var_bound(var_index, lconstraint_kind, const mpq&);
+    constraint_index add_var_bound_check_on_equal(var_index, lconstraint_kind, const mpq&, var_index&);
+
     var_index add_var(unsigned ext_j, bool is_integer);
     void set_cut_strategy(unsigned cut_frequency);
     inline unsigned column_count() const { return A_r().column_count(); }
     inline var_index local_to_external(var_index idx) const {
-        return tv::is_term(idx)?
-            m_term_register.local_to_external(idx) : m_var_register.local_to_external(idx);
+        return tv::is_term(idx) ? m_term_register.local_to_external(idx) : m_var_register.local_to_external(idx);
     }
-   
     bool column_corresponds_to_term(unsigned) const;
+    const lar_term& column_to_term(unsigned j) const {
+        SASSERT(column_corresponds_to_term(j));
+        return get_term(column2tv(to_column_index(j)));
+    }
+
     inline unsigned row_count() const { return A_r().row_count(); }
     bool var_is_registered(var_index vj) const;
     void clear_inf_heap() {
-        m_mpq_lar_core_solver.m_r_solver.inf_heap().clear();        
+        m_mpq_lar_core_solver.m_r_solver.inf_heap().clear();
     }
-    inline void remove_column_from_inf_heap(unsigned j) {
+
+    inline void remove_column_from_inf_set(unsigned j) {
         m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_heap(j);
     }
+
+    void pivot(int entering, int leaving) {
+        m_mpq_lar_core_solver.pivot(entering, leaving);
+    }
+
     template <typename ChangeReport>
     void change_basic_columns_dependend_on_a_given_nb_column_report(unsigned j,
-                                                                    const numeric_pair<mpq> & delta,
+                                                                    const numeric_pair<mpq>& delta,
                                                                     const ChangeReport& after) {
-        
-      for (const auto & c : A_r().m_columns[j]) {
-                unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()];
-                if (tableau_with_costs()) {
-                    m_basic_columns_with_changed_cost.insert(bj);
-                }
-                m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, - A_r().get_val(c) * delta);
-                after(bj);
-                TRACE("change_x_del",
-                      tout << "changed basis column " << bj << ", it is " <<
-                      ( m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj)?  "feas":"inf") << std::endl;);
-            }
-   }
+        for (const auto& c : A_r().m_columns[j]) {
+            unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()];
+            if (tableau_with_costs())
+                m_basic_columns_with_changed_cost.insert(bj);
+            m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, -A_r().get_val(c) * delta);
+            after(bj);
+            TRACE("change_x_del",
+                  tout << "changed basis column " << bj << ", it is " << (m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj) ? "feas" : "inf") << std::endl;);
+        }
+    }
 
     template <typename ChangeReport>
     void set_value_for_nbasic_column_report(unsigned j,
-                                            const impq & new_val,
+                                            const impq& new_val,
                                             const ChangeReport& after) {
-
         lp_assert(!is_base(j));
-        auto & x = m_mpq_lar_core_solver.m_r_x[j];        
+        auto& x = m_mpq_lar_core_solver.m_r_x[j];
         auto delta = new_val - x;
         x = new_val;
         after(j);
         change_basic_columns_dependend_on_a_given_nb_column_report(j, delta, after);
     }
-    
+
     template <typename Blocker, typename ChangeReport>
     bool try_to_patch(lpvar j, const mpq& val,
                       const Blocker& is_blocked,
@@ -445,8 +447,8 @@ public:
         impq delta = get_column_value(j) - ival;
         for (auto c : A_r().column(j)) {
             unsigned row_index = c.var();
-            const mpq & a = c.coeff();        
-            unsigned rj = m_mpq_lar_core_solver.m_r_basis[row_index];      
+            const mpq& a = c.coeff();
+            unsigned rj = m_mpq_lar_core_solver.m_r_basis[row_index];
             impq rj_new_val = a * delta + get_column_value(rj);
             // if (column_is_int(rj) && !rj_new_val.is_int())
             //     return false;
@@ -463,65 +465,62 @@ public:
 
     inline bool column_has_lower_bound(unsigned j) const {
         return m_mpq_lar_core_solver.m_r_solver.column_has_lower_bound(j);
-    }    
+    }
 
-    inline
-    constraint_index get_column_upper_bound_witness(unsigned j) const {
+    inline constraint_index get_column_upper_bound_witness(unsigned j) const {
         if (tv::is_term(j)) {
             j = m_var_register.external_to_local(j);
         }
         return m_columns_to_ul_pairs()[j].upper_bound_witness();
     }
 
-    inline
-    const impq& get_upper_bound(column_index j) const {
+    inline const impq& get_upper_bound(column_index j) const {
         return m_mpq_lar_core_solver.m_r_solver.m_upper_bounds[j];
     }
 
-    inline
-    const impq& get_lower_bound(column_index j) const {
+    inline const impq& get_lower_bound(column_index j) const {
         return m_mpq_lar_core_solver.m_r_solver.m_lower_bounds[j];
     }
-    bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const;    
+    bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const;
     bool has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const;
     bool has_value(var_index var, mpq& value) const;
-    bool fetch_normalized_term_column(const lar_term& t, std::pair<mpq, lpvar>& ) const;
+    bool fetch_normalized_term_column(const lar_term& t, std::pair<mpq, lpvar>&) const;
     unsigned map_term_index_to_column_index(unsigned j) const;
     bool column_is_fixed(unsigned j) const;
     bool column_is_free(unsigned j) const;
     unsigned column_to_reported_index(unsigned j) const;
-    lp_settings & settings();
-    lp_settings const & settings() const;
+    lp_settings& settings();
+    lp_settings const& settings() const;
     statistics& stats();
- 
+
     void updt_params(params_ref const& p);
     column_type get_column_type(unsigned j) const { return m_mpq_lar_core_solver.m_column_types()[j]; }
-    const impq & get_lower_bound(unsigned j) const { return m_mpq_lar_core_solver.m_r_lower_bounds()[j]; }
-    const impq & get_upper_bound(unsigned j) const { return m_mpq_lar_core_solver.m_r_upper_bounds()[j]; }
+    const impq& get_lower_bound(unsigned j) const { return m_mpq_lar_core_solver.m_r_lower_bounds()[j]; }
+    const impq& get_upper_bound(unsigned j) const { return m_mpq_lar_core_solver.m_r_upper_bounds()[j]; }
     std::ostream& print_terms(std::ostream& out) const;
-    std::ostream& print_term(lar_term const& term, std::ostream & out) const;
-    static std::ostream& print_term_as_indices(lar_term const& term, std::ostream & out);
-    std::ostream& print_constraint_indices_only(const lar_base_constraint * c, std::ostream & out) const;
-    std::ostream& print_implied_bound(const implied_bound& be, std::ostream & out) const;
+    std::ostream& print_term(lar_term const& term, std::ostream& out) const;
+    static std::ostream& print_term_as_indices(lar_term const& term, std::ostream& out);
+    std::ostream& print_constraint_indices_only(const lar_base_constraint* c, std::ostream& out) const;
+    std::ostream& print_implied_bound(const implied_bound& be, std::ostream& out) const;
     std::ostream& print_values(std::ostream& out) const;
     std::ostream& display(std::ostream& out) const;
 
     bool init_model() const;
     mpq get_value(column_index const& j) const;
     mpq get_tv_value(tv const& t) const;
-    const impq & get_tv_ivalue(tv const& t) const;
-    void get_model(std::unordered_map<var_index, mpq> & variable_values) const;
+    const impq& get_tv_ivalue(tv const& t) const;
+    void get_model(std::unordered_map<var_index, mpq>& variable_values) const;
     void get_rid_of_inf_eps();
-    void get_model_do_not_care_about_diff_vars(std::unordered_map<var_index, mpq> & variable_values) const;
+    void get_model_do_not_care_about_diff_vars(std::unordered_map<var_index, mpq>& variable_values) const;
     std::string get_variable_name(var_index vi) const override;
     void set_variable_name(var_index vi, std::string);
     inline unsigned number_of_vars() const { return m_var_register.size(); }
     inline bool is_base(unsigned j) const { return m_mpq_lar_core_solver.m_r_heading[j] >= 0; }
-    inline const impq & column_lower_bound(unsigned j) const {
+    inline const impq& column_lower_bound(unsigned j) const {
         return m_mpq_lar_core_solver.lower_bound(j);
     }
 
-    inline const impq & column_upper_bound(unsigned j) const {
+    inline const impq& column_upper_bound(unsigned j) const {
         return m_mpq_lar_core_solver.upper_bound(j);
     }
 
@@ -534,9 +533,9 @@ public:
     }
 
     std::pair<constraint_index, constraint_index> add_equality(lpvar j, lpvar k);
-    
-    inline void get_bound_constraint_witnesses_for_column(unsigned j, constraint_index & lc, constraint_index & uc) const {
-        const ul_pair & ul = m_columns_to_ul_pairs[j];
+
+    inline void get_bound_constraint_witnesses_for_column(unsigned j, constraint_index& lc, constraint_index& uc) const {
+        const ul_pair& ul = m_columns_to_ul_pairs[j];
         lc = ul.lower_bound_witness();
         uc = ul.upper_bound_witness();
     }
@@ -553,19 +552,19 @@ public:
     inline tv column2tv(column_index const& c) const {
         return tv::raw(column_to_reported_index(c));
     }
-    
+
     inline std::ostream& print_column_info(unsigned j, std::ostream& out) const {
         m_mpq_lar_core_solver.m_r_solver.print_column_info(j, out);
         if (tv::is_term(j)) {
             print_term_as_indices(get_term(j), out) << "\n";
-            
-        } else if (column_corresponds_to_term(j)) { 
+
+        } else if (column_corresponds_to_term(j)) {
             const lar_term& t = get_term(m_var_register.local_to_external(j));
             print_term_as_indices(t, out) << "\n";
         }
         return out;
     }
-    
+
     void subst_known_terms(lar_term*);
 
     inline std::ostream& print_column_bound_info(unsigned j, std::ostream& out) const {
@@ -576,67 +575,68 @@ public:
 
     inline bool has_inf_int() const {
         for (unsigned j = 0; j < column_count(); j++) {
-            if (column_is_int(j) && ! column_value_is_int(j))
+            if (column_is_int(j) && !column_value_is_int(j))
                 return true;
         }
         return false;
     }
-    inline const vector<lar_term*> & terms() const { return m_terms; }
+    inline const vector<lar_term*>& terms() const { return m_terms; }
     inline lar_term const& term(unsigned i) const { return *m_terms[i]; }
-    inline void set_int_solver(int_solver * int_slv) { m_int_solver = int_slv; }
-    inline int_solver * get_int_solver() { return m_int_solver; }
-    inline const int_solver * get_int_solver() const { return m_int_solver; }
-    inline const lar_term & get_term(tv const& t) const { lp_assert(t.is_term()); return *m_terms[t.id()]; }
-    lp_status find_feasible_solution();   
+    inline void set_int_solver(int_solver* int_slv) { m_int_solver = int_slv; }
+    inline int_solver* get_int_solver() { return m_int_solver; }
+    inline const int_solver* get_int_solver() const { return m_int_solver; }
+    inline const lar_term& get_term(tv const& t) const {
+        lp_assert(t.is_term());
+        return *m_terms[t.id()];
+    }
+    lp_status find_feasible_solution();
     void move_non_basic_columns_to_bounds(bool);
     bool move_non_basic_column_to_bounds(unsigned j, bool);
     inline bool r_basis_has_inf_int() const {
         for (unsigned j : r_basis()) {
-            if (column_is_int(j) && ! column_value_is_int(j))
+            if (column_is_int(j) && !column_value_is_int(j))
                 return true;
         }
         return false;
     }
     void round_to_integer_solution();
-    inline const row_strip<mpq> &  get_row(unsigned i) const { return A_r().m_rows[i]; }
-    inline const row_strip<mpq> &  basic2row(unsigned i) const { return A_r().m_rows[row_of_basic_column(i)]; }
-    inline const column_strip &  get_column(unsigned i) const { return A_r().m_columns[i]; }
+    inline const row_strip<mpq>& get_row(unsigned i) const { return A_r().m_rows[i]; }
+    inline const row_strip<mpq>& basic2row(unsigned i) const { return A_r().m_rows[row_of_basic_column(i)]; }
+    inline const column_strip& get_column(unsigned i) const { return A_r().m_columns[i]; }
     bool row_is_correct(unsigned i) const;
     bool ax_is_correct() const;
-    bool get_equality_and_right_side_for_term_on_current_x(tv const& t, mpq &rs, constraint_index& ci, bool &upper_bound) const;
+    bool get_equality_and_right_side_for_term_on_current_x(tv const& t, mpq& rs, constraint_index& ci, bool& upper_bound) const;
     bool var_is_int(var_index v) const;
-    inline const vector<int> & r_heading() const { return m_mpq_lar_core_solver.m_r_heading; }
-    inline const vector<unsigned> & r_basis() const { return m_mpq_lar_core_solver.r_basis(); }
-    inline const vector<unsigned> & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); }
-    inline bool column_is_real(unsigned j) const { return !column_is_int(j); }	
+    inline const vector<int>& r_heading() const { return m_mpq_lar_core_solver.m_r_heading; }
+    inline const vector<unsigned>& r_basis() const { return m_mpq_lar_core_solver.r_basis(); }
+    inline const vector<unsigned>& r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); }
+    inline bool column_is_real(unsigned j) const { return !column_is_int(j); }
     lp_status get_status() const;
-    bool has_changed_columns() const { return !m_columns_with_changed_bounds.empty();  }
+    bool has_changed_columns() const { return !m_columns_with_changed_bounds.empty(); }
     void set_status(lp_status s);
     lp_status solve();
-    void fill_explanation_from_crossed_bounds_column(explanation & evidence) const;
+    void fill_explanation_from_crossed_bounds_column(explanation& evidence) const;
     bool term_is_used_as_row(unsigned term) const;
     bool tighten_term_bounds_by_delta(tv const& t, const impq&);
     lar_solver();
     void set_track_pivoted_rows(bool v);
-    bool get_track_pivoted_rows() const;    
+    bool get_track_pivoted_rows() const;
     ~lar_solver() override;
     const vector<impq>& r_x() const { return m_mpq_lar_core_solver.m_r_x; }
     bool column_is_int(unsigned j) const;
     inline bool column_value_is_int(unsigned j) const { return m_mpq_lar_core_solver.m_r_x[j].is_int(); }
-    inline static_matrix<mpq, impq> & A_r() { return m_mpq_lar_core_solver.m_r_A; }
-    inline const static_matrix<mpq, impq> & A_r() const { return m_mpq_lar_core_solver.m_r_A; }
+    inline static_matrix<mpq, impq>& A_r() { return m_mpq_lar_core_solver.m_r_A; }
+    inline const static_matrix<mpq, impq>& A_r() const { return m_mpq_lar_core_solver.m_r_A; }
     // columns
     bool column_is_int(column_index const& j) const { return column_is_int((unsigned)j); }
-//    const impq& get_ivalue(column_index const& j) const { return get_column_value(j); }
+    //    const impq& get_ivalue(column_index const& j) const { return get_column_value(j); }
     const impq& get_column_value(column_index const& j) const { return m_mpq_lar_core_solver.m_r_x[j]; }
-    inline
-    var_index external_to_local(unsigned j) const {
+    inline var_index external_to_local(unsigned j) const {
         var_index local_j;
         if (m_var_register.external_is_used(j, local_j) ||
             m_term_register.external_is_used(j, local_j)) {
             return local_j;
-        }
-        else {
+        } else {
             return -1;
         }
     }
@@ -647,6 +647,5 @@ public:
     }
     friend int_solver;
     friend int_branch;
-    
 };
-}
+}  // namespace lp
diff --git a/src/math/lp/lp_bound_propagator.h b/src/math/lp/lp_bound_propagator.h
index dba93398e..a909e8472 100644
--- a/src/math/lp/lp_bound_propagator.h
+++ b/src/math/lp/lp_bound_propagator.h
@@ -4,39 +4,40 @@
   Nikolaj Bjorner (nbjorner)
   Lev Nachmanson   (levnach)
 */
+//clang-format off
 #pragma once
-#include "math/lp/lp_settings.h"
 #include <utility>
+
+#include "math/lp/lp_settings.h"
+#include "util/uint_set.h"
 namespace lp {
 template <typename T>
 class lp_bound_propagator {
-    class edge; // forward definition
+    class edge;  // forward definition
     // vertex represents a column
-    // The set of vertices is organised in a tree.
+    // The set of vertices is organized in a tree.
     // The edges of the tree are rows,
     // Vertices with m_neg set to false grow with the same rate as the root.
     // Vertices with m_neq set to true diminish with the same rate as the roow grows.
     // When two vertices with the same m_neg have the same value of columns
-    // then we have an equality betweet the columns.
-    class vertex {        
-        unsigned           m_column;
-        vector<edge>       m_edges;
-        edge               m_edge_from_parent;
-        unsigned           m_level; // the distance in hops to the root;
-                                    // it is handy to find the common ancestor
-    public:
+    // then we have an equality between the columns.
+    class vertex {
+        unsigned m_column;
+        vector<edge> m_edges;
+        edge m_edge_from_parent;
+        unsigned m_level;  // the distance in hops to the root;
+                           // it is handy to find the common ancestor
+       public:
         vertex() {}
-        vertex(unsigned column) :
-            m_column(column),
-            m_level(0)
-        {}
+        vertex(unsigned column) : m_column(column),
+                                  m_level(0) {}
         unsigned column() const { return m_column; }
         const vertex* parent() const { return m_edge_from_parent.source(); }
         vertex* parent() { return m_edge_from_parent.source(); }
         unsigned level() const { return m_level; }
-        void set_edge_from_parent(edge &e) { m_edge_from_parent = e; }
+        void set_edge_from_parent(edge& e) { m_edge_from_parent = e; }
         const edge& edge_from_parent() const { return m_edge_from_parent; }
-        
+
         void add_child(int row, vertex* child) {
             SASSERT(*this != *child);
             SASSERT(child->parent() == nullptr);
@@ -45,7 +46,7 @@ class lp_bound_propagator {
             child->set_edge_from_parent(e);
             child->m_level = m_level + 1;
         }
-        const vector<edge> & edges() const { return m_edges; }
+        const vector<edge>& edges() const { return m_edges; }
         bool operator==(const vertex& o) const {
             return m_column == o.m_column;
         }
@@ -58,7 +59,8 @@ class lp_bound_propagator {
         vertex* m_source;
         vertex* m_target;
         int m_row;
-    public:
+
+       public:
         edge(vertex* source, vertex* target, int row) : m_source(source), m_target(target), m_row(row) {}
         edge() : m_source(nullptr), m_target(nullptr), m_row(-1) {}
         const vertex* source() const { return m_source; }
@@ -69,57 +71,58 @@ class lp_bound_propagator {
         edge reverse() const { return edge(m_target, m_source, m_row); }
     };
 
-    static int other(int x, int y, int z) { SASSERT(x == z || y == z); return x == z ? y : x; }
-    std::ostream& print_vert(std::ostream & out, const vertex* v) const {
+    static int other(int x, int y, int z) {
+        SASSERT(x == z || y == z);
+        return x == z ? y : x;
+    }
+    std::ostream& print_vert(std::ostream& out, const vertex* v) const {
         out << "(c = " << v->column() << ", parent = {";
         if (v->parent())
             out << "(" << v->parent()->column() << ")";
         else
-            out << "null"; 
-        out <<  "} , lvl = " << v->level();
-        if (m_pol.contains(v->column())) 
-            out << (pol(v) == -1? " -":" +");
+            out << "null";
+        out << "} , lvl = " << v->level();
+        if (m_pol.contains(v->column()))
+            out << (pol(v) == -1 ? " -" : " +");
         else
             out << " not in m_pol";
         out << ')';
         return out;
     }
 
-    hashtable<unsigned, u_hash, u_eq>         m_visited_rows;
-    hashtable<unsigned, u_hash, u_eq>         m_visited_columns;
-    u_map<vertex*>                            m_vertices;
-    vertex*                                   m_root = nullptr;
+    hashtable<unsigned, u_hash, u_eq> m_visited_rows;
+    hashtable<unsigned, u_hash, u_eq> m_visited_columns;
+    u_map<vertex*> m_vertices;
+    vertex* m_root = nullptr;
     // At some point we can find a row with a single vertex non fixed vertex
     // then we can fix the whole tree,
     // by adjusting the vertices offsets, so they become absolute.
     // If the tree is fixed then in addition to checking with the m_vals_to_verts
     // we are going to check with the m_fixed_var_tables.
-    const vertex*                             m_fixed_vertex = nullptr;
-    explanation                               m_fixed_vertex_explanation;
+    const vertex* m_fixed_vertex = nullptr;
+    explanation m_fixed_vertex_explanation;
     // a pair (o, j) belongs to m_vals_to_verts iff x[j] = x[m_root->column()] + o
-    map<mpq, const vertex*, obj_hash<mpq>, default_eq<mpq>>  m_vals_to_verts;
+    map<mpq, const vertex*, obj_hash<mpq>, default_eq<mpq>> m_vals_to_verts;
     // a pair (o, j) belongs to m_vals_to_verts_neg iff -x[j] = x[m_root->column()] + o
-    map<mpq, const vertex*, obj_hash<mpq>, default_eq<mpq>>  m_vals_to_verts_neg;
+    map<mpq, const vertex*, obj_hash<mpq>, default_eq<mpq>> m_vals_to_verts_neg;
     // x[m_root->column()] - m_pol[j].pol()*x[j] == const;
     // to bind polarity and the vertex in the table
-    u_map<int>                                m_pol; 
+    u_map<int> m_pol;
     // if m_pos.contains(j) then  x[j] = x[m_root->column()] + o
-    uint_set                                  m_pos;
-    
+    uint_set m_pos;
+
     // these maps map a column index to the corresponding index in ibounds
-    std::unordered_map<unsigned, unsigned>    m_improved_lower_bounds;
-    std::unordered_map<unsigned, unsigned>    m_improved_upper_bounds;
-    
-    T&                                        m_imp;
-    vector<implied_bound>                     m_ibounds;
+    std::unordered_map<unsigned, unsigned> m_improved_lower_bounds;
+    std::unordered_map<unsigned, unsigned> m_improved_upper_bounds;
 
+    T& m_imp;
+    vector<implied_bound> m_ibounds;
 
+    map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>> m_val2fixed_row;
 
-    map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>  m_val2fixed_row;
-
-    bool is_fixed_row(unsigned r, unsigned & x) {
+    bool is_fixed_row(unsigned r, unsigned& x) {
         x = UINT_MAX;
-        const auto & row = lp().get_row(r);
+        const auto& row = lp().get_row(r);
         for (unsigned k = 0; k < row.size(); k++) {
             const auto& c = row[k];
             if (column_is_fixed(c.var()))
@@ -130,7 +133,7 @@ class lp_bound_propagator {
         }
         return x != UINT_MAX;
     }
-    
+
     void try_add_equation_with_internal_fixed_tables(unsigned r1) {
         SASSERT(m_fixed_vertex);
         unsigned v1, v2;
@@ -154,8 +157,8 @@ class lp_bound_propagator {
         TRACE("eq", print_row(tout, r1); print_row(tout, r2); tout << v1 << " == " << v2 << " = " << val(v1) << "\n");
         add_eq_on_columns(ex, v1, v2, true);
     }
-    
-    void try_add_equation_with_lp_fixed_tables(unsigned row_index, const vertex *v) {
+
+    void try_add_equation_with_lp_fixed_tables(unsigned row_index, const vertex* v) {
         SASSERT(m_fixed_vertex);
         unsigned v_j = v->column();
         unsigned j = null_lpvar;
@@ -163,23 +166,24 @@ class lp_bound_propagator {
             try_add_equation_with_internal_fixed_tables(row_index);
             return;
         }
-       
+
         TRACE("cheap_eq",
-              tout << "v_j = "; lp().print_column_info(v_j, tout) << std::endl;
+              tout << "v_j = ";
+              lp().print_column_info(v_j, tout) << std::endl;
               tout << "v = "; print_vert(tout, v) << std::endl;
-              tout << "found j " << j << std::endl; lp().print_column_info(j, tout)<< std::endl;
+              tout << "found j " << j << std::endl; lp().print_column_info(j, tout) << std::endl;
               tout << "found j = " << j << std::endl;);
         vector<edge> path = connect_in_tree(v, m_fixed_vertex);
         explanation ex = get_explanation_from_path(path);
         ex.add_expl(m_fixed_vertex_explanation);
         explain_fixed_column(j, ex);
-        add_eq_on_columns(ex, j, v_j, true);         
+        add_eq_on_columns(ex, j, v_j, true);
     }
 
-    void try_add_equation_with_val_table(const vertex *v) {
+    void try_add_equation_with_val_table(const vertex* v) {
         SASSERT(m_fixed_vertex);
         unsigned v_j = v->column();
-        const vertex *u = nullptr;
+        const vertex* u = nullptr;
         if (!m_vals_to_verts.find(val(v_j), u)) {
             m_vals_to_verts.insert(val(v_j), v);
             return;
@@ -187,7 +191,7 @@ class lp_bound_propagator {
         unsigned j = u->column();
         if (j == v_j || is_int(j) != is_int(v_j))
             return;
-    
+
         TRACE("cheap_eq", tout << "found j=" << j << " for v=";
               print_vert(tout, v) << "\n in m_vals_to_verts\n";);
         vector<edge> path = connect_in_tree(u, v);
@@ -198,7 +202,7 @@ class lp_bound_propagator {
 
     static bool not_set(unsigned j) { return j == UINT_MAX; }
     static bool is_set(unsigned j) { return j != UINT_MAX; }
-    
+
     void create_root(unsigned row_index) {
         SASSERT(!m_root && !m_fixed_vertex);
         unsigned x, y;
@@ -210,31 +214,30 @@ class lp_bound_propagator {
         }
         TRACE("cheap_eq", print_row(tout, row_index););
         m_root = alloc_v(x);
-        set_polarity(m_root, 1); // keep m_root in the positive table
+        set_polarity(m_root, 1);  // keep m_root in the positive table
         if (not_set(y)) {
             set_fixed_vertex(m_root);
             explain_fixed_in_row(row_index, m_fixed_vertex_explanation);
-        }
-        else {            
-            vertex *v = add_child_with_check(row_index, y, m_root, polarity);
+        } else {
+            vertex* v = add_child_with_check(row_index, y, m_root, polarity);
             if (v)
                 explore_under(v);
         }
         explore_under(m_root);
     }
 
-        void explore_under(vertex * v) {
+    void explore_under(vertex* v) {
         check_for_eq_and_add_to_val_tables(v);
         go_over_vertex_column(v);
     }
 
     // In case of only one non fixed column, and the function returns true,
     // this column would be represened by x.
-    bool is_tree_offset_row(unsigned row_index, unsigned & x, unsigned & y, int & polarity) const {
-        x = y =  UINT_MAX;
+    bool is_tree_offset_row(unsigned row_index, unsigned& x, unsigned& y, int& polarity) const {
+        x = y = UINT_MAX;
         const row_cell<mpq>* x_cell = nullptr;
         const row_cell<mpq>* y_cell = nullptr;
-        const auto & row = lp().get_row(row_index);
+        const auto& row = lp().get_row(row_index);
         for (unsigned k = 0; k < row.size(); k++) {
             const auto& c = row[k];
             if (column_is_fixed(c.var()))
@@ -242,25 +245,21 @@ class lp_bound_propagator {
             if (not_set(x)) {
                 if (c.coeff().is_one() || c.coeff().is_minus_one()) {
                     x = c.var();
-                    x_cell = & c;
-                }
-                else 
+                    x_cell = &c;
+                } else
                     return false;
-            }
-            else if (not_set(y)) {
+            } else if (not_set(y)) {
                 if (c.coeff().is_one() || c.coeff().is_minus_one()) {
                     y = c.var();
-                    y_cell = & c;
-                }
-                else
+                    y_cell = &c;
+                } else
                     return false;
-            }
-            else
+            } else
                 return false;
         }
         if (is_set(x)) {
             if (is_set(y))
-                polarity = x_cell->coeff().is_pos() == y_cell->coeff().is_pos()? -1 : 1;
+                polarity = x_cell->coeff().is_pos() == y_cell->coeff().is_pos() ? -1 : 1;
             else
                 polarity = 1;
             return true;
@@ -268,18 +267,18 @@ class lp_bound_propagator {
         return false;
     }
 
-    void go_over_vertex_column(vertex * v) {
+    void go_over_vertex_column(vertex* v) {
         lpvar j = v->column();
         if (!check_insert(m_visited_columns, j))
             return;
-        
-        for (const auto & c : lp().get_column(j)) {
+
+        for (const auto& c : lp().get_column(j)) {
             unsigned row_index = c.var();
             if (!check_insert(m_visited_rows, row_index))
                 continue;
             vertex* u = get_child_from_row(row_index, v);
-            if (u) 
-                explore_under(u);            
+            if (u)
+                explore_under(u);
         }
     }
 
@@ -295,21 +294,18 @@ class lp_bound_propagator {
         m_pol.reset();
         m_vertices.reset();
     }
-    
+
     struct reset_cheap_eq {
         lp_bound_propagator& p;
-        reset_cheap_eq(lp_bound_propagator& p):p(p) {}
+        reset_cheap_eq(lp_bound_propagator& p) : p(p) {}
         ~reset_cheap_eq() { p.reset_cheap_eq_eh(); }
     };
 
-
-public:
-    
-    lp_bound_propagator(T& imp):
-        m_imp(imp) {}
+   public:
+    lp_bound_propagator(T& imp) : m_imp(imp) {}
 
     const vector<implied_bound>& ibounds() const { return m_ibounds; }
-    
+
     void init() {
         m_improved_upper_bounds.clear();
         m_improved_lower_bounds.clear();
@@ -318,24 +314,24 @@ public:
 
     const lar_solver& lp() const { return m_imp.lp(); }
     lar_solver& lp() { return m_imp.lp(); }
-    
+
     column_type get_column_type(unsigned j) const {
         return m_imp.lp().get_column_type(j);
     }
-    
-    const impq & get_lower_bound(unsigned j) const {
+
+    const impq& get_lower_bound(unsigned j) const {
         return m_imp.lp().get_lower_bound(j);
     }
 
-    const mpq & get_lower_bound_rational(unsigned j) const {
+    const mpq& get_lower_bound_rational(unsigned j) const {
         return m_imp.lp().get_lower_bound(j).x;
     }
-    
-    const impq & get_upper_bound(unsigned j) const {
+
+    const impq& get_upper_bound(unsigned j) const {
         return m_imp.lp().get_upper_bound(j);
     }
 
-    const mpq & get_upper_bound_rational(unsigned j) const {
+    const mpq& get_upper_bound_rational(unsigned j) const {
         return m_imp.lp().get_upper_bound(j).x;
     }
 
@@ -343,40 +339,37 @@ public:
     bool column_is_fixed(lpvar j) const {
         return lp().column_is_fixed(j) && get_lower_bound(j).y.is_zero();
     }
-    
-    void try_add_bound(mpq const& v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
-        j = m_imp.lp().column_to_reported_index(j);    
 
-        lconstraint_kind kind = is_low? GE : LE;
+    void try_add_bound(mpq const& v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
+        j = m_imp.lp().column_to_reported_index(j);
+
+        lconstraint_kind kind = is_low ? GE : LE;
         if (strict)
             kind = static_cast<lconstraint_kind>(kind / 2);
-    
+
         if (!m_imp.bound_is_interesting(j, kind, v))
             return;
-        unsigned k; // index to ibounds
+        unsigned k;  // index to ibounds
         if (is_low) {
             if (try_get_value(m_improved_lower_bounds, j, k)) {
-                auto & found_bound = m_ibounds[k];
+                auto& found_bound = m_ibounds[k];
                 if (v > found_bound.m_bound || (v == found_bound.m_bound && !found_bound.m_strict && strict)) {
                     found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
                     TRACE("try_add_bound", m_imp.lp().print_implied_bound(found_bound, tout););
                 }
-            }
-            else {
+            } else {
                 m_improved_lower_bounds[j] = m_ibounds.size();
                 m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
                 TRACE("try_add_bound", m_imp.lp().print_implied_bound(m_ibounds.back(), tout););
             }
-        }
-        else { // the upper bound case
+        } else {  // the upper bound case
             if (try_get_value(m_improved_upper_bounds, j, k)) {
-                auto & found_bound = m_ibounds[k];
+                auto& found_bound = m_ibounds[k];
                 if (v < found_bound.m_bound || (v == found_bound.m_bound && !found_bound.m_strict && strict)) {
                     found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
                     TRACE("try_add_bound", m_imp.lp().print_implied_bound(found_bound, tout););
                 }
-            }
-            else {
+            } else {
                 m_improved_upper_bounds[j] = m_ibounds.size();
                 m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
                 TRACE("try_add_bound", m_imp.lp().print_implied_bound(m_ibounds.back(), tout););
@@ -388,18 +381,18 @@ public:
         m_imp.consume(a, ci);
     }
 
-    const mpq& val(unsigned j) const { 
-        return  lp().get_column_value(j).x;
+    const mpq& val(unsigned j) const {
+        return lp().get_column_value(j).x;
     }
-    
+
     const mpq& val(const vertex* v) const {
         return val(v->column());
     }
-    
-    bool tree_contains_r(vertex* root, vertex *v) const {
+
+    bool tree_contains_r(vertex* root, vertex* v) const {
         if (*root == *v)
             return true;
-        for (auto e : root->edges()) 
+        for (auto e : root->edges())
             if (tree_contains_r(e.target(), v))
                 return true;
         return false;
@@ -415,7 +408,7 @@ public:
         m_pol.insert(j, p);
     }
 
-    void check_and_set_polarity(vertex* v, int polarity, unsigned row_index, vertex*v_parent) {
+    void check_and_set_polarity(vertex* v, int polarity, unsigned row_index, vertex* v_parent) {
         int prev_pol;
         if (!m_pol.find(v->column(), prev_pol)) {
             set_polarity(v, polarity);
@@ -429,106 +422,107 @@ public:
         m_fixed_vertex_explanation = get_explanation_from_path(path);
         explain_fixed_in_row(row_index, m_fixed_vertex_explanation);
         set_fixed_vertex(v);
-        TRACE("cheap_eq", 
-              tout << "polarity switch: " << polarity << "\nv = "; print_vert(tout , v) << "\nu = ";        tout << "fixed vertex explanation\n";
-              for (auto p : m_fixed_vertex_explanation) 
-                  lp().constraints().display(tout, [this](lpvar j) { return lp().get_variable_name(j);}, p.ci()););
-
+        TRACE("cheap_eq",
+              tout << "polarity switch: " << polarity << "\nv = ";
+              print_vert(tout, v) << "\nu = "; tout << "fixed vertex explanation\n";
+              for (auto p
+                   : m_fixed_vertex_explanation)
+                  lp()
+                      .constraints()
+                      .display(
+                          tout, [this](lpvar j) { return lp().get_variable_name(j); }, p.ci()););
     }
-    
-    bool tree_contains(vertex *v) const {
+
+    bool tree_contains(vertex* v) const {
         if (!m_root)
             return false;
         return tree_contains_r(m_root, v);
     }
-    
-    vertex * alloc_v(unsigned column) {
-        vertex * v = alloc(vertex, column);
+
+    vertex* alloc_v(unsigned column) {
+        vertex* v = alloc(vertex, column);
         m_vertices.insert(column, v);
         SASSERT(!tree_contains(v));
         return v;
     }
 
-
     unsigned column(unsigned row, unsigned index) {
         return lp().get_row(row)[index].var();
     }
 
     bool fixed_phase() const { return m_fixed_vertex; }
-    
-    
+
     // Returns the vertex to start exploration from, or nullptr.
     // It is assumed that parent->column() is present in the row
     vertex* get_child_from_row(unsigned row_index, vertex* parent) {
         TRACE("cheap_eq_det", print_row(tout, row_index););
-        unsigned x, y; int row_polarity;
+        unsigned x, y;
+        int row_polarity;
         if (!is_tree_offset_row(row_index, x, y, row_polarity)) {
-            TRACE("cheap_eq_det", tout << "not an offset row\n"; );
+            TRACE("cheap_eq_det", tout << "not an offset row\n";);
             return nullptr;
         }
-        if (not_set(y)) { // there is only one fixed variable in the row
+        if (not_set(y)) {  // there is only one fixed variable in the row
             if (!fixed_phase()) {
                 set_fixed_vertex(parent);
                 explain_fixed_in_row(row_index, m_fixed_vertex_explanation);
             }
             return nullptr;
         }
-        
+
         SASSERT(is_set(x) && is_set(y));
         unsigned col = other(x, y, parent->column());
         return add_child_with_check(row_index, col, parent, row_polarity);
     }
 
-    vertex * add_child_with_check(unsigned row_index, unsigned col, vertex* parent, int row_polarity) {
+    vertex* add_child_with_check(unsigned row_index, unsigned col, vertex* parent, int row_polarity) {
         vertex* vy;
         if (m_vertices.find(col, vy)) {
             SASSERT(vy != nullptr);
             if (!fixed_phase()) {
-                check_and_set_polarity(vy, pol(parent) * row_polarity, row_index, parent); 
+                check_and_set_polarity(vy, pol(parent) * row_polarity, row_index, parent);
             }
-            return nullptr; // it is not a new vertex
-        } 
+            return nullptr;  // it is not a new vertex
+        }
         vy = alloc_v(col);
         parent->add_child(row_index, vy);
         if (!fixed_phase())
             check_and_set_polarity(vy, row_polarity * pol(parent), row_index, parent);
-        return vy;        
+        return vy;
     }
 
-    bool is_equal(lpvar j, lpvar k) const {        
+    bool is_equal(lpvar j, lpvar k) const {
         return m_imp.is_equal(col_to_imp(j), col_to_imp(k));
     }
 
-    void check_for_eq_and_add_to_val_table(vertex* v,  map<mpq, const vertex*, obj_hash<mpq>, default_eq<mpq>>& table) {
+    void check_for_eq_and_add_to_val_table(vertex* v, map<mpq, const vertex*, obj_hash<mpq>, default_eq<mpq>>& table) {
         TRACE("cheap_eq", tout << "v = "; print_vert(tout, v) << "\n";);
-        const vertex *k; // the other vertex
+        const vertex* k;  // the other vertex
         if (table.find(val(v), k)) {
-            TRACE("cheap_eq", tout << "found k " ; print_vert(tout, k) << "\n";);
+            TRACE("cheap_eq", tout << "found k "; print_vert(tout, k) << "\n";);
             if (k->column() != v->column() &&
                 is_int(k->column()) == is_int(v->column()) &&
                 !is_equal(k->column(), v->column())) {
                 report_eq(k, v);
-            }
-            else {
+            } else {
                 TRACE("cheap_eq", tout << "no report\n";);
             }
-        }
-        else {
+        } else {
             TRACE("cheap_eq", tout << "registered: " << val(v) << " -> { "; print_vert(tout, v) << "} \n";);
             table.insert(val(v), v);
         }
     }
-    
+
     void check_for_eq_and_add_to_val_tables(vertex* v) {
         TRACE("cheap_eq_det", print_vert(tout, v) << "\n";);
         if (!fixed_phase()) {
             if (pol(v->column()) == -1)
                 check_for_eq_and_add_to_val_table(v, m_vals_to_verts_neg);
-            else 
+            else
                 check_for_eq_and_add_to_val_table(v, m_vals_to_verts);
         }
     }
-        
+
     void clear_for_eq() {
         m_visited_rows.reset();
         m_visited_columns.reset();
@@ -542,41 +536,40 @@ public:
 
     std::ostream& print_path(const vector<edge>& path, std::ostream& out) const {
         out << "path = \n";
-        for (const edge& k : path) 
+        for (const edge& k : path)
             print_edge(k, out) << "\n";
         return out;
     }
-    
+
     // we have v_i and v_j, indices of vertices at the same offsets
     void report_eq(const vertex* v_i, const vertex* v_j) {
         SASSERT(v_i != v_j);
         SASSERT(lp().get_column_value(v_i->column()) == lp().get_column_value(v_j->column()));
         TRACE("cheap_eq", tout << v_i->column() << " = " << v_j->column() << "\nu = ";
-              print_vert(tout, v_i) << "\nv = "; print_vert(tout, v_j) <<"\n");
-        
+              print_vert(tout, v_i) << "\nv = "; print_vert(tout, v_j) << "\n");
+
         vector<edge> path = connect_in_tree(v_i, v_j);
         lp::explanation exp = get_explanation_from_path(path);
         add_eq_on_columns(exp, v_i->column(), v_j->column(), false);
-        
     }
 
-    std::ostream& print_expl(std::ostream & out, const explanation& exp) const {
-        for (auto p : exp) 
-            lp().constraints().display(out, [this](lpvar j) { return lp().get_variable_name(j);}, p.ci());
+    std::ostream& print_expl(std::ostream& out, const explanation& exp) const {
+        for (auto p : exp)
+            lp().constraints().display(
+                out, [this](lpvar j) { return lp().get_variable_name(j); }, p.ci());
         return out;
     }
-    
+
     bool add_eq_on_columns(const explanation& exp, lpvar j, lpvar k, bool is_fixed) {
         SASSERT(j != k);
         unsigned je = lp().column_to_reported_index(j);
         unsigned ke = lp().column_to_reported_index(k);
-        TRACE("cheap_eq", 
-              tout << "reporting eq " << j  << ", " << k << "\n";
+        TRACE("cheap_eq",
+              tout << "reporting eq " << j << ", " << k << "\n";
               tout << "reported idx " << je << ", " << ke << "\n";
               print_expl(tout, exp);
-              tout << "theory_vars v" << lp().local_to_external(je) << " == v" << lp().local_to_external(ke) << "\n";
-              );
-        
+              tout << "theory_vars v" << lp().local_to_external(je) << " == v" << lp().local_to_external(ke) << "\n";);
+
         bool added = m_imp.add_eq(je, ke, exp, is_fixed);
         if (added) {
             if (is_fixed)
@@ -600,75 +593,75 @@ public:
     bool is_int(lpvar j) const {
         return lp().column_is_int(j);
     }
-    
+
     explanation get_explanation_from_path(vector<edge>& path) const {
         explanation ex;
-        for (edge &e : path)
+        for (edge& e : path)
             explain_fixed_in_row(e.row(), ex);
         return ex;
     }
 
     void explain_fixed_in_row(unsigned row, explanation& ex) const {
         TRACE("cheap_eq", tout << lp().get_row(row) << std::endl);
-        for (const auto & c : lp().get_row(row)) 
-            if (lp().is_fixed(c.var())) 
+        for (const auto& c : lp().get_row(row))
+            if (lp().is_fixed(c.var()))
                 explain_fixed_column(c.var(), ex);
     }
 
-    void explain_fixed_column(unsigned j, explanation & ex) const {
+    void explain_fixed_column(unsigned j, explanation& ex) const {
         SASSERT(column_is_fixed(j));
-        constraint_index lc, uc;            
+        constraint_index lc, uc;
         lp().get_bound_constraint_witnesses_for_column(j, lc, uc);
         ex.push_back(lc);
-        ex.push_back(uc);        
+        ex.push_back(uc);
     }
-    
+
     vector<edge> connect_in_tree(const vertex* u, const vertex* v) const {
         vector<edge> path;
-        TRACE("cheap_eq_details", tout << "u = " ; print_vert(tout, u); tout << "\nv = ";print_vert(tout, v) << "\n";); 
+        TRACE("cheap_eq_details", tout << "u = "; print_vert(tout, u); tout << "\nv = "; print_vert(tout, v) << "\n";);
         vector<edge> v_branch;
         // equalize the levels
         while (u->level() > v->level()) {
             path.push_back(u->edge_from_parent().reverse());
             u = u->parent();
         }
-        
+
         while (u->level() < v->level()) {
             v_branch.push_back(v->edge_from_parent());
             v = v->parent();
         }
         SASSERT(u->level() == v->level());
-        TRACE("cheap_eq_details", tout << "u = " ; print_vert(tout, u); tout << "\nv = "; print_vert(tout, v) << "\n";); 
+        TRACE("cheap_eq_details", tout << "u = "; print_vert(tout, u); tout << "\nv = "; print_vert(tout, v) << "\n";);
         while (u != v) {
             path.push_back(u->edge_from_parent().reverse());
             v_branch.push_back(v->edge_from_parent());
             u = u->parent();
             v = v->parent();
         }
-        for (unsigned i = v_branch.size(); i--; ) {
-           path.push_back(v_branch[i]);
+        for (unsigned i = v_branch.size(); i--;) {
+            path.push_back(v_branch[i]);
         }
         TRACE("cheap_eq", print_path(path, tout););
         return path;
     }
-    
+
     bool tree_is_correct() const {
         std::unordered_set<int> vs;
         return tree_is_correct(m_root, vs);
     }
-    
+
     bool tree_is_correct(vertex* v, std::unordered_set<int>& visited_verts) const {
         if (fixed_phase())
             return true;
-        if (visited_verts.find(v->column()) != visited_verts.end()) 
+        if (visited_verts.find(v->column()) != visited_verts.end())
             return false;
         visited_verts.insert(v->column());
-        for (auto e : v->edges()) 
+        for (auto e : v->edges())
             if (!tree_is_correct(e.target(), visited_verts))
                 return false;
         return true;
     }
-    std::ostream& print_tree(std::ostream & out, vertex* v) const {
+    std::ostream& print_tree(std::ostream& out, vertex* v) const {
         print_vert(out, v);
         out << "\nchildren :\n";
         for (auto c : v->edges()) {
@@ -683,75 +676,74 @@ public:
         try_add_equation_with_lp_fixed_tables(row_index, v);
         try_add_equation_with_val_table(v);
     }
-    
+
     void handle_fixed_phase(unsigned row_index) {
         if (!fixed_phase())
             return;
         const vertex* v = m_root;
         try_add_equation_with_fixed_tables(row_index, v);
-        for (auto e: v->edges())
+        for (auto e : v->edges())
             try_add_equation_with_fixed_tables(row_index, e.target());
     }
 
-    
     void cheap_eq_tree(unsigned row_index) {
         reset_cheap_eq _reset(*this);
         TRACE("cheap_eq_det", tout << "row_index = " << row_index << "\n";);
-        if (!check_insert(m_visited_rows, row_index)) 
+        if (!check_insert(m_visited_rows, row_index))
             return;
         create_root(row_index);
         if (!m_root)
             return;
-        
-        TRACE("cheap_eq", tout << "tree = "; print_tree(tout, m_root) << "\n";);        
+
+        TRACE("cheap_eq", tout << "tree = "; print_tree(tout, m_root) << "\n";);
         SASSERT(tree_is_correct());
         handle_fixed_phase(row_index);
-        
+
         TRACE("cheap_eq",
               tout << "done for row_index " << row_index << "\n";
               tout << "tree size = " << verts_size(););
     }
 
-    std::ostream& print_row(std::ostream & out, unsigned row_index) const  {
-        unsigned x, y; int polarity;
+    std::ostream& print_row(std::ostream& out, unsigned row_index) const {
+        unsigned x, y;
+        int polarity;
         if (true || !is_tree_offset_row(row_index, x, y, polarity))
             return lp().get_int_solver()->display_row_info(out, row_index);
-        
+
         bool first = true;
-        for (const auto &c: lp().A_r().m_rows[row_index]) {
+        for (const auto& c : lp().A_r().m_rows[row_index]) {
             if (lp().column_is_fixed(c.var()))
                 continue;
             if (c.coeff().is_one()) {
-                    if (!first)
-                        out << "+";
-                }
-                else if (c.coeff().is_minus_one())
-                    out << "-";                             
+                if (!first)
+                    out << "+";
+            } else if (c.coeff().is_minus_one())
+                out << "-";
             out << lp().get_variable_name(c.var()) << " ";
-            first = false;        
+            first = false;
         }
         out << "\n";
         return out;
     }
-    
-    void set_fixed_vertex(vertex *v) {        
+
+    void set_fixed_vertex(vertex* v) {
         TRACE("cheap_eq", if (v) print_vert(tout, v); else tout << "set m_fixed_vertex to nullptr"; tout << "\n";);
         SASSERT(!m_fixed_vertex || v == nullptr);
         m_fixed_vertex = v;
     }
-    
+
     unsigned verts_size() const {
         return subtree_size(m_root);
     }
 
     unsigned subtree_size(vertex* v) const {
-        unsigned r = 1; // 1 for v
+        unsigned r = 1;  // 1 for v
         for (auto e : v->edges())
             r += subtree_size(e.target());
         return r;
     }
-            
-    void delete_tree(vertex * v) {
+
+    void delete_tree(vertex* v) {
         for (auto p : v->edges())
             delete_tree(p.target());
         dealloc(v);
@@ -763,7 +755,6 @@ public:
             return false;
         table.insert(j);
         return true;
-    }        
-
+    }
 };
-}
+}  // namespace lp
diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h
index b2e402c08..470405581 100644
--- a/src/math/lp/lp_primal_core_solver.h
+++ b/src/math/lp/lp_primal_core_solver.h
@@ -17,7 +17,7 @@ Revision History:
 
 
 --*/
-
+// clang-format off
 #pragma once
 #include "math/lp/core_solver_pretty_printer.h"
 #include "math/lp/lp_core_solver_base.h"
@@ -37,123 +37,108 @@ Revision History:
 
 namespace lp {
 
-// This core solver solves (Ax=b, lower_bound_values \leq x \leq
-// upper_bound_values, maximize costs*x ) The right side b is given implicitly
-// by x and the basis
-template <typename T, typename X>
-class lp_primal_core_solver : public lp_core_solver_base<T, X> {
-public:
-  int m_sign_of_entering_delta;
-  vector<T> m_costs_backup;
-  unsigned m_inf_row_index_for_tableau;
-  bool m_bland_mode_tableau;
-  u_set m_left_basis_tableau;
-  unsigned m_bland_mode_threshold;
-  unsigned m_left_basis_repeated;
-  vector<unsigned> m_leaving_candidates;
+    // This core solver solves (Ax=b, lower_bound_values \leq x \leq
+    // upper_bound_values, maximize costs*x ) The right side b is given implicitly
+    // by x and the basis
+    template <typename T, typename X>
+    class lp_primal_core_solver : public lp_core_solver_base<T, X> {
+    public:
+        int m_sign_of_entering_delta;
+        vector<T> m_costs_backup;
+        unsigned m_inf_row_index_for_tableau;
+        bool m_bland_mode_tableau;
+        u_set m_left_basis_tableau;
+        unsigned m_bland_mode_threshold;
+        unsigned m_left_basis_repeated;
+        vector<unsigned> m_leaving_candidates;
 
-  std::list<unsigned> m_non_basis_list;
-  void sort_non_basis();
-  int choose_entering_column_tableau();
+        std::list<unsigned> m_non_basis_list;
+        void sort_non_basis();
+        int choose_entering_column_tableau();
 
-  bool needs_to_grow(unsigned bj) const {
-    lp_assert(!this->column_is_feasible(bj));
-    switch (this->m_column_types[bj]) {
-    case column_type::free_column:
-      return false;
-    case column_type::fixed:
-    case column_type::lower_bound:
-    case column_type::boxed:
-      return this->x_below_low_bound(bj);
-    default:
-      return false;
-    }
-    UNREACHABLE(); // unreachable
-    return false;
-  }
+        bool needs_to_grow(unsigned bj) const {
+            lp_assert(!this->column_is_feasible(bj));
+            switch (this->m_column_types[bj]) {
+            case column_type::free_column:
+                return false;
+            case column_type::fixed:
+            case column_type::lower_bound:
+            case column_type::boxed:
+                return this->x_below_low_bound(bj);
+            default:
+                return false;
+            }
+            UNREACHABLE(); // unreachable
+            return false;
+        }
 
-  int inf_sign_of_column(unsigned bj) const {
-    lp_assert(!this->column_is_feasible(bj));
-    switch (this->m_column_types[bj]) {
-    case column_type::free_column:
-      return 0;
-    case column_type::lower_bound:
-      return 1;
-    case column_type::fixed:
-    case column_type::boxed:
-      return this->x_above_upper_bound(bj) ? -1 : 1;
-    default:
-      return -1;
-    }
-    UNREACHABLE(); // unreachable
-    return 0;
-  }
+        int inf_sign_of_column(unsigned bj) const {
+            lp_assert(!this->column_is_feasible(bj));
+            switch (this->m_column_types[bj]) {
+            case column_type::free_column:
+                return 0;
+            case column_type::lower_bound:
+                return 1;
+            case column_type::fixed:
+            case column_type::boxed:
+                return this->x_above_upper_bound(bj) ? -1 : 1;
+            default:
+                return -1;
+            }
+            UNREACHABLE(); // unreachable
+            return 0;
+        }
 
-  bool monoid_can_decrease(const row_cell<T> &rc) const {
-    unsigned j = rc.var();
-    lp_assert(this->column_is_feasible(j));
-    switch (this->m_column_types[j]) {
-    case column_type::free_column:
-      return true;
-    case column_type::fixed:
-      return false;
-    case column_type::lower_bound:
-      if (is_pos(rc.coeff())) {
-        return this->x_above_lower_bound(j);
-      }
+        bool monoid_can_decrease(const row_cell<T> &rc) const {
+            unsigned j = rc.var();
+            lp_assert(this->column_is_feasible(j));
+            switch (this->m_column_types[j]) {
+            case column_type::free_column:
+                return true;
+            case column_type::fixed:
+                return false;
+            case column_type::lower_bound:
+                return !is_pos(rc.coeff()) || this->x_above_lower_bound(j);
+            case column_type::upper_bound:
+                return is_pos(rc.coeff()) || this->x_below_upper_bound(j);
+            case column_type::boxed:
+                if (is_pos(rc.coeff())) 
+                    return this->x_above_lower_bound(j);
+                return this->x_below_upper_bound(j);
+            default:
+                return false;
+            }
+            UNREACHABLE(); // unreachable
+            return false;
+        }
 
-      return true;
-    case column_type::upper_bound:
-      if (is_pos(rc.coeff())) {
-        return true;
-      }
+        bool monoid_can_increase(const row_cell<T> &rc) const {
+            unsigned j = rc.var();
+            lp_assert(this->column_is_feasible(j));
+            switch (this->m_column_types[j]) {
+            case column_type::free_column:
+                return true;
+            case column_type::fixed:
+                return false;
+            case column_type::lower_bound:
+                if (is_neg(rc.coeff())) 
+                    return this->x_above_lower_bound(j);
+                return true;
+            case column_type::upper_bound:
+                if (is_neg(rc.coeff())) 
+                    return true;
+                return this->x_below_upper_bound(j);
+            case column_type::boxed:
+                if (is_neg(rc.coeff())) 
+                    return this->x_above_lower_bound(j);
+                return this->x_below_upper_bound(j);
+            default:
+                return false;
+            }
+            UNREACHABLE(); // unreachable
+            return false;
+        }
 
-      return this->x_below_upper_bound(j);
-    case column_type::boxed:
-      if (is_pos(rc.coeff())) {
-        return this->x_above_lower_bound(j);
-      }
-
-      return this->x_below_upper_bound(j);
-    default:
-      return false;
-    }
-    UNREACHABLE(); // unreachable
-    return false;
-  }
-
-  bool monoid_can_increase(const row_cell<T> &rc) const {
-    unsigned j = rc.var();
-    lp_assert(this->column_is_feasible(j));
-    switch (this->m_column_types[j]) {
-    case column_type::free_column:
-      return true;
-    case column_type::fixed:
-      return false;
-    case column_type::lower_bound:
-      if (is_neg(rc.coeff())) {
-        return this->x_above_lower_bound(j);
-      }
-
-      return true;
-    case column_type::upper_bound:
-      if (is_neg(rc.coeff())) {
-        return true;
-      }
-
-      return this->x_below_upper_bound(j);
-    case column_type::boxed:
-      if (is_neg(rc.coeff())) {
-        return this->x_above_lower_bound(j);
-      }
-
-      return this->x_below_upper_bound(j);
-    default:
-      return false;
-    }
-    UNREACHABLE(); // unreachable
-    return false;
-  }
 /**
  * Return the number of base non-free variables depending on the column j,
  * different from bj,
@@ -174,503 +159,495 @@ public:
     }
     return r;
   }
-
-  int find_beneficial_entering_in_row_tableau_rows_bland_mode(int i, T &a_ent) {
-    int j = -1;
-    unsigned bj = this->m_basis[i];
-    bool bj_needs_to_grow = needs_to_grow(bj);
-    for (const row_cell<T> &rc : this->m_A.m_rows[i]) {
-      if (rc.var() == bj)
-        continue;
-      if (bj_needs_to_grow) {
-        if (!monoid_can_decrease(rc))
-          continue;
-      } else {
-        if (!monoid_can_increase(rc))
-          continue;
-      }
-      if (rc.var() < static_cast<unsigned>(j)) {
-        j = rc.var();
-        a_ent = rc.coeff();
-      }
-    }
-    if (j == -1) {
-      m_inf_row_index_for_tableau = i;
-    }
-
-    return j;
-  }
-
-  int find_beneficial_entering_tableau_rows(int i, T &a_ent) {
-    if (m_bland_mode_tableau)
-      return find_beneficial_entering_in_row_tableau_rows_bland_mode(i, a_ent);
-    // a short row produces short infeasibility explanation and benefits at
-    // least one pivot operation
-    int choice = -1;
-    int nchoices = 0;
-    unsigned min_non_free_so_far = -1;
-    unsigned best_col_sz = -1;
-    unsigned bj = this->m_basis[i];
-    bool bj_needs_to_grow = needs_to_grow(bj);
-    for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) {
-      const row_cell<T> &rc = this->m_A.m_rows[i][k];
-      unsigned j = rc.var();
-      if (j == bj)
-        continue;
-      if (bj_needs_to_grow) {
-        if (!monoid_can_decrease(rc))
-          continue;
-      } else {
-        if (!monoid_can_increase(rc))
-          continue;
-      }
-      unsigned not_free = get_num_of_not_free_basic_dependent_vars(j, min_non_free_so_far, bj);
-      unsigned col_sz = this->m_A.m_columns[j].size();
-      if (not_free < min_non_free_so_far || (not_free == min_non_free_so_far && col_sz < best_col_sz)) {
-        min_non_free_so_far = not_free;
-        best_col_sz = this->m_A.m_columns[j].size();
-        choice = k;
-        nchoices = 1;
-      } else if (not_free == min_non_free_so_far &&
-                 col_sz == best_col_sz) {
-        if (this->m_settings.random_next(++nchoices) == 0){         
-          choice = k;          
+    // clang-format on
+    int find_beneficial_entering_in_row_tableau_rows_bland_mode(int i, T &a_ent) {
+        int j = -1;
+        unsigned bj = this->m_basis[i];
+        bool bj_needs_to_grow = needs_to_grow(bj);
+        for (const row_cell<T> &rc : this->m_A.m_rows[i]) {
+            if (rc.var() == bj)
+                continue;
+            if (bj_needs_to_grow) {
+                if (!monoid_can_decrease(rc))
+                    continue;
+            } else {
+                if (!monoid_can_increase(rc))
+                    continue;
+            }
+            if (rc.var() < static_cast<unsigned>(j)) {
+                j = rc.var();
+                a_ent = rc.coeff();
+            }
         }
-      }
+        if (j == -1)
+            m_inf_row_index_for_tableau = i;
+        return j;
+    }
+    //clang-format off
+    int find_beneficial_entering_tableau_rows(int i, T &a_ent) {
+        if (m_bland_mode_tableau)
+            return find_beneficial_entering_in_row_tableau_rows_bland_mode(i, a_ent);
+        // a short row produces short infeasibility explanation and benefits at
+        // least one pivot operation
+        int choice = -1;
+        int nchoices = 0;
+        unsigned min_non_free_so_far = -1;
+        unsigned best_col_sz = -1;
+        unsigned bj = this->m_basis[i];
+        bool bj_needs_to_grow = needs_to_grow(bj);
+        for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) {
+            const row_cell<T> &rc = this->m_A.m_rows[i][k];
+            unsigned j = rc.var();
+            if (j == bj)
+                continue;
+            if (bj_needs_to_grow) {
+                if (!monoid_can_decrease(rc))
+                    continue;
+            } else {
+                if (!monoid_can_increase(rc))
+                    continue;
+            }
+            unsigned not_free = get_num_of_not_free_basic_dependent_vars(j, min_non_free_so_far, bj);
+            unsigned col_sz = this->m_A.m_columns[j].size();
+            if (not_free < min_non_free_so_far || (not_free == min_non_free_so_far && col_sz < best_col_sz)) {
+                min_non_free_so_far = not_free;
+                best_col_sz = this->m_A.m_columns[j].size();
+                choice = k;
+                nchoices = 1;
+            } else if (not_free == min_non_free_so_far &&
+                       col_sz == best_col_sz) {
+                if (this->m_settings.random_next(++nchoices) == 0) {
+                    choice = k;
+                }
+            }
+        }
+
+        if (choice == -1) {
+            m_inf_row_index_for_tableau = i;
+            return -1;
+        }
+        const row_cell<T> &rc = this->m_A.m_rows[i][choice];
+        a_ent = rc.coeff();
+        return rc.var();
     }
 
-    if (choice == -1) {
-      m_inf_row_index_for_tableau = i;
-      return -1;
+    bool try_jump_to_another_bound_on_entering(unsigned entering, const X &theta, X &t, bool &unlimited);
+
+    bool try_jump_to_another_bound_on_entering_unlimited(unsigned entering, X &t);
+
+    int find_leaving_and_t_tableau(unsigned entering, X &t);
+
+    void limit_theta(const X &lim, X &theta, bool &unlimited) {
+        if (unlimited) {
+            theta = lim;
+            unlimited = false;
+        } else
+            theta = std::min(lim, theta);
     }
-    const row_cell<T> &rc = this->m_A.m_rows[i][choice];
-    a_ent = rc.coeff();
-    return rc.var();
-  }
 
-  bool try_jump_to_another_bound_on_entering(unsigned entering, const X &theta,
-                                             X &t, bool &unlimited);
-  bool try_jump_to_another_bound_on_entering_unlimited(unsigned entering, X &t);
-  int find_leaving_and_t_tableau(unsigned entering, X &t);
-
-  void limit_theta(const X &lim, X &theta, bool &unlimited) {
-    if (unlimited) {
-      theta = lim;
-      unlimited = false;
-    } else {
-      theta = std::min(lim, theta);
+    void limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound(
+        unsigned j, const T &m, X &theta, bool &unlimited) {
+        lp_assert(m < 0 && this->m_column_types[j] == column_type::upper_bound);
+        limit_inf_on_upper_bound_m_neg(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited);
     }
-  }
 
-  void limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound(
-      unsigned j, const T &m, X &theta, bool &unlimited) {
-    lp_assert(m < 0 && this->m_column_types[j] == column_type::upper_bound);
-    limit_inf_on_upper_bound_m_neg(m, this->m_x[j], this->m_upper_bounds[j],
-                                   theta, unlimited);
-  }
+    void limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound(
+        unsigned j, const T &m, X &theta, bool &unlimited) {
+        lp_assert(m < 0 && this->m_column_types[j] == column_type::lower_bound);
+        limit_inf_on_bound_m_neg(m, this->m_x[j], this->m_lower_bounds[j], theta, unlimited);
+    }
 
-  void limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound(
-      unsigned j, const T &m, X &theta, bool &unlimited) {
-    lp_assert(m < 0 && this->m_column_types[j] == column_type::lower_bound);
-    limit_inf_on_bound_m_neg(m, this->m_x[j], this->m_lower_bounds[j], theta,
-                             unlimited);
-  }
+    void limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound(
+        unsigned j, const T &m, X &theta, bool &unlimited) {
+        lp_assert(m > 0 && this->m_column_types[j] == column_type::lower_bound);
+        limit_inf_on_lower_bound_m_pos(m, this->m_x[j], this->m_lower_bounds[j], theta, unlimited);
+    }
 
-  void limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound(
-      unsigned j, const T &m, X &theta, bool &unlimited) {
-    lp_assert(m > 0 && this->m_column_types[j] == column_type::lower_bound);
-    limit_inf_on_lower_bound_m_pos(m, this->m_x[j], this->m_lower_bounds[j],
-                                   theta, unlimited);
-  }
+    void limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound(
+        unsigned j, const T &m, X &theta, bool &unlimited) {
+        lp_assert(m > 0 && this->m_column_types[j] == column_type::upper_bound);
+        limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited);
+    };
 
-  void limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound(
-      unsigned j, const T &m, X &theta, bool &unlimited) {
-    lp_assert(m > 0 && this->m_column_types[j] == column_type::upper_bound);
-    limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta,
-                             unlimited);
-  };
-
-  void get_bound_on_variable_and_update_leaving_precisely(
-      unsigned j, vector<unsigned> &leavings, T m, X &t,
-      T &abs_of_d_of_leaving);
-
- 
+    void get_bound_on_variable_and_update_leaving_precisely(
+        unsigned j, vector<unsigned> &leavings, T m, X &t,
+        T &abs_of_d_of_leaving);
 
 #ifdef Z3DEBUG
-  void check_Ax_equal_b();
-  void check_the_bounds();
-  void check_bound(unsigned i);
-  void check_correctness();
+    void check_Ax_equal_b();
+    void check_the_bounds();
+    void check_bound(unsigned i);
+    void check_correctness();
 #endif
 
+    void backup_and_normalize_costs();
 
-  void backup_and_normalize_costs();
+    void advance_on_entering_and_leaving_tableau(int entering, int leaving, X &t);
+    void advance_on_entering_equal_leaving_tableau(int entering, X &t);
 
-  void advance_on_entering_and_leaving_tableau(int entering, int leaving, X &t);
-  void advance_on_entering_equal_leaving_tableau(int entering, X &t);
-
-  bool need_to_switch_costs() const {
-    if (this->m_settings.simplex_strategy() ==
-        simplex_strategy_enum::tableau_rows)
-      return false;
-    //        lp_assert(calc_current_x_is_feasible() ==
-    //        current_x_is_feasible());
-    return this->current_x_is_feasible() == this->using_infeas_costs();
-  }
-
-  void advance_on_entering_tableau(int entering);
-
-  void push_forward_offset_in_non_basis(unsigned &offset_in_nb);
-
-  unsigned get_number_of_non_basic_column_to_try_for_enter();
-
-  // returns the number of iterations
-  unsigned solve();
-
-  void find_feasible_solution();
-
-  // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;}
-
-  void one_iteration_tableau();
-
-  // this version assumes that the leaving already has the right value, and does
-  // not update it
-  void update_x_tableau_rows(unsigned entering, unsigned leaving,
-                             const X &delta) {
-    this->add_delta_to_x(entering, delta);
-    for (const auto &c : this->m_A.m_columns[entering]) {
-       if (leaving != this->m_basis[c.var()]) {
-         this->add_delta_to_x_and_track_feasibility(
-         this->m_basis[c.var()], -delta * this->m_A.get_val(c));
-       }
-    }    
-  }
-
-  void update_basis_and_x_tableau_rows(int entering, int leaving, X const &tt) {
-    lp_assert(entering != leaving);
-    update_x_tableau_rows(entering, leaving, tt);
-    this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
-    this->change_basis(entering, leaving);
-  }
-
-  void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving,
-                                                    const X &theta) {
-    update_basis_and_x_tableau_rows(entering, leaving, theta);
-    this->track_column_feasibility(entering);
-  }
-
-  int find_smallest_inf_column() {
-    if (this->inf_heap().empty())
-      return -1;
-    return this->inf_heap().min_value();
-  }
-
-  
-
-  const X &get_val_for_leaving(unsigned j) const {
-    lp_assert(!this->column_is_feasible(j));
-    switch (this->m_column_types[j]) {
-    case column_type::fixed:
-    case column_type::upper_bound:
-      return this->m_upper_bounds[j];
-    case column_type::lower_bound:
-      return this->m_lower_bounds[j];
-      break;
-    case column_type::boxed:
-      if (this->x_above_upper_bound(j))
-        return this->m_upper_bounds[j];
-      else
-        return this->m_lower_bounds[j];
-      break;
-    default:
-      UNREACHABLE();
-      return this->m_lower_bounds[j];
-    }
-  }
-
-  void one_iteration_tableau_rows() {
-    int leaving = find_smallest_inf_column();
-    if (leaving == -1) {
-      this->set_status(lp_status::OPTIMAL);
-      return;
+    void pivot(int entering, int leaving) {
+        this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
+        this->change_basis(entering, leaving);
     }
 
-    SASSERT(this->column_is_base(leaving));
+    bool need_to_switch_costs() const {
+        if (this->m_settings.simplex_strategy() ==
+            simplex_strategy_enum::tableau_rows)
+            return false;
+        //        lp_assert(calc_current_x_is_feasible() ==
+        //        current_x_is_feasible());
+        return this->current_x_is_feasible() == this->using_infeas_costs();
+    }
 
-    if (!m_bland_mode_tableau) {
-      if (m_left_basis_tableau.contains(leaving)) {
-        if (++m_left_basis_repeated > m_bland_mode_threshold) {
-          m_bland_mode_tableau = true;
+    void advance_on_entering_tableau(int entering);
+
+    void push_forward_offset_in_non_basis(unsigned &offset_in_nb);
+
+    unsigned get_number_of_non_basic_column_to_try_for_enter();
+
+    // returns the number of iterations
+    unsigned solve();
+
+    void find_feasible_solution();
+
+    // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;}
+
+    void one_iteration_tableau();
+
+    // this version assumes that the leaving already has the right value, and does
+    // not update it
+    void update_x_tableau_rows(unsigned entering, unsigned leaving,
+                               const X &delta) {
+        this->add_delta_to_x(entering, delta);
+        for (const auto &c : this->m_A.m_columns[entering])
+            if (leaving != this->m_basis[c.var()])
+                this->add_delta_to_x_and_track_feasibility(
+                    this->m_basis[c.var()], -delta * this->m_A.get_val(c));
+    }
+
+    void update_basis_and_x_tableau_rows(int entering, int leaving, X const &tt) {
+        lp_assert(entering != leaving);
+        update_x_tableau_rows(entering, leaving, tt);
+        this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
+        this->change_basis(entering, leaving);
+    }
+
+    void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving,
+                                                      const X &theta) {
+        update_basis_and_x_tableau_rows(entering, leaving, theta);
+        this->track_column_feasibility(entering);
+    }
+
+    int find_smallest_inf_column() {
+        if (this->inf_heap().empty())
+            return -1;
+        return this->inf_heap().min_value();
+    }
+
+    const X &get_val_for_leaving(unsigned j) const {
+        lp_assert(!this->column_is_feasible(j));
+        switch (this->m_column_types[j]) {
+            case column_type::fixed:
+            case column_type::upper_bound:
+                return this->m_upper_bounds[j];
+            case column_type::lower_bound:
+                return this->m_lower_bounds[j];
+                break;
+            case column_type::boxed:
+                if (this->x_above_upper_bound(j))
+                    return this->m_upper_bounds[j];
+                else
+                    return this->m_lower_bounds[j];
+                break;
+            default:
+                UNREACHABLE();
+                return this->m_lower_bounds[j];
         }
-      } else {
-        m_left_basis_tableau.insert(leaving);
-      }
-    }
-    T a_ent;
-    int entering = find_beneficial_entering_tableau_rows(
-        this->m_basis_heading[leaving], a_ent);
-    if (entering == -1) {
-      this->set_status(lp_status::INFEASIBLE);
-      return;
-    }
-    const X &new_val_for_leaving = get_val_for_leaving(leaving);
-    X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
-    this->m_x[leaving] = new_val_for_leaving;
-    this->remove_column_from_inf_heap(leaving);
-    advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta);
-    if (this->current_x_is_feasible())
-      this->set_status(lp_status::OPTIMAL);
-  }
-
-  void decide_on_status_when_cannot_find_entering() {
-    this->set_status(this->current_x_is_feasible() ? lp_status::OPTIMAL
-                                                   : lp_status::INFEASIBLE);
-  }
-
-  void limit_theta_on_basis_column_for_feas_case_m_neg_no_check(
-      unsigned j, const T &m, X &theta, bool &unlimited) {
-    lp_assert(m < 0);
-    limit_theta((this->m_lower_bounds[j] - this->m_x[j]) / m, theta, unlimited);
-    if (theta < zero_of_type<X>())
-      theta = zero_of_type<X>();
-  }
-
-  bool limit_inf_on_bound_m_neg(const T &m, const X &x, const X &bound,
-                                X &theta, bool &unlimited) {
-    // x gets smaller
-    lp_assert(m < 0);
-    if (this->below_bound(x, bound))
-      return false;
-    if (this->above_bound(x, bound)) {
-      limit_theta((bound - x) / m, theta, unlimited);
-    } else {
-      theta = zero_of_type<X>();
-      unlimited = false;
-    }
-    return true;
-  }
-
-  bool limit_inf_on_bound_m_pos(const T &m, const X &x, const X &bound,
-                                X &theta, bool &unlimited) {
-    // x gets larger
-    lp_assert(m > 0);
-    if (this->above_bound(x, bound))
-      return false;
-    if (this->below_bound(x, bound)) {
-      limit_theta((bound - x) / m, theta, unlimited);
-    } else {
-      theta = zero_of_type<X>();
-      unlimited = false;
     }
 
-    return true;
-  }
+    void one_iteration_tableau_rows() {
+        int leaving = find_smallest_inf_column();
+        if (leaving == -1) {
+            this->set_status(lp_status::OPTIMAL);
+            return;
+        }
 
-  void limit_inf_on_lower_bound_m_pos(const T &m, const X &x, const X &bound,
-                                      X &theta, bool &unlimited) {
-    // x gets larger
-    lp_assert(m > 0);
-    if (this->below_bound(x, bound)) {
-      limit_theta((bound - x) / m, theta, unlimited);
+        SASSERT(this->column_is_base(leaving));
+
+        if (!m_bland_mode_tableau) {
+            if (m_left_basis_tableau.contains(leaving)) {
+                if (++m_left_basis_repeated > m_bland_mode_threshold) {
+                    m_bland_mode_tableau = true;
+                }
+            } else {
+                m_left_basis_tableau.insert(leaving);
+            }
+        }
+        T a_ent;
+        int entering = find_beneficial_entering_tableau_rows(
+            this->m_basis_heading[leaving], a_ent);
+        if (entering == -1) {
+            this->set_status(lp_status::INFEASIBLE);
+            return;
+        }
+        const X &new_val_for_leaving = get_val_for_leaving(leaving);
+        X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
+        this->m_x[leaving] = new_val_for_leaving;
+        this->remove_column_from_inf_heap(leaving);
+        advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta);
+        if (this->current_x_is_feasible())
+            this->set_status(lp_status::OPTIMAL);
     }
-  }
 
-  void limit_inf_on_upper_bound_m_neg(const T &m, const X &x, const X &bound,
-                                      X &theta, bool &unlimited) {
-    // x gets smaller
-    lp_assert(m < 0);
-    if (this->above_bound(x, bound)) {
-      limit_theta((bound - x) / m, theta, unlimited);
+    void decide_on_status_when_cannot_find_entering() {
+        this->set_status(this->current_x_is_feasible() ? lp_status::OPTIMAL
+                                                       : lp_status::INFEASIBLE);
     }
-  }
 
-  void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j,
-                                                            const T &m,
-                                                            X &theta,
-                                                            bool &unlimited) {
-    const X &x = this->m_x[j];
-    const X &lbound = this->m_lower_bounds[j];
-
-    if (this->below_bound(x, lbound)) {
-      limit_theta((lbound - x) / m, theta, unlimited);
-    } else {
-      const X &ubound = this->m_upper_bounds[j];
-      if (this->below_bound(x, ubound)) {
-        limit_theta((ubound - x) / m, theta, unlimited);
-      } else if (!this->above_bound(x, ubound)) {
-        theta = zero_of_type<X>();
-        unlimited = false;
-      }
+    void limit_theta_on_basis_column_for_feas_case_m_neg_no_check(
+        unsigned j, const T &m, X &theta, bool &unlimited) {
+        lp_assert(m < 0);
+        limit_theta((this->m_lower_bounds[j] - this->m_x[j]) / m, theta, unlimited);
+        if (theta < zero_of_type<X>())
+            theta = zero_of_type<X>();
     }
-  }
 
-  void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j,
-                                                            const T &m,
-                                                            X &theta,
-                                                            bool &unlimited) {
-    //  lp_assert(m < 0 && this->m_column_type[j] == column_type::boxed);
-    const X &x = this->m_x[j];
-    const X &ubound = this->m_upper_bounds[j];
-    if (this->above_bound(x, ubound)) {
-      limit_theta((ubound - x) / m, theta, unlimited);
-    } else {
-      const X &lbound = this->m_lower_bounds[j];
-      if (this->above_bound(x, lbound)) {
-        limit_theta((lbound - x) / m, theta, unlimited);
-      } else if (!this->below_bound(x, lbound)) {
-        theta = zero_of_type<X>();
-        unlimited = false;
-      }
-    }
-  }
-
-  void limit_theta_on_basis_column_for_feas_case_m_pos_no_check(
-      unsigned j, const T &m, X &theta, bool &unlimited) {
-    lp_assert(m > 0);
-    limit_theta((this->m_upper_bounds[j] - this->m_x[j]) / m, theta, unlimited);
-    if (theta < zero_of_type<X>()) {
-      theta = zero_of_type<X>();
-    }
-  }
-
-  // j is a basic column or the entering, in any case x[j] has to stay feasible.
-  // m is the multiplier. updating t in a way that holds the following
-  // x[j] + t * m >=  this->m_lower_bounds[j]( if m < 0 )
-  // or
-  // x[j] + t * m <= this->m_upper_bounds[j]  ( if m > 0)
-  void limit_theta_on_basis_column(unsigned j, T m, X &theta, bool &unlimited) {
-    switch (this->m_column_types[j]) {
-    case column_type::free_column:
-      break;
-    case column_type::upper_bound:
-      if (this->current_x_is_feasible()) {
-        if (m > 0)
-          limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta,
-                                                                   unlimited);
-      } else { // inside of feasibility_loop
-        if (m > 0)
-          limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound(
-              j, m, theta, unlimited);
-        else
-          limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound(
-              j, m, theta, unlimited);
-      }
-      break;
-    case column_type::lower_bound:
-      if (this->current_x_is_feasible()) {
-        if (m < 0)
-          limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta,
-                                                                   unlimited);
-      } else {
-        if (m < 0)
-          limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound(
-              j, m, theta, unlimited);
-        else
-          limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound(
-              j, m, theta, unlimited);
-      }
-      break;
-      // case fixed:
-      //     if (get_this->current_x_is_feasible()) {
-      //         theta = zero_of_type<X>();
-      //         break;
-      //     }
-      //     if (m < 0)
-      //         limit_theta_on_basis_column_for_inf_case_m_neg_fixed(j, m,
-      //         theta);
-      //     else
-      //         limit_theta_on_basis_column_for_inf_case_m_pos_fixed(j, m,
-      //         theta);
-      //     break;
-    case column_type::fixed:
-    case column_type::boxed:
-      if (this->current_x_is_feasible()) {
-        if (m > 0) {
-          limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta,
-                                                                   unlimited);
+    bool limit_inf_on_bound_m_neg(const T &m, const X &x, const X &bound,
+                                  X &theta, bool &unlimited) {
+        // x gets smaller
+        lp_assert(m < 0);
+        if (this->below_bound(x, bound))
+            return false;
+        if (this->above_bound(x, bound)) {
+            limit_theta((bound - x) / m, theta, unlimited);
         } else {
-          limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta,
-                                                                   unlimited);
+            theta = zero_of_type<X>();
+            unlimited = false;
         }
-      } else {
-        if (m > 0) {
-          limit_theta_on_basis_column_for_inf_case_m_pos_boxed(j, m, theta,
-                                                               unlimited);
+        return true;
+    }
+
+    bool limit_inf_on_bound_m_pos(const T &m, const X &x, const X &bound,
+                                  X &theta, bool &unlimited) {
+        // x gets larger
+        lp_assert(m > 0);
+        if (this->above_bound(x, bound))
+            return false;
+        if (this->below_bound(x, bound)) {
+            limit_theta((bound - x) / m, theta, unlimited);
         } else {
-          limit_theta_on_basis_column_for_inf_case_m_neg_boxed(j, m, theta,
-                                                               unlimited);
+            theta = zero_of_type<X>();
+            unlimited = false;
         }
-      }
 
-      break;
-    default:
-      UNREACHABLE();
+        return true;
     }
-    if (!unlimited && theta < zero_of_type<X>()) {
-      theta = zero_of_type<X>();
+
+    void limit_inf_on_lower_bound_m_pos(const T &m, const X &x, const X &bound,
+                                        X &theta, bool &unlimited) {
+        // x gets larger
+        lp_assert(m > 0);
+        if (this->below_bound(x, bound)) {
+            limit_theta((bound - x) / m, theta, unlimited);
+        }
     }
-  }
 
-  bool column_is_benefitial_for_entering_basis(unsigned j) const;
-  void init_infeasibility_costs();
-  void print_column(unsigned j, std::ostream &out);
-
-  void print_bound_info_and_x(unsigned j, std::ostream &out);
-
-  bool basis_column_is_set_correctly(unsigned j) const {
-    return this->m_A.m_columns[j].size() == 1;
-  }
-
-  bool basis_columns_are_set_correctly() const {
-    for (unsigned j : this->m_basis)
-      if (!basis_column_is_set_correctly(j))
-        return false;
-
-    return this->m_basis_heading.size() == this->m_A.column_count() &&
-           this->m_basis.size() == this->m_A.row_count();
-  }
-
-  void init_run_tableau();
-  void update_x_tableau(unsigned entering, const X &delta);
-  // the delta is between the old and the new cost (old - new)
-  void update_reduced_cost_for_basic_column_cost_change(const T &delta,
-                                                        unsigned j) {
-    lp_assert(this->m_basis_heading[j] >= 0);
-    unsigned i = static_cast<unsigned>(this->m_basis_heading[j]);
-    for (const row_cell<T> &rc : this->m_A.m_rows[i]) {
-      unsigned k = rc.var();
-      if (k == j)
-        continue;
-      this->m_d[k] += delta * rc.coeff();
+    void limit_inf_on_upper_bound_m_neg(const T &m, const X &x, const X &bound,
+                                        X &theta, bool &unlimited) {
+        // x gets smaller
+        lp_assert(m < 0);
+        if (this->above_bound(x, bound)) {
+            limit_theta((bound - x) / m, theta, unlimited);
+        }
     }
-  }
 
-  bool update_basis_and_x_tableau(int entering, int leaving, X const &tt);
-  void init_reduced_costs_tableau();
-  void init_tableau_rows() {
-    m_bland_mode_tableau = false;
-    m_left_basis_tableau.clear();
-    m_left_basis_tableau.resize(this->m_A.column_count());
-    m_left_basis_repeated = 0;
-  }
-  // stage1 constructor
-  lp_primal_core_solver(
-      static_matrix<T, X> &A,
-      vector<X> &b, // the right side vector
-      vector<X> &x, // the number of elements in x needs to be at least as large
-                    // as the number of columns in A
-      vector<unsigned> &basis, vector<unsigned> &nbasis, vector<int> &heading,
-      vector<T> &costs, const vector<column_type> &column_type_array,
-      const vector<X> &lower_bound_values, const vector<X> &upper_bound_values,
-      lp_settings &settings, const column_namer &column_names)
-      : lp_core_solver_base<T, X>(A, // b,
-                                  basis, nbasis, heading, x, costs, settings,
-                                  column_names, column_type_array,
-                                  lower_bound_values, upper_bound_values),
-        m_bland_mode_threshold(1000) {
-    this->set_status(lp_status::UNKNOWN);
-  }
+    void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j,
+                                                              const T &m,
+                                                              X &theta,
+                                                              bool &unlimited) {
+        const X &x = this->m_x[j];
+        const X &lbound = this->m_lower_bounds[j];
 
-  friend core_solver_pretty_printer<T, X>;
+        if (this->below_bound(x, lbound)) {
+            limit_theta((lbound - x) / m, theta, unlimited);
+        } else {
+            const X &ubound = this->m_upper_bounds[j];
+            if (this->below_bound(x, ubound)) {
+                limit_theta((ubound - x) / m, theta, unlimited);
+            } else if (!this->above_bound(x, ubound)) {
+                theta = zero_of_type<X>();
+                unlimited = false;
+            }
+        }
+    }
+
+    void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j,
+                                                              const T &m,
+                                                              X &theta,
+                                                              bool &unlimited) {
+        //  lp_assert(m < 0 && this->m_column_type[j] == column_type::boxed);
+        const X &x = this->m_x[j];
+        const X &ubound = this->m_upper_bounds[j];
+        if (this->above_bound(x, ubound)) {
+            limit_theta((ubound - x) / m, theta, unlimited);
+        } else {
+            const X &lbound = this->m_lower_bounds[j];
+            if (this->above_bound(x, lbound)) {
+                limit_theta((lbound - x) / m, theta, unlimited);
+            } else if (!this->below_bound(x, lbound)) {
+                theta = zero_of_type<X>();
+                unlimited = false;
+            }
+        }
+    }
+
+    void limit_theta_on_basis_column_for_feas_case_m_pos_no_check(
+        unsigned j, const T &m, X &theta, bool &unlimited) {
+        lp_assert(m > 0);
+        limit_theta((this->m_upper_bounds[j] - this->m_x[j]) / m, theta, unlimited);
+        if (theta < zero_of_type<X>()) {
+            theta = zero_of_type<X>();
+        }
+    }
+
+    // j is a basic column or the entering, in any case x[j] has to stay feasible.
+    // m is the multiplier. updating t in a way that holds the following
+    // x[j] + t * m >=  this->m_lower_bounds[j]( if m < 0 )
+    // or
+    // x[j] + t * m <= this->m_upper_bounds[j]  ( if m > 0)
+    void limit_theta_on_basis_column(unsigned j, T m, X &theta, bool &unlimited) {
+        switch (this->m_column_types[j]) {
+            case column_type::free_column:
+                break;
+            case column_type::upper_bound:
+                if (this->current_x_is_feasible()) {
+                    if (m > 0)
+                        limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta,
+                                                                                 unlimited);
+                } else {  // inside of feasibility_loop
+                    if (m > 0)
+                        limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound(
+                            j, m, theta, unlimited);
+                    else
+                        limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound(
+                            j, m, theta, unlimited);
+                }
+                break;
+            case column_type::lower_bound:
+                if (this->current_x_is_feasible()) {
+                    if (m < 0)
+                        limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta,
+                                                                                 unlimited);
+                } else {
+                    if (m < 0)
+                        limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound(
+                            j, m, theta, unlimited);
+                    else
+                        limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound(
+                            j, m, theta, unlimited);
+                }
+                break;
+                // case fixed:
+                //     if (get_this->current_x_is_feasible()) {
+                //         theta = zero_of_type<X>();
+                //         break;
+                //     }
+                //     if (m < 0)
+                //         limit_theta_on_basis_column_for_inf_case_m_neg_fixed(j, m,
+                //         theta);
+                //     else
+                //         limit_theta_on_basis_column_for_inf_case_m_pos_fixed(j, m,
+                //         theta);
+                //     break;
+            case column_type::fixed:
+            case column_type::boxed:
+                if (this->current_x_is_feasible()) {
+                    if (m > 0) {
+                        limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta,
+                                                                                 unlimited);
+                    } else {
+                        limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta,
+                                                                                 unlimited);
+                    }
+                } else {
+                    if (m > 0) {
+                        limit_theta_on_basis_column_for_inf_case_m_pos_boxed(j, m, theta,
+                                                                             unlimited);
+                    } else {
+                        limit_theta_on_basis_column_for_inf_case_m_neg_boxed(j, m, theta,
+                                                                             unlimited);
+                    }
+                }
+
+                break;
+            default:
+                UNREACHABLE();
+        }
+        if (!unlimited && theta < zero_of_type<X>()) {
+            theta = zero_of_type<X>();
+        }
+    }
+
+    bool column_is_benefitial_for_entering_basis(unsigned j) const;
+    void init_infeasibility_costs();
+    void print_column(unsigned j, std::ostream &out);
+
+    void print_bound_info_and_x(unsigned j, std::ostream &out);
+
+    bool basis_column_is_set_correctly(unsigned j) const {
+        return this->m_A.m_columns[j].size() == 1;
+    }
+
+    bool basis_columns_are_set_correctly() const {
+        for (unsigned j : this->m_basis)
+            if (!basis_column_is_set_correctly(j))
+                return false;
+
+        return this->m_basis_heading.size() == this->m_A.column_count() &&
+               this->m_basis.size() == this->m_A.row_count();
+    }
+
+    void init_run_tableau();
+    void update_x_tableau(unsigned entering, const X &delta);
+    // the delta is between the old and the new cost (old - new)
+    void update_reduced_cost_for_basic_column_cost_change(const T &delta,
+                                                          unsigned j) {
+        lp_assert(this->m_basis_heading[j] >= 0);
+        unsigned i = static_cast<unsigned>(this->m_basis_heading[j]);
+        for (const row_cell<T> &rc : this->m_A.m_rows[i]) {
+            unsigned k = rc.var();
+            if (k == j)
+                continue;
+            this->m_d[k] += delta * rc.coeff();
+        }
+    }
+
+    bool update_basis_and_x_tableau(int entering, int leaving, X const &tt);
+    void init_reduced_costs_tableau();
+    void init_tableau_rows() {
+        m_bland_mode_tableau = false;
+        m_left_basis_tableau.clear();
+        m_left_basis_tableau.resize(this->m_A.column_count());
+        m_left_basis_repeated = 0;
+    }
+    // stage1 constructor
+    lp_primal_core_solver(
+        static_matrix<T, X> &A,
+        vector<X> &b,  // the right side vector
+        vector<X> &x,  // the number of elements in x needs to be at least as large
+                       // as the number of columns in A
+        vector<unsigned> &basis, vector<unsigned> &nbasis, vector<int> &heading,
+        vector<T> &costs, const vector<column_type> &column_type_array,
+        const vector<X> &lower_bound_values, const vector<X> &upper_bound_values,
+        lp_settings &settings, const column_namer &column_names)
+        : lp_core_solver_base<T, X>(A,  // b,
+                                    basis, nbasis, heading, x, costs, settings,
+                                    column_names, column_type_array,
+                                    lower_bound_values, upper_bound_values),
+          m_bland_mode_threshold(1000) {
+        this->set_status(lp_status::UNKNOWN);
+    }
+
+    friend core_solver_pretty_printer<T, X>;
 };
-} // namespace lp
+}  // namespace lp
diff --git a/src/math/lp/lp_primal_core_solver_tableau_def.h b/src/math/lp/lp_primal_core_solver_tableau_def.h
index 0c4f6216e..dbcced4ab 100644
--- a/src/math/lp/lp_primal_core_solver_tableau_def.h
+++ b/src/math/lp/lp_primal_core_solver_tableau_def.h
@@ -59,6 +59,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
     }
     unsigned j_nz = this->m_m() + 1; // this number is greater than the max column size
     std::list<unsigned>::iterator entering_iter = m_non_basis_list.end();
+    unsigned n = 0;
     for (auto non_basis_iter = m_non_basis_list.begin(); number_of_benefitial_columns_to_go_over && non_basis_iter != m_non_basis_list.end(); ++non_basis_iter) {
         unsigned j = *non_basis_iter;
         if (!column_is_benefitial_for_entering_basis(j))
@@ -71,8 +72,9 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
             entering_iter = non_basis_iter;
             if (number_of_benefitial_columns_to_go_over)
                 number_of_benefitial_columns_to_go_over--;
+            n = 1;
         }
-        else if (t == j_nz && this->m_settings.random_next() % 2 == 0) {
+        else if (t == j_nz && this->m_settings.random_next(++n) == 0) {
             entering_iter = non_basis_iter;
         }
     }// while (number_of_benefitial_columns_to_go_over && initial_offset_in_non_basis != offset_in_nb);
@@ -166,7 +168,8 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
         }
         this->update_basis_and_x_tableau(entering, leaving, t);
         this->iters_with_no_cost_growing() = 0;
-    } else {
+    }
+    else {
         this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
         this->change_basis(entering, leaving);
     }
diff --git a/src/math/lp/lp_types.h b/src/math/lp/lp_types.h
index 3f9c107c5..e4ee535aa 100644
--- a/src/math/lp/lp_types.h
+++ b/src/math/lp/lp_types.h
@@ -20,8 +20,9 @@ Revision History:
 #pragma once
 
 #include <sstream>
-
-
+// clang-format off
+#include <limits.h>
+#include "util/debug.h"
 namespace nla {
     class core;
 }
diff --git a/src/test/lp/.clang-format b/src/test/lp/.clang-format
new file mode 100644
index 000000000..d70ca6b63
--- /dev/null
+++ b/src/test/lp/.clang-format
@@ -0,0 +1,3 @@
+BasedOnStyle: Google
+IndentWidth: 4
+ColumnLimit: 0
diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp
index 9120d64cf..750eaefb1 100644
--- a/src/test/lp/lp.cpp
+++ b/src/test/lp/lp.cpp
@@ -19,81 +19,81 @@
   --*/
 
 #include <limits>
+
+#include "util/rational.h"
 #ifndef _WINDOWS
 #include <dirent.h>
 #endif
-#include <algorithm>
-#include <string>
-#include <set>
-#include <iostream>
-#include <sys/types.h>
+#include <stdlib.h>
 #include <sys/stat.h>
+#include <sys/types.h>
+
+#include <algorithm>
 #include <cstdlib>
 #include <ctime>
-#include <stdlib.h>
+#include <iostream>
+#include <set>
+#include <string>
 #include <utility>
-#include "math/lp/lp_utils.h"
-#include "test/lp/smt_reader.h"
-#include "test/lp/argument_parser.h"
-#include "test/lp/test_file_reader.h"
-#include "math/lp/indexed_value.h"
-#include "math/lp/lar_solver.h"
-#include "math/lp/numeric_pair.h"
-#include "util/stacked_value.h"
-#include "math/lp/u_set.h"
-#include "util/stopwatch.h"
-#include <cstdlib>
-#include "test/lp/gomory_test.h"
-#include "math/lp/matrix.h"
-#include "math/lp/hnf.h"
-#include "math/lp/general_matrix.h"
-#include "math/lp/lp_bound_propagator.h"
-#include "math/lp/nla_solver.h"
-#include "math/lp/horner.h"
-#include "math/lp/cross_nested.h"
-#include "math/lp/int_cube.h"
-#include "math/lp/emonics.h"
-#include "math/lp/static_matrix.h"
 
-bool my_white_space(const char & a) {
-    return a == ' ' || a == '\t';
-}
-size_t number_of_whites(const std::string & s)  {
+#include "math/lp/cross_nested.h"
+#include "math/lp/emonics.h"
+#include "math/lp/general_matrix.h"
+#include "math/lp/hnf.h"
+#include "math/lp/horner.h"
+#include "math/lp/indexed_value.h"
+#include "math/lp/int_cube.h"
+#include "math/lp/lar_solver.h"
+#include "math/lp/lp_bound_propagator.h"
+#include "math/lp/lp_utils.h"
+#include "math/lp/matrix.h"
+#include "math/lp/nla_solver.h"
+#include "math/lp/numeric_pair.h"
+#include "math/lp/static_matrix.h"
+#include "math/lp/u_set.h"
+#include "test/lp/argument_parser.h"
+#include "test/lp/gomory_test.h"
+#include "test/lp/smt_reader.h"
+#include "test/lp/test_file_reader.h"
+#include "util/stacked_value.h"
+#include "util/stopwatch.h"
+void test_patching();
+bool my_white_space(const char &a) { return a == ' ' || a == '\t'; }
+size_t number_of_whites(const std::string &s) {
     size_t i = 0;
-    for(;i < s.size(); i++)
-        if (!my_white_space(s[i])) return i;
+    for (; i < s.size(); i++)
+        if (!my_white_space(s[i]))
+            return i;
     return i;
 }
-size_t number_of_whites_from_end(const std::string & s)  {
+size_t number_of_whites_from_end(const std::string &s) {
     size_t ret = 0;
-    for(int i = static_cast<int>(s.size()) - 1;i >= 0; i--)
-        if (my_white_space(s[i])) ret++;else break;
-    
+    for (int i = static_cast<int>(s.size()) - 1; i >= 0; i--)
+        if (my_white_space(s[i]))
+            ret++;
+        else
+            break;
+
     return ret;
 }
 
-
 std::string &ltrim(std::string &s) {
     s.erase(0, number_of_whites(s));
     return s;
 }
 
-
-
-
-    // trim from end
+// trim from end
 inline std::string &rtrim(std::string &s) {
-    //       s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun<int, int>(std::isspace))).base(), s.end());
+    //       s.erase(std::find_if(s.rbegin(), s.rend(),
+    //       std::not1(std::ptr_fun<int, int>(std::isspace))).base(), s.end());
     s.erase(s.end() - number_of_whites_from_end(s), s.end());
     return s;
 }
-    // trim from both ends
-inline std::string &trim(std::string &s) {
-    return ltrim(rtrim(s));
-}
+// trim from both ends
+inline std::string &trim(std::string &s) { return ltrim(rtrim(s)); }
 
-
-vector<std::string> string_split(const std::string &source, const char *delimiter, bool keep_empty)  {
+vector<std::string> string_split(const std::string &source,
+                                 const char *delimiter, bool keep_empty) {
     vector<std::string> results;
     size_t prev = 0;
     size_t next = 0;
@@ -118,7 +118,6 @@ vector<std::string> split_and_trim(const std::string &line) {
     return ret;
 }
 
-
 namespace nla {
 void test_horner();
 void test_monics();
@@ -131,7 +130,7 @@ void test_basic_lemma_for_mon_zero_from_factors_to_monomial();
 void test_basic_lemma_for_mon_neutral_from_monomial_to_factors();
 void test_basic_lemma_for_mon_neutral_from_factors_to_monomial();
 
-void test_cn_on_expr(nex_sum *t, cross_nested& cn) {
+void test_cn_on_expr(nex_sum *t, cross_nested &cn) {
     t = to_sum(cn.get_nex_creator().simplify(t));
     TRACE("nla_test", tout << "t=" << *t << '\n';);
     cn.run(t);
@@ -147,35 +146,34 @@ void test_nex_order() {
     r.set_number_of_vars(3);
     for (unsigned j = 0; j < r.get_number_of_vars(); j++)
         r.set_var_weight(j, 10 - j);
-    nex_var* a = r.mk_var(0);
-    nex_var* b = r.mk_var(1);
-    nex_var* c = r.mk_var(2);
+    nex_var *a = r.mk_var(0);
+    nex_var *b = r.mk_var(1);
+    nex_var *c = r.mk_var(2);
     ENSURE(r.gt(a, b));
     ENSURE(r.gt(b, c));
     ENSURE(r.gt(a, c));
-    
 
-    
-    nex* ab = r.mk_mul(a, b);
-    nex* ba = r.mk_mul(b, a);
-    nex* ac = r.mk_mul(a, c);
+    nex *ab = r.mk_mul(a, b);
+    nex *ba = r.mk_mul(b, a);
+    nex *ac = r.mk_mul(a, c);
     ENSURE(r.gt(ab, ac));
     ENSURE(!r.gt(ac, ab));
-    nex* _3ac = r.mk_mul(rational(3), a, c);
-    nex* _2ab = r.mk_mul(rational(2), a, b);
+    nex *_3ac = r.mk_mul(rational(3), a, c);
+    nex *_2ab = r.mk_mul(rational(2), a, b);
     ENSURE(r.gt(ab, _3ac));
     ENSURE(!r.gt(_3ac, ab));
     ENSURE(!r.gt(a, ab));
     ENSURE(r.gt(ab, a));
     ENSURE(r.gt(_2ab, _3ac));
     ENSURE(!r.gt(_3ac, _2ab));
-    nex* _2a = r.mk_mul(rational(2), a);
+    nex *_2a = r.mk_mul(rational(2), a);
     ENSURE(!r.gt(_2a, _2ab));
     ENSURE(r.gt(_2ab, _2a));
     ENSURE(nex_creator::equal(ab, ba));
-    nex_sum * five_a_pl_one = r.mk_sum(r.mk_mul(rational(5), a), r.mk_scalar(rational(1)));
-    nex_mul * poly = r.mk_mul(five_a_pl_one, b);
-    nex * p = r.simplify(poly);
+    nex_sum *five_a_pl_one =
+        r.mk_sum(r.mk_mul(rational(5), a), r.mk_scalar(rational(1)));
+    nex_mul *poly = r.mk_mul(five_a_pl_one, b);
+    nex *p = r.simplify(poly);
     std::cout << "poly = " << *poly << " , p = " << *p << "\n";
 #endif
 }
@@ -184,216 +182,208 @@ void test_simplify() {
 #ifdef Z3DEBUG
     nex_creator r;
     cross_nested cn(
-        [](const nex* n) {
-                           TRACE("nla_cn_test", tout << *n << "\n";);
-                           return false;
-                       } ,
-        [](unsigned) { return false; },
-        []() { return 1; }, // for random
-                    r);
+        [](const nex *n) {
+            TRACE("nla_cn_test", tout << *n << "\n";);
+            return false;
+        },
+        [](unsigned) { return false; }, []() { return 1; },  // for random
+        r);
     enable_trace("nla_cn");
     enable_trace("nla_cn_details");
     // enable_trace("nla_cn_details_");
     enable_trace("nla_test");
-    
+
     r.set_number_of_vars(3);
     for (unsigned j = 0; j < r.get_number_of_vars(); j++)
         r.set_var_weight(j, j);
-    nex_var* a = r.mk_var(0);
-    nex_var* b = r.mk_var(1);
-    nex_var* c = r.mk_var(2);
+    nex_var *a = r.mk_var(0);
+    nex_var *b = r.mk_var(1);
+    nex_var *c = r.mk_var(2);
     auto bc = r.mk_mul(b, c);
     auto a_plus_bc = r.mk_sum(a, bc);
     auto two_a_plus_bc = r.mk_mul(r.mk_scalar(rational(2)), a_plus_bc);
     auto simp_two_a_plus_bc = r.simplify(two_a_plus_bc);
-    TRACE("nla_test", tout << * simp_two_a_plus_bc << "\n";);
+    TRACE("nla_test", tout << *simp_two_a_plus_bc << "\n";);
     ENSURE(nex_creator::equal(simp_two_a_plus_bc, two_a_plus_bc));
     auto simp_a_plus_bc = r.simplify(a_plus_bc);
     ENSURE(to_sum(simp_a_plus_bc)->size() > 1);
 
     auto three_ab = r.mk_mul(r.mk_scalar(rational(3)), a, b);
     auto three_ab_square = r.mk_mul(three_ab, three_ab, three_ab);
-    
+
     TRACE("nla_test", tout << "before simplify " << *three_ab_square << "\n";);
     three_ab_square = to_mul(r.simplify(three_ab_square));
     TRACE("nla_test", tout << *three_ab_square << "\n";);
-    const rational& s = three_ab_square->coeff();
+    const rational &s = three_ab_square->coeff();
     ENSURE(s == rational(27));
     auto m = r.mk_mul(a, a);
     TRACE("nla_test_", tout << "m = " << *m << "\n";);
     /*
-    auto n = r.mk_mul(b, b, b, b, b, b, b);
-    n->add_child_in_power(b, 7);
-    n->add_child(r.mk_scalar(rational(3)));
-    n->add_child_in_power(r.mk_scalar(rational(2)), 2);
-    n->add_child(r.mk_scalar(rational(1)));
-    TRACE("nla_test_", tout << "n = " << *n << "\n";); 
-    m->add_child_in_power(n, 3);
-    n->add_child_in_power(r.mk_scalar(rational(1, 3)), 2);
-    TRACE("nla_test_", tout << "m = " << *m << "\n";); 
-    
-    nex_sum * e = r.mk_sum(a, r.mk_sum(b, m));
-    TRACE("nla_test", tout << "before simplify e = " << *e << "\n";);
-    e = to_sum(r.simplify(e));
-    TRACE("nla_test", tout << "simplified e = " << *e << "\n";);
-    ENSURE(e->children().size() > 2);
-    nex_sum * e_m = r.mk_sum();
-    for (const nex* ex: to_sum(e)->children()) {
-        nex* ce = r.mk_mul(r.clone(ex), r.mk_scalar(rational(3)));        
-        TRACE("nla_test", tout << "before simpl ce = " << *ce << "\n";);        
-        ce = r.simplify(ce);
-        TRACE("nla_test", tout << "simplified ce = " << *ce << "\n";);        
-        e_m->add_child(ce);
-    }
-    e->add_child(e_m);    
-    TRACE("nla_test", tout << "before simplify sum e = " << *e << "\n";);
-    e = to_sum(r.simplify(e));
-    TRACE("nla_test", tout << "simplified sum e = " << *e << "\n";);
+  auto n = r.mk_mul(b, b, b, b, b, b, b);
+  n->add_child_in_power(b, 7);
+  n->add_child(r.mk_scalar(rational(3)));
+  n->add_child_in_power(r.mk_scalar(rational(2)), 2);
+  n->add_child(r.mk_scalar(rational(1)));
+  TRACE("nla_test_", tout << "n = " << *n << "\n";);
+  m->add_child_in_power(n, 3);
+  n->add_child_in_power(r.mk_scalar(rational(1, 3)), 2);
+  TRACE("nla_test_", tout << "m = " << *m << "\n";);
 
-    nex * pr = r.mk_mul(a, b, b);
-    TRACE("nla_test", tout << "before simplify pr = " << *pr << "\n";);
-    r.simplify(pr);
-    TRACE("nla_test", tout << "simplified sum e = " << *pr << "\n";);
-    */
+  nex_sum * e = r.mk_sum(a, r.mk_sum(b, m));
+  TRACE("nla_test", tout << "before simplify e = " << *e << "\n";);
+  e = to_sum(r.simplify(e));
+  TRACE("nla_test", tout << "simplified e = " << *e << "\n";);
+  ENSURE(e->children().size() > 2);
+  nex_sum * e_m = r.mk_sum();
+  for (const nex* ex: to_sum(e)->children()) {
+      nex* ce = r.mk_mul(r.clone(ex), r.mk_scalar(rational(3)));
+      TRACE("nla_test", tout << "before simpl ce = " << *ce << "\n";);
+      ce = r.simplify(ce);
+      TRACE("nla_test", tout << "simplified ce = " << *ce << "\n";);
+      e_m->add_child(ce);
+  }
+  e->add_child(e_m);
+  TRACE("nla_test", tout << "before simplify sum e = " << *e << "\n";);
+  e = to_sum(r.simplify(e));
+  TRACE("nla_test", tout << "simplified sum e = " << *e << "\n";);
+
+  nex * pr = r.mk_mul(a, b, b);
+  TRACE("nla_test", tout << "before simplify pr = " << *pr << "\n";);
+  r.simplify(pr);
+  TRACE("nla_test", tout << "simplified sum e = " << *pr << "\n";);
+  */
 #endif
 }
 
 void test_cn_shorter() {
-//     nex_sum *clone;
-//     nex_creator cr;
-//     cross_nested cn(
-//         [](const nex* n) {
-//             TRACE("nla_test", tout <<"cn form = " <<  *n << "\n";
-                  
-// );
-//             return false;
-//         } ,
-//         [](unsigned) { return false; },
-//         []{ return 1; }, cr);
-//     enable_trace("nla_test");
-//     enable_trace("nla_cn");
-//     enable_trace("nla_cn_test");
-//     enable_trace("nla_cn_details");
-//     //    enable_trace("nla_cn_details_");
-//     enable_trace("nla_test_details");
-//     cr.set_number_of_vars(20);
-//     for (unsigned j = 0; j < cr.get_number_of_vars(); j++)
-//         cr.set_var_weight(j,j);
-        
-//     nex_var* a = cr.mk_var(0);
-//     nex_var* b = cr.mk_var(1);
-//     nex_var* c = cr.mk_var(2);
-//     nex_var* d = cr.mk_var(3);
-//     nex_var* e = cr.mk_var(4);
-//     nex_var* g = cr.mk_var(6);
+    //     nex_sum *clone;
+    //     nex_creator cr;
+    //     cross_nested cn(
+    //         [](const nex* n) {
+    //             TRACE("nla_test", tout <<"cn form = " <<  *n << "\n";
 
-//     nex* min_1 = cr.mk_scalar(rational(-1));
-//     // test_cn_on_expr(min_1*c*e + min_1*b*d + min_1*a*b + a*c);
-//     nex_mul* bcg = cr.mk_mul(b, c, g);
-//     /*
-//     bcg->add_child(min_1);
-//     nex* abcd = cr.mk_mul(a, b, c, d);
-//     nex* eae = cr.mk_mul(e, a, e);
-//     nex* three_eac = cr.mk_mul(e, a, c); to_mul(three_eac)->coeff() = rational(3);
-//     nex* _6aad = cr.mk_mul(cr.mk_scalar(rational(6)), a, a, d);
-//     clone = to_sum(cr.clone(cr.mk_sum(_6aad, abcd, eae, three_eac)));
-//     clone = to_sum(cr.simplify(clone));
-//     TRACE("nla_test", tout << "clone = " << *clone << "\n";);
-//     //    test_cn_on_expr(cr.mk_sum(aad,  abcd, aaccd, add, eae, eac, ed), cn);
-//     test_cn_on_expr(clone, cn);
-//     */
+    // );
+    //             return false;
+    //         } ,
+    //         [](unsigned) { return false; },
+    //         []{ return 1; }, cr);
+    //     enable_trace("nla_test");
+    //     enable_trace("nla_cn");
+    //     enable_trace("nla_cn_test");
+    //     enable_trace("nla_cn_details");
+    //     //    enable_trace("nla_cn_details_");
+    //     enable_trace("nla_test_details");
+    //     cr.set_number_of_vars(20);
+    //     for (unsigned j = 0; j < cr.get_number_of_vars(); j++)
+    //         cr.set_var_weight(j,j);
+
+    //     nex_var* a = cr.mk_var(0);
+    //     nex_var* b = cr.mk_var(1);
+    //     nex_var* c = cr.mk_var(2);
+    //     nex_var* d = cr.mk_var(3);
+    //     nex_var* e = cr.mk_var(4);
+    //     nex_var* g = cr.mk_var(6);
+
+    //     nex* min_1 = cr.mk_scalar(rational(-1));
+    //     // test_cn_on_expr(min_1*c*e + min_1*b*d + min_1*a*b + a*c);
+    //     nex_mul* bcg = cr.mk_mul(b, c, g);
+    //     /*
+    //     bcg->add_child(min_1);
+    //     nex* abcd = cr.mk_mul(a, b, c, d);
+    //     nex* eae = cr.mk_mul(e, a, e);
+    //     nex* three_eac = cr.mk_mul(e, a, c); to_mul(three_eac)->coeff() =
+    //     rational(3); nex* _6aad = cr.mk_mul(cr.mk_scalar(rational(6)), a, a,
+    //     d); clone = to_sum(cr.clone(cr.mk_sum(_6aad, abcd, eae, three_eac)));
+    //     clone = to_sum(cr.simplify(clone));
+    //     TRACE("nla_test", tout << "clone = " << *clone << "\n";);
+    //     //    test_cn_on_expr(cr.mk_sum(aad,  abcd, aaccd, add, eae, eac, ed),
+    //     cn); test_cn_on_expr(clone, cn);
+    //     */
 }
 
 void test_cn() {
-// #ifdef Z3DEBUG
-//     test_cn_shorter();
-//     nex_creator cr;
-//     cross_nested cn(
-//         [](const nex* n) {
-//             TRACE("nla_test", tout <<"cn form = " <<  *n << "\n";);
-//             return false;
-//         } ,
-//         [](unsigned) { return false; },
-//         []{ return 1; }, cr);
-//     enable_trace("nla_test");
-//     enable_trace("nla_cn_test");
-//     //    enable_trace("nla_cn");
-//     //   enable_trace("nla_test_details");
-//     cr.set_number_of_vars(20);
-//     for (unsigned j = 0; j < cr.get_number_of_vars(); j++)
-//         cr.set_var_weight(j, j);
-        
-//     nex_var* a = cr.mk_var(0);
-//     nex_var* b = cr.mk_var(1);
-//     nex_var* c = cr.mk_var(2);
-//     nex_var* d = cr.mk_var(3);
-//     nex_var* e = cr.mk_var(4);
-//     nex_var* g = cr.mk_var(6);
-//     nex_sum * a_p_ae_sq = cr.mk_sum(a, cr.mk_mul(a, e, e));
-//     a_p_ae_sq = to_sum(cr.simplify(a_p_ae_sq));
-//     test_cn_on_expr(a_p_ae_sq, cn);
+    // #ifdef Z3DEBUG
+    //     test_cn_shorter();
+    //     nex_creator cr;
+    //     cross_nested cn(
+    //         [](const nex* n) {
+    //             TRACE("nla_test", tout <<"cn form = " <<  *n << "\n";);
+    //             return false;
+    //         } ,
+    //         [](unsigned) { return false; },
+    //         []{ return 1; }, cr);
+    //     enable_trace("nla_test");
+    //     enable_trace("nla_cn_test");
+    //     //    enable_trace("nla_cn");
+    //     //   enable_trace("nla_test_details");
+    //     cr.set_number_of_vars(20);
+    //     for (unsigned j = 0; j < cr.get_number_of_vars(); j++)
+    //         cr.set_var_weight(j, j);
 
-//     nex* min_1 = cr.mk_scalar(rational(-1));
-//     // test_cn_on_expr(min_1*c*e + min_1*b*d + min_1*a*b + a*c);
-//     nex* bcd = cr.mk_mul(b, c, d);
-//     nex_mul* bcg = cr.mk_mul(b, c, g);
-//     /*
-//     bcg->add_child(min_1);
-//     nex_sum* t = cr.mk_sum(bcd, bcg);
-//     test_cn_on_expr(t, cn);
-//     nex* abd = cr.mk_mul(a, b, d);
-//     nex* abc = cr.mk_mul(a, b, c);
-//     nex* abcd = cr.mk_mul(a, b, c, d);
-//     nex* aaccd = cr.mk_mul(a, a, c, c, d);
-//     nex* add = cr.mk_mul(a, d, d);
-//     nex* eae = cr.mk_mul(e, a, e);
-//     nex* eac = cr.mk_mul(e, a, c);
-//     nex* ed = cr.mk_mul(e, d);
-//     nex* cbd = cr.mk_mul(c, b, d);
-//     nex* acd = cr.mk_mul(a, c, d);
-        
-//     nex* _6aad = cr.mk_mul(cr.mk_scalar(rational(6)), a, a, d);
-//     nex * clone = cr.clone(cr.mk_sum(_6aad, abcd, aaccd, add, eae, eac, ed));
-//     clone = cr.simplify(clone);
-//     ENSURE(cr.is_simplified(clone));
-//     TRACE("nla_test", tout << "clone = " << *clone << "\n";);
-//     //    test_cn_on_expr(cr.mk_sum(aad,  abcd, aaccd, add, eae, eac, ed), cn);
-//     test_cn_on_expr(to_sum(clone), cn);
-//     TRACE("nla_test", tout << "done\n";);
-//     test_cn_on_expr(cr.mk_sum(abd,  abc, cbd, acd), cn);
-//     TRACE("nla_test", tout << "done\n";);*/
-// #endif
-//     // test_cn_on_expr(a*b*b*d*d + a*b*b*c*d + c*b*b*d);
-//     // TRACE("nla_test", tout << "done\n";);
-//     // test_cn_on_expr(a*b*d + a*b*c + c*b*d);
+    //     nex_var* a = cr.mk_var(0);
+    //     nex_var* b = cr.mk_var(1);
+    //     nex_var* c = cr.mk_var(2);
+    //     nex_var* d = cr.mk_var(3);
+    //     nex_var* e = cr.mk_var(4);
+    //     nex_var* g = cr.mk_var(6);
+    //     nex_sum * a_p_ae_sq = cr.mk_sum(a, cr.mk_mul(a, e, e));
+    //     a_p_ae_sq = to_sum(cr.simplify(a_p_ae_sq));
+    //     test_cn_on_expr(a_p_ae_sq, cn);
+
+    //     nex* min_1 = cr.mk_scalar(rational(-1));
+    //     // test_cn_on_expr(min_1*c*e + min_1*b*d + min_1*a*b + a*c);
+    //     nex* bcd = cr.mk_mul(b, c, d);
+    //     nex_mul* bcg = cr.mk_mul(b, c, g);
+    //     /*
+    //     bcg->add_child(min_1);
+    //     nex_sum* t = cr.mk_sum(bcd, bcg);
+    //     test_cn_on_expr(t, cn);
+    //     nex* abd = cr.mk_mul(a, b, d);
+    //     nex* abc = cr.mk_mul(a, b, c);
+    //     nex* abcd = cr.mk_mul(a, b, c, d);
+    //     nex* aaccd = cr.mk_mul(a, a, c, c, d);
+    //     nex* add = cr.mk_mul(a, d, d);
+    //     nex* eae = cr.mk_mul(e, a, e);
+    //     nex* eac = cr.mk_mul(e, a, c);
+    //     nex* ed = cr.mk_mul(e, d);
+    //     nex* cbd = cr.mk_mul(c, b, d);
+    //     nex* acd = cr.mk_mul(a, c, d);
+
+    //     nex* _6aad = cr.mk_mul(cr.mk_scalar(rational(6)), a, a, d);
+    //     nex * clone = cr.clone(cr.mk_sum(_6aad, abcd, aaccd, add, eae, eac,
+    //     ed)); clone = cr.simplify(clone); ENSURE(cr.is_simplified(clone));
+    //     TRACE("nla_test", tout << "clone = " << *clone << "\n";);
+    //     //    test_cn_on_expr(cr.mk_sum(aad,  abcd, aaccd, add, eae, eac, ed),
+    //     cn); test_cn_on_expr(to_sum(clone), cn); TRACE("nla_test", tout <<
+    //     "done\n";); test_cn_on_expr(cr.mk_sum(abd,  abc, cbd, acd), cn);
+    //     TRACE("nla_test", tout << "done\n";);*/
+    // #endif
+    //     // test_cn_on_expr(a*b*b*d*d + a*b*b*c*d + c*b*b*d);
+    //     // TRACE("nla_test", tout << "done\n";);
+    //     // test_cn_on_expr(a*b*d + a*b*c + c*b*d);
 }
 
-} // end of namespace nla
+}  // end of namespace nla
 
 namespace lp {
 unsigned seed = 1;
 
-
 random_gen g_rand;
-static unsigned my_random() {
-    return g_rand();
-}
-struct simple_column_namer:public column_namer
-{
+static unsigned my_random() { return g_rand(); }
+struct simple_column_namer : public column_namer {
     std::string get_variable_name(unsigned j) const override {
-        return std::string("x") + T_to_string(j); 
+        return std::string("x") + T_to_string(j);
     }
 };
 
-
-
-vector<int> allocate_basis_heading(unsigned count) { // the rest of initialization will be handled by lu_QR
+vector<int> allocate_basis_heading(
+    unsigned count) {  // the rest of initialization will be handled by lu_QR
     vector<int> basis_heading(count, -1);
     return basis_heading;
 }
 
-
-void init_basic_part_of_basis_heading(vector<unsigned> & basis, vector<int> & basis_heading) {
+void init_basic_part_of_basis_heading(vector<unsigned> &basis,
+                                      vector<int> &basis_heading) {
     lp_assert(basis_heading.size() >= basis.size());
     unsigned m = basis.size();
     for (unsigned i = 0; i < m; i++) {
@@ -402,26 +392,28 @@ void init_basic_part_of_basis_heading(vector<unsigned> & basis, vector<int> & ba
     }
 }
 
-void init_non_basic_part_of_basis_heading(vector<int> & basis_heading, vector<unsigned> & non_basic_columns) {
+void init_non_basic_part_of_basis_heading(vector<int> &basis_heading,
+                                          vector<unsigned> &non_basic_columns) {
     non_basic_columns.clear();
-    for (int j = basis_heading.size(); j--;){
+    for (int j = basis_heading.size(); j--;) {
         if (basis_heading[j] < 0) {
             non_basic_columns.push_back(j);
             // the index of column j in m_nbasis is (- basis_heading[j] - 1)
-            basis_heading[j] = - static_cast<int>(non_basic_columns.size());
+            basis_heading[j] = -static_cast<int>(non_basic_columns.size());
         }
     }
 }
-void init_basis_heading_and_non_basic_columns_vector(vector<unsigned> & basis,
-                                                     vector<int> & basis_heading,
-                                                     vector<unsigned> & non_basic_columns) {
+void init_basis_heading_and_non_basic_columns_vector(
+    vector<unsigned> &basis, vector<int> &basis_heading,
+    vector<unsigned> &non_basic_columns) {
     init_basic_part_of_basis_heading(basis, basis_heading);
     init_non_basic_part_of_basis_heading(basis_heading, non_basic_columns);
 }
 
-void change_basis(unsigned entering, unsigned leaving, vector<unsigned>& basis, vector<unsigned>& nbasis, vector<int> & basis_heading) {
-    int place_in_basis =  basis_heading[leaving];
-    int place_in_non_basis = - basis_heading[entering] - 1;
+void change_basis(unsigned entering, unsigned leaving, vector<unsigned> &basis,
+                  vector<unsigned> &nbasis, vector<int> &basis_heading) {
+    int place_in_basis = basis_heading[leaving];
+    int place_in_non_basis = -basis_heading[entering] - 1;
     basis_heading[entering] = place_in_basis;
     basis_heading[leaving] = -place_in_non_basis - 1;
     basis[place_in_basis] = entering;
@@ -430,7 +422,8 @@ void change_basis(unsigned entering, unsigned leaving, vector<unsigned>& basis,
 
 int perm_id = 0;
 
-bool get_int_from_args_parser(const char * option, argument_parser & args_parser, unsigned & n) {
+bool get_int_from_args_parser(const char *option, argument_parser &args_parser,
+                              unsigned &n) {
     std::string s = args_parser.get_option_value(option);
     if (!s.empty()) {
         n = atoi(s.c_str());
@@ -439,7 +432,8 @@ bool get_int_from_args_parser(const char * option, argument_parser & args_parser
     return false;
 }
 
-bool get_double_from_args_parser(const char * option, argument_parser & args_parser, double & n) {
+bool get_double_from_args_parser(const char *option,
+                                 argument_parser &args_parser, double &n) {
     std::string s = args_parser.get_option_value(option);
     if (!s.empty()) {
         n = atof(s.c_str());
@@ -448,196 +442,18 @@ bool get_double_from_args_parser(const char * option, argument_parser & args_par
     return false;
 }
 
+void get_time_limit_and_max_iters_from_parser(
+    argument_parser &args_parser, unsigned &time_limit);  // forward definition
 
+int get_random_rows() { return 5 + my_random() % 2; }
 
-
-
-void get_time_limit_and_max_iters_from_parser(argument_parser & args_parser, unsigned & time_limit); // forward definition
-
-
-
-
-
-
-int get_random_rows() {
-    return 5 + my_random() % 2;
-}
-
-int get_random_columns() {
-    return 5 + my_random() % 3;
-}
+int get_random_columns() { return 5 + my_random() % 3; }
 
 int get_random_int() {
-    return -1 + my_random() % 2; // (1.0 + RAND_MAX);
+    return -1 + my_random() % 2;  // (1.0 + RAND_MAX);
 }
 
-#ifndef _WINDOWS
-void fill_file_names(vector<std::string> &file_names,  std::set<std::string> & minimums) {
-    char *home_dir = getenv("HOME");
-    if (home_dir == nullptr) {
-        std::cout << "cannot find home directory, don't know how to find the files";
-        return;
-    }
-    std::string home_dir_str(home_dir);
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/l0redund.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/l1.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/l2.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/l3.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/l4.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/l4fix.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/plan.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/samp2.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/murtagh.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/l0.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/AFIRO.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SC50B.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SC50A.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/KB2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SC105.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/STOCFOR1.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/ADLITTLE.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BLEND.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCAGR7.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SC205.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHARE2B.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/RECIPELP.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/LOTFI.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/VTP-BASE.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHARE1B.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BOEING2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BORE3D.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCORPION.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/CAPRI.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BRANDY.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCAGR25.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCTAP1.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/ISRAEL.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCFXM1.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BANDM.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/E226.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/AGG.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/GROW7.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/ETAMACRO.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/FINNIS.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCSD1.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/STANDATA.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/STANDGUB.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BEACONFD.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/STAIR.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/STANDMPS.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/GFRD-PNC.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCRS8.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BOEING1.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/MODSZK1.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/DEGEN2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/FORPLAN.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/AGG2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/AGG3.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCFXM2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHELL.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/PILOT4.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCSD6.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHIP04S.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SEBA.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/GROW15.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/FFFFF800.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BNL1.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/PEROLD.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/QAP8.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCFXM3.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHIP04L.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/GANGES.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCTAP2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/GROW22.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHIP08S.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/PILOT-WE.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/MAROS.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/STOCFOR2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/25FV47.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHIP12S.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCSD8.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/FIT1P.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SCTAP3.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SIERRA.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/PILOTNOV.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/CZPROB.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/FIT1D.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/PILOT-JA.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHIP08L.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/BNL2.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/NESM.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/CYCLE.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/acc-tight5.mps");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/SHIP12L.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/DEGEN3.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/GREENBEA.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/GREENBEB.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/80BAU3B.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/TRUSS.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/D2Q06C.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/WOODW.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/QAP12.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/D6CUBE.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/PILOT.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/DFL001.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/WOOD1P.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/FIT2P.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/PILOT87.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/STOCFOR3.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/QAP15.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/FIT2D.SIF");
-    file_names.push_back(home_dir_str + "/projects/lp/src/tests/math/lp/test_files/netlib/MAROS-R7.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/FIT2P.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/DFL001.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/D2Q06C.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/80BAU3B.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/GREENBEB.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/GREENBEA.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/BNL2.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/SHIP08L.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/FIT1D.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/SCTAP3.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/SCSD8.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/SCSD6.SIF");
-    minimums.insert("/projects/lp/src/tests/math/lp/test_files/netlib/MAROS-R7.SIF");
-}
-
-void test_out_dir(std::string out_dir) {
-    auto *out_dir_p = opendir(out_dir.c_str());
-    if (out_dir_p == nullptr) {
-        std::cout << "creating directory " << out_dir << std::endl;
-#ifdef LEAN_WINDOWS
-        int res = mkdir(out_dir.c_str());
-#else
-        int res = mkdir(out_dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
-#endif
-        if (res) {
-            std::cout << "Cannot open output directory \"" << out_dir << "\"" << std::endl;
-        }
-        return;
-    }
-    closedir(out_dir_p);
-}
-
-void find_dir_and_file_name(std::string a, std::string & dir, std::string& fn) {
-    // todo: make it system independent
-    size_t last_slash_pos = a.find_last_of('/');
-    if (last_slash_pos >= a.size()) {
-        std::cout << "cannot find file name in " << a << std::endl;
-        throw;
-    }
-    dir = a.substr(0, last_slash_pos);
-    // std::cout << "dir = " << dir << std::endl;
-    fn = a.substr(last_slash_pos + 1);
-    //    std::cout << "fn = " << fn << std::endl;
-}
-
-
-#endif
-
-
-
-std::string read_line(bool & end, std::ifstream & file) {
+std::string read_line(bool &end, std::ifstream &file) {
     std::string s;
     if (!getline(file, s)) {
         end = true;
@@ -647,67 +463,108 @@ std::string read_line(bool & end, std::ifstream & file) {
     return s;
 }
 
-bool contains(std::string const & s, char const * pattern) {
+bool contains(std::string const &s, char const *pattern) {
     return s.find(pattern) != std::string::npos;
 }
 
-
-
-void setup_args_parser(argument_parser & parser) {
+void setup_args_parser(argument_parser &parser) {
     parser.add_option_with_help_string("-monics", "test emonics");
     parser.add_option_with_help_string("-nex_order", "test nex order");
     parser.add_option_with_help_string("-nla_cn", "test cross nornmal form");
     parser.add_option_with_help_string("-nla_sim", "test nex simplify");
-    parser.add_option_with_help_string("-nla_blfmz_mf", "test_basic_lemma_for_mon_zero_from_factor_to_monomial");
-    parser.add_option_with_help_string("-nla_blfmz_fm", "test_basic_lemma_for_mon_zero_from_monomials_to_factor");
-    parser.add_option_with_help_string("-nla_order", "test nla_solver order lemma");
-    parser.add_option_with_help_string("-nla_monot", "test nla_solver order lemma");
+    parser.add_option_with_help_string(
+        "-nla_blfmz_mf", "test_basic_lemma_for_mon_zero_from_factor_to_monomial");
+    parser.add_option_with_help_string(
+        "-nla_blfmz_fm",
+        "test_basic_lemma_for_mon_zero_from_monomials_to_factor");
+    parser.add_option_with_help_string("-nla_order",
+                                       "test nla_solver order lemma");
+    parser.add_option_with_help_string("-nla_monot",
+                                       "test nla_solver order lemma");
     parser.add_option_with_help_string("-nla_tan", "test_tangent_lemma");
     parser.add_option_with_help_string("-nla_bsl", "test_basic_sign_lemma");
     parser.add_option_with_help_string("-horner", "test horner's heuristic");
-    parser.add_option_with_help_string("-nla_blnt_mf", "test_basic_lemma_for_mon_neutral_from_monomial_to_factors");
-    parser.add_option_with_help_string("-nla_blnt_fm", "test_basic_lemma_for_mon_neutral_from_factors_to_monomial");
+    parser.add_option_with_help_string(
+        "-nla_blnt_mf",
+        "test_basic_lemma_for_mon_neutral_from_monomial_to_factors");
+    parser.add_option_with_help_string(
+        "-nla_blnt_fm",
+        "test_basic_lemma_for_mon_neutral_from_factors_to_monomial");
     parser.add_option_with_help_string("-hnf", "test hermite normal form");
     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("--percent_for_enter", "which percent of columns check for entering column");
-    parser.add_option_with_help_string("--totalinf", "minimizes the total infeasibility  instead of diminishing infeasibility of the rows");
-    parser.add_option_with_after_string_with_help("--rep_frq", "the report frequency, in how many iterations print the cost and other info ");
+    parser.add_option_with_help_string("-xyz_sample",
+                                       "run a small interactive scenario");
+    parser.add_option_with_after_string_with_help(
+        "--percent_for_enter",
+        "which percent of columns check for entering column");
+    parser.add_option_with_help_string(
+        "--totalinf",
+        "minimizes the total infeasibility  instead of diminishing "
+        "infeasibility of the rows");
+    parser.add_option_with_after_string_with_help(
+        "--rep_frq",
+        "the report frequency, in how many iterations print the "
+        "cost and other info ");
     parser.add_option_with_help_string("--smt", "smt file format");
-    parser.add_option_with_after_string_with_help("--filelist", "the file containing the list of files");
-    parser.add_option_with_after_string_with_help("--file", "the input file name"); 
+    parser.add_option_with_after_string_with_help(
+        "--filelist", "the file containing the list of files");
+    parser.add_option_with_after_string_with_help("--file",
+                                                  "the input file name");
     parser.add_option_with_after_string_with_help("--random_seed", "random seed");
     parser.add_option_with_help_string("--bp", "bound propagation");
-    parser.add_option_with_help_string("--min", "will look for the minimum for the given file if --file is used; the default is looking for the max");
-    parser.add_option_with_help_string("--max", "will look for the maximum for the given file if --file is used; it is the default behavior");
-    parser.add_option_with_after_string_with_help("--max_iters", "maximum total iterations in a core solver stage");
-    parser.add_option_with_after_string_with_help("--time_limit", "time limit in seconds");
+    parser.add_option_with_help_string(
+        "--min",
+        "will look for the minimum for the given file if --file is "
+        "used; the default is looking for the max");
+    parser.add_option_with_help_string(
+        "--max",
+        "will look for the maximum for the given file if --file is "
+        "used; it is the default behavior");
+    parser.add_option_with_after_string_with_help(
+        "--max_iters", "maximum total iterations in a core solver stage");
+    parser.add_option_with_after_string_with_help("--time_limit",
+                                                  "time limit in seconds");
     parser.add_option_with_help_string("--mpq", "solve for rational numbers");
-    parser.add_option_with_after_string_with_help("--simplex_strategy", "sets simplex strategy for rational number");
+    parser.add_option_with_after_string_with_help(
+        "--simplex_strategy", "sets simplex strategy for rational number");
     parser.add_option_with_help_string("--test_lp_0", "solve a small lp");
-    parser.add_option_with_help_string("--solve_some_mps", "solves a list of mps problems");
-    parser.add_option_with_after_string_with_help("--test_file_directory", "loads files from the directory for testing");
-    parser.add_option_with_after_string_with_help("--out_dir", "setting the output directory for tests, if not set /tmp is used");
+    parser.add_option_with_help_string("--solve_some_mps",
+                                       "solves a list of mps problems");
+    parser.add_option_with_after_string_with_help(
+        "--test_file_directory", "loads files from the directory for testing");
+    parser.add_option_with_after_string_with_help(
+        "--out_dir",
+        "setting the output directory for tests, if not set /tmp is used");
     parser.add_option_with_help_string("--dual", "using the dual simplex solver");
-    parser.add_option_with_help_string("--compare_with_primal", "using the primal simplex solver for comparison");
+    parser.add_option_with_help_string(
+        "--compare_with_primal",
+        "using the primal simplex solver for comparison");
     parser.add_option_with_help_string("--lar", "test lar_solver");
-    parser.add_option_with_after_string_with_help("--maxng", "max iterations without progress");
-    parser.add_option_with_help_string("--randomize_lar", "test randomize functionality");
+    parser.add_option_with_after_string_with_help(
+        "--maxng", "max iterations without progress");
+    parser.add_option_with_help_string("--randomize_lar",
+                                       "test randomize functionality");
     parser.add_option_with_help_string("--smap", "test stacked_map");
     parser.add_option_with_help_string("--term", "simple term test");
-    parser.add_option_with_help_string("--eti"," run a small evidence test for total infeasibility scenario");
-    parser.add_option_with_help_string("--row_inf", "forces row infeasibility search");
+    parser.add_option_with_help_string(
+        "--eti", " run a small evidence test for total infeasibility scenario");
+    parser.add_option_with_help_string("--row_inf",
+                                       "forces row infeasibility search");
     parser.add_option_with_help_string("-pd", "presolve with double solver");
     parser.add_option_with_help_string("--test_int_set", "test int_set");
     parser.add_option_with_help_string("--test_mpq", "test rationals");
     parser.add_option_with_help_string("--test_mpq_np", "test rationals");
-    parser.add_option_with_help_string("--test_mpq_np_plus", "test rationals using plus instead of +=");
+    parser.add_option_with_help_string("--test_mpq_np_plus",
+                                       "test rationals using plus instead of +=");
     parser.add_option_with_help_string("--maximize_term", "test maximize_term()");
+    parser.add_option_with_help_string("--patching", "test patching");
 }
 
-struct fff { int a; int b;};
-
+struct fff {
+    int a;
+    int b;
+};
 
 void test_stacked_unsigned() {
     std::cout << "test stacked unsigned" << std::endl;
@@ -719,24 +576,21 @@ void test_stacked_unsigned() {
     v = 4;
     v.pop();
     lp_assert(v == 2);
-    v ++;
+    v++;
     v++;
     std::cout << "before push v=" << v << std::endl;
     v.push();
     v++;
     v.push();
-    v+=1;
+    v += 1;
     std::cout << "v = " << v << std::endl;
     v.pop(2);
     lp_assert(v == 4);
-    const unsigned & rr = v;
-    std::cout << rr << std:: endl;
-    
+    const unsigned &rr = v;
+    std::cout << rr << std::endl;
 }
 
-void test_stacked_value() {
-    test_stacked_unsigned();
-}
+void test_stacked_value() { test_stacked_unsigned(); }
 
 void test_stacked_vector() {
     std::cout << "test_stacked_vector" << std::endl;
@@ -751,31 +605,29 @@ void test_stacked_vector() {
     v.push_back(3);
     v.push_back(34);
     v.push();
-    v[1]=3;
+    v[1] = 3;
     v[2] = 3;
     v.push();
-    v[0]= 7;
+    v[0] = 7;
     v[1] = 9;
     v.pop(2);
     if (v.size())
-        v[v.size() -1 ] = 7;
+        v[v.size() - 1] = 7;
     v.push();
     v.push_back(33);
     v[0] = 13;
     v.pop();
-        
 }
 
-
 void test_stacked() {
     test_stacked_value();
     test_stacked_vector();
 }
 
-char * find_home_dir() {
+char *find_home_dir() {
 #ifdef _WINDOWS
 #else
-    char * home_dir =   getenv("HOME");
+    char *home_dir = getenv("HOME");
     if (home_dir == nullptr) {
         std::cout << "cannot find home directory" << std::endl;
         return nullptr;
@@ -784,9 +636,8 @@ char * find_home_dir() {
     return nullptr;
 }
 
-
 template <typename T>
-void print_chunk(T * arr, unsigned len) {
+void print_chunk(T *arr, unsigned len) {
     for (unsigned i = 0; i < len; i++) {
         std::cout << arr[i] << ", ";
     }
@@ -794,7 +645,7 @@ void print_chunk(T * arr, unsigned len) {
 }
 
 struct mem_cpy_place_holder {
-    static void mem_copy_hook(int * destination, unsigned num) {
+    static void mem_copy_hook(int *destination, unsigned num) {
         if (destination == nullptr || num == 0) {
             throw "bad parameters";
         }
@@ -803,13 +654,14 @@ struct mem_cpy_place_holder {
 
 void finalize(unsigned ret) {
     /*
-      finalize_util_module();
-      finalize_numerics_module();
-    */
+    finalize_util_module();
+    finalize_numerics_module();
+  */
     //    return ret;
 }
 
-void get_time_limit_and_max_iters_from_parser(argument_parser & args_parser, unsigned & time_limit) {
+void get_time_limit_and_max_iters_from_parser(argument_parser &args_parser,
+                                              unsigned &time_limit) {
     std::string time_limit_string = args_parser.get_option_value("--time_limit");
     if (!time_limit_string.empty()) {
         time_limit = atoi(time_limit_string.c_str());
@@ -818,21 +670,28 @@ void get_time_limit_and_max_iters_from_parser(argument_parser & args_parser, uns
     }
 }
 
-
-std::string create_output_file_name(bool minimize, std::string file_name, bool use_mpq) {
-    std::string ret = file_name + "_lp_tst_" +  (minimize?"min":"max");
-    if (use_mpq) return ret + "_mpq.out";
+std::string create_output_file_name(bool minimize, std::string file_name,
+                                    bool use_mpq) {
+    std::string ret = file_name + "_lp_tst_" + (minimize ? "min" : "max");
+    if (use_mpq)
+        return ret + "_mpq.out";
     return ret + ".out";
 }
 
-std::string create_output_file_name_for_glpsol(bool minimize, std::string file_name){
-    return file_name + (minimize?"_min":"_max") + "_glpk_out";
+std::string create_output_file_name_for_glpsol(bool minimize,
+                                               std::string file_name) {
+    return file_name + (minimize ? "_min" : "_max") + "_glpk_out";
 }
 
-int run_glpk(std::string file_name, std::string glpk_out_file_name, bool minimize, unsigned time_limit) {
-    std::string minmax(minimize?"--min":"--max");
-    std::string tmlim = time_limit > 0 ? std::string(" --tmlim ")  + std::to_string(time_limit)+ " ":std::string();
-    std::string command_line =  std::string("glpsol --nointopt --nomip ") +  minmax + tmlim +  + " -o " + glpk_out_file_name +" " + file_name + " > /dev/null";
+int run_glpk(std::string file_name, std::string glpk_out_file_name,
+             bool minimize, unsigned time_limit) {
+    std::string minmax(minimize ? "--min" : "--max");
+    std::string tmlim = time_limit > 0 ? std::string(" --tmlim ") +
+                                             std::to_string(time_limit) + " "
+                                       : std::string();
+    std::string command_line = std::string("glpsol --nointopt --nomip ") +
+                               minmax + tmlim + +" -o " + glpk_out_file_name +
+                               " " + file_name + " > /dev/null";
     return system(command_line.c_str());
 }
 
@@ -857,22 +716,13 @@ std::string get_status(std::string file_name) {
     throw 0;
 }
 
-
-
-
 struct sort_pred {
-    bool operator()(const std::pair<std::string, int> &left, const std::pair<std::string, int> &right) {
+    bool operator()(const std::pair<std::string, int> &left,
+                    const std::pair<std::string, int> &right) {
         return left.second < right.second;
     }
 };
 
-
-
-
-
-
-
-
 vector<std::string> get_file_names_from_file_list(std::string filelist) {
     std::ifstream file(filelist);
     if (!file.is_open()) {
@@ -892,17 +742,16 @@ vector<std::string> get_file_names_from_file_list(std::string filelist) {
     return ret;
 }
 
-
 void test_numeric_pair() {
     numeric_pair<lp::mpq> a;
     numeric_pair<lp::mpq> b(2, lp::mpq(6, 2));
     a = b;
     numeric_pair<lp::mpq> c(0.1, 0.5);
-    a += 2*c;
+    a += 2 * c;
     a -= c;
-    lp_assert (a == b + c);
+    lp_assert(a == b + c);
     numeric_pair<lp::mpq> d = a * 2;
-    std::cout << a  << std::endl;
+    std::cout << a << std::endl;
     lp_assert(b == b);
     lp_assert(b < a);
     lp_assert(b <= a);
@@ -913,12 +762,12 @@ void test_numeric_pair() {
     lp_assert(a < 2 * b);
     lp_assert(b + b > a);
     lp_assert(lp::mpq(2.1) * b + b > a);
-    lp_assert(-b * lp::mpq(2.1) - b < lp::mpq(0.99)  * a);
-    std::cout << - b * lp::mpq(2.1) - b << std::endl;
-    lp_assert(-b *(lp::mpq(2.1) + 1) == - b * lp::mpq(2.1) - b);
+    lp_assert(-b * lp::mpq(2.1) - b < lp::mpq(0.99) * a);
+    std::cout << -b * lp::mpq(2.1) - b << std::endl;
+    lp_assert(-b * (lp::mpq(2.1) + 1) == -b * lp::mpq(2.1) - b);
 }
 
-void get_matrix_dimensions(std::ifstream & f, unsigned & m, unsigned & n) {
+void get_matrix_dimensions(std::ifstream &f, unsigned &m, unsigned &n) {
     std::string line;
     getline(f, line);
     getline(f, line);
@@ -929,13 +778,10 @@ void get_matrix_dimensions(std::ifstream & f, unsigned & m, unsigned & n) {
     n = atoi(r[1].c_str());
 }
 
-
 void print_st(lp_status status) {
     std::cout << lp_status_to_string(status) << std::endl;
 }
 
-
-
 void test_term() {
     lar_solver solver;
     unsigned _x = 0;
@@ -953,10 +799,10 @@ void test_term() {
     solver.add_var_bound(x_plus_y, lconstraint_kind::LE, mpq(14, 3));
     pairs.pop_back();
     pairs.push_back(std::pair<mpq, var_index>(mpq(-1), y));
-    unsigned x_minus_y =  solver.add_term(pairs, ti++);
+    unsigned x_minus_y = solver.add_term(pairs, ti++);
     solver.add_var_bound(x_minus_y, lconstraint_kind::GE, mpq(5, 3));
     solver.add_var_bound(x_minus_y, lconstraint_kind::LE, mpq(14, 3));
-    auto status = solver.solve();   
+    auto status = solver.solve();
     std::cout << lp_status_to_string(status) << std::endl;
     std::unordered_map<var_index, mpq> model;
     if (status != lp_status::OPTIMAL) {
@@ -971,35 +817,36 @@ void test_term() {
     solver.set_int_solver(&i_s);
     int_cube cuber(i_s);
     lia_move m = cuber();
-    
-    std::cout <<"\n" << lia_move_to_string(m) << std::endl;
+
+    std::cout << "\n"
+              << lia_move_to_string(m) << std::endl;
     model.clear();
     solver.get_model(model);
-    for (auto & t : model) {
-        std::cout << solver.get_variable_name(t.first) << " = " << t.second.get_double() << ",";
+    for (auto &t : model) {
+        std::cout << solver.get_variable_name(t.first) << " = "
+                  << t.second.get_double() << ",";
     }
 
     std::cout << "\ntableu after cube\n";
     solver.pp(std::cout).print();
     std::cout << "Ax_is_correct = " << solver.ax_is_correct() << "\n";
-    
 }
 
-void test_evidence_for_total_inf_simple(argument_parser & args_parser) {
+void test_evidence_for_total_inf_simple(argument_parser &args_parser) {
     lar_solver solver;
     var_index x = solver.add_var(0, false);
     var_index y = solver.add_var(1, false);
     solver.add_var_bound(x, LE, mpq(-1));
     solver.add_var_bound(y, GE, mpq(0));
     vector<std::pair<mpq, var_index>> ls;
-    
+
     ls.push_back(std::pair<mpq, var_index>(mpq(1), x));
     ls.push_back(std::pair<mpq, var_index>(mpq(1), y));
 
     unsigned j = solver.add_term(ls, 1);
     solver.add_var_bound(j, GE, mpq(1));
     ls.pop_back();
-    ls.push_back(std::pair<mpq, var_index>(- mpq(1), y));
+    ls.push_back(std::pair<mpq, var_index>(-mpq(1), y));
     j = solver.add_term(ls, 2);
     solver.add_var_bound(j, GE, mpq(0));
     auto status = solver.solve();
@@ -1009,21 +856,19 @@ void test_evidence_for_total_inf_simple(argument_parser & args_parser) {
 }
 void test_bound_propagation_one_small_sample1() {
     /*
-      (<= (+ a (* (- 1.0) b)) 0.0)
-      (<= (+ b (* (- 1.0) x_13)) 0.0)
-      --> (<= (+ a (* (- 1.0) c)) 0.0)
+    (<= (+ a (* (- 1.0) b)) 0.0)
+    (<= (+ b (* (- 1.0) x_13)) 0.0)
+    --> (<= (+ a (* (- 1.0) c)) 0.0)
 
-      the inequality on (<= a c) is obtained from a triangle inequality (<= a b) (<= b c).
-      If b becomes basic variable, then it is likely the old solver ends up with a row that implies (<= a c).
-      a - b <= 0.0
-      b  - c <= 0.0
+    the inequality on (<= a c) is obtained from a triangle inequality (<= a b)
+    (<= b c). If b becomes basic variable, then it is likely the old solver ends
+    up with a row that implies (<= a c). a - b <= 0.0 b  - c <= 0.0
 
-      got to get a <= c
-    */
-    std::function<bool (unsigned, bool, bool, const mpq & )> bound_is_relevant =
-        [&](unsigned j, bool is_lower_bound, bool strict, const rational& bound_val) {
-        return true; 
-    };   
+    got to get a <= c
+  */
+    std::function<bool(unsigned, bool, bool, const mpq &)> bound_is_relevant =
+        [&](unsigned j, bool is_lower_bound, bool strict,
+            const rational &bound_val) { return true; };
     lar_solver ls;
     unsigned a = ls.add_var(0, false);
     unsigned b = ls.add_var(1, false);
@@ -1048,8 +893,8 @@ void test_bound_propagation_one_small_sample1() {
     // ls.solve();
     // my_bound_propagator bp(ls);
     // ls.propagate_bounds_for_touched_rows(bp);
-    // std::cout << " bound ev from test_bound_propagation_one_small_sample1" << std::endl;
-    // for (auto & be : bp.m_ibounds)  {
+    // std::cout << " bound ev from test_bound_propagation_one_small_sample1" <<
+    // std::endl; for (auto & be : bp.m_ibounds)  {
     //     std::cout << "bound\n";
     //     ls.print_implied_bound(be, std::cout);
     // } // todo: restore test
@@ -1058,35 +903,40 @@ void test_bound_propagation_one_small_sample1() {
 void test_bound_propagation_one_small_samples() {
     test_bound_propagation_one_small_sample1();
     /*
-      (>= x_46 0.0)
-      (<= x_29 0.0)
-      (not (<= x_68 0.0))
-      (<= (+ (* (/ 1001.0 1998.0) x_10) (* (- 1.0) x_151) x_68) (- (/ 1001.0 999.0)))
-      (<= (+ (* (/ 1001.0 999.0) x_9)
-      (* (- 1.0) x_152)
-      (* (/ 1001.0 999.0) x_151)
-      (* (/ 1001.0 999.0) x_68))
-      (- (/ 1502501.0 999000.0)))
-      (not (<= (+ (* (/ 999.0 2.0) x_10) (* (- 1.0) x_152) (* (- (/ 999.0 2.0)) x_151))
-      (/ 1001.0 2.0)))
-      (not (<= x_153 0.0))z
-      (>= (+ x_9 (* (- (/ 1001.0 999.0)) x_10) (* (- 1.0) x_153) (* (- 1.0) x_68))
-      (/ 5003.0 1998.0))
-      --> (not (<= (+ x_10 x_46 (* (- 1.0) x_29)) 0.0))
+    (>= x_46 0.0)
+    (<= x_29 0.0)
+    (not (<= x_68 0.0))
+    (<= (+ (* (/ 1001.0 1998.0) x_10) (* (- 1.0) x_151) x_68) (- (/ 1001.0
+    999.0)))
+    (<= (+ (* (/ 1001.0 999.0) x_9)
+    (* (- 1.0) x_152)
+    (* (/ 1001.0 999.0) x_151)
+    (* (/ 1001.0 999.0) x_68))
+    (- (/ 1502501.0 999000.0)))
+    (not (<= (+ (* (/ 999.0 2.0) x_10) (* (- 1.0) x_152) (* (- (/ 999.0 2.0))
+    x_151))
+    (/ 1001.0 2.0)))
+    (not (<= x_153 0.0))z
+    (>= (+ x_9 (* (- (/ 1001.0 999.0)) x_10) (* (- 1.0) x_153) (* (- 1.0) x_68))
+    (/ 5003.0 1998.0))
+    --> (not (<= (+ x_10 x_46 (* (- 1.0) x_29)) 0.0))
 
-      and 
+    and
 
-      (<= (+ a (* (- 1.0) b)) 0.0)
-      (<= (+ b (* (- 1.0) x_13)) 0.0)
-      --> (<= (+ a (* (- 1.0) x_13)) 0.0)
+    (<= (+ a (* (- 1.0) b)) 0.0)
+    (<= (+ b (* (- 1.0) x_13)) 0.0)
+    --> (<= (+ a (* (- 1.0) x_13)) 0.0)
 
-      In the first case, there typically are no atomic formulas for bounding x_10. So there is never some
-      basic lemma of the form (>= x46 0), (<= x29 0), (>= x10 0) -> (not (<= (+ x10 x46 (- x29)) 0)).
-      Instead the bound on x_10 falls out from a bigger blob of constraints. 
+    In the first case, there typically are no atomic formulas for bounding x_10.
+    So there is never some basic lemma of the form (>= x46 0), (<= x29 0), (>=
+    x10 0) -> (not (<= (+ x10 x46 (- x29)) 0)). Instead the bound on x_10 falls
+    out from a bigger blob of constraints.
 
-      In the second case, the inequality on (<= x19 x13) is obtained from a triangle inequality (<= x19 x9) (<= x9 x13).
-      If x9 becomes basic variable, then it is likely the old solver ends up with a row that implies (<= x19 x13).
-    */
+    In the second case, the inequality on (<= x19 x13) is obtained from a
+    triangle inequality (<= x19 x9) (<= x9 x13). If x9 becomes basic variable,
+    then it is likely the old solver ends up with a row that implies (<= x19
+    x13).
+  */
 }
 void test_bound_propagation_one_row() {
     lar_solver ls;
@@ -1102,7 +952,7 @@ void test_bound_propagation_one_row() {
     // ls.solve();
     // my_bound_propagator bp(ls);
     // ls.propagate_bounds_for_touched_rows(bp);
-} 
+}
 void test_bound_propagation_one_row_with_bounded_vars() {
     lar_solver ls;
     unsigned x0 = ls.add_var(0, false);
@@ -1134,7 +984,7 @@ void test_bound_propagation_one_row_mixed() {
     // ls.solve();
     // my_bound_propagator bp(ls);
     // ls.propagate_bounds_for_touched_rows(bp);
-} 
+}
 
 void test_bound_propagation_two_rows() {
     lar_solver ls;
@@ -1158,7 +1008,7 @@ void test_bound_propagation_two_rows() {
     // ls.solve();
     // my_bound_propagator bp(ls);
     // ls.propagate_bounds_for_touched_rows(bp);
-} 
+}
 
 void test_total_case_u() {
     std::cout << "test_total_case_u\n";
@@ -1180,14 +1030,15 @@ void test_total_case_u() {
     // my_bound_propagator bp(ls);
     // ls.propagate_bounds_for_touched_rows(bp);
 }
-bool contains_j_kind(unsigned j, lconstraint_kind kind, const mpq & rs, const vector<implied_bound> & ev) {
-    for (auto & e : ev) {
+bool contains_j_kind(unsigned j, lconstraint_kind kind, const mpq &rs,
+                     const vector<implied_bound> &ev) {
+    for (auto &e : ev) {
         if (e.m_j == j && e.m_bound == rs && e.kind() == kind)
             return true;
     }
     return false;
 }
-void test_total_case_l(){
+void test_total_case_l() {
     std::cout << "test_total_case_l\n";
     lar_solver ls;
     unsigned x = ls.add_var(0, false);
@@ -1218,7 +1069,6 @@ void test_bound_propagation() {
     test_bound_propagation_two_rows();
     test_bound_propagation_one_row_mixed();
     test_total_case_l();
-    
 }
 
 void test_int_set() {
@@ -1236,36 +1086,35 @@ void test_int_set() {
     s.insert(2);
     s.clear();
     lp_assert(s.size() == 0);
-    
-    
 }
 
 void test_rationals_no_numeric_pairs() {
     stopwatch sw;
 
     vector<mpq> c;
-    for (unsigned j = 0; j < 10; j ++)
-        c.push_back(mpq(my_random()%100, 1 + my_random()%100 )); 
-    
+    for (unsigned j = 0; j < 10; j++)
+        c.push_back(mpq(my_random() % 100, 1 + my_random() % 100));
+
     vector<mpq> x;
-    for (unsigned j = 0; j < 10; j ++)
-        x.push_back(mpq(my_random()%100, 1 + my_random()%100 )); 
+    for (unsigned j = 0; j < 10; j++)
+        x.push_back(mpq(my_random() % 100, 1 + my_random() % 100));
 
     unsigned k = 500000;
-    mpq r=zero_of_type<mpq>();
+    mpq r = zero_of_type<mpq>();
     sw.start();
-    
-    for (unsigned j = 0; j < k; j++){
+
+    for (unsigned j = 0; j < k; j++) {
         mpq val = zero_of_type<mpq>();
-        for (unsigned j=0;j< c.size(); j++){
-            val += c[j]*x[j];
+        for (unsigned j = 0; j < c.size(); j++) {
+            val += c[j] * x[j];
         }
-        
+
         r += val;
     }
-            
+
     sw.stop();
-    std::cout << "operation with rationals no pairs " << sw.get_seconds() << std::endl;
+    std::cout << "operation with rationals no pairs " << sw.get_seconds()
+              << std::endl;
     std::cout << T_to_string(r) << std::endl;
 }
 
@@ -1273,64 +1122,62 @@ void test_rationals_no_numeric_pairs_plus() {
     stopwatch sw;
 
     vector<mpq> c;
-    for (unsigned j = 0; j < 10; j ++)
-        c.push_back(mpq(my_random()%100, 1 + my_random()%100 )); 
-    
+    for (unsigned j = 0; j < 10; j++)
+        c.push_back(mpq(my_random() % 100, 1 + my_random() % 100));
+
     vector<mpq> x;
-    for (unsigned j = 0; j < 10; j ++)
-        x.push_back(mpq(my_random()%100, 1 + my_random()%100 )); 
+    for (unsigned j = 0; j < 10; j++)
+        x.push_back(mpq(my_random() % 100, 1 + my_random() % 100));
 
     unsigned k = 500000;
-    mpq r=zero_of_type<mpq>();
+    mpq r = zero_of_type<mpq>();
     sw.start();
-    
-    for (unsigned j = 0; j < k; j++){
+
+    for (unsigned j = 0; j < k; j++) {
         mpq val = zero_of_type<mpq>();
-        for (unsigned j=0;j< c.size(); j++){
-            val = val + c[j]*x[j];
+        for (unsigned j = 0; j < c.size(); j++) {
+            val = val + c[j] * x[j];
         }
-        
+
         r = r + val;
     }
-            
+
     sw.stop();
-    std::cout << "operation with rationals no pairs " << sw.get_seconds() << std::endl;
+    std::cout << "operation with rationals no pairs " << sw.get_seconds()
+              << std::endl;
     std::cout << T_to_string(r) << std::endl;
 }
 
-
-
 void test_rationals() {
     stopwatch sw;
 
     vector<rational> c;
-    for (unsigned j = 0; j < 10; j ++)
-        c.push_back(rational(my_random()%100, 1 + my_random()%100)); 
+    for (unsigned j = 0; j < 10; j++)
+        c.push_back(rational(my_random() % 100, 1 + my_random() % 100));
 
-    
-    
     vector<numeric_pair<rational>> x;
-    for (unsigned j = 0; j < 10; j ++)
-        x.push_back(numeric_pair<rational>(rational(my_random()%100, 1 + my_random()%100 ))); 
+    for (unsigned j = 0; j < 10; j++)
+        x.push_back(numeric_pair<rational>(
+            rational(my_random() % 100, 1 + my_random() % 100)));
 
     std::cout << "x = ";
     print_vector(x, std::cout);
-    
+
     unsigned k = 1000000;
-    numeric_pair<rational> r=zero_of_type<numeric_pair<rational>>();
+    numeric_pair<rational> r = zero_of_type<numeric_pair<rational>>();
     sw.start();
-    
+
     for (unsigned j = 0; j < k; j++) {
         for (unsigned i = 0; i < c.size(); i++) {
-            r+= c[i] * x[i];
+            r += c[i] * x[i];
         }
-    }   
+    }
     sw.stop();
     std::cout << "operation with rationals " << sw.get_seconds() << std::endl;
     std::cout << T_to_string(r) << std::endl;
 }
 
-void get_random_interval(bool& neg_inf, bool& pos_inf, int& x, int &y) {
+void get_random_interval(bool &neg_inf, bool &pos_inf, int &x, int &y) {
     int i = my_random() % 10;
     if (i == 0) {
         neg_inf = true;
@@ -1346,29 +1193,29 @@ void get_random_interval(bool& neg_inf, bool& pos_inf, int& x, int &y) {
         if (!neg_inf) {
             y = x + my_random() % (101 - x);
             lp_assert(y >= x);
-        }
-        else {
+        } else {
             y = my_random() % 100;
         }
     }
-    lp_assert((neg_inf || (0 <= x && x <= 100)) && (pos_inf || (0 <= y && y <= 100)));
+    lp_assert((neg_inf || (0 <= x && x <= 100)) &&
+              (pos_inf || (0 <= y && y <= 100)));
 }
 
-
 void test_gomory_cut_0() {
     gomory_test g(
-        [](unsigned j) { return "v" + T_to_string(j);} // name_function_p
+        [](unsigned j) { return "v" + T_to_string(j); }  // name_function_p
         ,
-        [](unsigned j) { //get_value_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);
+            if (j == 3)
+                return mpq(3);
             UNREACHABLE();
             return zero_of_type<mpq>();
         },
-        [](unsigned j) { // at_low_p
+        [](unsigned j) {  // at_low_p
             if (j == 1)
                 return false;
             if (j == 2)
@@ -1378,7 +1225,7 @@ void test_gomory_cut_0() {
             UNREACHABLE();
             return false;
         },
-        [](unsigned j) { // at_upper
+        [](unsigned j) {  // at_upper
             if (j == 1)
                 return false;
             if (j == 2)
@@ -1388,9 +1235,9 @@ void test_gomory_cut_0() {
             UNREACHABLE();
             return false;
         },
-        [](unsigned j) { // lower_bound
+        [](unsigned j) {  // lower_bound
             if (j == 1) {
-                UNREACHABLE(); //unlimited from below
+                UNREACHABLE();  // unlimited from below
                 return impq(0);
             }
             if (j == 2)
@@ -1400,9 +1247,9 @@ void test_gomory_cut_0() {
             UNREACHABLE();
             return impq(0);
         },
-        [](unsigned j) { // upper
+        [](unsigned j) {  // upper
             if (j == 1) {
-                UNREACHABLE(); //unlimited from above
+                UNREACHABLE();  // unlimited from above
                 return impq(0);
             }
             if (j == 2)
@@ -1412,9 +1259,7 @@ void test_gomory_cut_0() {
             UNREACHABLE();
             return impq(0);
         },
-        [] (unsigned) { return 0; },
-        [] (unsigned) { return 0; }
-                  );
+        [](unsigned) { return 0; }, [](unsigned) { return 0; });
     lar_term t;
     mpq k;
     explanation expl;
@@ -1423,14 +1268,14 @@ void test_gomory_cut_0() {
     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);
+    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) { return "v" + T_to_string(j); }  // name_function_p
         ,
-        [](unsigned j) { //get_value_p
+        [](unsigned j) {  // get_value_p
             if (j == 1)
                 return mpq(-2);
             if (j == 2)
@@ -1440,7 +1285,7 @@ void test_gomory_cut_1() {
             UNREACHABLE();
             return zero_of_type<mpq>();
         },
-        [](unsigned j) { // at_low_p
+        [](unsigned j) {  // at_low_p
             if (j == 1)
                 return false;
             if (j == 2)
@@ -1450,7 +1295,7 @@ void test_gomory_cut_1() {
             UNREACHABLE();
             return false;
         },
-        [](unsigned j) { // at_upper
+        [](unsigned j) {  // at_upper
             if (j == 1)
                 return true;
             if (j == 2)
@@ -1460,9 +1305,9 @@ void test_gomory_cut_1() {
             UNREACHABLE();
             return false;
         },
-        [](unsigned j) { // lower_bound
+        [](unsigned j) {  // lower_bound
             if (j == 1) {
-                UNREACHABLE(); //unlimited from below
+                UNREACHABLE();  // unlimited from below
                 return impq(0);
             }
             if (j == 2)
@@ -1472,7 +1317,7 @@ void test_gomory_cut_1() {
             UNREACHABLE();
             return impq(0);
         },
-        [](unsigned j) { // upper
+        [](unsigned j) {  // upper
             if (j == 1) {
                 return impq(-2);
             }
@@ -1483,9 +1328,7 @@ void test_gomory_cut_1() {
             UNREACHABLE();
             return impq(0);
         },
-        [] (unsigned) { return 0; },
-        [] (unsigned) { return 0; }
-                  );
+        [](unsigned) { return 0; }, [](unsigned) { return 0; });
     lar_term t;
     mpq k;
     explanation expl;
@@ -1494,10 +1337,10 @@ void test_gomory_cut_1() {
     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);
+    g.mk_gomory_cut(t, k, expl, inf_col, row);
 }
 
-void call_hnf(general_matrix & A);
+void call_hnf(general_matrix &A);
 
 void test_hnf_m_less_than_n() {
 #ifdef Z3DEBUG
@@ -1547,7 +1390,6 @@ void test_hnf_m_greater_than_n() {
 #endif
 }
 
-
 void cutting_the_mix_example_1() {
     mpq sev(7);
     mpq nine(9);
@@ -1560,9 +1402,9 @@ void cutting_the_mix_example_1() {
     hnf_calc::extended_gcd_minimal_uv(-nine, -nine, d, u, vv);
     std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl;
 
-    hnf_calc::extended_gcd_minimal_uv(-sev*2, sev, d, u, vv);
+    hnf_calc::extended_gcd_minimal_uv(-sev * 2, sev, d, u, vv);
     std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl;
-    
+
     hnf_calc::extended_gcd_minimal_uv(mpq(24), mpq(-7), d, u, vv);
     std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl;
     hnf_calc::extended_gcd_minimal_uv(-mpq(24), mpq(7), d, u, vv);
@@ -1578,7 +1420,7 @@ void cutting_the_mix_example_1() {
 
 #ifdef Z3DEBUG
 
-void fill_general_matrix(general_matrix & M) {
+void fill_general_matrix(general_matrix &M) {
     unsigned m = M.row_count();
     unsigned n = M.column_count();
     for (unsigned i = 0; i < m; i++)
@@ -1586,14 +1428,14 @@ void fill_general_matrix(general_matrix & M) {
             M[i][j] = mpq(static_cast<int>(my_random() % 13) - 6);
 }
 
-void call_hnf(general_matrix& A) {
+void call_hnf(general_matrix &A) {
     svector<unsigned> r;
-    mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r, mpq((int)1000000000));
+    mpq d =
+        hnf_calc::determinant_of_rectangular_matrix(A, r, mpq((int)1000000000));
     A.shrink_to_rank(r);
     hnf<general_matrix> h(A, d);
 }
 
-
 void test_hnf_for_dim(int m) {
     general_matrix M(m, m + my_random() % m);
     fill_general_matrix(M);
@@ -1643,7 +1485,7 @@ void test_hnf_3_3() {
     v.push_back(mpq(-4));
     v.push_back(mpq(-3));
     A.push_row(v);
-    
+
     call_hnf(A);
     std::cout << "test_hnf_3_3 passed" << std::endl;
 }
@@ -1768,7 +1610,7 @@ void test_larger_generated_hnf() {
 void test_maximize_term() {
     std::cout << "test_maximize_term\n";
     lar_solver solver;
-    int_solver i_solver(solver); // have to create it too
+    int_solver i_solver(solver);  // have to create it too
     unsigned _x = 0;
     unsigned _y = 1;
     var_index x = solver.add_var(_x, false);
@@ -1780,9 +1622,9 @@ void test_maximize_term() {
     term_ls.clear();
     term_ls.push_back(std::pair<mpq, var_index>(mpq(2), x));
     term_ls.push_back(std::pair<mpq, var_index>(mpq(2), y));
-    
+
     unsigned term_2x_pl_2y = solver.add_term(term_ls, -1);
-    solver.add_var_bound(term_x_min_y,  LE, zero_of_type<mpq>());
+    solver.add_var_bound(term_x_min_y, LE, zero_of_type<mpq>());
     solver.add_var_bound(term_2x_pl_2y, LE, mpq(5));
     solver.find_feasible_solution();
     lp_assert(solver.get_status() == lp_status::OPTIMAL);
@@ -1790,7 +1632,7 @@ void test_maximize_term() {
     std::unordered_map<var_index, mpq> model;
     solver.get_model(model);
     for (auto p : model) {
-        std::cout<< "v[" << p.first << "] = " << p.second << std::endl;
+        std::cout << "v[" << p.first << "] = " << p.second << std::endl;
     }
     std::cout << "calling int_solver\n";
     explanation ex;
@@ -1798,14 +1640,13 @@ void test_maximize_term() {
     VERIFY(lm == lia_move::sat);
     impq term_max;
     lp_status st = solver.maximize_term(term_2x_pl_2y, term_max);
-    
+
     std::cout << "status = " << lp_status_to_string(st) << std::endl;
     std::cout << "term_max = " << term_max << std::endl;
     solver.get_model(model);
     for (auto p : model) {
-        std::cout<< "v[" << p.first << "] = " << p.second << std::endl;
+        std::cout << "v[" << p.first << "] = " << p.second << std::endl;
     }
-    
 }
 #ifdef Z3DEBUG
 void test_hnf() {
@@ -1816,7 +1657,7 @@ void test_hnf() {
     test_hnf_4_4();
     test_hnf_5_5();
     test_hnf_2_2();
-    for (unsigned k=1000; k>0; k--)
+    for (unsigned k = 1000; k > 0; k--)
         for (int i = 1; i < 8; i++)
             test_hnf_for_dim(i);
     cutting_the_mix_example_1();
@@ -1829,13 +1670,9 @@ void test_gomory_cut() {
     test_gomory_cut_1();
 }
 
-void test_nla_order_lemma() {
-    nla::test_order_lemma();
-}
+void test_nla_order_lemma() { nla::test_order_lemma(); }
 
-
-void test_lp_local(int argn, char**argv) {
-    
+void test_lp_local(int argn, char **argv) {
     // initialize_util_module();
     // initialize_numerics_module();
     int ret;
@@ -1854,7 +1691,10 @@ void test_lp_local(int argn, char**argv) {
         return finalize(0);
     }
 
-    
+    if (args_parser.option_is_used("--patching")) {
+        test_patching();
+        return finalize(0);
+    }
     if (args_parser.option_is_used("-nla_cn")) {
 #ifdef Z3DEBUG
         nla::test_cn();
@@ -1874,7 +1714,6 @@ void test_lp_local(int argn, char**argv) {
         return finalize(0);
     }
 
-    
     if (args_parser.option_is_used("-nla_order")) {
 #ifdef Z3DEBUG
         test_nla_order_lemma();
@@ -1882,7 +1721,6 @@ void test_lp_local(int argn, char**argv) {
         return finalize(0);
     }
 
-    
     if (args_parser.option_is_used("-nla_monot")) {
 #ifdef Z3DEBUG
         nla::test_monotone_lemma();
@@ -1890,68 +1728,67 @@ void test_lp_local(int argn, char**argv) {
         return finalize(0);
     }
 
-    if (args_parser.option_is_used("-nla_bsl")) { 
+    if (args_parser.option_is_used("-nla_bsl")) {
 #ifdef Z3DEBUG
         nla::test_basic_sign_lemma();
 #endif
         return finalize(0);
     }
 
-    if (args_parser.option_is_used("-nla_horner")) { 
+    if (args_parser.option_is_used("-nla_horner")) {
 #ifdef Z3DEBUG
         nla::test_horner();
 #endif
         return finalize(0);
     }
 
-    if (args_parser.option_is_used("-nla_tan")) { 
+    if (args_parser.option_is_used("-nla_tan")) {
 #ifdef Z3DEBUG
         nla::test_tangent_lemma();
 #endif
         return finalize(0);
     }
 
-    if (args_parser.option_is_used("-nla_blfmz_mf")) { 
+    if (args_parser.option_is_used("-nla_blfmz_mf")) {
 #ifdef Z3DEBUG
         nla::test_basic_lemma_for_mon_zero_from_monomial_to_factors();
 #endif
         return finalize(0);
     }
 
-    if (args_parser.option_is_used("-nla_blfmz_fm")) { 
+    if (args_parser.option_is_used("-nla_blfmz_fm")) {
 #ifdef Z3DEBUG
         nla::test_basic_lemma_for_mon_zero_from_factors_to_monomial();
 #endif
         return finalize(0);
     }
 
-    if (args_parser.option_is_used("-nla_blnt_mf")) { 
+    if (args_parser.option_is_used("-nla_blnt_mf")) {
 #ifdef Z3DEBUG
         nla::test_basic_lemma_for_mon_neutral_from_monomial_to_factors();
 #endif
         return finalize(0);
     }
 
-    if (args_parser.option_is_used("-nla_blnt_fm")) { 
+    if (args_parser.option_is_used("-nla_blnt_fm")) {
 #ifdef Z3DEBUG
         nla::test_basic_lemma_for_mon_neutral_from_factors_to_monomial();
 #endif
         return finalize(0);
     }
-   
+
     if (args_parser.option_is_used("-hnf")) {
 #ifdef Z3DEBUG
         test_hnf();
 #endif
         return finalize(0);
     }
-    
+
     if (args_parser.option_is_used("-gomory")) {
         test_gomory_cut();
         return finalize(0);
     }
 
-    
     if (args_parser.option_is_used("--test_int_set")) {
         test_int_set();
         return finalize(0);
@@ -1960,10 +1797,120 @@ void test_lp_local(int argn, char**argv) {
         test_bound_propagation();
         return finalize(0);
     }
-    
-    return finalize(0); // has_violations() ? 1 : 0);
+
+    return finalize(0);  // has_violations() ? 1 : 0);
 }
-}
-void tst_lp(char ** argv, int argc, int& i) {
+}  // namespace lp
+void tst_lp(char **argv, int argc, int &i) {
     lp::test_lp_local(argc - 2, argv + 2);
 }
+// clang-format on
+bool coprime(int a, int b) {
+    return gcd(rational(a), rational(b)).is_one();
+}
+bool coprime(rational &a, rational &b) {
+    return gcd(a, b).is_one();
+}
+void asserts_on_patching(const rational &x, const rational &alpha) {
+    auto a1 = numerator(alpha);
+    auto a2 = denominator(alpha);
+    auto x1 = numerator(x);
+    auto x2 = denominator(x);
+    lp_assert(!a1.is_zero());
+    lp_assert(abs(a1) < abs(a2));
+    lp_assert(coprime(a1, a2));
+    lp_assert(!x1.is_zero());
+    lp_assert(abs(x1) < abs(x2));
+    lp_assert(coprime(x1, x2));
+    lp_assert((a2 / x2).is_int());
+}
+bool get_patching_deltas(const rational &x, const rational &alpha, rational &delta_0, rational &delta_1) {
+    std::cout << "get_patching_deltas(" << x << ", " << alpha << ")" << std::endl;
+    auto a1 = numerator(alpha);
+    auto a2 = denominator(alpha);
+    auto x1 = numerator(x);
+    auto x2 = denominator(x);
+    // delta has to be integral.
+    // We need to find delta such that x1/x2 + (a1/a2)*delta is integral.
+    // Then a2*x1/x2 + a1*delta is integral, that means that t = a2/x2 is integral.
+    // We established that a2 = x2*t
+    // Then x1 + a1*delta*(x2/a2)  = x1 + a1*(delta/t)  is integral. Taking into account
+    // that t and a1 are coprime we have delta = t*k, where k is an integer.
+    rational t = a2 / x2;
+    std::cout << "t = " << t << std::endl;
+    // Now we have x1/x2 + (a1/x2)*k is integral, or (x1 + a1*k)/x2 is integral.
+    // It is equivalent to x1 + a1*k = x2*m, where m is an integer
+    // We know that a2 and a1 are coprime, and x2 divides a2, so x2 and a1 are coprime.
+    rational u, v;
+    auto g = gcd(a1, x2, u, v);
+    lp_assert(g.is_one() && u.is_int() && v.is_int() && g == u * a1 + v * x2);
+    std::cout << "u = " << u << ", v = " << v << std::endl;
+    std::cout << "x= " << (x1 / x2) << std::endl;
+    std::cout << "x + (a1 / a2) * (-u * t) * x1 = " << x + (a1 / a2) * (-u * t) * x1 << std::endl;
+    lp_assert((x + (a1 / a2) * (-u * t) * x1).is_int());
+    // 1 = (u- l*x2 ) * a1 + (v + l*a1)*x2, for every integer l.
+    rational l_low, l_high;
+    auto sign = u.is_pos() ? 1 : -1;
+    auto p = sign * floor(abs(u / x2));
+    auto p_ = p + sign;
+    lp_assert(p * x2 <= u && u <= p_ * x2 || p * x2 >= u && u >= p_ * x2);
+    std::cout << "u = " << u << ", v = " << v << std::endl;
+    std::cout << "p = " << p << ", p_ = " << p_ << std::endl;
+    std::cout << "u - p*x2 = " << u - p * x2 << ", u - p_*x2 = " << u - p_ * x2 << std::endl;
+    delta_0 = (u - p * x2) * t * x1;
+    delta_1 = (u - p_ * x2) * t * x1;
+
+    std::cout << "delta_0 = " << delta_0 << std::endl;
+    std::cout << "delta_1 = " << delta_1 << std::endl;
+
+    return true;
+}
+void test_patching_alpha(const rational &x, const rational &alpha) {
+    std::cout << "\nstart patching x = " << x << ", alpha = " << alpha << "\n";
+    asserts_on_patching(x, alpha);
+    rational delta_0, delta_1;
+    bool r = get_patching_deltas(x, alpha, delta_0, delta_1);
+    if (r) {
+        lp_assert(delta_0 * delta_1 <= 0);
+
+        lp_assert((x - alpha * delta_0).is_int());
+        lp_assert((x - alpha * delta_1).is_int());
+
+        std::cout << "success\n";
+        // std::cout << "delta_minus = " << delta_minus << ", delta_1 = " << delta_1 << "\n";
+        // std::cout << "x + alpha*delta_minus = " << x + alpha * delta_minus << "\n";
+        // std::cout << "x + alpha*delta_1 = " << x + alpha * delta_1 << "\n";
+    }
+}
+
+void find_a1_x1_x2_and_fix_a2(int &x1, int &x2, int &a1, int &a2) {
+    x2 = (rand() % a2) + (int)(a2 / 3);
+    auto g = gcd(rational(a2), rational(x2));
+    a2 *= (x2 / numerator(g).get_uint64());
+    lp_assert(rational(a2, x2).is_int());
+    do {
+        x1 = rand() % (unsigned)x2 + 1;
+    } while (!coprime(x1, x2));
+
+    do {
+        a1 = rand() % (unsigned)a2 + 1;
+    } while (!coprime(a1, a2));
+    x1 *= (rand() % 2 == 0 ? 1 : -1);
+    a1 *= (rand() % 2 == 0 ? 1 : -1);
+}
+
+void test_patching() {
+    srand(1);
+    // repeat the test 100 times
+
+    int range = 40;
+    for (int i = 0; i < 100; i++) {
+        int a1;
+        int a2 = std::max((int)rand() % range, (int)range / 3);
+
+        int x1, x2;
+        find_a1_x1_x2_and_fix_a2(x1, x2, a1, a2);
+
+        test_patching_alpha(rational(x1, x2), rational(a1, a2));
+    }
+}
\ No newline at end of file