diff --git a/src/sat/smt/CMakeLists.txt b/src/sat/smt/CMakeLists.txt
index 2a6fb9e66..e27849391 100644
--- a/src/sat/smt/CMakeLists.txt
+++ b/src/sat/smt/CMakeLists.txt
@@ -39,6 +39,7 @@ z3_add_component(sat_smt
     polysat_internalize.cpp
     polysat_model.cpp
     polysat_solver.cpp
+    polysat_ule.cpp
     q_clause.cpp
     q_ematch.cpp
     q_eval.cpp
diff --git a/src/sat/smt/polysat_constraints.cpp b/src/sat/smt/polysat_constraints.cpp
index 1c9de327c..101d60e62 100644
--- a/src/sat/smt/polysat_constraints.cpp
+++ b/src/sat/smt/polysat_constraints.cpp
@@ -11,15 +11,21 @@ Author:
     Jakob Rath 2021-04-06
 
 --*/
+
 #include "sat/smt/polysat_core.h"
 #include "sat/smt/polysat_solver.h"
 #include "sat/smt/polysat_constraints.h"
+#include "sat/smt/polysat_ule.h"
 
 namespace polysat {
 
     signed_constraint constraints::ule(pdd const& p, pdd const& q) {
+        pdd lhs = p, rhs = q;
+        bool is_positive = true;
+        ule_constraint::simplify(is_positive, lhs, rhs);
         auto* c = alloc(ule_constraint, p, q);
         m_trail.push(new_obj_trail(c));
-        return signed_constraint(ckind_t::ule_t, c);
+        auto sc = signed_constraint(ckind_t::ule_t, c);
+        return is_positive ? sc : ~sc;
     }
 }
diff --git a/src/sat/smt/polysat_constraints.h b/src/sat/smt/polysat_constraints.h
index b8068c68a..da82431c4 100644
--- a/src/sat/smt/polysat_constraints.h
+++ b/src/sat/smt/polysat_constraints.h
@@ -14,11 +14,14 @@ Author:
 
 
 #pragma once 
+#include "util/trail.h"
 #include "sat/smt/polysat_types.h"
 
 namespace polysat {
 
     class core;
+    class ule_constraint;
+    class assignment;
 
     using pdd = dd::pdd;
     using pvar = unsigned;
@@ -33,15 +36,12 @@ namespace polysat {
         unsigned_vector const& vars() const { return m_vars; }
         unsigned var(unsigned idx) const { return m_vars[idx]; }
         bool contains_var(pvar v) const { return m_vars.contains(v); }
+        virtual std::ostream& display(std::ostream& out, lbool status) const = 0;
+        virtual std::ostream& display(std::ostream& out) const = 0;
+        virtual lbool eval() const = 0;
+        virtual lbool eval(assignment const& a) const = 0;
     };
 
-    class ule_constraint : public constraint {
-        pdd m_lhs, m_rhs;
-    public:
-        ule_constraint(pdd const& lhs, pdd const& rhs) : m_lhs(lhs), m_rhs(rhs) {}
-        pdd const& lhs() const { return m_lhs; }
-        pdd const& rhs() const { return m_rhs; }
-    };
 
     class umul_ovfl_constraint : public constraint {
         pdd m_lhs, m_rhs;
diff --git a/src/sat/smt/polysat_solver.cpp b/src/sat/smt/polysat_solver.cpp
index 65183578d..57098c447 100644
--- a/src/sat/smt/polysat_solver.cpp
+++ b/src/sat/smt/polysat_solver.cpp
@@ -27,6 +27,7 @@ The result of polysat::core::check is one of:
 #include "ast/euf/euf_bv_plugin.h"
 #include "sat/smt/polysat_solver.h"
 #include "sat/smt/euf_solver.h"
+#include "sat/smt/polysat_ule.h"
 
 
 namespace polysat {
@@ -170,11 +171,12 @@ namespace polysat {
     // Core uses the propagate callback to add unit propagations to the trail.
     // The polysat::solver takes care of translating signed constraints into expressions, which translate into literals.
     // Everything goes over expressions/literals. polysat::core is not responsible for replaying expressions. 
-    void solver::propagate(signed_constraint sc, dependency_vector const& deps) {
+    dependency solver::propagate(signed_constraint sc, dependency_vector const& deps) {
         sat::literal lit = ctx.mk_literal(constraint2expr(sc));
         auto [core, eqs] = explain_deps(deps);
         auto ex = euf::th_explain::propagate(*this, core, eqs, lit, nullptr);
         ctx.propagate(lit, ex);
+        return dependency(lit, s().lvl(lit));
     }
 
     void solver::propagate(dependency const& d, bool sign, dependency_vector const& deps) {
diff --git a/src/sat/smt/polysat_solver.h b/src/sat/smt/polysat_solver.h
index 408d0756c..76923f88f 100644
--- a/src/sat/smt/polysat_solver.h
+++ b/src/sat/smt/polysat_solver.h
@@ -131,7 +131,7 @@ namespace polysat {
         void add_eq_literal(pvar v, rational const& val);
         void set_conflict(dependency_vector const& core);
         void set_lemma(vector<signed_constraint> const& lemma, unsigned level, dependency_vector const& core);
-        void propagate(signed_constraint sc, dependency_vector const& deps);
+        dependency propagate(signed_constraint sc, dependency_vector const& deps);
         void propagate(dependency const& d, bool sign, dependency_vector const& deps);
         
         void add_lemma(vector<signed_constraint> const& lemma);
diff --git a/src/sat/smt/polysat_ule.cpp b/src/sat/smt/polysat_ule.cpp
new file mode 100644
index 000000000..08448b34d
--- /dev/null
+++ b/src/sat/smt/polysat_ule.cpp
@@ -0,0 +1,346 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    polysat unsigned <= constraints
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2021-03-19
+    Jakob Rath 2021-04-06
+
+Notes:
+
+   Canonical representation of equation p == 0 is the constraint p <= 0.
+   The alternatives p < 1, -1 <= q, q > -2 are eliminated.
+
+   Rewrite rules to simplify expressions.
+   In the following let k, k1, k2 be values.
+
+   - k1 <= k2       ==>   0 <= 0 if k1 <= k2
+   - k1 <= k2       ==>   1 <= 0 if k1 >  k2
+   - 0 <= p         ==>   0 <= 0
+   - p <= 0         ==>   1 <= 0 if p is never zero due to parity
+   - p <= -1        ==>   0 <= 0
+   - k <= p         ==>   p - k <= - k - 1
+   - k*2^n*p <= 0   ==>   2^n*p <= 0 if k is odd, leading coeffient is always a power of 2.
+
+   Note: the rules will rewrite alternative formulations of equations:
+   - -1 <= p        ==>   p + 1 <= 0
+   -  1 <= p        ==>   p - 1 <= -2
+
+   Rewrite rules on signed constraints:
+
+   - p > -2         ==>   p + 1 <= 0
+   - p <= -2        ==>   p + 1 > 0
+
+   At this point, all equations are in canonical form.
+
+TODO: clause simplifications:
+
+   - p + k <= p  ==> p + k <= k & p != 0   for k != 0
+   - p*q = 0     ==> p = 0 or q = 0        applies to any factoring
+   - 2*p <= 2*q  ==> (p >= 2^n-1 & q < 2^n-1) or (p >= 2^n-1 = q >= 2^n-1 & p <= q)
+                 ==> (p >= 2^n-1 => q < 2^n-1 or p <= q) &
+                     (p < 2^n-1 => p <= q) &
+                     (p < 2^n-1 => q < 2^n-1)
+
+   - 3*p <= 3*q ==> ?
+
+Note:
+     case p <= p + k is already covered because we test (lhs - rhs).is_val()
+
+     It can be seen as an instance of lemma 5.2 of Supratik and John.
+
+The following forms are equivalent:
+
+    p <= q
+    p <= p - q - 1
+    q - p <= q
+    q - p <= -p - 1
+    -q - 1 <= -p - 1
+    -q - 1 <= p - q - 1
+
+Useful lemmas:
+
+    p <= q  &&  q+1 != 0    ==>  p+1 <= q+1
+
+    p <= q  &&  p != 0      ==>   -q <= -p
+
+--*/
+
+#include "sat/smt/polysat_constraints.h"
+#include "sat/smt/polysat_ule.h"
+
+#define LOG(_msg_) verbose_stream() << _msg_ << "\n"
+
+namespace polysat {
+
+    // Simplify lhs <= rhs.
+    //
+    // NOTE: the result should not depend on the initial value of is_positive;
+    //       the purpose of is_positive is to allow flipping the sign as part of a rewrite rule.
+    static void simplify_impl(bool& is_positive, pdd& lhs, pdd& rhs) {
+        
+        SASSERT_EQ(lhs.power_of_2(), rhs.power_of_2());
+        unsigned const N = lhs.power_of_2();
+
+        // 0 <= p   -->   0 <= 0
+        if (lhs.is_zero()) {
+            rhs = 0;
+            return;
+        }
+        // p <= -1   -->   0 <= 0
+        if (rhs.is_max()) {
+            lhs = 0, rhs = 0;
+            return;
+        }
+        // p <= p   -->   0 <= 0
+        if (lhs == rhs) {
+            lhs = 0, rhs = 0;
+            return;
+        }
+        // Evaluate constants
+        if (lhs.is_val() && rhs.is_val()) {
+            if (lhs.val() <= rhs.val())
+                lhs = 0, rhs = 0;
+            else
+                lhs = 0, rhs = 0, is_positive = !is_positive;
+            return;
+        }
+
+        // Try to reduce the number of variables on one side using one of these rules:
+        //
+        //      p <= q  -->  p <= p - q - 1
+        //      p <= q  -->  q - p <= q
+        //
+        // Possible alternative to 1:
+        //      p <= q  -->  q - p <= -p - 1
+        // Possible alternative to 2:
+        //      p <= q  -->  -q-1 <= p - q - 1
+        //
+        // Example:
+        //
+        //      x <= x + y  -->  x <= - y - 1
+        //      x + y <= x  -->  -y <= x
+        if (!lhs.is_val() && !rhs.is_val()) {
+            unsigned const lhs_vars = lhs.free_vars().size();
+            unsigned const rhs_vars = rhs.free_vars().size();
+            unsigned const diff_vars = (lhs - rhs).free_vars().size();
+            if (diff_vars < lhs_vars || diff_vars < rhs_vars) {
+                LOG("reduce number of varables");
+                // verbose_stream() << "IN:  " << ule_pp(to_lbool(is_positive), lhs, rhs) << "\n";
+                if (lhs_vars <= rhs_vars)
+                    rhs = lhs - rhs - 1;
+                else
+                    lhs = rhs - lhs;
+                // verbose_stream() << "OUT: " << ule_pp(to_lbool(is_positive), lhs, rhs) << "\n";
+            }
+        }
+
+        // -p + k <= k      -->  p <= k
+        if (rhs.is_val() && !rhs.is_zero() && lhs.offset() == rhs.val()) {
+            LOG("-p + k <= k      -->  p <= k");
+            lhs = rhs - lhs;
+        }
+
+        // k <= p + k       -->  p <= -k-1
+        if (lhs.is_val() && !lhs.is_zero() && lhs.val() == rhs.offset()) {
+            LOG("k <= p + k       -->  p <= -k-1");
+            pdd k = lhs;
+            lhs = rhs - lhs;
+            rhs = -k - 1;
+        }
+
+        // k <= -p          -->  p-1 <= -k-1
+        if (lhs.is_val() && rhs.leading_coefficient().get_bit(N - 1) && !rhs.offset().is_zero()) {
+            LOG("k <= -p          -->  p-1 <= -k-1");
+            pdd k = lhs;
+            lhs = -(rhs + 1);
+            rhs = -k - 1;
+        }
+
+        // -p <= k          -->  -k-1 <= p-1
+        // if (rhs.is_val() && lhs.leading_coefficient() > rational::power_of_two(N - 1) && !lhs.offset().is_zero()) {
+        if (rhs.is_val() && lhs.leading_coefficient().get_bit(N - 1) && !lhs.offset().is_zero()) {
+            LOG("-p <= k          -->  -k-1 <= p-1");
+            pdd k = rhs;
+            rhs = -(lhs + 1);
+            lhs = -k - 1;
+        }
+
+        // NOTE: do not use pdd operations in conditions when comparing pdd values.
+        //       e.g.: "lhs.offset() == (rhs + 1).val()" is problematic with the following evaluation:
+        //          1. return reference into pdd_manager::m_values from lhs.offset()
+        //          2. compute rhs+1, which may reallocate pdd_manager::m_values
+        //          3. now the reference returned from lhs.offset() may be invalid
+        pdd const rhs_plus_one = rhs + 1;
+
+        // p - k <= -k - 1  -->  k <= p
+        // TODO: potential bug here: first call offset(), then rhs+1 has to reallocate pdd_manager::m_values, then the reference to offset is broken.
+        if (rhs.is_val() && !rhs.is_zero() && lhs.offset() == rhs_plus_one.val()) {
+            LOG("p - k <= -k - 1  -->  k <= p");
+            pdd k = -(rhs + 1);
+            rhs = lhs + k;
+            lhs = k;
+        }
+
+        pdd const lhs_minus_one = lhs - 1;
+
+        // k <= 2^(N-1)*p + q + k-1  -->  k <= 2^(N-1)*p - q
+        if (lhs.is_val() && rhs.leading_coefficient() == rational::power_of_two(N-1) && rhs.offset() == lhs_minus_one.val()) {
+            LOG("k <= 2^(N-1)*p + q + k-1  -->  k <= 2^(N-1)*p - q");
+            rhs = (lhs - 1) - rhs;
+        }
+
+        // -1 <= p  -->   p + 1 <= 0
+        if (lhs.is_max()) {
+            lhs = rhs + 1;
+            rhs = 0;
+        }
+
+        // 1 <= p   -->   p > 0
+        if (lhs.is_one()) {
+            lhs = rhs;
+            rhs = 0;
+            is_positive = !is_positive;
+        }
+
+        // p >  -2   -->   p + 1 <= 0
+        // p <= -2   -->   p + 1 >  0
+        if (rhs.is_val() && !rhs.is_zero() && (rhs + 2).is_zero()) {
+            // Note: rhs.is_zero() iff rhs.manager().power_of_2() == 1 (the rewrite is not wrong for M=2, but useless)
+            lhs = lhs + 1;
+            rhs = 0;
+            is_positive = !is_positive;
+        }
+        // 2p + 1 <= 0   -->   0 < 0
+        if (rhs.is_zero() && lhs.is_never_zero()) {
+            lhs = 0;
+            is_positive = !is_positive;
+            return;
+        }
+        // a*p + q <= 0   -->   p + a^-1*q <= 0   for a odd
+        if (rhs.is_zero() && !lhs.leading_coefficient().is_power_of_two()) {
+            rational lc = lhs.leading_coefficient();
+            rational x, y;
+            gcd(lc, lhs.manager().two_to_N(), x, y);
+            if (x.is_neg())
+                x = mod(x, lhs.manager().two_to_N());
+            lhs *= x;
+            SASSERT(lhs.leading_coefficient().is_power_of_two());
+        }
+    } // simplify_impl
+}
+
+
+namespace polysat {
+
+    ule_constraint::ule_constraint(pdd const& l, pdd const& r) :
+        m_lhs(l), m_rhs(r) {
+
+        SASSERT_EQ(m_lhs.power_of_2(), m_rhs.power_of_2());
+
+        vars().append(m_lhs.free_vars());
+        for (auto v : m_rhs.free_vars())
+            if (!vars().contains(v))
+                vars().push_back(v);
+    }
+
+    void ule_constraint::simplify(bool& is_positive, pdd& lhs, pdd& rhs) {
+        SASSERT_EQ(lhs.power_of_2(), rhs.power_of_2());
+#ifndef NDEBUG
+        bool const old_is_positive = is_positive;
+        pdd const old_lhs = lhs;
+        pdd const old_rhs = rhs;
+#endif
+        simplify_impl(is_positive, lhs, rhs);
+#ifndef NDEBUG
+        if (old_is_positive != is_positive || old_lhs != lhs || old_rhs != rhs) {
+            ule_pp const old_ule(to_lbool(old_is_positive), old_lhs, old_rhs);
+            ule_pp const new_ule(to_lbool(is_positive), lhs, rhs);
+            // always-false and always-true constraints should be rewritten to 0 != 0 and 0 == 0, respectively.
+            if (is_always_false(old_is_positive, old_lhs, old_rhs)) {
+                SASSERT(!is_positive);
+                SASSERT(lhs.is_zero());
+                SASSERT(rhs.is_zero());
+            }
+            if (is_always_true(old_is_positive, old_lhs, old_rhs)) {
+                SASSERT(is_positive);
+                SASSERT(lhs.is_zero());
+                SASSERT(rhs.is_zero());
+            }
+        }
+        SASSERT(is_simplified(lhs, rhs));  // rewriting should be idempotent
+#endif
+    }
+
+    bool ule_constraint::is_simplified(pdd const& lhs0, pdd const& rhs0) {
+        bool const pos0 = true;
+        bool pos1 = pos0;
+        pdd lhs1 = lhs0;
+        pdd rhs1 = rhs0;
+        simplify_impl(pos1, lhs1, rhs1);
+        bool const is_simplified = (pos1 == pos0 && lhs1 == lhs0 && rhs1 == rhs0);
+        DEBUG_CODE({
+            // check that simplification doesn't depend on initial sign
+            bool pos2 = !pos0;
+            pdd lhs2 = lhs0;
+            pdd rhs2 = rhs0;
+            simplify_impl(pos2, lhs2, rhs2);
+            SASSERT_EQ(pos2, !pos1);
+            SASSERT_EQ(lhs2, lhs1);
+            SASSERT_EQ(rhs2, rhs1);
+        });
+        return is_simplified;
+    }
+
+    std::ostream& ule_constraint::display(std::ostream& out, lbool status, pdd const& lhs, pdd const& rhs) {
+        out << lhs;
+        if (rhs.is_zero() && status == l_true) out << " == ";
+        else if (rhs.is_zero() && status == l_false) out << " != ";
+        else if (status == l_true) out << " <= ";
+        else if (status == l_false) out << " > ";
+        else out << " <=/> ";
+        return out << rhs;
+    }
+
+    std::ostream& ule_constraint::display(std::ostream& out, lbool status) const {
+        return display(out, status, m_lhs, m_rhs);
+    }
+
+    std::ostream& ule_constraint::display(std::ostream& out) const {
+        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
+        if (lhs.is_zero())
+            return l_true;      // 0 <= p
+        if (lhs == rhs)
+            return l_true;      // p <= p
+        if (rhs.is_max())
+            return l_true;      // p <= -1
+        if (rhs.is_zero() && lhs.is_never_zero())
+            return l_false;     // p <= 0 implies p == 0
+        if (lhs.is_one() && rhs.is_never_zero())
+            return l_true;      // 1 <= p implies p != 0
+        if (lhs.is_val() && rhs.is_val())
+            return to_lbool(lhs.val() <= rhs.val());
+        return l_undef;
+
+   }
+
+    lbool ule_constraint::eval() const {
+        return eval(lhs(), rhs());
+    }
+
+    lbool ule_constraint::eval(assignment const& a) const {
+        return eval(a.apply_to(lhs()), a.apply_to(rhs()));
+    }
+
+}
diff --git a/src/sat/smt/polysat_ule.h b/src/sat/smt/polysat_ule.h
new file mode 100644
index 000000000..12efe506a
--- /dev/null
+++ b/src/sat/smt/polysat_ule.h
@@ -0,0 +1,56 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    polysat unsigned <= constraint
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2021-03-19
+    Jakob Rath 2021-04-06
+
+--*/
+#pragma once
+#include "sat/smt/polysat_ule.h"
+#include "sat/smt/polysat_assignment.h"
+#include "sat/smt/polysat_constraints.h"
+
+
+namespace polysat {
+
+    class ule_constraint final : public constraint {
+        pdd m_lhs;
+        pdd m_rhs;        
+        static bool is_always_true(bool is_positive, pdd const& lhs, pdd const& rhs) { return eval(lhs, rhs) == to_lbool(is_positive); }
+        static bool is_always_false(bool is_positive, pdd const& lhs, pdd const& rhs) { return is_always_true(!is_positive, lhs, rhs); }
+        static lbool eval(pdd const& lhs, pdd const& rhs);
+
+    public:
+        ule_constraint(pdd const& l, pdd const& r);
+        ~ule_constraint() override {}
+        pdd const& lhs() const { return m_lhs; }
+        pdd const& rhs() const { return m_rhs; }
+        static std::ostream& display(std::ostream& out, lbool status, pdd const& lhs, pdd const& rhs);
+        std::ostream& display(std::ostream& out, lbool status) const override;
+        std::ostream& display(std::ostream& out) const override;
+        lbool eval() const override;
+        lbool eval(assignment const& a) const override;
+        bool is_eq() const { return m_rhs.is_zero(); }
+        unsigned power_of_2() const { return m_lhs.power_of_2(); }
+
+        static void simplify(bool& is_positive, pdd& lhs, pdd& rhs);
+        static bool is_simplified(pdd const& lhs, pdd const& rhs);  // return true if lhs <= rhs is not simplified further. this is meant to be used in assertions.
+    };
+
+    struct ule_pp {
+        lbool status;
+        pdd lhs;
+        pdd rhs;
+        ule_pp(lbool status, pdd const& lhs, pdd const& rhs): status(status), lhs(lhs), rhs(rhs) {}
+        ule_pp(lbool status, ule_constraint const& ule): status(status), lhs(ule.lhs()), rhs(ule.rhs()) {}
+    };
+
+    inline std::ostream& operator<<(std::ostream& out, ule_pp const& u) { return ule_constraint::display(out, u.status, u.lhs, u.rhs); }
+
+}