diff --git a/.clang-format b/.clang-format
index 7ef241e3d..05b4679ad 100644
--- a/.clang-format
+++ b/.clang-format
@@ -8,6 +8,7 @@ BasedOnStyle: LLVM
 IndentWidth: 4
 TabWidth: 4
 UseTab: Never
+IndentCaseLabels: false
 
 # Column width
 ColumnLimit: 120
@@ -55,6 +56,7 @@ BinPackArguments: true
 BinPackParameters: true
 BreakBeforeBinaryOperators: None
 BreakBeforeTernaryOperators: true
+# BreakBeforeElse: true
 
 # Includes
 SortIncludes: false  # Z3 has specific include ordering conventions
diff --git a/src/math/lp/.clang-format b/src/math/lp/.clang-format
deleted file mode 100644
index e9420012c..000000000
--- a/src/math/lp/.clang-format
+++ /dev/null
@@ -1,5 +0,0 @@
-BasedOnStyle: Google
-IndentWidth: 4
-ColumnLimit: 0
-NamespaceIndentation: All
-BreakBeforeBraces: Attach
\ No newline at end of file
diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt
index 7dea022b2..0b6f03b6d 100644
--- a/src/math/lp/CMakeLists.txt
+++ b/src/math/lp/CMakeLists.txt
@@ -24,6 +24,7 @@ z3_add_component(lp
     monomial_bounds.cpp
     nex_creator.cpp
     nla_basics_lemmas.cpp
+    nla_coi.cpp
     nla_common.cpp
     nla_core.cpp
     nla_divisions.cpp
diff --git a/src/math/lp/nla_coi.cpp b/src/math/lp/nla_coi.cpp
new file mode 100644
index 000000000..e4e2e6017
--- /dev/null
+++ b/src/math/lp/nla_coi.cpp
@@ -0,0 +1,78 @@
+
+#include "math/lp/nla_core.h"
+#include "math/lp/nla_coi.h"
+
+namespace nla {
+
+  void coi::init() {
+      indexed_uint_set visited;
+      unsigned_vector todo;
+      vector<occurs> var2occurs;
+      m_term_set.reset();
+      m_mon_set.reset();
+      m_constraint_set.reset();
+      auto& lra = c.lra_solver();
+
+      for (auto ci : lra.constraints().indices()) {
+          auto const& c = lra.constraints()[ci];
+          if (c.is_auxiliary())
+              continue;
+          for (auto const& [coeff, v] : c.coeffs()) {
+              var2occurs.reserve(v + 1);
+              var2occurs[v].constraints.push_back(ci);
+          }
+      }
+
+      for (auto const& m : c.emons()) {
+          for (auto v : m.vars()) {
+              var2occurs.reserve(v + 1);
+              var2occurs[v].monics.push_back(m.var());
+          }
+      }
+
+      for (const auto *t :  lra.terms() ) {
+          for (auto const iv : *t) {
+              auto v = iv.j();
+              var2occurs.reserve(v + 1);
+              var2occurs[v].terms.push_back(t->j());
+          }
+      }
+
+      for (auto const& m : c.to_refine())
+          todo.push_back(m);
+
+      for (unsigned i = 0; i < todo.size(); ++i) {
+          auto v = todo[i];
+          if (visited.contains(v))
+              continue;
+          visited.insert(v);
+          var2occurs.reserve(v + 1);
+          for (auto ci : var2occurs[v].constraints) {
+              m_constraint_set.insert(ci);
+              auto const& c = lra.constraints()[ci];
+              for (auto const& [coeff, w] : c.coeffs())
+                  todo.push_back(w);
+          }
+          for (auto w : var2occurs[v].monics)
+              todo.push_back(w);
+
+          for (auto ti : var2occurs[v].terms) {
+              for (auto iv : lra.get_term(ti))
+                  todo.push_back(iv.j());
+              todo.push_back(ti);
+          }
+
+          if (lra.column_has_term(v)) {
+              m_term_set.insert(v);
+              for (auto kv : lra.get_term(v))
+                  todo.push_back(kv.j());
+          }            
+
+          if (c.is_monic_var(v)) {
+              m_mon_set.insert(v);
+              for (auto w : c.emons()[v])
+                  todo.push_back(w);
+          }
+      }
+  }
+}
\ No newline at end of file
diff --git a/src/math/lp/nla_coi.h b/src/math/lp/nla_coi.h
new file mode 100644
index 000000000..7e12da50b
--- /dev/null
+++ b/src/math/lp/nla_coi.h
@@ -0,0 +1,29 @@
+
+#pragma once
+
+namespace nla {
+
+    class core;
+
+    class coi {
+        core& c;
+        indexed_uint_set m_mon_set, m_constraint_set;
+        indexed_uint_set m_term_set;
+
+        struct occurs {
+            unsigned_vector constraints;
+            unsigned_vector monics;
+            unsigned_vector terms;
+        };
+
+    public:
+        coi(core& c) : c(c) {}
+
+        void init();        
+
+        indexed_uint_set const& mons() const { return m_mon_set; }
+        indexed_uint_set const& constraints() const { return m_constraint_set; }
+        indexed_uint_set& terms() { return m_term_set; }
+
+    };
+}
\ No newline at end of file
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index b74282045..da3f60f7f 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -1332,6 +1332,9 @@ lbool core::check() {
         if (!m_lemmas.empty() || !m_literals.empty() || m_check_feasible)
             return l_false;
     }
+
+    if (false && no_effect())
+        ret = m_mul_saturate.saturate();
     
     if (no_effect() && should_run_bounded_nlsat()) 
         ret = bounded_nlsat();
@@ -1348,8 +1351,7 @@ lbool core::check() {
     if (no_effect())
         m_order.order_lemma();
 
-    if (false && no_effect())
-        ret = m_mul_saturate.saturate();
+
 
     if (no_effect()) {
         unsigned num_calls = lp_settings().stats().m_nla_calls;
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index a2b6327fa..ab02c5d17 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -70,7 +70,7 @@ class core {
     lbool bounded_nlsat();
 
     var_eqs<emonics>         m_evars;
-
+    
     lp::lar_solver&          lra;
     reslimit&                m_reslim;
     smt_params_helper        m_params;
@@ -127,6 +127,9 @@ 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; }
+    lp::lar_solver& lra_solver() { return lra; }
+    indexed_uint_set const& to_refine() const { return m_to_refine; }
+
     void insert_to_refine(lpvar j);
     void erase_from_to_refine(lpvar j);
 
diff --git a/src/math/lp/nla_mul_saturate.cpp b/src/math/lp/nla_mul_saturate.cpp
index f0e18f777..e41f4491d 100644
--- a/src/math/lp/nla_mul_saturate.cpp
+++ b/src/math/lp/nla_mul_saturate.cpp
@@ -19,13 +19,13 @@
   --*/
 
 #include "math/lp/nla_core.h"
+#include "math/lp/nla_coi.h"
 #include "math/lp/nla_mul_saturate.h"
 
 
 namespace nla {
 
-    mul_saturate::mul_saturate(core* core) : 
-        common(core) {}
+    mul_saturate::mul_saturate(core* core) : common(core), m_coi(*core) {}
 
     lbool mul_saturate::saturate() {
         lp::explanation ex;
@@ -39,26 +39,81 @@ namespace nla {
 
     void mul_saturate::init_solver() {
         local_solver = alloc(lp::lar_solver);
+        m_vars2mon.reset();
+        m_mon2vars.reset();
+        m_values.reset();
+        m_coi.init();
+        init_vars();
+    }
+
+    void mul_saturate::init_vars() {
+        auto const& lra = c().lra_solver();
+        auto sz = lra.number_of_vars();
+        for (unsigned v = 0; v < sz; ++v) {
+            unsigned w;
+            if (m_coi.mons().contains(v))
+                init_monomial(v);
+            else 
+                m_values.push_back(c().val(v));
+            if (m_coi.terms().contains(v)) {
+                auto const& t = lra.get_term(v);
+                // Assumption: variables in coefficients are always declared before term variable.
+                SASSERT(all_of(t, [&](auto p) { return p.j() < v; }));
+                w = local_solver->add_term(t.coeffs_as_vector(), v);
+            } 
+            else
+                w = local_solver->add_var(v, lra.var_is_int(v));
+
+            VERIFY(w == v);
+            if (lra.column_has_lower_bound(v)) {
+                auto lo_dep = lra.get_column_lower_bound_witness(v);
+                auto lo_bound = lra.get_lower_bound(v);
+                auto k = lo_bound.y > 0 ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE;
+                auto ci = local_solver->add_var_bound(v, k, lo_bound.x);
+            }
+            if (lra.column_has_upper_bound(v)) {
+                auto hi_dep = lra.get_column_upper_bound_witness(v);
+                auto hi_bound = lra.get_upper_bound(v);
+                auto k = hi_bound.y < 0 ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE;
+                auto ci = local_solver->add_var_bound(v, k, hi_bound.x);
+                m_ci2dep.setx(ci, hi_dep, nullptr);
+            }
+        }
+    }
+
+    void mul_saturate::init_monomial(unsigned mon_var) {
+        auto& mon = c().emons()[mon_var];
+        svector<lpvar> vars(mon.vars());
+        std::sort(vars.begin(), vars.end());
+        m_vars2mon.insert(vars, mon_var);
+        m_mon2vars.insert(mon_var, vars);
+        rational p(1);
+        for (auto v : vars)
+            p *= m_values[v];
+        m_values.push_back(p);
     }
 
     void mul_saturate::add_lemma(lp::explanation const& ex1) {
+        auto& lra = c().lra_solver();
         lp::explanation ex2;
         for (auto p : ex1) {
             lp::constraint_index src_ci;
-            if (m_new_mul_constraints.find(p.ci(), src_ci))
-                ex2.add_pair(src_ci, mpq(1));
-            else
-                ex2.add_pair(p.ci(), p.coeff());
+            if (!m_new_mul_constraints.find(p.ci(), src_ci))
+                src_ci = p.ci();
+            auto dep = m_ci2dep.get(src_ci, nullptr);
+            local_solver->push_explanation(dep, ex2);
+        }
+        for (auto [v, sign, dep] : m_var_signs) {            
+            if (!dep) {
+                verbose_stream() << "TODO explain non-implied bound\n";
+                continue;
+            }
+            local_solver->push_explanation(dep, ex2);
         }
         lemma_builder new_lemma(c(), "stellensatz");
         new_lemma &= ex2;
-        for (auto [v, sign] : m_var_signs) {
-            if (sign)
-                new_lemma.explain_existing_upper_bound(v);
-            else
-                new_lemma.explain_existing_lower_bound(v);
-        }
         IF_VERBOSE(1, verbose_stream() << "unsat \n" << new_lemma << "\n");
+        TRACE(arith, tout << "unsat\n" << new_lemma << "\n");
     }
 
     lbool mul_saturate::solve(lp::explanation& ex) {
@@ -68,49 +123,74 @@ namespace nla {
         lbool r = l_undef;
         if (st == lp::lp_status::INFEASIBLE) {
             local_solver->get_infeasibility_explanation(ex);
-            IF_VERBOSE(0, c().print_explanation(ex, verbose_stream()) << "\n";);
             r = l_false;
         }
         if (st == lp::lp_status::OPTIMAL || st == lp::lp_status::FEASIBLE) {
             // TODO: check model just in case it got lucky.
             IF_VERBOSE(1, verbose_stream() << "saturation is LP feasible, maybe it is a model of the NLA problem\n");
         }
-        IF_VERBOSE(0, local_solver->display(verbose_stream()); c().display(verbose_stream()));
+        IF_VERBOSE(0, display(verbose_stream()));
         return r;
     }
 
     // record new monomials that are created and recursively down-saturate with respect to these.
+    // this is a simplistic pass
     void mul_saturate::add_multiply_constraints() {
         m_new_mul_constraints.reset();
         m_seen_vars.reset();
         m_var_signs.reset();
-        for (auto j : c().m_to_refine) {
-            for (auto con_id : local_solver->constraints().indices()) {
-                unsigned num_vars = c().emon(j).vars().size();
-                for (unsigned i = 0; i < num_vars; ++i) {
-                    auto v = c().emon(j).vars()[i];
-                    for (auto [coeff, u] : local_solver->constraints()[con_id].coeffs())
+        m_to_refine.reset();
+        vector<svector<lp::constraint_index>> var2cs;
+
+        for (auto ci : local_solver->constraints().indices()) {
+            auto const& con = local_solver->constraints()[ci];
+            for (auto [coeff, v] : con.coeffs()) {
+                if (v >= var2cs.size())
+                    var2cs.resize(v + 1);
+                var2cs[v].push_back(ci);
+            }
+
+            // insert monomials to be refined
+            insert_monomials_from_constraint(ci);            
+        }
+
+        for (unsigned it = 0; it < m_to_refine.size(); ++it) {
+            auto j = m_to_refine[it];
+            verbose_stream() << "refining " << j << " := " << m_mon2vars[j] << "\n";
+            auto vars = m_mon2vars[j];
+            for (auto v : vars) {
+                if (v >= var2cs.size())
+                    continue;
+                auto cs = var2cs[v];
+                for (auto ci : cs) {
+                    for (auto [coeff, u] : local_solver->constraints()[ci].coeffs()) {
                         if (u == v)
-                            add_multiply_constraint(con_id, j, v);
+                            add_multiply_constraint(ci, j, v);
+                    }                    
                 }
             }
         }
+        IF_VERBOSE(0,
+                   c().lra_solver().display(verbose_stream() << "original\n");
+                   c().display(verbose_stream());
+                   display(verbose_stream() << "saturated\n"));
     }
 
     // multiply by remaining vars
     void mul_saturate::add_multiply_constraint(lp::constraint_index old_ci, lp::lpvar mi, lpvar x) {
         lp::lar_base_constraint const& con = local_solver->constraints()[old_ci];
+        auto &lra = c().lra_solver();
         auto const& lhs = con.coeffs();
         auto const& rhs = con.rhs();
-        auto k = con.kind();        
+        auto k = con.kind();
         if (k == lp::lconstraint_kind::NE || k == lp::lconstraint_kind::EQ)
             return;  // not supported
         auto sign = false;
         svector<lpvar> vars;
         bool first = true;
         for (auto v : c().emon(mi).vars()) {
-            if (v != x || !first) 
-                vars.push_back(v);            
+            if (v != x || !first)
+                vars.push_back(v);
             else
                 first = false;
         }
@@ -118,64 +198,186 @@ namespace nla {
         for (auto v : vars) {
             if (m_seen_vars.contains(v))
                 continue;
+            m_seen_vars.insert(v);
             // retrieve bounds of v
-            // if v has non-negative lower bound add as positive
-            // if v has non-positive upper bound add as negative
-            // otherwise, fail
-            if (local_solver->column_has_lower_bound(v) && !local_solver->get_lower_bound(v).is_neg()) {
-                m_var_signs.push_back({v, false});
-                m_seen_vars.insert(v);
+            // if v has positive lower bound add as positive
+            // if v has negative upper bound add as negative
+            // otherwise, soft-fail (for now unsound)
+            // proper signs of variables from old tableau should be extracted using lra_solver()
+            // instead of local_solver.
+            // TODO is to also add case where lower or upper bound is zero and then the sign
+            // of the inequality is relaxed if it is strict.
+            if (lra.number_of_vars() > v && lra.column_has_lower_bound(v) && lra.get_lower_bound(v).is_pos()) {
+                m_var_signs.push_back({v, false, lra.get_column_lower_bound_witness(v)});
             } 
-            else if (local_solver->column_has_upper_bound(v) && !local_solver->get_upper_bound(v).is_pos()) {
-                m_var_signs.push_back({v, true});
-                m_seen_vars.insert(v);
+            else if (lra.number_of_vars() > v && lra.column_has_upper_bound(v) && lra.get_upper_bound(v).is_neg()) {
+                m_var_signs.push_back({v, true, lra.get_column_upper_bound_witness(v)});
                 sign = !sign;
-            }
-            else 
-                return;
-        }
-        lp::lar_term new_lhs;        
-        rational new_rhs(rhs);
-        for (auto [coeff, v] : lhs) {
-            #if 0
-            vars.push_back(v);
-            lpvar new_monic_var = c().m_add_monomial(vars);
-            auto const& new_m = c().emons()[new_monic_var];            
-            verbose_stream() << vars << " v " << new_m.var() << " coeff " << coeff << "\n";
-            new_lhs.add_monomial(coeff, new_m.var());
-            vars.pop_back();
-            #endif
-        }
-        if (rhs != 0) {
-            if (vars.size() == 1) {
-                new_lhs.add_monomial(-rhs, vars[0]);
-                verbose_stream() << "rhs mul " << -rhs << " * j" << vars[0] << "\n";
+            } 
+            else if (m_values[v].is_neg()) {
+                m_var_signs.push_back({v, true, nullptr});
+                sign = !sign;
+            } 
+            else if (m_values[v].is_pos()) {
+                m_var_signs.push_back({v, false, nullptr});
             } 
             else {
-#if 0
-                lpvar new_monic_var = c().m_add_monomial(vars);
-                auto const& new_m = c().emons()[new_monic_var];
-                verbose_stream() << vars << " v " << new_m.var() << " coeff " << coeff << "\n";
-                new_lhs.add_monomial(-rhs, new_m.var());
-                verbose_stream() << "rhs mul " << -rhs << " * j" << new_m.var() << "\n";
-#endif
+                IF_VERBOSE(0, verbose_stream() << "not separated from 0\n");
+                return;
             }
+        }
+        lp::lar_term new_lhs;
+        rational new_rhs(rhs), term_value(0);
+        for (auto const & [coeff, v] : lhs) {
+            unsigned old_sz = vars.size();
+            if (m_mon2vars.contains(v))
+                vars.append(m_mon2vars[v]);
+            else
+                vars.push_back(v);
+            lpvar new_monic_var = add_monomial(vars);
+            new_lhs.add_monomial(coeff, new_monic_var);
+            term_value += coeff * m_values[new_monic_var];
+            vars.shrink(old_sz);
+        }
+        if (rhs != 0) {
+            lpvar new_monic_var = add_monomial(vars);
+            new_lhs.add_monomial(-rhs, new_monic_var);
+            term_value -= rhs * m_values[new_monic_var];
             new_rhs = 0;
         }
 
         if (sign) {
             switch (k) {
-            case lp::lconstraint_kind::LE: k = lp::lconstraint_kind::GE; break;
-            case lp::lconstraint_kind::LT: k = lp::lconstraint_kind::GT; break;
-            case lp::lconstraint_kind::GE: k = lp::lconstraint_kind::LE; break;
-            case lp::lconstraint_kind::GT: k = lp::lconstraint_kind::LT; break;
-            default: break;
+            case lp::lconstraint_kind::LE:
+                k = lp::lconstraint_kind::GE;
+                break;
+            case lp::lconstraint_kind::LT:
+                k = lp::lconstraint_kind::GT;
+                break;
+            case lp::lconstraint_kind::GE:
+                k = lp::lconstraint_kind::LE;
+                break;
+            case lp::lconstraint_kind::GT:
+                k = lp::lconstraint_kind::LT;
+                break;
+            default:
+                break;
             }
-        }   
-        c().display_constraint(verbose_stream(), old_ci) << " -> ";
-        c().display_constraint(verbose_stream(), new_lhs, k, new_rhs) << "\n";
-        // TODO: 
-        // auto new_ci = lra.m_add_constraint(new_lhs, k, new_rhs);
-        // m_new_mul_constraints.insert(new_ci, old_ci);        
+        }
+        display_constraint(verbose_stream(), old_ci) << " -> ";
+        display_constraint(verbose_stream(), new_lhs.coeffs_as_vector(), k, new_rhs) << "\n";
+        auto ti = local_solver->add_term(new_lhs.coeffs_as_vector(), local_solver->number_of_vars());
+        auto new_ci = local_solver->add_var_bound(ti, k, new_rhs);
+        insert_monomials_from_constraint(new_ci);
+        m_values.push_back(term_value);
+        SASSERT(m_values.size() - 1 == ti);
+        m_new_mul_constraints.insert(new_ci, old_ci);
+    }
+
+    lpvar mul_saturate::add_monomial(svector<lpvar> const& vars) {
+        lpvar v;
+        if (vars.size() == 1)
+            return vars[0];
+        svector<lpvar> _vars(vars);
+        std::sort(_vars.begin(), _vars.end());
+        if (m_vars2mon.find(_vars, v))
+            return v;
+
+        v = add_var(is_int(vars));
+        m_vars2mon.insert(_vars, v);
+        m_mon2vars.insert(v, _vars);
+        rational p(1);
+        for (auto v : vars)
+            p *= m_values[v];
+        m_values.push_back(p);
+        SASSERT(m_values.size() - 1 == v);
+        return v;
+    }
+
+    bool mul_saturate::is_int(svector<lp::lpvar> const& vars) const {
+        auto const& lra = m_core.lra;
+        return all_of(vars, [&](lpvar v) { return lra.var_is_int(v); });
+    }
+
+    lpvar mul_saturate::add_var(bool is_int) {
+        auto v = local_solver->number_of_vars();
+        auto w = local_solver->add_var(v, is_int);
+        VERIFY(v == w);
+        return w;
+    }
+
+    void mul_saturate::insert_monomials_from_constraint(lp::constraint_index ci) {
+        if (constraint_is_true(ci))
+            return;
+        auto const& con = local_solver->constraints()[ci];
+        for (auto [coeff, v] : con.coeffs()) 
+            if (c().is_monic_var(v))
+                m_to_refine.insert(v);        
+    }
+
+    bool mul_saturate::constraint_is_true(lp::constraint_index ci) {
+        auto const& con = local_solver->constraints()[ci];
+        auto lhs = -con.rhs();
+        for (auto const& [coeff, v] : con.coeffs())
+            lhs += coeff * m_values[v];
+        switch (con.kind()) {
+            case lp::lconstraint_kind::GT:
+                return lhs > 0;
+            case lp::lconstraint_kind::GE:
+                return lhs >= 0;
+            case lp::lconstraint_kind::LE:
+                return lhs <= 0;
+            case lp::lconstraint_kind::LT:
+                return lhs < 0;
+            case lp::lconstraint_kind::EQ:
+                return lhs == 0;
+            case lp::lconstraint_kind::NE:
+                return lhs != 0;
+            default:
+                UNREACHABLE();
+                return false;        
+        }
+    }
+
+    std::ostream& mul_saturate::display(std::ostream& out) const {
+        local_solver->display(out);
+        for (auto const& [vars, v] : m_vars2mon) {
+            out << "j" << v << " := ";
+            display_product(out, vars);
+            out << "\n";
+        }
+        return out;
+    }
+
+    std::ostream& mul_saturate::display_product(std::ostream& out, svector<lpvar> const& vars) const {
+        bool first = true;
+        for (auto v : vars) {
+            if (first)
+                first = false;
+            else
+                out << " * ";
+            out << "j" << v;
+        }
+        return out;
+    }
+
+    std::ostream& mul_saturate::display_constraint(std::ostream& out, lp::constraint_index ci) const {
+        auto const& con = local_solver->constraints()[ci];
+        return display_constraint(out, con.coeffs(), con.kind(), con.rhs());
+    }
+
+    std::ostream& mul_saturate::display_constraint(std::ostream& out, 
+        vector<std::pair<rational, lpvar>> const& lhs,
+        lp::lconstraint_kind k, rational const& rhs) const {
+        bool first = true;
+        for (auto [coeff, v] : lhs) {
+            c().display_coeff(out, first, coeff);
+            first = false;
+            if (m_mon2vars.contains(v))
+                display_product(out, m_mon2vars[v]);
+            else
+                out << "j" << v;
+        }
+        return out << " " << k << " " << rhs;
     }
 }
diff --git a/src/math/lp/nla_mul_saturate.h b/src/math/lp/nla_mul_saturate.h
index bcfa42277..0c5e09e1f 100644
--- a/src/math/lp/nla_mul_saturate.h
+++ b/src/math/lp/nla_mul_saturate.h
@@ -4,24 +4,65 @@
   --*/
 #pragma once
 
+#include "math/lp/nla_coi.h"
+
 namespace nla {
 
     class core;
-    class lar_solver;
+    class lar_solver;    
     class mul_saturate : common {
+
+        struct var_sign {
+            lpvar v = lp::null_lpvar;
+            bool is_neg = false;
+            u_dependency* dep = nullptr;
+        };
+        coi m_coi;
         // source of multiplication constraint
         u_map<lp::constraint_index> m_new_mul_constraints;
-        svector<std::pair<lpvar, bool>> m_var_signs;
+        svector<var_sign> m_var_signs;
         tracked_uint_set m_seen_vars;
+        indexed_uint_set m_to_refine;
         scoped_ptr<lp::lar_solver> local_solver;
+        ptr_vector<u_dependency> m_ci2dep;
+        vector<rational> m_values;
+        struct eq {
+            bool operator()(unsigned_vector const& a, unsigned_vector const& b) const {
+                return a == b;
+            }
+        };
+        map<unsigned_vector, unsigned, svector_hash<unsigned_hash>, eq> m_vars2mon;
+        u_map<unsigned_vector> m_mon2vars;
+
+        // initialization
         void init_solver();
+        void init_vars();
+        void init_monomial(unsigned mon_var);
+
+        bool constraint_is_true(lp::constraint_index ci);
+        void insert_monomials_from_constraint(lp::constraint_index ci);
+
+        // additional variables and monomials and constraints
+        lpvar add_monomial(svector<lp::lpvar> const& vars);
+        bool is_int(svector<lp::lpvar> const& vars) const;
+        lpvar add_var(bool is_int);
         void add_multiply_constraints();
         void add_multiply_constraint(lp::constraint_index con_id, lp::lpvar mi, lpvar x);
+
+        // solving
         lbool solve(lp::explanation& ex);
-        void add_lemma(lp::explanation const& ex1);
+
+        // lemmas
+        void add_lemma(lp::explanation const& ex);
+
+        std::ostream& display(std::ostream& out) const;
+        std::ostream& display_product(std::ostream& out, svector<lpvar> const& vars) const;
+        std::ostream& display_constraint(std::ostream& out, lp::constraint_index ci) const;
+        std::ostream& display_constraint(std::ostream& out, vector<std::pair<rational, lpvar>> const& lhs,
+                                         lp::lconstraint_kind k, rational const& rhs) const;
+
     public:
         mul_saturate(core* core);
-
         lbool saturate();
     };
 
diff --git a/src/math/lp/nla_pp.cpp b/src/math/lp/nla_pp.cpp
index 66e930f8c..199e81665 100644
--- a/src/math/lp/nla_pp.cpp
+++ b/src/math/lp/nla_pp.cpp
@@ -250,15 +250,17 @@ std::ostream& core::display_coeff(std::ostream& out, bool first, lp::mpq const&
     if (first && p == 1)
         return out;
     if (first && p > 0)
-        out << p;
+        out << p << " * ";
+    else if (first && p == -1)
+        out << "-";
+    else if (first)
+        out << p << " * ";
     else if (p == 1)
         out << " + ";
     else if (p > 0)
         out << " + " << p << " * ";
     else if (p == -1)
         out << " - ";
-    else if (first)
-        out << p << " * ";
     else
         out << " - " << -p << " * ";
     return out;
diff --git a/src/math/lp/nra_solver.cpp b/src/math/lp/nra_solver.cpp
index bc498f9e4..3a7f40093 100644
--- a/src/math/lp/nra_solver.cpp
+++ b/src/math/lp/nra_solver.cpp
@@ -9,6 +9,7 @@
 #include <fstream>
 #include "math/lp/lar_solver.h"
 #include "math/lp/nra_solver.h"
+#include "math/lp/nla_coi.h"
 #include "nlsat/nlsat_solver.h"
 #include "math/polynomial/polynomial.h"
 #include "math/polynomial/algebraic_numbers.h"
@@ -29,107 +30,30 @@ struct solver::imp {
     reslimit&                 m_limit;  
     params_ref                m_params; 
     u_map<polynomial::var>    m_lp2nl;  // map from lar_solver variables to nlsat::solver variables        
-    indexed_uint_set          m_term_set;
     scoped_ptr<nlsat::solver> m_nlsat;
     scoped_ptr<scoped_anum_vector>   m_values; // values provided by LRA solver
     scoped_ptr<scoped_anum> m_tmp1, m_tmp2;
+    nla::coi                  m_coi;
     nla::core&                m_nla_core;
     
     imp(lp::lar_solver& s, reslimit& lim, params_ref const& p, nla::core& nla_core): 
         lra(s), 
         m_limit(lim),
         m_params(p),
+        m_coi(nla_core),
         m_nla_core(nla_core) {}
 
     bool need_check() {
         return m_nla_core.m_to_refine.size() != 0;
     }
 
-    indexed_uint_set m_mon_set, m_constraint_set;
-
-    struct occurs {
-        unsigned_vector constraints;
-        unsigned_vector monics;
-        unsigned_vector terms;
-    };
-
-    void init_cone_of_influence() {
-        indexed_uint_set visited;
-        unsigned_vector todo;
-        vector<occurs> var2occurs;
-        m_term_set.reset();
-        m_mon_set.reset();
-        m_constraint_set.reset();
-
-        for (auto ci : lra.constraints().indices()) {
-            auto const& c = lra.constraints()[ci];
-            if (c.is_auxiliary())
-                continue;
-            for (auto const& [coeff, v] : c.coeffs()) {
-                var2occurs.reserve(v + 1);
-                var2occurs[v].constraints.push_back(ci);
-            }
-        }
-
-        for (auto const& m : m_nla_core.emons()) {
-            for (auto v : m.vars()) {
-                var2occurs.reserve(v + 1);
-                var2occurs[v].monics.push_back(m.var());
-            }
-        }
-
-        for (const auto *t :  lra.terms() ) {
-            for (auto const iv : *t) {
-                auto v = iv.j();
-                var2occurs.reserve(v + 1);
-                var2occurs[v].terms.push_back(t->j());
-            }
-        }
-
-        for (auto const& m : m_nla_core.m_to_refine)
-            todo.push_back(m);
-
-        for (unsigned i = 0; i < todo.size(); ++i) {
-            auto v = todo[i];
-            if (visited.contains(v))
-                continue;
-            visited.insert(v);
-            var2occurs.reserve(v + 1);
-            for (auto ci : var2occurs[v].constraints) {
-                m_constraint_set.insert(ci);
-                auto const& c = lra.constraints()[ci];
-                for (auto const& [coeff, w] : c.coeffs())
-                    todo.push_back(w);
-            }
-            for (auto w : var2occurs[v].monics)
-                todo.push_back(w);
-
-            for (auto ti : var2occurs[v].terms) {
-                for (auto iv : lra.get_term(ti))
-                    todo.push_back(iv.j());
-                todo.push_back(ti);
-            }
-
-            if (lra.column_has_term(v)) {
-                m_term_set.insert(v);
-                for (auto kv : lra.get_term(v))
-                    todo.push_back(kv.j());
-            }            
-
-            if (m_nla_core.is_monic_var(v)) {
-                m_mon_set.insert(v);
-                for (auto w : m_nla_core.emons()[v])
-                    todo.push_back(w);
-            }
-        }
-    }
+ 
 
     void reset() {
         m_values = nullptr;
         m_tmp1 = nullptr; m_tmp2 = nullptr;
         m_nlsat = alloc(nlsat::solver, m_limit, m_params, false);
         m_values = alloc(scoped_anum_vector, am());
-        m_term_set.reset();
         m_lp2nl.reset();
     }
 
@@ -149,17 +73,17 @@ struct solver::imp {
         reset();
         vector<nlsat::assumption, false> core;
 
-        init_cone_of_influence();
+        m_coi.init();
         // add linear inequalities from lra_solver
-        for (auto ci : m_constraint_set)
+        for (auto ci : m_coi.constraints())
             add_constraint(ci);
         
         // add polynomial definitions.
-        for (auto const& m : m_mon_set)
+        for (auto const& m : m_coi.mons())
             add_monic_eq(m_nla_core.emons()[m]);
 
         // add term definitions.
-        for (unsigned i : m_term_set)
+        for (unsigned i : m_coi.terms())
             add_term(i);
 
         TRACE(nra, m_nlsat->display(tout));
@@ -370,7 +294,7 @@ struct solver::imp {
         for (auto const& m : m_nla_core.emons())
             if (any_of(m.vars(), [&](lp::lpvar v) { return m_lp2nl.contains(v); }))
                 add_monic_eq_bound(m);
-        for (unsigned i : m_term_set)
+        for (unsigned i : m_coi.terms())
             add_term(i);
         for (auto const& [v, w] : m_lp2nl) {
             if (lra.column_has_lower_bound(v))
@@ -554,8 +478,8 @@ struct solver::imp {
         if (!m_lp2nl.find(v, r)) {
             r = m_nlsat->mk_var(is_int(v));
             m_lp2nl.insert(v, r);
-            if (!m_term_set.contains(v) && lra.column_has_term(v)) {
-                m_term_set.insert(v);
+            if (!m_coi.terms().contains(v) && lra.column_has_term(v)) {
+                m_coi.terms().insert(v);
             }
         }
         return r;
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 2155ee015..fd13d1575 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -2022,6 +2022,7 @@ public:
     }    
 
     void false_case_of_check_nla(const nla::lemma & l) {
+        TRACE(arith, tout << "nla false case\n";);
         m_lemma = l; //todo avoid the copy
         m_explanation = l.expl();
         literal_vector core;
@@ -3472,13 +3473,19 @@ public:
         // lp().shrink_explanation_to_minimum(m_explanation); // todo, enable when perf is fixed
         ++m_num_conflicts;
         ++m_stats.m_conflicts;
+        for (auto ev : m_explanation)
+            set_evidence(ev.ci(), m_core, m_eqs);
+        if (m_eqs.empty() && all_of(m_core, [&](literal l) { return ctx().get_assignment(l) == l_false; }))
+            is_conflict = true;
+        for (auto l : m_core) {
+            verbose_stream() << l << " " << ctx().get_assignment(l) << "\n";
+        }
         TRACE(arith_conflict,
               tout << "@" << ctx().get_scope_level() << (is_conflict ? " conflict":" lemma");
               for (auto const& p : m_params) tout << " " << p;
               tout << "\n";
               display_evidence(tout << core << " ", m_explanation););
-        for (auto ev : m_explanation) 
-            set_evidence(ev.ci(), m_core, m_eqs);
+
 
 
         if (params().m_arith_validate)