diff --git a/src/math/lp/bound_analyzer_on_row.h b/src/math/lp/bound_analyzer_on_row.h
index 0008a0ee9..576f14599 100644
--- a/src/math/lp/bound_analyzer_on_row.h
+++ b/src/math/lp/bound_analyzer_on_row.h
@@ -26,6 +26,8 @@ Revision History:
 #include "math/lp/test_bound_analyzer.h"
 
 namespace lp {
+  
+ 
 template <typename C, typename B> // C plays a role of a container, B - lp_bound_propagator
 class bound_analyzer_on_row {
     const C&                           m_row;
@@ -91,39 +93,17 @@ private:
     }
 
     bool bound_is_available(unsigned j, bool lower_bound) {
-        return (lower_bound && lower_bound_is_available(j)) ||
-            (!lower_bound && upper_bound_is_available(j));
-    }
-
-    bool upper_bound_is_available(unsigned j) const {
-        switch (m_bp.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 (m_bp.get_column_type(j)) {
-        case column_type::fixed:
-        case column_type::boxed:
-        case column_type::lower_bound:
-            return true;
-        default:
-            return false;
-        }
+        return (lower_bound && m_bp.lower_bound_is_available(j)) ||
+            (!lower_bound && m_bp.upper_bound_is_available(j));
     }
 
     const impq & ub(unsigned j) const {
-        lp_assert(upper_bound_is_available(j));
+        lp_assert(m_bp.upper_bound_is_available(j));
         return m_bp.get_upper_bound(j);
     }
 
     const impq & lb(unsigned j) const {
-        lp_assert(lower_bound_is_available(j));
+        lp_assert(m_bp.lower_bound_is_available(j));
         return m_bp.get_lower_bound(j);
     }
 
@@ -301,10 +281,32 @@ private:
     //     */
     // }
 
-    void limit_j(unsigned j, const mpq& u, bool coeff_before_j_is_pos, bool is_lower_bound, bool strict){
-        m_bp.try_add_bound(u, j, is_lower_bound, coeff_before_j_is_pos, m_row_index, strict);
+    void limit_j(unsigned bound_j, const mpq& u, bool coeff_before_j_is_pos, bool is_lower_bound, bool strict)
+    {
+        unsigned row_index = this->m_row_index;
+        auto* lar = &m_bp.lp();
+        auto explain = [bound_j, coeff_before_j_is_pos, is_lower_bound, strict, row_index, lar]() {
+            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 << ", row_index = " << row_index << "\n";);
+            int bound_sign = (is_lower_bound ? 1 : -1);
+            int j_sign = (coeff_before_j_is_pos ? 1 : -1) * bound_sign;
+
+            SASSERT(!tv::is_term(bound_j));
+            u_dependency* ret = nullptr;
+            for (auto const& r : lar->get_row(row_index)) {
+                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) {
         m_column_of_u = (m_column_of_u == -1) ? j : -2;
     }
@@ -335,6 +337,9 @@ private:
             break;
         }
     }   
+    
 };
+   
+
 }
 
diff --git a/src/math/lp/core_solver_pretty_printer_def.h b/src/math/lp/core_solver_pretty_printer_def.h
index b8048b241..cbe67ea36 100644
--- a/src/math/lp/core_solver_pretty_printer_def.h
+++ b/src/math/lp/core_solver_pretty_printer_def.h
@@ -279,10 +279,12 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print()
         print_row(i);
     }
     m_out << std::endl;
-    if (m_core_solver.inf_heap().size()) {
-        m_out << "inf columns: ";
+    if (!m_core_solver.inf_heap().empty()) {
+        m_out << "inf columns: size() = " << m_core_solver.inf_heap().size() << std::endl;
         print_vector(m_core_solver.inf_heap(), m_out);
         m_out << std::endl;
+    } else {
+        m_out << "inf columns: none\n";
     }
 }
 
diff --git a/src/math/lp/gomory.cpp b/src/math/lp/gomory.cpp
index 775025018..e4267cbd4 100644
--- a/src/math/lp/gomory.cpp
+++ b/src/math/lp/gomory.cpp
@@ -417,7 +417,7 @@ int gomory::find_basic_var() {
 }
     
 lia_move gomory::operator()() {
-    lra.move_non_basic_columns_to_bounds(true);
+    lra.move_non_basic_columns_to_bounds();
     int j = find_basic_var();
     if (j == -1)
         return lia_move::undef;
diff --git a/src/math/lp/implied_bound.h b/src/math/lp/implied_bound.h
index 9435edcdc..195ec0359 100644
--- a/src/math/lp/implied_bound.h
+++ b/src/math/lp/implied_bound.h
@@ -21,37 +21,40 @@ Revision History:
 #include "math/lp/lp_settings.h"
 #include "math/lp/lar_constraints.h"
 namespace lp {
-struct implied_bound {
+class implied_bound {
+    public:
     mpq m_bound;
-    unsigned m_j; // the column for which the bound has been found
+    // It is either the column for which the bound has been found, or,
+    // in the case the column was created as
+    // the slack variable to a term, it is the term index.
+    // It is the same index that was returned by lar_solver::add_var(), or
+    // by lar_solver::add_term()
+    unsigned m_j; 
     bool m_is_lower_bound;
-    bool m_coeff_before_j_is_pos;
-    unsigned m_row_or_term_index;
     bool m_strict;
+    private:
+    std::function<u_dependency*()> m_explain_bound = nullptr;
+    public:
+    // s is expected to be the pointer to lp_bound_propagator.
+    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)
             k = static_cast<lconstraint_kind>(k / 2);
         return k;
     }
-    bool operator==(const implied_bound & o) const {
-        return m_j == o.m_j && m_is_lower_bound == o.m_is_lower_bound && m_bound == o.m_bound &&
-            m_coeff_before_j_is_pos == o.m_coeff_before_j_is_pos &&
-            m_row_or_term_index == o.m_row_or_term_index && m_strict == o.m_strict;
-    }
     implied_bound(){}
     implied_bound(const mpq & a,
                   unsigned j,
-                  bool lower_bound,
-                  bool coeff_before_j_is_pos,
-                  unsigned row_or_term_index,
-                  bool strict):
+                  bool is_lower_bound,
+                  bool is_strict,
+                  std::function<u_dependency*()> get_dep):
         m_bound(a),
         m_j(j),
-        m_is_lower_bound(lower_bound),
-        m_coeff_before_j_is_pos(coeff_before_j_is_pos),
-        m_row_or_term_index(row_or_term_index),
-        m_strict(strict) {
+        m_is_lower_bound(is_lower_bound),
+        m_strict(is_strict),
+        m_explain_bound(get_dep) {
     }
 };
 }
diff --git a/src/math/lp/int_branch.cpp b/src/math/lp/int_branch.cpp
index 4bfe2b827..d6262a4d0 100644
--- a/src/math/lp/int_branch.cpp
+++ b/src/math/lp/int_branch.cpp
@@ -24,7 +24,7 @@ namespace lp {
 int_branch::int_branch(int_solver& lia):lia(lia), lra(lia.lra) {}
 
 lia_move int_branch::operator()() {
-    lra.move_non_basic_columns_to_bounds(true);
+    lra.move_non_basic_columns_to_bounds();
     int j = find_inf_int_base_column();
     return j == -1? lia_move::sat : create_branch_on_column(j);        
 }
diff --git a/src/math/lp/int_cube.cpp b/src/math/lp/int_cube.cpp
index da724a543..9ef9aa341 100644
--- a/src/math/lp/int_cube.cpp
+++ b/src/math/lp/int_cube.cpp
@@ -43,7 +43,7 @@ namespace lp {
         if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
             TRACE("cube", tout << "cannot find a feasible solution";);
             lra.pop();
-            lra.move_non_basic_columns_to_bounds(false);
+            lra.move_non_basic_columns_to_bounds();
             // it can happen that we found an integer solution here
             return !lra.r_basis_has_inf_int()? lia_move::sat: lia_move::undef;
         }
diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index 1bd55b9bb..ebe14e407 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -22,7 +22,8 @@ namespace lp {
     }
 
     lar_solver::lar_solver() :
-        m_crossed_bounds_column(-1),
+        m_crossed_bounds_column(null_lpvar),
+        m_crossed_bounds_deps(nullptr),
         m_mpq_lar_core_solver(m_settings, *this),
         m_var_register(false),
         m_term_register(true),
@@ -180,7 +181,10 @@ namespace lp {
 
     lp_status lar_solver::get_status() const { return m_status; }
 
-    void lar_solver::set_status(lp_status s) { m_status = s; }
+    void lar_solver::set_status(lp_status s) { 
+        TRACE("lar_solver", tout << "setting status to " << s << "\n";);
+        m_status = s; 
+    }
 
     lp_status lar_solver::find_feasible_solution() {
         stats().m_make_feasible++;
@@ -199,9 +203,9 @@ namespace lp {
     }
 
     lp_status lar_solver::solve() {
-        if (m_status == lp_status::INFEASIBLE) {
+        if (m_status == lp_status::INFEASIBLE) 
             return m_status;
-        }
+        
         solve_with_core_solver();
         if (m_status != lp_status::INFEASIBLE) {
             if (m_settings.bound_propagation())
@@ -213,17 +217,10 @@ namespace lp {
     }
 
     void lar_solver::fill_explanation_from_crossed_bounds_column(explanation& evidence) const {
-        lp_assert(static_cast<int>(get_column_type(m_crossed_bounds_column)) >= static_cast<int>(column_type::boxed));
-        lp_assert(!column_is_feasible(m_crossed_bounds_column));
-
         // this is the case when the lower bound is in conflict with the upper one
-        const ul_pair& ul = m_columns_to_ul_pairs[m_crossed_bounds_column];
         svector<constraint_index> deps;
-        m_dependencies.linearize(ul.upper_bound_witness(), deps);
-        for (auto d : deps)
-            evidence.add_pair(d, numeric_traits<mpq>::one());
-        deps.reset();
-        m_dependencies.linearize(ul.lower_bound_witness(), deps);
+        SASSERT(m_crossed_bounds_deps != nullptr);
+        m_dependencies.linearize(m_crossed_bounds_deps, deps);
         for (auto d : deps)
             evidence.add_pair(d, -numeric_traits<mpq>::one());
     }
@@ -232,7 +229,8 @@ namespace lp {
         m_simplex_strategy = m_settings.simplex_strategy();
         m_simplex_strategy.push();
         m_columns_to_ul_pairs.push();
-        m_crossed_bounds_column.push();
+        m_crossed_bounds_column = null_lpvar;
+        m_crossed_bounds_deps = nullptr;
         m_mpq_lar_core_solver.push();
         m_term_count = m_terms.size();
         m_term_count.push();
@@ -262,7 +260,8 @@ namespace lp {
 
     void lar_solver::pop(unsigned k) {
         TRACE("lar_solver", tout << "k = " << k << std::endl;);
-        m_crossed_bounds_column.pop(k);
+        m_crossed_bounds_column = null_lpvar;
+        m_crossed_bounds_deps = nullptr;
         unsigned n = m_columns_to_ul_pairs.peek_size(k);
         m_var_register.shrink(n);
         pop_tableau(n);
@@ -403,7 +402,7 @@ namespace lp {
 
     void lar_solver::prepare_costs_for_r_solver(const lar_term& term) {
         TRACE("lar_solver", print_term(term, tout << "prepare: ") << "\n";);
-        move_non_basic_columns_to_bounds(false);
+        move_non_basic_columns_to_bounds();
         auto& rslv = m_mpq_lar_core_solver.m_r_solver;
         lp_assert(costs_are_zeros_for_r_solver());
         lp_assert(reduced_costs_are_zeroes_for_r_solver());
@@ -420,11 +419,11 @@ namespace lp {
         lp_assert(rslv.reduced_costs_are_correct_tableau());
     }
 
-    void lar_solver::move_non_basic_columns_to_bounds(bool shift_randomly) {
+    void lar_solver::move_non_basic_columns_to_bounds() {
         auto& lcs = m_mpq_lar_core_solver;
         bool change = false;
         for (unsigned j : lcs.m_r_nbasis) {
-            if (move_non_basic_column_to_bounds(j, shift_randomly))
+            if (move_non_basic_column_to_bounds(j))
                 change = true;
         }
         if (!change)
@@ -435,46 +434,40 @@ namespace lp {
         find_feasible_solution();
     }
 
-    bool lar_solver::move_non_basic_column_to_bounds(unsigned j, bool force_change) {
+    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];
         switch (lcs.m_column_types()[j]) {
         case column_type::boxed: {
-            bool at_l = val == lcs.m_r_lower_bounds()[j];
-            bool at_u = !at_l && (val == lcs.m_r_upper_bounds()[j]);
-            if (!at_l && !at_u) {
-                if (m_settings.random_next() % 2)
-                    set_value_for_nbasic_column(j, lcs.m_r_lower_bounds()[j]);
-                else
-                    set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
+            const auto& l = lcs.m_r_lower_bounds()[j];
+            if (val == l || val == lcs.m_r_upper_bounds()[j]) return false;
+            set_value_for_nbasic_column(j, l);
+            return true;
+        }
+                                        
+        case column_type::lower_bound: {
+            const auto& l = lcs.m_r_lower_bounds()[j];
+            if (val != l) {
+                set_value_for_nbasic_column(j, l);
                 return true;
             }
-            else if (force_change && m_settings.random_next() % 3 == 0) {
-                set_value_for_nbasic_column(j,
-                    at_l ? lcs.m_r_upper_bounds()[j] : lcs.m_r_lower_bounds()[j]);
-                return true;
-            }
-            break;
-        }                               
-        case column_type::lower_bound:
-            if (val != lcs.m_r_lower_bounds()[j]) {
-                set_value_for_nbasic_column(j, lcs.m_r_lower_bounds()[j]);
-                return true;
-            }
-            break;
+            return false;
+        }
         case column_type::fixed:
-        case column_type::upper_bound:
-            if (val != lcs.m_r_upper_bounds()[j]) {
-                set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
+        case column_type::upper_bound: {
+            const auto & u = lcs.m_r_upper_bounds()[j];
+            if (val != u) {
+                set_value_for_nbasic_column(j, u);
                 return true;
             }
-            break;
+            return false;
+        }
         case column_type::free_column:
             if (column_is_int(j) && !val.is_int()) {
                 set_value_for_nbasic_column(j, impq(floor(val)));
                 return true;
             }
-            break;
+            return false;
         default:
             SASSERT(false);
         }
@@ -486,6 +479,8 @@ namespace lp {
         auto& x = m_mpq_lar_core_solver.m_r_x[j];
         auto delta = new_val - x;
         x = new_val;
+        TRACE("lar_solver_feas", tout << "setting " << j << " to "
+         << new_val << (column_is_feasible(j)?"feas":"non-feas") << "\n";);
         change_basic_columns_dependend_on_a_given_nb_column(j, delta);
     }
 
@@ -793,6 +788,8 @@ namespace lp {
     void lar_solver::detect_rows_with_changed_bounds() {
         for (auto j : m_columns_with_changed_bounds)
             detect_rows_with_changed_bounds_for_column(j);
+        if (m_find_monics_with_changed_bounds_func)
+            m_find_monics_with_changed_bounds_func(m_columns_with_changed_bounds);
     }
 
     void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau() {
@@ -1021,7 +1018,7 @@ namespace lp {
 
     void lar_solver::get_infeasibility_explanation(explanation& exp) const {
         exp.clear();
-        if (m_crossed_bounds_column != -1) {
+        if (m_crossed_bounds_column != null_lpvar) {
             fill_explanation_from_crossed_bounds_column(exp);
             return;
         }
@@ -1073,12 +1070,16 @@ namespace lp {
     }
 
     bool lar_solver::init_model() const {
-        if (get_status() != lp_status::OPTIMAL && get_status() != lp_status::FEASIBLE)
+        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());
+        if (status != lp_status::OPTIMAL && status != lp_status::FEASIBLE)
             return false;
         if (!m_columns_with_changed_bounds.empty())
             return false;
 
-        lp_assert(m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis());
         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();
@@ -1625,10 +1626,9 @@ namespace lp {
         SASSERT(m_terms.size() == m_term_register.size());
         unsigned adjusted_term_index = m_terms.size() - 1;
         var_index ret = tv::mask_term(adjusted_term_index);
-        if (!coeffs.empty()) {
+        if (!coeffs.empty())
             add_row_from_term_no_constraint(m_terms.back(), ret);
-            add_touched_row(A_r().row_count() - 1);
-        }
+        
         lp_assert(m_var_register.size() == A_r().column_count());
         if (m_need_register_terms) 
             register_normalized_term(*t, A_r().column_count() - 1);
@@ -1828,7 +1828,6 @@ namespace lp {
 
         if (is_base(j) && column_is_fixed(j))
             m_fixed_base_var_set.insert(j);
-
         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;);    
     }
 
@@ -1924,7 +1923,6 @@ namespace lp {
         }
     }
 
-    // clang-format on
     void lar_solver::update_bound_with_ub_lb(var_index 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 ||
@@ -1937,12 +1935,14 @@ namespace lp {
             case LE: {
                 auto up = numeric_pair<mpq>(right_side, y_of_bound);
                 if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
-                    set_infeasible_column(j);
+                    set_crossed_bounds_column_and_deps(j, true, dep);
+                }
+                else {
+                    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);
                 }
-                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);
                 break;
             }
             case GT:
@@ -1950,25 +1950,33 @@ namespace lp {
             case GE: {
                 auto low = numeric_pair<mpq>(right_side, y_of_bound);
                 if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                    set_infeasible_column(j);
+                    set_crossed_bounds_column_and_deps(j, false, dep);
+                } else {
+                    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);
+                    insert_to_columns_with_changed_bounds(j);
                 }
-                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);
-                insert_to_columns_with_changed_bounds(j);
                 break;
+                
             } 
             case EQ: {
                 auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
-                if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j] || v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
-                    set_infeasible_column(j);
+                if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]){
+                    set_crossed_bounds_column_and_deps(j, false, dep);
+                }  
+                else if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
+                    set_crossed_bounds_column_and_deps(j, true, dep);
+                }
+                else {
+                    set_upper_bound_witness(j, dep);
+                    set_lower_bound_witness(j, dep);
+                    m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
+                    insert_to_columns_with_changed_bounds(j);
                 }
-                set_upper_bound_witness(j, dep);
-                set_lower_bound_witness(j, dep);
-                m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
                 break;
             }
 
@@ -1979,7 +1987,7 @@ namespace lp {
             m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
         }
     }
-    // clang-format off
+    
     void lar_solver::update_bound_with_no_ub_lb(var_index 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);
@@ -1991,12 +1999,14 @@ namespace lp {
             case LE: {
                 auto up = numeric_pair<mpq>(right_side, y_of_bound);
                 if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
-                    set_infeasible_column(j);
+                    set_crossed_bounds_column_and_deps(j, true, dep);
+                }
+                else {
+                    m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+                    set_upper_bound_witness(j, dep);
+                    m_mpq_lar_core_solver.m_column_types[j] = (up == m_mpq_lar_core_solver.m_r_lower_bounds[j] ? column_type::fixed : column_type::boxed);
+                    insert_to_columns_with_changed_bounds(j);
                 }
-                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
-                set_upper_bound_witness(j, dep);
-                m_mpq_lar_core_solver.m_column_types[j] = (up == m_mpq_lar_core_solver.m_r_lower_bounds[j] ? column_type::fixed : column_type::boxed);
-                insert_to_columns_with_changed_bounds(j);
                 break;
             } 
             case GT:
@@ -2014,13 +2024,15 @@ namespace lp {
             case EQ: {
                 auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
                 if (v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
-                    set_infeasible_column(j);
+                    set_crossed_bounds_column_and_deps(j, true, dep);
+                } 
+                else {
+                    set_upper_bound_witness(j, dep);
+                    set_lower_bound_witness(j, dep);
+                    m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
+                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+                    insert_to_columns_with_changed_bounds(j);
                 }
-
-                set_upper_bound_witness(j, dep);
-                set_lower_bound_witness(j, dep);
-                m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
-                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
                 break;
             }
 
@@ -2051,26 +2063,29 @@ namespace lp {
         {
             auto low = numeric_pair<mpq>(right_side, y_of_bound);
             if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                set_infeasible_column(j);
+                set_crossed_bounds_column_and_deps(j, false, dep);
+            }
+            else {
+                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);
+                insert_to_columns_with_changed_bounds(j);
             }
-            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);
-            insert_to_columns_with_changed_bounds(j);
-            
         }
         break;
         case EQ:
         {
             auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
             if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
-                set_infeasible_column(j);
+                set_crossed_bounds_column_and_deps(j, false, dep);
+            }
+            else {
+                set_upper_bound_witness(j, dep);
+                set_lower_bound_witness(j, dep);
+                m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+                insert_to_columns_with_changed_bounds(j);
             }
-
-            set_upper_bound_witness(j, dep);
-            set_lower_bound_witness(j, dep);
-            m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v;
-            m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
             break;
         }
 
@@ -2345,6 +2360,27 @@ namespace lp {
             return false;
         return true;
     }
+    // If lower_bound is true than the new asserted upper bound is less than the existing lower bound.
+    // Otherwise the new asserted lower bound is is greater than the existing upper bound.
+    // dep is the reason for the new bound
+
+    void lar_solver::set_crossed_bounds_column_and_deps(unsigned j, bool lower_bound, u_dependency* dep) {
+        if (m_crossed_bounds_column != null_lpvar) return; // already set
+        SASSERT(m_crossed_bounds_deps == nullptr);
+        set_status(lp_status::INFEASIBLE);
+        m_crossed_bounds_column = j;
+        const auto& ul = this->m_columns_to_ul_pairs()[j];
+        u_dependency* bdep = lower_bound? ul.lower_bound_witness() : ul.upper_bound_witness();
+        SASSERT(bdep != nullptr);
+        m_crossed_bounds_deps = m_dependencies.mk_join(bdep, dep);
+        insert_to_columns_with_changed_bounds(j);
+    }
+
+    void lar_solver::collect_more_rows_for_lp_propagation(){
+        for (auto j : m_columns_with_changed_bounds)
+            detect_rows_with_changed_bounds_for_column(j); 
+    }
+
 
 } // namespace lp
 
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index 275f7446d..d902f9a3a 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -78,7 +78,8 @@ class lar_solver : public column_namer {
     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;
+    lpvar m_crossed_bounds_column;
+    u_dependency* m_crossed_bounds_deps;
     lar_core_solver m_mpq_lar_core_solver;
     int_solver* m_int_solver = nullptr;
     bool m_need_register_terms = false;
@@ -139,12 +140,22 @@ 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(); }
+ public:
+    const auto& columns_with_changed_bounds() const { return m_columns_with_changed_bounds; }
     void insert_to_columns_with_changed_bounds(unsigned j);
+    const u_dependency* crossed_bounds_deps() const { return m_crossed_bounds_deps;}
+    u_dependency*& crossed_bounds_deps() { return m_crossed_bounds_deps;}
+
+    lpvar crossed_bounds_column() const { return m_crossed_bounds_column; }
+    lpvar& crossed_bounds_column() { return m_crossed_bounds_column; } 
+        
+
+ private:   
     void update_column_type_and_bound_check_on_equal(unsigned j, const mpq& right_side, constraint_index ci, unsigned&);
     void update_column_type_and_bound(unsigned j, const mpq& right_side, constraint_index ci);
-public:
+ public:   
     void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
-private:
+ private:   
     void update_column_type_and_bound_with_ub(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
     void update_column_type_and_bound_with_no_ub(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
     void update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
@@ -154,10 +165,7 @@ private:
     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);
-    inline void set_infeasible_column(unsigned j) {
-        set_status(lp_status::INFEASIBLE);
-        m_crossed_bounds_column = j;
-    }
+    void set_crossed_bounds_column_and_deps(unsigned j, bool lower_bound, u_dependency* dep);
     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);
     unsigned row_of_basic_column(unsigned) const;
@@ -304,26 +312,34 @@ private:
 
     template <typename T>
     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))
-            bound_j = m_var_register.external_to_local(bound_j);
+        u_dependency* dep = ib.explain_implied();
+        for (auto ci : flatten(dep))
+            bp.consume(mpq(1), ci); // TODO: flatten should provide the coefficients
+        /*
+        if (ib.m_is_monic) {
+            NOT_IMPLEMENTED_YET();
+        } else {
+            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))
+                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)
-                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];
-            auto* witness = sign > 0 ? ul.upper_bound_witness() : ul.lower_bound_witness();
-            lp_assert(witness);
-            for (auto ci : flatten(witness))
-                bp.consume(a, ci);
-        }
+            for (auto const& r : get_row(i)) {
+                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;
+                const ul_pair& ul = m_columns_to_ul_pairs[j];
+                auto* witness = sign > 0 ? ul.upper_bound_witness() : ul.lower_bound_witness();
+                lp_assert(witness);
+                for (auto ci : flatten(witness))
+                    bp.consume(a, ci);
+            }
+            }*/
     }
 
     void set_value_for_nbasic_column(unsigned j, const impq& new_val);
@@ -363,7 +379,7 @@ private:
         }
         m_touched_rows.reset();
     }
-
+    void collect_more_rows_for_lp_propagation();
     template <typename T>
     void check_missed_propagations(lp_bound_propagator<T>& bp) {
         for (unsigned i = 0; i < A_r().row_count(); i++)
@@ -557,6 +573,16 @@ private:
         const ul_pair& ul = m_columns_to_ul_pairs[j];
         return m_dependencies.mk_join(ul.lower_bound_witness(), ul.upper_bound_witness());
     }
+    template <typename T>
+    u_dependency* get_bound_constraint_witnesses_for_columns(const T& collection) {
+        u_dependency* dep = nullptr;
+        for (auto j : collection) {
+            u_dependency* d = get_bound_constraint_witnesses_for_column(j);
+            dep = m_dependencies.mk_join(dep, d);
+        }
+        return dep;
+    }
+    u_dependency* join_deps(u_dependency* a, u_dependency *b) { return m_dependencies.mk_join(a, b); }
     inline constraint_set const& constraints() const { return m_constraints; }
     void push();
     void pop();
@@ -609,8 +635,8 @@ private:
         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);
+    void move_non_basic_columns_to_bounds();
+    bool move_non_basic_column_to_bounds(unsigned j);
     inline bool r_basis_has_inf_int() const {
         for (unsigned j : r_basis()) {
             if (column_is_int(j) && !column_value_is_int(j))
@@ -663,6 +689,7 @@ private:
             return 0;
         return m_usage_in_terms[j];
     }
+    std::function<void (const indexed_uint_set& columns_with_changed_bound)> m_find_monics_with_changed_bounds_func = nullptr;
     friend int_solver;
     friend int_branch;
 };
diff --git a/src/math/lp/lp_bound_propagator.h b/src/math/lp/lp_bound_propagator.h
index 9bd485b1a..c6fa5796d 100644
--- a/src/math/lp/lp_bound_propagator.h
+++ b/src/math/lp/lp_bound_propagator.h
@@ -6,19 +6,21 @@
 */
 #pragma once
 #include <utility>
-
 #include "math/lp/lp_settings.h"
 #include "util/uint_set.h"
+#include "math/lp/implied_bound.h"
+#include "util/vector.h"
 namespace lp {
+    
 template <typename T>
 class lp_bound_propagator {
-	uint_set m_visited_rows;
+    uint_set m_visited_rows;
     // 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;
+    u_map<unsigned> m_improved_lower_bounds;
+    u_map<unsigned> m_improved_upper_bounds;
 
     T& m_imp;
-    vector<implied_bound> m_ibounds;
+    std_vector<implied_bound>& m_ibounds;
 
     map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>> m_val2fixed_row;
     // works for rows of the form x + y + sum of fixed = 0
@@ -39,7 +41,30 @@ class lp_bound_propagator {
         }
         return x != UINT_MAX;
     }
-
+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;
         if (!only_one_nfixed(r1, v1))
@@ -83,21 +108,18 @@ class lp_bound_propagator {
         ~reset_cheap_eq() { p.reset_cheap_eq_eh(); }
     };
 
-   public:
-    lp_bound_propagator(T& imp) : m_imp(imp) {}
+public:
+    lp_bound_propagator(T& imp, std_vector<implied_bound> & ibounds) : m_imp(imp), m_ibounds(ibounds) {}
 
-    const vector<implied_bound>& ibounds() const { return m_ibounds; }
+    const std_vector<implied_bound>& ibounds() const { return m_ibounds; }
 
     void init() {
-        m_improved_upper_bounds.clear();
-        m_improved_lower_bounds.clear();
-        m_ibounds.reset();
+        m_improved_upper_bounds.reset();
+        m_improved_lower_bounds.reset();
+        m_ibounds.clear();
         m_column_types = &lp().get_column_types();
     }
-
-    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_column_types)[j];
     }
@@ -123,7 +145,8 @@ class lp_bound_propagator {
         return (*m_column_types)[j] == column_type::fixed && 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) {
+
+    void add_bound(mpq const& v, unsigned j, bool is_low, bool strict, std::function<u_dependency* ()> explain_bound) {
         j = lp().column_to_reported_index(j);
 
         lconstraint_kind kind = is_low ? GE : LE;
@@ -132,30 +155,37 @@ class lp_bound_propagator {
 
         if (!m_imp.bound_is_interesting(j, kind, v))
             return;
-        unsigned k;  // index to ibounds
         if (is_low) {
-            if (try_get_value(m_improved_lower_bounds, j, k)) {
+            unsigned k;
+            if (m_improved_lower_bounds.find(j, 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", lp().print_implied_bound(found_bound, tout););
+
+                    found_bound.m_bound = v;
+                    found_bound.m_strict = strict;
+                    found_bound.set_explain(explain_bound);
+                    TRACE("add_bound", lp().print_implied_bound(found_bound, tout););
                 }
             } 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", lp().print_implied_bound(m_ibounds.back(), tout););
+                m_improved_lower_bounds.insert(j, static_cast<unsigned>(m_ibounds.size()));
+                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
-            if (try_get_value(m_improved_upper_bounds, j, k)) {
+            unsigned k;
+            if (m_improved_upper_bounds.find(j, 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", lp().print_implied_bound(found_bound, tout););
+
+                    found_bound.m_bound = v;
+                    found_bound.m_strict = strict;
+                    found_bound.set_explain(explain_bound);
+                    TRACE("add_bound", lp().print_implied_bound(found_bound, tout););
                 }
             } 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", lp().print_implied_bound(m_ibounds.back(), tout););
+                m_improved_upper_bounds.insert(j, static_cast<unsigned>(m_ibounds.size()));
+                m_ibounds.push_back(implied_bound(v, j, is_low, strict, explain_bound));
+                TRACE("add_bound", lp().print_implied_bound(m_ibounds.back(), tout););
             }
         }
     }
@@ -383,7 +413,8 @@ class lp_bound_propagator {
             lp_assert(y_sign == 1 || y_sign == -1);
             auto& table = y_sign == 1 ? m_row2index_pos : m_row2index_neg;
             const auto& v = val(x);
-            unsigned found_i;
+            unsigned found_i;;
+            
             if (!table.find(v, found_i)) {
                 table.insert(v, i);
             } else {
diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h
index e7570d5dd..2c0993d71 100644
--- a/src/math/lp/lp_core_solver_base.h
+++ b/src/math/lp/lp_core_solver_base.h
@@ -540,9 +540,9 @@ public:
     }
 
     void add_delta_to_x_and_track_feasibility(unsigned j, const X & del) {
-        TRACE("lar_solver_feas_bug", tout << "del = " << del << ", was x[" << j << "] = " << m_x[j] << "\n";);
+        TRACE("lar_solver_feas", tout << "del = " << del << ", was x[" << j << "] = " << m_x[j] << "\n";);
         m_x[j] += del;
-        TRACE("lar_solver_feas_bug", tout << "became x[" << j << "] = " << m_x[j] << "\n";);
+        TRACE("lar_solver_feas", tout << "became x[" << j << "] = " << m_x[j] << "\n";);
         track_column_feasibility(j);
     }
 
@@ -564,6 +564,7 @@ public:
     }
     void insert_column_into_inf_heap(unsigned j) {        
 		if (!m_inf_heap.contains(j)) {
+            m_inf_heap.reserve(j+1);
 	        m_inf_heap.insert(j);
             TRACE("lar_solver_inf_heap", tout << "insert into inf_heap j = " << j << "\n";);
         }
@@ -571,7 +572,7 @@ public:
     }
     void remove_column_from_inf_heap(unsigned j) {
 		if (m_inf_heap.contains(j)) {
-            TRACE("lar_solver_inf_heap", tout << "insert into heap j = " << j << "\n";);
+            TRACE("lar_solver_inf_heap", tout << "erase from heap j = " << j << "\n";);
         	m_inf_heap.erase(j);
         }
         lp_assert(column_is_feasible(j));
diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h
index 0552b8e82..b4a53e872 100644
--- a/src/math/lp/lp_core_solver_base_def.h
+++ b/src/math/lp/lp_core_solver_base_def.h
@@ -115,10 +115,12 @@ pretty_print(std::ostream & out) {
 
 template <typename T, typename X> void lp_core_solver_base<T, X>::
 add_delta_to_entering(unsigned entering, const X& delta) {
-    m_x[entering] += delta;     
+    m_x[entering] += delta;  
+    TRACE("lar_solver_feas", tout << "not tracking feas entering = " << entering << " = " << m_x[entering] << (column_is_feasible(entering) ? " feas" : " non-feas") << "\n";); 
     for (const auto & c : m_A.m_columns[entering]) {
         unsigned i = c.var();
         m_x[m_basis[i]] -= delta * m_A.get_val(c);
+        TRACE("lar_solver_feas", tout << "not tracking feas m_basis[i] = " << m_basis[i] << " = " << m_x[m_basis[i]] << (column_is_feasible(m_basis[i]) ? " feas" : " non-feas") << "\n";);
     }
 }
 
diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h
index 9871ec691..f5d7314db 100644
--- a/src/math/lp/lp_primal_core_solver.h
+++ b/src/math/lp/lp_primal_core_solver.h
@@ -394,9 +394,10 @@ namespace lp {
         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 will remove the leaving from the heap
-        TRACE("lar_solver_inf_heap", tout << "leaving = " << leaving
+        TRACE("lar_solver_feas", tout << "entering = " << entering << ", leaving = " << leaving << ", new_val_for_leaving = " << new_val_for_leaving << ", theta = " << theta << "\n";);
+        TRACE("lar_solver_feas", tout << "leaving = " << leaving
                                  << " removed from inf_heap()\n";);
+        // this will remove the leaving from the heap
         this->inf_heap().erase_min();
         advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta);
         if (this->current_x_is_feasible())
diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp
index 9bc49fc43..81c00d2f0 100644
--- a/src/math/lp/monomial_bounds.cpp
+++ b/src/math/lp/monomial_bounds.cpp
@@ -260,20 +260,31 @@ namespace nla {
     }
 
     void monomial_bounds::unit_propagate() {        
-        for (auto const& m : c().m_emons) 
-            unit_propagate(m);
+        for (lpvar v : c().m_monics_with_changed_bounds) {
+            if (!c().is_monic_var(v)) continue;
+            unit_propagate(c().emons()[v]);
+            if (c().lra.get_status() == lp::lp_status::INFEASIBLE) {
+                lp::explanation exp;
+                c().lra.get_infeasibility_explanation(exp);
+                new_lemma lemma(c(), "propagate fixed - infeasible lra");
+                lemma &= exp;
+                break;
+            }
+            if (c().m_conflicts > 0 ) {
+                break;
+            }
+        }   
     }
 
-
-    void monomial_bounds::unit_propagate(monic const& m) {
+    void monomial_bounds::unit_propagate(monic & m) {
         if (m.is_propagated())
             return;
 
         if (!is_linear(m))
             return;
 
-        c().m_emons.set_propagated(m);
-        
+        c().emons().set_propagated(m);
+
         rational k = fixed_var_product(m);
         lpvar w = non_fixed_var(m);
         if (w == null_lpvar || k == 0)
@@ -297,10 +308,12 @@ namespace nla {
             lp::impq val(k);
             c().lra.set_value_for_nbasic_column(m.var(), val);
         }
+        TRACE("nla_solver", tout << "propagate fixed " << m << " = " << k << "\n";);
         c().lra.update_column_type_and_bound(m.var(), lp::lconstraint_kind::EQ, k, dep);
+        
         // propagate fixed equality
-        auto exp = get_explanation(dep);
-        c().add_fixed_equality(m.var(), k, exp);
+        auto exp = get_explanation(dep);        
+        c().add_fixed_equality(c().lra.column_to_reported_index(m.var()), k, exp);
     }
 
     void monomial_bounds::propagate_nonfixed(monic const& m, rational const& k, lpvar w) {
@@ -311,11 +324,12 @@ namespace nla {
         lp::lpvar term_index = c().lra.add_term(coeffs, UINT_MAX);
         auto* dep = explain_fixed(m, k);
         term_index = c().lra.map_term_index_to_column_index(term_index);
+        TRACE("nla_solver", tout << "propagate nonfixed " << m << " = " << k << " " << w << "\n";);
         c().lra.update_column_type_and_bound(term_index, lp::lconstraint_kind::EQ, mpq(0), dep);
 
         if (k == 1) {
             lp::explanation exp = get_explanation(dep);
-            c().add_equality(m.var(), w, exp);
+            c().add_equality(c().lra.column_to_reported_index(m.var()), c().lra.column_to_reported_index(w), exp);
         }
     }
 
@@ -356,8 +370,9 @@ namespace nla {
     rational monomial_bounds::fixed_var_product(monic const& m) {
         rational r(1);
         for (lpvar v : m) {
+            //  we have to use the column bounds here, because the column value may be outside the bounds
             if (c().var_is_fixed(v))
-                r *= c().lra.get_column_value(v).x;
+                r *= c().lra.get_lower_bound(v).x;            
         }
         return r;
     }
diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h
index 747aca9a2..b2079b4f1 100644
--- a/src/math/lp/monomial_bounds.h
+++ b/src/math/lp/monomial_bounds.h
@@ -35,11 +35,10 @@ namespace nla {
         bool is_zero(lpvar v) const;
 
         // monomial propagation
-        void unit_propagate(monic const& m);
+        void unit_propagate(monic & m);
         bool is_linear(monic const& m);
         rational fixed_var_product(monic const& m);
         lpvar non_fixed_var(monic const& m);
-
     public:
         monomial_bounds(core* core);
         void propagate();
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index e65829ed7..2a23ba47b 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -41,6 +41,17 @@ core::core(lp::lar_solver& s, params_ref const& p, reslimit & lim) :
     m_nra(s, m_nra_lim, *this) 
 {
     m_nlsat_delay = lp_settings().nlsat_delay();
+    lra.m_find_monics_with_changed_bounds_func = [&](const indexed_uint_set& columns_with_changed_bounds) {
+        for (lpvar j : columns_with_changed_bounds) {
+            if (is_monic_var(j))
+                m_monics_with_changed_bounds.insert(j);
+            else {
+                for (const auto & m: m_emons.get_use_list(j)) {
+                    m_monics_with_changed_bounds.insert(m.var());
+                }
+            }    
+        }
+    };
 }
     
 bool core::compare_holds(const rational& ls, llc cmp, const rational& rs) const {
@@ -137,6 +148,7 @@ void core::add_monic(lpvar v, unsigned sz, lpvar const* vs) {
         m_add_buffer[i] = j;
     }
     m_emons.add(v, m_add_buffer);
+    m_monics_with_changed_bounds.insert(v);
 }
     
 void core::push() {
@@ -812,6 +824,7 @@ void core::clear() {
     m_literals.clear();
     m_fixed_equalities.clear();
     m_equalities.clear();
+    m_conflicts = 0;
 }
     
 void core::init_search() {
@@ -1065,6 +1078,9 @@ new_lemma::~new_lemma() {
     (void)i;
     (void)name;
     // code for checking lemma can be added here
+    if (current().is_conflict()) {
+        c.m_conflicts++;
+    }
     TRACE("nla_solver", tout << name << " " << (++i) << "\n" << *this; );
 }
 
@@ -1813,6 +1829,7 @@ bool core::improve_bounds() {
 void core::propagate() {
     clear();
     m_monomial_bounds.unit_propagate();
+    m_monics_with_changed_bounds.reset();
 }
 
 
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index 17aa304ac..ae76f6d5a 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -89,6 +89,7 @@ class core {
     vector<equality>         m_equalities;
     vector<fixed_equality>   m_fixed_equalities;
     indexed_uint_set         m_to_refine;
+    indexed_uint_set         m_monics_with_changed_bounds;
     tangents                 m_tangents;
     basics                   m_basics;
     order                    m_order;
@@ -97,7 +98,7 @@ class core {
     divisions                m_divisions;
     intervals                m_intervals; 
     monomial_bounds          m_monomial_bounds;
-    
+    unsigned                 m_conflicts;
     horner                   m_horner;
     grobner                  m_grobner;
     emonics                  m_emons;
@@ -120,7 +121,7 @@ class core {
 public:    
     // constructor
     core(lp::lar_solver& s, params_ref const& p, reslimit&);
-
+    const auto& monics_with_changed_bounds() const { return m_monics_with_changed_bounds; }
     void insert_to_refine(lpvar j);
     void erase_from_to_refine(lpvar j);
     
diff --git a/src/math/lp/nla_solver.h b/src/math/lp/nla_solver.h
index fec27c32b..c508e68d0 100644
--- a/src/math/lp/nla_solver.h
+++ b/src/math/lp/nla_solver.h
@@ -26,7 +26,7 @@ namespace nla {
 
         solver(lp::lar_solver& s, params_ref const& p, reslimit& limit);
         ~solver();
-
+        const auto& monics_with_changed_bounds() const { return m_core->monics_with_changed_bounds(); }
         void add_monic(lpvar v, unsigned sz, lpvar const* vs);
         void add_idivision(lpvar q, lpvar x, lpvar y);
         void add_rdivision(lpvar q, lpvar x, lpvar y);
diff --git a/src/sat/smt/arith_solver.cpp b/src/sat/smt/arith_solver.cpp
index d64af9bb9..17ab03ee6 100644
--- a/src/sat/smt/arith_solver.cpp
+++ b/src/sat/smt/arith_solver.cpp
@@ -26,7 +26,7 @@ namespace arith {
         m_model_eqs(DEFAULT_HASHTABLE_INITIAL_CAPACITY, var_value_hash(*this), var_value_eq(*this)),
         m_local_search(*this),
         m_resource_limit(*this),
-        m_bp(*this),
+        m_bp(*this, m_implied_bounds),
         a(m),
         m_bound_terms(m),
         m_bound_predicate(m)
diff --git a/src/sat/smt/arith_solver.h b/src/sat/smt/arith_solver.h
index 6a577e753..801eb474c 100644
--- a/src/sat/smt/arith_solver.h
+++ b/src/sat/smt/arith_solver.h
@@ -243,6 +243,7 @@ namespace arith {
         resource_limit               m_resource_limit;
         lp_bounds                    m_new_bounds;
         symbol                       m_farkas;
+        std_vector<lp::implied_bound> m_implied_bounds;
         lp::lp_bound_propagator<solver> m_bp;
         mutable vector<std::pair<lp::tv, rational>> m_todo_terms;
 
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index b33688209..db4dba137 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -225,6 +225,7 @@ class theory_lra::imp {
     lp_bounds                    m_new_bounds;
     symbol                       m_farkas;
     vector<parameter>            m_bound_params;
+    std_vector<lp::implied_bound>   m_implied_bounds;
     lp::lp_bound_propagator<imp> m_bp;
 
     context& ctx() const { return th.get_context(); }
@@ -873,7 +874,7 @@ public:
         m_solver(nullptr),
         m_resource_limit(*this),
         m_farkas("farkas"),
-        m_bp(*this),
+        m_bp(*this, m_implied_bounds),
         m_bounded_range_idx(0),
         m_bounded_range_lit(null_literal),
         m_bound_terms(m),
@@ -2113,12 +2114,13 @@ public:
         m_model_is_initialized = false;
         flush_bound_axioms();
         // disabled in master:
-        // propagate_nla();
+        propagate_nla(); 
         if (ctx().inconsistent())
             return true;
         
-        if (!can_propagate_core())
+        if (!can_propagate_core()) 
             return false;
+        
         m_new_def = false;        
         while (m_asserted_qhead < m_asserted_atoms.size() && !ctx().inconsistent() && m.inc()) {
             auto [bv, is_true] = m_asserted_atoms[m_asserted_qhead];
@@ -2151,6 +2153,7 @@ public:
             propagate_bounds_with_lp_solver();
             break;
         case l_undef:
+            UNREACHABLE();
             break;
         }
         return true;            
@@ -2161,6 +2164,7 @@ public:
             m_nla->propagate();
             add_lemmas();
             add_equalities();
+            lp().collect_more_rows_for_lp_propagation();
         }
     }
 
@@ -2225,13 +2229,11 @@ public:
             // verbose_stream() << "unsat\n";
         }
         else {
-            unsigned count = 0, prop = 0;
             for (auto& ib : m_bp.ibounds()) {
                 m.inc();
                 if (ctx().inconsistent())
                     break;
-                ++prop;
-                count += propagate_lp_solver_bound(ib);
+                propagate_lp_solver_bound(ib);
             }
         }
     }
diff --git a/src/test/lp/nla_solver_test.cpp b/src/test/lp/nla_solver_test.cpp
index 054c6583f..fa5a29e99 100644
--- a/src/test/lp/nla_solver_test.cpp
+++ b/src/test/lp/nla_solver_test.cpp
@@ -193,9 +193,9 @@ void test_basic_lemma_for_mon_neutral_from_factors_to_monomial_0() {
 
     
     VERIFY(nla.get_core().test_check() == l_false);
-    
-    auto const& lv = nla.get_core().lemmas();
-    nla.get_core().print_lemma(lv.back(), std::cout);
+    auto const& lemmas = nla.get_core().lemmas();
+    nla.get_core().print_lemma(lemmas.back(), std::cout);
+
 
     ineq i0(lp_ac, llc::NE, 1);
     lp::lar_term t1, t2;
@@ -208,7 +208,7 @@ void test_basic_lemma_for_mon_neutral_from_factors_to_monomial_0() {
     bool found0 = false;
     bool found1 = false;
     bool found2 = false;
-    for (const auto& k : lv[0].ineqs()){
+    for (const auto& k : lemmas[0].ineqs()){
         if (k == i0) {
             found0 = true;
         } else if (k == i1) {
@@ -250,6 +250,7 @@ void test_basic_lemma_for_mon_neutral_from_factors_to_monomial_1() {
     svector<lpvar> v; v.push_back(lp_b);v.push_back(lp_d);v.push_back(lp_e);
     nla.add_monic(lp_bde, v.size(), v.begin());
 
+
     s_set_column_value_test(s, lp_a, rational(1));
     s_set_column_value_test(s, lp_b, rational(1));
     s_set_column_value_test(s, lp_c, rational(1));
@@ -395,6 +396,7 @@ void test_basic_lemma_for_mon_zero_from_monomial_to_factors() {
 
     VERIFY(nla.get_core().test_check() == l_false);
     auto const& lemma = nla.get_core().lemmas();
+
     nla.get_core().print_lemma(lemma.back(), std::cout);
 
     ineq i0(lp_a, llc::EQ, 0);
@@ -511,7 +513,8 @@ void test_horner() {
     
     reslimit l;
     params_ref p;
-    solver nla(s, p, l);
+    std_vector<lp::implied_bound> ib;
+    solver nla(s, p, l, ib);
     vector<lpvar> v;
     v.push_back(a); v.push_back(b);
     nla.add_monic(lp_ab, v.size(), v.begin());
@@ -583,6 +586,7 @@ void test_basic_sign_lemma() {
    
     VERIFY(nla.get_core().test_check() == l_false);
     auto const& lemmas = nla.get_core().lemmas();
+
     lp::lar_term t;
     t.add_var(lp_bde);
     t.add_var(lp_acd);
@@ -622,7 +626,8 @@ void test_order_lemma_params(bool var_equiv, int sign) {
     
     reslimit l;
     params_ref p;
-    solver nla(s,p,l);
+    std_vector<lp::implied_bound> ib;
+    solver nla(s,p,l,ib);
     // create monomial ab
     vector<unsigned> vec;
     vec.push_back(lp_a);
@@ -753,7 +758,8 @@ void test_monotone_lemma() {
     
     reslimit l;
     params_ref p;
-    solver nla(s, p, l);
+    std_vector<lp::implied_bound> ib;
+    solver nla(s, p, l, ib);
     // create monomial ab
     vector<unsigned> vec;
     vec.push_back(lp_a);
@@ -880,7 +886,8 @@ void test_tangent_lemma_equiv() {
     s_set_column_value_test(s, lp_a, - s.get_column_value(lp_k));
     reslimit l;
     params_ref p;
-    solver nla(s, p, l);
+    std_vector<lp::implied_bound> ib;
+    solver nla(s, p, l, ib);
     // create monomial ab
     vector<unsigned> vec;
     vec.push_back(lp_a);
diff --git a/src/util/heap.h b/src/util/heap.h
index df67b31be..332f4e53e 100644
--- a/src/util/heap.h
+++ b/src/util/heap.h
@@ -159,7 +159,7 @@ public:
     }
 
     unsigned size() const {
-        return m_value2indices.size();
+        return m_values.size() - 1;
     }
 
     void reserve(int s) {