diff --git a/src/ast/arith_decl_plugin.h b/src/ast/arith_decl_plugin.h
index 4b24ea5e6..4850b52a8 100644
--- a/src/ast/arith_decl_plugin.h
+++ b/src/ast/arith_decl_plugin.h
@@ -370,7 +370,7 @@ public:
     app * mk_lt(expr * arg1, expr * arg2) const { return m_manager.mk_app(m_afid, OP_LT, arg1, arg2); }
     app * mk_gt(expr * arg1, expr * arg2) const { return m_manager.mk_app(m_afid, OP_GT, arg1, arg2); }
 
-    app * mk_add(unsigned num_args, expr * const * args) const { return m_manager.mk_app(m_afid, OP_ADD, num_args, args); }
+    app * mk_add(unsigned num_args, expr * const * args) const { return num_args == 1 && is_app(args[0]) ? to_app(args[0]) : m_manager.mk_app(m_afid, OP_ADD, num_args, args); }
     app * mk_add(expr * arg1, expr * arg2) const { return m_manager.mk_app(m_afid, OP_ADD, arg1, arg2); }
     app * mk_add(expr * arg1, expr * arg2, expr* arg3) const { return m_manager.mk_app(m_afid, OP_ADD, arg1, arg2, arg3); }
 
@@ -378,7 +378,7 @@ public:
     app * mk_sub(unsigned num_args, expr * const * args) const { return m_manager.mk_app(m_afid, OP_SUB, num_args, args); }
     app * mk_mul(expr * arg1, expr * arg2) const { return m_manager.mk_app(m_afid, OP_MUL, arg1, arg2); }
     app * mk_mul(expr * arg1, expr * arg2, expr* arg3) const { return m_manager.mk_app(m_afid, OP_MUL, arg1, arg2, arg3); }
-    app * mk_mul(unsigned num_args, expr * const * args) const { return m_manager.mk_app(m_afid, OP_MUL, num_args, args); }
+    app * mk_mul(unsigned num_args, expr * const * args) const { return num_args == 1 && is_app(args[0]) ? to_app(args[0]) : m_manager.mk_app(m_afid, OP_MUL, num_args, args); }
     app * mk_uminus(expr * arg) const { return m_manager.mk_app(m_afid, OP_UMINUS, arg); }
     app * mk_div(expr * arg1, expr * arg2) { return m_manager.mk_app(m_afid, OP_DIV, arg1, arg2); }
     app * mk_idiv(expr * arg1, expr * arg2) { return m_manager.mk_app(m_afid, OP_IDIV, arg1, arg2); }
diff --git a/src/math/polynomial/polynomial.cpp b/src/math/polynomial/polynomial.cpp
index d6d392148..51d8489e3 100644
--- a/src/math/polynomial/polynomial.cpp
+++ b/src/math/polynomial/polynomial.cpp
@@ -150,7 +150,6 @@ namespace polynomial {
         return r;
     }
 
-
     /**
        \brief Monomials (power products)
     */
@@ -192,7 +191,7 @@ namespace polynomial {
         };
 
         static unsigned get_obj_size(unsigned sz) { return sizeof(monomial) + sz * sizeof(power); }
-
+        
         monomial(unsigned id, unsigned sz, power const * pws, unsigned h):
             m_ref_count(0),
             m_id(id),
@@ -257,9 +256,7 @@ namespace polynomial {
             if (m_size < SMALL_MONOMIAL) {
                 // use linear search for small monomials
                 // search backwards since we usually ask for the degree of "big" variables
-                unsigned i = last;
-                while (i > 0) {
-                    --i;
+                for (unsigned i = last; i-- > 0; ) {
                     if (get_var(i) == x)
                         return i;
                 }
@@ -798,9 +795,8 @@ namespace polynomial {
             dec_ref(m_unit);
             CTRACE("polynomial", !m_monomials.empty(),
                    tout << "monomials leaked\n";
-                   monomial_table::iterator it = m_monomials.begin(); monomial_table::iterator end = m_monomials.end();
-                   for (; it != end; ++it) {
-                       (*it)->display(tout); tout << "\n";
+                   for (auto * m : m_monomials) {
+                       m->display(tout); tout << "\n";
                    });
             SASSERT(m_monomials.empty());
             if (m_own_allocator)
@@ -1510,6 +1506,8 @@ namespace polynomial {
         unsigned id() const { return m_id; }
         unsigned size() const { return m_size; }
         monomial * m(unsigned idx) const { SASSERT(idx < size()); return m_ms[idx]; }
+        monomial *const* begin() const { return m_ms; }
+        monomial *const* end() const { return m_ms + size(); }
         numeral const & a(unsigned idx) const { SASSERT(idx < size()); return m_as[idx]; }
         numeral & a(unsigned idx) { SASSERT(idx < size()); return m_as[idx]; }
         numeral const * as() const { return m_as; }
@@ -1773,11 +1771,9 @@ namespace polynomial {
     }
 
     bool manager::is_linear(polynomial const * p) {
-        unsigned sz = p->size();
-        for (unsigned i = 0; i < sz; i++) {
-            if (!is_linear(p->m(0)))
+        for (monomial* m : *p) 
+            if (!is_linear(m)) 
                 return false;
-        }
         return true;
     }
 
@@ -2396,6 +2392,7 @@ namespace polynomial {
             return mm().is_valid(x);
         }
 
+
         void add_del_eh(del_eh * eh) {
             eh->m_next = m_del_eh;
             m_del_eh = eh;
@@ -6101,6 +6098,33 @@ namespace polynomial {
                   });
         }
 
+        lbool sign(monomial* m, numeral const& c, svector<lbool> const& sign_of_vars) {
+            unsigned sz = size(m);
+            lbool sign1 = m_manager.is_pos(c) ? l_true : l_false;
+            for (unsigned i = 0; i < sz; ++i) {
+                var v = get_var(m, i);
+                unsigned d = degree(m, i);
+                lbool sign2 = sign_of_vars.get(v, l_undef);
+                if (sign2 == l_undef) 
+                    return l_undef;
+                else if (1 == (d % 2) && sign2 == l_false) {
+                    sign1 = sign1 == l_true ? l_false : l_true;
+                }
+            }
+            return sign1;
+        }
+
+        lbool sign(polynomial const * p, svector<lbool> const& sign_of_vars) {
+            unsigned sz = size(p);
+            if (sz == 0) return l_undef;
+            lbool sign1 = sign(p->m(0), p->a(0), sign_of_vars);
+            for (unsigned i = 1; sign1 != l_undef && i < sz; ++i) {
+                if (sign(p->m(i), p->a(i), sign_of_vars) != sign1) 
+                    return l_undef;
+            }
+            return sign1;
+        }
+
         bool is_pos(polynomial const * p) {
             bool found_unit = false;
             unsigned sz = p->size();
@@ -6372,6 +6396,31 @@ namespace polynomial {
                 R.add(new_a, mk_monomial(new_m));
             }
             return R.mk();
+        }        
+
+        void substitute(polynomial const* r, var x, polynomial const* p, polynomial const* q, polynomial_ref& result) {
+            unsigned md = degree(r, x);
+            if (md == 0) {
+                result = const_cast<polynomial*>(r);
+                return;
+            }
+            result = 0;
+            polynomial_ref p1(pm()), q1(pm());
+            polynomial_ref_buffer ps(pm());
+            unsigned sz = r->size();
+            for (unsigned i = 0; i < sz; i++) {
+                monomial * m0 = r->m(i);
+                unsigned dm = m0->degree_of(x);
+                SASSERT(md >= dm);
+                monomial_ref m1(div_x(m0, x), pm());
+                pw(p, dm, p1);
+                pw(q, md - dm, q1);
+                p1 = mul(r->a(i), m1, p1 * q1);
+                if (result) 
+                    result = add(result, p1);
+                else 
+                    result = p1;
+            }
         }
 
 
@@ -6918,6 +6967,18 @@ namespace polynomial {
         return m_imp->m().set_zp(p);
     }
 
+    bool manager::is_var(polynomial const* p, var& v) {
+        return p->size() == 1 && is_var(p->m(0), v) && m_imp->m().is_one(p->a(0));
+    }
+
+    bool manager::is_var(monomial const* m, var& v) {
+        return m->size() == 1 && m->degree(0) == 1 && (v = m->get_var(0), true);
+    }
+
+    bool manager::is_var_num(polynomial const* p, var& v, scoped_numeral& n) {
+        return p->size() == 2 && m_imp->m().is_one(p->a(0)) && is_var(p->m(0), v) && is_unit(p->m(1)) && (n = p->a(1), true);
+    }
+
     small_object_allocator & manager::allocator() const {
         return m_imp->mm().allocator();
     }
@@ -7271,6 +7332,10 @@ namespace polynomial {
     void manager::psc_chain(polynomial const * p, polynomial const * q, var x, polynomial_ref_vector & S) {
         m_imp->psc_chain(p, q, x, S);
     }
+    
+    lbool manager::sign(polynomial const * p, svector<lbool> const& sign_of_vars) {
+        return m_imp->sign(p, sign_of_vars);
+    }
 
     bool manager::is_pos(polynomial const * p) {
         return m_imp->is_pos(p);
@@ -7307,6 +7372,10 @@ namespace polynomial {
     polynomial * manager::substitute(polynomial const * p, unsigned xs_sz, var const * xs, numeral const * vs) {
         return m_imp->substitute(p, xs_sz, xs, vs);
     }
+    
+    void manager::substitute(polynomial const* r, var x, polynomial const* p, polynomial const* q, polynomial_ref& result) {
+        m_imp->substitute(r, x, p, q, result);
+    }
 
     void manager::factor(polynomial const * p, factors & r, factor_params const & params) {
         m_imp->factor(p, r, params);
diff --git a/src/math/polynomial/polynomial.h b/src/math/polynomial/polynomial.h
index 43b3c1138..b09af94bc 100644
--- a/src/math/polynomial/polynomial.h
+++ b/src/math/polynomial/polynomial.h
@@ -29,6 +29,7 @@ Notes:
 #include "util/params.h"
 #include "util/mpbqi.h"
 #include "util/rlimit.h"
+#include "util/lbool.h"
 
 class small_object_allocator;
 
@@ -98,7 +99,7 @@ namespace polynomial {
     };
 
     struct display_var_proc {
-        virtual void operator()(std::ostream & out, var x) const { out << "x" << x; }
+        virtual std::ostream& operator()(std::ostream & out, var x) const { return out << "x" << x; }
     };
 
     class polynomial;
@@ -306,12 +307,27 @@ namespace polynomial {
            \brief Return true if m is linear (i.e., it is of the form 1 or x).
         */
         static bool is_linear(monomial const * m);
-        
+       
         /**
            \brief Return true if all monomials in p are linear.
         */
         static bool is_linear(polynomial const * p);
 
+        /**
+           \brief Return true if the monomial is a variable.
+        */
+        static bool is_var(monomial const* p, var& v);
+
+        /**
+           \brief Return true if the polynomial is a variable.
+        */
+        bool is_var(polynomial const* p, var& v);
+
+        /**
+           \brief Return true if the polynomial is of the form x + k
+        */
+        bool is_var_num(polynomial const* p, var& v, scoped_numeral& n);
+
         /**
            \brief Return the degree of variable x in p.
         */
@@ -860,7 +876,13 @@ namespace polynomial {
            \brief Return true if p is a square, and store its square root in r.
         */
         bool sqrt(polynomial const * p, polynomial_ref & r);
-        
+       
+
+        /**
+           \brief obtain the sign of the polynomial given sign of variables.
+        */
+        lbool sign(polynomial const* p, svector<lbool> const& sign_of_vars);
+ 
         /**
            \brief Return true if p is always positive for any assignment of its variables.
            
@@ -936,6 +958,13 @@ namespace polynomial {
             return substitute(p, 1, &x, &v);
         }
 
+        /**
+           \brief Apply substiution [x -> p/q] in r.
+           That is, given r \in Z[x, y_1, .., y_m] return
+           polynomial q^k * r(p/q, y_1, .., y_m), where k is the maximal degree of x in r.
+        */
+        void substitute(polynomial const* r, var x, polynomial const* p, polynomial const* q, polynomial_ref& result);
+
         /**
            \brief Factorize the given polynomial p and store its factors in r.
         */
diff --git a/src/nlsat/nlsat_clause.h b/src/nlsat/nlsat_clause.h
index 898c32449..93dfb0e80 100644
--- a/src/nlsat/nlsat_clause.h
+++ b/src/nlsat/nlsat_clause.h
@@ -44,6 +44,8 @@ namespace nlsat {
         bool is_learned() const { return m_learned; }
         literal * begin() { return m_lits; }
         literal * end() { return m_lits + m_size; }
+        literal const * begin() const { return m_lits; }
+        literal const * end() const { return m_lits + m_size; }
         literal const * c_ptr() const { return m_lits; }
         void inc_activity() { m_activity++; }
         void set_activity(unsigned v) { m_activity = v; }
diff --git a/src/nlsat/nlsat_solver.cpp b/src/nlsat/nlsat_solver.cpp
index e1dadf12b..4cba6d123 100644
--- a/src/nlsat/nlsat_solver.cpp
+++ b/src/nlsat/nlsat_solver.cpp
@@ -18,19 +18,20 @@ Author:
 Revision History:
 
 --*/
+#include "util/z3_exception.h"
+#include "util/chashtable.h"
+#include "util/id_gen.h"
+#include "util/map.h"
+#include "util/dependency.h"
+#include "util/permutation.h"
+#include "math/polynomial/algebraic_numbers.h"
+#include "math/polynomial/polynomial_cache.h"
 #include "nlsat/nlsat_solver.h"
 #include "nlsat/nlsat_clause.h"
 #include "nlsat/nlsat_assignment.h"
 #include "nlsat/nlsat_justification.h"
 #include "nlsat/nlsat_evaluator.h"
 #include "nlsat/nlsat_explain.h"
-#include "math/polynomial/algebraic_numbers.h"
-#include "util/z3_exception.h"
-#include "util/chashtable.h"
-#include "util/id_gen.h"
-#include "util/dependency.h"
-#include "math/polynomial/polynomial_cache.h"
-#include "util/permutation.h"
 #include "nlsat/nlsat_params.hpp"
 
 #define NLSAT_EXTRA_VERBOSE
@@ -68,6 +69,7 @@ namespace nlsat {
 
         reslimit&              m_rlimit;
         small_object_allocator m_allocator;
+        bool                   m_incremental;
         unsynch_mpq_manager    m_qm;
         pmanager               m_pm;
         cache                  m_cache;
@@ -78,6 +80,8 @@ namespace nlsat {
         interval_set_manager & m_ism;
         ineq_atom_table        m_ineq_atoms;
         root_atom_table        m_root_atoms;
+        svector<bool_var>      m_patch_var;
+        polynomial_ref_vector  m_patch_num, m_patch_denom;
         
         id_gen                 m_cid_gen;
         clause_vector          m_clauses; // set of clauses
@@ -108,11 +112,12 @@ namespace nlsat {
                 m_perm(perm),
                 m_proc(0) {
             }
-            virtual void operator()(std::ostream & out, var x) const { 
+            std::ostream& operator()(std::ostream & out, var x) const { 
                 if (m_proc == 0)
                     m_default_display_var(out, x);
                 else
                     (*m_proc)(out, m_perm[x]);
+                return out;
             }
         };
         perm_display_var_proc  m_display_var;
@@ -158,9 +163,10 @@ namespace nlsat {
         unsigned               m_stages;
         unsigned               m_irrational_assignments; // number of irrational witnesses
 
-        imp(solver& s, reslimit& rlim, params_ref const & p):
+        imp(solver& s, reslimit& rlim, params_ref const & p, bool incremental):
             m_rlimit(rlim),
             m_allocator("nlsat"),
+            m_incremental(incremental),
             m_pm(rlim, m_qm, &m_allocator),
             m_cache(m_pm),
             m_am(rlim, m_qm, p, &m_allocator),
@@ -168,6 +174,8 @@ namespace nlsat {
             m_assignment(m_am),
             m_evaluator(s, m_assignment, m_pm, m_allocator), 
             m_ism(m_evaluator.ism()),
+            m_patch_num(m_pm),
+            m_patch_denom(m_pm),
             m_num_bool_vars(0),
             m_display_var(m_perm),
             m_explain(s, m_assignment, m_cache, m_atoms, m_var2eq, m_evaluator),
@@ -330,9 +338,7 @@ namespace nlsat {
         */
         bool_var max_bvar(clause const & cls) const {
             bool_var b = null_bool_var;
-            unsigned sz = cls.size();
-            for (unsigned i = 0; i < sz; i++) {
-                literal l = cls[i];
+            for (literal l : cls) {
                 if (b == null_bool_var || l.var() > b)
                     b = l.var();
             }
@@ -367,9 +373,7 @@ namespace nlsat {
             if (x == null_var)
                 return 0;
             unsigned max = 0;
-            unsigned sz  = c.size();
-            for (unsigned i = 0; i < sz; i++) {
-                literal l = c[i];
+            for (literal l : c) {
                 atom const * a  = m_atoms[l.var()];
                 if (a == 0)
                     continue;
@@ -496,13 +500,12 @@ namespace nlsat {
         
         // Delete atoms with ref_count == 0
         void del_unref_atoms() {
-            unsigned sz = m_atoms.size();
-            for (unsigned i = 0; i < sz; i++) {
-                del(m_atoms[i]);
+            for (auto* a : m_atoms) {
+                del(a);
             }
         }
 
-        bool_var mk_ineq_atom(atom::kind k, unsigned sz, poly * const * ps, bool const * is_even) {
+        ineq_atom* mk_ineq_atom(atom::kind k, unsigned sz, poly * const * ps, bool const * is_even, bool& is_new) {
             SASSERT(sz >= 1);
             SASSERT(k == atom::LT || k == atom::GT || k == atom::EQ);
             int sign = 1;
@@ -522,23 +525,36 @@ namespace nlsat {
             void * mem = m_allocator.allocate(ineq_atom::get_obj_size(sz));
             if (sign < 0)
                 k = atom::flip(k);
-            ineq_atom * new_atom = new (mem) ineq_atom(k, sz, uniq_ps.c_ptr(), is_even, max);
+            ineq_atom * tmp_atom = new (mem) ineq_atom(k, sz, uniq_ps.c_ptr(), is_even, max);
             TRACE("nlsat_table_bug", ineq_atom::hash_proc h; 
-                  tout << "mk_ineq_atom hash: " << h(new_atom) << "\n"; display(tout, *new_atom, m_display_var); tout << "\n";);
-            ineq_atom * old_atom = m_ineq_atoms.insert_if_not_there(new_atom);
-            CTRACE("nlsat_table_bug", old_atom->max_var() != max, display(tout, *old_atom, m_display_var); tout << "\n";);
-            SASSERT(old_atom->max_var() == max);
-            if (old_atom != new_atom) {
-                deallocate(new_atom);
-                return old_atom->bvar();
+                  tout << "mk_ineq_atom hash: " << h(tmp_atom) << "\n"; display(tout, *tmp_atom, m_display_var); tout << "\n";);
+            ineq_atom * atom = m_ineq_atoms.insert_if_not_there(tmp_atom);
+            CTRACE("nlsat_table_bug", atom->max_var() != max, display(tout, *atom, m_display_var); tout << "\n";);
+            SASSERT(atom->max_var() == max);
+            is_new = (atom == tmp_atom);
+            if (is_new) {
+                for (unsigned i = 0; i < sz; i++) {
+                    m_pm.inc_ref(atom->p(i));
+                }
             }
-            bool_var b = mk_bool_var_core();
-            m_atoms[b] = new_atom;
-            new_atom->m_bool_var = b;
-            for (unsigned i = 0; i < sz; i++) {
-                m_pm.inc_ref(new_atom->p(i));
+            else {
+                deallocate(tmp_atom);
+            }
+            return atom;
+        }
+
+        bool_var mk_ineq_atom(atom::kind k, unsigned sz, poly * const * ps, bool const * is_even) {
+            bool is_new = false;
+            ineq_atom* atom = mk_ineq_atom(k, sz, ps, is_even, is_new);
+            if (!is_new) {
+                return atom->bvar();
+            }
+            else {
+                bool_var b = mk_bool_var_core();
+                m_atoms[b] = atom;
+                atom->m_bool_var = b;
+                return b;
             }
-            return b;
         }
 
         literal mk_ineq_literal(atom::kind k, unsigned sz, poly * const * ps, bool const * is_even) {
@@ -619,10 +635,14 @@ namespace nlsat {
             deallocate(cls);
         }
 
+        void del_clause(clause * cls, clause_vector& clauses) {
+            clauses.erase(cls);
+            del_clause(cls);
+        }
+
         void del_clauses(ptr_vector<clause> & cs) {
-            unsigned sz = cs.size();
-            for (unsigned i = 0; i < sz; i++)
-                del_clause(cs[i]);
+            for (clause* cp : cs) 
+                del_clause(cp);
         }
 
         void del_clauses() {
@@ -927,10 +947,9 @@ namespace nlsat {
            \brief Return true if the given clause is already satisfied in the current partial interpretation.
         */
         bool is_satisfied(clause const & cls) const {
-            unsigned sz = cls.size();
-            for (unsigned i = 0; i < sz; i++) {
-                if (const_cast<imp*>(this)->value(cls[i]) == l_true) {
-                    TRACE("value_bug:", tout << cls[i] << " := true\n";);
+            for (literal l : cls) {
+                if (const_cast<imp*>(this)->value(l) == l_true) {
+                    TRACE("value_bug:", tout << l << " := true\n";);
                     return true;
                 }
             }
@@ -1023,7 +1042,7 @@ namespace nlsat {
             SASSERT(x != null_var);
             if (m_var2eq[x] != 0 && degree(m_var2eq[x]) <= degree(a))
                 return; // we only update m_var2eq if the new equality has smaller degree
-            TRACE("simplify_core", tout << "Saving equality for "; m_display_var(tout, x); tout << " (x" << x << ") ";
+            TRACE("simplify_core", tout << "Saving equality for "; m_display_var(tout, x) << " (x" << x << ") ";
                   tout << "scope-lvl: " << scope_lvl() << "\n"; display(tout, literal(b, false)); tout << "\n";);
             save_updt_eq_trail(m_var2eq[x]);
             m_var2eq[x] = a;
@@ -1045,10 +1064,9 @@ namespace nlsat {
             interval_set_ref first_undef_set(m_ism); // infeasible region of the first undefined literal
             interval_set * xk_set = m_infeasible[m_xk]; // current set of infeasible interval for current variable
             SASSERT(!m_ism.is_full(xk_set));
-            unsigned sz = cls.size();
-            for (unsigned i = 0; i < sz; i++) {
+            for (unsigned idx = 0; idx < cls.size(); ++idx) {
+                literal l = cls[idx];
                 checkpoint();
-                literal l = cls[i];
                 if (value(l) == l_false)
                     continue;
                 SASSERT(value(l) == l_undef);
@@ -1085,7 +1103,7 @@ namespace nlsat {
                 }
                 num_undef++;
                 if (first_undef == UINT_MAX) {
-                    first_undef = i;
+                    first_undef = idx;
                     first_undef_set = curr_set;
                 }
             }
@@ -1131,9 +1149,7 @@ namespace nlsat {
            Return 0, if the set was satisfied, or the violating clause otherwise
         */
         clause * process_clauses(clause_vector const & cs) {
-            unsigned sz = cs.size();
-            for (unsigned i = 0; i < sz; i++) {
-                clause * c = cs[i];
+            for (clause* c : cs) {
                 if (!process_clause(*c, false))
                     return c;
             }
@@ -1175,7 +1191,7 @@ namespace nlsat {
             m_ism.peek_in_complement(m_infeasible[m_xk], m_is_int[m_xk], w, m_randomize);
             TRACE("nlsat", 
                   tout << "infeasible intervals: "; m_ism.display(tout, m_infeasible[m_xk]); tout << "\n";
-                  tout << "assigning "; m_display_var(tout, m_xk); tout << "(x" << m_xk << ") -> " << w << "\n";);
+                  tout << "assigning "; m_display_var(tout, m_xk) << "(x" << m_xk << ") -> " << w << "\n";);
             TRACE("nlsat_root", tout << "value as root object: "; m_am.display_root(tout, w); tout << "\n";);
             if (!m_am.is_rational(w))
                 m_irrational_assignments++;
@@ -1205,6 +1221,7 @@ namespace nlsat {
                 TRACE("nlsat_bug", tout << "xk: x" << m_xk << " bk: b" << m_bk << "\n";);
                 if (m_bk == null_bool_var && m_xk >= num_vars()) {
                     TRACE("nlsat", tout << "found model\n"; display_assignment(tout););
+                    fix_patch();
                     return l_true; // all variables were assigned, and all clauses were satisfied.
                 }
                 TRACE("nlsat", tout << "processing variable "; 
@@ -1295,6 +1312,13 @@ namespace nlsat {
             init_search();
             m_explain.set_full_dimensional(is_full_dimensional());
             bool reordered = false;
+
+#if 0
+            // disabled
+            if (!m_incremental) 
+                simplify();
+#endif
+
             if (!can_reorder()) {
 
             }
@@ -1497,9 +1521,9 @@ namespace nlsat {
             TRACE("nlsat_mathematica", display_mathematica_lemma(tout, m_lazy_clause.size(), m_lazy_clause.c_ptr()););            
             TRACE("nlsat_proof_sk", tout << "theory lemma\n"; display_abst(tout, m_lazy_clause.size(), m_lazy_clause.c_ptr()); tout << "\n";); 
             TRACE("nlsat_resolve", 
-                  tout << "m_xk: " << m_xk << ", "; m_display_var(tout, m_xk); tout << "\n";
+                  tout << "m_xk: " << m_xk << ", "; m_display_var(tout, m_xk) << "\n";
                   tout << "new valid clause:\n";
-                  display(tout, m_lazy_clause.size(), m_lazy_clause.c_ptr()); tout << "\n";);
+                  display(tout, m_lazy_clause.size(), m_lazy_clause.c_ptr()) << "\n";);
 
             DEBUG_CODE({
                 unsigned sz = m_lazy_clause.size();
@@ -1633,15 +1657,10 @@ namespace nlsat {
             TRACE("nlsat_proof", tout << "STARTING RESOLUTION\n";);
             TRACE("nlsat_proof_sk", tout << "STARTING RESOLUTION\n";);
             m_conflicts++;
-            TRACE("nlsat", tout << "resolve, conflicting clause:\n"; display(tout, *conflict_clause); tout << "\n";
+            TRACE("nlsat", tout << "resolve, conflicting clause:\n"; display(tout, *conflict_clause) << "\n";
                   tout << "xk: "; if (m_xk != null_var) m_display_var(tout, m_xk); else tout << "<null>"; tout << "\n";
                   tout << "scope_lvl: " << scope_lvl() << "\n";
                   tout << "current assignment\n"; display_assignment(tout););
-
-            // static unsigned counter = 0;
-            // counter++;
-            // if (counter > 6)
-            //    exit(0);
             
             m_num_marks = 0;
             m_lemma.reset();
@@ -1750,7 +1769,7 @@ namespace nlsat {
                 SASSERT(m_xk == new_max_var);
                 new_cls = mk_clause(sz, m_lemma.c_ptr(), true, m_lemma_assumptions.get());
                 TRACE("nlsat", tout << "new_level: " << scope_lvl() << "\nnew_stage: " << new_max_var << " "; 
-                      if (new_max_var != null_var) m_display_var(tout, new_max_var); tout << "\n";);
+                      if (new_max_var != null_var) m_display_var(tout, new_max_var) << "\n";);
             }
             else {
                 SASSERT(scope_lvl() >= 1);
@@ -1931,11 +1950,12 @@ namespace nlsat {
                     collect(*(cs[i]));
             }
 
-            void display(std::ostream & out, display_var_proc const & proc) {
+            std::ostream& display(std::ostream & out, display_var_proc const & proc) {
                 unsigned sz = m_num_occs.size();
                 for (unsigned i = 0; i < sz; i++) {
                     proc(out, i); out << " -> " << m_max_degree[i] << " : " << m_num_occs[i] << "\n";
                 }
+                return out;
             }
         };
         
@@ -2169,12 +2189,10 @@ namespace nlsat {
         }
 
         void reattach_arith_clauses(clause_vector const & cs) {
-            unsigned sz = cs.size();
-            for (unsigned i = 0; i < sz; i++) {
-                clause & c = *cs[i];
-                var x      = max_var(c);
+            for (clause* cp : cs) {
+                var x = max_var(*cp);
                 if (x != null_var)
-                    m_watches[x].push_back(&c);
+                    m_watches[x].push_back(cp);
             }
         }
 
@@ -2251,18 +2269,16 @@ namespace nlsat {
         }
 
         bool is_full_dimensional(clause const & c) const {
-            unsigned sz = c.size();
-            for (unsigned i = 0; i < sz; i++) {
-                if (!is_full_dimensional(c[i]))
+            for (literal l : c) {
+                if (!is_full_dimensional(l))
                     return false;
             }
             return true;
         }
 
         bool is_full_dimensional(clause_vector const & cs) const {
-            unsigned sz = cs.size();
-            for (unsigned i = 0; i < sz; i++) {
-                if (!is_full_dimensional(*(cs[i])))
+            for (clause* c : cs) {
+                if (!is_full_dimensional(*c)) 
                     return false;
             }
             return true;
@@ -2272,13 +2288,271 @@ namespace nlsat {
             return is_full_dimensional(m_clauses);
         }
 
+
+        // -----------------------
+        //
+        // Simplification
+        //
+        // -----------------------
+
+        // solve simple equalities
+        // TBD WU-Reit decomposition?
+
+        /**
+           \brief isolate variables in unit equalities.
+           Assume a clause is c == v*p + q
+           and the context implies p > 0
+           
+           replace v by -q/p
+           remove clause c, 
+           The for other occurrences of v, 
+              replace v*r + v*v*r' > 0 by 
+              by p*p*v*r + p*p*v*v*r' > 0
+              by p*q*r + q*q*r' > 0        
+
+           The method ignores lemmas and assumes constraints don't use roots.
+        */
+
+        void simplify() {
+            polynomial_ref p(m_pm), q(m_pm);
+            var v;
+            init_var_signs();
+            SASSERT(m_learned.empty());
+            bool change = true;
+            while (change) {
+                change = false;
+                for (clause* c : m_clauses) {
+                    if (solve_var(*c, v, p, q)) {
+                        q = -q;
+                        m_patch_var.push_back(v);
+                        m_patch_num.push_back(q);
+                        m_patch_denom.push_back(p);
+                        del_clause(c, m_clauses);
+                        substitute_var(v, p, q);
+                        change = true;
+                        break;
+                    }
+                }
+            }
+        }
+
+        void fix_patch() {
+            for (unsigned i = m_patch_var.size(); i-- > 0; ) {
+                var v = m_patch_var[i];
+                poly* q = m_patch_num.get(i);
+                poly* p = m_patch_denom.get(i);
+                scoped_anum pv(m_am), qv(m_am), val(m_am);
+                m_pm.eval(p, m_assignment, pv);
+                m_pm.eval(q, m_assignment, qv);
+                val = qv / pv;
+                TRACE("nlsat", 
+                      m_display_var(tout << "patch ", v) << "\n";
+                      if (m_assignment.is_assigned(v)) m_am.display(tout << "previous value: ", m_assignment.value(v)); tout << "\n";
+                      m_am.display(tout << "updated value: ", val); tout << "\n";
+                      );
+                m_assignment.set_core(v, val);
+            }
+        }
+
+        void substitute_var(var x, poly* p, poly* q) {
+            polynomial_ref pr(m_pm);
+            polynomial_ref_vector ps(m_pm);
+            u_map<literal> b2l;
+            svector<bool> even;
+            unsigned num_atoms = m_atoms.size();
+            for (unsigned j = 0; j < num_atoms; ++j) {
+                atom* a = m_atoms[j];
+                if (a && a->is_ineq_atom()) {
+                    ineq_atom const& a1 = *to_ineq_atom(a);
+                    unsigned sz = a1.size();
+                    ps.reset();
+                    even.reset();
+                    bool change = false;
+                    for (unsigned i = 0; i < sz; ++i) {
+                        poly * po = a1.p(i);
+                        m_pm.substitute(po, x, q, p, pr);
+                        ps.push_back(pr);                                
+                        even.push_back(a1.is_even(i));
+                        change |= pr != po;
+                        if (m_pm.is_zero(pr)) {
+                            ps.reset();
+                            even.reset();
+                            change = true;
+                            break;
+                        }
+                    }               
+                    if (!change) continue;
+                    literal l = mk_ineq_literal(a1.get_kind(), ps.size(), ps.c_ptr(), even.c_ptr());                            
+                    if (a1.m_bool_var != l.var()) {                        
+                        b2l.insert(a1.m_bool_var, l);
+                        inc_ref(l);
+                    }
+                }
+            }
+            update_clauses(b2l);
+            for (auto const& kv : b2l) {
+                dec_ref(kv.m_value);
+            }
+        }
+
+
+        void update_clauses(u_map<literal> const& b2l) {
+            literal_vector lits;
+            clause_vector to_delete;
+            unsigned n = m_clauses.size();
+            for (unsigned i = 0; i < n; ++i) {
+                clause* c = m_clauses[i];
+                lits.reset();
+                bool changed = false;
+                bool is_tautology = false;
+                for (literal l : *c) {               
+                    literal lit = null_literal;
+                    if (b2l.find(l.var(), lit)) {
+                        lit = l.sign() ? ~lit : lit;
+                        if (lit == true_literal) {
+                            is_tautology = true;
+                        }
+                        else if (lit != false_literal) {                            
+                            lits.push_back(lit);
+                        }
+                        changed = true;
+                    }
+                    else {
+                        lits.push_back(l);
+                    }
+                }
+                if (changed) {
+                    to_delete.push_back(c);
+                    if (!is_tautology) {
+                        mk_clause(lits.size(), lits.c_ptr(), c->is_learned(), static_cast<_assumption_set>(c->assumptions()));
+                    }
+                }                        
+            }        
+            for (clause* c : to_delete) {
+                del_clause(c, m_clauses);
+            }
+        }
+
+        bool is_unit_ineq(clause const& c) const {
+            return 
+                c.size() == 1 &&
+                m_atoms[c[0].var()] && 
+                m_atoms[c[0].var()]->is_ineq_atom();
+        }
+
+        bool is_unit_eq(clause const& c) const {
+            return 
+                is_unit_ineq(c) && 
+                !c[0].sign() && 
+                m_atoms[c[0].var()]->is_eq();
+        }
+
+        /**
+           \brief determine whether the clause is a comparison v > k or v < k', where k >= 0 or k' <= 0.
+         */
+        lbool is_cmp0(clause const& c, var& v) {
+            if (!is_unit_ineq(c)) return l_undef;
+            literal lit = c[0];
+            ineq_atom const& a = *to_ineq_atom(m_atoms[lit.var()]);
+            bool sign = lit.sign();
+            poly * p0;
+            if (!is_single_poly(a, p0)) return l_undef;
+            if (m_pm.is_var(p0, v)) {
+                if (!sign && a.get_kind() == atom::GT) {
+                    return l_true;
+                }
+                if (!sign && a.get_kind() == atom::LT) {
+                    return l_false;
+                }
+                return l_undef;
+            }
+            polynomial::scoped_numeral n(m_pm.m());
+            if (m_pm.is_var_num(p0, v, n)) {
+                // x - k > 0
+                if (!sign && a.get_kind() == atom::GT && m_pm.m().is_nonneg(n)) {
+                    return l_true;
+                }
+                // x + k < 0
+                if (!sign && a.get_kind() == atom::LT && m_pm.m().is_nonpos(n)) {
+                    return l_false;
+                }
+                // !(x + k > 0)
+                if (sign && a.get_kind() == atom::GT && m_pm.m().is_pos(n)) {
+                    return l_false;
+                }
+                // !(x - k < 0)
+                if (sign && a.get_kind() == atom::LT && m_pm.m().is_neg(n)) {
+                    return l_true;
+                }
+            }
+            return l_undef;
+        }
+
+        bool is_single_poly(ineq_atom const& a, poly*& p) {
+            unsigned sz = a.size();
+            return sz == 1 && a.is_odd(0) && (p = a.p(0), true);
+        }
+
+        svector<lbool> m_var_signs;
+
+        void init_var_signs() {
+            m_var_signs.reset();
+            for (clause* cp : m_clauses) {
+                clause& c = *cp;
+                var x = 0;
+                lbool cmp = is_cmp0(c, x);
+                switch (cmp) {
+                case l_true:
+                    m_var_signs.setx(x, l_true, l_undef);
+                    break;
+                case l_false:
+                    m_var_signs.setx(x, l_false, l_undef);
+                    break;
+                default:
+                    break;                
+                }
+            }
+        }
+
+        /**
+           \brief returns true then c is an equality that is equivalent to v*p + q, 
+           and p > 0, v does not occur in p, q.
+        */
+        bool solve_var(clause& c, var& v, polynomial_ref& p, polynomial_ref& q) {
+            poly* p0;
+            if (!is_unit_eq(c)) return false;
+            ineq_atom & a = *to_ineq_atom(m_atoms[c[0].var()]);
+            if (!is_single_poly(a, p0)) return false;
+            var mx = max_var(p0);
+            if (mx >= m_is_int.size()) return false;
+            for (var x = 0; x <= mx; ++x) {
+                if (m_is_int[x]) continue;
+                if (1 == m_pm.degree(p0, x)) {                    
+                    p = m_pm.coeff(p0, x, 1, q);
+                    switch (m_pm.sign(p, m_var_signs)) {
+                    case l_true:
+                        v = x;
+                        return true;
+                    case l_false:
+                        v = x;
+                        p = -p;
+                        q = -q;
+                        return true;
+                    default:
+                        break;
+                    }
+                }
+            }
+            return false;
+        }        
+
         // -----------------------
         //
         // Pretty printing
         //
         // -----------------------
         
-        void display_num_assignment(std::ostream & out, display_var_proc const & proc) const {
+        std::ostream& display_num_assignment(std::ostream & out, display_var_proc const & proc) const {
             for (var x = 0; x < num_vars(); x++) {
                 if (m_assignment.is_assigned(x)) {
                     proc(out, x);
@@ -2287,9 +2561,10 @@ namespace nlsat {
                     out << "\n";
                 }
             }
+            return out;
         }
 
-        void display_bool_assignment(std::ostream & out) const {
+        std::ostream& display_bool_assignment(std::ostream & out) const {
             unsigned sz = m_atoms.size();
             for (bool_var b = 0; b < sz; b++) {
                 if (m_atoms[b] == 0 && m_bvalues[b] != l_undef) {
@@ -2300,6 +2575,7 @@ namespace nlsat {
                   for (bool_var b = 0; b < sz; b++) {
                       out << "b" << b << " -> " << m_bvalues[b] << " " << m_atoms[b] << "\n";
                   });
+            return out;
         }
 
         bool display_mathematica_assignment(std::ostream & out) const {
@@ -2317,16 +2593,17 @@ namespace nlsat {
             return !first;
         }
 
-        void display_num_assignment(std::ostream & out) const { 
-            display_num_assignment(out, m_display_var);
+        std::ostream& display_num_assignment(std::ostream & out) const { 
+            return display_num_assignment(out, m_display_var);
         }
 
-        void display_assignment(std::ostream& out) const {
+        std::ostream& display_assignment(std::ostream& out) const {
             display_bool_assignment(out);
             display_num_assignment(out);
+            return out;
         }
        
-        void display(std::ostream & out, ineq_atom const & a, display_var_proc const & proc, bool use_star = false) const {
+        std::ostream& display(std::ostream & out, ineq_atom const & a, display_var_proc const & proc, bool use_star = false) const {
             unsigned sz = a.size();
             for (unsigned i = 0; i < sz; i++) {
                 if (use_star && i > 0)
@@ -2346,9 +2623,10 @@ namespace nlsat {
             case atom::EQ: out << " = 0"; break;
             default: UNREACHABLE(); break;
             }
+            return out;
         }
-
-        void display_mathematica(std::ostream & out, ineq_atom const & a) const {
+        
+        std::ostream& display_mathematica(std::ostream & out, ineq_atom const & a) const {
             unsigned sz = a.size();
             for (unsigned i = 0; i < sz; i++) {
                 if (i > 0)
@@ -2370,9 +2648,10 @@ namespace nlsat {
             case atom::EQ: out << " == 0"; break;
             default: UNREACHABLE(); break;
             }
+            return out;
         }
 
-        void display_smt2(std::ostream & out, ineq_atom const & a, display_var_proc const & proc) const {
+        std::ostream& display_smt2(std::ostream & out, ineq_atom const & a, display_var_proc const & proc) const {
             switch (a.get_kind()) {
             case atom::LT: out << "(< "; break;
             case atom::GT: out << "(> "; break;
@@ -2398,9 +2677,10 @@ namespace nlsat {
             if (sz > 1)
                 out << ")";
             out << " 0)";
+            return out;
         }
         
-        void display(std::ostream & out, root_atom const & a, display_var_proc const & proc) const {
+        std::ostream& display(std::ostream & out, root_atom const & a, display_var_proc const & proc) const {
             proc(out, a.x());
             switch (a.get_kind()) {
             case atom::ROOT_LT: out << " < "; break;
@@ -2413,21 +2693,22 @@ namespace nlsat {
             out << "root[" << a.i() << "](";
             m_pm.display(out, a.p(), proc);
             out << ")";
+            return out;
         }
 
         struct mathematica_var_proc : public display_var_proc {
             var m_x;
         public:
             mathematica_var_proc(var x):m_x(x) {}
-            virtual void operator()(std::ostream & out, var x) const { 
+            virtual std::ostream& operator()(std::ostream & out, var x) const { 
                 if (m_x == x)
-                    out << "#1";
+                    return out << "#1";
                 else
-                    out << "x" << x; 
+                    return out << "x" << x; 
             }
         };
 
-        void display_mathematica(std::ostream & out, root_atom const & a) const {
+        std::ostream& display_mathematica(std::ostream & out, root_atom const & a) const {
             out << "x" << a.x();
             switch (a.get_kind()) {
             case atom::ROOT_LT: out << " < "; break;
@@ -2440,65 +2721,74 @@ namespace nlsat {
             out << "Root[";
             m_pm.display(out, a.p(), mathematica_var_proc(a.x()), true);
             out << " &, " << a.i() << "]";
+            return out;
         }
 
-        void display_smt2(std::ostream & out, root_atom const & a) const {
+        std::ostream& display_smt2(std::ostream & out, root_atom const & a) const {
             NOT_IMPLEMENTED_YET();
+            return out;
         }
         
-        void display(std::ostream & out, atom const & a, display_var_proc const & proc) const {
+        std::ostream& display(std::ostream & out, atom const & a, display_var_proc const & proc) const {
             if (a.is_ineq_atom())
-                display(out, static_cast<ineq_atom const &>(a), proc);
+                return display(out, static_cast<ineq_atom const &>(a), proc);
             else
-                display(out, static_cast<root_atom const &>(a), proc);
+                return display(out, static_cast<root_atom const &>(a), proc);
         }
 
-        void display_mathematica(std::ostream & out, atom const & a) const {
-            if (a.is_ineq_atom())
-                display_mathematica(out, static_cast<ineq_atom const &>(a));
-            else
-                display_mathematica(out, static_cast<root_atom const &>(a));
+        std::ostream& display(std::ostream & out, atom const & a) const {
+            return display(out, a, m_display_var);
         }
 
-        void display_smt2(std::ostream & out, atom const & a, display_var_proc const & proc) const {
+        std::ostream& display_mathematica(std::ostream & out, atom const & a) const {
             if (a.is_ineq_atom())
-                display_smt2(out, static_cast<ineq_atom const &>(a), proc);
+                return display_mathematica(out, static_cast<ineq_atom const &>(a));
             else
-                display_smt2(out, static_cast<root_atom const &>(a), proc);
+                return display_mathematica(out, static_cast<root_atom const &>(a));
         }
 
-        void display_atom(std::ostream & out, bool_var b, display_var_proc const & proc) const {
+        std::ostream& display_smt2(std::ostream & out, atom const & a, display_var_proc const & proc) const {
+            if (a.is_ineq_atom())
+                return display_smt2(out, static_cast<ineq_atom const &>(a), proc);
+            else
+                return display_smt2(out, static_cast<root_atom const &>(a), proc);
+        }
+
+        std::ostream& display_atom(std::ostream & out, bool_var b, display_var_proc const & proc) const {
             if (b == 0)
                 out << "true";
             else if (m_atoms[b] == 0)
                 out << "b" << b;
             else
                 display(out, *(m_atoms[b]), proc);
+            return out;
         }
 
-        void display_atom(std::ostream & out, bool_var b) const {
-            display_atom(out, b, m_display_var);
+        std::ostream& display_atom(std::ostream & out, bool_var b) const {
+            return display_atom(out, b, m_display_var);
         }
 
-        void display_mathematica_atom(std::ostream & out, bool_var b) const {
+        std::ostream& display_mathematica_atom(std::ostream & out, bool_var b) const {
             if (b == 0)
                 out << "(0 < 1)";
             else if (m_atoms[b] == 0)
                 out << "b" << b;
             else
                 display_mathematica(out, *(m_atoms[b]));
+            return out;
         }
 
-        void display_smt2_atom(std::ostream & out, bool_var b, display_var_proc const & proc) const {
+        std::ostream& display_smt2_atom(std::ostream & out, bool_var b, display_var_proc const & proc) const {
             if (b == 0)
                 out << "true";
             else if (m_atoms[b] == 0)
                 out << "b" << b;
             else
                 display_smt2(out, *(m_atoms[b]), proc);
+            return out;
         }
 
-        void display(std::ostream & out, literal l, display_var_proc const & proc) const {
+        std::ostream& display(std::ostream & out, literal l, display_var_proc const & proc) const {
             if (l.sign()) {
                 bool_var b = l.var();
                 out << "!";
@@ -2511,13 +2801,14 @@ namespace nlsat {
             else {
                 display_atom(out, l.var(), proc);
             }
+            return out;
         }
 
-        void display(std::ostream & out, literal l) const {
-            display(out, l, m_display_var);
+        std::ostream& display(std::ostream & out, literal l) const {
+            return display(out, l, m_display_var);
         }
 
-        void display_mathematica(std::ostream & out, literal l) const {
+        std::ostream& display_mathematica(std::ostream & out, literal l) const {
             if (l.sign()) {
                 bool_var b = l.var();
                 out << "!";
@@ -2530,9 +2821,10 @@ namespace nlsat {
             else {
                 display_mathematica_atom(out, l.var());
             }
+            return out;
         }
 
-        void display_smt2(std::ostream & out, literal l, display_var_proc const & proc) const {
+        std::ostream& display_smt2(std::ostream & out, literal l, display_var_proc const & proc) const {
             if (l.sign()) {
                 bool_var b = l.var();
                 out << "(not ";
@@ -2542,41 +2834,43 @@ namespace nlsat {
             else {
                 display_smt2_atom(out, l.var(), proc);
             }
+            return out;
         }
             
-        void display_assumptions(std::ostream & out, _assumption_set s) const {
-            
+        std::ostream& display_assumptions(std::ostream & out, _assumption_set s) const {
+            return out;
         }
         
-        void display(std::ostream & out, unsigned num, literal const * ls, display_var_proc const & proc) const {
+        std::ostream& display(std::ostream & out, unsigned num, literal const * ls, display_var_proc const & proc) const {
             for (unsigned i = 0; i < num; i++) {
                 if (i > 0)
                     out << " or ";
                 display(out, ls[i], proc);
             }
+            return out;
         }
 
-        void display(std::ostream & out, unsigned num, literal const * ls) {
-            display(out, num, ls, m_display_var);
+        std::ostream& display(std::ostream & out, unsigned num, literal const * ls) {
+            return display(out, num, ls, m_display_var);
+        }
+        
+        std::ostream& display(std::ostream & out, scoped_literal_vector const & cs) {
+            return display(out, cs.size(), cs.c_ptr(), m_display_var);
         }
 
-        void display(std::ostream & out, scoped_literal_vector const & cs) {
-            display(out, cs.size(), cs.c_ptr(), m_display_var);
-        }
-
-        void display(std::ostream & out, clause const & c, display_var_proc const & proc) const {
+        std::ostream& display(std::ostream & out, clause const & c, display_var_proc const & proc) const {
             if (c.assumptions() != 0) {
                 display_assumptions(out, static_cast<_assumption_set>(c.assumptions()));
                 out << " |- ";
             }
-            display(out, c.size(), c.c_ptr(), proc);
+            return display(out, c.size(), c.c_ptr(), proc);
         }
 
-        void display(std::ostream & out, clause const & c) const {
-            display(out, c, m_display_var);
+        std::ostream& display(std::ostream & out, clause const & c) const {
+            return display(out, c, m_display_var);
         }
 
-        void display_smt2(std::ostream & out, unsigned num, literal const * ls, display_var_proc const & proc) const {
+        std::ostream& display_smt2(std::ostream & out, unsigned num, literal const * ls, display_var_proc const & proc) const {
             if (num == 0) {
                 out << "false";
             }
@@ -2591,13 +2885,14 @@ namespace nlsat {
                 }
                 out << ")";
             }
+            return out;
         }
 
-        void display_smt2(std::ostream & out, clause const & c, display_var_proc const & proc = display_var_proc()) const {
-            display_smt2(out, c.size(), c.c_ptr(), proc);
+        std::ostream& display_smt2(std::ostream & out, clause const & c, display_var_proc const & proc = display_var_proc()) const {
+            return display_smt2(out, c.size(), c.c_ptr(), proc);
         }
 
-        void display_abst(std::ostream & out, literal l) const {
+        std::ostream& display_abst(std::ostream & out, literal l) const {
             if (l.sign()) {
                 bool_var b = l.var();
                 out << "!";
@@ -2609,25 +2904,27 @@ namespace nlsat {
             else {
                 out << "b" << l.var();
             }
+            return out;
         }
 
-        void display_abst(std::ostream & out, unsigned num, literal const * ls) const {
+        std::ostream& display_abst(std::ostream & out, unsigned num, literal const * ls) const {
             for (unsigned i = 0; i < num; i++) {
                 if (i > 0)
                     out << " or ";
                 display_abst(out, ls[i]);
             }
+            return out;
         }
 
-        void display_abst(std::ostream & out, scoped_literal_vector const & cs) const {
-            display_abst(out, cs.size(), cs.c_ptr());
+        std::ostream& display_abst(std::ostream & out, scoped_literal_vector const & cs) const {
+            return display_abst(out, cs.size(), cs.c_ptr());
         }
 
-        void display_abst(std::ostream & out, clause const & c) const {
-            display_abst(out, c.size(), c.c_ptr());
+        std::ostream& display_abst(std::ostream & out, clause const & c) const {
+            return display_abst(out, c.size(), c.c_ptr());
         }
 
-        void display_mathematica(std::ostream & out, clause const & c) const {
+        std::ostream& display_mathematica(std::ostream & out, clause const & c) const {
             out << "(";
             unsigned sz = c.size();
             for (unsigned i = 0; i < sz; i++) {
@@ -2636,12 +2933,13 @@ namespace nlsat {
                 display_mathematica(out, c[i]);
             }
             out << ")";
+            return out;
         }
 
         // Debugging support:
         // Display generated lemma in Mathematica format.
         // Mathematica must reduce lemma to True (modulo resource constraints).
-        void display_mathematica_lemma(std::ostream & out, unsigned num, literal const * ls, bool include_assignment = false) const {
+        std::ostream& display_mathematica_lemma(std::ostream & out, unsigned num, literal const * ls, bool include_assignment = false) const {
             out << "Resolve[ForAll[{";
             // var definition
             for (unsigned i = 0; i < num_vars(); i++) {
@@ -2662,71 +2960,69 @@ namespace nlsat {
                 display_mathematica(out, ls[i]);
             }
             out << "], Reals]\n"; // end of exists
+            return out;
         }
         
-        void display(std::ostream & out, clause_vector const & cs, display_var_proc const & proc) const {
-            unsigned sz = cs.size();
-            for (unsigned i = 0; i < sz; i++) {
-                display(out, *(cs[i]), proc);
-                out << "\n";
+        std::ostream& display(std::ostream & out, clause_vector const & cs, display_var_proc const & proc) const {
+            for (clause* c : cs) {
+                display(out, *c, proc) << "\n";
             }
+            return out;
         }
 
-        void display(std::ostream & out, clause_vector const & cs) const {
-            display(out, cs, m_display_var);
+        std::ostream& display(std::ostream & out, clause_vector const & cs) const {
+            return display(out, cs, m_display_var);
         }
 
-        void display_mathematica(std::ostream & out, clause_vector const & cs) const {
+        std::ostream& display_mathematica(std::ostream & out, clause_vector const & cs) const {            
             unsigned sz = cs.size();
             for (unsigned i = 0; i < sz; i++) {
                 if (i > 0) out << ",\n";
-                out << "  ";
-                display_mathematica(out, *(cs[i]));
+                display_mathematica(out << " ", *(cs[i]));
             }
+            return out;
         }
 
-        void display_abst(std::ostream & out, clause_vector const & cs) const {
-            unsigned sz = cs.size();
-            for (unsigned i = 0; i < sz; i++) {
-                display_abst(out, *(cs[i]));
-                out << "\n";
+        std::ostream& display_abst(std::ostream & out, clause_vector const & cs) const {
+            for (clause* c : cs) {
+                display_abst(out, *c) << "\n";
             }
+            return out;
         }
 
-        void display(std::ostream & out, display_var_proc const & proc) const {
+        std::ostream& display(std::ostream & out, display_var_proc const & proc) const {
             display(out, m_clauses, proc);
             if (!m_learned.empty()) {
-                out << "Lemmas:\n";
-                display(out, m_learned, proc);
+                display(out << "Lemmas:\n", m_learned, proc);
             }
+            return out;
         }
 
-        void display_mathematica(std::ostream & out) const {
-            out << "{\n";
-            display_mathematica(out, m_clauses);
-            out << "}\n";
+        std::ostream& display_mathematica(std::ostream & out) const {
+            return display_mathematica(out << "{\n", m_clauses) << "}\n";
         }
 
-        void display_abst(std::ostream & out) const {
+        std::ostream& display_abst(std::ostream & out) const {
             display_abst(out, m_clauses);
             if (!m_learned.empty()) {
-                out << "Lemmas:\n";
-                display_abst(out, m_learned);
+                display_abst(out << "Lemmas:\n", m_learned);
             }
+            return out;
         }
 
-        void display(std::ostream & out) const {
+        std::ostream& display(std::ostream & out) const {
             display(out, m_display_var);
-            display_assignment(out);
+            return display_assignment(out);
         }
 
-        void display_vars(std::ostream & out) const {
+        std::ostream& display_vars(std::ostream & out) const {
             for (unsigned i = 0; i < num_vars(); i++) {
                 out << i << " -> "; m_display_var(out, i); out << "\n";
             }
+            return out;
         }
 
-        void display_smt2_arith_decls(std::ostream & out) const {
+        std::ostream& display_smt2_arith_decls(std::ostream & out) const {
             unsigned sz = m_is_int.size();
             for (unsigned i = 0; i < sz; i++) {
                 if (m_is_int[i])
@@ -2734,31 +3030,32 @@ namespace nlsat {
                 else
                     out << "(declare-fun x" << i << " () Real)\n";
             }
+            return out;
         }
 
-        void display_smt2_bool_decls(std::ostream & out) const {
+        std::ostream& display_smt2_bool_decls(std::ostream & out) const {
             unsigned sz = m_atoms.size();
             for (unsigned i = 0; i < sz; i++) {
                 if (m_atoms[i] == 0)
                     out << "(declare-fun b" << i << " () Bool)\n";
             }
+            return out;
         }
 
-        void display_smt2(std::ostream & out) const {
+        std::ostream& display_smt2(std::ostream & out) const {
             display_smt2_bool_decls(out);
             display_smt2_arith_decls(out);
             out << "(assert (and true\n";
-            unsigned sz = m_clauses.size();
-            for (unsigned i = 0; i < sz; i++) {
-                display_smt2(out, *(m_clauses[i]));
-                out << "\n";
+            for (clause* c : m_clauses) {
+                display_smt2(out, *c) << "\n";
             }
             out << "))\n(check-sat)" << std::endl;
+            return out;
         }
     };
     
-    solver::solver(reslimit& rlim, params_ref const & p) {
-        m_imp = alloc(imp, *this, rlim, p);
+    solver::solver(reslimit& rlim, params_ref const & p, bool incremental) {
+        m_imp = alloc(imp, *this, rlim, p, incremental);
     }
         
     solver::~solver() {
diff --git a/src/nlsat/nlsat_solver.h b/src/nlsat/nlsat_solver.h
index ac503c603..c632a0334 100644
--- a/src/nlsat/nlsat_solver.h
+++ b/src/nlsat/nlsat_solver.h
@@ -35,7 +35,7 @@ namespace nlsat {
         struct imp;
         imp * m_imp;
     public:
-        solver(reslimit& rlim, params_ref const & p);
+        solver(reslimit& rlim, params_ref const & p, bool incremental);
         ~solver();
 
         /**
diff --git a/src/nlsat/nlsat_types.h b/src/nlsat/nlsat_types.h
index db3b6ffa5..647a5e3ee 100644
--- a/src/nlsat/nlsat_types.h
+++ b/src/nlsat/nlsat_types.h
@@ -47,6 +47,8 @@ namespace nlsat {
     typedef polynomial::var_vector     var_vector;
     typedef polynomial::manager        pmanager;
     typedef polynomial::polynomial     poly;
+    typedef polynomial::monomial       monomial;
+    typedef polynomial::numeral        numeral;
     const var null_var = polynomial::null_var;
     
     const var true_bool_var = 0;
diff --git a/src/nlsat/tactic/nlsat_tactic.cpp b/src/nlsat/tactic/nlsat_tactic.cpp
index 510f503e7..9de7215e3 100644
--- a/src/nlsat/tactic/nlsat_tactic.cpp
+++ b/src/nlsat/tactic/nlsat_tactic.cpp
@@ -32,11 +32,11 @@ class nlsat_tactic : public tactic {
         ast_manager & m;
         expr_ref_vector m_var2expr;
         expr_display_var_proc(ast_manager & _m):m(_m), m_var2expr(_m) {}
-        virtual void operator()(std::ostream & out, nlsat::var x) const { 
+        virtual std::ostream& operator()(std::ostream & out, nlsat::var x) const { 
             if (x < m_var2expr.size())
-                out << mk_ismt2_pp(m_var2expr.get(x), m); 
+                return out << mk_ismt2_pp(m_var2expr.get(x), m); 
             else
-                out << "x!" << x;
+                return out << "x!" << x;
         }
     };
 
@@ -51,7 +51,7 @@ class nlsat_tactic : public tactic {
             m(_m),
             m_params(p),
             m_display_var(_m),
-            m_solver(m.limit(), p) {
+            m_solver(m.limit(), p, false) {
         }
         
         void updt_params(params_ref const & p) {
diff --git a/src/nlsat/tactic/qfnra_nlsat_tactic.cpp b/src/nlsat/tactic/qfnra_nlsat_tactic.cpp
index 22f64af47..47b9e0505 100644
--- a/src/nlsat/tactic/qfnra_nlsat_tactic.cpp
+++ b/src/nlsat/tactic/qfnra_nlsat_tactic.cpp
@@ -48,11 +48,15 @@ tactic * mk_qfnra_nlsat_tactic(ast_manager & m, params_ref const & p) {
                                           purify_p),
                              mk_propagate_values_tactic(m, p),
                              mk_solve_eqs_tactic(m, p),
+                             using_params(mk_purify_arith_tactic(m, p),
+                                          purify_p),
                              mk_elim_uncnstr_tactic(m, p),
                              mk_elim_term_ite_tactic(m, p)),
                     and_then(/* mk_degree_shift_tactic(m, p), */ // may affect full dimensionality detection
                              factor,
                              mk_solve_eqs_tactic(m, p),
+                             using_params(mk_purify_arith_tactic(m, p),
+                                          purify_p),
                              using_params(mk_simplify_tactic(m, p),
                                           main_p),
                              mk_tseitin_cnf_core_tactic(m, p),
diff --git a/src/qe/nlqsat.cpp b/src/qe/nlqsat.cpp
index ab5095328..e28f6ae4f 100644
--- a/src/qe/nlqsat.cpp
+++ b/src/qe/nlqsat.cpp
@@ -782,7 +782,7 @@ namespace qe {
             m(m),
             m_mode(mode),
             m_params(p),
-            m_solver(m.limit(), p),
+            m_solver(m.limit(), p, true),
             m_nftactic(0),
             m_rmodel(m_solver.am()),
             m_rmodel0(m_solver.am()),
diff --git a/src/tactic/core/solve_eqs_tactic.cpp b/src/tactic/core/solve_eqs_tactic.cpp
index 65d474182..5df523ff7 100644
--- a/src/tactic/core/solve_eqs_tactic.cpp
+++ b/src/tactic/core/solve_eqs_tactic.cpp
@@ -44,6 +44,7 @@ class solve_eqs_tactic : public tactic {
         expr_sparse_mark              m_candidate_set;
         ptr_vector<expr>              m_candidates;
         ptr_vector<app>               m_vars;
+        expr_sparse_mark              m_nonzero;
         ptr_vector<app>               m_ordered_vars;
         bool                          m_produce_proofs;
         bool                          m_produce_unsat_cores;
@@ -55,8 +56,7 @@ class solve_eqs_tactic : public tactic {
             m_r_owner(r == 0 || owner),
             m_a_util(m),
             m_num_steps(0),
-            m_num_eliminated_vars(0)
-             {
+            m_num_eliminated_vars(0) {
             updt_params(p);
             if (m_r == 0)
                 m_r       = mk_default_expr_replacer(m);
@@ -78,7 +78,7 @@ class solve_eqs_tactic : public tactic {
         void checkpoint() {
             if (m().canceled())
                 throw tactic_exception(m().limit().get_cancel_msg());
-        cooperate("solve-eqs");
+            cooperate("solve-eqs");
         }
         
         // Check if the number of occurrences of t is below the specified threshold :solve-eqs-max-occs
@@ -106,7 +106,8 @@ class solve_eqs_tactic : public tactic {
             }
         }
         bool trivial_solve(expr * lhs, expr * rhs, app_ref & var, expr_ref & def, proof_ref & pr) {
-            if (trivial_solve1(lhs, rhs, var, def, pr)) return true;
+            if (trivial_solve1(lhs, rhs, var, def, pr)) 
+                return true;
             if (trivial_solve1(rhs, lhs, var, def, pr)) {
                 if (m_produce_proofs) {
                     pr = m().mk_commutativity(m().mk_eq(lhs, rhs));
@@ -187,6 +188,77 @@ class solve_eqs_tactic : public tactic {
             }
             return false;
         }
+
+        void add_pos(expr* f) {
+            expr* lhs = nullptr, *rhs = nullptr;
+            rational val;
+            if (m_a_util.is_le(f, lhs, rhs) && m_a_util.is_numeral(rhs, val) && val.is_neg()) {
+                m_nonzero.mark(lhs);
+            }
+            else if (m_a_util.is_ge(f, lhs, rhs) && m_a_util.is_numeral(rhs, val) && val.is_pos()) {
+                m_nonzero.mark(lhs);
+            }
+            else if (m().is_not(f, f)) {
+                if (m_a_util.is_le(f, lhs, rhs) && m_a_util.is_numeral(rhs, val) && !val.is_neg()) {
+                    m_nonzero.mark(lhs);
+                }
+                else if (m_a_util.is_ge(f, lhs, rhs) && m_a_util.is_numeral(rhs, val) && !val.is_pos()) {
+                    m_nonzero.mark(lhs);
+                }
+                else if (m().is_eq(f, lhs, rhs) && m_a_util.is_numeral(rhs, val) && val.is_zero()) {
+                    m_nonzero.mark(lhs);                    
+                }
+            }            
+        }
+
+        bool is_nonzero(expr* e) {
+            return m_nonzero.is_marked(e);
+        }
+
+        bool isolate_var(app* arg, app_ref& var, expr_ref& div, unsigned i, app* lhs, expr* rhs) {
+            if (!m_a_util.is_mul(arg)) return false;
+            unsigned n = arg->get_num_args();
+            for (unsigned j = 0; j < n; ++j) {
+                expr* e = arg->get_arg(j);
+                bool ok = is_uninterp_const(e) && check_occs(e) && !occurs(e, rhs) && !occurs_except(e, lhs, i);
+                if (!ok) continue;
+                var = to_app(e);
+                for (unsigned k = 0; ok && k < n; ++k) {
+                    expr* arg_k = arg->get_arg(k);
+                    ok = k == j || (!occurs(var, arg_k) && is_nonzero(arg_k));
+                }
+                if (!ok) continue;
+                ptr_vector<expr> args;
+                for (unsigned k = 0; k < n; ++k) {
+                    if (k != j) args.push_back(arg->get_arg(k));
+                }
+                div = m_a_util.mk_mul(args.size(), args.c_ptr());
+                return true;
+            }
+            return false;
+        }
+
+        bool solve_nl(app * lhs, expr * rhs, expr* eq, app_ref& var, expr_ref & def, proof_ref & pr) {
+            SASSERT(m_a_util.is_add(lhs));
+            if (m_a_util.is_int(lhs)) return false;
+            unsigned num = lhs->get_num_args();
+            expr_ref div(m());
+            for (unsigned i = 0; i < num; i++) {
+                expr * arg = lhs->get_arg(i);
+                if (is_app(arg) && isolate_var(to_app(arg), var, div, i, lhs, rhs)) {
+                    ptr_vector<expr> args;
+                    for (unsigned k = 0; k < num; ++k) {
+                        if (k != i) args.push_back(lhs->get_arg(k));
+                    }
+                    def = m_a_util.mk_sub(rhs, m_a_util.mk_add(args.size(), args.c_ptr()));
+                    def = m_a_util.mk_div(def, div);
+                    if (m_produce_proofs)
+                        pr = m().mk_rewrite(eq, m().mk_eq(var, def));
+                    return true;
+                }
+            }
+            return false;
+        }
         
         bool solve_arith_core(app * lhs, expr * rhs, expr * eq, app_ref & var, expr_ref & def, proof_ref & pr) {
             SASSERT(m_a_util.is_add(lhs));
@@ -204,7 +276,8 @@ class solve_eqs_tactic : public tactic {
                     break;
                 }
                 else if (m_a_util.is_mul(arg, a, v) && 
-                         is_uninterp_const(v) && !m_candidate_vars.is_marked(v) &&
+                         is_uninterp_const(v) && 
+                         !m_candidate_vars.is_marked(v) &&
                          m_a_util.is_numeral(a, a_val) &&
                          !a_val.is_zero() &&
                          (!is_int || a_val.is_minus_one()) &&
@@ -252,16 +325,20 @@ class solve_eqs_tactic : public tactic {
             return 
                 (m_a_util.is_add(lhs) && solve_arith_core(to_app(lhs), rhs, eq, var, def, pr)) ||
                 (m_a_util.is_add(rhs) && solve_arith_core(to_app(rhs), lhs, eq, var, def, pr));
+#if 0
+            // better done inside of nlsat
+                (m_a_util.is_add(lhs) && solve_nl(to_app(lhs), rhs, eq, var, def, pr)) ||
+                (m_a_util.is_add(rhs) && solve_nl(to_app(rhs), lhs, eq, var, def, pr));
+#endif
         }
         
         bool solve(expr * f, app_ref & var, expr_ref & def, proof_ref & pr) {
-            if (m().is_eq(f)) {
-                if (trivial_solve(to_app(f)->get_arg(0), to_app(f)->get_arg(1), var, def, pr))
+            expr* arg1 = 0, *arg2 = 0;
+            if (m().is_eq(f, arg1, arg2)) {
+                if (trivial_solve(arg1, arg2, var, def, pr))
                     return true;
                 if (m_theory_solver) {
-                    expr * lhs = to_app(f)->get_arg(0);
-                    expr * rhs = to_app(f)->get_arg(1);
-                    if (solve_arith(lhs, rhs, f, var, def, pr))
+                    if (solve_arith(arg1, arg2, f, var, def, pr))
                         return true;
                 }
                 return false;
@@ -321,11 +398,14 @@ class solve_eqs_tactic : public tactic {
             m_candidate_set.reset();
             m_candidates.reset();
             m_vars.reset();
-            
+            m_nonzero.reset();
             app_ref  var(m());
             expr_ref  def(m());
             proof_ref pr(m());
             unsigned size = g.size();
+            for (unsigned idx = 0; idx < size; idx++) {
+                add_pos(g.form(idx));
+            }
             for (unsigned idx = 0; idx < size; idx++) {
                 checkpoint();
                 expr * f = g.form(idx);
@@ -347,10 +427,8 @@ class solve_eqs_tactic : public tactic {
             
             TRACE("solve_eqs", 
                   tout << "candidate vars:\n";
-                  ptr_vector<app>::iterator it = m_vars.begin();
-                  ptr_vector<app>::iterator end = m_vars.end();
-                  for (; it != end; ++it) {
-                      tout << mk_ismt2_pp(*it, m()) << " ";
+                  for (app* v : m_vars) {
+                      tout << mk_ismt2_pp(v, m()) << " ";
                   }
                   tout << "\n";);
         }
@@ -374,12 +452,9 @@ class solve_eqs_tactic : public tactic {
             
             typedef std::pair<expr *, unsigned> frame;
             svector<frame> todo;
-            ptr_vector<app>::const_iterator it  = m_vars.begin();
-            ptr_vector<app>::const_iterator end = m_vars.end();
-            unsigned num;
-            for (; it != end; ++it) {
+            unsigned num = 0;
+            for (app* v : m_vars) {
                 checkpoint();
-                app * v = *it;
                 if (!m_candidate_vars.is_marked(v))
                     continue;
                 todo.push_back(frame(v, 0));
@@ -483,20 +558,19 @@ class solve_eqs_tactic : public tactic {
             }
             
             // cleanup
-            it  = m_vars.begin();
-            for (unsigned idx = 0; it != end; ++it, ++idx) {
-                if (!m_candidate_vars.is_marked(*it)) {
+            unsigned idx = 0;
+            for (expr* v : m_vars) {
+                if (!m_candidate_vars.is_marked(v)) {
                     m_candidate_set.mark(m_candidates[idx], false);
                 }
+                ++idx;
             }
             
             TRACE("solve_eqs", 
                   tout << "ordered vars:\n";
-                  ptr_vector<app>::iterator it = m_ordered_vars.begin();
-                  ptr_vector<app>::iterator end = m_ordered_vars.end();
-                  for (; it != end; ++it) {
-                      SASSERT(m_candidate_vars.is_marked(*it));
-                      tout << mk_ismt2_pp(*it, m()) << " ";
+                  for (app* v : m_ordered_vars) {
+                      SASSERT(m_candidate_vars.is_marked(v));
+                      tout << mk_ismt2_pp(v, m()) << " ";
                   }
                   tout << "\n";);
             m_candidate_vars.reset();
@@ -609,10 +683,7 @@ class solve_eqs_tactic : public tactic {
             if (m_produce_models) {
                 if (mc.get() == 0)
                     mc = alloc(gmc, m());
-                ptr_vector<app>::iterator it  = m_ordered_vars.begin();
-                ptr_vector<app>::iterator end = m_ordered_vars.end();
-                for (; it != end; ++it) {
-                    app * v    = *it;
+                for (app* v : m_ordered_vars) {
                     expr * def = 0;
                     proof * pr;
                     expr_dependency * dep;