diff --git a/src/math/interval/dep_intervals.cpp b/src/math/interval/dep_intervals.cpp
index 860647883..56f3e1a0c 100644
--- a/src/math/interval/dep_intervals.cpp
+++ b/src/math/interval/dep_intervals.cpp
@@ -89,21 +89,23 @@ bool dep_intervals::separated_from_zero_on_upper(const interval& i) const {
 std::ostream& dep_intervals::display(std::ostream& out, const interval& i) const {
     if (m_imanager.lower_is_inf(i)) {
         out << "(-oo";
-    } else {
+    } 
+    else {
         out << (m_imanager.lower_is_open(i)? "(":"[") << rational(m_imanager.lower(i));        
     }
     out << ",";
     if (m_imanager.upper_is_inf(i)) {
         out << "oo)";
-    } else {
+    } 
+    else {
         out << rational(m_imanager.upper(i)) << (m_imanager.upper_is_open(i)? ")":"]");         
     }
     if (i.m_lower_dep) {
-        out << "\nlower deps\n";
+        // out << "\nlower deps\n";
         // TBD: print_dependencies(i.m_lower_dep, out);
     }
     if (i.m_upper_dep) {
-        out << "\nupper deps\n";
+        // out << "\nupper deps\n";
         // TBD: print_dependencies(i.m_upper_dep, out);   
     }
     return out;
@@ -125,21 +127,21 @@ bool dep_intervals::is_empty(interval const& a) const {
 bool dep_intervals::is_above(const interval& i, const rational& r) const {
     if (lower_is_inf(i))
         return false;
-    if (m_num_manager.lt(lower(i), r.to_mpq()))
-        return false;
-    if (m_num_manager.eq(lower(i), r.to_mpq()) && !m_config.lower_is_open(i))
-        return false;
-    return true;
+    if (m_num_manager.lt(r.to_mpq(), lower(i)))
+        return true;
+    if (m_num_manager.eq(lower(i), r.to_mpq()) && m_config.lower_is_open(i))
+        return true;
+    return false;
 }
 
 bool dep_intervals::is_below(const interval& i, const rational& r) const {
     if (upper_is_inf(i))
         return false;
     if (m_num_manager.lt(upper(i), r.to_mpq()))
-        return false;
-    if (m_num_manager.eq(upper(i), r.to_mpq()) && !m_config.upper_is_open(i))
-        return false;
-    return true;
+        return true;
+    if (m_num_manager.eq(upper(i), r.to_mpq()) && m_config.upper_is_open(i))
+        return true;
+    return false;
 }
 
 
diff --git a/src/math/interval/dep_intervals.h b/src/math/interval/dep_intervals.h
index b0b752e36..e316c8095 100644
--- a/src/math/interval/dep_intervals.h
+++ b/src/math/interval/dep_intervals.h
@@ -411,12 +411,10 @@ private:
     void add(const interval& a, const interval& b, interval& c, interval_deps_combine_rule& deps) { m_imanager.add(a, b, c, deps); }
    
     void combine_deps(interval const& a, interval const& b, interval_deps_combine_rule const& deps, interval& i) const {
-        SASSERT(&a != &i && &b != &i);
         m_config.add_deps(a, b, deps, i);
     }
 
     void combine_deps(interval const& a, interval_deps_combine_rule const& deps, interval& i) const {
-        SASSERT(&a != &i);
         m_config.add_deps(a, deps, i);
     }
 
diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp
index 09bed9554..69cdb3364 100644
--- a/src/math/lp/monomial_bounds.cpp
+++ b/src/math/lp/monomial_bounds.cpp
@@ -13,7 +13,9 @@
 
 namespace nla {
 
-    monomial_bounds::monomial_bounds(core* c):common(c) {}
+    monomial_bounds::monomial_bounds(core* c):
+        common(c), 
+        dep(c->m_intervals.get_dep_intervals()) {}
 
     bool monomial_bounds::operator()() {
         bool propagated = false;
@@ -29,9 +31,7 @@ namespace nla {
      * Accumulate product of variables in monomial starting at position 'start'
      */
     void monomial_bounds::compute_product(unsigned start, monic const& m, scoped_dep_interval& product) {
-        auto & intervals = c().m_intervals;
-        auto & dep_intervals = intervals.get_dep_intervals();
-        scoped_dep_interval vi(dep_intervals);
+        scoped_dep_interval vi(dep);
         for (unsigned i = start; i < m.size(); ) {
             lpvar v = m.vars()[i];
             unsigned power = 1;
@@ -40,8 +40,8 @@ namespace nla {
             for (; i < m.size() && m.vars()[i] == v; ++i) {
                 ++power;
             }
-            dep_intervals.power<dep_intervals::with_deps>(vi, power, vi);            
-            dep_intervals.mul<dep_intervals::with_deps>(product, vi, product);
+            dep.power<dep_intervals::with_deps>(vi, power, vi);            
+            dep.mul<dep_intervals::with_deps>(product, vi, product);
         }
     }
 
@@ -51,26 +51,27 @@ namespace nla {
      * a bounds axiom.
      */
     bool monomial_bounds::propagate_value(dep_interval& range, lpvar v) {
-        auto & intervals = c().m_intervals;
-        auto & dep_intervals = intervals.get_dep_intervals();
-        if (dep_intervals.is_below(range, c().val(v))) {
+        auto val = c().val(v);
+        if (dep.is_below(range, val)) {
             lp::explanation ex;
-            dep_intervals.get_upper_dep(range, ex);
-            auto const& upper = dep_intervals.upper(range);
-            auto cmp = dep_intervals.upper_is_open(range) ? llc::LT : llc::LE;
+            dep.get_upper_dep(range, ex);
+            auto const& upper = dep.upper(range);
+            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";);
             return true;
         }
-        else if (dep_intervals.is_above(range, c().val(v))) {
+        else if (dep.is_above(range, val)) {
             lp::explanation ex;
-            dep_intervals.get_lower_dep(range, ex);
-            auto const& lower = dep_intervals.lower(range);
-            auto cmp = dep_intervals.lower_is_open(range) ? llc::GT : llc::GE;
+            dep.get_lower_dep(range, ex);
+            auto const& lower = dep.lower(range);
+            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";);
             return true;
         }
         else {
@@ -79,26 +80,26 @@ namespace nla {
     }
 
     void monomial_bounds::var2interval(lpvar v, scoped_dep_interval& i) {
-        auto & intervals = c().m_intervals;
-        auto & dep_intervals = intervals.get_dep_intervals();
         lp::constraint_index ci;
         rational bound;
         bool is_strict;
         if (c().has_lower_bound(v, ci, bound, is_strict)) {
-            dep_intervals.set_lower_is_open(i, is_strict);
-            dep_intervals.set_lower(i, bound);
-            dep_intervals.set_lower_dep(i, dep_intervals.mk_leaf(ci));
+            dep.set_lower_is_open(i, is_strict);
+            dep.set_lower(i, bound);
+            dep.set_lower_dep(i, dep.mk_leaf(ci));
+            dep.set_lower_is_inf(i, false);
         }
         else {
-            dep_intervals.set_lower_is_inf(i, true);
+            dep.set_lower_is_inf(i, true);
         }
         if (c().has_upper_bound(v, ci, bound, is_strict)) {
-            dep_intervals.set_upper_is_open(i, is_strict);
-            dep_intervals.set_upper(i, bound);
-            dep_intervals.set_upper_dep(i, dep_intervals.mk_leaf(ci));            
+            dep.set_upper_is_open(i, is_strict);
+            dep.set_upper(i, bound);
+            dep.set_upper_dep(i, dep.mk_leaf(ci));            
+            dep.set_upper_is_inf(i, false);
         }
         else {
-            dep_intervals.set_upper_is_inf(i, true);
+            dep.set_upper_is_inf(i, true);
         }
     }
 
@@ -110,15 +111,22 @@ namespace nla {
      * If the value of m.var() is outside of product_of_all_vars, add a bounds lemma.
      */
     bool monomial_bounds::propagate(monic const& m) {
-        auto & intervals = c().m_intervals;
-        auto & dep_intervals = intervals.get_dep_intervals();
-        scoped_dep_interval product(dep_intervals);
-        scoped_dep_interval vi(dep_intervals), mi(dep_intervals);
-        scoped_dep_interval other_product(dep_intervals);
-        var2interval(m.var(), mi);
-        if (dep_intervals.lower_is_inf(mi) && dep_intervals.upper_is_inf(mi))
+        unsigned num_free, power;
+        lpvar free_var;
+        analyze_monomial(m, num_free, free_var, power);
+        bool m_is_free = is_free(m.var());
+        if (num_free >= 2)
             return false;
-        dep_intervals.set_value(product, rational::one());
+        if (num_free >= 1 && m_is_free)
+            return false;
+        SASSERT(num_free == 0 || !m_is_free);
+        bool do_propagate_up   = num_free == 0;
+        bool do_propagate_down = !m_is_free;
+        scoped_dep_interval product(dep);
+        scoped_dep_interval vi(dep), mi(dep);
+        scoped_dep_interval other_product(dep);
+        var2interval(m.var(), mi);
+        dep.set_value(product, rational::one());
         for (unsigned i = 0; i < m.size(); ) {
             lpvar v = m.vars()[i];
             ++i;
@@ -126,28 +134,60 @@ namespace nla {
             for (; i < m.size() && v == m.vars()[i]; ++i) 
                 ++power;
             var2interval(v, vi);
-            dep_intervals.power<dep_intervals::with_deps>(vi, power, vi);            
-            if (power == 1) {
-                dep_intervals.set<dep_intervals::with_deps>(other_product, product);
+            dep.power<dep_intervals::with_deps>(vi, power, vi);            
+
+            if (power == 1 && do_propagate_down && (num_free == 0 || free_var == v)) {
+                dep.set<dep_intervals::with_deps>(other_product, product);
                 compute_product(i, m, other_product);
                 if (propagate_down(m, mi, v, other_product))
                     return true;
             }
-            dep_intervals.mul<dep_intervals::with_deps>(product, vi, product);
+            dep.mul<dep_intervals::with_deps>(product, vi, product);
         }
-        return propagate_value(product, m.var());
+        return do_propagate_up && propagate_value(product, m.var());
     }
 
     bool monomial_bounds::propagate_down(monic const& m, dep_interval& mi, lpvar v, dep_interval& product) {
-        auto & intervals = c().m_intervals;
-        auto & dep_intervals = intervals.get_dep_intervals();
-        if (!dep_intervals.separated_from_zero(product)) 
+        if (!dep.separated_from_zero(product)) 
             return false;
-        scoped_dep_interval range(dep_intervals);
-        dep_intervals.set<dep_intervals::with_deps>(range, mi);
-        dep_intervals.div<dep_intervals::with_deps>(range, product, range);
+        scoped_dep_interval range(dep);
+        dep.div<dep_intervals::with_deps>(mi, product, range);
         return propagate_value(range, v);
     }
 
+    bool monomial_bounds::is_free(lpvar v) const {
+        return !c().has_lower_bound(v) && !c().has_upper_bound(v);
+    }    
+
+    bool monomial_bounds::is_zero(lpvar v) const {
+        return 
+            c().has_lower_bound(v) && 
+            c().has_upper_bound(v) &&
+            c().get_lower_bound(v).is_zero() && 
+            c().get_lower_bound(v) == c().get_upper_bound(v);
+    }    
+
+    void monomial_bounds::analyze_monomial(monic const& m, unsigned& num_free, lpvar& fv, unsigned& fv_power) const {
+        unsigned power = 0;
+        num_free = 0;
+        fv = null_lpvar;
+        fv_power = 0;
+        for (unsigned i = 0; i < m.vars().size(); ) {
+            lpvar v = m.vars()[i];
+            unsigned power = 1;
+            ++i;
+            for (; i < m.vars().size() && m.vars()[i] == v; ++i, ++power);
+            if (is_zero(v)) {
+                num_free = 0;
+                return;
+            }
+            if (power % 2 == 1 && is_free(v)) {
+                ++num_free;
+                fv_power = power;
+                fv = v;
+            }
+        }
+    }
+
 }
 
diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h
index 714bcbd51..aaf39d175 100644
--- a/src/math/lp/monomial_bounds.h
+++ b/src/math/lp/monomial_bounds.h
@@ -10,20 +10,21 @@
 
 #include "math/lp/nla_common.h"
 #include "math/lp/nla_intervals.h"
-#include "math/lp/nex.h"
-#include "math/lp/cross_nested.h"
 #include "math/lp/u_set.h"
 
 namespace nla {
     class core;
-
     class monomial_bounds : common {
+        dep_intervals& dep;
         void var2interval(lpvar v, scoped_dep_interval& i);
         bool propagate_down(monic const& m, lpvar u);
         bool propagate_value(dep_interval& range, lpvar v);
         void compute_product(unsigned start, monic const& m, scoped_dep_interval& i);
         bool propagate(monic const& m);
         bool propagate_down(monic const& m, dep_interval& mi, lpvar v, dep_interval& product);
+        void analyze_monomial(monic const& m, unsigned& num_free, lpvar& free_v, unsigned& power) const;
+        bool is_free(lpvar v) const;
+        bool is_zero(lpvar v) const;
     public:
         monomial_bounds(core* core);
         bool operator()();
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index 08f3841d9..d83bc89c5 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -29,6 +29,7 @@ core::core(lp::lar_solver& s, reslimit & lim) :
     m_monotone(this),
     m_intervals(this, lim),
     m_horner(this),
+    m_monomial_bounds(this),
     m_pdd_manager(s.number_of_vars()),
     m_pdd_grobner(lim, m_pdd_manager),
     m_emons(m_evars),
@@ -137,13 +138,13 @@ void core::add_monic(lpvar v, unsigned sz, lpvar const* vs) {
 }
     
 void core::push() {
-    TRACE("nla_solver",);
+    TRACE("nla_solver_verbose", tout << "\n";);
     m_emons.push();
 }
 
      
 void core::pop(unsigned n) {
-    TRACE("nla_solver", tout << "n = " << n << "\n";);
+    TRACE("nla_solver_verbose", tout << "n = " << n << "\n";);
     m_emons.pop(n);
     SASSERT(elists_are_consistent(false));
 }
@@ -403,10 +404,8 @@ bool core::explain_by_equiv(const lp::lar_term& t, lp::explanation& e) const {
             
     m_evars.explain(signed_var(i, false), signed_var(j, sign), e);
     TRACE("nla_solver", tout << "explained :"; m_lar_solver.print_term_as_indices(t, tout););
-    return true;
-            
+    return true;            
 }
-    
 
 void core::mk_ineq_no_expl_check(new_lemma& lemma, lp::lar_term& t, llc cmp, const rational& rs) {
     TRACE("nla_solver_details", m_lar_solver.print_term_as_indices(t, tout << "t = "););
@@ -424,8 +423,7 @@ llc apply_minus(llc cmp) {
     default: break;
     }
     return cmp;
-}
-    
+}   
     
 // the monics should be equal by modulo sign but this is not so in the model
 void core::fill_explanation_and_lemma_sign(new_lemma& lemma, const monic& a, const monic & b, rational const& sign) {
@@ -475,9 +473,7 @@ int core::vars_sign(const svector<lpvar>& v) {
     }
     return sign;
 }
-
-    
-
+   
 bool core::has_upper_bound(lpvar j) const {
     return m_lar_solver.column_has_upper_bound(j);
 } 
@@ -1214,10 +1210,8 @@ std::ostream& new_lemma::display(std::ostream & out) const {
     for (lpvar j : c.collect_vars(lemma)) {
         c.print_var(j, out);
     }
-
     return out;
 }
-
     
 void core::negate_relation(new_lemma& lemma, unsigned j, const rational& a) {
     SASSERT(val(j) != a);
@@ -1238,21 +1232,6 @@ bool core::done() const {
         lp_settings().get_cancel_flag();
 }
 
-void core::incremental_linearization(bool constraint_derived) {
-    TRACE("nla_solver_details", print_terms(tout); tout << m_lar_solver.constraints(););
-    m_basics.basic_lemma(constraint_derived);
-    if (!m_lemma_vec->empty() || constraint_derived || done())
-        return;
-    TRACE("nla_solver", tout << "passed constraint_derived and basic lemmas\n";);
-    SASSERT(elists_are_consistent(true));
-    if (!done())
-        m_order.order_lemma();
-    if (!done()) 
-        m_monotone.monotonicity_lemma();
-    if (!done())
-        m_tangents.tangent_lemma();    
-}
-
 bool core::elist_is_consistent(const std::unordered_set<lpvar> & list) const {
     bool first = true;
     bool p;
@@ -1359,12 +1338,9 @@ bool core::patch_blocker(lpvar u, const monic& m) const {
 }
 
 bool core::try_to_patch(lpvar k, const rational& v, const monic & m) {
-    return m_lar_solver.try_to_patch(k, v,
-                                     [this, k, m](lpvar u) {
-                                         if (u == k)
-                                             return false; // ok to patch
-                                         return patch_blocker(u, m); },
-                                     [this](lpvar u) { update_to_refine_of_var(u); });
+    auto blocker = [this, k, m](lpvar u) { return u != k && patch_blocker(u, m); };
+    auto change_report = [this](lpvar u) { update_to_refine_of_var(u); };
+    return m_lar_solver.try_to_patch(k, v, blocker,  change_report);
 }
 
 bool in_power(const svector<lpvar>& vs, unsigned l) {
@@ -1465,22 +1441,25 @@ lbool core::check(vector<lemma>& l_vec) {
     patch_monomials_with_real_vars();
     if (m_to_refine.is_empty()) { return l_true; }   
     init_search();
-  
+
+    lbool ret = l_undef;
+
     set_use_nra_model(false);    
 
-    if (need_to_call_algebraic_methods()) {
-        if (!m_horner.horner_lemmas() && m_nla_settings.run_grobner() && !done()) {
-            clear_and_resize_active_var_set(); 
-            find_nl_cluster();
-            run_grobner();                
-        }
-    }
-    TRACE("nla_solver_details", print_terms(tout); tout << m_lar_solver.constraints(););
-    if (!done()) 
-        m_basics.basic_lemma(true);    
+    if (false && l_vec.empty() && !done()) 
+        m_monomial_bounds();
+    
+    if (l_vec.empty() && !done () && need_to_call_algebraic_methods()) 
+        m_horner.horner_lemmas();
 
-    TRACE("nla_solver", tout << "passed constraint_derived and basic lemmas\n";);
-    SASSERT(!l_vec.empty() || elists_are_consistent(true));
+    if (l_vec.empty() && !done() && m_nla_settings.run_grobner()) {
+        clear_and_resize_active_var_set(); 
+        find_nl_cluster();
+        run_grobner();                
+    }
+
+    if (l_vec.empty() && !done()) 
+        m_basics.basic_lemma(true);    
 
     if (l_vec.empty() && !done()) 
         m_basics.basic_lemma(false);
@@ -1495,16 +1474,13 @@ lbool core::check(vector<lemma>& l_vec) {
             m_tangents.tangent_lemma();
     }
 
-    if (!m_reslim.inc())
-        return l_undef;
-  
-    lbool ret = l_vec.empty() ? l_undef : l_false;
-
 #if 0
-    if (l_vec.empty()) 
+    if (l_vec.empty() && !done()) 
         ret = m_nra.check();
-    }        
 #endif
+
+    if (ret == l_undef && !l_vec.empty() && m_reslim.inc()) 
+        ret = l_false;
     
     TRACE("nla_solver", tout << "ret = " << ret << ", lemmas count = " << l_vec.size() << "\n";);
     IF_VERBOSE(2, if(ret == l_undef) {verbose_stream() << "Monomials\n"; print_monics(verbose_stream());});
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index 669f26114..fc303d8b1 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -22,6 +22,7 @@
 #include "math/lp/nla_settings.h"
 #include "math/lp/nex.h"
 #include "math/lp/horner.h"
+#include "math/lp/monomial_bounds.h"
 #include "math/lp/nla_intervals.h"
 #include "math/grobner/pdd_solver.h"
 #include "nlsat/nlsat_solver.h"
@@ -148,7 +149,8 @@ public:
     basics                   m_basics;
     order                    m_order;
     monotone                 m_monotone;
-    intervals                m_intervals;                
+    intervals                m_intervals; 
+    monomial_bounds          m_monomial_bounds;
     horner                   m_horner;
     nla_settings             m_nla_settings;    
     dd::pdd_manager          m_pdd_manager;