diff --git a/src/ast/bv_decl_plugin.cpp b/src/ast/bv_decl_plugin.cpp
index 8f0ebe3d1..243e4273a 100644
--- a/src/ast/bv_decl_plugin.cpp
+++ b/src/ast/bv_decl_plugin.cpp
@@ -841,28 +841,6 @@ bool bv_recognizers::is_bit2bool(expr* e, expr*& bv, unsigned& idx) const {
     return true;
 }
 
-bool bv_recognizers::mult_inverse(rational const & n, unsigned bv_size, rational & result) {
-    if (n.is_one()) {
-        result = n;
-        return true;
-    }
-
-    if (!mod(n, rational(2)).is_one()) {
-        return false;
-    }
-
-    rational g;
-    rational x;
-    rational y;
-    g = gcd(n, rational::power_of_two(bv_size), x, y);
-    if (x.is_neg()) {
-        x = mod(x, rational::power_of_two(bv_size));
-    }
-    SASSERT(x.is_pos());
-    SASSERT(mod(x * n, rational::power_of_two(bv_size)).is_one());
-    result = x;
-    return true;
-}
 
 bv_util::bv_util(ast_manager & m):
     bv_recognizers(m.mk_family_id(symbol("bv"))),
diff --git a/src/ast/bv_decl_plugin.h b/src/ast/bv_decl_plugin.h
index cf1ef1cdd..be8c3f55c 100644
--- a/src/ast/bv_decl_plugin.h
+++ b/src/ast/bv_decl_plugin.h
@@ -379,7 +379,6 @@ public:
     rational norm(rational const & val, unsigned bv_size, bool is_signed) const ;
     rational norm(rational const & val, unsigned bv_size) const { return norm(val, bv_size, false); }
     bool has_sign_bit(rational const & n, unsigned bv_size) const;
-    bool mult_inverse(rational const & n, unsigned bv_size, rational & result);
 };
 
 class bv_util : public bv_recognizers {
diff --git a/src/ast/rewriter/bv_rewriter.cpp b/src/ast/rewriter/bv_rewriter.cpp
index 1f9395c63..1543ab1a2 100644
--- a/src/ast/rewriter/bv_rewriter.cpp
+++ b/src/ast/rewriter/bv_rewriter.cpp
@@ -2520,7 +2520,7 @@ br_status bv_rewriter::mk_mul_eq(expr * lhs, expr * rhs, expr_ref & result) {
     unsigned sz;
     if (m_util.is_bv_mul(lhs, c, x) &&
         m_util.is_numeral(c, c_val, sz) &&
-        m_util.mult_inverse(c_val, sz, c_inv_val)) {
+        c_val.mult_inverse(sz, c_inv_val)) {
 
         SASSERT(m_util.norm(c_val * c_inv_val, sz).is_one());
 
@@ -2562,7 +2562,7 @@ br_status bv_rewriter::mk_mul_eq(expr * lhs, expr * rhs, expr_ref & result) {
         for (; !found && i < num_args; ++i) {
             expr* arg = to_app(rhs)->get_arg(i);
             if (m_util.is_bv_mul(arg, c2, x2) && m_util.is_numeral(c2, c2_val, sz) &&
-                m_util.mult_inverse(c2_val, sz, c2_inv_val)) {
+                c2_val.mult_inverse(sz, c2_inv_val)) {
                 found = true;
             }
         }
diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp
index 224de5ee2..7fe80b214 100644
--- a/src/math/dd/dd_pdd.cpp
+++ b/src/math/dd/dd_pdd.cpp
@@ -23,7 +23,7 @@ Revision History:
 
 namespace dd {
 
-    pdd_manager::pdd_manager(unsigned num_vars, semantics s) {
+    pdd_manager::pdd_manager(unsigned num_vars, semantics s, unsigned power_of_2) {
         m_spare_entry = nullptr;
         m_max_num_nodes = 1 << 24; // up to 16M nodes
         m_mark_level = 0;
@@ -31,6 +31,8 @@ namespace dd {
         m_disable_gc = false;
         m_is_new_node = false;
         m_semantics = s;
+        m_mod2N = rational::power_of_two(power_of_2);
+        m_power_of_2 = power_of_2;
         unsigned_vector l2v;
         for (unsigned i = 0; i < num_vars; ++i) l2v.push_back(i);
         init_nodes(l2v);
@@ -109,6 +111,56 @@ namespace dd {
     pdd pdd_manager::mk_xor(pdd const& p, unsigned x) { pdd q(mk_val(x)); if (m_semantics == mod2_e) return p + q; return (p*q*2) - p - q; }
     pdd pdd_manager::mk_not(pdd const& p) { return 1 - p; }
 
+    pdd pdd_manager::subst_val(pdd const& p, unsigned v, rational const& val) {
+        pdd r = mk_var(v) + val;
+        return pdd(apply(p.root, r.root, pdd_subst_val_op), this);
+    }
+
+
+    /**
+     * A polynomial is non-zero if the constant coefficient
+     * is a power of two such that none of the coefficients to
+     * non-constant monomials divide it.
+     * Example: 2^4*x + 2 is non-zero for every x.
+     */
+
+    bool pdd_manager::is_non_zero(PDD p) {
+        if (is_val(p))
+            return !is_zero(p);
+        if (m_semantics != mod2N_e)
+            return false;
+        PDD q = p;
+        while (!is_val(q))
+            q = lo(q);
+        auto const& v = val(q);
+        if (v.is_zero())
+            return false;
+        unsigned p2 = v.trailing_zeros();
+        init_mark();
+        SASSERT(m_todo.empty());
+        m_todo.push_back(hi(p));
+        while (!is_val(lo(p))) {
+            p = lo(p);
+            m_todo.push_back(hi(p));
+        }
+        while (!m_todo.empty()) {
+            PDD r = m_todo.back();
+            m_todo.pop_back();
+            if (is_marked(r)) 
+                continue;
+            set_mark(r);
+            if (!is_val(r)) {
+                m_todo.push_back(lo(r));
+                m_todo.push_back(hi(r));
+            }
+            else if (val(r).trailing_zeros() <= p2) {
+                m_todo.reset();
+                return false;
+            }
+        }
+        return true;
+    }
+
     pdd pdd_manager::subst_val(pdd const& p, vector<std::pair<unsigned, rational>> const& _s) {
         typedef std::pair<unsigned, rational> pr;
         vector<pr> s(_s);        
@@ -116,9 +168,8 @@ namespace dd {
             [&](pr const& a, pr const& b) { return m_var2level[a.first] < m_var2level[b.first]; };
         std::sort(s.begin(), s.end(), compare_level);
         pdd r(one());
-        for (auto const& q : s) {
+        for (auto const& q : s) 
             r = (r*mk_var(q.first)) + q.second;            
-        }
         return pdd(apply(p.root, r.root, pdd_subst_val_op), this);
     }
 
@@ -256,7 +307,7 @@ namespace dd {
                 r = make_node(level_p, read(2), read(1));
             }
             else if (level_p == level_q) {
-                if (m_semantics != free_e) {
+                if (m_semantics != free_e && m_semantics != mod2N_e) {
                     //
                     // (xa+b)*(xc+d)  == x(ac+bc+ad) + bd
                     //                == x((a+b)(c+d)-bd) + bd
@@ -334,9 +385,9 @@ namespace dd {
         case pdd_subst_val_op:
             SASSERT(!is_val(p));
             SASSERT(!is_val(q));
-            SASSERT(level_p = level_q);
-            push(apply_rec(lo(p), hi(q), pdd_subst_val_op));   // lo := subst(lo(p), s)
-            push(apply_rec(hi(p), hi(q), pdd_subst_val_op));   // hi := subst(hi(p), s)
+            SASSERT(level_p == level_q);
+            push(apply_rec(lo(p), q, pdd_subst_val_op));   // lo := subst(lo(p), s)
+            push(apply_rec(hi(p), q, pdd_subst_val_op));   // hi := subst(hi(p), s)
             push(apply_rec(lo(q), read(1), pdd_mul_op));       // hi := hi*s[var(p)]
             r = apply_rec(read(1), read(3), pdd_add_op);       // r := hi + lo := subst(lo(p),s) + s[var(p)]*subst(hi(p),s)
             npop = 3;
@@ -721,9 +772,14 @@ namespace dd {
     }
 
     pdd_manager::PDD pdd_manager::imk_val(rational const& r) {
-        if (r.is_zero()) return zero_pdd;
-        if (r.is_one()) return one_pdd;
-        if (m_semantics == mod2_e) return imk_val(mod(r, rational(2)));
+        if (r.is_zero()) 
+            return zero_pdd;
+        if (r.is_one()) 
+            return one_pdd;
+        if (m_semantics == mod2_e) 
+            return imk_val(mod(r, rational(2)));
+        if (m_semantics == mod2N_e && (r < 0 || r >= m_mod2N)) 
+            return imk_val(mod(r, m_mod2N));
         const_info info;
         if (!m_mpq_table.find(r, info)) {
             init_value(info, r);
@@ -1073,13 +1129,12 @@ namespace dd {
         auto mons = to_monomials(b);
         bool first = true;
         for (auto& m : mons) {
-            if (!first) {
-                if (m.first.is_neg()) out << " - ";
-                else out << " + ";
-            }
-            else {
-                if (m.first.is_neg()) out << "- ";
-            }
+            if (!first)
+                out << " ";
+            if (m.first.is_neg()) 
+                out << "- ";
+            else if (!first)
+                out << "+ ";
             first = false;
             rational c = abs(m.first);
             m.second.reverse();
diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h
index 20b04fa39..a8e65c0bf 100644
--- a/src/math/dd/dd_pdd.h
+++ b/src/math/dd/dd_pdd.h
@@ -43,7 +43,7 @@ namespace dd {
 
     class pdd_manager {
     public:
-        enum semantics { free_e, mod2_e, zero_one_vars_e };
+        enum semantics { free_e, mod2_e, zero_one_vars_e, mod2N_e };
     private:
         friend test;
         friend pdd;
@@ -140,7 +140,7 @@ namespace dd {
 
         typedef ptr_hashtable<op_entry, hash_entry, eq_entry> op_table;
 
-        svector<node>          m_nodes;
+        svector<node>              m_nodes;
         vector<rational>           m_values;
         op_table                   m_op_cache;
         node_table                 m_node_table;
@@ -163,6 +163,8 @@ namespace dd {
         unsigned_vector            m_free_vars;
         unsigned_vector            m_free_values;
         rational                   m_freeze_value;
+        rational                   m_mod2N;
+        unsigned                   m_power_of_2 { 0 };
 
         void reset_op_cache();
         void init_nodes(unsigned_vector const& l2v);
@@ -204,6 +206,7 @@ namespace dd {
         inline bool is_one(PDD p) const { return p == one_pdd; } 
         inline bool is_val(PDD p) const { return m_nodes[p].is_val(); }
         inline bool is_internal(PDD p) const { return m_nodes[p].is_internal(); }
+        bool is_non_zero(PDD p);
         inline unsigned level(PDD p) const { return m_nodes[p].m_level; }
         inline unsigned var(PDD p) const { return m_level2var[level(p)]; }
         inline PDD lo(PDD p) const { return m_nodes[p].m_lo; }
@@ -252,7 +255,7 @@ namespace dd {
 
         struct mem_out {};
 
-        pdd_manager(unsigned nodes, semantics s = free_e);
+        pdd_manager(unsigned nodes, semantics s = free_e, unsigned power_of_2 = 0);
         ~pdd_manager();
 
         semantics get_semantics() const { return m_semantics; }
@@ -278,6 +281,7 @@ namespace dd {
         pdd mk_not(pdd const& p);
         pdd reduce(pdd const& a, pdd const& b);
         pdd subst_val(pdd const& a, vector<std::pair<unsigned, rational>> const& s);
+        pdd subst_val(pdd const& a, unsigned v, rational const& val);
 
         bool is_linear(PDD p) { return degree(p) == 1; }
         bool is_linear(pdd const& p);
@@ -300,6 +304,8 @@ namespace dd {
         double tree_size(pdd const& p);
         unsigned num_vars() const { return m_var2pdd.size(); }
 
+        unsigned power_of_2() const { return m_power_of_2; }
+
         unsigned_vector const& free_vars(pdd const& p);
 
         std::ostream& display(std::ostream& out);
@@ -334,6 +340,7 @@ namespace dd {
         bool is_unary() const { return !is_val() && lo().is_zero() && hi().is_val(); } 
         bool is_binary() const { return m.is_binary(root); }
         bool is_monomial() const { return m.is_monomial(root); }
+        bool is_non_zero() const { return m.is_non_zero(root); }
         bool var_is_leaf(unsigned v) const { return m.var_is_leaf(root, v); }
 
         pdd operator-() const { return m.minus(*this); }
@@ -353,11 +360,14 @@ namespace dd {
         bool different_leading_term(pdd const& other) const { return m.different_leading_term(*this, other); }
 
         pdd subst_val(vector<std::pair<unsigned, rational>> const& s) const { return m.subst_val(*this, s); }
+        pdd subst_val(unsigned v, rational const& val) const { return m.subst_val(*this, v, val); }
 
         std::ostream& display(std::ostream& out) const { return m.display(out, *this); }
         bool operator==(pdd const& other) const { return root == other.root; }
         bool operator!=(pdd const& other) const { return root != other.root; }
 
+        unsigned power_of_2() const { return m.power_of_2(); }
+
         unsigned dag_size() const { return m.dag_size(*this); }
         double tree_size() const { return m.tree_size(*this); }
         unsigned degree() const { return m.degree(*this); }
diff --git a/src/math/polysat/polysat.cpp b/src/math/polysat/polysat.cpp
index 884888964..73eddc545 100644
--- a/src/math/polysat/polysat.cpp
+++ b/src/math/polysat/polysat.cpp
@@ -19,28 +19,30 @@ Author:
 
 namespace polysat {
 
-    std::ostream& poly::display(std::ostream& out) const {
-        return out;
-    }
-
     std::ostream& constraint::display(std::ostream& out) const {
-        return out;
+        return out << "constraint";
     }
     
     std::ostream& linear::display(std::ostream& out) const {
-        return out;
+        return out << "linear";
     }
 
     std::ostream& mono::display(std::ostream& out) const {
-        return out;
+        return out << "mono";
     }
 
-    unsigned solver::poly2size(poly const& p) const {
-        return 0;
+    dd::pdd_manager& solver::sz2pdd(unsigned sz) {
+        m_pdd.reserve(sz + 1);
+        if (!m_pdd[sz])
+            m_pdd.set(sz, alloc(dd::pdd_manager, 1000));
+        return *m_pdd[sz];
     }
 
-    bool solver::is_viable(unsigned var, rational const& val) const {
-        return false;
+    bool solver::is_viable(unsigned var, rational const& val) {
+        bdd b = m_viable[var];
+        for (unsigned k = size(var); k-- > 0 && !b.is_false(); ) 
+            b &= val.get_bit(k) ? m_bdd.mk_var(k) : m_bdd.mk_nvar(k);
+        return !b.is_false();
     }
 
     struct solver::del_var : public trail {
@@ -60,11 +62,10 @@ namespace polysat {
         var_unassign(solver& s): s(s) {}
         void undo() override { s.do_var_unassign(); }
     };
+
     
     solver::solver(trail_stack& s): 
         m_trail(s),
-        m_region(s.get_region()),
-        m_pdd(1000),
         m_bdd(1000),
         m_free_vars(m_activity) {
     }
@@ -77,20 +78,21 @@ namespace polysat {
         
     unsigned solver::add_var(unsigned sz) {
         unsigned v = m_viable.size();
-        // TODO: set var size sz into m_pdd.
         m_value.push_back(rational::zero());
         m_justification.push_back(justification::unassigned());
         m_viable.push_back(m_bdd.mk_true());
         m_vdeps.push_back(m_dep_manager.mk_empty());
-        m_pdeps.push_back(vector<poly>());
-        m_watch.push_back(unsigned_vector());
+        m_pdeps.push_back(vector<pdd>());
+        m_watch.push_back(ptr_vector<constraint>());
         m_activity.push_back(0);
-        m_vars.push_back(m_pdd.mk_var(v));
+        m_vars.push_back(sz2pdd(sz).mk_var(v));
+        m_size.push_back(sz);
         m_trail.push(del_var(*this));
         return v;
     }
 
     void solver::do_del_var() {
+        // TODO also remove v from all learned constraints.
         unsigned v = m_viable.size() - 1;
         m_free_vars.del_var_eh(v);
         m_viable.pop_back();
@@ -101,34 +103,42 @@ namespace polysat {
         m_watch.pop_back();
         m_activity.pop_back();
         m_vars.pop_back();
+        m_size.pop_back();
     }
 
     void solver::do_del_constraint() {
-        // TODO remove variables from watch lists
+        // TODO rewrite to allow for learned constraints
+        // so have to gc these.
+        constraint& c = *m_constraints.back();
+        if (c.vars().size() > 0)
+            erase_watch(c.vars()[0], c);
+        if (c.vars().size() > 1)
+            erase_watch(c.vars()[1], c);
         m_constraints.pop_back();
-        m_cdeps.pop_back();
     }
 
     void solver::do_var_unassign() {
-
+        unsigned v = m_search.back();
+        m_justification[v] = justification::unassigned();
+        m_free_vars.unassign_var_eh(v);
     }
 
-    poly solver::var(unsigned v) {
-        return poly(*this, m_vars[v]);
-    }
-
-    vector<mono> solver::poly2monos(poly const& p) const {
-        return vector<mono>();
-    }
-
-    void solver::add_eq(poly const& p, unsigned dep) {
-        m_constraints.push_back(constraint::eq(p));
-        m_cdeps.push_back(m_dep_manager.mk_leaf(dep));
-        // TODO init watch list
+    void solver::add_eq(pdd const& p, unsigned dep) {
+        //
+        // TODO reduce p using assignment (at current level, 
+        // assuming constraint is removed also at current level).
+        // 
+        constraint* c = constraint::eq(p, m_dep_manager.mk_leaf(dep));
+        m_constraints.push_back(c);
+        auto const& vars = c->vars();
+        if (vars.size() > 0) 
+            m_watch[vars[0]].push_back(c);
+        if (vars.size() > 1) 
+            m_watch[vars[1]].push_back(c);
         m_trail.push(del_constraint(*this));
     }
 
-    void solver::add_diseq(poly const& p, unsigned dep) {
+    void solver::add_diseq(pdd const& p, unsigned dep) {
 #if 0
         // Basically: 
         auto slack = add_var(p.size());
@@ -138,40 +148,175 @@ namespace polysat {
 #endif
     }
 
-    void solver::add_ule(poly const& p, poly const& q, unsigned dep) {
+    void solver::add_ule(pdd const& p, pdd const& q, unsigned dep) {
         // save for later
     }
 
-    void solver::add_sle(poly const& p, poly const& q, unsigned dep) {
+    void solver::add_sle(pdd const& p, pdd const& q, unsigned dep) {
         // save for later
     }
 
     void solver::assign(unsigned var, unsigned index, bool value, unsigned dep) {
-
+        m_viable[var] &= value ? m_bdd.mk_var(index) : m_bdd.mk_nvar(index);
+        m_trail.push(vector_value_trail<u_dependency*, false>(m_vdeps, var));
+        m_vdeps[var] = m_dep_manager.mk_join(m_vdeps[var], m_dep_manager.mk_leaf(dep));
+        if (m_viable[var].is_false()) {
+            // TBD: set conflict
+        }
     }
 
-    bool  solver::can_propagate() {
+    bool solver::can_propagate() {
+        return m_qhead < m_search.size() && !is_conflict();
+    }
+
+    void solver::propagate() {
+        m_trail.push(value_trail(m_qhead));
+        while (can_propagate()) 
+            propagate(m_search[m_qhead++]);
+    }
+
+    void solver::propagate(unsigned v) {
+        auto& wlist = m_watch[v];
+        unsigned i = 0, j = 0, sz = wlist.size();
+        for (; i < sz && !is_conflict(); ++i) 
+            if (!propagate(v, *wlist[i]))
+                wlist[j++] = wlist[i];
+        for (; i < sz; ++i)
+            wlist[j++] = wlist[i];
+        wlist.shrink(j);
+    }
+
+    bool solver::propagate(unsigned v, constraint& c) {
+        switch (c.kind()) {
+        case ckind_t::eq_t:
+            return propagate_eq(v, c);
+        case ckind_t::ule_t:
+        case ckind_t::sle_t:
+            NOT_IMPLEMENTED_YET();
+            return false;
+        }
+        UNREACHABLE();
         return false;
     }
 
-    lbool solver::propagate() {
-        return l_false;
+    bool solver::propagate_eq(unsigned v, constraint& c) {
+        SASSERT(c.kind() == ckind_t::eq_t);
+        SSSERT(!c.vars().empty());
+        auto var = m_vars[v].var();
+        auto& vars = c.vars();
+        unsigned idx = 0;
+        if (vars[idx] != v)
+            idx = 1;
+        SASSERT(v == vars[idx]);
+        // find other watch variable.
+        for (unsigned i = var.size(); i-- > 2; ) {
+            if (!is_assigned(vars[i])) {
+                std::swap(vars[idx], vars[i]);
+                return true;
+            }
+        }        
+
+        vector<std::pair<unsigned, rational>> sub;
+        for (auto w : vars) 
+            if (is_assigned(w))
+                sub.push_back(std::make_pair(w, m_value[w]));
+
+        auto p = c.p().subst_val(sub);
+        if (p.is_zero()) 
+            return false;
+        if (p.is_non_zero()) {
+            // we could tag constraint to allow early substitution before 
+            // swapping watch variable in case we can detect conflict earlier.
+            set_conflict(v, c);
+            return false;
+        }
+
+        // one variable remains unassigned.
+        auto other_var = vars[1 - idx];
+        SASSERT(!is_assigned(other_var));
+
+        // Detect and apply unit propagation.
+            
+        if (!p.is_linear())
+            return false;
+
+        // a*x + b == 0
+        rational a = p.hi().val();
+        rational b = p.lo().val();
+        rational inv_a;
+        if (p.lo().val().is_odd()) {
+            // v1 = -b * inverse(a)
+            unsigned sz = p.power_of_2();
+            VERIFY(a.mult_inverse(sz, inv_a)); 
+            rational val = mod(inv_a * -b, rational::power_of_two(sz));
+            assign(other_var, val, justification::propagation(m_level, &c));
+            return false;
+        }
+
+        // TBD
+        // constrain viable using condition on x
+        // 2*x + 2 == 0 mod 4 => x is odd
+        
+        return false;
+    }
+
+    void solver::inc_level() {
+        m_trail.push(value_trail(m_level));
+        ++m_level;
+    }
+
+    void solver::erase_watch(unsigned v, constraint& c) {
+        if (v == UINT_MAX)
+            return;
+        auto& wlist = m_watch[v];
+        unsigned sz = wlist.size();
+        for (unsigned i = 0; i < sz; ++i) {
+            if (&c == wlist[i]) {
+                wlist[i] = wlist.back();
+                wlist.pop_back();
+                return;
+            }
+        }
     }
 
     bool solver::can_decide() {
-        return false;
+        return !m_free_vars.empty();
     }
 
     void solver::decide(rational & val, unsigned& var) {
-
+        SASSERT(can_decide());
+        inc_level();
+        var = m_free_vars.next_var();
+        auto viable = m_viable[var];
+        SASSERT(!viable.is_false());
+        // TBD, choose some value from viable and construct val.
+        assign_core(var, val, justification::decision(m_level));
     }
 
-    void solver::assign(unsigned var, rational const& val) {
+    void solver::assign(unsigned var, rational const& val, justification const& j) {
+        if (is_viable(var, val)) 
+            assign_core(var, val, j);        
+        else {
+            SASSERT(!j.is_decision());
+            // TBD: set conflict, assumes justification is a propagation.
+        }
+    }
 
+    void solver::assign_core(unsigned var, rational const& val, justification const& j) {
+        SASSERT(is_viable(var, val));
+        m_trail.push(var_unassign(*this));
+        m_search.push_back(var);
+        m_value[var] = val;
+        m_justification[var] = j; 
     }
         
     bool solver::is_conflict() {
-        return false;
+        return m_conflict != nullptr;
+    }
+
+    void solver::set_conflict(unsigned v, constraint& c) {
+        m_conflict = &c;
+        m_conflict_var = v;
     }
 
     unsigned resolve_conflict(unsigned_vector& deps) {
diff --git a/src/math/polysat/polysat.h b/src/math/polysat/polysat.h
index d495077bf..50e8edaea 100644
--- a/src/math/polysat/polysat.h
+++ b/src/math/polysat/polysat.h
@@ -19,6 +19,7 @@ Author:
 #include "util/dependency.h"
 #include "util/trail.h"
 #include "util/lbool.h"
+#include "util/scoped_ptr_vector.h"
 #include "util/var_queue.h"
 #include "math/dd/dd_pdd.h"
 #include "math/dd/dd_bdd.h"
@@ -29,38 +30,33 @@ namespace polysat {
     typedef dd::pdd pdd;
     typedef dd::bdd bdd;
 
-    class poly {
-        friend class solver;
-        solver&  s;
-        pdd      m_pdd;
-    public:
-        poly(solver& s, pdd const& p): s(s), m_pdd(p) {}
-        poly(solver& s, rational const& r, unsigned sz);
-        std::ostream& display(std::ostream& out) const;
-        unsigned size() const { throw default_exception("nyi query pdd for size"); }
-        
-        poly operator*(rational const& r);
-        poly operator+(poly const& other) { return poly(s, m_pdd + other.m_pdd); }
-        poly operator*(poly const& other) { return poly(s, m_pdd * other.m_pdd); }
-    };
-
-    inline std::ostream& operator<<(std::ostream& out, poly const& p) { return p.display(out); }
-
     enum ckind_t { eq_t, ule_t, sle_t };
 
     class constraint {
         ckind_t m_kind;
-        poly    m_poly;
-        poly    m_other;
-        constraint(poly const& p, poly const& q, ckind_t k): m_kind(k), m_poly(p), m_other(q) {}
+        unsigned m_v1, m_v2;
+        pdd    m_poly;
+        pdd    m_other;
+        u_dependency* m_dep;
+        unsigned_vector m_vars;
+        constraint(pdd const& p, pdd const& q, u_dependency* dep, ckind_t k): 
+            m_kind(k), m_v1(UINT_MAX), m_v2(UINT_MAX), m_poly(p), m_other(q), m_dep(dep) {
+            m_vars.append(p.free_vars());
+            if (q != p) 
+                for (auto v : q.free_vars())
+                    m_vars.insert(v);                
+        }
     public:
-        static constraint eq(poly const& p) { return constraint(p, p, ckind_t::eq_t); }
-        static constraint ule(poly const& p, poly const& q) { return constraint(p, q, ckind_t::ule_t); }
+        static constraint* eq(pdd const& p, u_dependency* d) { return alloc(constraint, p, p, d, ckind_t::eq_t); }
+        static constraint* ule(pdd const& p, pdd const& q, u_dependency* d) { return alloc(constraint, p, q, d, ckind_t::ule_t); }
         ckind_t kind() const { return m_kind; }
-        poly const &  p() const { return m_poly; }
-        poly const &  lhs() const { return m_poly; }
-        poly const &  rhs() const { return m_other; }
+        pdd const &  p() const { return m_poly; }
+        pdd const &  lhs() const { return m_poly; }
+        pdd const &  rhs() const { return m_other; }
         std::ostream& display(std::ostream& out) const;
+        void set_poly(pdd const& p) { m_poly = p; }
+        u_dependency* dep() const { return m_dep; }
+        unsigned_vector& vars() { return m_vars; }
     };
 
     inline std::ostream& operator<<(std::ostream& out, constraint const& c) { return c.display(out); }
@@ -98,63 +94,61 @@ namespace polysat {
 
     class justification {
         justification_k m_kind;
-        unsigned        m_idx;
-        justification(justification_k k, unsigned constraint_idx): m_kind(k), m_idx(constraint_idx) {}
+        constraint*     m_c { nullptr };
+        unsigned        m_level;
+        justification(justification_k k, unsigned lvl, constraint* c): m_kind(k), m_c(c), m_level(lvl) {}
     public:
-        justification(): m_kind(justification_k::unassigned), m_idx(0) {}
-        static justification unassigned() { return justification(justification_k::unassigned, 0); }
-        static justification decision() { return justification(justification_k::decision, 0); }
-        static justification propagation(unsigned idx) { return justification(justification_k::propagation, idx); }
+        justification(): m_kind(justification_k::unassigned), m_c(nullptr) {}
+        static justification unassigned() { return justification(justification_k::unassigned, 0, nullptr); }
+        static justification decision(unsigned lvl) { return justification(justification_k::decision, lvl, nullptr); }
+        static justification propagation(unsigned lvl, constraint* c) { return justification(justification_k::propagation, lvl, c); }
+        bool is_decision() const { return m_kind == justification_k::decision; }
+        bool is_unassigned() const { return m_kind == justification_k::unassigned; }
         justification_k kind() const { return m_kind; }
-        unsigned constraint_index() const { return m_idx; }
+        constraint& get_constraint() const { return *m_c; }
+        unsigned level() const { return m_level; }
         std::ostream& display(std::ostream& out) const;
     };
 
     inline std::ostream& operator<<(std::ostream& out, justification const& j) { return j.display(out); }
 
     class solver {
-        friend class poly;
 
         trail_stack&             m_trail;
-        region&                  m_region;
-        dd::pdd_manager          m_pdd;
+        scoped_ptr_vector<dd::pdd_manager> m_pdd;
         dd::bdd_manager          m_bdd;
         u_dependency_manager     m_dep_manager;
         var_queue                m_free_vars;
 
-        /**
-         * store of linear polynomials. The l_idx points to linear monomials.
-         * could also just use pdds.
-         */
-        vector<linear>           m_linear; 
-
         // Per constraint state
-        ptr_vector<u_dependency> m_cdeps;   // each constraint has set of dependencies
-        vector<constraint>       m_constraints;
+        scoped_ptr_vector<constraint>   m_constraints;
+        // TODO: vector<constraint> m_redundant; // learned constraints
 
         // Per variable information
         vector<bdd>              m_viable;   // set of viable values.
         ptr_vector<u_dependency> m_vdeps;    // dependencies for viable values
-        vector<vector<poly>>     m_pdeps;    // dependencies in polynomial form
+        vector<vector<pdd>>      m_pdeps;    // dependencies in polynomial form
         vector<rational>         m_value;    // assigned value
         vector<justification>    m_justification; // justification for variable assignment
-        vector<unsigned_vector>  m_watch;    // watch list datastructure into constraints.
+        vector<ptr_vector<constraint>>  m_watch;    // watch list datastructure into constraints.
         unsigned_vector          m_activity; 
         vector<pdd>              m_vars;
+        unsigned_vector          m_size;     // store size of variables.
 
         // search state that lists assigned variables
         unsigned_vector          m_search;
         unsigned                 m_qhead { 0 };
+        unsigned                 m_level { 0 };
 
-        /**
-         * retrieve bit-size associated with polynomial.
-         */
-        unsigned poly2size(poly const& p) const;
+        // conflict state
+        unsigned    m_conflict_var { UINT_MAX };
+        constraint* m_conflict { nullptr };
 
+        unsigned size(unsigned var) const { return m_size[var]; }
         /**
          * check if value is viable according to m_viable.
          */
-        bool is_viable(unsigned var, rational const& val) const;
+        bool is_viable(unsigned var, rational const& val);
 
         /**
          * undo trail operations for backtracking.
@@ -168,6 +162,23 @@ namespace polysat {
         void do_del_constraint();
         void do_var_unassign();
 
+        dd::pdd_manager& sz2pdd(unsigned sz);
+
+        void inc_level();
+
+        void assign(unsigned var, rational const& val, justification const& j);
+        void assign_core(unsigned var, rational const& val, justification const& j);
+
+        bool is_assigned(unsigned var) const { return !m_justification[var].is_unassigned(); }
+
+        void propagate(unsigned v);
+        bool propagate(unsigned v, constraint& c);
+        bool propagate_eq(unsigned v, constraint& c);
+        void erase_watch(unsigned v, constraint& c);
+
+        void set_conflict(unsigned v, constraint& c);
+        void clear_conflict() { m_conflict = nullptr; m_conflict_var = UINT_MAX; }
+
         /**
          * push / pop are used only in self-contained mode from check_sat.
          */
@@ -198,33 +209,33 @@ namespace polysat {
         /**
          * Create polynomial terms
          */
-        poly var(unsigned v);
+        pdd var(unsigned v) { return m_vars[v]; }
 
-        /**
-         * deconstruct polynomials into sum of monomials.
-         */
-        vector<mono> poly2monos(poly const& p) const;        
 
         /**
          * Add polynomial constraints
          * Each constraint is tracked by a dependency.
          * assign sets the 'index'th bit of var.
          */
-        void add_eq(poly const& p, unsigned dep);
-        void add_diseq(poly const& p, unsigned dep);
-        void add_ule(poly const& p, poly const& q, unsigned dep);
-        void add_sle(poly const& p, poly const& q, unsigned dep);
+        void add_eq(pdd const& p, unsigned dep);
+        void add_diseq(pdd const& p, unsigned dep);
+        void add_ule(pdd const& p, pdd const& q, unsigned dep);
+        void add_sle(pdd const& p, pdd const& q, unsigned dep);
         void assign(unsigned var, unsigned index, bool value, unsigned dep);        
 
         /**
          * main state transitions.
          */
-        bool  can_propagate();
-        lbool propagate();
+        bool can_propagate();
+        void propagate();
 
         bool can_decide();
         void decide(rational & val, unsigned& var);
-        void assign(unsigned var, rational const& val);
+        
+        /**
+         * external decision
+         */
+        void assign(unsigned var, rational const& val) { inc_level(); assign(var, val, justification::decision(m_level)); }
         
         bool is_conflict();
         /**
diff --git a/src/tactic/core/elim_uncnstr_tactic.cpp b/src/tactic/core/elim_uncnstr_tactic.cpp
index ab0487a83..35e3e6cc4 100644
--- a/src/tactic/core/elim_uncnstr_tactic.cpp
+++ b/src/tactic/core/elim_uncnstr_tactic.cpp
@@ -470,7 +470,7 @@ class elim_uncnstr_tactic : public tactic {
             if (num == 2 && 
                 uncnstr(args[1]) && 
                 m_bv_util.is_numeral(args[0], val, bv_size) &&
-                m_bv_util.mult_inverse(val, bv_size, inv)) {
+                val.mult_inverse(bv_size, inv)) {
                 app * r;
                 if (!mk_fresh_uncnstr_var_for(f, num, args, r))
                     return r;
diff --git a/src/tactic/core/reduce_invertible_tactic.cpp b/src/tactic/core/reduce_invertible_tactic.cpp
index 06fde2757..cfee5bde7 100644
--- a/src/tactic/core/reduce_invertible_tactic.cpp
+++ b/src/tactic/core/reduce_invertible_tactic.cpp
@@ -342,7 +342,7 @@ private:
                 ensure_mc(mc);
                 expr_ref def(m);
                 rational inv_r;
-                VERIFY(m_bv.mult_inverse(r, sz, inv_r));
+                VERIFY(r.mult_inverse(sz, inv_r));
                 def = m_bv.mk_bv_mul(m_bv.mk_numeral(inv_r, sz), v);
                 (*mc)->add(v, def);
                 TRACE("invertible_tactic", tout << def << "\n";);
diff --git a/src/test/pdd.cpp b/src/test/pdd.cpp
index 216b3d48c..b10abb891 100644
--- a/src/test/pdd.cpp
+++ b/src/test/pdd.cpp
@@ -284,6 +284,44 @@ public :
         
     }
 
+    /**
+     * Test polynomials mod 2^2
+     */
+    static void mod4_operations() {
+        std::cout << "operations mod4\n";
+        pdd_manager m(4, pdd_manager::mod2N_e, 2);
+        unsigned va = 0;
+        unsigned vb = 1;
+        unsigned vc = 2;
+        unsigned vd = 3;
+        pdd a = m.mk_var(va);
+        pdd b = m.mk_var(vb);
+        pdd c = m.mk_var(vc);
+        pdd d = m.mk_var(vd);
+        pdd p = a - b;
+        std::cout << p << "\n";
+
+        // should be reduced to 0:
+        p = a*a*(a*a - 1);
+        std::cout << p << "\n"; 
+        vector<std::pair<unsigned, rational>> sub0, sub1, sub2, sub3;
+        sub0.push_back(std::make_pair(va, rational(0)));
+        sub1.push_back(std::make_pair(va, rational(1)));
+        sub2.push_back(std::make_pair(va, rational(2)));
+        sub3.push_back(std::make_pair(va, rational(3)));
+        std::cout << "sub 0 " << p.subst_val(sub0) << "\n";
+        std::cout << "sub 1 " << p.subst_val(sub1) << "\n";
+        std::cout << "sub 2 " << p.subst_val(sub2) << "\n";
+        std::cout << "sub 3 " << p.subst_val(sub3) << "\n";
+
+        std::cout << "expect 1 " << (2*a + 1).is_non_zero() << "\n";
+        std::cout << "expect 1 " << (2*a*b + 2*b + 1).is_non_zero() << "\n";
+        std::cout << "expect 0 " << (2*a + 2).is_non_zero() << "\n";
+        SASSERT((2*a + 1).is_non_zero());
+        SASSERT((2*a + 2*b + 1).is_non_zero());
+        SASSERT(!(2*a*b + 3*b + 2).is_non_zero());
+    }
+
 };
 
 }
@@ -297,4 +335,5 @@ void tst_pdd() {
     dd::test::iterator();
     dd::test::order();
     dd::test::order_lm();
+    dd::test::mod4_operations();
 }
diff --git a/src/test/rational.cpp b/src/test/rational.cpp
index ac477a02f..711958de5 100644
--- a/src/test/rational.cpp
+++ b/src/test/rational.cpp
@@ -453,6 +453,19 @@ static void tst11(bool use_ints) {
     std::cout << "\n";
 }
 
+static void tst12() {
+    std::cout << "test12\n";
+    rational r;
+    r = 5;
+    SASSERT(r.get_bit(0));
+    SASSERT(!r.get_bit(1));
+    SASSERT(r.get_bit(2));
+    SASSERT(!r.get_bit(3));
+    r = rational("10000000000000000000000000000000001");
+    for (unsigned i = 0; i < r.get_num_bits(); ++i)
+        std::cout << i << ": " << r.get_bit(i) << "\n";
+}
+
 
 void tst_rational() {
     TRACE("rational", tout << "starting rational test...\n";);
@@ -478,4 +491,5 @@ void tst_rational() {
     tst11(true);
     tst10(true);
     tst10(false);
+    tst12();
 }
diff --git a/src/util/mpq.h b/src/util/mpq.h
index 22147af63..1b97b6ece 100644
--- a/src/util/mpq.h
+++ b/src/util/mpq.h
@@ -800,6 +800,8 @@ public:
     unsigned storage_size(mpz const & a) { return mpz_manager<SYNCH>::size_info(a); }
     unsigned storage_size(mpq const & a) { return mpz_manager<SYNCH>::size_info(a.m_num) + mpz_manager<SYNCH>::size_info(a.m_den); }
 
+    bool get_bit(mpq const& a, unsigned index) { SASSERT(is_int(a) && !is_neg(a)); return mpz_manager<SYNCH>::get_bit(a.m_num, index); }
+
     /**
        \brief Return true if the number is a perfect square, and 
        store the square root in 'root'.
diff --git a/src/util/mpz.cpp b/src/util/mpz.cpp
index 8e08560d7..5e426eb18 100644
--- a/src/util/mpz.cpp
+++ b/src/util/mpz.cpp
@@ -2469,6 +2469,33 @@ bool mpz_manager<SYNCH>::decompose(mpz const & a, svector<digit_t> & digits) {
     }
 }
 
+template<bool SYNCH>
+bool mpz_manager<SYNCH>::get_bit(mpz const & a, unsigned index) {
+    if (is_small(a)) {
+        SASSERT(a.m_val >= 0);
+        if (index >= 8*sizeof(digit_t))
+            return false;
+        return 0 != (a.m_val & (1ull << (digit_t)index));
+    }
+    unsigned i = index / (sizeof(digit_t)*8);
+    unsigned o = index % (sizeof(digit_t)*8);
+
+#ifndef _MP_GMP
+    mpz_cell * cell_a = a.m_ptr;
+    unsigned sz = cell_a->m_size;
+    if (sz*sizeof(digit_t)*8 <= index)
+        return false;
+    return 0 != (cell_a->m_digits[i] & (1ull << (digit_t)o));
+#else
+    SASSERT(!is_neg(a));
+    svector<digit_t> digits;
+    decompose(a, digits);
+    if (digits.size()*sizeof(digit_t)*8 <= index)
+        return false;
+    return 0 != (digits[i] & (1ull << (digit_t)o));    
+#endif
+}
+
 template<bool SYNCH>
 bool mpz_manager<SYNCH>::divides(mpz const & a, mpz const & b) {
     _scoped_numeral<mpz_manager<SYNCH> > tmp(*this);
diff --git a/src/util/mpz.h b/src/util/mpz.h
index ebfe98704..f3c3b19b9 100644
--- a/src/util/mpz.h
+++ b/src/util/mpz.h
@@ -621,6 +621,7 @@ public:
     */
     void display_bin(std::ostream & out, mpz const & a, unsigned num_bits) const;
 
+
     static unsigned hash(mpz const & a);
 
     static bool is_one(mpz const & a) {
@@ -719,6 +720,8 @@ public:
     // Store the digits of n into digits, and return the sign.
     bool decompose(mpz const & n, svector<digit_t> & digits);
 
+    bool get_bit(mpz const& a, unsigned bit);
+
 };
 
 #ifndef SINGLE_THREAD
diff --git a/src/util/rational.cpp b/src/util/rational.cpp
index dcbf09969..bb443935a 100644
--- a/src/util/rational.cpp
+++ b/src/util/rational.cpp
@@ -130,3 +130,26 @@ bool rational::limit_denominator(rational &num, rational const& limit) {
     return false;
 }
 
+bool rational::mult_inverse(unsigned num_bits, rational & result) {
+    rational const& n = *this;
+    if (n.is_one()) {
+        result = n;
+        return true;
+    }
+    
+    if (n.is_even())
+        return false;
+    
+    rational g;
+    rational x;
+    rational y;
+    g = gcd(n, rational::power_of_two(num_bits), x, y);
+    if (x.is_neg()) {
+        x = mod(x, rational::power_of_two(num_bits));
+    }
+    SASSERT(x.is_pos());
+    SASSERT(mod(x * n, rational::power_of_two(num_bits)).is_one());
+    result = x;
+    return true;
+}
+
diff --git a/src/util/rational.h b/src/util/rational.h
index 6585fe5a7..3b8ee5649 100644
--- a/src/util/rational.h
+++ b/src/util/rational.h
@@ -308,6 +308,10 @@ public:
     bool is_even() const {
         return m().is_even(m_val);
     }
+
+    bool is_odd() const { 
+        return !is_even(); 
+    }
     
     friend inline rational floor(rational const & r) {
         rational f;
@@ -333,6 +337,8 @@ public:
         return m().is_power_of_two(m_val, shift);
     }
 
+    bool mult_inverse(unsigned num_bits, rational & result);
+
     static rational const & zero() {
         return m_zero;
     }
@@ -459,6 +465,18 @@ public:
         return get_num_digits(rational(10));
     }
 
+    bool get_bit(unsigned index) const {
+        return m().get_bit(m_val, index);
+    }
+
+    unsigned trailing_zeros() const {
+        if (is_zero())
+            return 0;
+        unsigned k = 0;
+        for (; !get_bit(k); ++k); 
+        return k;
+    }
+
     static bool limit_denominator(rational &num, rational const& limit);
 };