diff --git a/src/math/interval/dep_intervals.h b/src/math/interval/dep_intervals.h
index f9768b6f4..f4774c186 100644
--- a/src/math/interval/dep_intervals.h
+++ b/src/math/interval/dep_intervals.h
@@ -175,6 +175,8 @@ public:
     void set_upper_is_inf(interval& a, bool inf) const { m_config.set_upper_is_inf(a, inf); }
     void set_lower_dep(interval& a, u_dependency* d) const { m_config.set_lower_dep(a, d); }
     void set_upper_dep(interval& a, u_dependency* d) const { m_config.set_upper_dep(a, d); }
+    u_dependency* get_lower_dep(interval const& a) const { return a.m_lower_dep; }
+    u_dependency* get_upper_dep(interval const& a) const { return a.m_upper_dep; }
     void reset(interval& a) const { set_lower_is_inf(a, true); set_upper_is_inf(a, true); }
     void set_value(interval& a, rational const& n) const { 
         set_lower(a, n); 
diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h
index 7ffffe5c5..362d8966a 100644
--- a/src/math/lp/lp_settings.h
+++ b/src/math/lp/lp_settings.h
@@ -115,7 +115,12 @@ struct statistics {
     unsigned m_hnf_cutter_calls;
     unsigned m_hnf_cuts;
     unsigned m_nla_calls;
-    unsigned m_nla_bounds;
+    unsigned m_nla_add_bounds;
+    unsigned m_nla_propagate_bounds;
+    unsigned m_nla_propagate_eq;
+    unsigned m_nla_lemmas;
+    unsigned m_nra_calls;
+    unsigned m_nla_bounds_improvements;
     unsigned m_horner_calls;
     unsigned m_horner_conflicts;
     unsigned m_cross_nested_forms;
@@ -145,7 +150,13 @@ struct statistics {
         st.update("arith-grobner-conflicts", m_grobner_conflicts);
         st.update("arith-offset-eqs", m_offset_eqs);
         st.update("arith-fixed-eqs", m_fixed_eqs);
-        st.update("arith-nla-bounds", m_nla_bounds);
+        st.update("arith-nla-add-bounds", m_nla_add_bounds);
+        st.update("arith-nla-propagate-bounds", m_nla_propagate_bounds);
+        st.update("arith-nla-propagate-eq", m_nla_propagate_eq);
+        st.update("arith-nla-lemmas", m_nla_lemmas);
+        st.update("arith-nra-calls", m_nra_calls);   
+        st.update("arith-bounds-improvements", m_nla_bounds_improvements);
+
     }
 };
 
diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp
index 81c00d2f0..46fa6db81 100644
--- a/src/math/lp/monomial_bounds.cpp
+++ b/src/math/lp/monomial_bounds.cpp
@@ -12,6 +12,8 @@
 #include "math/lp/nla_intervals.h"
 #include "math/lp/numeric_pair.h"
 
+#define UNIT_PROPAGATE_BOUNDS 0
+
 namespace nla {
 
     monomial_bounds::monomial_bounds(core* c):
@@ -19,10 +21,8 @@ namespace nla {
         dep(c->m_intervals.get_dep_intervals()) {}
 
     void monomial_bounds::propagate() {
-        for (lpvar v : c().m_to_refine) {
-            monic const& m = c().emons()[v];
-            propagate(m);
-        }
+        for (lpvar v : c().m_to_refine) 
+            propagate(c().emon(v));
     }
 
 
@@ -55,29 +55,41 @@ namespace nla {
     bool monomial_bounds::propagate_value(dep_interval& range, lpvar v) {
         auto val = c().val(v);
         if (dep.is_below(range, val)) {
+            auto const& upper = dep.upper(range);
+            auto cmp = dep.upper_is_open(range) ? llc::LT : llc::LE;
+            ++c().lra.settings().stats().m_nla_propagate_bounds;
+#if UNIT_PROPAGATE_BOUNDS
+            auto* d = dep.get_upper_dep(range);
+            c().lra.update_column_type_and_bound(v, cmp, upper, d);
+#else
             lp::explanation ex;
             dep.get_upper_dep(range, ex);
-            auto const& upper = dep.upper(range);
             if (is_too_big(upper))
                 return false;
-            auto cmp = dep.upper_is_open(range) ? llc::LT : llc::LE;
             new_lemma lemma(c(), "propagate value - upper bound of range is below value");
             lemma &= ex;
             lemma |= ineq(v, cmp, upper); 
-            TRACE("nla_solver", dep.display(tout << val << " > ", range) << "\n" << lemma << "\n";);
+            TRACE("nla_solver", dep.display(tout << c().val(v) << " > ", range) << "\n" << lemma << "\n";);
+#endif
             return true;
         }
         else if (dep.is_above(range, val)) {
+            auto const& lower = dep.lower(range);
+            auto cmp = dep.lower_is_open(range) ? llc::GT : llc::GE;
+            ++c().lra.settings().stats().m_nla_propagate_bounds;
+#if UNIT_PROPAGATE_BOUNDS
+            auto* d = dep.get_lower_dep(range);
+            c().lra.update_column_type_and_bound(v, cmp, lower, d);            
+#else
             lp::explanation ex;
             dep.get_lower_dep(range, ex);
-            auto const& lower = dep.lower(range);
             if (is_too_big(lower))
                 return false;
-            auto cmp = dep.lower_is_open(range) ? llc::GT : llc::GE;
             new_lemma lemma(c(), "propagate value - lower bound of range is above value");
             lemma &= ex;
             lemma |= ineq(v, cmp, lower); 
-            TRACE("nla_solver", dep.display(tout << val << " < ", range) << "\n" << lemma << "\n";);
+            TRACE("nla_solver", dep.display(tout << c().val(v) << " < ", range) << "\n" << lemma << "\n";);
+#endif
             return true;
         }
         else {
@@ -106,6 +118,7 @@ namespace nla {
             lp::explanation ex;
             dep.get_upper_dep(range, ex);
             if (p % 2 == 0 && rational(dep.upper(range)).is_neg()) {
+                ++c().lra.settings().stats().m_nla_propagate_bounds;
                 new_lemma lemma(c(), "range requires a non-negative upper bound");
                 lemma &= ex;
                 return true;
@@ -114,18 +127,30 @@ namespace nla {
                 // v = -2, [-4,-3]^3 < v^3 -> add bound v <= -3
                 // v = -2, [-1,+1]^2 < v^2 -> add bound v >= -1                
                 if ((p % 2 == 1) || val_v.is_pos()) {
+                    ++c().lra.settings().stats().m_nla_propagate_bounds;
                     auto le = dep.upper_is_open(range) ? llc::LT : llc::LE;
+#if UNIT_PROPAGATE_BOUNDS
+                    auto* d = dep.get_upper_dep();
+                    c().lra.update_column_type_and_bound(v, le, r, d);
+#else
                     new_lemma lemma(c(), "propagate value - root case - upper bound of range is below value");
                     lemma &= ex;
                     lemma |= ineq(v, le, r);
+#endif
                     return true;
                 }
                 if (p % 2 == 0 && val_v.is_neg()) {
+                    ++c().lra.settings().stats().m_nla_propagate_bounds;
                     SASSERT(!r.is_neg());
                     auto ge = dep.upper_is_open(range) ? llc::GT : llc::GE;
+#if UNIT_PROPAGATE_BOUNDS
+                    auto* d = dep.get_upper_dep();
+                    c().lra.update_column_type_and_bound(v, ge, -r, d);
+#else
                     new_lemma lemma(c(), "propagate value - root case - upper bound of range is below negative value");
                     lemma &= ex;
                     lemma |= ineq(v, ge, -r);
+#endif
                     return true;
                 }
             }
@@ -133,6 +158,7 @@ namespace nla {
         }
         else if (dep.is_above(range, val)) {
             if (rational(dep.lower(range)).root(p, r)) {
+                ++c().lra.settings().stats().m_nla_propagate_bounds;
                 lp::explanation ex;
                 dep.get_lower_dep(range, ex);
                 auto ge = dep.lower_is_open(range) ? llc::GT : llc::GE;
@@ -140,9 +166,8 @@ namespace nla {
                 new_lemma lemma(c(), "propagate value - root case - lower bound of range is above value");
                 lemma &= ex;
                 lemma |= ineq(v, ge, r); 
-                if (p % 2 == 0) {
+                if (p % 2 == 0) 
                     lemma |= ineq(v, le, -r);                 
-                }
                 return true;
             }
             // TBD: add bounds as long as difference to val is above some epsilon.
@@ -261,8 +286,10 @@ namespace nla {
 
     void monomial_bounds::unit_propagate() {        
         for (lpvar v : c().m_monics_with_changed_bounds) {
-            if (!c().is_monic_var(v)) continue;
-            unit_propagate(c().emons()[v]);
+            if (!c().is_monic_var(v))
+                continue;
+            monic& m = c().emon(v);
+            unit_propagate(m);
             if (c().lra.get_status() == lp::lp_status::INFEASIBLE) {
                 lp::explanation exp;
                 c().lra.get_infeasibility_explanation(exp);
@@ -270,9 +297,8 @@ namespace nla {
                 lemma &= exp;
                 break;
             }
-            if (c().m_conflicts > 0 ) {
+            if (c().m_conflicts > 0) 
                 break;
-            }
         }   
     }
 
@@ -280,8 +306,12 @@ namespace nla {
         if (m.is_propagated())
             return;
 
-        if (!is_linear(m))
+        if (!is_linear(m)) {
+#if UNIT_PROPAGATE_BOUNDS
+            propagate(m);
+#endif            
             return;
+        }
 
         c().emons().set_propagated(m);
 
@@ -291,6 +321,7 @@ namespace nla {
             propagate_fixed(m, k);
         else
             propagate_nonfixed(m, k, w);
+        ++c().lra.settings().stats().m_nla_propagate_eq;
     }
 
     lp::explanation monomial_bounds::get_explanation(u_dependency* dep) {
@@ -317,7 +348,6 @@ namespace nla {
     }
 
     void monomial_bounds::propagate_nonfixed(monic const& m, rational const& k, lpvar w) {
-        VERIFY(k != 0);
         vector<std::pair<lp::mpq, unsigned>> coeffs;        
         coeffs.push_back(std::make_pair(-k, w));
         coeffs.push_back(std::make_pair(rational::one(), m.var()));
diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h
index b2079b4f1..6fa439e01 100644
--- a/src/math/lp/monomial_bounds.h
+++ b/src/math/lp/monomial_bounds.h
@@ -17,7 +17,6 @@ namespace nla {
     class monomial_bounds : common {
         dep_intervals& dep;
 
-
         void var2interval(lpvar v, scoped_dep_interval& i);
         bool is_too_big(mpq const& q) const;
         bool propagate_down(monic const& m, lpvar u);
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index 2a23ba47b..056df3b69 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -45,11 +45,8 @@ core::core(lp::lar_solver& s, params_ref const& p, reslimit & lim) :
         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());
-                }
-            }    
+            for (const auto & m: m_emons.get_use_list(j)) 
+                m_monics_with_changed_bounds.insert(m.var());
         }
     };
 }
@@ -1257,7 +1254,7 @@ bool core::var_breaks_correct_monic_as_factor(lpvar j, const monic& m) const {
 
 bool core::var_breaks_correct_monic(lpvar j) const {
     if (is_monic_var(j) && !m_to_refine.contains(j)) {
-        TRACE("nla_solver", tout << "j = " << j << ", m  = "; print_monic(emons()[j], tout) << "\n";);
+        TRACE("nla_solver", tout << "j = " << j << ", m  = "; print_monic(emon(j), tout) << "\n";);
         return true; // changing the value of a correct monic
     }
     
@@ -1279,7 +1276,7 @@ void core::update_to_refine_of_var(lpvar j) {
             insert_to_refine(var(m));
     }
     if (is_monic_var(j)) {
-        const monic& m = emons()[j];
+        const monic& m = emon(j);
         if (var_val(m) == mul_val(m))
             erase_from_to_refine(j);
         else
@@ -1362,10 +1359,10 @@ bool in_power(const svector<lpvar>& vs, unsigned l) {
 bool core::to_refine_is_correct() const {
     for (unsigned j = 0; j < lra.number_of_vars(); j++) {
         if (!is_monic_var(j)) continue;
-        bool valid = check_monic(emons()[j]);
+        bool valid = check_monic(emon(j));
         if (valid == m_to_refine.contains(j)) {
             TRACE("nla_solver", tout << "inconstency in m_to_refine : ";
-                  print_monic(emons()[j], tout) << "\n";
+                  print_monic(emon(j), tout) << "\n";
                   if (valid) tout << "should NOT be in to_refine\n";
                   else tout << "should be in to_refine\n";);
             return false;
@@ -1375,7 +1372,7 @@ bool core::to_refine_is_correct() const {
 }
 
 void core::patch_monomial(lpvar j) {    
-    m_patched_monic =& (emons()[j]);
+    m_patched_monic =& (emon(j));
     m_patched_var = j;
     TRACE("nla_solver", tout << "m = "; print_monic(*m_patched_monic, tout) << "\n";);
     rational v = mul_val(*m_patched_monic);
@@ -1530,7 +1527,7 @@ void core::add_bounds() {
             if (!var_is_free(j)) continue;
             // split the free variable (j <= 0, or j > 0), and return
             m_literals.push_back(ineq(j, lp::lconstraint_kind::EQ, rational::zero()));  
-            ++lp_settings().stats().m_nla_bounds;
+            ++lp_settings().stats().m_nla_add_bounds;
             return;
         }
     }    
@@ -1611,15 +1608,13 @@ lbool core::check() {
 
     if (no_effect() && params().arith_nl_nra() && ret == l_undef) {
         ret = m_nra.check();
-        m_stats.m_nra_calls++;
+        lp_settings().stats().m_nra_calls++;
     }
     
     if (ret == l_undef && !m_lemmas.empty() && m_reslim.inc()) 
         ret = l_false;
 
-    m_stats.m_nla_lemmas += m_lemmas.size();
-    for (const auto& l : m_lemmas)
-        m_stats.m_nla_explanations += static_cast<unsigned>(l.expl().size());
+    lp_settings().stats().m_nla_lemmas += m_lemmas.size();
 
     
     TRACE("nla_solver", tout << "ret = " << ret << ", lemmas count = " << m_lemmas.size() << "\n";);
@@ -1649,7 +1644,7 @@ lbool core::bounded_nlsat() {
     }
     p.set_uint("max_conflicts", UINT_MAX);           
     m_nra.updt_params(p);
-    m_stats.m_nra_calls++;
+    lp_settings().stats().m_nra_calls++;
     if (ret == l_undef) 
         ++m_nlsat_delay;    
     else { 
@@ -1678,7 +1673,7 @@ lbool core::test_check() {
 }
 
 std::ostream& core::print_terms(std::ostream& out) const {
-    for (unsigned i = 0; i< lra.terms().size(); i++) {
+    for (unsigned i = 0; i < lra.terms().size(); i++) {
         unsigned ext = lp::tv::mask_term(i);
         if (!lra.var_is_registered(ext)) {
             out << "term is not registered\n";
@@ -1798,10 +1793,6 @@ void core::set_use_nra_model(bool m) {
 }
     
 void core::collect_statistics(::statistics & st) {
-    st.update("arith-nla-explanations", m_stats.m_nla_explanations);
-    st.update("arith-nla-lemmas", m_stats.m_nla_lemmas);
-    st.update("arith-nra-calls", m_stats.m_nra_calls);   
-    st.update("arith-bounds-improvements", m_stats.m_bounds_improvements);
 }
 
 bool core::improve_bounds() {
@@ -1814,9 +1805,9 @@ bool core::improve_bounds() {
             return;
         seen.insert(v);
         if (lra.improve_bound(v, false))
-            bounds_improved = true, m_stats.m_bounds_improvements++;
+            bounds_improved = true, lp_settings().stats().m_nla_bounds_improvements++;
         if (lra.improve_bound(v, true))
-            bounds_improved = true, m_stats.m_bounds_improvements++;
+            bounds_improved = true, lp_settings().stats().m_nla_bounds_improvements++;
     };
     for (auto & m : m_emons) {
         insert(m.var());
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index ae76f6d5a..3ff33a879 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -60,18 +60,6 @@ class core {
     friend class nra::solver;
     friend class divisions;
 
-    struct stats {
-        unsigned m_nla_explanations;
-        unsigned m_nla_lemmas;
-        unsigned m_nra_calls;
-        unsigned m_bounds_improvements;
-        stats() { reset(); }
-        void reset() {
-            memset(this, 0, sizeof(*this));
-        }
-    };
-
-    stats    m_stats;
     unsigned m_nlsat_delay = 50;
     unsigned m_nlsat_fails = 0;
 
@@ -139,6 +127,8 @@ public:
     reslimit& reslim() { return m_reslim; }  
     emonics& emons() { return m_emons; }
     const emonics& emons() const { return m_emons; }
+    monic& emon(unsigned i) { return m_emons[i]; }
+    monic const& emon(unsigned i) const { return m_emons[i]; }
 
     bool has_relevant_monomial() const;