diff --git a/src/sat/smt/intblast_solver.cpp b/src/sat/smt/intblast_solver.cpp
index 5c12d4061..da3526bb4 100644
--- a/src/sat/smt/intblast_solver.cpp
+++ b/src/sat/smt/intblast_solver.cpp
@@ -71,8 +71,6 @@ namespace intblast {
         return n && n->is_attached_to(get_id());
     }
 
-
-
     bool solver::post_visit(expr* e, bool sign, bool root) {
         euf::enode* n = expr2enode(e);
         app* a = to_app(e);
@@ -236,7 +234,7 @@ namespace intblast {
 
         m_core.reset();
         m_vars.reset();
-        m_trail.reset();
+        m_translate.reset();
         m_solver = mk_smt2_solver(m, s.params(), symbol::null);
 
         expr_ref_vector es(m);
@@ -275,7 +273,6 @@ namespace intblast {
                     m_core.push_back(ctx.mk_literal(e));
             }
         }
-
         return r;
     };
 
@@ -446,7 +443,6 @@ namespace intblast {
                 sorts.push_back(a.mk_int());
             }
             else
-
                 sorts.push_back(s);
         }
         b = translated(b);
@@ -855,7 +851,6 @@ namespace intblast {
         }
         set_translated(e, r);
     }
-
     void solver::translate_basic(app* e) {
         if (m.is_eq(e)) {
             bool has_bv_arg = any_of(*e, [&](expr* arg) { return bv.is_bv(arg); });
diff --git a/src/sat/smt/intblast_solver.h b/src/sat/smt/intblast_solver.h
index f4eb24b92..35790c6aa 100644
--- a/src/sat/smt/intblast_solver.h
+++ b/src/sat/smt/intblast_solver.h
@@ -8,7 +8,8 @@ Module Name:
 Abstract:
 
     Int-blast solver.
-    It assumes a full assignemnt to literals in 
+
+    check_solver_state assumes a full assignment to literals in 
     irredundant clauses. 
     It picks a satisfying Boolean assignment and 
     checks if it is feasible for bit-vectors using
@@ -45,6 +46,14 @@ namespace euf {
 namespace intblast {
 
     class solver : public euf::th_euf_solver {
+<<<<<<< HEAD
+=======
+        struct var_info {
+            expr* dst;
+            rational sz;
+        };
+
+>>>>>>> 4cadf6d9f (preparing intblaster as self-contained solver.)
         euf::solver& ctx;
         sat::solver& s;
         ast_manager& m;
@@ -94,6 +103,27 @@ namespace intblast {
         void add_value_plugin(euf::enode* n, model& mdl, expr_ref_vector& values);
         void add_value_solver(euf::enode* n, model& mdl, expr_ref_vector& values);
 
+        expr* translated(expr* e) { expr* r = m_translate.get(e->get_id(), nullptr); SASSERT(r); return r; }
+        void set_translated(expr* e, expr* r) { m_translate.setx(e->get_id(), r); }
+        expr* arg(unsigned i) { return m_args.get(i); }
+
+        expr* mk_mod(expr* x);
+        expr* mk_smod(expr* x);
+        expr* bv_expr = nullptr;
+        rational bv_size();
+
+        void translate_expr(expr* e);
+        void translate_bv(app* e);
+        void translate_basic(app* e);
+        void translate_app(app* e);
+        void translate_quantifier(quantifier* q);
+        void translate_var(var* v);
+
+        void ensure_args(app* e);
+        void internalize_bv(app* e);
+
+        euf::theory_var mk_var(euf::enode* n) override;
+
     public:
         solver(euf::solver& ctx);
         
@@ -137,7 +167,7 @@ namespace intblast {
 
         sat::literal_vector const& unsat_core();
 
-        void add_value(euf::enode* n, model& mdl, expr_ref_vector& values);
+        void add_value(euf::enode* n, model& mdl, expr_ref_vector& values) override;
 
 
     };
diff --git a/src/sat/smt/polysat/constraints.h b/src/sat/smt/polysat/constraints.h
index 15d8dfa09..a9ec63165 100644
--- a/src/sat/smt/polysat/constraints.h
+++ b/src/sat/smt/polysat/constraints.h
@@ -41,6 +41,8 @@ namespace polysat {
         virtual std::ostream& display(std::ostream& out) const = 0;
         virtual lbool eval() const = 0;
         virtual lbool eval(assignment const& a) const = 0;
+        virtual void activate(core& c, bool sign, dependency const& d) = 0;
+        virtual void propagate(core& c, lbool value, dependency const& d) = 0;
     };
 
     inline std::ostream& operator<<(std::ostream& out, constraint const& c) { return c.display(out); }
@@ -61,6 +63,8 @@ namespace polysat {
         unsigned_vector const& vars() const { return m_constraint->vars(); }
         unsigned var(unsigned idx) const { return m_constraint->var(idx); }
         bool contains_var(pvar v) const { return m_constraint->contains_var(v); }
+        void activate(core& c, dependency const& d) { m_constraint->activate(c, m_sign, d); }
+        void propagate(core& c, lbool value, dependency const& d) { m_constraint->propagate(c, value, d); }
         bool is_always_true() const { return eval() == l_true; }
         bool is_always_false() const { return eval() == l_false; }
         lbool eval(assignment& a) const;
@@ -84,6 +88,8 @@ namespace polysat {
 
         signed_constraint eq(pdd const& p) { return ule(p, p.manager().mk_val(0)); }
         signed_constraint eq(pdd const& p, rational const& v) { return eq(p - p.manager().mk_val(v)); }
+        signed_constraint eq(pdd const& p, unsigned v) { return eq(p - p.manager().mk_val(v)); }
+        signed_constraint eq(pdd const& p, pdd const& q) { return eq(p - q); }
         signed_constraint ule(pdd const& p, pdd const& q);
         signed_constraint sle(pdd const& p, pdd const& q) { auto sh = rational::power_of_two(p.power_of_2() - 1); return ule(p + sh, q + sh); }
         signed_constraint ult(pdd const& p, pdd const& q) { return ~ule(q, p); }
diff --git a/src/sat/smt/polysat/core.cpp b/src/sat/smt/polysat/core.cpp
index a552bb9ab..dd30e8227 100644
--- a/src/sat/smt/polysat/core.cpp
+++ b/src/sat/smt/polysat/core.cpp
@@ -187,12 +187,13 @@ namespace polysat {
         return sc;
     }
 
-
     void core::propagate_assignment(prop_item& dc) { 
         auto [idx, sign, dep] = dc;
         auto sc = get_constraint(idx, sign);
         if (sc.is_eq(m_var, m_value))
             propagate_assignment(m_var, m_value, dep);
+        else 
+            sc.activate(*this, dep);        
     }
 
     void core::add_watch(unsigned idx, unsigned var) {
@@ -216,7 +217,7 @@ namespace polysat {
         // for entries where there is only one free variable left add to viable set 
         unsigned j = 0;
         for (auto idx : m_watch[v]) {
-            auto [sc, as, value] = m_constraint_index[idx];
+            auto [sc, dep, value] = m_constraint_index[idx];
             auto& vars = sc.vars();
             if (vars[0] != v)
                 std::swap(vars[0], vars[1]);
@@ -231,6 +232,8 @@ namespace polysat {
                 }
             }
 
+            sc.propagate(*this, value, dep);
+
             SASSERT(!swapped || vars.size() <= 1 || (!is_assigned(vars[0]) && !is_assigned(vars[1])));
             if (swapped)
                 continue;
@@ -262,6 +265,11 @@ namespace polysat {
         default:
             break;
         }
+        // propagate current assignment for sc
+        sc.propagate(*this, to_lbool(!sign), dep);
+        if (s.inconsistent())
+            return;
+
         // if sc is v == value, then check the watch list for v to propagate truth assignments
         if (sc.is_eq(m_var, m_value)) {
             for (auto idx1 : m_watch[m_var]) {
@@ -360,4 +368,13 @@ namespace polysat {
 
     }
 
+    void core::add_axiom(signed_constraint sc) {
+        auto idx = register_constraint(sc, dependency::axiom());
+        assign_eh(idx, false, dependency::axiom());
+    }
+
+    void core::add_clause(char const* name, core_vector const& cs, bool is_redundant) {
+        s.add_polysat_clause(name, cs, is_redundant);
+    }
+
 }
diff --git a/src/sat/smt/polysat/core.h b/src/sat/smt/polysat/core.h
index c3dddfece..6297e567e 100644
--- a/src/sat/smt/polysat/core.h
+++ b/src/sat/smt/polysat/core.h
@@ -84,7 +84,7 @@ namespace polysat {
         void get_bitvector_prefixes(pvar v, pvar_vector& out);
         void get_fixed_bits(pvar v, svector<justified_fixed_bits>& fixed_bits);
         bool inconsistent() const;
-        void add_clause(char const* name, std::initializer_list<signed_constraint> cs, bool is_redundant);
+
 
 
         void add_watch(unsigned idx, unsigned var);
@@ -94,6 +94,8 @@ namespace polysat {
         lbool eval(signed_constraint const& sc);
         dependency_vector explain_eval(signed_constraint const& sc);
 
+        void add_axiom(signed_constraint sc);
+
     public:
         core(solver_interface& s);
 
@@ -118,13 +120,20 @@ namespace polysat {
         signed_constraint bit(pdd const& p, unsigned i) { return m_constraints.bit(p, i); }
 
 
-        signed_constraint lshr(pdd const& a, pdd const& b, pdd const& r) { return m_constraints.lshr(a, b, r); }
-        signed_constraint ashr(pdd const& a, pdd const& b, pdd const& r) { return m_constraints.ashr(a, b, r); }
-        signed_constraint shl(pdd const& a, pdd const& b, pdd const& r) { return m_constraints.shl(a, b, r); }
-        signed_constraint band(pdd const& a, pdd const& b, pdd const& r) { return m_constraints.band(a, b, r); }
+        void lshr(pdd const& a, pdd const& b, pdd const& r) { add_axiom(m_constraints.lshr(a, b, r)); }
+        void ashr(pdd const& a, pdd const& b, pdd const& r) { add_axiom(m_constraints.ashr(a, b, r)); }
+        void shl(pdd const& a, pdd const& b, pdd const& r) { add_axiom(m_constraints.shl(a, b, r)); }
+        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; }
 
+
+        /*
+        * Add a named clause. Dependencies are assumed, signed constraints are guaranteeed.
+        * In other words, the clause represents the formula /\ d_i -> \/ sc_j
+        * Where d_i are logical interpretations of dependencies and sc_j are signed constraints.
+        */
+        void add_clause(char const* name, core_vector const& cs, bool is_redundant);
         
         pvar add_var(unsigned sz);
         pdd var(pvar p) { return m_vars[p]; }
diff --git a/src/sat/smt/polysat/op_constraint.cpp b/src/sat/smt/polysat/op_constraint.cpp
index 401b1ca52..234ddea04 100644
--- a/src/sat/smt/polysat/op_constraint.cpp
+++ b/src/sat/smt/polysat/op_constraint.cpp
@@ -13,13 +13,13 @@ Notes:
 
 Additional possible functionality on constraints:
 
-- activate - when operation is first activated. It may be created and only activated later.
 - bit-wise assignments - narrow based on bit assignment, not entire word assignment.
 - integration with congruence tables
 - integration with conflict resolution
 
 --*/
 
+#include "util/log.h"
 #include "sat/smt/polysat/op_constraint.h"
 #include "sat/smt/polysat/core.h"
 
@@ -157,7 +157,6 @@ namespace polysat {
             return out << "&";
         case op_constraint::code::inv_op:
             return out << "inv";
-
         default:
             UNREACHABLE();
             return out;
@@ -176,96 +175,95 @@ namespace polysat {
         return out << r() << " " << eq << " " << p() << " " << m_op << " " << q();
     }
 
-#if 0
-    /**
-     * Produce lemmas that contradict the given assignment.
-     *
-     * We can assume that op_constraint is only asserted positive.
-     */
-    clause_ref op_constraint::produce_lemma(solver& s, assignment const& a, bool is_positive) {
-        SASSERT(is_positive);
-
-        if (is_currently_true(a, is_positive))
-            return {};
-
-        return produce_lemma(s, a);
-    }
-
-    clause_ref op_constraint::produce_lemma(solver& s, assignment const& a) {
+    void op_constraint::activate(core& c, bool sign, dependency const& dep) {
+        SASSERT(!sign);
         switch (m_op) {
-        case code::lshr_op:
-            return lemma_lshr(s, a);
-        case code::shl_op:
-            return lemma_shl(s, a);
         case code::and_op:
-            return lemma_and(s, a);
-        case code::inv_op:
-            return lemma_inv(s, a);
+            activate_and(c, dep);
+            break;
         default:
-            NOT_IMPLEMENTED_YET();
-            return {};
+            break;
         }
     }
 
+    void op_constraint::propagate(core& c, lbool value, dependency const& dep) {
+        SASSERT(value == l_true);
+        switch (m_op) {
+        case code::lshr_op:
+            propagate_lshr(c, dep);
+            break;
+        case code::shl_op:
+            propagate_shl(c, dep);
+            break;
+        case code::and_op:
+            propagate_and(c, dep);
+            break;
+        case code::inv_op:
+            propagate_inv(c, dep);
+            break;
+        default:
+            NOT_IMPLEMENTED_YET();
+            break;
+        }
+    }
+
+    void op_constraint::propagate_inv(core& s, dependency const& dep) {
+
+    }
+
     /**
-     * Enforce basic axioms for r == p >> q:
-     *
-     *      q >= N  ->  r = 0
-     *      q >= k  ->  r[i] = 0 for N - k <= i < N     (bit indices range from 0 to N-1, inclusive)
-     *      q >= k  ->  r <= 2^{N-k} - 1
-     *      q = k   ->  r[i] = p[i+k] for 0 <= i < N - k
-     *      r <= p
-     *      q != 0  ->  r <= p      (subsumed by previous axiom)
-     *      q != 0  /\  p > 0 ->  r < p
-     *      q = 0   ->  r = p
-     *      p = q   ->  r = 0
-     *
-     * when q is a constant, several axioms can be enforced at activation time.
-     *
-     * Enforce also inferences and bounds
-     *
-     * TODO: use also
-     * s.m_viable.min_viable();
-     * s.m_viable.max_viable()
-     * when r, q are variables.
-     */
-    clause_ref op_constraint::lemma_lshr(solver& s, assignment const& a) {
+  * Enforce basic axioms for r == p >> q:
+  *
+  *      q >= N  ->  r = 0
+  *      q >= k  ->  r[i] = 0 for N - k <= i < N     (bit indices range from 0 to N-1, inclusive)
+  *      q >= k  ->  r <= 2^{N-k} - 1
+  *      q = k   ->  r[i] = p[i+k] for 0 <= i < N - k
+  *      r <= p
+  *      q != 0  ->  r <= p      (subsumed by previous axiom)
+  *      q != 0  /\  p > 0 ->  r < p
+  *      q = 0   ->  r = p
+  *      p = q   ->  r = 0
+  *
+  * when q is a constant, several axioms can be enforced at activation time.
+  *
+  * Enforce also inferences and bounds
+  *
+  * TODO: use also
+  * s.m_viable.min_viable();
+  * s.m_viable.max_viable()
+  * when r, q are variables.
+  */
+    void op_constraint::propagate_lshr(core& c, dependency const& d) {
         auto& m = p().manager();
-        auto const pv = a.apply_to(p());
-        auto const qv = a.apply_to(q());
-        auto const rv = a.apply_to(r());
+        auto const pv = c.subst(p());
+        auto const qv = c.subst(q());
+        auto const rv = c.subst(r());
         unsigned const N = m.power_of_2();
 
-        signed_constraint const lshr(this, true);
+        
+        signed_constraint const lshr(polysat::ckind_t::op_t, this);
+        auto& C = c.cs();
 
         if (pv.is_val() && rv.is_val() && rv.val() > pv.val())
-            // r <= p
-            return s.mk_clause(~lshr, s.ule(r(), p()), true);
+            c.add_clause("lshr 1", { d, C.ule(r(), p()) }, false);
+
         else if (qv.is_val() && qv.val() >= N && rv.is_val() && !rv.is_zero())
             // TODO: instead of rv.is_val() && !rv.is_zero(), we should use !is_forced_zero(r) which checks whether eval(r) = 0 or bvalue(r=0) = true; see saturation.cpp
-            // q >= N -> r = 0
-            return s.mk_clause(~lshr, ~s.ule(N, q()), s.eq(r()), true);
+            c.add_clause("q >= N -> r = 0", { d, ~C.ule(N, q()), C.eq(r()) }, true);
         else if (qv.is_zero() && pv.is_val() && rv.is_val() && pv != rv)
-            // q = 0 -> p = r
-            return s.mk_clause(~lshr, ~s.eq(q()), s.eq(p(), r()), true);
+            c.add_clause("q = 0 -> p = r", { d, ~C.eq(q()), C.eq(p(), r()) } , true);
         else if (qv.is_val() && !qv.is_zero() && pv.is_val() && rv.is_val() && !pv.is_zero() && rv.val() >= pv.val())
-            // q != 0 & p > 0 -> r < p
-            return s.mk_clause(~lshr, s.eq(q()), s.ule(p(), 0), s.ult(r(), p()), true);
+            c.add_clause("q != 0 & p > 0 -> r < p", { d, C.eq(q()), C.ule(p(), 0), C.ult(r(), p()) }, true);
         else if (qv.is_val() && !qv.is_zero() && qv.val() < N && rv.is_val() && rv.val() > rational::power_of_two(N - qv.val().get_unsigned()) - 1)
-            // q >= k -> r <= 2^{N-k} - 1
-            return s.mk_clause(~lshr, ~s.ule(qv.val(), q()), s.ule(r(), rational::power_of_two(N - qv.val().get_unsigned()) - 1), true);
-        // else if (pv == qv && !rv.is_zero())
-        //     return s.mk_clause(~lshr, ~s.eq(p(), q()), s.eq(r()), true);
+            c.add_clause("q >= k -> r <= 2^{N-k} - 1", { d, ~C.ule(qv.val(), q()), C.ule(r(), rational::power_of_two(N - qv.val().get_unsigned()) - 1)}, true);
         else if (pv.is_val() && rv.is_val() && qv.is_val() && !qv.is_zero()) {
             unsigned k = qv.val().get_unsigned();
-            // q = k  ->  r[i] = p[i+k] for 0 <= i < N - k
             for (unsigned i = 0; i < N - k; ++i) {
-                if (rv.val().get_bit(i) && !pv.val().get_bit(i + k)) {
-                    return s.mk_clause(~lshr, ~s.eq(q(), k), ~s.bit(r(), i), s.bit(p(), i + k), true);
-                }
-                if (!rv.val().get_bit(i) && pv.val().get_bit(i + k)) {
-                    return s.mk_clause(~lshr, ~s.eq(q(), k), s.bit(r(), i), ~s.bit(p(), i + k), true);
-                }
+                if (rv.val().get_bit(i) && !pv.val().get_bit(i + k)) 
+                    c.add_clause("q = k  ->  r[i] = p[i+k] for 0 <= i < N - k", { d, ~C.eq(q(), k), ~C.bit(r(), i), C.bit(p(), i + k) }, true);
+                
+                if (!rv.val().get_bit(i) && pv.val().get_bit(i + k)) 
+                    c.add_clause("q = k  ->  r[i] = p[i+k] for 0 <= i < N - k", { d, ~C.eq(q(), k), C.bit(r(), i), ~C.bit(p(), i + k) }, true);                
             }
         }
         else {
@@ -276,19 +274,44 @@ namespace polysat {
                 rational const& q_val = qv.val();
                 if (q_val >= N)
                     // q >= N ==> r = 0
-                    return s.mk_clause(~lshr, ~s.ule(N, q()), s.eq(r()), true);
-                if (pv.is_val()) {
+                    c.add_clause("q >= N ==> r = 0", { d, ~C.ule(N, q()), C.eq(r()) }, true);
+                else if (pv.is_val()) {
                     SASSERT(q_val.is_unsigned());
-                    // p = p_val & q = q_val ==> r = p_val / 2^q_val
+                    // 
                     rational const r_val = machine_div2k(pv.val(), q_val.get_unsigned());
-                    return s.mk_clause(~lshr, ~s.eq(p(), pv), ~s.eq(q(), qv), s.eq(r(), r_val), true);
+                    c.add_clause("p = p_val & q = q_val ==> r = p_val / 2^q_val", { d, ~C.eq(p(), pv), ~C.eq(q(), qv), C.eq(r(), r_val) }, true);
                 }
             }
         }
-        return {};
     }
 
 
+    void op_constraint::activate_and(core& c, dependency const& d) {
+        auto x = p(), y = q();
+        auto& C = c.cs();
+        if (x.is_val())
+            std::swap(x, y);
+        if (!y.is_val())
+            return;
+        auto& m = x.manager();
+        auto yv = y.val();
+        if (!(yv + 1).is_power_of_two())
+            return;
+        if (yv == m.max_value())
+            c.add_clause("band-mask-true", { d, C.eq(x, r()) }, false);
+        else if (yv == 0)
+            c.add_clause("band-mask-false", { d, C.eq(r()) }, false);
+        else {
+            unsigned N = m.power_of_2();
+            unsigned k = yv.get_num_bits();
+            SASSERT(k < N);
+            rational exp = rational::power_of_two(N - k);
+            c.add_clause("band-mask 1", { d, C.eq(x * exp, r() * exp) }, false);
+            c.add_clause("band-mask 2", { d, C.ule(r(), y) }, false);  // maybe always activate these constraints regardless?
+        }
+    }   
+
+
     /**
      * Enforce axioms for constraint: r == p << q
      *
@@ -298,35 +321,33 @@ namespace polysat {
      *      q = k   ->  r[i+k] = p[i] for 0 <= i < N - k
      *      q = 0   ->  r = p
      */
-    clause_ref op_constraint::lemma_shl(solver& s, assignment const& a) {
+    void op_constraint::propagate_shl(core& c, dependency const& d) {
         auto& m = p().manager();
-        auto const pv = a.apply_to(p());
-        auto const qv = a.apply_to(q());
-        auto const rv = a.apply_to(r());
+        auto const pv = c.subst(p());
+        auto const qv = c.subst(q());
+        auto const rv = c.subst(r());
         unsigned const N = m.power_of_2();
+        auto& C = c.cs();
 
-        signed_constraint const shl(this, true);
-
-        if (qv.is_val() && qv.val() >= N && rv.is_val() && !rv.is_zero())
-            // q >= N  ->  r = 0
-            return s.mk_clause(~shl, ~s.ule(N, q()), s.eq(r()), true);
+        if (qv.is_val() && qv.val() >= N && rv.is_val() && !rv.is_zero())            
+            c.add_clause("q >= N  ->  r = 0", { d, ~C.ule(N, q()), C.eq(r()) }, true);
         else if (qv.is_zero() && pv.is_val() && rv.is_val() && rv != pv)
-            // q = 0  ->  r = p
-            return s.mk_clause(~shl, ~s.eq(q()), s.eq(r(), p()), true);
+            // 
+            c.add_clause("q = 0  ->  r = p", { d, ~C.eq(q()), C.eq(r(), p()) }, true);
         else if (qv.is_val() && !qv.is_zero() && qv.val() < N && rv.is_val() &&
             !rv.is_zero() && rv.val() < rational::power_of_two(qv.val().get_unsigned()))
             // q >= k  ->  r = 0  \/  r >= 2^k  (intuitive version)
             // q >= k  ->  r - 1 >= 2^k - 1     (equivalent unit constraint to better support narrowing)
-            return s.mk_clause(~shl, ~s.ule(qv.val(), q()), s.ule(rational::power_of_two(qv.val().get_unsigned()) - 1, r() - 1), true);
+            c.add_clause("q >= k  ->  r - 1 >= 2^k - 1", { d, ~C.ule(qv.val(), q()), C.ule(rational::power_of_two(qv.val().get_unsigned()) - 1, r() - 1) }, true);
         else if (pv.is_val() && rv.is_val() && qv.is_val() && !qv.is_zero()) {
             unsigned k = qv.val().get_unsigned();
             // q = k  ->  r[i+k] = p[i] for 0 <= i < N - k
             for (unsigned i = 0; i < N - k; ++i) {
                 if (rv.val().get_bit(i + k) && !pv.val().get_bit(i)) {
-                    return s.mk_clause(~shl, ~s.eq(q(), k), ~s.bit(r(), i + k), s.bit(p(), i), true);
+                    c.add_clause("q = k  ->  r[i+k] = p[i] for 0 <= i < N - k", { d, ~C.eq(q(), k), ~C.bit(r(), i + k), C.bit(p(), i) }, true);
                 }
                 if (!rv.val().get_bit(i + k) && pv.val().get_bit(i)) {
-                    return s.mk_clause(~shl, ~s.eq(q(), k), s.bit(r(), i + k), ~s.bit(p(), i), true);
+                    c.add_clause("q = k  ->  r[i+k] = p[i] for 0 <= i < N - k", { d, ~C.eq(q(), k), C.bit(r(), i + k), ~C.bit(p(), i) }, true);
                 }
             }
         }
@@ -338,43 +359,15 @@ namespace polysat {
                 rational const& q_val = qv.val();
                 if (q_val >= N)
                     // q >= N ==> r = 0
-                    return s.mk_clause("shl forward 1", {~shl, ~s.ule(N, q()), s.eq(r())}, true);
+                    c.add_clause("shl forward 1", {d, ~C.ule(N, q()), C.eq(r())}, true);
                 if (pv.is_val()) {
                     SASSERT(q_val.is_unsigned());
                     // p = p_val & q = q_val ==> r = p_val * 2^q_val
                     rational const r_val = pv.val() * rational::power_of_two(q_val.get_unsigned());
-                    return s.mk_clause("shl forward 2", {~shl, ~s.eq(p(), pv), ~s.eq(q(), qv), s.eq(r(), r_val)}, true);
+                    c.add_clause("shl forward 2", {d, ~C.eq(p(), pv), ~C.eq(q(), qv), C.eq(r(), r_val)}, true);
                 }
             }
         }
-        return {};
-    }
-
-
-
-    void op_constraint::activate_and(solver& s) {
-        auto x = p(), y = q();
-        if (x.is_val())
-            std::swap(x, y);
-        if (!y.is_val())
-            return;
-        auto& m = x.manager();
-        auto yv = y.val();
-        if (!(yv + 1).is_power_of_two())
-            return;
-        signed_constraint const andc(this, true);
-        if (yv == m.max_value())
-            s.add_clause(~andc, s.eq(x, r()), false);
-        else if (yv == 0)
-            s.add_clause(~andc, s.eq(r()), false);
-        else {
-            unsigned N = m.power_of_2();
-            unsigned k = yv.get_num_bits();
-            SASSERT(k < N);
-            rational exp = rational::power_of_two(N - k);
-            s.add_clause(~andc, s.eq(x * exp, r() * exp), false);
-            s.add_clause(~andc, s.ule(r(), y), false);  // maybe always activate these constraints regardless?
-        }
     }
 
     /**
@@ -390,48 +383,39 @@ namespace polysat {
      * p = 2^k - 1 && r = 0 && q != 0 => q >= 2^k
      * q = 2^k - 1 && r = 0 && p != 0 => p >= 2^k
      */
-    clause_ref op_constraint::lemma_and(solver& s, assignment const& a) {
+    void op_constraint::propagate_and(core& c, dependency const& d) {
         auto& m = p().manager();
-        auto pv = a.apply_to(p());
-        auto qv = a.apply_to(q());
-        auto rv = a.apply_to(r());
+        auto pv = c.subst(p());
+        auto qv = c.subst(q());
+        auto rv = c.subst(r());
+        auto& C = c.cs();
 
-        signed_constraint const andc(this, true); // op_constraints are always true
-
-        // r <= p
         if (pv.is_val() && rv.is_val() && rv.val() > pv.val())
-            return s.mk_clause(~andc, s.ule(r(), p()), true);
-        // r <= q
-        if (qv.is_val() && rv.is_val() && rv.val() > qv.val())
-            return s.mk_clause(~andc, s.ule(r(), q()), true);
-        // p = q => r = p
-        if (pv.is_val() && qv.is_val() && rv.is_val() && pv == qv && rv != pv)
-            return s.mk_clause(~andc, ~s.eq(p(), q()), s.eq(r(), p()), true);
-        if (pv.is_val() && qv.is_val() && rv.is_val()) {
-            // p = -1 => r = q
+            c.add_clause("p&q <= p", { d, C.ule(r(), p()) }, true);
+        else if (qv.is_val() && rv.is_val() && rv.val() > qv.val())
+            c.add_clause("p&q <= q", { d, C.ule(r(), q()) }, true);
+        else if (pv.is_val() && qv.is_val() && rv.is_val() && pv == qv && rv != pv)
+            c.add_clause("p = q => r = p", { d, ~C.eq(p(), q()), C.eq(r(), p()) }, true);
+        else if (pv.is_val() && qv.is_val() && rv.is_val()) {
             if (pv.is_max() && qv != rv)
-                return s.mk_clause(~andc, ~s.eq(p(), m.max_value()), s.eq(q(), r()), true);
-            // q = -1 => r = p
+                c.add_clause("p = -1 => r = q", { d, ~C.eq(p(), m.max_value()), C.eq(q(), r()) }, true);
             if (qv.is_max() && pv != rv)
-                return s.mk_clause(~andc, ~s.eq(q(), m.max_value()), s.eq(p(), r()), true);
+                c.add_clause("q = -1 => r = p", { d, ~C.eq(q(), m.max_value()), C.eq(p(), r()) }, true);
 
             unsigned const N = m.power_of_2();
             unsigned pow;
             if ((pv.val() + 1).is_power_of_two(pow)) {
-                // p = 2^k - 1 && r = 0 && q != 0 => q >= 2^k
                 if (rv.is_zero() && !qv.is_zero() && qv.val() <= pv.val())
-                    return s.mk_clause(~andc, ~s.eq(p(), pv), ~s.eq(r()), s.eq(q()), s.ule(pv + 1, q()), true);
-                // p = 2^k - 1  ==>  r*2^{N - k} = q*2^{N - k}
+                    c.add_clause("p = 2^k - 1 && r = 0 && q != 0 => q >= 2^k", { d, ~C.eq(p(), pv), ~C.eq(r()), C.eq(q()), C.ule(pv + 1, q()) }, true);
                 if (rv != qv)
-                    return s.mk_clause(~andc, ~s.eq(p(), pv), s.eq(r() * rational::power_of_two(N - pow), q() * rational::power_of_two(N - pow)), true);
+                    c.add_clause("p = 2^k - 1  ==>  r*2^{N - k} = q*2^{N - k}", { d, ~C.eq(p(), pv), C.eq(r() * rational::power_of_two(N - pow), q() * rational::power_of_two(N - pow)) }, true);
             }
             if ((qv.val() + 1).is_power_of_two(pow)) {
-                // q = 2^k - 1 && r = 0 && p != 0  ==>  p >= 2^k
                 if (rv.is_zero() && !pv.is_zero() && pv.val() <= qv.val())
-                    return s.mk_clause(~andc, ~s.eq(q(), qv), ~s.eq(r()), s.eq(p()), s.ule(qv + 1, p()), true);
-                // q = 2^k - 1  ==>  r*2^{N - k} = p*2^{N - k}
+                    c.add_clause("q = 2^k - 1 && r = 0 && p != 0  ==>  p >= 2^k", { d, ~C.eq(q(), qv), ~C.eq(r()), C.eq(p()), C.ule(qv + 1, p()) }, true);
+                // 
                 if (rv != pv)
-                    return s.mk_clause(~andc, ~s.eq(q(), qv), s.eq(r() * rational::power_of_two(N - pow), p() * rational::power_of_two(N - pow)), true);
+                    c.add_clause("q = 2^k - 1  ==>  r*2^{N - k} = p*2^{N - k}", { d, ~C.eq(q(), qv), C.eq(r() * rational::power_of_two(N - pow), p() * rational::power_of_two(N - pow)) }, true);
             }
 
             for (unsigned i = 0; i < N; ++i) {
@@ -441,33 +425,31 @@ namespace polysat {
                 if (rb == (pb && qb))
                     continue;
                 if (pb && qb && !rb)
-                    return s.mk_clause(~andc, ~s.bit(p(), i), ~s.bit(q(), i), s.bit(r(), i), true);
+                    c.add_clause("p&q[i] = p[i]&q[i]", { d, ~C.bit(p(), i), ~C.bit(q(), i), C.bit(r(), i) }, true);
                 else if (!pb && rb)
-                    return s.mk_clause(~andc, s.bit(p(), i), ~s.bit(r(), i), true);
+                    c.add_clause("p&q[i] = p[i]&q[i]", { d, C.bit(p(), i), ~C.bit(r(), i) }, true);
                 else if (!qb && rb)
-                    return s.mk_clause(~andc, s.bit(q(), i), ~s.bit(r(), i), true);
+                    c.add_clause("p&q[i] = p[i]&q[i]", { d, C.bit(q(), i), ~C.bit(r(), i) }, true);
                 else
                     UNREACHABLE();
             }
-            return {};
+            return;
         }
 
         // Propagate r if p or q are 0
-        if (pv.is_zero() && !rv.is_zero())  // rv not necessarily fully evaluated
-            return s.mk_clause(~andc, s.ule(r(), p()), true);
-        if (qv.is_zero() && !rv.is_zero())  // rv not necessarily fully evaluated
-            return s.mk_clause(~andc, s.ule(r(), q()), true);
+        else if (pv.is_zero() && !rv.is_zero())  // rv not necessarily fully evaluated
+            c.add_clause("p = 0 -> p&q = 0", { d, C.ule(r(), p()) }, true);
+        else if (qv.is_zero() && !rv.is_zero())  // rv not necessarily fully evaluated
+            c.add_clause("q = 0 -> p&q = 0", { d, C.ule(r(), q()) }, true);
         // p = a && q = b ==> r = a & b
-        if (pv.is_val() && qv.is_val() && !rv.is_val()) {
+        else if (pv.is_val() && qv.is_val() && !rv.is_val()) {
             // Just assign by this very weak justification. It will be strengthened in saturation in case of a conflict
             LOG(p() << " = " << pv << " and " << q() << " = " << qv << " yields [band] " << r() << " = " << bitwise_and(pv.val(), qv.val()));
-            return s.mk_clause(~andc, ~s.eq(p(), pv), ~s.eq(q(), qv), s.eq(r(), bitwise_and(pv.val(), qv.val())), true);
+            c.add_clause("p = a & q = b => r = a&b", { d, ~C.eq(p(), pv), ~C.eq(q(), qv), C.eq(r(), bitwise_and(pv.val(), qv.val())) }, true);
         }
-
-        return {};
     }
 
-
+#if 0
 
     /**
      * Produce lemmas for constraint: r == inv p
@@ -490,15 +472,15 @@ namespace polysat {
 
         // p = 0  ==>  r = 0
         if (pv.is_zero())
-            return s.mk_clause(~invc, ~s.eq(p()), s.eq(r()), true);
+            c.add_clause(~invc, ~C.eq(p()), C.eq(r()), true);
         // r = 0  ==>  p = 0
         if (rv.is_zero())
-            return s.mk_clause(~invc, ~s.eq(r()), s.eq(p()), true);
+            c.add_clause(~invc, ~C.eq(r()), C.eq(p()), true);
 
         // forward propagation: p assigned  ==>  r = pseudo_inverse(eval(p))
         // TODO: (later) this should be propagated instead of adding a clause
         /*if (pv.is_val() && !rv.is_val())
-            return s.mk_clause(~invc, ~s.eq(p(), pv), s.eq(r(), pv.val().pseudo_inverse(m.power_of_2())), true);*/
+            c.add_clause(~invc, ~C.eq(p(), pv), C.eq(r(), pv.val().pseudo_inverse(m.power_of_2())), true);*/
 
         if (!pv.is_val() || !rv.is_val())
             return {};
@@ -511,7 +493,7 @@ namespace polysat {
 
         // p != 0  ==>  odd(r)
         if (parity_rv != 0)
-            return s.mk_clause("r = inv p  &  p != 0  ==>  odd(r)", {~invc, s.eq(p()), s.odd(r())}, true);
+            c.add_clause("r = inv p  &  p != 0  ==>  odd(r)", {~invc, C.eq(p()), s.odd(r())}, true);
 
         pdd prod = p() * r();
         rational prodv = (pv * rv).val();
@@ -527,13 +509,13 @@ namespace polysat {
                 LOG("Its in [" << lower << "; " << upper << ")");
                 // parity(p) >= k  ==>  p * r >= 2^k
                 if (prodv < rational::power_of_two(middle))
-                    return s.mk_clause("r = inv p  &  parity(p) >= k  ==>  p*r >= 2^k",
+                    c.add_clause("r = inv p  &  parity(p) >= k  ==>  p*r >= 2^k",
                         {~invc, ~s.parity_at_least(p(), middle), s.uge(prod, rational::power_of_two(middle))}, false);
                 // parity(p) >= k  ==>  r <= 2^(N - k) - 1     (because r is the smallest pseudo-inverse)
                 rational const max_rv = rational::power_of_two(m.power_of_2() - middle) - 1;
                 if (rv.val() > max_rv)
-                    return s.mk_clause("r = inv p  &  parity(p) >= k  ==>  r <= 2^(N - k) - 1",
-                        {~invc, ~s.parity_at_least(p(), middle), s.ule(r(), max_rv)}, false);
+                    c.add_clause("r = inv p  &  parity(p) >= k  ==>  r <= 2^(N - k) - 1",
+                        {~invc, ~s.parity_at_least(p(), middle), C.ule(r(), max_rv)}, false);
             }
             else { // parity less than middle
                 SASSERT(parity_pv < middle);
@@ -541,8 +523,8 @@ namespace polysat {
                 LOG("Its in [" << lower << "; " << upper << ")");
                 // parity(p) < k   ==>  p * r <= 2^k - 1
                 if (prodv > rational::power_of_two(middle))
-                    return s.mk_clause("r = inv p  &  parity(p) < k  ==>  p*r <= 2^k - 1",
-                        {~invc, s.parity_at_least(p(), middle), s.ule(prod, rational::power_of_two(middle) - 1)}, false);
+                    c.add_clause("r = inv p  &  parity(p) < k  ==>  p*r <= 2^k - 1",
+                        {~invc, s.parity_at_least(p(), middle), C.ule(prod, rational::power_of_two(middle) - 1)}, false);
             }
         }
          // Why did it evaluate to false in this case?
@@ -550,114 +532,5 @@ namespace polysat {
         return {};
     }
 
-
-    
-    void op_constraint::activate_udiv(solver& s) {
-        // signed_constraint const udivc(this, true); Do we really need this premiss? We anyway assert these constraints as unit clauses
-
-        pdd const& quot = r();
-        pdd const& rem = m_linked->r();
-        
-        // Axioms for quotient/remainder:
-        //      a = b*q + r
-        //      multiplication does not overflow in b*q
-        //      addition does not overflow in (b*q) + r; for now expressed as: r <= bq+r
-        //      b ≠ 0  ==>  r < b
-        //      b = 0  ==>  q = -1
-        // TODO: when a,b become evaluable, can we actually propagate q,r? doesn't seem like it.
-        //       Maybe we need something like an op_constraint for better propagation.
-        s.add_clause(s.eq(q() * quot + rem - p()), false);
-        s.add_clause(~s.umul_ovfl(q(), quot), false);
-        // r <= b*q+r
-        //  { apply equivalence:  p <= q  <=>  q-p <= -p-1 }
-        // b*q <= -r-1
-        s.add_clause(s.ule(q() * quot, -rem - 1), false);
-
-        auto c_eq = s.eq(q());
-        s.add_clause(c_eq, s.ult(rem, q()), false);
-        s.add_clause(~c_eq, s.eq(quot + 1), false);
-    }
-    
-    /**
-     * Produce lemmas for constraint: r == p / q
-     * q = 0   ==>  r = max_value
-     * p = 0   ==>  r = 0 || r = max_value
-     * q = 1   ==>  r = p
-     */
-    clause_ref op_constraint::lemma_udiv(solver& s, assignment const& a) {
-        auto pv = a.apply_to(p());
-        auto qv = a.apply_to(q());
-        auto rv = a.apply_to(r());
-
-        if (eval_udiv(pv, qv, rv) == l_true)
-            return {};
-
-        signed_constraint const udivc(this, true);
-
-        if (qv.is_zero() && !rv.is_val())
-            return s.mk_clause(~udivc, ~s.eq(q()), s.eq(r(), r().manager().max_value()), true);
-        if (pv.is_zero() && !rv.is_val())
-            return s.mk_clause(~udivc, ~s.eq(p()), s.eq(r()), s.eq(r(), r().manager().max_value()), true);
-        if (qv.is_one()) 
-            return s.mk_clause(~udivc, ~s.eq(q(), 1), s.eq(r(), p()), true);
-        
-        if (pv.is_val() && qv.is_val() && !rv.is_val()) {
-            SASSERT(!qv.is_zero());
-            // TODO: We could actually propagate an interval. Instead of p = 9 & q = 4 => r = 2 we could do p >= 8 && p < 12 && q = 4 => r = 2
-            return s.mk_clause(~udivc, ~s.eq(p(), pv.val()), ~s.eq(q(), qv.val()), s.eq(r(), div(pv.val(), qv.val())), true);
-        }
-        
-        return {};
-    }
-
-    
-    /**
-     * Produce lemmas for constraint: r == p % q
-     * p = 0   ==>  r = 0
-     * q = 1   ==>  r = 0
-     * q = 0   ==>  r = p
-     */
-    clause_ref op_constraint::lemma_urem(solver& s, assignment const& a) {
-        auto pv = a.apply_to(p());
-        auto qv = a.apply_to(q());
-        auto rv = a.apply_to(r());
-
-        if (eval_urem(pv, qv, rv) == l_true)
-            return {};
-
-        signed_constraint const urem(this, true);
-
-        if (pv.is_zero() && !rv.is_val())
-            return s.mk_clause(~urem, ~s.eq(p()), s.eq(r()), true);
-        if (qv.is_one() && !rv.is_val())
-            return s.mk_clause(~urem, ~s.eq(q(), 1), s.eq(r()), true);
-        if (qv.is_zero()) 
-            return s.mk_clause(~urem, ~s.eq(q()), s.eq(r(), p()), true);
-        
-        if (pv.is_val() && qv.is_val() && !rv.is_val()) {
-            SASSERT(!qv.is_zero());
-            return s.mk_clause(~urem, ~s.eq(p(), pv.val()), ~s.eq(q(), qv.val()), s.eq(r(), mod(pv.val(), qv.val())), true);
-        }
-        
-        return {};
-    }
-
-    /** Evaluate constraint: r == p % q */
-    lbool op_constraint::eval_urem(pdd const& p, pdd const& q, pdd const& r) {
-        
-        if (q.is_one() && r.is_val()) {
-            return r.val().is_zero() ? l_true : l_false;
-        }
-        if (q.is_zero()) {
-            if (r == p)
-                return l_true;
-        }
-        
-        if (!p.is_val() || !q.is_val() || !r.is_val())
-            return l_undef;
-
-        return r.val() == mod(p.val(), q.val()) ? l_true : l_false; // mod == rem as we know hat q > 0
-    }
-
 #endif
 }
diff --git a/src/sat/smt/polysat/op_constraint.h b/src/sat/smt/polysat/op_constraint.h
index a33f1b705..d7b6be392 100644
--- a/src/sat/smt/polysat/op_constraint.h
+++ b/src/sat/smt/polysat/op_constraint.h
@@ -50,32 +50,22 @@ namespace polysat {
         
         op_constraint(code c, pdd const& r, pdd const& p, pdd const& q);
         lbool eval(pdd const& r, pdd const& p, pdd const& q) const;
-//        clause_ref produce_lemma(core& s, assignment const& a);
 
-  //      clause_ref lemma_lshr(core& s, assignment const& a);
         static lbool eval_lshr(pdd const& p, pdd const& q, pdd const& r);
-
-    //    clause_ref lemma_shl(core& s, assignment const& a);
         static lbool eval_shl(pdd const& p, pdd const& q, pdd const& r);
-
-      //  clause_ref lemma_and(core& s, assignment const& a);
         static lbool eval_and(pdd const& p, pdd const& q, pdd const& r);
-
-       // clause_ref lemma_inv(core& s, assignment const& a);
         static lbool eval_inv(pdd const& p, pdd const& r);
+
+        void propagate_lshr(core& s, dependency const& dep);
+        void propagate_shl(core& s, dependency const& dep);
+        void propagate_and(core& s, dependency const& dep);
+        void propagate_inv(core& s, dependency const& dep);
+
         
-       // clause_ref lemma_udiv(core& s, assignment const& a);
-        static lbool eval_udiv(pdd const& p, pdd const& q, pdd const& r);
-        
-       // clause_ref lemma_urem(core& s, assignment const& a);
-        static lbool eval_urem(pdd const& p, pdd const& q, pdd const& r);
 
         std::ostream& display(std::ostream& out, char const* eq) const;
 
-        void activate(core& s);
-
-        void activate_and(core& s);
-        void activate_udiv(core& s);
+        void activate_and(core& s, dependency const& d);
 
     public:
         ~op_constraint() override {}
@@ -89,6 +79,8 @@ namespace polysat {
         lbool eval(assignment const& a) const override;
         bool is_always_true() const { return false; }
         bool is_always_false() const { return false; }
+        void activate(core& c, bool sign, dependency const& dep) override;
+        void propagate(core& c, lbool value, dependency const& dep) override;
     };
 
 }
diff --git a/src/sat/smt/polysat/types.h b/src/sat/smt/polysat/types.h
index 207ea091e..6f855f98b 100644
--- a/src/sat/smt/polysat/types.h
+++ b/src/sat/smt/polysat/types.h
@@ -29,12 +29,17 @@ namespace polysat {
     class signed_constraint;
 
     class dependency {
-        std::variant<sat::literal, std::pair<theory_var, theory_var>> m_data;
+        struct axiom_t {};
+        std::variant<axiom_t, sat::literal, std::pair<theory_var, theory_var>> m_data;
         unsigned m_level;
+        dependency(): m_data(axiom_t()), m_level(0) {}
     public:
         dependency(sat::literal lit, unsigned level) : m_data(lit), m_level(level) {}
         dependency(theory_var v1, theory_var v2, unsigned level) : m_data(std::make_pair(v1, v2)), m_level(level) {} 
+        static dependency axiom() { return dependency(); } 
         bool is_null() const { return is_literal() && *std::get_if<sat::literal>(&m_data) == sat::null_literal; }
+        bool is_axiom() const { return std::holds_alternative<axiom_t>(m_data); }
+        bool is_eq() const { return std::holds_alternative<std::pair<theory_var, theory_var>>(m_data); }
         bool is_literal() const { return std::holds_alternative<sat::literal>(m_data); }
         sat::literal literal() const { SASSERT(is_literal()); return *std::get_if<sat::literal>(&m_data); }
         std::pair<theory_var, theory_var> eq() const { SASSERT(!is_literal()); return *std::get_if<std::pair<theory_var, theory_var>>(&m_data); }
@@ -46,6 +51,8 @@ namespace polysat {
     inline std::ostream& operator<<(std::ostream& out, dependency d) {
         if (d.is_null())
             return out << "null";
+        else if (d.is_axiom())
+            return out << "axiom@" << d.level();
         else if (d.is_literal())
             return out << d.literal() << "@" << d.level();
         else
@@ -87,7 +94,7 @@ namespace polysat {
 
     using dependency_vector = vector<dependency>;
 
-    using core_vector = vector<std::variant<signed_constraint, dependency>>;
+    using core_vector = std::initializer_list<std::variant<signed_constraint, dependency>>;
 
 
 
@@ -101,6 +108,7 @@ namespace polysat {
         virtual void add_eq_literal(pvar v, rational const& val) = 0;
         virtual void set_conflict(dependency_vector const& core) = 0;
         virtual void set_lemma(core_vector const& aux_core, unsigned level, dependency_vector const& core) = 0;
+        virtual void add_polysat_clause(char const* name, core_vector cs, bool redundant) = 0;
         virtual dependency propagate(signed_constraint sc, dependency_vector const& deps) = 0;
         virtual void propagate(dependency const& d, bool sign, dependency_vector const& deps) = 0;
         virtual trail_stack& trail() = 0;
diff --git a/src/sat/smt/polysat/ule_constraint.cpp b/src/sat/smt/polysat/ule_constraint.cpp
index 185dad0ee..bdfcb7c5f 100644
--- a/src/sat/smt/polysat/ule_constraint.cpp
+++ b/src/sat/smt/polysat/ule_constraint.cpp
@@ -70,6 +70,8 @@ Useful lemmas:
 
 --*/
 
+#include "util/log.h"
+#include "sat/smt/polysat/core.h"
 #include "sat/smt/polysat/constraints.h"
 #include "sat/smt/polysat/ule_constraint.h"
 
@@ -314,8 +316,6 @@ namespace polysat {
         return display(out, l_true, m_lhs, m_rhs);
     }
 
-
-
     // Evaluate lhs <= rhs
     lbool ule_constraint::eval(pdd const& lhs, pdd const& rhs) {
         // NOTE: don't assume simplifications here because we also call this on partially substituted constraints
@@ -343,4 +343,15 @@ namespace polysat {
         return eval(a.apply_to(lhs()), a.apply_to(rhs()));
     }
 
+    void ule_constraint::activate(core& c, bool sign, dependency const& d) {
+        auto p = c.subst(lhs());
+        auto q = c.subst(rhs());
+        auto& C = c.cs();
+        if (sign && !lhs().is_val() && !rhs().is_val()) {
+            c.add_clause("lhs > rhs  ==>  -1 > rhs", { d, C.ult(rhs(), -1) }, false);
+            c.add_clause("lhs > rhs  ==>  lhs > 0", { d, C.ult(0, lhs()) }, false);
+        }
+    }
+
+
 }
diff --git a/src/sat/smt/polysat/ule_constraint.h b/src/sat/smt/polysat/ule_constraint.h
index aa53e6a4f..81a0b64c5 100644
--- a/src/sat/smt/polysat/ule_constraint.h
+++ b/src/sat/smt/polysat/ule_constraint.h
@@ -35,6 +35,8 @@ namespace polysat {
         std::ostream& display(std::ostream& out) const override;
         lbool eval() const override;
         lbool eval(assignment const& a) const override;
+        void activate(core& c, bool sign, dependency const& dep);
+        void propagate(core& c, lbool value, dependency const& dep) {}
         bool is_eq() const { return m_rhs.is_zero(); }
         unsigned power_of_2() const { return m_lhs.power_of_2(); }
 
diff --git a/src/sat/smt/polysat/umul_ovfl_constraint.cpp b/src/sat/smt/polysat/umul_ovfl_constraint.cpp
index e7dc5801c..5d185e7ee 100644
--- a/src/sat/smt/polysat/umul_ovfl_constraint.cpp
+++ b/src/sat/smt/polysat/umul_ovfl_constraint.cpp
@@ -10,6 +10,8 @@ Author:
     Jakob Rath, Nikolaj Bjorner (nbjorner) 2021-12-09
 
 --*/
+#include "util/log.h"
+#include "sat/smt/polysat/core.h"
 #include "sat/smt/polysat/constraints.h"
 #include "sat/smt/polysat/assignment.h"
 #include "sat/smt/polysat/umul_ovfl_constraint.h"
@@ -70,4 +72,84 @@ namespace polysat {
         return eval(a.apply_to(p()), a.apply_to(q()));
     }
 
+    void umul_ovfl_constraint::activate(core& c, bool sign, dependency const& dep) {
+
+    }
+
+    void umul_ovfl_constraint::propagate(core& c, lbool value, dependency const& dep) {
+        auto& C = c.cs();
+        auto p1 = c.subst(p());
+        auto q1 = c.subst(q());
+        if (narrow_bound(c, value == l_true, p(), q(), p1, q1, dep))
+            return;
+        if (narrow_bound(c, value == l_true, q(), p(), q1, p1, dep))
+            return;
+    }
+
+    /**
+     * if p constant, q, propagate inequality
+     */
+    bool umul_ovfl_constraint::narrow_bound(core& c, bool is_positive, pdd const& p0, pdd const& q0, pdd const& p, pdd const& q, dependency const& d) {
+        LOG("p: " << p0 << " := " << p);
+        LOG("q: " << q0 << " := " << q);
+
+        if (!p.is_val())
+            return false;
+        VERIFY(!p.is_zero() && !p.is_one());  // evaluation should catch this case
+
+        rational const& M = p.manager().two_to_N();
+        auto& C = c.cs();
+
+        // q_bound
+        // = min q . Ovfl(p_val, q)
+        // = min q . p_val * q >= M
+        // = min q . q >= M / p_val
+        // = ceil(M / p_val)
+        rational const q_bound = ceil(M / p.val());
+        SASSERT(2 <= q_bound && q_bound <= M / 2);
+        SASSERT(p.val() * q_bound >= M);
+        SASSERT(p.val() * (q_bound - 1) < M);
+        // LOG("q_bound: " << q.manager().mk_val(q_bound));
+
+        // We need the following properties for the bounds:
+        //
+        //      p_bound * (q_bound - 1) < M
+        //      p_bound * q_bound >= M
+        //
+        // With these properties we get:
+        //
+        //      p <= p_bound & q  < q_bound ==> ~Ovfl(p, q)
+        //      p >= p_bound & q >= q_bound ==> Ovfl(p, q)
+        //
+        // Written as lemmas:
+        //
+        //       Ovfl(p, q) & p <= p_bound ==> q >= q_bound
+        //      ~Ovfl(p, q) & p >= p_bound ==> q <  q_bound
+        //
+        if (is_positive) {
+            // Find largest bound for p such that q_bound is still correct.
+            // p_bound = max p . (q_bound - 1)*p < M
+            //         = max p . p < M / (q_bound - 1)
+            //         = ceil(M / (q_bound - 1)) - 1
+            rational const p_bound = ceil(M / (q_bound - 1)) - 1;
+            SASSERT(p.val() <= p_bound);
+            SASSERT(p_bound * q_bound >= M);
+            SASSERT(p_bound * (q_bound - 1) < M);
+            // LOG("p_bound: " << p.manager().mk_val(p_bound));
+            c.add_clause("~Ovfl(p, q) & p <= p_bound ==> q < q_bound", { d, ~C.ule(p0, p_bound), C.ule(q_bound, q0) }, false);
+        }
+        else {
+            // Find lowest bound for p such that q_bound is still correct.
+            // p_bound = min p . Ovfl(p, q_bound) = ceil(M / q_bound)
+            rational const p_bound = ceil(M / q_bound);
+            SASSERT(p_bound <= p.val());
+            SASSERT(p_bound * q_bound >= M);
+            SASSERT(p_bound * (q_bound - 1) < M);
+            // LOG("p_bound: " << p.manager().mk_val(p_bound));
+            c.add_clause("~Ovfl(p, q) & p >= p_bound ==> q < q_bound", { d, ~C.ule(p_bound, p0), C.ult(q0, q_bound) }, false);
+        }
+        return true;
+    }
+
+
 }
diff --git a/src/sat/smt/polysat/umul_ovfl_constraint.h b/src/sat/smt/polysat/umul_ovfl_constraint.h
index c9d03fb01..374d346f7 100644
--- a/src/sat/smt/polysat/umul_ovfl_constraint.h
+++ b/src/sat/smt/polysat/umul_ovfl_constraint.h
@@ -25,6 +25,8 @@ namespace polysat {
         static bool is_always_false(bool is_positive, pdd const& p, pdd const& q) { return is_always_true(!is_positive, p, q); }
         static lbool eval(pdd const& p, pdd const& q);
 
+        bool narrow_bound(core& c, bool is_positive, pdd const& p0, pdd const& q0, pdd const& p, pdd const& q, dependency const& d);
+
     public:
         umul_ovfl_constraint(pdd const& p, pdd const& q);
         ~umul_ovfl_constraint() override {}
@@ -34,6 +36,8 @@ namespace polysat {
         std::ostream& display(std::ostream& out) const override;
         lbool eval() const override;
         lbool eval(assignment const& a) const override;
+        void activate(core& c, bool sign, dependency const& dep) override;
+        void propagate(core& c, lbool value, dependency const& dep) override;
     };
 
 }
diff --git a/src/sat/smt/polysat_internalize.cpp b/src/sat/smt/polysat_internalize.cpp
index 5e5647bd3..ef469fe6f 100644
--- a/src/sat/smt/polysat_internalize.cpp
+++ b/src/sat/smt/polysat_internalize.cpp
@@ -235,9 +235,7 @@ namespace polysat {
         if (n->get_num_args() == 2) {
             expr* x, * y;
             VERIFY(bv.is_bv_and(n, x, y));
-            auto sc = m_core.band(expr2pdd(x), expr2pdd(y), expr2pdd(n));
-            // auto index = m_core.register_constraint(sc, dependency::axiom());
-            // 
+            m_core.band(expr2pdd(x), expr2pdd(y), expr2pdd(n));
         }
         else {
             expr_ref z(n->get_arg(0), m);
@@ -252,19 +250,19 @@ namespace polysat {
     void solver::internalize_lshr(app* n) {
         expr* x, * y;
         VERIFY(bv.is_bv_lshr(n, x, y));
-        auto sc = m_core.lshr(expr2pdd(x), expr2pdd(y), expr2pdd(n));
+        m_core.lshr(expr2pdd(x), expr2pdd(y), expr2pdd(n));
     }
 
     void solver::internalize_ashr(app* n) {
         expr* x, * y;
         VERIFY(bv.is_bv_ashr(n, x, y));
-        auto sc = m_core.ashr(expr2pdd(x), expr2pdd(y), expr2pdd(n));
+        m_core.ashr(expr2pdd(x), expr2pdd(y), expr2pdd(n));
     }
 
     void solver::internalize_shl(app* n) {
         expr* x, * y;
         VERIFY(bv.is_bv_shl(n, x, y));
-        auto sc = m_core.shl(expr2pdd(x), expr2pdd(y), expr2pdd(n));
+        m_core.shl(expr2pdd(x), expr2pdd(y), expr2pdd(n));
     }  
 
     void solver::internalize_urem_i(app* rem) {
diff --git a/src/sat/smt/polysat_solver.cpp b/src/sat/smt/polysat_solver.cpp
index 9f185b22d..690548aaa 100644
--- a/src/sat/smt/polysat_solver.cpp
+++ b/src/sat/smt/polysat_solver.cpp
@@ -64,7 +64,7 @@ namespace polysat {
         case sat::check_result::CR_GIVEUP: {
             if (!m.inc())
                 return sat::check_result::CR_GIVEUP;
-            switch (m_intblast.check()) {
+            switch (m_intblast.check_solver_state()) {
             case l_true:
                 trail().push(value_trail(m_use_intblast_model));
                 m_use_intblast_model = true;
@@ -254,10 +254,25 @@ namespace polysat {
         return ctx.get_trail_stack();
     }
 
-    void solver::add_polysat_clause(char const* name, std::initializer_list<signed_constraint> cs, bool is_redundant) {
+    void solver::add_polysat_clause(char const* name, core_vector cs, bool is_redundant) {
         sat::literal_vector lits;
-        for (auto sc : cs)
-            lits.push_back(ctx.mk_literal(constraint2expr(sc)));
+        for (auto e : cs) {
+            if (std::holds_alternative<dependency>(e)) {
+                auto d = *std::get_if<dependency>(&e);
+                SASSERT(!d.is_null());
+                if (d.is_literal())
+                    lits.push_back(~d.literal());
+                else if (d.is_eq()) {
+                    auto [v1, v2] = d.eq();
+                    lits.push_back(~eq_internalize(var2enode(v1), var2enode(v2)));
+                }
+                else {
+                    SASSERT(d.is_axiom());
+                }
+            }
+            else if (std::holds_alternative<signed_constraint>(e))
+                lits.push_back(ctx.mk_literal(constraint2expr(*std::get_if<signed_constraint>(&e))));
+        }
         s().add_clause(lits.size(), lits.data(), sat::status::th(is_redundant, get_id(), nullptr));
     }
 
diff --git a/src/sat/smt/polysat_solver.h b/src/sat/smt/polysat_solver.h
index f54bafb1c..a04c76618 100644
--- a/src/sat/smt/polysat_solver.h
+++ b/src/sat/smt/polysat_solver.h
@@ -151,7 +151,7 @@ namespace polysat {
         bool inconsistent() const override;
         void get_bitvector_prefixes(pvar v, pvar_vector& out) override;
         void get_fixed_bits(pvar v, svector<justified_fixed_bits>& fixed_bits) override;
-        void add_polysat_clause(char const* name, std::initializer_list<signed_constraint> cs, bool is_redundant);
+        void add_polysat_clause(char const* name, core_vector cs, bool redundant) override;
 
         std::pair<sat::literal_vector, euf::enode_pair_vector> explain_deps(dependency_vector const& deps);