diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index 43acce576..d47447b55 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -204,6 +204,12 @@ namespace lp {
             return m_status;
         
         solve_with_core_solver();
+        if (m_status != lp_status::INFEASIBLE) {
+            if (m_settings.bound_propagation())
+                detect_rows_with_changed_bounds();
+        }
+
+        clear_columns_with_changed_bounds();
         return m_status;
     }
 
@@ -478,7 +484,6 @@ namespace lp {
         x = new_val;
         TRACE("lar_solver_feas", tout << "setting " << j << " to "
          << new_val << (column_is_feasible(j)?"feas":"non-feas") << "\n";);
-        m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
         change_basic_columns_dependend_on_a_given_nb_column(j, delta);
     }
 
@@ -786,7 +791,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() {
@@ -1106,6 +1112,7 @@ namespace lp {
 
     mpq lar_solver::get_value(column_index const& j) const {
         SASSERT(get_status() == lp_status::OPTIMAL || get_status() == lp_status::FEASIBLE);
+        SASSERT(m_columns_with_changed_bounds.empty());
         numeric_pair<mpq> const& rp = get_column_value(j);
         return from_model_in_impq_to_mpq(rp);        
     }
@@ -1822,8 +1829,7 @@ namespace lp {
 
         if (is_base(j) && column_is_fixed(j))
             m_fixed_base_var_set.insert(j);
-        track_column_feasibility(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") << " val = " << get_column_value(j) << std::endl;);   
+        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;);    
     }
 
     void lar_solver::insert_to_columns_with_changed_bounds(unsigned j) {
@@ -2370,10 +2376,12 @@ namespace lp {
         m_crossed_bounds_deps = m_dependencies.mk_join(bdep, dep);
     }
 
-    void lar_solver::track_column_feasibility(lpvar j) {
-        m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(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 9b0dd8381..298d1f879 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -89,13 +89,6 @@ class lar_solver : public column_namer {
     constraint_set m_constraints;
     // the set of column indices j such that bounds have changed for j
     indexed_uint_set m_columns_with_changed_bounds;
-public: 
-    const indexed_uint_set& columns_with_changed_bounds() const { return m_columns_with_changed_bounds; }
-    inline void clear_columns_with_changed_bounds() { m_columns_with_changed_bounds.reset(); }
-    void track_column_feasibility(lpvar j);
-
-private:
-    // m_touched_rows contains rows that have been changed by a pivoting or have a column with changed bounds
     indexed_uint_set m_touched_rows;
     unsigned_vector m_row_bounds_to_replay;
     u_dependency_manager m_dependencies;
@@ -145,22 +138,23 @@ private:
     void add_row_from_term_no_constraint(const lar_term* term, unsigned term_ext_index);
     void add_basic_var_to_core_fields();
     bool compare_values(impq const& lhs, lconstraint_kind k, const mpq& rhs);
-public:   
+
+    inline void clear_columns_with_changed_bounds() { m_columns_with_changed_bounds.reset(); }
+ public:   
     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; } 
-    bool current_x_is_feasible() const { return m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible(); }     
-    
+        
 
-private:   
+ 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);
@@ -364,8 +358,6 @@ private:
     void add_column_rows_to_touched_rows(lpvar j);
     template <typename T>
     void propagate_bounds_for_touched_rows(lp_bound_propagator<T>& bp) {
-        detect_rows_with_changed_bounds();
-        clear_columns_with_changed_bounds();
         if (settings().propagate_eqs()) {
             if (settings().random_next() % 10 == 0) 
                 remove_fixed_vars_from_base();
@@ -386,7 +378,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++)
@@ -696,6 +688,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/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp
index 0e4a5041e..95f4586d8 100644
--- a/src/math/lp/monomial_bounds.cpp
+++ b/src/math/lp/monomial_bounds.cpp
@@ -260,8 +260,9 @@ 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);
@@ -276,14 +277,15 @@ namespace nla {
     }
 
 
-    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().emons().set_propagated(m);
+
         rational k = fixed_var_product(m);
         lpvar w = non_fixed_var(m);
         if (w == null_lpvar || k == 0) {
diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h
index 747aca9a2..a65ca7f4a 100644
--- a/src/math/lp/monomial_bounds.h
+++ b/src/math/lp/monomial_bounds.h
@@ -35,7 +35,7 @@ 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);
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index 6a241e951..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() {
@@ -1817,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 5a597ae67..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;
@@ -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/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 7f7951b20..e854db802 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -2161,7 +2161,7 @@ public:
             m_nla->propagate();
             add_lemmas();
             add_equalities();
-            propagate_bounds_with_lp_solver();
+            lp().collect_more_rows_for_lp_propagation();
         }
     }
 
@@ -2212,10 +2212,6 @@ public:
     }
     
     void propagate_bounds_with_lp_solver() {
-        if (!lp().current_x_is_feasible()) {  
-            lp().clear_columns_with_changed_bounds();
-            return;
-        }
         m_bp.init();
         lp().propagate_bounds_for_touched_rows(m_bp);