diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp
index 7e4c00873..ca0f22abf 100644
--- a/src/math/dd/dd_pdd.cpp
+++ b/src/math/dd/dd_pdd.cpp
@@ -1711,6 +1711,10 @@ namespace dd {
         return p.val();
     }
 
+    pdd pdd::shl(unsigned n) const {
+        return (*this) * rational::power_of_two(n);
+    }
+
     std::ostream& operator<<(std::ostream& out, pdd const& b) { return b.display(out); }
 
     void pdd_iterator::next() {
diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h
index b81d5c915..421790e7f 100644
--- a/src/math/dd/dd_pdd.h
+++ b/src/math/dd/dd_pdd.h
@@ -428,6 +428,7 @@ namespace dd {
         pdd operator*(rational const& other) const { return m.mul(other, *this); }
         pdd operator+(rational const& other) const { return m.add(other, *this); }
         pdd operator~() const { return m.mk_not(*this); }
+        pdd shl(unsigned n) const; 
         pdd rev_sub(rational const& r) const { return m.sub(m.mk_val(r), *this); }
         pdd div(rational const& other) const { return m.div(*this, other); }
         bool try_div(rational const& other, pdd& out_result) const { return m.try_div(*this, other, out_result); }
diff --git a/src/math/polysat/CMakeLists.txt b/src/math/polysat/CMakeLists.txt
index 6c62db8ef..87633140a 100644
--- a/src/math/polysat/CMakeLists.txt
+++ b/src/math/polysat/CMakeLists.txt
@@ -5,19 +5,20 @@ z3_add_component(polysat
     clause_builder.cpp
     conflict.cpp
     constraint.cpp
+    eq_explain.cpp
     explain.cpp
     forbidden_intervals.cpp
     justification.cpp
     linear_solver.cpp
     log.cpp
     mul_ovfl_constraint.cpp
+    op_constraint.cpp
     restart.cpp
     saturation.cpp
     search_state.cpp
     simplify.cpp
-    op_constraint.cpp
+    smul_ovfl_constraint.cpp
     solver.cpp
-    solve_explain.cpp
     ule_constraint.cpp
     viable.cpp
     variable_elimination.cpp
diff --git a/src/math/polysat/boolean.cpp b/src/math/polysat/boolean.cpp
index 3373cae07..51724a319 100644
--- a/src/math/polysat/boolean.cpp
+++ b/src/math/polysat/boolean.cpp
@@ -54,7 +54,8 @@ namespace polysat {
         m_deps[var] = null_dependency;
         m_watch[lit.index()].reset();
         m_watch[(~lit).index()].reset();
-        m_free_vars.del_var_eh(var);
+        if (m_tracked.get(var, false))
+            m_free_vars.del_var_eh(var);
         // TODO: this is disabled for now, since re-using variables for different constraints may be confusing during debugging. Should be enabled later.
         // m_unused.push_back(var);
     }
@@ -98,7 +99,8 @@ namespace polysat {
         m_kind[lit.var()] = k;
         m_clause[lit.var()] = reason;
         m_deps[lit.var()] = dep;
-        m_free_vars.del_var_eh(lit.var());
+        if (m_tracked.get(lit.var(), false))
+            m_free_vars.del_var_eh(lit.var());
     }
 
     void bool_var_manager::unassign(sat::literal lit) {
@@ -109,7 +111,8 @@ namespace polysat {
         m_kind[lit.var()] = kind_t::unassigned;
         m_clause[lit.var()] = nullptr;
         m_deps[lit.var()] = null_dependency;
-        m_free_vars.unassign_var_eh(lit.var());
+        if (m_tracked.get(lit.var(), false))
+            m_free_vars.unassign_var_eh(lit.var());
     }
 
     std::ostream& bool_var_manager::display(std::ostream& out) const {
diff --git a/src/math/polysat/boolean.h b/src/math/polysat/boolean.h
index 7b85dcd70..0f210e121 100644
--- a/src/math/polysat/boolean.h
+++ b/src/math/polysat/boolean.h
@@ -32,6 +32,7 @@ namespace polysat {
         svector<lbool>          m_value;    // current value (indexed by literal)
         unsigned_vector         m_level;    // level of assignment (indexed by variable)
         dependency_vector       m_deps;     // dependencies of external asserts
+        bool_vector             m_tracked;  // is the variable tracked
         unsigned_vector         m_activity; //
         svector<kind_t>         m_kind;     // decision or propagation?
         // Clause associated with the assignment (indexed by variable):
@@ -80,7 +81,7 @@ namespace polysat {
         unsigned_vector& activity() { return m_activity; }
         bool can_decide() const { return !m_free_vars.empty(); }
         sat::bool_var next_var() { return m_free_vars.next_var(); }
-        void track_var(sat::literal lit) { m_free_vars.mk_var_eh(lit.var()); }
+        void track_var(sat::literal lit) { m_tracked.setx(lit.var(), true, false);  m_free_vars.mk_var_eh(lit.var()); }
 
         // TODO connect activity updates with solver
         void inc_activity(sat::literal lit) { m_activity[lit.var()]++; }
diff --git a/src/math/polysat/conflict.cpp b/src/math/polysat/conflict.cpp
index db75f46b3..a261b4024 100644
--- a/src/math/polysat/conflict.cpp
+++ b/src/math/polysat/conflict.cpp
@@ -26,7 +26,7 @@ Notes:
 #include "math/polysat/log.h"
 #include "math/polysat/log_helper.h"
 #include "math/polysat/explain.h"
-#include "math/polysat/solve_explain.h"
+#include "math/polysat/eq_explain.h"
 #include "math/polysat/forbidden_intervals.h"
 #include "math/polysat/saturation.h"
 #include "math/polysat/variable_elimination.h"
@@ -36,7 +36,7 @@ namespace polysat {
 
     conflict::conflict(solver& s):s(s) {
         ex_engines.push_back(alloc(ex_polynomial_superposition, s));
-        ex_engines.push_back(alloc(solve_explain, s));
+        ex_engines.push_back(alloc(eq_explain, s));
         ve_engines.push_back(alloc(ve_reduction));
         inf_engines.push_back(alloc(inf_saturate, s));
     }
diff --git a/src/math/polysat/constraint.cpp b/src/math/polysat/constraint.cpp
index 67a357fd2..a3ae7a1d5 100644
--- a/src/math/polysat/constraint.cpp
+++ b/src/math/polysat/constraint.cpp
@@ -19,6 +19,7 @@ Author:
 #include "math/polysat/log_helper.h"
 #include "math/polysat/ule_constraint.h"
 #include "math/polysat/mul_ovfl_constraint.h"
+#include "math/polysat/smul_ovfl_constraint.h"
 #include "math/polysat/op_constraint.h"
 
 namespace polysat {
@@ -237,6 +238,10 @@ namespace polysat {
         return { dedup(alloc(mul_ovfl_constraint, *this, a, b)), true };
     }
 
+    signed_constraint constraint_manager::smul_ovfl(pdd const& a, pdd const& b) {
+        return { dedup(alloc(smul_ovfl_constraint, *this, a, b)), true };
+    }
+
     signed_constraint constraint_manager::lshr(pdd const& p, pdd const& q, pdd const& r) {
         return { dedup(alloc(op_constraint, *this, op_constraint::code::lshr_op, p, q, r)), true };
     }
@@ -291,6 +296,14 @@ namespace polysat {
         return *dynamic_cast<mul_ovfl_constraint const*>(this);
     }
 
+    smul_ovfl_constraint& constraint::to_smul_ovfl() {
+        return *dynamic_cast<smul_ovfl_constraint*>(this);
+    }
+
+    smul_ovfl_constraint const& constraint::to_smul_ovfl() const {
+        return *dynamic_cast<smul_ovfl_constraint const*>(this);
+    }
+
     op_constraint& constraint::to_op() {
         return *dynamic_cast<op_constraint*>(this);
     }
@@ -320,7 +333,6 @@ namespace polysat {
             if (!s.is_assigned(other_v)) {
                 s.add_watch({this, is_positive}, other_v);
                 std::swap(vars()[idx], vars()[i]);
-               // narrow(s, is_positive);
                 return true;
             }
         }
@@ -333,7 +345,7 @@ namespace polysat {
     void constraint::propagate_core(solver& s, bool is_positive, pvar v, pvar other_v) {
         (void)v;
         (void)other_v;
-        narrow(s, is_positive);
+        narrow(s, is_positive, false);
     }
 
 
diff --git a/src/math/polysat/constraint.h b/src/math/polysat/constraint.h
index db131acba..ec508a55c 100644
--- a/src/math/polysat/constraint.h
+++ b/src/math/polysat/constraint.h
@@ -20,11 +20,12 @@ Author:
 
 namespace polysat {
 
-    enum ckind_t { ule_t, mul_ovfl_t, op_t };
+    enum ckind_t { ule_t, mul_ovfl_t, smul_ovfl_t, op_t };
 
     class constraint;
     class ule_constraint;
     class mul_ovfl_constraint;
+    class smul_ovfl_constraint;
     class op_constraint;
     class signed_constraint;
 
@@ -98,6 +99,7 @@ namespace polysat {
         signed_constraint sle(pdd const& a, pdd const& b);
         signed_constraint slt(pdd const& a, pdd const& b);
         signed_constraint mul_ovfl(pdd const& p, pdd const& q);
+        signed_constraint smul_ovfl(pdd const& p, pdd const& q);
         signed_constraint bit(pdd const& p, unsigned i);
         signed_constraint lshr(pdd const& p, pdd const& q, pdd const& r);
         signed_constraint band(pdd const& p, pdd const& q, pdd const& r);
@@ -139,6 +141,7 @@ namespace polysat {
         friend class clause;
         friend class ule_constraint;
         friend class mul_ovfl_constraint;
+        friend class smul_ovfl_constraint;
         friend class op_constraint;
 
         // constraint_manager* m_manager;
@@ -167,6 +170,7 @@ namespace polysat {
         virtual bool is_diseq() const { return false; }
         bool is_ule() const { return m_kind == ckind_t::ule_t; }
         bool is_mul_ovfl() const { return m_kind == ckind_t::mul_ovfl_t; }
+        bool is_smul_ovfl() const { return m_kind == ckind_t::smul_ovfl_t; }
         bool is_op() const { return m_kind == ckind_t::op_t; }
         ckind_t kind() const { return m_kind; }
         virtual std::ostream& display(std::ostream& out, lbool status) const = 0;
@@ -179,13 +183,15 @@ namespace polysat {
         virtual bool is_currently_true(solver& s, bool is_positive) const = 0;
         virtual bool is_currently_false(solver& s, assignment_t const& sub, bool is_positive) const = 0;
         virtual bool is_currently_true(solver& s, assignment_t const& sub, bool is_positive) const = 0;
-        virtual void narrow(solver& s, bool is_positive) = 0;
+        virtual void narrow(solver& s, bool is_positive, bool first) = 0;
         virtual inequality as_inequality(bool is_positive) const = 0;
 
         ule_constraint& to_ule();
         ule_constraint const& to_ule() const;
         mul_ovfl_constraint& to_mul_ovfl();
         mul_ovfl_constraint const& to_mul_ovfl() const;
+        smul_ovfl_constraint& to_smul_ovfl();
+        smul_ovfl_constraint const& to_smul_ovfl() const;
         op_constraint& to_op();
         op_constraint const& to_op() const;
         unsigned_vector& vars() { return m_vars; }
@@ -246,7 +252,7 @@ namespace polysat {
         bool is_currently_false(solver& s, assignment_t const& sub) const { return get()->is_currently_false(s, sub, is_positive()); }
         lbool bvalue(solver& s) const;
         unsigned level(solver& s) const { return get()->level(s); }
-        void narrow(solver& s) { get()->narrow(s, is_positive()); }
+        void narrow(solver& s, bool first) { get()->narrow(s, is_positive(), first); }
         inequality as_inequality() const { return get()->as_inequality(is_positive()); }
 
         sat::bool_var bvar() const { return m_constraint->bvar(); }
diff --git a/src/math/polysat/forbidden_intervals.cpp b/src/math/polysat/forbidden_intervals.cpp
index 0007abf08..cd5db069e 100644
--- a/src/math/polysat/forbidden_intervals.cpp
+++ b/src/math/polysat/forbidden_intervals.cpp
@@ -17,6 +17,7 @@ Author:
 #include "math/polysat/interval.h"
 #include "math/polysat/solver.h"
 #include "math/polysat/log.h"
+#include "math/polysat/mul_ovfl_constraint.h"
 
 namespace polysat {
 
@@ -28,20 +29,112 @@ namespace polysat {
      * \returns True iff a forbidden interval exists and the output parameters were set.
      */
     bool forbidden_intervals::get_interval(signed_constraint const& c, pvar v, fi_record& fi) {
-        if (!c->is_ule())
-            return false;
-        
-        struct backtrack {
-            bool released = false;
-            vector<signed_constraint>& side_cond;
-            unsigned sz;
-            backtrack(vector<signed_constraint>& s):side_cond(s), sz(s.size()) {}
-            ~backtrack() {
-                if (!released)
-                    side_cond.shrink(sz);
-            }
-        };
+        if (c->is_ule())
+            return get_interval_ule(c, v, fi);
+        if (c->is_mul_ovfl())
+            return get_interval_mul_ovfl(c, v, fi);
+        return false;
+    }
 
+    
+    
+    bool forbidden_intervals::get_interval_mul_ovfl(signed_constraint const& c, pvar v, fi_record& fi) {
+
+
+        std::cout << "FORBIDDEN v" << v << "\n";
+
+        backtrack _backtrack(fi.side_cond);     
+
+        fi.coeff = 1;
+        fi.src = c;
+
+        // eval(lhs) = a1*v + eval(e1) = a1*v + b1
+        // eval(rhs) = a2*v + eval(e2) = a2*v + b2
+        // We keep the e1, e2 around in case we need side conditions such as e1=b1, e2=b2.
+        auto [ok1, a1, e1, b1] = linear_decompose(v, c->to_mul_ovfl().p(), fi.side_cond);
+        auto [ok2, a2, e2, b2] = linear_decompose(v, c->to_mul_ovfl().q(), fi.side_cond);
+
+        auto& m = e1.manager();
+        rational bound = m.max_value();
+
+        if (ok2 && !ok1) {
+            std::swap(a1, a2);
+            std::swap(e1, e2);
+            std::swap(b1, b2);
+            std::swap(ok1, ok2);            
+        }
+        if (ok1 && !ok2 && a1.is_one() && b1.is_zero()) {
+            if (c.is_positive()) {
+                _backtrack.released = true;
+                rational lo_val(0);
+                rational hi_val(2);
+                pdd lo = m.mk_val(lo_val);
+                pdd hi = m.mk_val(hi_val);
+                fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val);
+                return true;
+            }
+        }
+        
+        if (!ok1 || !ok2) 
+            return false;
+
+
+        if (a2.is_one() && a1.is_zero()) {
+            std::swap(a1, a2);
+            std::swap(e1, e2);
+            std::swap(b1, b2);
+        }
+
+        if (!a1.is_one() || !a2.is_zero())
+            return false;
+
+        if (!b1.is_zero())
+            return false;
+
+        _backtrack.released = true;
+
+        // Ovfl(v, e2)
+
+
+        if (c.is_positive()) {
+            if (b2.val() <= 1) {
+                fi.interval = eval_interval::full();
+                fi.side_cond.push_back(s.ule(e2, 1));
+            }
+            else {
+                // [0, div(bound, b2.val()) + 1[
+                rational lo_val(0);
+                rational hi_val(div(bound, b2.val()) + 1);
+                pdd lo = m.mk_val(lo_val);
+                pdd hi = m.mk_val(hi_val);
+                fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val);
+                fi.side_cond.push_back(s.ule(e2, b2.val()));
+            }
+
+        }
+        else {
+            if (b2.val() <= 1) {
+                _backtrack.released = false;
+                return false;
+            }
+            else {
+                // [div(bound, b2.val()) + 1, 0[
+                rational lo_val(div(bound, b2.val()) + 1);
+                rational hi_val(0);
+                pdd lo = m.mk_val(lo_val);
+                pdd hi = m.mk_val(hi_val);
+                fi.interval = eval_interval::proper(lo, lo_val, hi, hi_val);
+                fi.side_cond.push_back(s.ule(b2.val(), e2));
+            }
+        }
+
+        LOG("overflow interval " << fi.interval);
+
+        return true;
+    }
+    
+    bool forbidden_intervals::get_interval_ule(signed_constraint const& c, pvar v, fi_record& fi) {
+        
         backtrack _backtrack(fi.side_cond);     
 
         fi.coeff = 1;
diff --git a/src/math/polysat/forbidden_intervals.h b/src/math/polysat/forbidden_intervals.h
index 2d88216d3..17e84c370 100644
--- a/src/math/polysat/forbidden_intervals.h
+++ b/src/math/polysat/forbidden_intervals.h
@@ -70,6 +70,22 @@ namespace polysat {
             rational const & a2, pdd const& b2, pdd const& e2,
             fi_record& fi);
 
+        bool get_interval_ule(signed_constraint const& c, pvar v, fi_record& fi);
+        
+        bool get_interval_mul_ovfl(signed_constraint const& c, pvar v, fi_record& fi);
+
+        struct backtrack {
+            bool released = false;
+            vector<signed_constraint>& side_cond;
+            unsigned sz;
+            backtrack(vector<signed_constraint>& s):side_cond(s), sz(s.size()) {}
+            ~backtrack() {
+                if (!released)
+                    side_cond.shrink(sz);
+            }
+        };
+
+
     public:
         forbidden_intervals(solver& s) :s(s) {}
         bool get_interval(signed_constraint const& c, pvar v, fi_record& fi);
diff --git a/src/math/polysat/mul_ovfl_constraint.cpp b/src/math/polysat/mul_ovfl_constraint.cpp
index b07f1efa6..c52c999de 100644
--- a/src/math/polysat/mul_ovfl_constraint.cpp
+++ b/src/math/polysat/mul_ovfl_constraint.cpp
@@ -90,7 +90,7 @@ namespace polysat {
         return is_always_true(is_positive, s.subst(p()), s.subst(q()));
     }
 
-    void mul_ovfl_constraint::narrow(solver& s, bool is_positive) {    
+    void mul_ovfl_constraint::narrow(solver& s, bool is_positive, bool first) {    
         auto p1 = s.subst(p());
         auto q1 = s.subst(q());
         
@@ -100,10 +100,16 @@ namespace polysat {
         }
         if (is_always_true(is_positive, p1, q1))
             return;
+
+        if (try_viable(s, is_positive, p(), q(), p1, q1))
+            return;
+        
         if (narrow_bound(s, is_positive, p(), q(), p1, q1))
             return;
         if (narrow_bound(s, is_positive, q(), p(), q1, p1))
             return;
+
+        
     }
 
     /**
@@ -145,6 +151,14 @@ namespace polysat {
         return true;
     }
 
+    bool mul_ovfl_constraint::try_viable(
+        solver& s, bool is_positive,
+        pdd const& p0, pdd const& q0, pdd const& p, pdd const& q) {
+        signed_constraint sc(this, is_positive);
+        return s.m_viable.intersect(p0, q0, sc);
+    }
+
+
     unsigned mul_ovfl_constraint::hash() const {
     	return mk_mix(p().hash(), q().hash(), kind());
     }
diff --git a/src/math/polysat/mul_ovfl_constraint.h b/src/math/polysat/mul_ovfl_constraint.h
index 8932c10d6..22ea39273 100644
--- a/src/math/polysat/mul_ovfl_constraint.h
+++ b/src/math/polysat/mul_ovfl_constraint.h
@@ -28,6 +28,7 @@ namespace polysat {
         bool is_always_false(bool is_positive, pdd const& p, pdd const& q) const;
         bool is_always_true(bool is_positive, pdd const& p, pdd const& q) const;
         bool narrow_bound(solver& s, bool is_positive, pdd const& p0, pdd const& q0, pdd const& p, pdd const& q);
+        bool try_viable(solver& s, bool is_positive, pdd const& p0, pdd const& q0, pdd const& p, pdd const& q);
         lbool eval(pdd const& p, pdd const& q) const;
         
     public:
@@ -37,7 +38,7 @@ namespace polysat {
         std::ostream& display(std::ostream& out, lbool status) const override;
         std::ostream& display(std::ostream& out) const override;
         bool is_always_false(bool is_positive) const override;
-        void narrow(solver& s, bool is_positive) override;
+        void narrow(solver& s, bool is_positive, bool first) override;
         bool is_currently_false(solver & s, bool is_positive) const override;
         bool is_currently_true(solver& s, bool is_positive) const override;
         bool is_currently_false(solver& s, assignment_t const& sub, bool is_positive) const override { return false; }
diff --git a/src/math/polysat/op_constraint.cpp b/src/math/polysat/op_constraint.cpp
index 71fd9c684..3c8b866d1 100644
--- a/src/math/polysat/op_constraint.cpp
+++ b/src/math/polysat/op_constraint.cpp
@@ -127,7 +127,7 @@ namespace polysat {
     *
     * We can assume that op_constraint is only asserted positive.
     */
-    void op_constraint::narrow(solver& s, bool is_positive) {
+    void op_constraint::narrow(solver& s, bool is_positive, bool first) {
         SASSERT(is_positive);
 
         if (is_currently_true(s, is_positive))
diff --git a/src/math/polysat/op_constraint.h b/src/math/polysat/op_constraint.h
index f4d7a118a..e3dfec82d 100644
--- a/src/math/polysat/op_constraint.h
+++ b/src/math/polysat/op_constraint.h
@@ -60,7 +60,7 @@ namespace polysat {
         bool is_currently_true(solver& s, bool is_positive) const override;
         bool is_currently_false(solver& s, assignment_t const& sub, bool is_positive) const override { return false; }
         bool is_currently_true(solver& s, assignment_t const& sub, bool is_positive) const override { return false; }
-        void narrow(solver& s, bool is_positive) override;
+        void narrow(solver& s, bool is_positive, bool first) override;
         inequality as_inequality(bool is_positive) const override { throw default_exception("is not an inequality"); }
         unsigned hash() const override;
         bool operator==(constraint const& other) const override;
diff --git a/src/math/polysat/smul_ovfl_constraint.cpp b/src/math/polysat/smul_ovfl_constraint.cpp
new file mode 100644
index 000000000..a0209c393
--- /dev/null
+++ b/src/math/polysat/smul_ovfl_constraint.cpp
@@ -0,0 +1,79 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    polysat multiplication overflow constraint
+
+Author:
+
+    Jakob Rath, Nikolaj Bjorner (nbjorner) 2021-12-09
+
+--*/
+#include "math/polysat/smul_ovfl_constraint.h"
+#include "math/polysat/solver.h"
+
+namespace polysat {
+
+    smul_ovfl_constraint::smul_ovfl_constraint(constraint_manager& m, pdd const& p, pdd const& q):
+        constraint(m, ckind_t::smul_ovfl_t), m_p(p), m_q(q) {
+        simplify();
+        m_vars.append(m_p.free_vars());
+        for (auto v : m_q.free_vars())
+            if (!m_vars.contains(v))
+                m_vars.push_back(v);
+
+    }
+    void smul_ovfl_constraint::simplify() {
+        if (m_p.is_zero() || m_q.is_zero() ||
+            m_p.is_one() || m_q.is_one()) {
+            m_q = 0;
+            m_p = 0;
+            return;
+        }
+        if (m_p.index() > m_q.index())
+            std::swap(m_p, m_q);
+    }
+
+    std::ostream& smul_ovfl_constraint::display(std::ostream& out, lbool status) const {
+        switch (status) {
+        case l_true: return display(out);
+        case l_false: return display(out << "~");
+        case l_undef: return display(out << "?");
+        }       
+        return out;
+    }
+
+    std::ostream& smul_ovfl_constraint::display(std::ostream& out) const {
+        return out << "sovfl*(" << m_p << ", " << m_q << ")";       
+    }
+
+    void smul_ovfl_constraint::narrow(solver& s, bool is_positive, bool first) {
+        if (!first)
+            return;
+        signed_constraint sc(this, is_positive);
+        if (is_positive) {
+            s.add_clause(~sc, s.ule(2, p()), false);
+            s.add_clause(~sc, s.ule(2, q()), false);
+            s.add_clause(~sc, ~s.sgt(p(), 0), s.sgt(q(), 0), false);
+            s.add_clause(~sc, ~s.sgt(q(), 0), s.sgt(p(), 0), false);
+            s.add_clause(~sc, ~s.sgt(p(), 0), s.slt(p()*q(), 0), s.mul_ovfl(p(), q()), false);
+            s.add_clause(~sc, s.sgt(p(), 0),  s.slt(p()*q(), 0), s.mul_ovfl(-p(), -q()), false);
+        }
+        else {
+            // smul_noovfl(p,q) => sign(p) != sign(q) or p'*q' < 2^{K-1}
+            s.add_clause(~sc, ~s.sgt(p(), 0), ~s.sgt(q(), 0), ~s.mul_ovfl(p(), q()), false);
+            s.add_clause(~sc, ~s.sgt(p(), 0), ~s.sgt(q(), 0), ~s.slt(p()*q(), 0), false);
+            s.add_clause(~sc, ~s.slt(p(), 0), ~s.slt(q(), 0), ~s.mul_ovfl(-p(), -q()), false);
+            s.add_clause(~sc, ~s.slt(p(), 0), ~s.slt(q(), 0), ~s.slt((-p())*(-q()), 0), false);
+        }
+    }
+
+    unsigned smul_ovfl_constraint::hash() const {
+    	return mk_mix(p().hash(), q().hash(), kind());
+    }
+
+    bool smul_ovfl_constraint::operator==(constraint const& other) const {
+        return other.is_smul_ovfl() && p() == other.to_smul_ovfl().p() && q() == other.to_smul_ovfl().q();
+    }
+}
diff --git a/src/math/polysat/smul_ovfl_constraint.h b/src/math/polysat/smul_ovfl_constraint.h
new file mode 100644
index 000000000..4c07d7ac1
--- /dev/null
+++ b/src/math/polysat/smul_ovfl_constraint.h
@@ -0,0 +1,48 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    polysat multiplication overflow constraint
+
+Author:
+
+    Jakob Rath, Nikolaj Bjorner (nbjorner) 2021-12-09
+
+--*/
+#pragma once
+#include "math/polysat/constraint.h"
+
+namespace polysat {
+
+    class solver;
+
+    class smul_ovfl_constraint final : public constraint {
+        friend class constraint_manager;
+
+        pdd m_p;
+        pdd m_q;
+
+        void simplify();
+        smul_ovfl_constraint(constraint_manager& m, pdd const& p, pdd const& q);
+        
+    public:
+        ~smul_ovfl_constraint() override {}
+        pdd const& p() const { return m_p; }
+        pdd const& q() const { return m_q; }
+        std::ostream& display(std::ostream& out, lbool status) const override;
+        std::ostream& display(std::ostream& out) const override;
+        bool is_always_false(bool is_positive) const override { return false; }
+        void narrow(solver& s, bool is_positive, bool first) override;
+        bool is_currently_false(solver & s, bool is_positive) const override { return false; }
+        bool is_currently_true(solver& s, bool is_positive) const override { return false; }
+        bool is_currently_false(solver& s, assignment_t const& sub, bool is_positive) const override { return false; }
+        bool is_currently_true(solver& s, assignment_t const& sub, bool is_positive) const override { return false; }
+
+        inequality as_inequality(bool is_positive) const override { throw default_exception("is not an inequality"); }
+        unsigned hash() const override;
+        bool operator==(constraint const& other) const override;
+        bool is_eq() const override { return false; }
+    };
+
+}
diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp
index 25f094e33..6fc935cce 100644
--- a/src/math/polysat/solver.cpp
+++ b/src/math/polysat/solver.cpp
@@ -819,7 +819,7 @@ namespace polysat {
         SASSERT(!c->is_active());
         c->set_active(true);
         add_watch(c);
-        c.narrow(*this);
+        c.narrow(*this, true);
 #if ENABLE_LINEAR_SOLVER
         m_linear_solver.activate_constraint(c);
 #endif
diff --git a/src/math/polysat/solver.h b/src/math/polysat/solver.h
index 88cb90914..e77266940 100644
--- a/src/math/polysat/solver.h
+++ b/src/math/polysat/solver.h
@@ -56,6 +56,7 @@ namespace polysat {
         friend class constraint;
         friend class ule_constraint;
         friend class mul_ovfl_constraint;
+        friend class smul_ovfl_constraint;
         friend class op_constraint;
         friend class signed_constraint;
         friend class clause;
@@ -115,7 +116,6 @@ namespace polysat {
 
         search_state             m_search;
         assignment_t const& assignment() const { return m_search.assignment(); }
-        pdd subst(pdd const& p) const; 
         pdd subst(assignment_t const& sub, pdd const& p) const;
 
         unsigned                 m_qhead = 0; // next item to propagate (index into m_search)
@@ -309,6 +309,12 @@ namespace polysat {
          */
         bool try_eval(pdd const& p, rational& out_value) const;
 
+        /**
+         * Apply current substitution to p.
+         */
+        pdd subst(pdd const& p) const; 
+
+
         /** Create constraints */
         signed_constraint eq(pdd const& p) { return m_constraints.eq(p); }
         signed_constraint diseq(pdd const& p) { return ~m_constraints.eq(p); }
@@ -328,8 +334,16 @@ namespace polysat {
         signed_constraint ult(rational const& p, pdd const& q) { return ult(q.manager().mk_val(p), q); }
         signed_constraint sle(pdd const& p, pdd const& q) { return m_constraints.sle(p, q); }
         signed_constraint slt(pdd const& p, pdd const& q) { return m_constraints.slt(p, q); }
+        signed_constraint slt(pdd const& p, rational const& q) { return slt(p, p.manager().mk_val(q)); }
+        signed_constraint slt(rational const& p, pdd const& q) { return slt(q.manager().mk_val(p), q); }
+        signed_constraint slt(pdd const& p, int n) { return slt(p, rational(n)); }
+        signed_constraint slt(int n, pdd const& p) { return slt(rational(n), p); }
+        signed_constraint sgt(pdd const& p, pdd const& q) { return slt(q, p); }
+        signed_constraint sgt(pdd const& p, int n) { return slt(n, p); }
+        signed_constraint sgt(int n, pdd const& p) { return slt(p, n); }
         signed_constraint mul_ovfl(pdd const& p, pdd const& q) { return m_constraints.mul_ovfl(p, q); }
         signed_constraint mul_ovfl(rational const& p, pdd const& q) { return mul_ovfl(q.manager().mk_val(p), q); }
+        signed_constraint smul_ovfl(pdd const& p, pdd const& q) { return m_constraints.smul_ovfl(p, q); }
         signed_constraint bit(pdd const& p, unsigned i) { return m_constraints.bit(p, i); }
 
         /** Create and activate polynomial constraints. */
diff --git a/src/math/polysat/ule_constraint.cpp b/src/math/polysat/ule_constraint.cpp
index 8045cd42c..974d83eeb 100644
--- a/src/math/polysat/ule_constraint.cpp
+++ b/src/math/polysat/ule_constraint.cpp
@@ -121,7 +121,7 @@ namespace polysat {
         return out << m_lhs << (is_eq() ? " == " : " <= ") << m_rhs;
     }
 
-    void ule_constraint::narrow(solver& s, bool is_positive) {
+    void ule_constraint::narrow(solver& s, bool is_positive, bool first) {
         auto p = s.subst(lhs());
         auto q = s.subst(rhs());
         
@@ -142,34 +142,7 @@ namespace polysat {
             return;
         }
 
-        pvar v = null_var;
-        bool first = true;
-        if (p.is_unilinear())
-            v = p.var();
-        else if (q.is_unilinear())
-            v = q.var(), first = false;
-        else
-            return;
-
-    try_viable:
-        if (s.m_viable.intersect(v, sc)) {
-            rational val;
-            switch (s.m_viable.find_viable(v, val)) {
-            case dd::find_t::singleton:
-                s.propagate(v, val, sc); // TBD why is sc used as justification? It should be all of viable            
-                break;
-            case dd::find_t::empty:
-                s.set_conflict(v);
-                return;
-            default:
-                break;
-            }
-        }
-        if (first && q.is_unilinear() && q.var() != v) {
-            v = q.var();
-            first = false;
-            goto try_viable;
-        }
+        s.m_viable.intersect(p, q, sc);
     }
 
     bool ule_constraint::is_always_false(bool is_positive, pdd const& lhs, pdd const& rhs) const {
diff --git a/src/math/polysat/ule_constraint.h b/src/math/polysat/ule_constraint.h
index 25ac4681f..57fab0a1e 100644
--- a/src/math/polysat/ule_constraint.h
+++ b/src/math/polysat/ule_constraint.h
@@ -40,7 +40,7 @@ namespace polysat {
         bool is_currently_true(solver& s, bool is_positive) const override;
         bool is_currently_false(solver& s, assignment_t const& sub, bool is_positive) const override;
         bool is_currently_true(solver& s, assignment_t const& sub, bool is_positive) const override;
-        void narrow(solver& s, bool is_positive) override;
+        void narrow(solver& s, bool is_positive, bool first) override;
         inequality as_inequality(bool is_positive) const override;
         unsigned hash() const override;
         bool operator==(constraint const& other) const override;
diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp
index b6ca77b3c..2a04d8386 100644
--- a/src/math/polysat/viable.cpp
+++ b/src/math/polysat/viable.cpp
@@ -80,6 +80,40 @@ namespace polysat {
         SASSERT(well_formed(m_units[v]));
         m_trail.pop_back();
     }
+    
+    bool viable::intersect(pdd const & p, pdd const & q, signed_constraint const& sc) {
+        pvar v = null_var;
+        bool first = true;
+        bool prop = false;
+        if (p.is_unilinear())
+            v = p.var();
+        else if (q.is_unilinear())
+            v = q.var(), first = false;
+        else
+            return prop;
+
+    try_viable:
+        if (s.m_viable.intersect(v, sc)) {
+            rational val;
+            switch (s.m_viable.find_viable(v, val)) {
+            case dd::find_t::singleton:
+                s.propagate(v, val, sc); // TBD why is sc used as justification? It should be all of viable
+                prop = true;
+                break;
+            case dd::find_t::empty:
+                s.set_conflict(v);
+                return true;
+            default:
+                break;
+            }
+        }
+        if (first && q.is_unilinear() && q.var() != v) {
+            v = q.var();
+            first = false;
+            goto try_viable;
+        }
+        return prop;
+    }    
 
     bool viable::intersect(pvar v, signed_constraint const& c) {
         auto& fi = s.m_forbidden_intervals;
diff --git a/src/math/polysat/viable.h b/src/math/polysat/viable.h
index 1a5224d13..b4b1880c5 100644
--- a/src/math/polysat/viable.h
+++ b/src/math/polysat/viable.h
@@ -81,6 +81,8 @@ namespace polysat {
          */
         bool intersect(pvar v, signed_constraint const& c);
 
+        bool intersect(pdd const & p, pdd const & q, signed_constraint const& c);
+
         /**
          * Check whether variable v has any viable values left according to m_viable.
          */
diff --git a/src/sat/smt/bv_polysat.cpp b/src/sat/smt/bv_polysat.cpp
index 4fbc59cae..00a929da4 100644
--- a/src/sat/smt/bv_polysat.cpp
+++ b/src/sat/smt/bv_polysat.cpp
@@ -57,6 +57,8 @@ namespace bv {
         case OP_SGT:              polysat_le<true,  false, true>(a); break;
 
         case OP_BUMUL_NO_OVFL:    polysat_umul_noovfl(a); break;
+        case OP_BSMUL_NO_OVFL:    polysat_smul_noovfl(a); break;
+                        
         case OP_BUDIV_I:          polysat_div_rem_i(a, true); break;       
         case OP_BUREM_I:          polysat_div_rem_i(a, false); break;
 
@@ -74,7 +76,6 @@ namespace bv {
         case OP_BSDIV:            
         case OP_BSREM:            
         case OP_BSMOD:                        
-        case OP_BSMUL_NO_OVFL:    
         case OP_BSMUL_NO_UDFL:    
         case OP_BSDIV_I:            
         case OP_BSREM_I:                        
@@ -110,6 +111,15 @@ namespace bv {
         a->m_sc = sc;
     }
 
+    void solver::polysat_smul_noovfl(app* e) {
+        auto p = expr2pdd(e->get_arg(0));
+        auto q = expr2pdd(e->get_arg(1));
+        auto sc = ~m_polysat.smul_ovfl(p, q);
+        sat::literal lit = expr2literal(e);
+        atom* a = mk_atom(lit.var());
+        a->m_sc = sc;
+    }
+
     void solver::polysat_div_rem_i(app* e, bool is_div) {
         auto p = expr2pdd(e->get_arg(0));
         auto q = expr2pdd(e->get_arg(1));
diff --git a/src/sat/smt/bv_solver.h b/src/sat/smt/bv_solver.h
index 85d3502ba..39fbd6d18 100644
--- a/src/sat/smt/bv_solver.h
+++ b/src/sat/smt/bv_solver.h
@@ -281,6 +281,7 @@ namespace bv {
         void polysat_num(app* a);
         void polysat_mkbv(app* a);
         void polysat_umul_noovfl(app* e);
+        void polysat_smul_noovfl(app* e);
         void polysat_div_rem_i(app* e, bool is_div);
         void polysat_div_rem(app* e, bool is_div);
         void polysat_bit2bool(atom* a, expr* e, unsigned idx);