diff --git a/src/math/polysat/saturation.cpp b/src/math/polysat/saturation.cpp
index f8da38a4b..10f8e5182 100644
--- a/src/math/polysat/saturation.cpp
+++ b/src/math/polysat/saturation.cpp
@@ -1140,7 +1140,8 @@ namespace polysat {
      */
     bool saturation::try_add_overflow_bound(pvar x, conflict& core, inequality const& axb_l_y) {
         set_rule("[x] x >= x + y & x <= n => y = 0 or y >= 2^N - n");
-        signed_constraint y_eq_0, x_ge_bound;
+        signed_constraint y_eq_0;
+        vector<signed_constraint> x_ge_bound;
         auto& m = s.var2pdd(x);
         pdd y = m.zero();
         if (!is_add_overflow(x, axb_l_y, y))
@@ -1154,7 +1155,8 @@ namespace polysat {
         m_lemma.reset();
         if (!axb_l_y.is_strict())
             m_lemma.insert_eval(y_eq_0);        
-        m_lemma.insert_eval(~x_ge_bound);
+        for (auto c : x_ge_bound)
+            m_lemma.insert_eval(~c);
         return propagate(x, core, axb_l_y, s.uge(y, m.two_to_N() - bound));
     }
 
@@ -1186,11 +1188,11 @@ namespace polysat {
         return false;
     }
 
-    bool saturation::has_upper_bound(pvar x, conflict& core, rational& bound, signed_constraint& x_le_bound) {
+    bool saturation::has_upper_bound(pvar x, conflict& core, rational& bound, vector<signed_constraint>& x_le_bound) {
         return s.m_viable.has_upper_bound(x, bound, x_le_bound);
     }
 
-    bool saturation::has_lower_bound(pvar x, conflict& core, rational& bound, signed_constraint& x_ge_bound) {
+    bool saturation::has_lower_bound(pvar x, conflict& core, rational& bound, vector<signed_constraint>& x_ge_bound) {
         return s.m_viable.has_lower_bound(x, bound, x_ge_bound);
     }
 
@@ -1250,8 +1252,16 @@ namespace polysat {
      * -1195: v85 + 1 > 2^128+1
      * v25 := 353
      * v81 := -1
+     *
+     * -1*v85*v25 + v81 < v25    
+     * a -1*v25 := -315  b v81 := -1   y v25 := 315
+     * & v25 <= 315
+     * & -v81 <= 1 
+     * 
+     * The example illustrates that fixing y_val produces a weaker bound.
+     * The result should be a forbidden interval around v25 based on bounds for
+     * v85 and v81.
      * 
-     * Note that lower bound for v85 should use forbidden intervals which subsumes pattern matching against -1195.
      */
     
     bool saturation::try_add_mul_bound(pvar x, conflict& core, inequality const& a_l_b) {
@@ -1260,17 +1270,18 @@ namespace polysat {
         pdd const X = s.var(x);
         pdd a = s.var(x);
         pdd b = a, c = a, y = a;
-        rational b_val, c_val, y_val, x_bound;
-        signed_constraint x_ge_bound, x_le_bound, b_bound, ax_bound;
-        if (is_AxB_l_Y(x, a_l_b, a, b, y) && s.try_eval(y, y_val) && s.try_eval(b, b_val) && !y_val.is_zero() && !a.is_val()) {
-            verbose_stream() << a_l_b << "   a " << a << "   b " << b << " := " << dd::val_pp(m, b_val, false) << "   y " << y << "\n";
+        rational a_val, b_val, c_val, y_val, x_bound;
+        vector<signed_constraint> x_le_bound, x_ge_bound;
+        signed_constraint b_bound;
+        if (is_AxB_l_Y(x, a_l_b, a, b, y) && !a.is_val() && s.try_eval(y, y_val) && s.try_eval(b, b_val) && s.try_eval(a, a_val) && !y_val.is_zero()) {
+            IF_VERBOSE(2, verbose_stream() << "v" << x << ": " << a_l_b << "   a " << a << " := " << dd::val_pp(m, a_val, false) << "  b " << b << " := " << dd::val_pp(m, b_val, false) << "   y " << y << " := " << dd::val_pp(m, y_val, false) << "\n");
             SASSERT(!a.is_zero());
 
             // define c := -b
             c = -b;
             VERIFY(s.try_eval(c, c_val));
 
-            if (has_upper_bound(x, core, x_bound, x_le_bound) && x_le_bound != a_l_b.as_signed_constraint()) {
+            if (has_upper_bound(x, core, x_bound, x_le_bound) && !x_le_bound.contains(a_l_b.as_signed_constraint())) {
                 // ax - c <= y
                 // ==>  ax <= y + c                              if int(y) + int(c) <= 2^N, y <= int(y), c <= int(c)
                 // ==>  a not in [-floor(-int(y+c) / int(x), 0[
@@ -1278,18 +1289,23 @@ namespace polysat {
                 if (c_val + y_val <= m.max_value()) {
                     auto bound = floor((m.two_to_N() - y_val - c_val) / x_bound);
                     m_lemma.reset();
-                    m_lemma.insert_eval(~x_le_bound);           // x <= x_bound
+                    for (auto c : x_le_bound)
+                        m_lemma.insert_eval(~c);               // x <= x_bound
                     m_lemma.insert_eval(~s.ule(c, c_val));      // c <= c_val
                     m_lemma.insert_eval(~s.ule(y, y_val));      // y <= y_val
                     auto conclusion = s.uge(-a, bound);         // ==>  -a >= bound
-                    verbose_stream() << "XX: " << conclusion << "\n";
-                    if (propagate(x, core, a_l_b, conclusion)) {
+                    IF_VERBOSE(2, 
+                               verbose_stream() << core << "\n";
+                               verbose_stream() << "& " << X << " <= " << dd::val_pp(m, x_bound, false) << " := " << x_le_bound << "\n";
+                               verbose_stream() << "& " << s.ule(c, c_val) << "\n";
+                               verbose_stream() << "& " << s.ule(y, y_val) << "\n";
+                               verbose_stream() << "==> " << -a << " >= " << dd::val_pp(m, bound, false) << "\n");
+                    if (propagate(x, core, a_l_b, conclusion)) 
                         return true;
-                    }
                 }
                 // verbose_stream() << "TODO bound 1 " << a_l_b << " " << x_ge_bound << " " << b << " " << b_val << " " << y << "\n";
             }
-            if (has_lower_bound(x, core, x_bound, x_le_bound) && x_le_bound != a_l_b.as_signed_constraint()) {
+            if (has_lower_bound(x, core, x_bound, x_le_bound) && !x_le_bound.contains(a_l_b.as_signed_constraint())) {
                 
                 // verbose_stream() << "TODO bound 2 " << a_l_b << " " << x_le_bound << "\n";
             }
diff --git a/src/math/polysat/saturation.h b/src/math/polysat/saturation.h
index 2fe2119f7..86c2f0190 100644
--- a/src/math/polysat/saturation.h
+++ b/src/math/polysat/saturation.h
@@ -120,8 +120,9 @@ namespace polysat {
         // i := x + y >= x or x + y > x
         bool is_add_overflow(pvar x, inequality const& i, pdd& y);
 
-        bool has_upper_bound(pvar x, conflict& core, rational& bound, signed_constraint& x_ge_bound);
-        bool has_lower_bound(pvar x, conflict& core, rational& bound, signed_constraint& x_le_bound);
+        bool has_upper_bound(pvar x, conflict& core, rational& bound, vector<signed_constraint>& x_ge_bound);
+
+        bool has_lower_bound(pvar x, conflict& core, rational& bound, vector<signed_constraint>& x_le_bound);
 
         // determine min/max parity of polynomial
         unsigned min_parity(pdd const& p);
diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp
index f4d960352..09367b088 100644
--- a/src/math/polysat/viable.cpp
+++ b/src/math/polysat/viable.cpp
@@ -643,46 +643,69 @@ namespace polysat {
         return query<query_t::max_viable>(v, hi);
     }
 
-    // TBD: generalize to multiple intervals?
-    // Multiple intervals could be used to constrain upper/lower bounds, not just the last one.
-    bool viable::has_upper_bound(pvar v, rational& out_hi, signed_constraint& out_c) {
+    bool viable::has_upper_bound(pvar v, rational& out_hi, vector<signed_constraint>& out_c) {
         entry const* first = m_units[v];
         entry const* e = first;
+        bool found = false;
+        out_c.reset();
         do {
-            if (!e->refined) {
-                auto const& lo = e->interval.lo();
-                auto const& hi = e->interval.hi();
-                if (lo.is_val() && hi.is_val() && lo.val() > hi.val()) {
-                    out_c = e->src;
-                    out_hi = lo.val() - 1;
-                    // verbose_stream() << "Upper bound v" << v << " " << out_hi << " " << out_c << "\n";
-                    return true;
+            found = false;
+            do {
+                if (!e->refined) {
+                    auto const& lo = e->interval.lo();
+                    auto const& hi = e->interval.hi();
+                    if (lo.is_val() && hi.is_val()) {
+                        if (out_c.empty() && lo.val() > hi.val()) {
+                            out_c.push_back(e->src);
+                            out_hi = lo.val() - 1;
+                            found = true;
+                        }
+                        else if (!out_c.empty() && lo.val() <= out_hi && out_hi < hi.val()) {
+                            out_c.push_back(e->src);
+                            out_hi = lo.val() - 1;
+                            found = true;                        
+                        }
+                    }
                 }
+                e = e->next();
             }
-            e = e->next();
+            while (e != first);
         }
-        while (e != first);
-        return false;
+        while (found);
+        return !out_c.empty();
     }
 
-    bool viable::has_lower_bound(pvar v, rational& out_lo, signed_constraint& out_c) {
+    // TBD: generalize to multiple intervals similar to upper bound
+    bool viable::has_lower_bound(pvar v, rational& out_lo, vector<signed_constraint>& out_c) {
         entry const* first = m_units[v];
         entry const* e = first;
+        bool found = false;
+        out_c.reset();
         do {
-            if (!e->refined) {
-                auto const& lo = e->interval.lo();
-                auto const& hi = e->interval.hi();
-                if (lo.is_val() && hi.is_val() && (lo.val() == 0 || lo.val() > hi.val())) {
-                    out_c = e->src;
-                    out_lo = hi.val();
-                    // verbose_stream() << "Lower bound " << v << " " << out_lo << " " << out_c << "\n";
-                    return true;
+            found = false;
+            do {
+                if (!e->refined) {
+                    auto const& lo = e->interval.lo();
+                    auto const& hi = e->interval.hi();
+                    if (lo.is_val() && hi.is_val()) {
+                        if (out_c.empty() && (lo.val() == 0 || lo.val() > hi.val())) {
+                            out_c.push_back(e->src);
+                            out_lo = hi.val();
+                            found = true;
+                        }
+                        else if (!out_c.empty() && lo.val() <= out_lo && out_lo < hi.val()) {
+                            out_c.push_back(e->src);
+                            out_lo = hi.val();
+                            found = true;
+                        }
+                    }
                 }
+                e = e->next();
             }
-            e = e->next();
+            while (e != first);
         }
-        while (e != first);
-        return false;
+        while (found);
+        return !out_c.empty();
     }
 
     template <query_t mode>
diff --git a/src/math/polysat/viable.h b/src/math/polysat/viable.h
index a4f3f6b24..e677bb550 100644
--- a/src/math/polysat/viable.h
+++ b/src/math/polysat/viable.h
@@ -189,13 +189,13 @@ namespace polysat {
          * Query for an upper bound literal for v together with justification.
          * @return true if a non-trivial upper bound is found, return justifying constraint.
          */
-        bool has_upper_bound(pvar v, rational& out_hi, signed_constraint& out_c);
+        bool has_upper_bound(pvar v, rational& out_hi, vector<signed_constraint>& out_c);
 
         /**
          * Query for an lower bound literal for v together with justification.
          * @return true if a non-trivial lower bound is found, return justifying constraint.
          */
-        bool has_lower_bound(pvar v, rational& out_lo, signed_constraint& out_c);
+        bool has_lower_bound(pvar v, rational& out_lo, vector<signed_constraint>& out_c);
 
         /**
          * Find a next viable value for variable.