diff --git a/src/sat/smt/polysat/CMakeLists.txt b/src/sat/smt/polysat/CMakeLists.txt
index 80168e713..597f21e2e 100644
--- a/src/sat/smt/polysat/CMakeLists.txt
+++ b/src/sat/smt/polysat/CMakeLists.txt
@@ -6,6 +6,7 @@ z3_add_component(polysat
     fixed_bits.cpp
     forbidden_intervals.cpp
     inequality.cpp
+    monomials.cpp
     op_constraint.cpp
     refine.cpp
     saturation.cpp
diff --git a/src/sat/smt/polysat/constraints.cpp b/src/sat/smt/polysat/constraints.cpp
index 5bf8df60b..90b241a08 100644
--- a/src/sat/smt/polysat/constraints.cpp
+++ b/src/sat/smt/polysat/constraints.cpp
@@ -25,7 +25,9 @@ namespace polysat {
         pdd lhs = p, rhs = q;
         bool is_positive = true;
         ule_constraint::simplify(is_positive, lhs, rhs);
-        auto* cnstr = alloc(ule_constraint, lhs, rhs);
+        pdd unfold_lhs = c.ms().subst(lhs);
+        pdd unfold_rhs = c.ms().subst(rhs);
+        auto* cnstr = alloc(ule_constraint, lhs, rhs, unfold_lhs, unfold_rhs);
         c.trail().push(new_obj_trail(cnstr));
         auto sc = signed_constraint(ckind_t::ule_t, cnstr);
         return is_positive ? sc : ~sc;
@@ -61,6 +63,25 @@ namespace polysat {
         return signed_constraint(ckind_t::op_t, cnstr);
     }
 
+    // parity p >= k if low order k bits of p are 0
+    signed_constraint constraints::parity_at_least(pdd const& p, unsigned k) {
+        if (k == 0)
+            return uge(p, 0);
+        unsigned N = p.manager().power_of_2();
+        // parity(p) >= k
+        // <=> p * 2^(N - k) == 0
+        if (k > N) 
+            // parity(p) > N is never true
+            return ~eq(p.manager().zero());        
+        else if (k == 0) 
+            // parity(p) >= 0 is a tautology
+            return eq(p.manager().zero());        
+        else if (k == N)
+            return eq(p);
+        else
+            return eq(p * rational::power_of_two(N - k));
+    }
+
     bool signed_constraint::is_eq(pvar& v, rational& val) {
         if (m_sign)
             return false;
@@ -87,6 +108,11 @@ namespace polysat {
         return m_sign ? ~r : r;
     }
 
+    lbool signed_constraint::eval_unfold(assignment& a) const {
+        lbool r = m_constraint->eval_unfold(a);
+        return m_sign ? ~r : r;
+    }
+
     std::ostream& signed_constraint::display(std::ostream& out) const {
         if (m_sign) out << "~";
         return out << *m_constraint;        
diff --git a/src/sat/smt/polysat/constraints.h b/src/sat/smt/polysat/constraints.h
index 13a077eee..0851f1505 100644
--- a/src/sat/smt/polysat/constraints.h
+++ b/src/sat/smt/polysat/constraints.h
@@ -40,10 +40,12 @@ namespace polysat {
         bool contains_var(pvar v) const { return m_vars.contains(v); }
         unsigned num_watch() const { return m_num_watch;  }
         void set_num_watch(unsigned n) { SASSERT(n <= 2);  m_num_watch = n; }
+        virtual unsigned_vector const& unfold_vars() const { return m_vars; }
         virtual std::ostream& display(std::ostream& out, lbool status) const = 0;
         virtual std::ostream& display(std::ostream& out) const = 0;
         virtual lbool eval() const = 0;
         virtual lbool eval(assignment const& a) const = 0;
+        virtual lbool eval_unfold(assignment const& a) const { return eval(a); }
         virtual void activate(core& c, bool sign, dependency const& d) = 0;
         virtual void propagate(core& c, lbool value, dependency const& d) = 0;
         virtual bool is_linear() const { return false; }
@@ -65,6 +67,7 @@ namespace polysat {
         bool is_negative() const { return m_sign; }
         unsigned_vector& vars() { return m_constraint->vars(); }
         unsigned_vector const& vars() const { return m_constraint->vars(); }
+        unsigned_vector const& unfold_vars() const { return m_constraint->unfold_vars(); }
         unsigned var(unsigned idx) const { return m_constraint->var(idx); }
         bool contains_var(pvar v) const { return m_constraint->contains_var(v); }
         unsigned num_watch() const { return m_constraint->num_watch(); }
@@ -77,6 +80,7 @@ namespace polysat {
         bool is_currently_false(core& c) const;
         bool is_linear() const { return m_constraint->is_linear(); }
         lbool eval(assignment& a) const;
+        lbool eval_unfold(assignment& a) const;
         lbool eval() const { return m_sign ? ~m_constraint->eval() : m_constraint->eval();}
         ckind_t op() const { return m_op; }
         bool is_ule() const { return m_op == ule_t; }
@@ -132,6 +136,14 @@ namespace polysat {
         signed_constraint ult(pdd const& p, int q) { return ult(p, rational(q)); }
         signed_constraint ult(pdd const& p, unsigned q) { return ult(p, rational(q)); }
 
+        signed_constraint ugt(pdd const& p, pdd const& q) { return ult(q, p); }
+        signed_constraint ugt(pdd const& p, rational const& q) { return ugt(p, p.manager().mk_val(q)); }
+        signed_constraint ugt(rational const& p, pdd const& q) { return ugt(q.manager().mk_val(p), q); }
+        signed_constraint ugt(int             p, pdd const& q) { return ugt(rational(p), q); }
+        signed_constraint ugt(unsigned        p, pdd const& q) { return ugt(rational(p), q); }
+        signed_constraint ugt(pdd const& p, int q) { return ugt(p, rational(q)); }
+        signed_constraint ugt(pdd const& p, unsigned q) { return ugt(p, rational(q)); }
+
         signed_constraint slt(pdd const& p, rational const& q) { return slt(p, p.manager().mk_val(q)); }
         signed_constraint slt(rational const& p, pdd const& q) { return slt(q.manager().mk_val(p), q); }
         signed_constraint slt(pdd const& p, int             q) { return slt(p, rational(q)); }
@@ -156,6 +168,8 @@ namespace polysat {
         signed_constraint umul_ovfl(int             p, pdd const& q) { return umul_ovfl(rational(p), q); }
         signed_constraint umul_ovfl(unsigned        p, pdd const& q) { return umul_ovfl(rational(p), q); }
 
+        signed_constraint parity_at_least(pdd const& p, unsigned k);
+
         signed_constraint lshr(pdd const& a, pdd const& b, pdd const& r);
         signed_constraint ashr(pdd const& a, pdd const& b, pdd const& r);
         signed_constraint shl(pdd const& a, pdd const& b, pdd const& r);
diff --git a/src/sat/smt/polysat/core.cpp b/src/sat/smt/polysat/core.cpp
index d557fa789..0fb1e1635 100644
--- a/src/sat/smt/polysat/core.cpp
+++ b/src/sat/smt/polysat/core.cpp
@@ -88,6 +88,7 @@ namespace polysat {
         m_viable(*this),
         m_constraints(*this),
         m_assignment(*this),
+        m_monomials(*this),
         m_var_queue(m_activity)
     {}
 
@@ -200,12 +201,21 @@ namespace polysat {
     sat::check_result core::final_check() {
         unsigned n = 0;
         constraint_id conflict_idx = { UINT_MAX };
-        // verbose_stream() << "final-check\n";
+
+        switch (m_monomials.refine()) {
+        case l_true:
+            break;
+        case l_false:
+            return sat::check_result::CR_CONTINUE;
+        case l_undef:
+            break;
+        }
+        
         for (auto idx : m_prop_queue) {
             auto [sc, d, value] = m_constraint_index[idx.id];
             SASSERT(value != l_undef);
             // verbose_stream() << "constraint " << (value == l_true ? sc : ~sc) << "\n";
-            lbool eval_value = eval(sc);
+            lbool eval_value = eval_unfold(sc);
             CTRACE("bv", eval_value == l_undef, sc.display(tout << "eval: ") << " evaluates to " << eval_value << "\n"; display(tout););
             SASSERT(eval_value != l_undef);
             if (eval_value == value)
@@ -224,11 +234,21 @@ namespace polysat {
         // If all constraints evaluate to true, we are done.
         if (conflict_idx.is_null())
             return sat::check_result::CR_DONE;
+
+        switch (m_monomials.bit_blast()) {
+        case l_true:
+            break;
+        case l_false:
+            return sat::check_result::CR_CONTINUE;
+        case l_undef:
+            break;
+        }
         
         // If no saturation propagation was possible, explain the conflict using the variable assignment.
         m_unsat_core = explain_eval(get_constraint(conflict_idx));
         m_unsat_core.push_back(get_dependency(conflict_idx));
         s.set_conflict(m_unsat_core, "polysat-bail-out-conflict");
+        decay_activity();
         return sat::check_result::CR_CONTINUE;
     }
 
@@ -239,7 +259,7 @@ namespace polysat {
         svector<pvar> result;
         auto [sc, d, value] = m_constraint_index[idx.id];
         unsigned lvl = s.level(d);
-        for (auto v : sc.vars()) {
+        for (auto v : sc.unfold_vars()) {
             if (!is_assigned(v))
                 continue;
             auto new_level = s.level(m_justification[v]);
@@ -265,8 +285,8 @@ namespace polysat {
 
     void core::viable_conflict(pvar v) {
         TRACE("bv", tout << "viable-conflict v" << v << "\n");
-        m_var_queue.activity_increased_eh(v);
         s.set_conflict(m_viable.explain(), "viable-conflict");
+        decay_activity();
     }
 
     void core::propagate_assignment(constraint_id idx) { 
@@ -278,6 +298,7 @@ namespace polysat {
         TRACE("bv", tout << "propagate " << sc << " using " << dep << " := " << value << "\n");
         if (sc.is_always_false()) {
             s.set_conflict({dep}, "infeasible assignment");
+            decay_activity();
             return;
         }
         rational var_value;
@@ -398,6 +419,7 @@ namespace polysat {
             m_unsat_core = explain_eval(sc);
             m_unsat_core.push_back(m_constraint_index[id.id].d);
             s.set_conflict(m_unsat_core, "polysat-constraint-core");
+            decay_activity();
         }                   
     }
 
@@ -436,9 +458,12 @@ namespace polysat {
 
     dependency_vector core::explain_eval(signed_constraint const& sc) {
         dependency_vector deps;
-        for (auto v : sc.vars()) 
-            if (is_assigned(v))
+        for (auto v : sc.vars()) {
+            if (is_assigned(v)) {
+                inc_activity(v);
                 deps.push_back(m_justification[v]);
+            }
+        }
         return deps;
     }
 
@@ -446,6 +471,10 @@ namespace polysat {
         return sc.eval(m_assignment);
     }
 
+    lbool core::eval_unfold(signed_constraint const& sc) {
+        return sc.eval_unfold(m_assignment);
+    }
+
     pdd core::subst(pdd const& p) { 
         return m_assignment.apply_to(p);
     }
@@ -462,6 +491,10 @@ namespace polysat {
         s.add_axiom(name, begin, end, is_redundant);
     }
 
+    void core::add_axiom(char const* name, constraint_or_dependency_vector const& cs, bool is_redundant) {
+        s.add_axiom(name, cs.begin(), cs.end(), is_redundant);
+    }
+
     std::ostream& core::display(std::ostream& out) const {
         if (m_constraint_index.empty() && m_vars.empty())
             return out;
@@ -504,4 +537,27 @@ namespace polysat {
         return get_constraint(id).eval(m_assignment);
     }
 
+    lbool core::eval_unfold(constraint_id id) {
+        return get_constraint(id).eval_unfold(m_assignment);
+    }
+
+    void core::inc_activity(pvar v) {
+        unsigned& act = m_activity[v].second;
+        act += m_activity_inc;
+        m_var_queue.activity_increased_eh(v);
+        if (act > (1 << 24))
+            rescale_activity();        
+    }
+
+    void core::rescale_activity() {
+        for (auto& act : m_activity) 
+            act.second >>= 14;        
+        m_activity_inc >>= 14;
+    }
+
+    void core::decay_activity() {
+        m_activity_inc *= 110;
+        m_activity_inc /= 100;
+    }
+
 }
diff --git a/src/sat/smt/polysat/core.h b/src/sat/smt/polysat/core.h
index d14748aa9..8f58aa118 100644
--- a/src/sat/smt/polysat/core.h
+++ b/src/sat/smt/polysat/core.h
@@ -25,6 +25,7 @@ Author:
 #include "sat/smt/polysat/constraints.h"
 #include "sat/smt/polysat/viable.h"
 #include "sat/smt/polysat/assignment.h"
+#include "sat/smt/polysat/monomials.h"
 
 namespace polysat {
 
@@ -53,6 +54,7 @@ namespace polysat {
         viable m_viable;
         constraints m_constraints;
         assignment m_assignment;
+        monomials  m_monomials;
         unsigned m_qhead = 0;
         constraint_id_vector m_prop_queue;
         svector<constraint_info> m_constraint_index;  // index of constraints
@@ -92,6 +94,11 @@ namespace polysat {
 
         void add_axiom(signed_constraint sc);
 
+        unsigned m_activity_inc = 128;
+        void inc_activity(pvar v);
+        void rescale_activity();
+        void decay_activity();
+
     public:
         core(solver_interface& s);
 
@@ -123,6 +130,7 @@ namespace polysat {
         void band(pdd const& a, pdd const& b, pdd const& r) { add_axiom(m_constraints.band(a, b, r)); }
 
         pdd bnot(pdd p) { return -p - 1; }
+        pdd mul(unsigned n, pdd const* args) { return m_monomials.mk(n, args); }
 
 
         /*
@@ -132,12 +140,14 @@ namespace polysat {
         */
         void add_axiom(char const* name, constraint_or_dependency_list const& cs, bool is_redundant);
         void add_axiom(char const* name, constraint_or_dependency const* begin, constraint_or_dependency const* end, bool is_redundant);
+        void add_axiom(char const* name, constraint_or_dependency_vector const& cs, bool is_redundant);
         
         pvar add_var(unsigned sz);
         pdd var(pvar p) { return m_vars[p]; }
         unsigned size(pvar v) const { return m_vars[v].power_of_2(); }
 
         constraints& cs() { return m_constraints; }
+        monomials& ms() { return m_monomials; }
         trail_stack& trail();
 
         std::ostream& display(std::ostream& out) const;
@@ -158,8 +168,10 @@ namespace polysat {
         dependency get_dependency(constraint_id idx) const;
         // dependency_vector get_dependencies(constraint_id_vector const& ids) const;
         lbool eval(constraint_id id);
+        lbool eval_unfold(constraint_id id);
         dependency propagate(signed_constraint const& sc, dependency_vector const& deps) { return s.propagate(sc, deps, nullptr); }
         lbool eval(signed_constraint const& sc);
+        lbool eval_unfold(signed_constraint const& sc);
         dependency_vector explain_eval(signed_constraint const& sc);
         bool inconsistent() const;
 
@@ -167,6 +179,8 @@ namespace polysat {
         * Constraints
         */
         assignment& get_assignment() { return m_assignment; }
+
+        random_gen& rand() { return m_rand; }
     };
 
 }
diff --git a/src/sat/smt/polysat/inequality.cpp b/src/sat/smt/polysat/inequality.cpp
index 9a6ba486b..177adec37 100644
--- a/src/sat/smt/polysat/inequality.cpp
+++ b/src/sat/smt/polysat/inequality.cpp
@@ -25,9 +25,9 @@ namespace polysat {
         auto src = c.get_constraint(id);
         ule_constraint const& ule = src.to_ule();
         if (src.is_positive())
-            return inequality(c, id, ule.lhs(), ule.rhs(), src);
+            return inequality(c, id, ule.unfold_lhs(), ule.unfold_rhs(), src);
         else
-            return inequality(c, id, ule.rhs(), ule.lhs(), src);
+            return inequality(c, id, ule.unfold_rhs(), ule.unfold_lhs(), src);
     }
 
 
@@ -36,108 +36,11 @@ namespace polysat {
     }
 
     bool inequality::is_l_v(pdd const& v, signed_constraint const& sc) {
-        return sc.is_ule() && v == (sc.sign() ? sc.to_ule().rhs() : sc.to_ule().lhs());
+        return sc.is_ule() && v == (sc.sign() ? sc.to_ule().unfold_rhs() : sc.to_ule().unfold_lhs());
     }
 
     bool inequality::is_g_v(pdd const& v, signed_constraint const& sc) {
-        return sc.is_ule() && v == (sc.sign() ? sc.to_ule().lhs() : sc.to_ule().rhs());
+        return sc.is_ule() && v == (sc.sign() ? sc.to_ule().unfold_lhs() : sc.to_ule().unfold_rhs());
     }
 
-#if 0
-
-
-    bool saturation::verify_Y_l_AxB(pvar x, inequality const& i, pdd const& y, pdd const& a, pdd& b) {
-        return i.lhs() == y && i.rhs() == a * c.var(x) + b;
-    }
-
-
-    /**
-     * Match [x] a*x + b <= y, val(y) = 0
-     */
-    bool saturation::is_AxB_eq_0(pvar x, inequality const& i, pdd& a, pdd& b, pdd& y) {
-        y.reset(i.rhs().manager());
-        y = i.rhs();
-        rational y_val;
-        if (!c.try_eval(y, y_val) || y_val != 0)
-            return false;
-        return i.lhs().degree(x) == 1 && (i.lhs().factor(x, 1, a, b), true);
-    }
-
-    bool saturation::verify_AxB_eq_0(pvar x, inequality const& i, pdd const& a, pdd const& b, pdd const& y) {
-        return y.is_val() && y.val() == 0 && i.rhs() == y && i.lhs() == a * c.var(x) + b;
-    }
-
-    bool saturation::is_AxB_diseq_0(pvar x, inequality const& i, pdd& a, pdd& b, pdd& y) {
-        if (!i.is_strict())
-            return false;
-        y.reset(i.lhs().manager());
-        y = i.lhs();
-        if (i.rhs().is_val() && i.rhs().val() == 1)
-            return false;
-        rational y_val;
-        if (!c.try_eval(y, y_val) || y_val != 0)
-            return false;
-        a.reset(i.rhs().manager());
-        b.reset(i.rhs().manager());
-        return i.rhs().degree(x) == 1 && (i.rhs().factor(x, 1, a, b), true);
-    }
-
-    /**
-     * Match [coeff*x] coeff*x*Y where x is a variable
-     */
-    bool saturation::is_coeffxY(pdd const& x, pdd const& p, pdd& y) {
-        pdd xy = x.manager().zero();
-        return x.is_unary() && p.try_div(x.hi().val(), xy) && xy.factor(x.var(), 1, y);
-    }
-
-    /**
- * Match [v] v*x <= z*x with x a variable
- */
-    bool saturation::is_Xy_l_XZ(pvar v, inequality const& i, pdd& x, pdd& z) {
-        return is_xY(v, i.lhs(), x) && is_coeffxY(x, i.rhs(), z);
-    }
-
-    bool saturation::verify_Xy_l_XZ(pvar v, inequality const& i, pdd const& x, pdd const& z) {
-        return i.lhs() == c.var(v) * x && i.rhs() == z * x;
-    }
-
-
-    /**
-     * Determine whether values of x * y is non-overflowing.
-     */
-    bool saturation::is_non_overflow(pdd const& x, pdd const& y) {
-        rational x_val, y_val;
-        rational bound = x.manager().two_to_N();
-        return c.try_eval(x, x_val) && c.try_eval(y, y_val) && x_val * y_val < bound;
-    }
-
-
-    /**
-     * Match [z] yx <= zx with x a variable
-     */
-    bool saturation::is_YX_l_zX(pvar z, inequality const& c, pdd& x, pdd& y) {
-        return is_xY(z, c.rhs(), x) && is_coeffxY(x, c.lhs(), y);
-    }
-
-    bool saturation::verify_YX_l_zX(pvar z, inequality const& i, pdd const& x, pdd const& y) {
-        return i.lhs() == y * x && i.rhs() == c.var(z) * x;
-    }
-
-    /**
-     * Match [x] xY <= xZ
-     */
-    bool saturation::is_xY_l_xZ(pvar x, inequality const& c, pdd& y, pdd& z) {
-        return is_xY(x, c.lhs(), y) && is_xY(x, c.rhs(), z);
-    }
-
-    /**
-     * Match xy = x * Y
-     */
-    bool saturation::is_xY(pvar x, pdd const& xy, pdd& y) {
-        return xy.degree(x) == 1 && xy.factor(x, 1, y);
-    }
-
-#endif
-
-
 }
diff --git a/src/sat/smt/polysat/monomials.cpp b/src/sat/smt/polysat/monomials.cpp
new file mode 100644
index 000000000..6fa890116
--- /dev/null
+++ b/src/sat/smt/polysat/monomials.cpp
@@ -0,0 +1,338 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    Polysat monomials
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2021-03-19
+    Jakob Rath 2021-04-6
+
+--*/
+
+#include "sat/smt/polysat/monomials.h"
+#include "sat/smt/polysat/core.h"
+#include "sat/smt/polysat/number.h"
+
+namespace polysat {
+
+    monomials::monomials(core& c):
+        c(c), 
+        C(c.cs()),
+        m_hash(*this),
+        m_eq(*this),
+        m_table(DEFAULT_HASHTABLE_INITIAL_CAPACITY, m_hash, m_eq) 
+    {}
+        
+    pdd monomials::mk(unsigned n, pdd const* args) {
+        SASSERT(n > 0);
+        if (n == 1)
+            return args[0];
+        if (n == 2 && args[0].is_val())
+            return args[0]*args[1];
+        if (n == 2 && args[1].is_val())
+            return args[0]*args[1];
+        auto& m = args[0].manager();
+        unsigned sz = m.power_of_2();
+        m_tmp.reset();
+        pdd offset = c.value(rational(1), sz);
+        for (unsigned i = 0; i < n; ++i) {
+            pdd const& p = args[i];
+            if (p.is_val())
+                offset *= p;
+            else 
+                m_tmp.push_back(p);
+        }
+        if (m_tmp.empty())
+            return offset;
+        if (m_tmp.size() == 1)
+            return offset * m_tmp[0];
+
+        std::stable_sort(m_tmp.begin(), m_tmp.end(), 
+                         [&](pdd const& a, pdd const& b) { return a.index() < b.index(); });
+
+        unsigned index = m_monomials.size();
+
+        m_monomials.push_back({m_tmp, offset, offset, {}, rational(0) });
+        unsigned j;
+        if (m_table.find(index, j)) {
+            m_monomials.pop_back();
+            return offset * m_monomials[j].var;
+        }
+
+        struct del_monomial : public trail {
+            monomials& m;
+            del_monomial(monomials& m):m(m) {}
+            void undo() override {
+                unsigned index = m.m_monomials.size()-1;
+                m.m_table.erase(index);
+                m.m_monomials.pop_back();
+            }
+        };
+
+        auto & mon = m_monomials.back();
+        mon.var = c.add_var(sz);
+        mon.def = c.value(rational(1), sz);
+        for (auto p : m_tmp) {
+            mon.arg_vals.push_back(rational(0));
+            mon.def *= p;
+        }
+        m_table.insert(index);
+        c.trail().push(del_monomial(*this));
+        return offset * mon.var;
+    }
+
+    pdd monomials::subst(pdd const& p) {
+        pdd r = p;
+        for (unsigned i = m_monomials.size(); i-- > 0;) 
+            r = r.subst_pdd(m_monomials[i].var.var(), m_monomials[i].def);        
+        return r;
+    }
+        
+    lbool monomials::refine() {
+        init_to_refine();
+        if (m_to_refine.empty())
+            return l_true;
+        shuffle(m_to_refine.size(), m_to_refine.data(), c.rand());
+        if (any_of(m_to_refine, [&](auto i) { return mul0(m_monomials[i]); }))
+            return l_false;
+        if (any_of(m_to_refine, [&](auto i) { return non_overflow_unit(m_monomials[i]); }))
+            return l_false;
+        if (any_of(m_to_refine, [&](auto i) { return parity0(m_monomials[i]); }))
+            return l_false;
+        if (any_of(m_to_refine, [&](auto i) { return prefix_overflow(m_monomials[i]); }))
+            return l_false;
+        if (any_of(m_to_refine, [&](auto i) { return non_overflow_monotone(m_monomials[i]); }))
+            return l_false;
+        if (any_of(m_to_refine, [&](auto i) { return mulp2(m_monomials[i]); }))
+            return l_false;
+        if (any_of(m_to_refine, [&](auto i) { return parity(m_monomials[i]); }))
+            return l_false;
+        return l_undef;
+    }
+
+    // bit blast a monomial definition
+    lbool monomials::bit_blast() {
+        init_to_refine();
+        if (m_to_refine.empty())
+            return l_true;
+        shuffle(m_to_refine.size(), m_to_refine.data(), c.rand());
+        if (any_of(m_to_refine, [&](auto i) { return bit_blast(m_monomials[i]); }))
+            return l_false;
+        return l_undef;
+    }
+
+    void monomials::init_to_refine() {
+        m_to_refine.reset();
+        for (unsigned i = 0; i < m_monomials.size(); ++i) 
+            if (eval_to_false(i))
+                m_to_refine.push_back(i);
+    }
+
+    bool monomials::eval_to_false(unsigned i) {
+        rational rhs, lhs, p_val;
+        auto& mon = m_monomials[i];
+        if (!c.try_eval(mon.var, mon.val))
+            return false;
+        if (!c.try_eval(mon.def, rhs))
+            return false;
+        if (rhs == mon.val)
+            return false;
+        for (unsigned j = mon.size(); j-- > 0; ) {
+            if (!c.try_eval(mon.args[j], p_val))
+                return false;                       
+            mon.arg_vals[j] = p_val;
+        }
+        return true;
+    }
+
+    // check p = 0 => p * q = 0
+    bool monomials::mul0(monomial const& mon) {
+        for (unsigned j = mon.size(); j-- > 0; ) {
+            if (mon.arg_vals[j] == 0) {
+                auto const& p = mon.args[j];
+                c.add_axiom("p = 0 => p * q = 0", { ~C.eq(p), C.eq(mon.var) }, true);
+                return true;
+            }                
+        }
+        return false;
+    }
+
+    // check p = k => p * q = k * q
+    bool monomials::mulp2(monomial const& mon) {
+        unsigned free_index = UINT_MAX;
+        auto& m = mon.args[0].manager();
+        for (unsigned j = mon.size(); j-- > 0; ) {
+            auto const& arg_val = mon.arg_vals[j];
+            if (arg_val == m.max_value() || arg_val.is_power_of_two())
+                continue;
+            if (free_index != UINT_MAX)
+                return false;
+            free_index = j;
+        }
+        constraint_or_dependency_vector cs;
+        pdd offset = c.value(rational(1), m.power_of_2());
+        for (unsigned j = mon.size(); j-- > 0; ) {
+            if (j != free_index) {
+                cs.push_back(~C.eq(mon.args[j], mon.arg_vals[j]));
+                offset *= mon.arg_vals[j];
+            }
+        }
+        if (free_index == UINT_MAX)
+            cs.push_back(C.eq(mon.var, offset));
+        else
+            cs.push_back(C.eq(mon.var, offset * mon.args[free_index]));
+        
+        c.add_axiom("p = k => p * q = k * q", cs, true);
+        return true;
+    }
+
+    // parity p >= i => parity p * q >= i
+    bool monomials::parity(monomial const& mon) {       
+        unsigned parity_val = get_parity(mon.val, mon.num_bits());
+        for (unsigned j = 0; j < mon.args.size(); ++j) {
+            unsigned k = get_parity(mon.arg_vals[j], mon.num_bits());
+            if (k > parity_val) {
+                auto const& p = mon.args[j];
+                c.add_axiom("parity p >= i => parity p * q >= i", { ~C.parity_at_least(p, k), C.parity_at_least(mon.var, k) }, true);
+                return true;
+            }
+        }
+        return false;
+    }
+
+    // ~ovfl*(p,q) & q != 0 => p <= p*q    
+    bool monomials::non_overflow_monotone(monomial const& mon) {
+        rational product(1);
+        unsigned big_index = UINT_MAX;
+        for (unsigned i = 0; i < mon.args.size(); ++i) {
+            auto const& val = mon.arg_vals[i];
+            product *= val;
+            if (val > mon.val)
+                big_index = i;
+        }
+        if (big_index == UINT_MAX)
+            return false;
+        if (product > mon.var.manager().max_value())
+            return false;
+        pdd p = mon.args[0];
+        constraint_or_dependency_vector clause;
+        for (unsigned i = 1; i < mon.args.size(); ++i) {
+            clause.push_back(~C.umul_ovfl(p, mon.args[i]));
+            p *= mon.args[i];
+        }
+        for (unsigned i = 0; i < mon.args.size(); ++i)
+            if (i != big_index)
+                clause.push_back(~C.eq(mon.args[i]));
+        clause.push_back(C.ule(mon.args[big_index], mon.var));
+        return false;
+    }
+
+    // ~ovfl*(p,q) & p*q = 1 => p = 1, q = 1
+    bool monomials::non_overflow_unit(monomial const& mon) {
+        if (mon.val != 1)
+            return false;
+        rational product(1);
+        for (auto const& val : mon.arg_vals)
+            product *= val;
+        if (product > mon.var.manager().max_value())
+            return false;
+        constraint_or_dependency_vector clause;
+        pdd p = mon.args[0];
+        clause.push_back(~C.eq(mon.var, 1));
+        for (unsigned i = 1; i < mon.args.size(); ++i) {
+            clause.push_back(~C.umul_ovfl(p, mon.args[i]));
+            p *= mon.args[i];
+        }        
+        for (auto const& q : mon.args) {
+            clause.push_back(C.eq(q, 1));
+            c.add_axiom("~ovfl*(p,q) & p*q = 1 => p = 1", clause, true);
+            clause.pop_back();
+        }
+        return true;
+    }
+
+    // p * q = 0 => p = 0 or even(q)
+    // p * q = 0 => p = 0 or q = 0 or even(q)
+    bool monomials::parity0(monomial const& mon) {
+        if (mon.val != 0)
+            return false;
+        constraint_or_dependency_vector clause;
+        clause.push_back(~C.eq(mon.var, 0));
+        for (auto const& val : mon.arg_vals)
+            if (!val.is_odd() || val == 0)
+                return false;
+        for (auto const& p : mon.args) 
+            clause.push_back(C.eq(p));        
+        for (auto const& p : mon.args) {
+            clause.push_back(C.parity_at_least(p, 1));
+            c.add_axiom("p * q = 0 => p = 0 or q = 0 or even(q)", clause, true);
+            clause.pop_back();
+        }
+        return false;
+    }
+
+    // 0p * 0q >= 2^k => ovfl(p,q), where |p| = |q| = k
+    bool monomials::prefix_overflow(monomial const& mon) {
+        if (mon.size() != 2)
+            return false;
+        if (!mon.args[0].is_var())
+            return false;
+        if (!mon.args[1].is_var())
+            return false;
+        if (mon.val <= mon.arg_vals[0])
+            return false;
+        if (mon.val <= mon.arg_vals[1])
+            return false;
+        auto x = mon.args[0].var();
+        auto y = mon.args[1].var();
+        offset_slices x_suffixes, y_suffixes;
+        c.get_bitvector_suffixes(x, x_suffixes);
+        c.get_bitvector_suffixes(y, y_suffixes);
+        rational x_val, y_val;
+        for (auto const& xslice : x_suffixes) {
+            if (c.size(xslice.v) == mon.num_bits())
+                continue;
+            auto const& xmax_value = c.var(xslice.v).manager().max_value();
+            if (mon.val <= xmax_value)
+                continue;
+            if (!c.try_eval(c.var(xslice.v), x_val) || x_val != mon.arg_vals[0])
+                continue;
+            for (auto const& yslice : y_suffixes) {
+                if (c.size(yslice.v) != c.size(xslice.v))
+                    continue;
+                if (!c.try_eval(c.var(yslice.v), y_val) || y_val != mon.arg_vals[1])
+                    continue;
+                c.add_axiom("0p * 0q >= 2^k => ovfl(p,q), where |p| = |q| = k",
+                    { dependency({x, xslice}), dependency({y, yslice}),
+                     ~C.ule(mon.args[0], xmax_value),
+                     ~C.ule(mon.args[1], xmax_value),
+                     ~C.ugt(mon.var, xmax_value), 
+                      C.umul_ovfl(c.var(xslice.v), c.var(yslice.v)) }, 
+                    true);
+                return true;
+            }
+        }       
+        return false;
+    }
+
+    bool monomials::bit_blast(monomial const& mon) {
+        return false;
+    }
+
+    std::ostream& monomials::display(std::ostream& out) const {
+        for (auto const& mon : m_monomials) {
+            out << mon.var << " := ";
+            char const* sep = "";
+            for (auto p : mon.args) 
+                if (p.is_var())
+                    out << sep << p, sep = " * ";
+                else 
+                    out << sep << "(" << p << ")", sep = " * ";
+            out << "\n";
+        }
+        return out;
+    }
+}
diff --git a/src/sat/smt/polysat/monomials.h b/src/sat/smt/polysat/monomials.h
new file mode 100644
index 000000000..65d16ee96
--- /dev/null
+++ b/src/sat/smt/polysat/monomials.h
@@ -0,0 +1,102 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    Polysat monomials
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2021-03-19
+    Jakob Rath 2021-04-6
+
+--*/
+#pragma once
+
+#include "math/dd/dd_pdd.h"
+#include "sat/smt/polysat/types.h"
+
+namespace polysat {
+
+    class core;
+    class constraints;
+
+    using pdd_vector = vector<pdd>;
+    using rational_vector = vector<rational>;
+    
+    class monomials {
+        core& c;
+        constraints& C;
+
+        struct monomial {
+            pdd_vector args;
+            pdd        var;
+            pdd        def;
+            rational_vector arg_vals;
+            rational   val;
+            unsigned size() const { return args.size(); }
+            unsigned num_bits() const { return args[0].manager().power_of_2(); }
+        };
+        vector<monomial>   m_monomials;
+        pdd_vector         m_tmp;
+        
+        struct mon_eq {
+            monomials& m;
+            mon_eq(monomials& m):m(m) {}
+            bool operator()(unsigned i, unsigned j) const {
+                auto const& a = m.m_monomials[i].args;
+                auto const& b = m.m_monomials[j].args;
+                return a == b;
+            };
+        };
+
+        struct pdd_hash {
+            typedef pdd data_t;
+            unsigned operator()(pdd const& p) const { return p.hash(); }
+        };
+        struct mon_hash {
+            monomials& m;
+            mon_hash(monomials& m):m(m) {}
+            unsigned operator()(unsigned i) const {
+                auto const& a = m.m_monomials[i].args;
+                return vector_hash<pdd_hash>()(a);
+            }
+        };
+        mon_hash m_hash;
+        mon_eq   m_eq;
+        hashtable<unsigned, mon_hash, mon_eq> m_table;
+
+        unsigned_vector m_to_refine;
+        void init_to_refine();
+        bool eval_to_false(unsigned i);
+
+        bool mul0(monomial const& mon);
+        bool mulp2(monomial const& mon);
+        bool parity(monomial const& mon);
+        bool non_overflow_monotone(monomial const& mon);
+        bool non_overflow_unit(monomial const& mon);
+        bool parity0(monomial const& mon);
+        bool prefix_overflow(monomial const& mon);
+
+        bool bit_blast(monomial const& mon);
+        
+    public:
+
+        monomials(core& c);
+
+        pdd mk(unsigned n, pdd const* args);
+
+        pdd subst(pdd const& p);
+        
+        lbool refine();
+
+        lbool bit_blast();
+
+        std::ostream& display(std::ostream& out) const;
+               
+    };
+
+    inline std::ostream& operator<<(std::ostream& out, monomials const& m) {
+        return m.display(out);
+    }
+}
diff --git a/src/sat/smt/polysat/saturation.cpp b/src/sat/smt/polysat/saturation.cpp
index db0bc8410..015d34c16 100644
--- a/src/sat/smt/polysat/saturation.cpp
+++ b/src/sat/smt/polysat/saturation.cpp
@@ -35,7 +35,7 @@ namespace polysat {
     saturation::saturation(core& c) : c(c), C(c.cs()) {}
 
     bool saturation::resolve(pvar v, constraint_id id) {
-        if (c.eval(id) == l_true)
+        if (c.eval_unfold(id) == l_true)
             return false;
         auto sc = c.get_constraint(id);
         if (!sc.vars().contains(v))
diff --git a/src/sat/smt/polysat/types.h b/src/sat/smt/polysat/types.h
index f9910deab..8fa891b66 100644
--- a/src/sat/smt/polysat/types.h
+++ b/src/sat/smt/polysat/types.h
@@ -124,6 +124,7 @@ namespace polysat {
     using constraint_id_or_dependency = std::variant<constraint_id, dependency>;
 
     using constraint_or_dependency_list = std::initializer_list<constraint_or_dependency>;
+    using constraint_or_dependency_vector = vector<constraint_or_dependency>;
     using constraint_id_vector = svector<constraint_id>;
     using constraint_id_list = std::initializer_list<constraint_id>;
     using offset_slices = vector<offset_slice>;
diff --git a/src/sat/smt/polysat/ule_constraint.cpp b/src/sat/smt/polysat/ule_constraint.cpp
index 60b7dcd61..802cea772 100644
--- a/src/sat/smt/polysat/ule_constraint.cpp
+++ b/src/sat/smt/polysat/ule_constraint.cpp
@@ -238,8 +238,8 @@ namespace polysat {
 
 namespace polysat {
 
-    ule_constraint::ule_constraint(pdd const& l, pdd const& r) :
-        m_lhs(l), m_rhs(r) {
+    ule_constraint::ule_constraint(pdd const& l, pdd const& r, pdd const& ul, pdd const& ur) :
+        m_lhs(l), m_rhs(r), m_unfold_lhs(ul), m_unfold_rhs(ur) {
 
         SASSERT_EQ(m_lhs.power_of_2(), m_rhs.power_of_2());
 
@@ -247,6 +247,11 @@ namespace polysat {
         for (auto v : m_rhs.free_vars())
             if (!vars().contains(v))
                 vars().push_back(v);
+
+        m_unfold_vars.append(m_unfold_lhs.free_vars());
+        for (auto v : m_unfold_rhs.free_vars())
+            if (!m_unfold_vars.contains(v))
+                m_unfold_vars.push_back(v);
     }
 
     void ule_constraint::simplify(bool& is_positive, pdd& lhs, pdd& rhs) {
@@ -312,7 +317,10 @@ namespace polysat {
     }
 
     std::ostream& ule_constraint::display(std::ostream& out) const {
-        return display(out, l_true, m_lhs, m_rhs);
+        display(out, l_true, m_lhs, m_rhs);
+        if (m_lhs != m_unfold_lhs || m_rhs != m_unfold_rhs)
+            display(out << " alias ( ", l_true, m_unfold_lhs, m_unfold_rhs) << ")";
+        return out;
     }
 
     // Evaluate lhs <= rhs
@@ -342,6 +350,10 @@ namespace polysat {
         return eval(a.apply_to(lhs()), a.apply_to(rhs()));
     }
 
+    lbool ule_constraint::eval_unfold(assignment const& a) const {
+        return eval(a.apply_to(unfold_lhs()), a.apply_to(unfold_rhs()));
+    }
+
     void ule_constraint::activate(core& c, bool sign, dependency const& d) {
         auto p = c.subst(lhs());
         auto q = c.subst(rhs());
diff --git a/src/sat/smt/polysat/ule_constraint.h b/src/sat/smt/polysat/ule_constraint.h
index 01ab9570f..8b838bcbf 100644
--- a/src/sat/smt/polysat/ule_constraint.h
+++ b/src/sat/smt/polysat/ule_constraint.h
@@ -19,22 +19,27 @@ Author:
 namespace polysat {
 
     class ule_constraint final : public constraint {
-        pdd m_lhs;
-        pdd m_rhs;        
+        pdd m_lhs, m_rhs;
+        unsigned_vector m_unfold_vars;
+        pdd m_unfold_lhs, m_unfold_rhs;        
         static bool is_always_true(bool is_positive, pdd const& lhs, pdd const& rhs) { return eval(lhs, rhs) == to_lbool(is_positive); }
         static bool is_always_false(bool is_positive, pdd const& lhs, pdd const& rhs) { return is_always_true(!is_positive, lhs, rhs); }
         static lbool eval(pdd const& lhs, pdd const& rhs);
 
     public:
-        ule_constraint(pdd const& l, pdd const& r);
+        ule_constraint(pdd const& l, pdd const& r, pdd const& ul, pdd const& ur);
         ~ule_constraint() override {}
         pdd const& lhs() const { return m_lhs; }
         pdd const& rhs() const { return m_rhs; }
+        pdd const& unfold_lhs() const { return m_unfold_lhs; }
+        pdd const& unfold_rhs() const { return m_unfold_rhs; }
         static std::ostream& display(std::ostream& out, lbool status, pdd const& lhs, pdd const& rhs);
+        unsigned_vector const& unfold_vars() const override { return m_unfold_vars; }
         std::ostream& display(std::ostream& out, lbool status) const override;
         std::ostream& display(std::ostream& out) const override;
         lbool eval() const override;
         lbool eval(assignment const& a) const override;
+        lbool eval_unfold(assignment const& a) const override;
         bool is_linear() const override { return lhs().is_linear() && rhs().is_linear(); }
         void activate(core& c, bool sign, dependency const& dep);
         void propagate(core& c, lbool value, dependency const& dep) {}