diff --git a/src/math/polysat/saturation.cpp b/src/math/polysat/saturation.cpp
index eedaff140..2735f9b53 100644
--- a/src/math/polysat/saturation.cpp
+++ b/src/math/polysat/saturation.cpp
@@ -1441,69 +1441,10 @@ namespace polysat {
      * The range is a "forbidden interval" for y and is implied. It is much stronger than resolving on y0.
      *
      */
-
-#if 0
-    // outline of what should be a more general approach
-    
-    pdd p = a_l_b.lhs(), q = a_l_b.rhs();
-    if (p.degree(x) > 1 || q.degree(x) > 1)
-        return false;
-    if (p.degree(x) == 0 && q.degree(x) == 0)
-        return false;
-    vector<signed_constraint> bounds;
-    if (!m_viable.has_max_forbidden(x, lo_x, hi_x, bounds))
-        return false;
-    SASSERT(lo_x != hi_x);
-    if (lo_x > hi_x)
-        lo_x += m.two_to_N();
-    SASSERT(lo_x < hi_x);
-
-    auto is_bounded = [&](pdd& p, rational& lo, rational& hi) {
-        if (!lower_bound(x, lo_x, p, lo))
-            return false;
-        if (!upper_bound(x, hi_x, p, hi))
-            return false;
-        SASSERT(0 <= lo && lo <= hi);
-        if (lo + m.two_to_N() < hi)
-            return false;
-        ratinoal offset = floor(lo / m.two_to_N()) * m.two_to_N();
-        lo -= offset;
-        hi -= offset;
-        return true;
-    };
-
-    rational lo_p, hi_p, offset_p;
-    rational lo_q, hi_q, offset_q;
-    rational lo_r, hi_r, offset_r;
-    pdd r = q - p;
-    if (!is_bounded(p, lo_p, hi_p, offset_p))
-        return false;
-    if (!is_bounded(q, lo_q, hi_q, offset_q))
-        return false;
-    if (!is_bounded(r, lo_r, hi_r, offset_r))
-        return false;
-    SASSERT(0 <= lo_r && lo_r <= hi_r);
-
-    
-    // for every value of x, p, q are bewteen lo_p, hi_p, lo_q, hi_q
-
-
-    // if a_l_b is non-strict, it is false under all assignments to x, so r > 0
-    // if a_l_b is strict, we bail also if r = 1
-    if (lo_r == 0)
-        return false;
-    if (a_l_b.is_strict() && lo_r == 1)
-        return false;
-
-    // isolate what is 'y'
-    
-    // then compute range around y that is admissible based on x_lo, x_hi
-    
-    
-#endif
     
     bool saturation::try_add_mul_bound(pvar x, conflict& core, inequality const& a_l_b) {
         set_rule("[x] ax + b <= y, ... => a >= u_a");
+        // try_add_mul_bound2(x, core, a_l_b);
         auto& m = s.var2pdd(x);        
         pdd const X = s.var(x);
         pdd a = s.var(x);
@@ -1556,6 +1497,126 @@ namespace polysat {
         return false;
     }
 
+
+    bool saturation::get_bound(pvar x, rational const& bound_x, pdd const& p, rational& bound_p) {
+        if (p.degree(x) == 0)
+            return s.try_eval(p, bound_p);
+        pdd a = p, b = p;
+        rational a_val, b_val;
+        p.factor(x, 1, a, b);
+        if (!get_bound(x, bound_x, a, a_val))
+            return false;
+        if (!get_bound(x, bound_x, b, b_val))
+            return false;
+        bound_p = bound_x * a_val + b_val;
+        return true;
+    }
+
+    // wip - outline of what should be a more general approach
+    bool saturation::try_add_mul_bound2(pvar x, conflict& core, inequality const& a_l_b) {
+        auto& m = s.var2pdd(x);    
+        pdd p = a_l_b.lhs(), q = a_l_b.rhs();
+        if (p.degree(x) > 1 || q.degree(x) > 1)
+            return false;
+        if (p.degree(x) == 0 && q.degree(x) == 0)
+            return false;
+        vector<signed_constraint> bounds;
+        rational x_min, x_max;
+        if (!s.m_viable.has_max_forbidden(x, x_min, x_max, bounds))
+            return false;
+
+        VERIFY(x_min != x_max);
+        // From forbidden interval [x_min, x_max[ compute
+        // allowed range: [x_max, x_min - 1]
+        SASSERT(0 <= x_min && x_min <= m.max_value());
+        SASSERT(0 <= x_max && x_max <= m.max_value());
+        rational hi = x_min == 0 ? m.max_value() : x_min - 1; 
+        x_min = x_max;
+        x_max = hi;
+        SASSERT(x_min != x_max);
+        if (x_min > x_max)
+            x_min += m.two_to_N();
+        SASSERT(x_min <= x_max);
+        if (x_max == 0) {
+            SASSERT(x_min == 0);
+            return false;
+        }
+        
+        auto isolate_y = [&](pvar x, pdd const& p, pdd& y, pdd& b) {
+            if (p.degree(x) != 1)
+                return false;
+            p.factor(x, 1, y, b);            
+            return !y.is_val();
+        };
+        pdd y = p, b = p;
+        rational b_val, y_val;
+        // handle just one-side with x for now
+        if (p.degree(x) == 1 && q.degree(x) == 1)
+            return false;
+        if (!isolate_y(x, p, y, b) && !isolate_y(x, q, y, b))
+            return false;
+        if (!s.try_eval(b, b_val))
+            return false;
+        if (!s.try_eval(y, y_val))
+            return false;
+        // at this point p = x*y + b or q = x*y + b
+        
+        auto is_bounded = [&](pdd& p, rational& lo, rational& hi) {
+            verbose_stream() << "is-bounded " << p << "\n";
+            if (!get_bound(x, x_min, p, lo))
+                return false;
+            if (!get_bound(x, x_max, p, hi))
+                return false;
+            SASSERT(0 <= lo && lo <= hi);
+            if (lo + m.two_to_N() < hi)
+                return false;
+            rational offset = floor(lo / m.two_to_N()) * m.two_to_N();
+            lo -= offset;
+            hi -= offset;
+            return true;
+        };
+        
+        rational lo_p, hi_p;
+        rational lo_q, hi_q;
+        rational lo_r, hi_r;
+        pdd r = q - p;
+        if (!is_bounded(p, lo_p, hi_p))
+            return false;
+        if (!is_bounded(q, lo_q, hi_q))
+            return false;
+        if (!is_bounded(r, lo_r, hi_r))
+            return false;
+        SASSERT(0 <= lo_r && lo_r <= hi_r);
+
+        verbose_stream() << "bounded ranges\n";
+        verbose_stream() << a_l_b << ": v" << x << " y: " << y << " := " << y_val << "\n";
+        verbose_stream() << p << " " << lo_p << " " << hi_p << "\n";
+        verbose_stream() << q << " " << lo_q << " " << hi_q << "\n";
+        verbose_stream() << r << " " << lo_r << " " << hi_r << "\n";
+        verbose_stream() << x_min << " " << x_max << "\n";
+
+        
+        // for every value of x, p, q are bewteen lo_p, hi_p, lo_q, hi_q
+        // if a_l_b is non-strict, it is false under all assignments to x, so r > 0
+        // if a_l_b is strict, we bail also if r = 1
+        if (lo_r == 0)
+            return false;
+        if (a_l_b.is_strict() && lo_r == 1)
+            return false;
+        
+        // then compute range around y that is admissible based on x_lo, x_hi
+        rational max_y, min_y = m.max_value();
+
+        if (q.degree(x) == 1) {
+            // q = x*y + b <= 2^N - 1
+            // y <= (2^N - 1 - b) / x_max
+            max_y = floor((m.max_value() - b_val)/x_max);
+            verbose_stream() << b << " " << b_val << " max-y " << max_y << "\n";
+        }
+
+        return false;
+    }
+
     
     /*
      * TODO
diff --git a/src/math/polysat/saturation.h b/src/math/polysat/saturation.h
index 4ce268021..6f2ad19f3 100644
--- a/src/math/polysat/saturation.h
+++ b/src/math/polysat/saturation.h
@@ -63,7 +63,10 @@ namespace polysat {
         bool try_tangent(pvar v, conflict& core, inequality const& c);
         bool try_add_overflow_bound(pvar x, conflict& core, inequality const& axb_l_y);
         bool try_add_mul_bound(pvar x, conflict& core, inequality const& axb_l_y);
+        bool try_add_mul_bound2(pvar x, conflict& core, inequality const& axb_l_y);
 
+        bool get_bound(pvar x, rational const& bound_x, pdd const& p, rational& bound_p);
+        
         // c := lhs ~ v
         //  where ~ is < or <=
         bool is_l_v(pvar v, inequality const& c);
diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp
index e0568328a..bd0d64e75 100644
--- a/src/math/polysat/viable.cpp
+++ b/src/math/polysat/viable.cpp
@@ -663,6 +663,8 @@ namespace polysat {
         entry const* e = first;
         bool found = false;
         out_c.reset();
+        if (!e)
+            return false;
         do {
             found = false;
             do {
@@ -695,6 +697,8 @@ namespace polysat {
         entry const* e = first;
         bool found = false;
         out_c.reset();
+        if (!e)
+            return false;
         do {
             found = false;
             do {
@@ -727,6 +731,20 @@ namespace polysat {
         entry const* first = m_units[v];
         entry const* e = first;
         bool found = false;
+        if (!e)
+            return false;
+
+        auto covers_all = [&](rational const& lo1, rational const& hi1, rational const lo2, rational const& hi2) {
+            SASSERT(lo1 != hi1);
+            if (lo1 < hi1) {
+                return lo2 <= hi1 && lo1 <= hi2;
+            }
+            else {
+                // hi1 < lo1
+                return hi1 <= hi2 && hi2 <= lo2 && lo2 <= lo1;
+            }
+        };
+        
         do {
             found = false;
             do {
@@ -737,13 +755,16 @@ namespace polysat {
                 auto const& hi = e->interval.hi();
                 if (!lo.is_val() || !hi.is_val())
                     goto next;
-
+                if (out_c.contains(e->src))
+                    goto next;
                 if (out_c.empty()) {
                     out_c.push_back(e->src);
                     out_lo = lo.val();
                     out_hi = hi.val();
                     found = true;
                 }
+                else if (covers_all(out_lo, out_hi, lo.val(), hi.val()))
+                    return false;
                 // [lo, hi0, hi[
                 // [lo, hi0, 0, hi[
                 else if (lo.val() <= out_hi && (out_hi < hi.val() || hi.val() < lo.val())) {