diff --git a/src/math/lp/bound_analyzer_on_row.h b/src/math/lp/bound_analyzer_on_row.h
index 341b18e67..cf0ff6436 100644
--- a/src/math/lp/bound_analyzer_on_row.h
+++ b/src/math/lp/bound_analyzer_on_row.h
@@ -275,22 +275,27 @@ namespace lp {
 
         void limit_j(unsigned bound_j, const mpq& u, bool coeff_before_j_is_pos, bool is_lower_bound, bool strict) {
             auto* lar = &m_bp.lp();
-            int bound_sign = (is_lower_bound ? 1 : -1);
-            int j_sign = (coeff_before_j_is_pos ? 1 : -1) * bound_sign;
-                
-            u_dependency* dep = nullptr;
-            for (auto const& r : m_row) {
-                unsigned j = r.var();
-                if (j == bound_j)
-                    continue;
-                mpq const& a = r.coeff();
-                int a_sign = is_pos(a) ? 1 : -1;
-                int sign = j_sign * a_sign;
-                u_dependency* witness = sign > 0 ? lar->get_column_upper_bound_witness(j) : lar->get_column_lower_bound_witness(j);
-                dep = lar->join_deps(dep, witness);
-            }
-                
-            m_bp.add_bound(u, bound_j, is_lower_bound, strict, dep);
+            const auto& row = this->m_row;
+            auto explain = [row, bound_j, coeff_before_j_is_pos, is_lower_bound, strict, lar]() {
+                (void) strict;
+                TRACE("bound_analyzer", tout << "explain_bound_on_var_on_coeff, bound_j = " << bound_j << ", coeff_before_j_is_pos = " << coeff_before_j_is_pos << ", is_lower_bound = " << is_lower_bound << ", strict = " << strict << "\n";);
+                int bound_sign = (is_lower_bound ? 1 : -1);
+                int j_sign = (coeff_before_j_is_pos ? 1 : -1) * bound_sign;
+
+                u_dependency* ret = nullptr;
+                for (auto const& r : row) {
+                    unsigned j = r.var();
+                    if (j == bound_j)
+                        continue;
+                    mpq const& a = r.coeff();
+                    int a_sign = is_pos(a) ? 1 : -1;
+                    int sign = j_sign * a_sign;
+                    u_dependency* witness = sign > 0 ? lar->get_column_upper_bound_witness(j) : lar->get_column_lower_bound_witness(j);
+                    ret = lar->join_deps(ret, witness);
+                }
+                return ret;
+            };
+            m_bp.add_bound(u, bound_j, is_lower_bound, strict, explain);
         }
 
         void advance_u(unsigned j) {
diff --git a/src/math/lp/dioph_eq.cpp b/src/math/lp/dioph_eq.cpp
index cbd5f5004..a840b44ac 100644
--- a/src/math/lp/dioph_eq.cpp
+++ b/src/math/lp/dioph_eq.cpp
@@ -176,8 +176,6 @@ namespace lp {
         }
     };
     class dioph_eq::imp {
-        int m_n_of_lemmas = 0;
-        
         // This class represents a term with an added constant number c, in form sum
         // {x_i*a_i} + c.
         class term_o : public lar_term {
@@ -495,15 +493,6 @@ namespace lp {
         unsigned m_conflict_index = -1;  // the row index of the conflict
         unsigned m_max_of_branching_iterations = 0;
         unsigned m_number_of_branching_calls;
-        struct prop_bound {
-            mpq m_bound;
-            unsigned m_j;
-            bool m_is_low;
-            bool m_strict;
-            u_dependency* m_dep;
-        };
-
-        std_vector<prop_bound> m_prop_bounds;
         struct branch {
             unsigned m_j = UINT_MAX;
             mpq m_rs;
@@ -1218,7 +1207,6 @@ namespace lp {
                 unsigned j = p.var();
                 if (j == k)
                     continue;
-                
                 m_espace.add(p.coeff() * coeff, j);
                 // do we need to add j to the queue?
                 if (m_espace.has(j) && can_substitute(j))
@@ -1346,7 +1334,6 @@ namespace lp {
             return ret;
         }
 
-        
         const unsigned sub_index(unsigned k) const {
             return m_k2s[k];
         }
@@ -1410,8 +1397,7 @@ namespace lp {
                 if (j >= lra.column_count() ||
                     !lra.column_has_term(j) ||
                     lra.column_is_free(j) ||
-                    !lia.column_is_int(j) ||
-                    !term_has_int_inf_vars(j)) {
+                    is_fixed(j) || !lia.column_is_int(j)) {
                     cleanup.push_back(j);
                     continue;
                 }
@@ -1433,7 +1419,7 @@ namespace lp {
             for (unsigned j : sorted_changed_terms) {
                 m_changed_terms.remove(j); 
                 
-                r = tighten_bounds_for_term_column_gcd(j);
+                r = tighten_bounds_for_term_column(j);
                 if (r != lia_move::undef) {
                     break;
                 }                
@@ -1492,62 +1478,6 @@ namespace lp {
             return m_var_register.external_to_local(j);
         }
 
-        std::ostream& print_prop_bound(const prop_bound & b, std::ostream& out) {
-            out << "low:" << b.m_is_low << ", bound:" << b.m_bound << "\n";
-            out << "deps:\n";
-            lra.print_explanation(out, lra.flatten(b.m_dep));
-            return out;
-        }
-        
-        template <typename T>        
-        lia_move propagate_bounds_on_espace(const mpq& sum_of_fixed, const T& meta_term) {
-            // change m_espace indices from local to lar_solver: have to restore for clear to work
-            if (has_fresh_var(m_espace)) return lia_move::undef;
-            for (auto & p: m_espace.m_data) {
-                p.m_j = local_to_lar_solver(p.m_j);        
-            }
-            m_prop_bounds.clear();
-            bound_analyzer_on_row<term_with_index, dioph_eq::imp>::analyze_row(m_espace, impq(-sum_of_fixed), *this);
-            //restore m_espace to local variables
-            for (auto & p: m_espace.m_data) p.m_j = lar_solver_to_local(p.m_j);
-            if (m_prop_bounds.size() == 0) {
-                TRACE("dio", tout <<"no new bounds\n";);
-                return lia_move::undef;
-            }
-            TRACE("dio", tout << "prop_bounds:\n"; for (auto& pb: m_prop_bounds) print_prop_bound(pb, tout););
-            u_dependency * dep  = explain_fixed_in_meta_term(meta_term);
-            bool change = false;
-            for (auto& pb: m_prop_bounds) {
-                pb.m_dep = lra.join_deps(pb.m_dep, dep);                
-                if (update_bound(pb)) {
-                    change = true;
-                    TRACE("dio",
-                          tout << "change after update_bound pb.m_j:" << pb.m_j << "\n";
-                          lra.print_explanation(tout, lra.flatten(pb.m_dep)));
-                }
-            }
-            if (!change)
-                return lia_move::undef;
-            auto st = lra.find_feasible_solution();
-            if (st == lp_status::INFEASIBLE) {
-                lra.get_infeasibility_explanation(m_infeas_explanation);
-                TRACE("dio", tout << "inf explanation:\n"; lra.print_explanation(tout, m_infeas_explanation););
-                return lia_move::conflict;
-            }
-            TRACE("dio", tout << "lra is feasible\n";);                        
-            return lia_move::undef;
-        }
-
-        // h is the entry that is ready to be moved to S
-        lia_move tighten_bounds_on_entry(unsigned h) {
-            SASSERT(entry_invariant(h));
-            protected_queue q;
-            copy_row_to_espace(h);
-            remove_fresh_from_espace();
-            return propagate_bounds_on_espace(m_sum_of_fixed[h], m_l_matrix[h]);
-        }
-
-
         void process_fixed_in_espace() {
             std_vector<unsigned> fixed_vars;
             for (const auto & p: m_espace) {
@@ -1566,10 +1496,11 @@ namespace lp {
   we have x+y = 7t - 1 <= 8, where t is a term. Then we have 7t <= 9, or t <= 9/7, or we can enforce t <= floor(9/7) = 1.
   Then x + y = 7*1 - 1 <= 6: the bound is strenthgened
 */
-        lia_move tighten_bounds_for_term_column_gcd(unsigned j) {
+        lia_move tighten_bounds_for_term_column(unsigned j) {
             term_o term_to_tighten = lra.get_term(j);  // copy the term aside
             
-            SASSERT(get_extended_term_value(term_to_tighten).is_zero() && all_vars_are_int(term_to_tighten));
+            if (!all_vars_are_int(term_to_tighten))
+                return lia_move::undef;
             // q is the queue of variables that can be substituted in term_to_tighten
             protected_queue q;
             TRACE("dio", tout << "j:" << j << " , intitial term t: "; print_lar_term_L(term_to_tighten, tout) << std::endl;
@@ -1612,15 +1543,7 @@ namespace lp {
             if (tighten_bounds_for_non_trivial_gcd(g, j, false) != lia_move::undef) {
                 return lia_move::conflict;
             }
-
-            // try more bound progagation
-            if (is_fixed(j)) {
-                m_c -= lra.get_lower_bound(j).x;
-            }
-            else {
-                m_espace.add(-mpq(1), lar_solver_to_local(j));
-            }
-            return propagate_bounds_on_espace(m_c, m_lspace);
+            return lia_move::undef;
         }
 
         bool should_report_branch() const {
@@ -1744,81 +1667,6 @@ namespace lp {
             }
         }
 
-        void add_bound(mpq const& bound, unsigned j, bool is_low, bool strict, u_dependency* dep) {
-            m_prop_bounds.push_back({bound, j, is_low, strict, dep});
-        }
-
-        lconstraint_kind get_bound_kind(const prop_bound& pb) const {
-            if (!pb.m_is_low) {
-                return pb.m_strict ? lp::LT : lp::LE; 
-            }
-            else {
-                return pb.m_strict ? lp::GT : lp::GE;
-            }
-        }
-
-        bool cut(const prop_bound& pb) const {            
-            const auto & v = lra.get_column_value(pb.m_j);
-            switch (get_bound_kind(pb)) {
-            case LT:
-                return v >= pb.m_bound;
-            case LE:
-                return v > pb.m_bound;
-            case GT:
-                return v <= pb.m_bound;
-            case GE:
-                return v < pb.m_bound;
-            default:
-                UNREACHABLE();
-            }
-            return false;
-        }
-
-        std::ostream& print_int_inf_vars(std::ostream& out) const {
-            out << "Integer infeasible variables:\n";
-            for (unsigned j = 0; j < lra.column_count(); j++) {
-                if (!lia.column_is_int_inf(j)) continue;
-                out << "x" << j << " = " << lra.get_column_value(j).x << " bounds:\n";
-                if (lra.column_has_lower_bound(j))
-                    out << lra.get_lower_bound(j).x << " <= ";
-                else 
-                    out << "-oo" << " <= ";
-
-                out << "x" << j;
-                if (lra.column_has_upper_bound(j)) 
-                    out << " <= " << lra.get_upper_bound(j).x;
-                else 
-                    out << " <= " << "oo";    
-                out << "\n";
-                
-            }
-            return out;
-        }
-
-        
-        // return true iff the column bound has been changed
-        bool update_bound(const prop_bound& pb) {            
-            TRACE("dio", tout << "pb: " << "x" << pb.m_j << ", low:" << pb.m_is_low << " , strict:" << pb.m_strict << " , bound:" << pb.m_bound << "\n"; lra.print_column_info(pb.m_j, tout, true););
-            bool r = lra.update_column_type_and_bound(pb.m_j, get_bound_kind(pb), pb.m_bound, pb.m_dep);
-            CTRACE("dio", r, tout << "change in lar_solver: "; tout << "was updating with " << (pb.m_is_low? "lower" : "upper") << " bound " << pb.m_bound << "\n";
-                   tout << "the column became:\n";
-                   lra.print_column_info(pb.m_j, tout);tout <<"return " << r << "\n";);
-            if (lra.settings().dump_bound_lemmas()) {
-                std::string lemma_name = "lemma" + std::to_string(m_n_of_lemmas++); 
-                lra.write_bound_lemma_to_file(pb.m_j, pb.m_is_low, lemma_name, std::string( __FILE__)+ ","+ std::to_string(__LINE__));
-            }
-            return r;
-        }
-        
-        bool row_has_int_inf(unsigned ei) {
-            for (const auto& p: m_e_matrix.m_rows[ei]) {
-                unsigned j = local_to_lar_solver(p.var());
-                if (lia.column_is_int_inf(j))
-                    return true;
-            }
-            return false;
-        }
-
         
         // m_espace contains the coefficients of the term
         // m_c contains the fixed part of the term
@@ -1845,10 +1693,9 @@ namespace lp {
                     TRACE("dio", tout << "no improvement in the bound\n";);
                 }
             }
-            
             return lia_move::undef;
         }
-
+        unsigned m_n_of_lemmas = 0;
         // returns true only on a conflict
         bool tighten_bound_kind(const mpq& g, unsigned j, const mpq& ub, bool upper) {
             // ub = (upper_bound(j) - m_c)/g.
@@ -1949,13 +1796,9 @@ namespace lp {
                 if (m_conflict_index != UINT_MAX) {
                     lra.stats().m_dio_rewrite_conflicts++;
                 } 
-                else {
-                    lra.stats().m_dio_bound_propagation_conflicts++;
-                }
                 return lia_move::conflict;                
             }
-            TRACE("dio", print_int_inf_vars(tout) << "\n"; 
-                  print_S(tout));
+            TRACE("dio", print_S(tout));
 
             return lia_move::undef;
         }
@@ -2482,7 +2325,7 @@ namespace lp {
             }
             return true;
         }
-        
+
         term_o term_to_lar_solver(const term_o& eterm) const {
             term_o ret;
             for (const auto& p : eterm) {
@@ -2785,14 +2628,6 @@ namespace lp {
             SASSERT(h == f_vector[ih]);
             if (min_ahk.is_one()) {
                 TRACE("dioph_eq", tout << "push to S:\n"; print_entry(h, tout););
-                SASSERT(m_new_fixed_columns.size() == 0);
-                if (tighten_bounds_on_entry(h) == lia_move::conflict){
-                    TRACE("dio", tout << "conflict\n";);
-                    return lia_move::conflict;
-                }
-                if (m_new_fixed_columns.size()) {
-                    return lia_move::undef;
-                }
                 move_entry_from_f_to_s(kh, h);
                 eliminate_var_in_f(h, kh, kh_sign);
                 if (ih != f_vector.size() - 1) {
diff --git a/src/math/lp/implied_bound.h b/src/math/lp/implied_bound.h
index ed6a81692..6d7fc54b9 100644
--- a/src/math/lp/implied_bound.h
+++ b/src/math/lp/implied_bound.h
@@ -33,11 +33,11 @@ class implied_bound {
     bool m_is_lower_bound;
     bool m_strict;
     private:
-    u_dependency* m_dep = nullptr;
+    std::function<u_dependency*()> m_explain_bound = nullptr;
     public:
     // s is expected to be the pointer to lp_bound_propagator.
-    u_dependency* dep() const { return m_dep; }
-    void set_dep(u_dependency* dep) { m_dep = dep; }
+    u_dependency* explain_implied() const { return m_explain_bound(); }
+    void set_explain(std::function<u_dependency*()> f) { m_explain_bound = f; }
     lconstraint_kind kind() const {
         lconstraint_kind k = m_is_lower_bound? GE : LE;
         if (m_strict)
@@ -49,12 +49,12 @@ class implied_bound {
                   unsigned j,
                   bool is_lower_bound,
                   bool is_strict,
-                  u_dependency* dep):
+                  std::function<u_dependency*()> get_dep):
         m_bound(a),
         m_j(j),
         m_is_lower_bound(is_lower_bound),
         m_strict(is_strict),
-        m_dep(dep) {
+        m_explain_bound(get_dep) {
     }
 };
 }
diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp
index 103b974ec..8c6d8065d 100644
--- a/src/math/lp/int_solver.cpp
+++ b/src/math/lp/int_solver.cpp
@@ -44,16 +44,7 @@ namespace lp {
         unsigned            m_dioph_eq_period;
         dioph_eq            m_dio;  
         int_gcd_test        m_gcd;
-        struct prop_bound {
-            mpq m_bound;
-            unsigned m_j;
-            bool m_is_low;
-            bool m_strict;
-            u_dependency* m_dep;
-        };
 
-        std_vector<prop_bound> m_prop_bounds;
-        
         bool column_is_int_inf(unsigned j) const {
             return lra.column_is_int(j) && (!lia.value_is_int(j));
         }
@@ -226,50 +217,74 @@ namespace lp {
         // consider multi-round bound tightnening as well and deal with divergence issues.
 
         unsigned m_bounds_refine_depth = 0;
-       
-        void add_bound(mpq const& bound, unsigned j, bool is_low, bool strict, u_dependency* dep) {
-            m_prop_bounds.push_back({bound, j, is_low, strict, dep});
-        }        
 
-        lconstraint_kind get_bound_kind(bool upper, bool is_strict) {
-            if (upper) {
-                return is_strict ? lp::LT : lp::LE; 
-            }
-            else {
-                return is_strict ? lp::GT : lp::GE;
-            }
-        }
-
-        void add_propagated_bounds() {
-            for (const auto& pb: m_prop_bounds) {
-                lconstraint_kind kind = get_bound_kind(!pb.m_is_low, pb.m_strict);            
-                lra.update_column_type_and_bound(pb.m_j, kind, pb.m_bound, pb.m_dep);
-            }
-        }
-        
         lia_move tighten_bounds() {
-            if (m_bounds_refine_depth > 1)
-                return lia_move::undef;
-            bool bounds_refined = false;
 
-            unsigned start = random();
-            for (unsigned k = 0; k < lra.m_touched_rows.size(); k++) {
-                m_prop_bounds.clear();
-                unsigned i = (start + k) % lra.m_touched_rows.size();
-                bound_analyzer_on_row<row_strip<mpq>, int_solver::imp>::analyze_row(lra.get_row(i), impq(0), *this);
-                if (m_prop_bounds.size() == 0) continue;
-                add_propagated_bounds();
+            if (m_bounds_refine_depth > 10)
+                return lia_move::undef;
+
+            struct bound_consumer {
+                imp& i;
+                bound_consumer(imp& i) : i(i) {}
+                lar_solver& lp() { return i.lra; }
+                bool bound_is_interesting(unsigned vi, lp::lconstraint_kind kind, const rational& bval) const { return true; }
+                bool add_eq(lpvar u, lpvar v, lp::explanation const& e, bool is_fixed) { return false; }
+            };
+            bound_consumer bc(*this);
+            std_vector<implied_bound> ibounds;
+            lp_bound_propagator<bound_consumer> bp(bc, ibounds);
+
+            auto set_conflict = [&](u_dependency * d1, u_dependency * d2) {
+                ++settings().stats().m_bounds_tightening_conflicts;
+                for (auto e : lra.flatten(d1))
+                    m_ex->push_back(e);
+                for (auto e : lra.flatten(d2))
+                    m_ex->push_back(e);
+            };
+
+            bool bounds_refined = false;
+            auto refine_bound = [&](implied_bound const& ib) {
+                unsigned j = ib.m_j;
+                rational bound = ib.m_bound;
+                if (ib.m_is_lower_bound) {
+                    if (lra.column_is_int(j))
+                        bound = ceil(bound);
+                    if (lra.column_has_lower_bound(j) && lra.column_lower_bound(j) >= bound)
+                        return l_undef;
+                    auto d = ib.explain_implied();
+                    if (lra.column_has_upper_bound(j) && lra.column_upper_bound(j) < bound) {
+                        set_conflict(d, lra.get_column_upper_bound_witness(j));
+                        return l_false;
+                    }
+                    lra.update_column_type_and_bound(j, ib.m_strict ? lconstraint_kind::GT : lconstraint_kind::GE, bound, d);                    
+                } 
+                else {
+                    if (lra.column_is_int(j))
+                        bound = floor(bound);
+                    if (lra.column_has_upper_bound(j) && lra.column_upper_bound(j) <= bound)
+                        return l_undef;
+                    auto d = ib.explain_implied();
+                    if (lra.column_has_lower_bound(j) && lra.column_lower_bound(j) > bound) {                        
+                        set_conflict(d, lra.get_column_lower_bound_witness(j));
+                        return l_false;
+                    }
+                    lra.update_column_type_and_bound(j, ib.m_strict ? lconstraint_kind::LT : lconstraint_kind::LE, bound, d);
+                }
+                ++settings().stats().m_bounds_tightenings;
                 bounds_refined = true;
-                lp_status st = lra.find_feasible_solution();
-                if ((int)st >= (int)lp::lp_status::FEASIBLE) {
-                    continue;
-                }
-                if (st == lp_status::CANCELLED) {
-                    return lia_move::undef;
-                }
-                SASSERT(st == lp_status::INFEASIBLE);
-                lra.get_infeasibility_explanation(*m_ex);
-                return lia_move::conflict;
+                return l_true;
+            };
+
+            for (unsigned row_index = 0; row_index < lra.row_count(); ++row_index) {
+                bp.init();
+                bound_analyzer_on_row<row_strip<mpq>, lp_bound_propagator<bound_consumer>>::analyze_row(
+                    lra.A_r().m_rows[row_index],
+                    zero_of_type<numeric_pair<mpq>>(),                    
+                    bp);
+
+                for (auto const& ib : ibounds)
+                    if (l_false == refine_bound(ib))
+                        return lia_move::conflict;
             }
             
             if (bounds_refined) {
@@ -277,14 +292,9 @@ namespace lp {
                 ++m_bounds_refine_depth;
             }
 
-            return lia_move::undef;
+            return bounds_refined  ? lia_move::continue_with_check : lia_move::undef;
         }
-        bool should_tighten_bounds() {
-            return false && m_number_of_calls % 4 == 0;
-        }
-     // needed for the template bound_analyzer_on_row.h
-        const lar_solver& lp() const { return lra; }
-        lar_solver& lp() {return lra;}
+
         
         lia_move check(lp::explanation * e) {
             SASSERT(lra.ax_is_correct());
@@ -311,10 +321,10 @@ namespace lp {
             if (r == lia_move::undef) r = patch_basic_columns();
             if (r == lia_move::undef && should_find_cube()) r = int_cube(lia)();
             if (r == lia_move::undef) lra.move_non_basic_columns_to_bounds();
+            // if (r == lia_move::undef) r = tighten_bounds();
             if (r == lia_move::undef && should_hnf_cut()) r = hnf_cut();
             if (r == lia_move::undef && should_gomory_cut()) r = gomory(lia).get_gomory_cuts(2);
             if (r == lia_move::undef && should_solve_dioph_eq()) r = solve_dioph_eq();
-            if (r == lia_move::undef && should_tighten_bounds()) r = tighten_bounds();
             if (r == lia_move::undef) r = int_branch(lia)();
             if (settings().get_cancel_flag()) r = lia_move::undef;        
             return r;
@@ -330,7 +340,7 @@ namespace lp {
 
         bool current_solution_is_inf_on_cut() const {
             SASSERT(cut_indices_are_columns());
-            const auto & x = lrac.m_r_x;
+            const auto & x = lrac.r_x();
             impq v = m_t.apply(x);
             mpq sign = m_upper ? one_of_type<mpq>()  : -one_of_type<mpq>();
             CTRACE("current_solution_is_inf_on_cut", v * sign <= impq(m_k) * sign,
@@ -672,7 +682,7 @@ namespace lp {
     }
 
     const impq & int_solver::get_value(unsigned j) const {
-        return lrac.m_r_x[j];
+        return lrac.r_x(j);
     }
 
     std::ostream& int_solver::display_column(std::ostream & out, unsigned j) const {
diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h
index 070fcd285..8e727bce9 100644
--- a/src/math/lp/lar_core_solver.h
+++ b/src/math/lp/lar_core_solver.h
@@ -22,19 +22,21 @@ class lar_core_solver  {
     vector<numeric_pair<mpq>> m_right_sides_dummy;
     vector<mpq> m_costs_dummy;
     stacked_value<simplex_strategy_enum> m_stacked_simplex_strategy;
+    vector<impq> m_r_x;  // the solution
+    vector<impq> m_backup_x;
 
 public:
     
     stacked_vector<column_type> m_column_types;
     // r - solver fields, for rational numbers
-    vector<numeric_pair<mpq>> m_r_x; // the solution
+
     stacked_vector<numeric_pair<mpq>> m_r_lower_bounds;
     stacked_vector<numeric_pair<mpq>> m_r_upper_bounds;
     static_matrix<mpq, numeric_pair<mpq>> m_r_A;
     stacked_vector<unsigned> m_r_pushed_basis;
     vector<unsigned>         m_r_basis;
     vector<unsigned>         m_r_nbasis;
-    std_vector<int>              m_r_heading;
+    std_vector<int>          m_r_heading;
     
 
     lp_primal_core_solver<mpq, numeric_pair<mpq>> m_r_solver; // solver in rational numbers
@@ -71,9 +73,23 @@ public:
         
         m_r_solver.print_column_bound_info(m_r_solver.m_basis[row_index], out);        
     }
-    
-    
+
+
     void prefix_r();
+    
+    // access to x:
+
+    void backup_x() { m_backup_x = m_r_x; }
+
+    void restore_x() {
+        m_r_x = m_backup_x;
+        m_r_x.reserve(m_m());
+    }
+
+    vector<impq> const& r_x() const { return m_r_x; }
+    impq& r_x(unsigned j) { return m_r_x[j]; }
+    impq const& r_x(unsigned j) const { return m_r_x[j]; }
+    void resize_x(unsigned n) { m_r_x.resize(n); }
 
     unsigned m_m() const { return m_r_A.row_count();  }
 
diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h
index e0ffed16e..8fc98db1d 100644
--- a/src/math/lp/lar_core_solver_def.h
+++ b/src/math/lp/lar_core_solver_def.h
@@ -76,7 +76,7 @@ void lar_core_solver::fill_not_improvable_zero_sum() {
 
 unsigned lar_core_solver::get_number_of_non_ints() const {
     unsigned n = 0;
-    for (auto & x : m_r_solver.m_x) 
+    for (auto & x : r_x()) 
         if (!x.is_int())
             n++;    
     return n;
diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index bff371d0d..c0eac9dc3 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -45,7 +45,7 @@ namespace lp {
     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());
         lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
-        lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_x.size());
+        lp_assert(A_r().column_count() == m_mpq_lar_core_solver.r_x().size());
         return true;
     }
 
@@ -65,8 +65,8 @@ namespace lp {
     }
 
     std::ostream& lar_solver::print_values(std::ostream& out) const {
-        for (unsigned i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++) {
-            const numeric_pair<mpq>& rp = m_mpq_lar_core_solver.m_r_x[i];
+        for (unsigned i = 0; i < m_mpq_lar_core_solver.r_x().size(); i++) {
+            const numeric_pair<mpq>& rp = m_mpq_lar_core_solver.r_x(i);
             out << this->get_variable_name(i) << " -> " << rp << "\n";
         }
         return out;
@@ -264,7 +264,7 @@ namespace lp {
             return false;
         }
         else {
-            term_max = term.apply(m_mpq_lar_core_solver.m_r_x);
+            term_max = term.apply(m_mpq_lar_core_solver.r_x());
             return true;
         }
     }
@@ -418,7 +418,7 @@ namespace lp {
 
     bool lar_solver::move_non_basic_column_to_bounds(unsigned j) {
         auto& lcs = m_mpq_lar_core_solver;
-        auto& val = lcs.m_r_x[j];
+        auto& val = lcs.r_x(j);
         switch (lcs.m_column_types()[j]) {
         case column_type::boxed: {
             const auto& l = lcs.m_r_lower_bounds()[j];
@@ -458,7 +458,7 @@ namespace lp {
 
     void lar_solver::set_value_for_nbasic_column(unsigned j, const impq& new_val) {
         lp_assert(!is_base(j));
-        auto& x = m_mpq_lar_core_solver.m_r_x[j];
+        auto& x = m_mpq_lar_core_solver.r_x(j);
         auto delta = new_val - x;
         x = new_val;
         TRACE("lar_solver_feas", tout << "setting " << j << " to "
@@ -505,7 +505,7 @@ namespace lp {
         if (column_has_term(j)) {
             return * m_columns[j].term();
         }
-        if (j < m_mpq_lar_core_solver.m_r_x.size()) {
+        if (j < m_mpq_lar_core_solver.r_x().size()) {
             lar_term r;
             r.add_monomial(one_of_type<mpq>(), j);
             return r;
@@ -519,41 +519,44 @@ namespace lp {
         SASSERT(m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis());
         lar_term term = get_term_to_maximize(j);
         if (term.is_empty()) return lp_status::UNBOUNDED;
-        impq prev_value = term.apply(m_mpq_lar_core_solver.m_r_x);
-        auto backup = m_mpq_lar_core_solver.m_r_x;
+        m_mpq_lar_core_solver.backup_x();
+        impq prev_value = term.apply(m_mpq_lar_core_solver.r_x());
+        auto restore = [&]() {
+            m_mpq_lar_core_solver.restore_x();
+        };
         if (!maximize_term_on_feasible_r_solver(term, term_max, nullptr)) {
-            m_mpq_lar_core_solver.m_r_x = backup;
+            restore();
             return lp_status::UNBOUNDED;
         }
 
         impq opt_val = term_max;
 
         bool change = false;
-        for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_x.size(); j++) {
+        for (unsigned j = 0; j < m_mpq_lar_core_solver.r_x().size(); j++) {
             if (!column_is_int(j))
                 continue;
             if (column_value_is_integer(j))
                 continue;
             if (m_int_solver->is_base(j)) {
                 if (!remove_from_basis(j)) { // consider a special version of remove_from_basis that would not remove inf_int columns
-                    m_mpq_lar_core_solver.m_r_x = backup;
+                    restore();
                     term_max = prev_value;
                     return lp_status::FEASIBLE; // it should not happen
                 }
             }
             if (!column_value_is_integer(j)) {
                 term_max = prev_value;
-                m_mpq_lar_core_solver.m_r_x = backup;
+                restore();
                 return lp_status::FEASIBLE;
             }
             change = true;
         }
         if (change) {
-            term_max = term.apply(m_mpq_lar_core_solver.m_r_x);
+            term_max = term.apply(m_mpq_lar_core_solver.r_x());
         }
         if (term_max < prev_value) {
             term_max = prev_value;
-            m_mpq_lar_core_solver.m_r_x = backup;
+            restore();
         }
         TRACE("lar_solver", print_values(tout););
         if (term_max == opt_val) {
@@ -823,7 +826,7 @@ namespace lp {
     bool lar_solver::row_is_correct(unsigned i) const {
         numeric_pair<mpq> r = zero_of_type<numeric_pair<mpq>>();
         for (const auto& c : A_r().m_rows[i]) {
-            r += c.coeff() * m_mpq_lar_core_solver.m_r_x[c.var()];
+            r += c.coeff() * m_mpq_lar_core_solver.r_x(c.var());
         }
         CTRACE("lar_solver", !is_zero(r), tout << "row = " << i << ", j = " << m_mpq_lar_core_solver.m_r_basis[i] << "\n";
                print_row(A_r().m_rows[i], tout);  tout << " = " << r << "\n");
@@ -932,7 +935,7 @@ namespace lp {
         unsigned bj = m_mpq_lar_core_solver.m_r_solver.m_basis[i];
         for (const auto& c : A_r().m_rows[i]) {
             if (c.var() == bj) continue;
-            const auto& x = m_mpq_lar_core_solver.m_r_x[c.var()];
+            const auto& x = m_mpq_lar_core_solver.r_x(c.var());
             if (!is_zero(x))
                 r -= c.coeff() * x;
         }
@@ -1170,7 +1173,7 @@ namespace lp {
         if (!init_model())
             return;
 
-        unsigned n = m_mpq_lar_core_solver.m_r_x.size();
+        unsigned n = m_mpq_lar_core_solver.r_x().size();
 
         for (unsigned j = 0; j < n; j++) 
             variable_values[j] = get_value(j);
@@ -1180,11 +1183,15 @@ namespace lp {
     }
 
     bool lar_solver::init_model() const {
+        auto& rslv = m_mpq_lar_core_solver.m_r_solver;
+        lp_assert(A_r().column_count() == rslv.m_costs.size());
+        lp_assert(A_r().column_count() == m_mpq_lar_core_solver.r_x().size());
+        lp_assert(A_r().column_count() == rslv.m_d.size());
         CTRACE("lar_solver_model",!m_columns_with_changed_bounds.empty(), tout << "non-empty changed bounds\n");
         TRACE("lar_solver_model", tout << get_status() << "\n");
         auto status = get_status();
         SASSERT((status != lp_status::OPTIMAL && status != lp_status::FEASIBLE) 
-            || m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis());
+            || rslv.calc_current_x_is_feasible_include_non_basis());
         if (status != lp_status::OPTIMAL && status != lp_status::FEASIBLE)
             return false;
         if (!m_columns_with_changed_bounds.empty())
@@ -1192,12 +1199,12 @@ namespace lp {
 
         m_delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1));
         unsigned j;
-        unsigned n = m_mpq_lar_core_solver.m_r_x.size();
+        unsigned n = m_mpq_lar_core_solver.r_x().size();
         do {
             m_set_of_different_pairs.clear();
             m_set_of_different_singles.clear();
             for (j = 0; j < n; j++) {
-                const numeric_pair<mpq>& rp = m_mpq_lar_core_solver.m_r_x[j];
+                const numeric_pair<mpq>& rp = m_mpq_lar_core_solver.r_x(j);
                 mpq x = rp.x + m_delta * rp.y;
                 m_set_of_different_pairs.insert(rp);
                 m_set_of_different_singles.insert(x);
@@ -1213,8 +1220,8 @@ namespace lp {
 
     void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map<lpvar, mpq>& variable_values) const {
         mpq delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1));
-        for (unsigned i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++) {
-            const impq& rp = m_mpq_lar_core_solver.m_r_x[i];
+        for (unsigned i = 0; i < m_mpq_lar_core_solver.r_x().size(); i++) {
+            const impq& rp = m_mpq_lar_core_solver.r_x(i);
             variable_values[i] = rp.x + delta * rp.y;
         }
     }
@@ -1229,7 +1236,7 @@ namespace lp {
     void lar_solver::get_rid_of_inf_eps() {
         bool y_is_zero = true;
         for (unsigned j = 0; j < number_of_vars(); j++) {
-            if (!m_mpq_lar_core_solver.m_r_x[j].y.is_zero()) {
+            if (!m_mpq_lar_core_solver.r_x(j).y.is_zero()) {
                 y_is_zero = false;
                 break;
             }
@@ -1238,7 +1245,7 @@ namespace lp {
             return;
         mpq delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1));
         for (unsigned j = 0; j < number_of_vars(); j++) {
-            auto& v = m_mpq_lar_core_solver.m_r_x[j];
+            auto& v = m_mpq_lar_core_solver.r_x(j);
             if (!v.y.is_zero()) {
                 v = impq(v.x + delta * v.y);
                 TRACE("lar_solver_feas", tout << "x[" << j << "] = " << v << "\n";);
@@ -1448,7 +1455,7 @@ namespace lp {
     void lar_solver::remove_last_column_from_tableau() {
         auto& rslv = m_mpq_lar_core_solver.m_r_solver;
         unsigned j = A_r().column_count() - 1;
-        lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
+        lp_assert(A_r().column_count() == rslv.m_costs.size());
         if (column_represents_row_in_tableau(j)) {
             remove_last_row_and_column_from_tableau(j);
             if (rslv.m_basis_heading[j] < 0)
@@ -1457,13 +1464,15 @@ namespace lp {
         else {
             remove_last_column_from_A();
         }
-        rslv.m_x.pop_back();
+        m_mpq_lar_core_solver.resize_x(A_r().column_count());
         rslv.m_d.pop_back();
         rslv.m_costs.pop_back();
 
         remove_last_column_from_basis_tableau(j);
         lp_assert(m_mpq_lar_core_solver.r_basis_is_OK());
-        lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
+        lp_assert(A_r().column_count() == rslv.m_costs.size());
+        lp_assert(A_r().column_count() == m_mpq_lar_core_solver.r_x().size());
+        lp_assert(A_r().column_count() == rslv.m_d.size());
     }
 
 
@@ -1611,14 +1620,15 @@ namespace lp {
         unsigned j = A_r().column_count();
         TRACE("add_var", tout << "j = " << j << std::endl;);
         A_r().add_column();
-        lp_assert(m_mpq_lar_core_solver.m_r_x.size() == j);
+        lp_assert(m_mpq_lar_core_solver.r_x().size() == j);
         //        lp_assert(m_mpq_lar_core_solver.m_r_lower_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j);  // restore later
-        m_mpq_lar_core_solver.m_r_x.resize(j + 1);
+        m_mpq_lar_core_solver.resize_x(j + 1);
+        auto& rslv = m_mpq_lar_core_solver.m_r_solver;
         m_mpq_lar_core_solver.m_r_lower_bounds.increase_size_by_one();
         m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
-        m_mpq_lar_core_solver.m_r_solver.inf_heap_increase_size_by_one();
-        m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
-        m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
+        rslv.inf_heap_increase_size_by_one();
+        rslv.m_costs.resize(j + 1);
+        rslv.m_d.resize(j + 1);
         lp_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method
         if (register_in_basis) {
             A_r().add_row();
@@ -1867,14 +1877,14 @@ namespace lp {
         }
     }
 
-    bool lar_solver::update_column_type_and_bound(unsigned j,
+    void lar_solver::update_column_type_and_bound(unsigned j,
         const mpq& right_side,
         constraint_index constr_index) {
         TRACE("lar_solver_feas", tout << "j = " << j << " was " << (this->column_is_feasible(j)?"feas":"non-feas") << std::endl;);
         m_constraints.activate(constr_index);
         lconstraint_kind kind = m_constraints[constr_index].kind();
         u_dependency* dep = m_constraints[constr_index].dep();
-        return update_column_type_and_bound(j, kind, right_side, dep);
+        update_column_type_and_bound(j, kind, right_side, dep);
     }
 
 
@@ -1956,8 +1966,7 @@ namespace lp {
         }
         ls.add_var_bound(tv, c.kind(), c.rhs());
     }
-    // return true iff a bound has been changed or or the criss-crossed bounds have been found
-    bool lar_solver::update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
+    void lar_solver::update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
         // SASSERT(validate_bound(j, kind, right_side, dep));
         TRACE(
             "lar_solver_feas",
@@ -1973,12 +1982,11 @@ namespace lp {
             });
         bool was_fixed = column_is_fixed(j);
         mpq rs = adjust_bound_for_int(j, kind, right_side);
-        bool r;
         if (column_has_upper_bound(j))
-            r = update_column_type_and_bound_with_ub(j, kind, rs, dep);
+            update_column_type_and_bound_with_ub(j, kind, rs, dep);
         else
-            r = update_column_type_and_bound_with_no_ub(j, kind, rs, dep);
- 
+            update_column_type_and_bound_with_no_ub(j, kind, rs, dep);
+
         if (!was_fixed && column_is_fixed(j) && m_fixed_var_eh)
             m_fixed_var_eh(j);
 
@@ -1989,7 +1997,6 @@ namespace lp {
         TRACE("lar_solver_feas", tout << "j = " << j << " became " << (this->column_is_feasible(j) ? "feas" : "non-feas") << ", and " << (this->column_is_bounded(j) ? "bounded" : "non-bounded") << std::endl;);
         if (m_update_column_bound_callback)
             m_update_column_bound_callback(j);
-        return r;
     }
 
     void lar_solver::insert_to_columns_with_changed_bounds(unsigned j) {
@@ -2014,37 +2021,37 @@ namespace lp {
         return m_constraints.add_term_constraint(j, m_columns[j].term(), kind, rs);
     }
 
-    struct scoped_backup {
+    struct lar_solver::scoped_backup {
         lar_solver& m_s;
         scoped_backup(lar_solver& s) : m_s(s) {
-            m_s.backup_x(); 
+            m_s.m_mpq_lar_core_solver.backup_x();
         }
         ~scoped_backup() {
-            m_s.restore_x();
+            m_s.m_mpq_lar_core_solver.restore_x();
         }
     };
    
-    bool lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
+    void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
         SASSERT(column_has_upper_bound(j));
         if (column_has_lower_bound(j)) {
-            return update_bound_with_ub_lb(j, kind, right_side, dep);
+            update_bound_with_ub_lb(j, kind, right_side, dep);
         }
         else {
-            return update_bound_with_ub_no_lb(j, kind, right_side, dep);
+            update_bound_with_ub_no_lb(j, kind, right_side, dep);
         }
     }
 
-    bool lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
+    void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
         SASSERT(!column_has_upper_bound(j));
         if (column_has_lower_bound(j)) {
-            return update_bound_with_no_ub_lb(j, kind, right_side, dep);
+            update_bound_with_no_ub_lb(j, kind, right_side, dep);
         }
         else {
-            return update_bound_with_no_ub_no_lb(j, kind, right_side, dep);
+            update_bound_with_no_ub_no_lb(j, kind, right_side, dep);
         }
     }
 
-    bool lar_solver::update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
+    void lar_solver::update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
         lp_assert(column_has_lower_bound(j) && column_has_upper_bound(j));
         lp_assert(m_mpq_lar_core_solver.m_column_types[j] == column_type::boxed ||
                   m_mpq_lar_core_solver.m_column_types[j] == column_type::fixed);
@@ -2060,7 +2067,7 @@ namespace lp {
                     set_crossed_bounds_column_and_deps(j, true, dep);
                 }
                 else {
-                    if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return false;
+                    if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return;
                     m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
                     set_upper_bound_witness(j, dep);
                     insert_to_columns_with_changed_bounds(j);
@@ -2075,7 +2082,9 @@ namespace lp {
                 if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
                     set_crossed_bounds_column_and_deps(j, false, dep);
                 } else {
-                    if (low < m_mpq_lar_core_solver.m_r_lower_bounds[j]) return false;
+                    if (low < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                        return;
+                    }
                     m_mpq_lar_core_solver.m_r_lower_bounds[j] = low;
                     set_lower_bound_witness(j, dep);
                     m_mpq_lar_core_solver.m_column_types[j] = (low == m_mpq_lar_core_solver.m_r_upper_bounds[j] ? column_type::fixed : column_type::boxed);
@@ -2108,10 +2117,9 @@ namespace lp {
         numeric_pair<mpq> const& hi = m_mpq_lar_core_solver.m_r_upper_bounds[j];
         if (lo == hi)
             m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
-        return true;
     }
     
-    bool lar_solver::update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
+    void lar_solver::update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
         lp_assert(column_has_lower_bound(j) && !column_has_upper_bound(j));
         lp_assert(m_mpq_lar_core_solver.m_column_types[j] == column_type::lower_bound);
 
@@ -2138,7 +2146,7 @@ namespace lp {
             case GE: {
                 auto low = numeric_pair<mpq>(right_side, y_of_bound);
                 if (low < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
-                    return false;
+                    return;
                 }
                 m_mpq_lar_core_solver.m_r_lower_bounds[j] = low;
                 set_lower_bound_witness(j, dep);
@@ -2163,10 +2171,9 @@ namespace lp {
             default:
                 UNREACHABLE();
         }
-        return true;
     }
    
-    bool lar_solver::update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
+    void lar_solver::update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
         lp_assert(!column_has_lower_bound(j) && column_has_upper_bound(j));
         lp_assert(m_mpq_lar_core_solver.m_column_types[j] == column_type::upper_bound);
         mpq y_of_bound(0);
@@ -2177,7 +2184,7 @@ namespace lp {
         case LE:
         {
             auto up = numeric_pair<mpq>(right_side, y_of_bound);
-            if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return false;
+            if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return;
             m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
             set_upper_bound_witness(j, dep);
             insert_to_columns_with_changed_bounds(j);
@@ -2219,10 +2226,9 @@ namespace lp {
         default:
             UNREACHABLE();
         }
-        return true;
     }
    
-    bool lar_solver::update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
+    void lar_solver::update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
         lp_assert(!column_has_lower_bound(j) && !column_has_upper_bound(j));
 
         mpq y_of_bound(0);
@@ -2259,7 +2265,6 @@ namespace lp {
                 UNREACHABLE();
         }
         insert_to_columns_with_changed_bounds(j);
-        return true;
     }    
 
     lpvar lar_solver::to_column(unsigned ext_j) const {
@@ -2316,7 +2321,7 @@ namespace lp {
         for (unsigned j = 0; j < column_count(); j++) {
             if (!column_is_int(j)) continue;
             if (column_has_term(j)) continue;
-            impq& v = m_mpq_lar_core_solver.m_r_x[j];
+            impq & v = m_mpq_lar_core_solver.r_x(j);
             if (v.is_int())
                 continue;
             TRACE("cube", m_int_solver->display_column(tout, j););
@@ -2353,7 +2358,7 @@ namespace lp {
                 }
             }
             if (need_to_fix) {
-                impq v = t->apply(m_mpq_lar_core_solver.m_r_x);
+                impq v = t->apply(m_mpq_lar_core_solver.r_x());
                 m_mpq_lar_core_solver.m_r_solver.update_x(j, v);
             }
         }
@@ -2365,7 +2370,7 @@ namespace lp {
     bool lar_solver::sum_first_coords(const lar_term& t, mpq& val) const {
         val = zero_of_type<mpq>();
         for (lar_term::ival c : t) {
-            const auto& x = m_mpq_lar_core_solver.m_r_x[c.j()];
+            const auto& x = m_mpq_lar_core_solver.r_x(c.j());
             if (!is_zero(x.y))
                 return false;
             val += x.x * c.coeff();
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index 8ca018072..9c6212c3c 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -105,7 +105,7 @@ class lar_solver : public column_namer {
     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;
     // ((x[j], is_int(j))->j)  for fixed j, used in equalities propagation
     // maps values to integral fixed vars
@@ -139,6 +139,8 @@ class lar_solver : public column_namer {
     bool compare_values(impq const& lhs, lconstraint_kind k, const mpq& rhs);
 
     inline void clear_columns_with_changed_bounds() { m_columns_with_changed_bounds.reset(); }
+
+    struct scoped_backup;
  public:
     const auto& columns_with_changed_bounds() const { return m_columns_with_changed_bounds; }
     void insert_to_columns_with_changed_bounds(unsigned j);
@@ -156,19 +158,19 @@ class lar_solver : public column_namer {
     void add_constraint_to_validate(lar_solver& ls, constraint_index ci);
     bool m_validate_blocker = false;
     void update_column_type_and_bound_check_on_equal(unsigned j, const mpq& right_side, constraint_index ci, unsigned&);
-    bool update_column_type_and_bound(unsigned j, const mpq& right_side, constraint_index ci);
+    void update_column_type_and_bound(unsigned j, const mpq& right_side, constraint_index ci);
  public:   
     bool validate_blocker() const { return m_validate_blocker; }
     bool & validate_blocker() { return m_validate_blocker; }   
-    bool update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
+	void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
  private:
     void require_nbasis_sort() { m_mpq_lar_core_solver.m_r_solver.m_nbasis_sort_counter = 0; }   
-    bool update_column_type_and_bound_with_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
-    bool update_column_type_and_bound_with_no_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
-    bool update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
-    bool update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
-    bool update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
-    bool update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
+    void update_column_type_and_bound_with_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
+    void update_column_type_and_bound_with_no_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
+    void update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
+    void update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
+    void update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
+    void update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
     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(lpvar j, lconstraint_kind kind, const mpq& right_side);
@@ -305,14 +307,12 @@ public:
 
     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; }
 
     std::function<void(lpvar)> m_fixed_var_eh;
     template <typename T>
     void explain_implied_bound(const implied_bound& ib, lp_bound_propagator<T>& bp) {
-        u_dependency* dep = ib.dep();
+        u_dependency* dep = ib.explain_implied();
         for (auto ci : flatten(dep))
             bp.consume(mpq(1), ci); // TODO: flatten should provide the coefficients
         /*
@@ -454,7 +454,7 @@ public:
                                             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.r_x(j);
         auto delta = new_val - x;
         x = new_val;
         after(j);
@@ -540,6 +540,9 @@ public:
     lp_settings const& settings() const;
     statistics& stats();
 
+    void backup_x() { m_mpq_lar_core_solver.backup_x(); }
+    void restore_x() { m_mpq_lar_core_solver.restore_x(); }
+
     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 vector<column_type>&  get_column_types() const { return m_mpq_lar_core_solver.m_column_types(); }
@@ -702,13 +705,13 @@ public:
     void track_touched_rows(bool v);
     bool touched_rows_are_tracked() const;
     ~lar_solver() override;
-    const vector<impq>& r_x() const { return m_mpq_lar_core_solver.m_r_x; }
+    const vector<impq>& r_x() const { return m_mpq_lar_core_solver.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 bool column_value_is_int(unsigned j) const { return m_mpq_lar_core_solver.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; }
     // columns
-    const impq& get_column_value(lpvar j) const { return m_mpq_lar_core_solver.m_r_x[j]; }
+    const impq& get_column_value(lpvar j) const { return m_mpq_lar_core_solver.r_x(j); }
     inline lpvar external_to_local(unsigned j) const {
         lpvar local_j;
         if (m_var_register.external_is_used(j, local_j)) {
diff --git a/src/math/lp/lp_bound_propagator.h b/src/math/lp/lp_bound_propagator.h
index 6f7bee7e5..9da65db04 100644
--- a/src/math/lp/lp_bound_propagator.h
+++ b/src/math/lp/lp_bound_propagator.h
@@ -44,6 +44,26 @@ class lp_bound_propagator {
 public:
     const lar_solver& lp() const { return m_imp.lp(); }
     lar_solver& lp() { return m_imp.lp(); }
+    bool upper_bound_is_available(unsigned j) const {
+        switch (get_column_type(j)) {
+        case column_type::fixed:
+        case column_type::boxed:
+        case column_type::upper_bound:
+            return true;
+        default:
+            return false;
+        }
+    }
+    bool lower_bound_is_available(unsigned j) const {
+        switch (get_column_type(j)) {
+        case column_type::fixed:
+        case column_type::boxed:
+        case column_type::lower_bound:
+            return true;
+        default:
+            return false;
+        }
+    }
 private:
     void try_add_equation_with_internal_fixed_tables(unsigned r1) {
         unsigned v1, v2;
@@ -125,7 +145,8 @@ public:
         return (*m_column_types)[j] == column_type::fixed && get_lower_bound(j).y.is_zero();
     }
 
-    void add_bound(mpq const& v, unsigned j, bool is_low, bool strict, u_dependency* dep) {
+
+    void add_bound(mpq const& v, unsigned j, bool is_low, bool strict, std::function<u_dependency* ()> explain_bound) {
         lconstraint_kind kind = is_low ? GE : LE;
         if (strict)
             kind = static_cast<lconstraint_kind>(kind / 2);
@@ -140,12 +161,12 @@ public:
 
                     found_bound.m_bound = v;
                     found_bound.m_strict = strict;
-                    found_bound.set_dep(dep);
+                    found_bound.set_explain(explain_bound);
                     TRACE("add_bound", lp().print_implied_bound(found_bound, tout););
                 }
             } else {
                 m_improved_lower_bounds.insert(j, static_cast<unsigned>(m_ibounds.size()));
-                m_ibounds.push_back(implied_bound(v, j, is_low, strict, dep));
+                m_ibounds.push_back(implied_bound(v, j, is_low, strict, explain_bound));
                 TRACE("add_bound", lp().print_implied_bound(m_ibounds.back(), tout););
             }
         } else {  // the upper bound case
@@ -156,12 +177,12 @@ public:
 
                     found_bound.m_bound = v;
                     found_bound.m_strict = strict;
-                    found_bound.set_dep(dep);
+                    found_bound.set_explain(explain_bound);
                     TRACE("add_bound", lp().print_implied_bound(found_bound, tout););
                 }
             } else {
                 m_improved_upper_bounds.insert(j, static_cast<unsigned>(m_ibounds.size()));
-                m_ibounds.push_back(implied_bound(v, j, is_low, strict, dep));
+                m_ibounds.push_back(implied_bound(v, j, is_low, strict, explain_bound));
                 TRACE("add_bound", lp().print_implied_bound(m_ibounds.back(), tout););
             }
         }
diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h
index 946b8bf4f..11ea140c0 100644
--- a/src/math/lp/lp_settings.h
+++ b/src/math/lp/lp_settings.h
@@ -138,7 +138,6 @@ struct statistics {
     unsigned m_dio_rewrite_conflicts = 0;
     unsigned m_dio_branching_sats = 0;
     unsigned m_dio_branching_conflicts = 0;
-    unsigned m_dio_bound_propagation_conflicts = 0;
     unsigned m_bounds_tightening_conflicts = 0;
     unsigned m_bounds_tightenings = 0;
     ::statistics m_st = {};
@@ -185,7 +184,6 @@ struct statistics {
         st.update("arith-dio-branching-depth", m_dio_branching_depth);
         st.update("arith-dio-branching-conflicts", m_dio_branching_conflicts);
         st.update("arith-bounds-tightening-conflicts", m_bounds_tightening_conflicts);
-        st.update("arith-dio-propagation-conflicts", m_dio_bound_propagation_conflicts);
         st.update("arith-bounds-tightenings", m_bounds_tightenings);
         st.copy(m_st);
     }