diff --git a/src/sat/smt/polysat/core.cpp b/src/sat/smt/polysat/core.cpp
index de9791bb1..fb49e59c8 100644
--- a/src/sat/smt/polysat/core.cpp
+++ b/src/sat/smt/polysat/core.cpp
@@ -534,7 +534,7 @@ namespace polysat {
 
     }
 
-    void core::add_axiom(signed_constraint sc) {
+    void core::add_opdef(signed_constraint sc) {
         auto idx = register_constraint(sc, dependency::axiom());
         assign_eh(idx, false);
     }
diff --git a/src/sat/smt/polysat/core.h b/src/sat/smt/polysat/core.h
index 645bdaef9..d51bcd839 100644
--- a/src/sat/smt/polysat/core.h
+++ b/src/sat/smt/polysat/core.h
@@ -92,7 +92,7 @@ namespace polysat {
 
         lbool assign_variable();
         
-        void add_axiom(signed_constraint sc);
+        void add_opdef(signed_constraint sc);
 
         unsigned m_activity_inc = 128;
         void inc_activity(pvar v);
@@ -124,10 +124,10 @@ namespace polysat {
         signed_constraint bit(pdd const& p, unsigned i) { return m_constraints.bit(p, i); }
 
 
-        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)); }
+        void lshr(pdd const& a, pdd const& b, pdd const& r) { add_opdef(m_constraints.lshr(a, b, r)); }
+        void ashr(pdd const& a, pdd const& b, pdd const& r) { add_opdef(m_constraints.ashr(a, b, r)); }
+        void shl(pdd const& a, pdd const& b, pdd const& r) { add_opdef(m_constraints.shl(a, b, r)); }
+        void band(pdd const& a, pdd const& b, pdd const& r) { add_opdef(m_constraints.band(a, b, r)); }
 
         pdd bnot(pdd p) { return -p - 1; }
         pdd mul(unsigned n, pdd const* args) { return m_monomials.mk(n, args); }
diff --git a/src/sat/smt/polysat/op_constraint.cpp b/src/sat/smt/polysat/op_constraint.cpp
index 80b58ca8b..bfab8da95 100644
--- a/src/sat/smt/polysat/op_constraint.cpp
+++ b/src/sat/smt/polysat/op_constraint.cpp
@@ -198,13 +198,28 @@ namespace polysat {
     }
 
     void op_constraint::activate(core& c, bool sign, dependency const& dep) {
+        auto& m = p.manager();
+        unsigned const N = m.power_of_2();
+        auto& C = c.cs();
         SASSERT(!sign);
         switch (m_op) {
         case code::and_op:
             activate_and(c, dep);
             break;
         case code::ashr_op:
-            activate_ashr(c, dep);
+            c.add_axiom("q >= N & p < 0 -> p <<a q = -1", { ~C.uge(q, N), ~C.slt(p, 0), C.eq(r, m.max_value()) }, false);
+            c.add_axiom("q >= N & p >= 0 -> p <<a q = 0", { ~C.uge(q, N), ~C.sge(p, 0), C.eq(r) }, false);
+            c.add_axiom("q = 0 -> p <<a q = p", { ~C.eq(q), C.eq(r, p) }, false);
+            break;
+        case code::lshr_op: 
+            c.add_axiom("q >= N  -> p <<l q = 0", { ~C.uge(q, N), C.eq(r) }, false);
+            c.add_axiom("q = 0 -> p <<l q = p", { ~C.eq(q), C.eq(r, p) }, false);
+            break;       
+        case code::shl_op:
+            c.add_axiom("q >= N -> p >> q = 0", { ~C.uge(q, N), C.eq(r) }, false);
+            c.add_axiom("q = 0 -> p >> q = p", { ~C.eq(q), C.eq(r, p) }, false);
+            break;
+        case code::inv_op:
             break;
         default:
             break;
@@ -312,22 +327,13 @@ namespace polysat {
         }
     }
 
-    void op_constraint::activate_ashr(core& c, dependency const& d) {
-        auto& m = p.manager();
-        unsigned const N = m.power_of_2();
-        auto& C = c.cs();
-        c.add_axiom("q >= N & p < 0 -> p << q = -1", {~C.uge(q, N), ~C.slt(p, 0), C.eq(r, m.max_value())}, false);
-        c.add_axiom("q >= N & p >= 0 -> p << q = 0", {~C.uge(q, N), ~C.sge(p, 0), C.eq(r)}, false);
-        c.add_axiom("q = 0 -> p << q = p", { ~C.eq(q), C.eq(r, p) }, false);
-    }
-
-
     void op_constraint::activate_and(core& c, dependency const& d) {
         auto x = p, y = q;
         auto& C = c.cs();
 
-        c.add_axiom("band-mask p&q <= p", { C.ule(r, p) }, false);
-        c.add_axiom("band-mask p&q <= q", { C.ule(r, q) }, false);
+        c.add_axiom("p & q <= p", { C.ule(r, p) }, false);
+        c.add_axiom("p & q <= q", { C.ule(r, q) }, false);
+        c.add_axiom("p = q -> p & q = p", { ~C.eq(p, q), C.eq(r, p) }, false);
 
         if (x.is_val())
             std::swap(x, y);
@@ -338,15 +344,15 @@ namespace polysat {
         if (!(yv + 1).is_power_of_two())
             return;
         if (yv == m.max_value())
-            c.add_axiom("band-mask-true", { C.eq(x, r) }, false);
+            c.add_axiom("p & 1 = p", { C.eq(x, r) }, false);
         else if (yv == 0)
-            c.add_axiom("band-mask-false", { C.eq(r) }, false);
+            c.add_axiom("p & 0 = 0", { 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_axiom("band-mask 1", { C.eq(x * exp, r * exp) }, false);
+            c.add_axiom("(p & 0011)*2^k = p*2^k", { C.eq(x * exp, r * exp) }, false);
         }
     }   
 
@@ -409,10 +415,10 @@ namespace polysat {
         auto& C = c.cs();
 
         if (qv.is_val() && qv.val() >= N && rv.is_val() && !rv.is_zero())            
-            c.add_axiom("q >= N  ->  r = 0", { ~C.ule(N, q), C.eq(r) }, true);
+            c.add_axiom("q >= N  ->  p >> q = 0", { ~C.ule(N, q), C.eq(r) }, true);
         else if (qv.is_zero() && pv.is_val() && rv.is_val() && rv != pv)
             // 
-            c.add_axiom("q = 0  ->  r = p", { ~C.eq(q), C.eq(r, p) }, true);
+            c.add_axiom("q = 0  ->  p >> q = p", { ~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)
@@ -423,27 +429,25 @@ namespace polysat {
             // 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)) {
-                    c.add_axiom("q = k  ->  r[i+k] = p[i] for 0 <= i < N - k", { ~C.eq(q, k), ~C.bit(r, i + k), C.bit(p, i) }, true);
+                    c.add_axiom("q = k  ->  p>>q[i+k] = p[i] for 0 <= i < N - k", { ~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)) {
-                    c.add_axiom("q = k  ->  r[i+k] = p[i] for 0 <= i < N - k", { ~C.eq(q, k), C.bit(r, i + k), ~C.bit(p, i) }, true);
+                    c.add_axiom("q = k  ->  p>>q[i+k] = p[i] for 0 <= i < N - k", { ~C.eq(q, k), C.bit(r, i + k), ~C.bit(p, i) }, true);
                 }
             }
         }
         else {
             // forward propagation
             SASSERT(!(pv.is_val() && qv.is_val() && rv.is_val()));
-            // LOG(p << " = " << pv << " and " << q << " = " << qv << " yields [<<] " << r << " = " << (qv.val().is_unsigned() ? rational::power_of_two(qv.val().get_unsigned()) * pv.val() : rational::zero()));
             if (qv.is_val() && !rv.is_val()) {
                 rational const& q_val = qv.val();
-                if (q_val >= N)
-                    // q >= N ==> r = 0
-                    c.add_axiom("shl forward 1", {~C.ule(N, q), C.eq(r)}, true);
+                if (q_val >= N)                    
+                    c.add_axiom("q >= N ==> p << q = 0", {~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());
-                    c.add_axiom("shl forward 2", {~C.eq(p, pv), ~C.eq(q, qv), C.eq(r, r_val)}, true);
+                    c.add_axiom("p = v1, q = v2, p << q -> v1 << v2", {~C.eq(p, pv), ~C.eq(q, qv), C.eq(r, r_val)}, true);
                 }
             }
         }
@@ -470,164 +474,54 @@ namespace polysat {
         auto& C = c.cs();
 
         if (pv.is_val() && rv.is_val() && rv.val() > pv.val())
-            c.add_axiom("p&q <= p", { C.ule(r, p) }, true);
+            c.add_axiom("p & q <= p", { C.ule(r, p) }, true);
         else if (qv.is_val() && rv.is_val() && rv.val() > qv.val())
-            c.add_axiom("p&q <= q", { C.ule(r, q) }, true);
+            c.add_axiom("p & q <= q", { C.ule(r, q) }, true);
         else if (pv.is_val() && qv.is_val() && rv.is_val() && pv == qv && rv != pv)
-            c.add_axiom("p = q => r = p", { ~C.eq(p, q), C.eq(r, p) }, true);
+            c.add_axiom("p = q => p & q = p", { ~C.eq(p, q), C.eq(r, p) }, true);
+        // p = a && q = b ==> r = a & b
+        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
+            c.add_axiom("p = a & q = b => r = a&b", { ~C.eq(p, pv), ~C.eq(q, qv), C.eq(r, bitwise_and(pv.val(), qv.val())) }, true);        
         else if (pv.is_val() && qv.is_val() && rv.is_val()) {
             if (pv.is_max() && qv != rv)
-                c.add_axiom("p = -1 => r = q", { ~C.eq(p, m.max_value()), C.eq(q, r) }, true);
-            if (qv.is_max() && pv != rv)
-                c.add_axiom("q = -1 => r = p", { ~C.eq(q, m.max_value()), C.eq(p, r) }, true);
+                c.add_axiom("p = -1 => p & q = q", { ~C.eq(p, m.max_value()), C.eq(q, r) }, true);
+            else if (qv.is_max() && pv != rv)
+                c.add_axiom("q = -1 => p & q = p", { ~C.eq(q, m.max_value()), C.eq(p, r) }, true);
+            else {
+                unsigned const N = m.power_of_2();
+                unsigned pow;
+                if ((pv.val() + 1).is_power_of_two(pow)) {
+                    if (rv.is_zero() && !qv.is_zero() && qv.val() <= pv.val())
+                        c.add_axiom("p = 2^k - 1 && p & q = 0 && q != 0 => q >= 2^k", { ~C.eq(p, pv), ~C.eq(r), C.eq(q), C.ule(pv + 1, q) }, true);
+                    else if (rv != qv)
+                        c.add_axiom("p = 2^k - 1  ==>  (p&q)*2^{N - k} = q*2^{N - k}", { ~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)) {
+                    if (rv.is_zero() && !pv.is_zero() && pv.val() <= qv.val())
+                        c.add_axiom("q = 2^k - 1 && p & q = 0 && p != 0  ==>  p >= 2^k", { ~C.eq(q, qv), ~C.eq(r), C.eq(p), C.ule(qv + 1, p) }, true);
+                    else if (rv != pv)
+                        c.add_axiom("q = 2^k - 1  ==>  (p&q)*2^{N - k} = p*2^{N - k}", { ~C.eq(q, qv), C.eq(r * rational::power_of_two(N - pow), p * rational::power_of_two(N - pow)) }, true);
+                }
 
-            unsigned const N = m.power_of_2();
-            unsigned pow;
-            if ((pv.val() + 1).is_power_of_two(pow)) {
-                if (rv.is_zero() && !qv.is_zero() && qv.val() <= pv.val())
-                    c.add_axiom("p = 2^k - 1 && r = 0 && q != 0 => q >= 2^k", { ~C.eq(p, pv), ~C.eq(r), C.eq(q), C.ule(pv + 1, q) }, true);
-                if (rv != qv)
-                    c.add_axiom("p = 2^k - 1  ==>  r*2^{N - k} = q*2^{N - k}", { ~C.eq(p, pv), C.eq(r * rational::power_of_two(N - pow), q * rational::power_of_two(N - pow)) }, true);
+                for (unsigned i = 0; i < N; ++i) {
+                    bool pb = pv.val().get_bit(i);
+                    bool qb = qv.val().get_bit(i);
+                    bool rb = rv.val().get_bit(i);
+                    if (rb == (pb && qb))
+                        continue;
+                    if (pb && qb && !rb)
+                        c.add_axiom("p&q[i] = p[i]&q[i]", { ~C.bit(p, i), ~C.bit(q, i), C.bit(r, i) }, true);
+                    else if (!pb && rb)
+                        c.add_axiom("p&q[i] = p[i]&q[i]", { C.bit(p, i), ~C.bit(r, i) }, true);
+                    else if (!qb && rb)
+                        c.add_axiom("p&q[i] = p[i]&q[i]", { C.bit(q, i), ~C.bit(r, i) }, true);
+                    else
+                        UNREACHABLE();
+                    return;
+                }
             }
-            if ((qv.val() + 1).is_power_of_two(pow)) {
-                if (rv.is_zero() && !pv.is_zero() && pv.val() <= qv.val())
-                    c.add_axiom("q = 2^k - 1 && r = 0 && p != 0  ==>  p >= 2^k", { ~C.eq(q, qv), ~C.eq(r), C.eq(p), C.ule(qv + 1, p) }, true);
-                // 
-                if (rv != pv)
-                    c.add_axiom("q = 2^k - 1  ==>  r*2^{N - k} = p*2^{N - k}", { ~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) {
-                bool pb = pv.val().get_bit(i);
-                bool qb = qv.val().get_bit(i);
-                bool rb = rv.val().get_bit(i);
-                if (rb == (pb && qb))
-                    continue;
-                if (pb && qb && !rb)
-                    c.add_axiom("p&q[i] = p[i]&q[i]", { ~C.bit(p, i), ~C.bit(q, i), C.bit(r, i) }, true);
-                else if (!pb && rb)
-                    c.add_axiom("p&q[i] = p[i]&q[i]", { C.bit(p, i), ~C.bit(r, i) }, true);
-                else if (!qb && rb)
-                    c.add_axiom("p&q[i] = p[i]&q[i]", { C.bit(q, i), ~C.bit(r, i) }, true);
-                else
-                    UNREACHABLE();
-            }
-            return;
         }
 
-        // Propagate r if p or q are 0
-        else if (pv.is_zero() && !rv.is_zero())  // rv not necessarily fully evaluated
-            c.add_axiom("p = 0 -> p&q = 0", { C.ule(r, p) }, true);
-        else if (qv.is_zero() && !rv.is_zero())  // rv not necessarily fully evaluated
-            c.add_axiom("q = 0 -> p&q = 0", { C.ule(r, q) }, true);
-        // p = a && q = b ==> r = a & b
-        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()));
-            c.add_axiom("p = a & q = b => r = a&b", { ~C.eq(p, pv), ~C.eq(q, qv), C.eq(r, bitwise_and(pv.val(), qv.val())) }, true);
-        }
     }
-
-
-#if 0
-
-    // introduce multiplication constraint and do away with non-linear polynomials in inequalities.
-    // 
-    // z = x * y
-    // x = 0 or y = 0 => z = 0
-    // x = 1 => z = y
-    // y = 1 => z = y
-    // ~ovfl(x, y) => z >= x & z >= y
-    // ~ovfl(x, y) & x > 1 & y > 1 => z > x, z > y
-    // parity(x) + parity(y) >= N => z = 0
-    // parity(x) + parity(y) < N => parity(z) = parity(x) + parity(y)
-    // blast:
-    // z = sum_i bit(x,i) ? y*2^i : 0
-
-#endif
-
-#if 0
-
-    /**
-     * Produce lemmas for constraint: r == inv p
-     * p = 0   ==>  r = 0
-     * r = 0   ==>  p = 0
-     * p != 0  ==>  odd(r)
-     * parity(p) >= k   ==>  p * r >= 2^k
-     * parity(p) < k    ==>  p * r <= 2^k - 1
-     * parity(p) < k    ==>  r <= 2^(N - k) - 1     (because r is the smallest pseudo-inverse)
-     */
-    clause_ref op_constraint::lemma_inv(solver& s, assignment const& a) {
-        auto& m = p.manager();
-        auto pv = a.apply_to(p);
-        auto rv = a.apply_to(r);
-
-        if (eval_inv(pv, rv) == l_true)
-            return {};
-
-        signed_constraint const invc(this, true);
-
-        // p = 0  ==>  r = 0
-        if (pv.is_zero())
-            c.add_axiom(~invc, ~C.eq(p), C.eq(r), true);
-        // r = 0  ==>  p = 0
-        if (rv.is_zero())
-            c.add_axiom(~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())
-            c.add_axiom(~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 {};
-
-        unsigned parity_pv = pv.val().trailing_zeros();
-        unsigned parity_rv = rv.val().trailing_zeros();
-
-        LOG("p: " << p << " := " << pv << "   parity " << parity_pv);
-        LOG("r: " << r << " := " << rv << "   parity " << parity_rv);
-
-        // p != 0  ==>  odd(r)
-        if (parity_rv != 0)
-            c.add_axiom("r = inv p  &  p != 0  ==>  odd(r)", {~invc, C.eq(p), s.odd(r)}, true);
-
-        pdd prod = p * r;
-        rational prodv = (pv * rv).val();
-//        if (prodv != rational::power_of_two(parity_pv))
-//            verbose_stream() << prodv << " " << rational::power_of_two(parity_pv) << " " << parity_pv << " " << pv << " " << rv << "\n";
-        unsigned lower = 0, upper = m.power_of_2();
-        // binary search for the parity (otw. we would have justifications like "parity_at_most(k) && parity_at_least(k)" for at most "k" widths
-        while (lower + 1 < upper) {
-            unsigned middle = (upper + lower) / 2;
-            LOG("Splitting on " << middle);
-            if (parity_pv >= middle) { // parity at least middle
-                lower = middle;
-                LOG("Its in [" << lower << "; " << upper << ")");
-                // parity(p) >= k  ==>  p * r >= 2^k
-                if (prodv < rational::power_of_two(middle))
-                    c.add_axiom("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)
-                    c.add_axiom("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);
-                upper = middle;
-                LOG("Its in [" << lower << "; " << upper << ")");
-                // parity(p) < k   ==>  p * r <= 2^k - 1
-                if (prodv > rational::power_of_two(middle))
-                    c.add_axiom("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?
-        UNREACHABLE();
-        return {};
-    }
-
-#endif
 }
diff --git a/src/sat/smt/polysat/op_constraint.h b/src/sat/smt/polysat/op_constraint.h
index 4b44b3009..c5400fbab 100644
--- a/src/sat/smt/polysat/op_constraint.h
+++ b/src/sat/smt/polysat/op_constraint.h
@@ -68,7 +68,6 @@ namespace polysat {
         std::ostream& display(std::ostream& out, char const* eq) const;
 
         void activate_and(core& s, dependency const& d);
-        void activate_ashr(core& s, dependency const& d);
 
     public:
         ~op_constraint() override {}