diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp
index 3c0eb9173..c102e447e 100644
--- a/src/math/dd/dd_pdd.cpp
+++ b/src/math/dd/dd_pdd.cpp
@@ -1908,6 +1908,36 @@ namespace dd {
             return m->mk_var(var())*h + l;
     }
 
+    void pdd::merge(unsigned_vector& lo, unsigned_vector& hi, unsigned_vector& common) {
+        std::sort(lo.begin(), lo.end());
+        std::sort(hi.begin(), hi.end());
+        unsigned ir = 0, jr = 0;
+        for (unsigned i = 0, j = 0; i < lo.size() || j < hi.size(); ) {
+            if (i == lo.size()) 
+                hi[jr++] = hi[j++];
+            else if (j == hi.size()) 
+                lo[ir++] = lo[i++];
+            else if (lo[i] == hi[j]) {
+                common.push_back(lo[i]);
+                ++i;
+                ++j;
+            }
+            else if (lo[i] > hi[j]) 
+                hi[jr++] = hi[j++];
+            else 
+                lo[ir++] = lo[i++];
+        }
+        lo.shrink(ir);
+        hi.shrink(jr);    
+    }
+
+    pdd pdd::mul(unsigned_vector const& vars) const {
+        pdd r = *this;
+        for (auto v : vars)
+            r = m->mk_var(v) * r;
+        return r;
+    }
+
     std::pair<unsigned_vector, pdd> pdd::var_factors() const {
         if (is_val())
             return { unsigned_vector(), *this };
@@ -1924,49 +1954,23 @@ namespace dd {
             return { unsigned_vector(), *this };
         
         unsigned_vector lo_and_hi;
-        auto merge = [&](unsigned_vector& lo_vars, unsigned_vector& hi_vars) {
-            unsigned ir = 0, jr = 0;
-            for (unsigned i = 0, j = 0; i < lo_vars.size() || j < hi_vars.size(); ) {
-                if (i == lo_vars.size()) 
-                    hi_vars[jr++] = hi_vars[j++];
-                else if (j == hi_vars.size()) 
-                    lo_vars[ir++] = lo_vars[i++];
-                else if (lo_vars[i] == hi_vars[j]) {
-                    lo_and_hi.push_back(lo_vars[i]);
-                    ++i;
-                    ++j;
-                }
-                else if (m->m_var2level[lo_vars[i]] > m->m_var2level[hi_vars[j]]) 
-                    hi_vars[jr++] = hi_vars[j++];
-                else 
-                    lo_vars[ir++] = lo_vars[i++];
-            }
-            lo_vars.shrink(ir);
-            hi_vars.shrink(jr);
-        };
-
-        auto mul = [&](unsigned_vector const& vars, pdd p) {
-            for (auto v : vars)
-                p *= m->mk_var(v);
-            return p;
-        };
 
         auto [hi_vars, p] = hi().var_factors();
         if (lo_vars.back() == v) {
             lo_vars.pop_back();
-            merge(lo_vars, hi_vars);
+            merge(lo_vars, hi_vars, lo_and_hi);
             lo_and_hi.push_back(v);
-            return { lo_and_hi, mul(lo_vars, q) + mul(hi_vars, p) };
+            return { lo_and_hi, q.mul(lo_vars) + p.mul(hi_vars) };
         }
         if (hi_vars.empty())
             return { unsigned_vector(), *this };
         
-        merge(lo_vars, hi_vars);
+        merge(lo_vars, hi_vars, lo_and_hi);
         hi_vars.push_back(v);
         if (lo_and_hi.empty())
             return { unsigned_vector(), *this };
         else 
-            return { lo_and_hi, mul(lo_vars, q) + mul(hi_vars, p) };                
+            return { lo_and_hi, q.mul(lo_vars) + p.mul(hi_vars) };                
     }
 
 
diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h
index b96e5f04b..896e6bccb 100644
--- a/src/math/dd/dd_pdd.h
+++ b/src/math/dd/dd_pdd.h
@@ -483,6 +483,8 @@ namespace dd {
          * \brief factor out variables
          */
         std::pair<unsigned_vector, pdd> var_factors() const;
+        static void merge(unsigned_vector &lo, unsigned_vector &hi, unsigned_vector& common);
+        pdd mul(unsigned_vector const &vars) const;
 
         pdd subst_val0(vector<std::pair<unsigned, rational>> const& s) const { return m->subst_val0(*this, s); }
         pdd subst_val(pdd const& s) const { VERIFY_EQ(m, s.m); return m->subst_val(*this, s); }
diff --git a/src/math/lp/nla_stellensatz.cpp b/src/math/lp/nla_stellensatz.cpp
index 87c0f3ab1..7d1313a20 100644
--- a/src/math/lp/nla_stellensatz.cpp
+++ b/src/math/lp/nla_stellensatz.cpp
@@ -6,6 +6,7 @@
   */
 
 #include "util/heap.h"
+#include "math/dd/pdd_eval.h"
 #include "math/lp/nla_core.h"
 #include "math/lp/nla_coi.h"
 #include "math/lp/nla_stellensatz.h"
@@ -29,16 +30,37 @@ namespace nla {
         return k1;
     }
 
+    lp::lconstraint_kind meet(lp::lconstraint_kind k1, lp::lconstraint_kind k2) {
+        if (k1 == k2)
+            return k1;
+        if (k1 == lp::lconstraint_kind::EQ)
+            return k1;
+        if (k2 == lp::lconstraint_kind::EQ)
+            return k2;
+        if ((k1 == lp::lconstraint_kind::LE && k2 == lp::lconstraint_kind::LT) ||
+            (k1 == lp::lconstraint_kind::LT && k2 == lp::lconstraint_kind::LE))
+            return lp::lconstraint_kind::LE;
+        if ((k1 == lp::lconstraint_kind::GE && k2 == lp::lconstraint_kind::GT) ||
+            (k1 == lp::lconstraint_kind::GT && k2 == lp::lconstraint_kind::GE))
+            return lp::lconstraint_kind::GE;
+        UNREACHABLE();
+        return k1;
+    }
+
     stellensatz::stellensatz(core* core) : 
         common(core), 
         m_solver(*this), 
-        m_coi(*core) 
+        m_coi(*core),
+        pddm(m_solver.lra().number_of_vars())
     {}
 
     lbool stellensatz::saturate() {
         init_solver();
-        TRACE(arith, display(tout << "stellensatz after saturation\n"));
+        TRACE(arith, display(tout << "stellensatz before saturation\n"));
         eliminate_variables();
+        // TODO: populate solver
+        TRACE(arith, display(tout << "stellensatz after saturation\n"));
+        return l_undef;
         lbool r = m_solver.solve();
         // IF_VERBOSE(0, verbose_stream() << "stellensatz " << r << "\n");
         if (r == l_false)
@@ -50,9 +72,7 @@ namespace nla {
         m_solver.init();
         m_vars2mon.reset();
         m_mon2vars.reset();
-        m_values.reset();
-        m_justifications.reset();
-        m_max_monomial_degree = 0;
+        m_constraints.reset();
         m_coi.init();
         init_vars();
         init_occurs();
@@ -60,34 +80,18 @@ namespace nla {
 
     void stellensatz::init_vars() {
         auto const& src = c().lra_solver();
-        auto &dst = m_solver.lra();
         auto sz = src.number_of_vars();
         for (unsigned v = 0; v < sz; ++v) {
-            // Declare v into m_solver.lra()
-            unsigned w;
             if (m_coi.mons().contains(v))
                 init_monomial(v);
-            else 
-                m_values.push_back(c().val(v));
-            if (m_coi.terms().contains(v)) {
-                auto const& t = src.get_term(v);
-                // justification: variables in coefficients are always declared before term variable.
-                SASSERT(all_of(t, [&](auto p) { return p.j() < v; }));
-                w = dst.add_term(t.coeffs_as_vector(), v);
-            } 
-            else
-                w = dst.add_var(v, src.var_is_int(v));
-
             // assert bounds on v in the new solver.
-            VERIFY(w == v);
             if (src.column_has_lower_bound(v)) {
                 auto lo_bound = src.get_lower_bound(v);
                 SASSERT(lo_bound.y >= 0);
                 auto k   = lo_bound.y > 0 ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE;
                 auto rhs = lo_bound.x;
                 auto dep = src.get_column_lower_bound_witness(v);
-                auto ci  = dst.add_var_bound(v, k, rhs);
-                add_justification(ci, external_justification(dep));
+                add_var_bound(v, k, rhs, external_justification(dep));
             }
             if (src.column_has_upper_bound(v)) {
                 auto hi_bound = src.get_upper_bound(v);
@@ -95,32 +99,97 @@ namespace nla {
                 auto k   = hi_bound.y < 0 ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE;
                 auto rhs = hi_bound.x;
                 auto dep = src.get_column_upper_bound_witness(v);
-                auto ci  = dst.add_var_bound(v, k, rhs);
-                add_justification(ci, external_justification(dep));
+                add_var_bound(v, k, rhs, external_justification(dep));
             }
         }
         m_occurs.reserve(sz);
-        for (auto const &[v, vars] : m_mon2vars) 
-            m_max_monomial_degree = std::max(m_max_monomial_degree, vars.size());        
+    }
+
+    dd::pdd stellensatz::to_pdd(lpvar v) {
+        if (m_coi.mons().contains(v)) {
+            auto& mon = c().emons()[v];
+            dd::pdd prod(pddm);
+            prod = 1;
+            for (auto w : mon.vars())
+                prod *= to_pdd(w);
+            return prod;
+        }
+        if (m_coi.terms().contains(v)) {
+            auto const& t = c().lra_solver().get_term(v);
+            dd::pdd sum(pddm);
+            sum = 0;
+            for (auto cv : t) 
+                sum *= to_pdd(cv.j())*cv.coeff();
+            return sum;
+        }
+        return pddm.mk_var(v);
+    }
+
+    stellensatz::term_offset stellensatz::to_term_offset(dd::pdd const &p) {
+        term_offset to;
+        for (auto const &[coeff, vars] : p) {
+            if (vars.empty())
+                to.second += coeff;
+            else
+                to.first.add_monomial(coeff, mk_monomial(vars));                
+        }        
+        return to;
+    }
+
+    lpvar stellensatz::mk_monomial(svector<lpvar> const &vars) {
+        lpvar v;
+        if (vars.size() == 1)
+            return vars[0];
+        svector<lpvar> _vars(vars);
+        std::sort(_vars.begin(), _vars.end());
+        if (m_vars2mon.find(_vars, v))
+            return v;
+
+        v = add_var(is_int(vars));
+        m_vars2mon.insert(_vars, v);
+        m_mon2vars.insert(v, _vars);
+        SASSERT(m_values.size() == v);
+        return v;
     }
 
     void stellensatz::init_monomial(unsigned mon_var) {
-        auto& mon = c().emons()[mon_var];
+        auto &mon = c().emons()[mon_var];
         svector<lpvar> vars(mon.vars());
         std::sort(vars.begin(), vars.end());
         m_vars2mon.insert(vars, mon_var);
         m_mon2vars.insert(mon_var, vars);
-        m_values.push_back(cvalue(vars));
+    }
+
+    lp::constraint_index stellensatz::add_var_bound(lp::lpvar v, lp::lconstraint_kind k, rational const& rhs, justification j) {
+        auto p = to_pdd(v) - rhs;
+        rational d(1);
+        for (auto const& [c, is_constant] : p.coefficients()) 
+            d = lcm(d, denominator(c));
+        if (d != 1)
+            p *= d;        
+        return gcd_normalize(add_constraint(p, k, j));
+    }
+
+    lp::constraint_index stellensatz::add_constraint(dd::pdd p, lp::lconstraint_kind k, justification j) {
+        if (k == lp::lconstraint_kind::LE) {
+            k = lp::lconstraint_kind::GE;
+            p = -p;
+        }
+        if (k == lp::lconstraint_kind::LT) {
+            k = lp::lconstraint_kind::GT;
+            p = -p;
+        }
+        m_constraints.push_back({p, k, j });
+        return m_constraints.size() - 1;
     }
 
     // initialize active set of constraints that evaluate to false in the current model
     // loop over active set to eliminate variables.
     void stellensatz::eliminate_variables() {
         vector<std::pair<lp::constraint_index, uint_set>> active;
-        for (auto ci : m_solver.lra().constraints().indices()) {
-            if (constraint_is_true(ci))
-                continue;
-            active.push_back({ci, uint_set()});
+        for (unsigned ci = 0; ci < m_constraints.size(); ++ci) {
+            if (!constraint_is_true(ci))
+                active.push_back({ci, uint_set()});
         } 
         while (!active.empty()) {
             // factor ci into x*p + q <= rhs, where degree(x, q) < degree(x, p) + 1
@@ -128,126 +197,126 @@ namespace nla {
             // otherwise
             // find complementary constraints that contain x^k for k = 0 .. degree -1
             // eliminate x (and other variables) by combining ci with complementary constraints. 
-            auto [ci, tabu] = active.back();      
+            auto [ci, tabu] = active.back();               
             active.pop_back();           
             auto x = select_variable_to_eliminate(ci);
             auto f = factor(x, ci);
             auto p_value = cvalue(f.p);
-            if (!f.p.is_empty() && p_value == 0) {
+            if (!f.p.is_zero() && p_value == 0) {
                 // add p = 0 as assumption and reduce to q.
-                auto p_is_0 = assume(f.p, lp::lconstraint_kind::EQ, rational(0));
-                if (!f.q.is_empty()) {
-                    auto const &con = m_solver.lra().constraints()[ci];
-                    // ci & -p_is_0*f.coeff*x^f.degree  => new_ci
-              
-                    if (f.coeff != -1)
-                        p_is_0 = multiply(p_is_0, rational(-f.coeff));
+                auto p_is_0 = assume(f.p, lp::lconstraint_kind::EQ);
+                if (!f.q.is_zero()) {
+                    // ci & -p_is_0*x^f.degree  => new_ci
+                    p_is_0 = multiply(p_is_0, rational(-1));             
                     for (unsigned i = 0; i < f.degree; ++i)
-                        p_is_0 = multiply_var(p_is_0, x);                    
+                        p_is_0 = multiply_var(p_is_0, x);                  
                     auto new_ci = add(ci, p_is_0);
+                    init_occurs(new_ci);
                     uint_set new_tabu(tabu);
                     uint_set q_vars;
-                    for (auto cv : f.q)
-                        q_vars.insert(cv.j());
-                    for (auto cv : f.p)
-                        if (!q_vars.contains(cv.j()))
-                            new_tabu.insert(cv.j());
+                    for (auto j : f.q.free_vars())
+                        q_vars.insert(j);
+                    for (auto j : f.p.free_vars())
+                        if (!q_vars.contains(j))
+                            new_tabu.insert(j);
                     active.push_back({new_ci, new_tabu});
                 }
                 continue;
             }
             for (auto other_ci : m_occurs[x]) {
-            
-            }
+                if (other_ci == ci)
+                    continue;
+                auto const &other_ineq = m_constraints[other_ci];
+                auto g = factor(x, other_ci);
+                auto p_value2 = cvalue(g.p);
+                // check that p_value and p_value2 have different signs
+                // check that other_ineq.lhs() is disjoint from tabu
+                // compute overlaps extending x
+                // 
+                SASSERT(g.degree > 0);
+                if (g.degree > f.degree) // future: could handle this too by considering tabu to be a map into degrees.
+                    continue;
+                if (p_value > 0 && p_value2 > 0)
+                    continue;
+                if (p_value < 0 && p_value2 < 0)
+                    continue;
+                if (any_of(other_ineq.p.free_vars(), [&](unsigned j) { return tabu.contains(j); }))
+                    continue;
 
-            // 
+                for (unsigned i = 0; i < f.degree; ++i)
+                    f.p *= to_pdd(x);
+                for (unsigned i = 0; i < g.degree; ++i)
+                    g.p *= to_pdd(x);
+
+                auto [m1, f_p] = f.p.var_factors();
+                auto [m2, g_p] = g.p.var_factors();
+                unsigned_vector m1m2;
+                dd::pdd::merge(m1, m2, m1m2);
+                SASSERT(m1m2.contains(x));
+                f_p = f_p.mul(m1);
+                g_p = g_p.mul(m2);
+                
+                auto sign_f = cvalue(f_p) < 0;
+                SASSERT(sign_f != cvalue(g_p) < 0);
+                SASSERT(cvalue(f_p) != 0);
+                SASSERT(cvalue(g_p) != 0);
+
+                // m1m2 * f_p + f.q >= 0
+                // m1m2 * g_p + g.q >= 0
+                // 
+                auto ci_a = ci;
+                auto ci_b = other_ci;
+                auto aj = assumption_justification();
+                if (f_p.is_val())
+                    ci_a = multiply(other_ci, f_p.val());
+                else if (sign_f)
+                    ci_a = multiply(other_ci, add_constraint(f_p, lp::lconstraint_kind::LT, aj));
+                else
+                    ci_a = multiply(other_ci, add_constraint(f_p, lp::lconstraint_kind::GT, aj));
+                
+                if (g_p.is_val())
+                    ci_b = multiply(ci, g_p.val());
+                else if (!sign_f)
+                    ci_b = multiply(ci, add_constraint(g_p, lp::lconstraint_kind::LT, aj));
+                else
+                    ci_b = multiply(ci, add_constraint(g_p, lp::lconstraint_kind::GT, aj));
+
+                auto new_ci = add(ci_a, ci_b);
+                init_occurs(new_ci);
+                if (!constraint_is_true(new_ci)) {
+                    auto const &p = m_constraints[ci].p;
+                    auto const &[new_p, new_k, new_j] = m_constraints[new_ci];
+                    uint_set new_tabu(tabu), fv;                 
+                    for (auto v : new_p.free_vars())
+                        fv.insert(v);
+                    for (auto v : p.free_vars())
+                        if (!fv.contains(v))
+                            new_tabu.insert(v);
+                    active.push_back({new_ci, new_tabu});
+                }                
+            }          
         }
     }
 
     lp::lpvar stellensatz::select_variable_to_eliminate(lp::constraint_index ci) {
-        auto &ineq = m_solver.lra().constraints()[ci];
+        auto const& [p, k, j] = m_constraints[ci];
         lpvar best_var = lp::null_lpvar;
-        for (auto cv : ineq.lhs()) {
-            auto v = cv.j();
-            if (best_var == lp::null_lpvar) {
+        for (auto v : p.free_vars())
+            if (best_var > v)
                 best_var = v;
-                continue;
-            }
-            if (is_mon_var(v)) {
-                for (auto w : m_mon2vars[v]) {
-                    if (best_var > v) {
-                        best_var = v;                     
-                    }
-                }
-            }
-            else if (best_var > v) {
-                best_var = v;
-            }
-        }
         return best_var;
     }
 
     unsigned stellensatz::degree_of_var_in_constraint(lpvar var, lp::constraint_index ci) const {
-        auto &ineq = m_solver.lra().constraints()[ci];
-        unsigned max_degree = 0;
-        for (auto cv : ineq.lhs()) {            
-            auto degree = degree_of_var_in_monomial(var, cv.j());
-            if (degree > max_degree)
-                max_degree = degree;
-        }
-        return max_degree;
-    }
-
-    unsigned stellensatz::degree_of_var_in_monomial(lpvar v, lpvar mi) const {
-        unsigned degree = 0;
-        if (is_mon_var(mi)) {
-            for (auto w : m_mon2vars[mi]) 
-                if (v == w)
-                    degree++;            
-        }
-        else if (mi == v)
-            ++degree; 
-        return degree;
+        return m_constraints[ci].p.degree(var);
     }
 
     stellensatz::factorization stellensatz::factor(lpvar v, lp::constraint_index ci) {
-        auto &ineq = m_solver.lra().constraints()[ci];
+        auto const& [p, k, j] = m_constraints[ci];
         auto degree = degree_of_var_in_constraint(v, ci);
-        rational coeff(0);
-        lp::lar_term p, q;
-        for (auto cv : ineq.lhs()) {
-            auto d = degree_of_var_in_monomial(v, cv.j());
-            if (d < degree) 
-                q.add_monomial(cv.coeff(), cv.j());
-            else {
-                auto w = divide(v, degree, cv.j());
-                if (w == lp::null_lpvar)
-                    coeff += cv.coeff();
-                else
-                    p.add_monomial(cv.coeff(), w);               
-            }
-                
-        }
-        return {coeff, degree, p, q};
-    }
-
-    lp::lpvar stellensatz::divide(lpvar v, unsigned degree, lpvar mi) {
-        SASSERT(degree_of_var(v, mi) >= degree);
-        if (is_mon_var(mi)) {
-            auto const &vars = m_mon2vars[mi];
-            if (degree == vars.size())
-                return lp::null_lpvar;
-            svector<lpvar> _vars;
-            for (auto w : vars) {
-                if (w == v && degree > 0) 
-                    --degree;                
-                else
-                    _vars.push_back(w);
-            }
-            return mk_monomial(_vars);
-        }
-        SASSERT(degree == 1 && mi == v);
-        return lp::null_lpvar;
+        dd::pdd lc(pddm), rest(pddm);
+        p.factor(v, degree, lc, rest);
+        return {degree, lc, rest};
     }
 
     //
@@ -282,10 +351,7 @@ namespace nla {
         if (m_constraints_in_conflict.contains(ci))
             return;
         m_constraints_in_conflict.insert(ci);
-        auto just = m_justifications.get(ci);
-        if (just == nullptr)
-            return;
-        auto &b = *just;
+        auto &[p, k, b] = m_constraints[ci];
         if (std::holds_alternative<external_justification>(b)) {
             auto dep = std::get<external_justification>(b);
             m_solver.lra().push_explanation(dep.dep, ex);
@@ -305,8 +371,8 @@ namespace nla {
             explain_constraint(new_lemma, m.ci, ex);
         } 
         else if (std::holds_alternative<assumption_justification>(b)) {
-            auto &con = m_solver.lra().constraints()[ci];
-            new_lemma |= ineq(con.lhs(), negate(con.kind()), con.rhs());
+            auto [t, coeff] = to_term_offset(p);
+            new_lemma |= ineq(t, negate(k), -coeff);
         }
         else if (std::holds_alternative<gcd_justification>(b)) {
             auto& m = std::get<gcd_justification>(b);
@@ -316,297 +382,98 @@ namespace nla {
             UNREACHABLE();
     }
  
-    rational stellensatz::cvalue(lp::lar_term const &t) const {
-        rational r(0);
-        for (auto cv : t)
-            r += c().val(cv.j()) * cv.coeff();
-        return r;
+    rational stellensatz::cvalue(dd::pdd const& p) const {
+        dd::pdd_eval eval;
+        // eval.var2val() = [&](unsigned v) -> void { return c().val(v); };
+        return eval(p);
     }
 
-    rational stellensatz::mvalue(lp::lar_term const &t) const {
-        rational r(0);
-        for (auto cv : t)
-            r += m_values[cv.j()] * cv.coeff();
-        return r;
-    }
-
-    rational stellensatz::cvalue(svector<lpvar> const &prod) const {
-        rational p(1);
-        for (auto v : prod)
-            p *= c().val(v);
-        return p;
-    }
-
-    rational stellensatz::mvalue(svector<lpvar> const &prod) const {
-        rational p(1);
-            for (auto v : prod)
-                p *= m_values[v];
-        return p;
-    }
-
-    void stellensatz::gcd_normalize(vector<std::pair<rational, lpvar>> &t, lp::lconstraint_kind k, rational &rhs) {
-        rational d(1);
-        for (auto &[c,v] : t)
-            d = lcm(d, denominator(c));
-        d = lcm(d, denominator(rhs));
-        if (d != 1) {
-            for (auto &[c, v] : t)
-                c *= d;
-            rhs *= d;
-        }
+    lp::constraint_index stellensatz::gcd_normalize(lp::constraint_index ci) {
+        auto [p, k, j] = m_constraints[ci];
         rational g(0);
-        for (auto &[c, v] : t)
-            g = gcd(g, c);
-        bool is_int = true;
-        for (auto &[c, v] : t)
-            is_int &= m_solver.lra().var_is_int(v);
-        if (!is_int)
-            g = gcd(g, rhs);
-        if (g == 1 || g == 0)
-            return;
+        bool is_int = all_of(p.free_vars(), [&](unsigned v) { return c().lra_solver().var_is_int(v); });
+        for (auto const& [c, is_constant] : p.coefficients()) 
+            if (!is_constant || !is_int)
+                g = gcd(g, c);
+        if (g == 0 || g == 1)
+            return ci;
         switch (k) {
-        case lp::lconstraint_kind::GE:
-            for (auto &[c, v] : t)
-                c /= g;
-            rhs /= g;
-            rhs = ceil(rhs);
-            break;
-        case lp::lconstraint_kind::GT:
-            for (auto &[c, v] : t)
-                c /= g;
-            if (is_int)
-                rhs += 1;
-            rhs /= g;
-            rhs = ceil(rhs);
-            break;
-        case lp::lconstraint_kind::LE:
-            for (auto &[c, v] : t)
-                c /= g;
-            rhs /= g;
-            rhs = floor(rhs);
-            break;
+        case lp::lconstraint_kind::GE: {
+            auto q = p * (1/ g);
+            q += (ceil(q.offset()) - q.offset());
+            return add_constraint(q, k, gcd_justification(ci));
+        }
+        case lp::lconstraint_kind::GT: {
+            auto q = p;
+            if (is_int) {
+                q += rational(1);
+                k = lp::lconstraint_kind::GE;                
+            }
+            q *= (1 / g);
+            q += (ceil(q.offset()) - q.offset());
+            return add_constraint(q, k, gcd_justification(ci));
+        }
         case lp::lconstraint_kind::LT:
-            for (auto &[c, v] : t)
-                c /= g;
-            if (is_int)
-                rhs -= 1;
-            rhs /= g;
-            rhs = floor(rhs);
-            break;
+        case lp::lconstraint_kind::LE:            
+            UNREACHABLE();
         case lp::lconstraint_kind::EQ:
         case lp::lconstraint_kind::NE:
-            if (!divides(g, rhs))
-                break;
-            for (auto &[c, v] : t)
-                c /= g;
-            rhs /= g;
-            break;
-        }    
+            if (!divides(g, p.offset()))
+                return ci;
+            return add_constraint(p * (1/g), k, gcd_justification(ci));
+        default:
+            UNREACHABLE(); 
+            return ci;
+        }        
     }
-
-    lp::constraint_index stellensatz::normalize_ineq(lp::constraint_index ci) {
-        auto const &con = m_solver.lra().constraints()[ci];
-        auto k = con.kind();
-        if (k == lp::lconstraint_kind::LE || k == lp::lconstraint_kind::LT) {
-            auto mjust = multiplication_const_justification({ci, rational(-1)});
-            k = swap_side(k);
-            rational value(0);
-            auto coeffs = con.coeffs();
-            for (auto &[c, v] : coeffs) {
-                c = -c;
-                value += c * m_values[v];
-            }
-            auto ti = m_solver.lra().add_term(coeffs, m_solver.lra().number_of_vars());
-            m_occurs.reserve(m_solver.lra().number_of_vars());
-            m_values.push_back(value);
-            ci = m_solver.lra().add_var_bound(ti, k, -con.rhs());
-        }
-        return ci;
-    }
-
-    lp::constraint_index stellensatz::add_ineq( 
-        justification const& just, lp::lar_term & t, lp::lconstraint_kind k,
-        rational const & rhs_) {
-        auto coeffs = t.coeffs_as_vector();
-        SASSERT(!coeffs.empty());
-        auto add = [&](justification const& just, vector<std::pair<rational, lpvar>> const& coeffs, rational const& rhs) {
-            auto ti = m_solver.lra().add_term(coeffs, m_solver.lra().number_of_vars());
-            m_occurs.reserve(m_solver.lra().number_of_vars());
-            m_values.push_back(mvalue(t));
-            auto new_ci = m_solver.lra().add_var_bound(ti, k, rhs);
-            TRACE(arith, display(tout, just) << "\n");
-            SASSERT(m_values.size() - 1 == ti);
-            add_justification(new_ci, just);
-            return new_ci;
-        };
-        auto new_ci = add(just, coeffs, rhs_);  
-        rational rhs(rhs_);
-        gcd_normalize(coeffs, k, rhs);
-        if (rhs != rhs_) 
-            new_ci = add(gcd_justification(new_ci), coeffs, rhs);        
-        init_occurs(new_ci);
-        return new_ci;
-    }
-
-    lp::constraint_index stellensatz::add_ineq(justification const& just, lpvar j, lp::lconstraint_kind k,
-        rational const& rhs) {
-        auto new_ci = m_solver.lra().add_var_bound(j, k, rhs);
-        add_justification(new_ci, just);
-        init_occurs(new_ci);
-        return new_ci;
-    }
-
-    lp::constraint_index stellensatz::assume(lp::lar_term &t, lp::lconstraint_kind k, rational const &rhs) {
-        return add_ineq(assumption_justification(), t, k, rhs);
+    
+    lp::constraint_index stellensatz::assume(dd::pdd const& p, lp::lconstraint_kind k) {
+        return add_constraint(p, k, assumption_justification());
     }
+    
     // p1 >= lo1, p2 >= lo2 => p1 + p2 >= lo1 + lo2
     lp::constraint_index stellensatz::add(lp::constraint_index left, lp::constraint_index right) {
-        auto const &lcon = m_solver.lra().constraints()[left];
-        auto const &rcon = m_solver.lra().constraints()[right];
-        SASSERT(lcon.kind() == lp::lconstraint_kind::GE || lcon.kind() == lp::lconstraint_kind::GT);
-        SASSERT(rcon.kind() == lp::lconstraint_kind::GE || rcon.kind() == lp::lconstraint_kind::GT);
-        lp::lar_term t;
-        for (auto const &[v, c] : lcon.lhs())
-            t.add_monomial(c, v);
-        for (auto const &[v, c] : rcon.lhs())
-            t.add_monomial(c, v);
-        lp::lconstraint_kind k = join(lcon.kind(), rcon.kind());
-        rational rhs = lcon.rhs() + rcon.rhs();
-        auto just = addition_justification{left, right};
-        return add_ineq(just, t, k, rhs);
+        auto const &[p, k1, j1] = m_constraints[left];
+        auto const &[q, k2, j2] = m_constraints[right];
+        lp::lconstraint_kind k = join(k1, k2);
+        return gcd_normalize(add_constraint(p + q, k, addition_justification{left, right}));
     }
 
     // p >= lo  =>  a * p >= a * lo if a > 0
     lp::constraint_index stellensatz::multiply(lp::constraint_index ci, rational const& mul) {
         SASSERT(mul != 0);
-        auto const &con = m_solver.lra().constraints()[ci];
-        lp::lar_term t;
-        for (auto const &[v, c] : con.lhs())
-            t.add_monomial(c * mul, v);
-        lp::lconstraint_kind k = con.kind();
-        if (mul < 0)
-            k = swap_side(k);
-        rational rhs = con.rhs() * mul;
-        auto just = multiplication_const_justification{ci, mul};
-        return add_ineq(just, t, k, rhs);        
+        auto const& [p, k, j] = m_constraints[ci];
+        return add_constraint(p * mul, mul > 0 ? k : swap_side(k), multiplication_const_justification{ci, mul});
     }
 
     lp::constraint_index stellensatz::multiply(lp::constraint_index left, lp::constraint_index right) {
-        auto const &lcon = m_solver.lra().constraints()[left];
-        auto const &rcon = m_solver.lra().constraints()[right];
-        SASSERT(lcon.kind() == lp::lconstraint_kind::GE || lcon.kind() == lp::lconstraint_kind::GT);
-        SASSERT(rcon.kind() == lp::lconstraint_kind::GE || rcon.kind() == lp::lconstraint_kind::GT);
-        lp::lar_term t;
-        for (auto const &[v1, c1] : lcon.lhs()) {
-            for (auto const &[v2, c2] : rcon.lhs()) {
-                auto j = mk_monomial(expand_monomial({v1, v2}));
-                t.add_monomial(c1 * c2, j);
-            }
-        }
-        if (rcon.rhs() != 0)
-            for (auto const &[v, c] : lcon.lhs()) 
-                t.add_monomial(-c * rcon.rhs(), v);
-        if (lcon.rhs() != 0)
-            for (auto const &[v, c] : rcon.lhs()) 
-                t.add_monomial(-c * lcon.rhs(), v);
-        lp::lconstraint_kind k = join(lcon.kind(), rcon.kind());
-        rational rhs{0};
-        return add_ineq(multiplication_justification{left, right}, t, k, rhs);
+        auto const &[p, k1, j1] = m_constraints[left];
+        auto const &[q, k2, j2] = m_constraints[right];
+        lp::lconstraint_kind k = meet(k1, k2);
+        return add_constraint(p*q, k, multiplication_justification{left, right});
     }
 
     lp::constraint_index stellensatz::multiply_var(lp::constraint_index ci, lpvar x) {
-        auto const &con = m_solver.lra().constraints()[ci];
-        lp::lconstraint_kind k = con.kind();
-        SASSERT(k == lp::lconstriant_kind::EQ);
-        SASSERT(con.rhs() == 0);
-        lp::lar_term t;
-        for (auto const &[v, c] : con.lhs())
-            t.add_monomial(c, mk_monomial(expand_monomial({x, v})));        
-        auto just = multiplication_const_justification{ci, rational(1)}; // multiply var
-        return add_ineq(just, t, k, rational(0));  
+        auto const& [p, k, j] = m_constraints[ci];
+        SASSERT(k == lp::lconstriant_kind::EQ);        
+        SASSERT(p.offset() == 0);
+        return add_constraint(to_pdd(x) * p, k, multiplication_var_justification{ci, x});          
     }
 
     void stellensatz::init_occurs() {
         m_occurs.reset();
         m_occurs.reserve(c().lra_solver().number_of_vars());
-        for (auto ci : m_solver.lra().constraints().indices())
+        for (lp::constraint_index ci = 0; ci < m_constraints.size(); ++ci) 
             init_occurs(ci);        
     }
 
     void stellensatz::init_occurs(lp::constraint_index ci) {
         if (ci == lp::null_ci)
             return;
-        ci = normalize_ineq(ci);
-        auto const &con = m_solver.lra().constraints()[ci];
-        for (auto cv : con.lhs()) {
-            auto v = cv.j();
-            if (is_mon_var(v)) {
-                for (auto w : m_mon2vars[v]) 
-                    m_occurs[w].push_back(ci);                
-            }
-            else 
-                m_occurs[v].push_back(ci);
-        }
-    }
-
-    // Insert a (new) monomial representing product of vars.
-    // if the length of vars is 1 then the new monomial is vars[0]
-    // if the same monomial was previous defined, return the previously defined monomial.
-    // otherwise create a fresh variable representing the monomial.
-    // todo: if _vars is a square, then add bound axiom.
-    lpvar stellensatz::mk_monomial(svector<lpvar> const& vars) {
-        lpvar v;
-        if (vars.empty())
-            return lp::null_lpvar;
-        if (vars.size() == 1)
-            return vars[0];
-        svector<lpvar> _vars(vars);
-        std::sort(_vars.begin(), _vars.end());
-        if (m_vars2mon.find(_vars, v))
-            return v;
-        v = add_var(is_int(vars));
-        m_vars2mon.insert(_vars, v);
-        m_mon2vars.insert(v, _vars);
-        SASSERT(m_values.size() == v);
-        m_values.push_back(cvalue(vars));
-        return v;
-    }
-
-    lpvar stellensatz::mk_monomial(svector<lp::lpvar> const &_vars, lp::lpvar j) {
-        svector<lpvar> vars(_vars);
-        if (is_mon_var(j))
-            vars.append(m_mon2vars[j]);
-        else
-            vars.push_back(j);
-        return mk_monomial(vars);
-    }
-
-    svector<lp::lpvar> stellensatz::expand_monomial(svector<lp::lpvar> const& vars) {
-        svector<lp::lpvar> result;
-        for (auto v : vars) {
-            if (is_mon_var(v)) 
-                result.append(m_mon2vars[v]);            
-            else 
-                result.push_back(v);            
-        }
-        return result;
-    }
-
-    lpvar stellensatz::mk_term(term_offset& t) {
-        if (t.second != 0) {            
-            auto w = add_var(t.second.is_int()); // or reuse the constant 1 that is already there.
-            m_values.push_back(t.second);
-            m_solver.lra().add_var_bound(w, lp::lconstraint_kind::GE, t.second);
-            m_solver.lra().add_var_bound(w, lp::lconstraint_kind::LE, t.second);
-            m_justifications.push_back(nullptr);
-            m_justifications.push_back(nullptr);
-            t.first.add_monomial(rational(1), w);
-            t.second = 0;
-        }
-        auto ti = m_solver.lra().add_term(t.first.coeffs_as_vector(), m_solver.lra().number_of_vars());
-        m_values.push_back(cvalue(t.first));
-        m_occurs.reserve(m_solver.lra().number_of_vars());
-        return ti;
+        auto const &con = m_constraints[ci];
+        for (auto v : con.p.free_vars())
+            m_occurs[v].push_back(ci);
+        
     }
 
     bool stellensatz::is_int(svector<lp::lpvar> const& vars) const {
@@ -622,18 +489,16 @@ namespace nla {
     }
 
     bool stellensatz::constraint_is_true(lp::constraint_index ci) const {
-        auto const& con = m_solver.lra().constraints()[ci];
-        auto lhs = -con.rhs();
-        for (auto const& cv : con.lhs())
-            lhs += cv.coeff() * m_values[cv.j()];
-        switch (con.kind()) {
+        auto const& [p, k, j] = m_constraints[ci];
+        auto lhs = cvalue(p);
+        switch (k) {
         case lp::lconstraint_kind::GT: return lhs > 0;
         case lp::lconstraint_kind::GE: return lhs >= 0;
         case lp::lconstraint_kind::LE: return lhs <= 0;
         case lp::lconstraint_kind::LT: return lhs < 0;
         case lp::lconstraint_kind::EQ: return lhs == 0;
         case lp::lconstraint_kind::NE: return lhs != 0;
-        default: UNREACHABLE(); return false;
+        default: UNREACHABLE();  return false;
         }
     }
 
@@ -644,53 +509,18 @@ namespace nla {
             display_product(out, vars);
             out << "\n";
         }
-        for (auto ci : m_solver.lra().constraints().indices()) {
+        for (unsigned ci = 0; ci < m_constraints.size(); ++ci) {      
             out << "(" << ci << ") ";
             display_constraint(out, ci);
             bool is_true = constraint_is_true(ci);
             out << (is_true ? " [true]" : " [false]") << "\n";
-            auto j = m_justifications.get(ci);
-            if (!j)
-                continue;
             out << "\t<- ";
-            display(out, *j);
+            display(out, m_constraints[ci].j);
             out << "\n";
         }
         return out;
     }
 
-    std::ostream& stellensatz::display(std::ostream& out, justification const& b) const {
-        if (std::holds_alternative<external_justification>(b)) {
-            auto &dep = std::get<external_justification>(b);
-            unsigned_vector cs;
-            c().lra_solver().dep_manager().linearize(dep.dep, cs);
-            for (auto c : cs)
-                out << "[o " << c << "] ";  // constraint index from c().lra_solver.
-        } 
-        else if (std::holds_alternative<multiplication_justification>(b)) {
-            auto &m = std::get<multiplication_justification>(b);
-            out << "mult (" << m.left << ") (" << m.right << ")";
-        }
-        else if (std::holds_alternative<addition_justification>(b)) {
-            auto &m = std::get<addition_justification>(b);
-            out << "(" << m.left << ") (" << m.right << ")";
-        }
-        else if (std::holds_alternative<gcd_justification>(b)) {
-            auto &m = std::get<gcd_justification>(b);
-            out << "gcd(" << m.ci << ")";
-        }
-        else if (std::holds_alternative<assumption_justification>(b)) {
-            out << "assumption";
-        }
-        else if (std::holds_alternative<multiplication_const_justification>(b)) {
-            auto &m = std::get<multiplication_const_justification>(b);
-            out << m.mul << " * (" << m.ci << ")";
-        }
-        else
-            UNREACHABLE();
-        return out;
-    }
-
     std::ostream& stellensatz::display_product(std::ostream& out, svector<lpvar> const& vars) const {
         bool first = true;
         for (auto v : vars) {
@@ -716,32 +546,56 @@ namespace nla {
     }
 
     std::ostream& stellensatz::display_var(std::ostream& out, lpvar j) const {
-        if (is_mon_var(j))
-            display_product(out, m_mon2vars[j]);
-        else
+        if (is_mon_var(j)) {
+//            display_product(out, c().emons()[j]);
+        }
+        else {
             out << "j" << j;
+        }
+
         return out;
     }
 
     std::ostream& stellensatz::display_constraint(std::ostream& out, lp::constraint_index ci) const {
-        auto const& con = m_solver.lra().constraints()[ci];
-        return display_constraint(out, con.lhs(), con.kind(), con.rhs());
+        return display_constraint(out, m_constraints[ci]);
     }
 
-    std::ostream& stellensatz::display_constraint(std::ostream& out, 
-        lp::lar_term const& lhs,
-        lp::lconstraint_kind k, rational const& rhs) const {
-        bool first = true;
-        for (auto cv : lhs) {
-            auto v = cv.j();
-            c().display_coeff(out, first, cv.coeff());
-            first = false;
-            if (is_mon_var(v))
-                display_product(out, m_mon2vars[v]);
-            else
-                out << "j" << v;
+    std::ostream& stellensatz::display_constraint(std::ostream& out, constraint const &c) const {
+        auto const &[p, k, j] = c;
+        p.display(out);
+        return out << " " << k << " 0";
+    }
+
+    std::ostream &stellensatz::display(std::ostream &out, justification const &j) const {
+        if (std::holds_alternative<external_justification>(j)) {
+            auto dep = std::get<external_justification>(j).dep;
+            unsigned_vector cs;
+            c().lra_solver().dep_manager().linearize(dep, cs);
+            for (auto c : cs)
+                out << "[o " << c << "] ";  // constraint index from c().lra_solver.
         }
-        return out << " " << k << " " << rhs;
+        else if (std::holds_alternative<multiplication_justification>(j)) {
+            auto m = std::get<multiplication_justification>(j);
+            out << "(" << m.left << ") * (" << m.right << ")";
+        }
+        else if (std::holds_alternative<addition_justification>(j)) {
+            auto m = std::get<addition_justification>(j);
+            out << "(" << m.left << ") + (" << m.right << ")";
+        }
+        else if (std::holds_alternative<multiplication_const_justification>(j)) {
+            auto m = std::get<multiplication_const_justification>(j);
+            out << m.mul << " * (" << m.ci << ")";
+        }
+        else if (std::holds_alternative<assumption_justification>(j)) {
+            out << "assumption";
+        }
+        else if (std::holds_alternative<gcd_justification>(j)) {
+            auto m = std::get<gcd_justification>(j);
+            out << " gcd (" << m.ci << ")";
+        }
+        else
+            UNREACHABLE();
+        return out;
     }
 
     std::ostream& stellensatz::display_lemma(std::ostream& out, lp::explanation const& ex,
@@ -757,32 +611,31 @@ namespace nla {
             m_constraints_in_conflict.insert(ci);
             out << "(" << ci << ") ";
             display_constraint(out, ci) << " ";
-            auto j = m_justifications.get(ci);
-            if (!j)
-                continue;
-            display(out, *j) << "\n";
-            if (std::holds_alternative<multiplication_justification>(*j)) {
-                auto m = std::get<multiplication_justification>(*j);
+            auto const& j = m_constraints[ci].j;
+
+            display(out, j) << "\n";
+            if (std::holds_alternative<multiplication_justification>(j)) {
+                auto m = std::get<multiplication_justification>(j);
                 todo.push_back(m.left);
                 todo.push_back(m.right);
             }
-            else if (std::holds_alternative<addition_justification>(*j)) {
-                auto m = std::get<addition_justification>(*j);
+            else if (std::holds_alternative<addition_justification>(j)) {
+                auto m = std::get<addition_justification>(j);
                 todo.push_back(m.left);
                 todo.push_back(m.right);
             }
-            else if (std::holds_alternative<multiplication_const_justification>(*j)) {
-                auto m = std::get<multiplication_const_justification>(*j);
+            else if (std::holds_alternative<multiplication_const_justification>(j)) {
+                auto m = std::get<multiplication_const_justification>(j);
                 todo.push_back(m.ci);
             }
-            else if (std::holds_alternative<external_justification>(*j)) {
+            else if (std::holds_alternative<external_justification>(j)) {
                 // do nothing
             }
-            else if (std::holds_alternative<assumption_justification>(*j)) {
+            else if (std::holds_alternative<assumption_justification>(j)) {
                 // do nothing
             }
-            else if (std::holds_alternative<gcd_justification>(*j)) {
-                auto m = std::get<gcd_justification>(*j);
+            else if (std::holds_alternative<gcd_justification>(j)) {
+                auto m = std::get<gcd_justification>(j);
                 todo.push_back(m.ci);
             }
             else 
@@ -860,18 +713,22 @@ namespace nla {
             auto const &value = values[i];
             bool is_int = lra_solver->var_is_int(i); 
             SASSERT(!is_int || value.is_int());
-            if (s.is_mon_var(i)) 
-                s.m_values[i] = s.mvalue(s.m_mon2vars[i]);             
-            else
-                s.m_values[i] = value;
+            if (!s.is_mon_var(i))
+                continue;
+            rational val(1);
+            for (auto w : s.c().emons()[i])
+                val *= values[w];
+            values[i] = val;
         }
         auto indices = lra_solver->constraints().indices();
         bool satisfies_products = all_of(indices, [&](auto ci) {
+            NOT_IMPLEMENTED_YET();
+            // todo: wrong, needs to be at level of lra constraint evaluation
             return s.constraint_is_true(ci);});
         if (satisfies_products) {
             TRACE(arith, tout << "found satisfying model\n");
             for (auto v : s.m_coi.vars()) 
-                s.c().lra_solver().set_column_value(v, lp::impq(s.m_values[v], rational(0)));            
+                s.c().lra_solver().set_column_value(v, lp::impq(values[v], rational(0)));            
         }
         return satisfies_products;
     }
diff --git a/src/math/lp/nla_stellensatz.h b/src/math/lp/nla_stellensatz.h
index d58a3de24..4449fe7ba 100644
--- a/src/math/lp/nla_stellensatz.h
+++ b/src/math/lp/nla_stellensatz.h
@@ -5,9 +5,9 @@
 #pragma once
 
 #include "util/scoped_ptr_vector.h"
+#include "math/dd/dd_pdd.h"
 #include "math/lp/nla_coi.h"
 #include "math/lp/int_solver.h"
-#include "math/polynomial/polynomial.h"
 #include <variant>
 
 namespace nla {
@@ -26,7 +26,6 @@ namespace nla {
             lbool solve_lia(); 
             bool update_values();
             vector<std::pair<lpvar, rational>> m_to_refine;
-            void saturate_basic_linearize();
         public:  
             solver(stellensatz &s) : s(s) {};                
             void init();
@@ -39,12 +38,11 @@ namespace nla {
 
         solver m_solver;
 
-        // factor t into coeff*x^degree*p + q, such that degree(x, q) < degree, 
+
+        // factor t into x^degree*p + q, such that degree(x, q) < degree, 
         struct factorization {
-            rational coeff;
             unsigned degree;
-            lp::lar_term p;
-            lp::lar_term q;
+            dd::pdd p, q;
         };
 
         struct external_justification {
@@ -60,6 +58,10 @@ namespace nla {
             lp::constraint_index ci;
             rational mul;
         };
+        struct multiplication_var_justification {
+            lp::constraint_index ci;
+            lpvar v;
+        };
         struct multiplication_justification {
             lp::constraint_index left, right;
         };
@@ -73,16 +75,32 @@ namespace nla {
             gcd_justification,
             addition_justification, 
             multiplication_const_justification, 
+            multiplication_var_justification, 
             multiplication_justification
         >;
 
+        using term_offset = std::pair<lp::lar_term, rational>;
+
+        struct constraint {
+            dd::pdd p;
+            lp::lconstraint_kind k;
+            justification j;
+        };
+
         coi m_coi;
-        scoped_ptr_vector<justification> m_justifications;
-        void add_justification(lp::constraint_index ci, justification const &j) {
-            VERIFY(ci == m_justifications.size());
-            m_justifications.push_back(alloc(justification, j));
-        }
-        vector<rational> m_values;
+        dd::pdd_manager pddm;
+        vector<constraint> m_constraints;
+
+
+        dd::pdd to_pdd(lpvar v);
+        lpvar mk_monomial(svector<lpvar> const &vars);
+        void init_monomial(unsigned mon_var); 
+        term_offset to_term_offset(dd::pdd const &p);
+
+        lp::constraint_index add_constraint(dd::pdd p, lp::lconstraint_kind k, justification j);
+        
+        lp::constraint_index add_var_bound(lp::lpvar v, lp::lconstraint_kind k, rational const &rhs, justification j);
+
         struct eq {
             bool operator()(unsigned_vector const& a, unsigned_vector const& b) const {
                 return a == b;
@@ -92,55 +110,30 @@ namespace nla {
         u_map<unsigned_vector> m_mon2vars;
         bool is_mon_var(lpvar v) const { return m_mon2vars.contains(v); }
         
-        unsigned m_max_monomial_degree = 0;
-
         vector<svector<lp::constraint_index>> m_occurs; // map from variable to constraints they occur. 
 
-        bool is_new_inequality(vector<std::pair<rational, lpvar>> lhs, lp::lconstraint_kind k, rational const &rhs);
-
-
         // initialization
         void init_solver();
         void init_vars();
-        void init_monomial(unsigned mon_var);
         void init_occurs();
         void init_occurs(lp::constraint_index ci);
 
         void eliminate_variables();
         lp::lpvar select_variable_to_eliminate(lp::constraint_index ci);
         unsigned degree_of_var_in_constraint(lpvar v, lp::constraint_index ci) const;
-        unsigned degree_of_var_in_monomial(lpvar v, lpvar mi) const;
-        factorization factor(lpvar v, lp::constraint_index ci);
-        lp::lpvar divide(lpvar v, unsigned degree, lpvar mi);
+        factorization factor(lpvar v, lp::constraint_index ci);    
 
         bool constraint_is_true(lp::constraint_index ci) const;
 
-        // additional variables and monomials and constraints
-        using term_offset = std::pair<lp::lar_term, rational>;  // term and its offset
-
-        svector<lp::lpvar> expand_monomial(svector<lp::lpvar> const & vars);
-        lpvar mk_monomial(svector<lp::lpvar> const& vars);
-        lpvar mk_monomial(svector<lp::lpvar> const &vars, lp::lpvar j);
-        lpvar mk_term(term_offset &t);
-
-        void gcd_normalize(vector<std::pair<rational, lpvar>> &t, lp::lconstraint_kind k, rational &rhs);
-        lp::constraint_index add_ineq(justification const& just, lp::lar_term &t, lp::lconstraint_kind k, rational const &rhs);
-
-        lp::constraint_index add_ineq(justification const &just, lpvar j, lp::lconstraint_kind k,
-                                      rational const &rhs);
-
-        lp::constraint_index assume(lp::lar_term &t, lp::lconstraint_kind k, rational const &rhs);
-        lp::constraint_index normalize_ineq(lp::constraint_index ci);
+        lp::constraint_index gcd_normalize(lp::constraint_index ci);
+        lp::constraint_index assume(dd::pdd const& p, lp::lconstraint_kind k);
         lp::constraint_index add(lp::constraint_index left, lp::constraint_index right);
         lp::constraint_index multiply(lp::constraint_index ci, rational const &mul);
         lp::constraint_index multiply(lp::constraint_index left, lp::constraint_index right);
         lp::constraint_index multiply_var(lp::constraint_index ci, lpvar x);
 
         bool is_int(svector<lp::lpvar> const& vars) const;
-        rational cvalue(lp::lar_term const &t) const;
-        rational mvalue(lp::lar_term const &t) const;
-        rational cvalue(svector<lpvar> const &prod) const;
-        rational mvalue(svector<lpvar> const &prod) const;
+        rational cvalue(dd::pdd const& p) const;
         lpvar add_var(bool is_int);
 
         // lemmas
@@ -162,12 +155,11 @@ namespace nla {
         std::ostream& display(std::ostream& out) const;
         std::ostream& display_product(std::ostream& out, svector<lpvar> const& vars) const;
         std::ostream& display_constraint(std::ostream& out, lp::constraint_index ci) const;
-        std::ostream& display_constraint(std::ostream& out, lp::lar_term const& lhs,
-                                         lp::lconstraint_kind k, rational const& rhs) const;
-        std::ostream& display(std::ostream &out, term_offset const &t) const;
-        std::ostream& display(std::ostream &out, justification const &bounds) const;
+        std::ostream& display_constraint(std::ostream& out, constraint const& c) const;
+        std::ostream &display(std::ostream &out, justification const &j) const;
         std::ostream &display_var(std::ostream &out, lpvar j) const;
         std::ostream &display_lemma(std::ostream &out, lp::explanation const &ex, vector<ineq> const &ineqs);
+        std::ostream &display(std::ostream &out, term_offset const &t) const;
 
     public:
         stellensatz(core* core);