diff --git a/src/math/polysat/constraint.cpp b/src/math/polysat/constraint.cpp
index d9e5a55db..1e19f5284 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/shr_constraint.h"
 
 namespace polysat {
 
@@ -232,6 +233,10 @@ namespace polysat {
         return { dedup(alloc(mul_ovfl_constraint, *this, a, b)), true };
     }
 
+    signed_constraint constraint_manager::shr(pdd const& p, pdd const& q) {
+        return { dedup(alloc(shr_constraint, *this, p, q)), true };
+    }
+
 
     // To do signed comparison of bitvectors, flip the msb and do unsigned comparison:
     //
@@ -279,6 +284,14 @@ namespace polysat {
         return *dynamic_cast<mul_ovfl_constraint const*>(this);
     }
 
+    shr_constraint& constraint::to_shr() {
+        return *dynamic_cast<shr_constraint*>(this);
+    }
+
+    shr_constraint const& constraint::to_shr() const {
+        return *dynamic_cast<shr_constraint const*>(this);
+    }
+
     std::string constraint::bvar2string() const {
         std::stringstream out;
         out << " (b";
diff --git a/src/math/polysat/constraint.h b/src/math/polysat/constraint.h
index b58efcc46..b62e3b3f2 100644
--- a/src/math/polysat/constraint.h
+++ b/src/math/polysat/constraint.h
@@ -97,6 +97,7 @@ namespace polysat {
         signed_constraint slt(pdd const& a, pdd const& b);
         signed_constraint mul_ovfl(pdd const& p, pdd const& q);
         signed_constraint bit(pdd const& p, unsigned i);
+        signed_constraint shr(pdd const& p, pdd const& q);
 
         constraint *const* begin() const { return m_constraints.data(); }
         constraint *const* end() const { return m_constraints.data() + m_constraints.size(); }
@@ -164,6 +165,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_shr() const { return m_kind == ckind_t::shr_t; }
         ckind_t kind() const { return m_kind; }
         virtual std::ostream& display(std::ostream& out, lbool status) const = 0;
         virtual std::ostream& display(std::ostream& out) const = 0;
@@ -180,6 +182,8 @@ namespace polysat {
         ule_constraint const& to_ule() const;
         mul_ovfl_constraint& to_mul_ovfl();
         mul_ovfl_constraint const& to_mul_ovfl() const;
+        shr_constraint& to_shr();
+        shr_constraint const& to_shr() const;
         unsigned_vector& vars() { return m_vars; }
         unsigned_vector const& vars() const { return m_vars; }
         unsigned var(unsigned idx) const { return m_vars[idx]; }
diff --git a/src/math/polysat/shr_constraint.cpp b/src/math/polysat/shr_constraint.cpp
index 17e336c5a..da736d4cc 100644
--- a/src/math/polysat/shr_constraint.cpp
+++ b/src/math/polysat/shr_constraint.cpp
@@ -19,6 +19,13 @@ namespace polysat {
 
     shr_constraint::shr_constraint(constraint_manager& m, pdd const& p, pdd const& q, pdd const& r):
         constraint(m, ckind_t::shr_t), m_p(p), m_q(q), m_r(r) {
+        m_vars.append(p.free_vars());
+        for (auto v : q.free_vars())
+            if (!m_vars.contains(v))
+                m_vars.push_back(v);
+        for (auto v : r.free_vars())
+            if (!m_vars.contains(v))
+                m_vars.push_back(v);
     }    
 
     lbool shr_constraint::eval(pdd const& p, pdd const& q, pdd const& r) const {
@@ -76,7 +83,12 @@ namespace polysat {
     * r = p >> q & q >= k -> r[i] = 0 for i > K - k
     * r = p >> q & q >= K -> r = 0
     * r = p >> q & q = k -> r[i] = p[i+k] for k + i < K
-    * r = p >> q & q >= k -> r <= 2^K-k-1
+    * r = p >> q & q >= k -> r <= 2^{K-k-1}
+    * r = p >> q => r <= p
+    * r = p >> q & q != 0 => r < p
+    * r = p >> q & q = 0 => r = p
+    * 
+    * when q is a constant, several axioms can be enforced at activation time.
     * 
     * Enforce also inferences and bounds
     * 
@@ -94,7 +106,7 @@ namespace polysat {
     bool shr_constraint::operator==(constraint const& other) const {
         if (other.kind() != ckind_t::shr_t)
             return false;
-        auto const& o = dynamic_cast<shr_constraint const&>(other);
+        auto const& o = other.to_shr();
         return p() == o.p() && q() == o.q() && r() == o.r();
     }
 
diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp
index 35af5aaa6..cadd8aba9 100644
--- a/src/math/polysat/solver.cpp
+++ b/src/math/polysat/solver.cpp
@@ -135,6 +135,14 @@ namespace polysat {
         return std::tuple<pdd, pdd>(quot, rem);
     }
 
+    pdd solver::shr(pdd const& p, pdd const& q) {
+        auto& m = p.manager();
+        unsigned sz = m.power_of_2();
+        pdd r = m.mk_var(add_var(sz));
+        assign_eh(m_constraints.shr(r, p, q));
+        return r;
+    }
+
 
     void solver::assign_eh(signed_constraint c, unsigned dep) {
         SASSERT(at_base_level());
diff --git a/src/math/polysat/solver.h b/src/math/polysat/solver.h
index a3f4d7254..05096dcd0 100644
--- a/src/math/polysat/solver.h
+++ b/src/math/polysat/solver.h
@@ -267,6 +267,11 @@ namespace polysat {
         */
         std::tuple<pdd, pdd> quot_rem(pdd const& a, pdd const& b);
 
+        /**
+        * Create expression for the logical right shift of p by q.
+        */
+        pdd shr(pdd const& p, pdd const& q);
+
         /**
          * Create polynomial constant.
          */