diff --git a/src/math/polysat/boolean.cpp b/src/math/polysat/boolean.cpp
index 2691c1214..3373cae07 100644
--- a/src/math/polysat/boolean.cpp
+++ b/src/math/polysat/boolean.cpp
@@ -83,9 +83,9 @@ namespace polysat {
         SASSERT(is_value_propagation(lit));
     }
 
-    void bool_var_manager::asserted(sat::literal lit, unsigned lvl, dependency dep) {
+    void bool_var_manager::assumption(sat::literal lit, unsigned lvl, dependency dep) {
         LOG("Asserted " << lit << " @ " << lvl);
-        assign(dep == null_dependency ? kind_t::decision : kind_t::assumption, lit, lvl, nullptr, dep);
+        assign(kind_t::assumption, lit, lvl, nullptr, dep);
         SASSERT(is_decision(lit) || is_assumption(lit));
     }
 
diff --git a/src/math/polysat/boolean.h b/src/math/polysat/boolean.h
index cb942da33..7b85dcd70 100644
--- a/src/math/polysat/boolean.h
+++ b/src/math/polysat/boolean.h
@@ -91,7 +91,7 @@ namespace polysat {
         void decide(sat::literal lit, unsigned lvl, clause& lemma);
         void decide(sat::literal lit, unsigned lvl);
         void eval(sat::literal lit, unsigned lvl);
-        void asserted(sat::literal lit, unsigned lvl, dependency dep);
+        void assumption(sat::literal lit, unsigned lvl, dependency dep);
         void unassign(sat::literal lit);
 
         std::ostream& display(std::ostream& out) const;
diff --git a/src/math/polysat/conflict.cpp b/src/math/polysat/conflict.cpp
index f7c4d26de..f740b39bf 100644
--- a/src/math/polysat/conflict.cpp
+++ b/src/math/polysat/conflict.cpp
@@ -50,8 +50,11 @@ namespace polysat {
         if (!m_vars.empty())
             out << " vars";
         for (auto v : m_vars)
-            if (is_pmarked(v))
-                out << " v" << v;
+            out << " v" << v;
+        if (!m_bail_vars.empty())
+            out << " bail vars";
+        for (auto v : m_bail_vars)
+            out << " v" << v;
         return out;
     }
 
@@ -61,6 +64,7 @@ namespace polysat {
         m_constraints.reset();
         m_literals.reset();
         m_vars.reset();
+        m_bail_vars.reset();
         m_conflict_var = null_var;
         m_bailout = false;
         SASSERT(empty());        
@@ -83,7 +87,7 @@ namespace polysat {
         else {
             SASSERT(c.is_currently_false(s));
             // TBD: fails with test_subst       SASSERT(c.bvalue(s) == l_true);
-            c->set_var_dependent();
+            insert_vars(c);
             insert(c);
         }
         SASSERT(!empty());
@@ -140,6 +144,25 @@ namespace polysat {
             m_constraints.push_back(c);
     }
 
+    void conflict::propagate(signed_constraint c) {
+        cm().ensure_bvar(c.get());
+        switch (c.bvalue(s)) {
+        case l_undef:
+            s.assign_eval(c.blit());
+            break;
+        case l_true:
+            break;
+        case l_false:            
+            break;
+        }
+        insert(c);
+    }
+
+    void conflict::insert_vars(signed_constraint c) {
+        for (pvar v : c->vars()) 
+            if (s.is_assigned(v)) 
+                m_vars.insert(v);  
+    }
 
     /**
      * Premises can either be justified by a Clause or by a value assignment.
@@ -210,7 +233,7 @@ namespace polysat {
         SASSERT(std::count(cl.begin(), cl.end(), ~lit) == 0);
         
         remove_literal(lit);        
-        unset_bmark(s.lit2cnstr(lit));         
+        unset_mark(s.lit2cnstr(lit));         
         for (sat::literal lit2 : cl)
             if (lit2 != lit)
                 insert(s.lit2cnstr(~lit2));
@@ -226,14 +249,17 @@ namespace polysat {
         SASSERT(contains_literal(lit));
         SASSERT(!contains_literal(~lit));
 
-        remove_literal(lit);
         signed_constraint c = s.lit2cnstr(lit);
-        unset_mark(c);
-        for (pvar v : c->vars()) {
-            if (s.is_assigned(v) && s.get_level(v) <= lvl) {
-                m_vars.insert(v);  
-                inc_pref(v);
-            }
+        bool has_decision = false;
+        for (pvar v : c->vars()) 
+            if (s.is_assigned(v) && s.m_justification[v].is_decision()) 
+                m_bail_vars.insert(v), has_decision = true;
+
+        if (!has_decision) {
+            remove(c);
+            for (pvar v : c->vars()) 
+                if (s.is_assigned(v)) 
+                    m_vars.insert(v);            
         }
     }
 
@@ -268,8 +294,6 @@ namespace polysat {
             lemma.push(~c);
 
         for (unsigned v : m_vars) {
-            if (!is_pmarked(v))
-                continue;
             auto eq = s.eq(s.var(v), s.get_value(v));
             cm().ensure_bvar(eq.get());
             if (eq.bvalue(s) == l_undef) 
@@ -314,42 +338,24 @@ namespace polysat {
     bool conflict::resolve_value(pvar v) {
         // NOTE:
         // In the "standard" case where "v = val" is on the stack:
-        //      - core contains both false and true constraints 
-        //        (originally only false ones, but additional true ones may come from saturation)
-
         // forbidden interval projection is performed at top level
-        SASSERT(v != conflict_var());
 
-        if (is_bailout()) {
-            if (!s.m_justification[v].is_decision())
-                m_vars.remove(v);
-            return false;
-        }
+        SASSERT(v != conflict_var());
 
         auto const& j = s.m_justification[v];
 
-        s.inc_activity(v); 
-
-#if 0
-        if (j.is_decision() && m_vars.contains(v)) {
-            set_bailout();
+        if (j.is_decision() && m_bail_vars.contains(v))
             return false;
-        }
-#endif
-   
+        
+        s.inc_activity(v);    
         m_vars.remove(v);
 
-        if (j.is_propagation()) {
-            for (auto const& c : s.m_viable.get_constraints(v)) {
-                if (!c->has_bvar()) {
-                    // it is a side-condition that was not propagated already.
-                    // it is true in the current variable assignment.
-                    cm().ensure_bvar(c.get());
-                    s.assign_eval(c.blit());
-                }
-                insert(c);
-            }        
-        }
+        if (is_bailout())
+            goto bailout;
+        
+        if (j.is_propagation()) 
+            for (auto const& c : s.m_viable.get_constraints(v)) 
+                propagate(c);
 
         LOG("try-explain v" << v);
         if (try_explain(v))
@@ -366,6 +372,7 @@ namespace polysat {
         }
         LOG("bailout v" << v);
         set_bailout();
+    bailout:
         if (s.is_assigned(v) && j.is_decision())
             m_vars.insert(v);
         return false;
@@ -406,19 +413,12 @@ namespace polysat {
         c->set_mark();
         if (c->has_bvar())
             set_bmark(c->bvar());    
-        if (c->is_var_dependent()) {
-            for (auto v : c->vars()) {
-                if (s.is_assigned(v))
-                    m_vars.insert(v);
-                inc_pref(v);
-            }
-        }
     }
 
     /**
      * unset marking on the constraint, but keep variable dependencies.
      */
-    void conflict::unset_bmark(signed_constraint c) {
+    void conflict::unset_mark(signed_constraint c) {
         if (!c->is_marked())
             return;
         c->unset_mark();
@@ -426,32 +426,6 @@ namespace polysat {
             unset_bmark(c->bvar());
     }
 
-    void conflict::unset_mark(signed_constraint c) {
-        if (!c->is_marked())
-            return;
-        unset_bmark(c);
-        if (!c->is_var_dependent()) 
-            return;
-        c->unset_var_dependent();
-        for (auto v : c->vars())
-            dec_pref(v);        
-    }
-
-    void conflict::inc_pref(pvar v) {
-        if (v >= m_pvar2count.size())
-            m_pvar2count.resize(v + 1);        
-        m_pvar2count[v]++;
-    }
-
-    void conflict::dec_pref(pvar v) {
-        SASSERT(m_pvar2count[v] > 0);
-        m_pvar2count[v]--;
-    }
-
-    bool conflict::is_pmarked(pvar v) const {
-        return m_pvar2count.get(v, 0) > 0;
-    }
-
     void conflict::set_bmark(sat::bool_var b) {
         if (b >= m_bvar2mark.size())
             m_bvar2mark.resize(b + 1);
diff --git a/src/math/polysat/conflict.h b/src/math/polysat/conflict.h
index 38a7a6c17..7fc013686 100644
--- a/src/math/polysat/conflict.h
+++ b/src/math/polysat/conflict.h
@@ -10,6 +10,64 @@ Author:
     Nikolaj Bjorner (nbjorner) 2021-03-19
     Jakob Rath 2021-04-6
 
+Notes:
+
+    A conflict state is of the form <Vars, Constraints>
+    Where Vars are shorthand for the constraints v = value(v) for v in Vars and value(v) is the assignent.
+    
+    The conflict state is unsatisfiable under background clauses F.
+    Dually, the negation is a consequence of F.
+
+    Conflict resolution resolves an assignment in the search stack against the conflict state.
+
+    Assignments are of the form:
+  
+    lit <- D => lit              - lit is propagated by the clause D => lit
+    lit <- ?                     - lit is a decision literal.
+    lit <- asserted              - lit is asserted
+    lit <- Vars                  - lit is propagated from variable evaluation.
+
+    v = value <- D               - v is assigned value by constraints D
+    v = value <- ?               - v is a decision literal.
+
+    - All literals should be assigned in the stack prior to their use.
+    
+    l <- D => l,       < Vars, { l } u C >   ===>   < Vars, C u D >
+    l <- ?,            < Vars, { l } u C >   ===>   ~l <- (C & Vars = value(Vars) => ~l)
+    l <- asserted,     < Vars, { l } u C >   ===>   < Vars, { l } u C >
+    l <- Vars',        < Vars, { l } u C >   ===>   < Vars u Vars', C >          if all Vars' are propagated
+    l <- Vars',        < Vars, { l } u C >   ===>   Mark < Vars, { l } u C > as bailout
+
+    v = value <- D,  < Vars u { v }, C >     ===>   < Vars, D u C >
+    v = value <- ?,  < Vars u { v }, C >     ===>   v != value <- (C & Vars = value(Vars) => v != value)
+    
+    
+Example derivations:
+
+Trail:       z <= y  <- asserted 
+             xz > xy <- asserted 
+             x = a   <- ? 
+             y = b   <- ? 
+             z = c   <- ?   
+Conflict:    < {x, y, z}, xz > xy > when ~O(a,b) and c <= b
+Append       x <= a <- { x }
+Append       y <= b <- { y }
+Conflict:    < {}, y >= z, xz > xy, x <= a, y <= b >
+Based on:    z <= y & x <= a & y <= b => xz <= xy
+Resolve:     y <= b <- { y }, y is a decision variable.
+Bailout:     lemma ~(y >= z & xz > xy & x <= a & y <= b) at decision level of lemma
+
+With overflow predicate:
+Append       ~O(x, y) <- { x, y }                  
+Conflict:    < {}, y >= z, xz > xy, ~O(x,y) >
+Based on     z <= y & ~O(x,y) => xz <= xy
+Resolve:     ~O(x, y) <- { x, y } both x, y are decision variables
+Lemma:       y < z or xz <= xy or O(x,y)
+
+
+
+ 
+    
 --*/
 #pragma once
 #include "math/polysat/constraint.h"
@@ -30,15 +88,12 @@ namespace polysat {
         signed_constraints m_constraints;   // new constraints used as premises
         indexed_uint_set m_literals;        // set of boolean literals in the conflict
         uint_set m_vars;                    // variable assignments used as premises
+        uint_set m_bail_vars;
 
         // If this is not null_var, the conflict was due to empty viable set for this variable.
         // Can be treated like "v = x" for any value x.
         pvar m_conflict_var = null_var;
 
-        unsigned_vector m_pvar2count;             // reference count of variables
-        void inc_pref(pvar v);
-        void dec_pref(pvar v);
-
         bool_vector m_bvar2mark;                  // mark of Boolean variables
         void set_bmark(sat::bool_var b);
         void unset_bmark(sat::bool_var b);
@@ -77,7 +132,7 @@ namespace polysat {
 
         void reset();
 
-        bool is_pmarked(pvar v) const;
+        bool contains_pvar(pvar v) const { return m_vars.contains(v) || m_bail_vars.contains(v); }
         bool is_bmarked(sat::bool_var b) const;
 
         /** conflict because the constraint c is false under current variable assignment */
@@ -87,7 +142,9 @@ namespace polysat {
         /** all literals in clause are false */
         void set(clause const& cl);
 
+        void propagate(signed_constraint c);
         void insert(signed_constraint c);
+        void insert_vars(signed_constraint c);
         void insert(signed_constraint c, vector<signed_constraint> const& premises);
         void remove(signed_constraint c);
         void replace(signed_constraint c_old, signed_constraint c_new, vector<signed_constraint> const& c_new_premises);
diff --git a/src/math/polysat/constraint.h b/src/math/polysat/constraint.h
index e1f385208..db131acba 100644
--- a/src/math/polysat/constraint.h
+++ b/src/math/polysat/constraint.h
@@ -146,7 +146,6 @@ namespace polysat {
         unsigned_vector     m_vars;
         lbool               m_external_sign = l_undef;
         bool                m_is_marked = false;
-        bool                m_is_var_dependent = false;
         bool                m_is_active = false;
         /** The boolean variable associated to this constraint, if any.
          *  If this is not null_bool_var, then the constraint corresponds to a literal on the assignment stack.
@@ -207,9 +206,6 @@ namespace polysat {
         void unset_mark() { m_is_marked = false; }
         bool is_marked() const { return m_is_marked; }
 
-        void set_var_dependent() { m_is_var_dependent = true; }
-        void unset_var_dependent() { m_is_var_dependent = false; }
-        bool is_var_dependent() const { return m_is_var_dependent;  }
 
         bool is_active() const { return m_is_active; }
         void set_active(bool f) { m_is_active = f; }
diff --git a/src/math/polysat/explain.cpp b/src/math/polysat/explain.cpp
index cbc7d4040..961ee6839 100644
--- a/src/math/polysat/explain.cpp
+++ b/src/math/polysat/explain.cpp
@@ -165,10 +165,9 @@ namespace polysat {
             auto c2 = s.ule(a, b);
             if (!c.is_positive())
                 c2 = ~c2;
-            LOG("try-reduce is false " << c2.is_currently_false(s));
             if (!c2.is_currently_false(s))
                 continue;
-            if (c2.bvalue(s) == l_false)
+            if (c2.is_always_false() || c2.bvalue(s) == l_false)
                 return false;
             if (!c2->has_bvar() || l_undef == c2.bvalue(s)) {
                 vector<signed_constraint> premises;
diff --git a/src/math/polysat/forbidden_intervals.cpp b/src/math/polysat/forbidden_intervals.cpp
index 63b97ed5b..0007abf08 100644
--- a/src/math/polysat/forbidden_intervals.cpp
+++ b/src/math/polysat/forbidden_intervals.cpp
@@ -73,11 +73,9 @@ namespace polysat {
         // a*v <= 0, a odd
         if (match_zero(c, a1, b1, e1, a2, b2, e2, fi))
             return true;
-
         // a*v + b > 0, a odd
         if (match_non_zero_linear(c, a1, b1, e1, a2, b2, e2, fi))
             return true;
-
         if (match_linear1(c, a1, b1, e1, a2, b2, e2, fi))
             return true;
         if (match_linear2(c, a1, b1, e1, a2, b2, e2, fi))
@@ -349,7 +347,7 @@ namespace polysat {
         signed_constraint const& c,
         rational const & a1, pdd const& b1, pdd const& e1,
         fi_record& fi) {
-        if (a1.is_odd() && b1.is_zero() && !c.is_positive()) {
+        if (a1.is_odd() && b1.is_zero() && c.is_negative()) {
             auto& m = e1.manager();
             rational lo_val(0);
             auto lo = m.zero();
@@ -361,7 +359,7 @@ namespace polysat {
                 fi.side_cond.push_back(s.eq(b1, e1));
             return true;
         }
-        if (a1.is_odd() && b1.is_val() && !c.is_positive()) {
+        if (a1.is_odd() && b1.is_val() && c.is_negative()) {
             auto& m = e1.manager();
             rational const& mod_value = m.two_to_N();
             rational a1_inv;
@@ -389,7 +387,7 @@ namespace polysat {
         signed_constraint const& c,
         rational const & a2, pdd const& b2, pdd const& e2,
         fi_record& fi) {
-        if (a2.is_one() && b2.is_val() && !c.is_positive()) {
+        if (a2.is_one() && b2.is_val() && c.is_negative()) {
             auto& m = e2.manager();
             rational const& mod_value = m.two_to_N();
             rational lo_val(mod(-b2.val() - 1, mod_value));
diff --git a/src/math/polysat/saturation.cpp b/src/math/polysat/saturation.cpp
index e68d87a75..a0e76297d 100644
--- a/src/math/polysat/saturation.cpp
+++ b/src/math/polysat/saturation.cpp
@@ -88,7 +88,7 @@ namespace polysat {
         for (auto d : m_new_constraints)
             core.insert(d);
         if (c.bvalue(s) != l_false)
-            c->set_var_dependent();
+            core.insert_vars(c);
         core.insert(~c);       
         LOG("Core " << core);
         return true;        
@@ -423,6 +423,22 @@ namespace polysat {
             return false;
         if (c.lhs.is_val() || c.rhs.is_val())
             return false;
+
+        pdd q_l(c.lhs), e_l(c.lhs), q_r(c.rhs), e_r(c.rhs);
+        bool is_linear = true;
+        is_linear &= c.lhs.degree(v) <= 1;
+        is_linear &= c.rhs.degree(v) <= 1;
+        if (c.lhs.degree(v) == 1) {
+            c.lhs.factor(v, 1, q_l, e_l); 
+            is_linear &= q_l.is_val();
+        }
+        if (c.rhs.degree(v) == 1) {
+            c.rhs.factor(v, 1, q_r, e_r);
+            is_linear &= q_r.is_val();
+        }
+        if (is_linear)
+            return false;
+
         if (!c.as_signed_constraint().is_currently_false(s))
             return false;
         rational l_val, r_val;
diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp
index f99fa41cd..25f094e33 100644
--- a/src/math/polysat/solver.cpp
+++ b/src/math/polysat/solver.cpp
@@ -190,14 +190,14 @@ namespace polysat {
         switch (m_bvars.value(lit)) {
         case l_false:
             set_conflict(c);
+            SASSERT(dep == null_dependency && "track dependencies is TODO");
             break;
         case l_true:
             // constraint c is already asserted
             SASSERT(m_bvars.level(lit) <= m_level);
-            // TODO: track additional dep?
             break;
         case l_undef:
-            m_bvars.asserted(lit, m_level, dep);
+            m_bvars.assumption(lit, m_level, dep);
             m_trail.push_back(trail_instr_t::assign_bool_i);
             m_search.push_boolean(lit);
             if (c.is_currently_false(*this))
@@ -535,7 +535,7 @@ namespace polysat {
                 // Resolve over variable assignment
                 pvar v = item.var();
                 LOG(m_justification[v]);
-                if (!m_conflict.is_pmarked(v) && !m_conflict.is_bailout()) {
+                if (!m_conflict.contains_pvar(v) && !m_conflict.is_bailout()) {
                     m_search.pop_assignment();
                     continue;
                 }
diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp
index 08ef1b9c5..b6ca77b3c 100644
--- a/src/math/polysat/viable.cpp
+++ b/src/math/polysat/viable.cpp
@@ -46,9 +46,11 @@ namespace polysat {
 
     void viable::pop_viable() {
         auto& [v, k, e] = m_trail.back();
+        SASSERT(well_formed(m_units[v]));
         switch (k) {
         case entry_kind::unit_e: 
             e->remove_from(m_units[v], e); 
+            SASSERT(well_formed(m_units[v]));
             break;
         case entry_kind::equal_e: 
             e->remove_from(m_equal_lin[v], e); 
@@ -67,13 +69,15 @@ namespace polysat {
         SASSERT(e->prev() != e || e->next() == e);
         SASSERT(k == entry_kind::unit_e);
         (void)k;
+        SASSERT(well_formed(m_units[v]));
         if (e->prev() != e) {
             e->prev()->insert_after(e);
-            if (e->interval.lo_val() < e->next()->interval.lo_val())
+            if (e->interval.lo_val() < m_units[v]->interval.lo_val())
                 m_units[v] = e;
         }
         else 
-            m_units[v] = e;        
+            m_units[v] = e;  
+        SASSERT(well_formed(m_units[v]));
         m_trail.pop_back();
     }
 
@@ -102,6 +106,7 @@ namespace polysat {
     }
 
     void viable::insert(entry* e, pvar v, ptr_vector<entry>& entries, entry_kind k) {
+        SASSERT(well_formed(m_units[v]));
         m_trail.push_back({ v, k, e });
         s.m_trail.push_back(trail_instr_t::viable_add_i);
         e->init(e);
@@ -109,6 +114,7 @@ namespace polysat {
             entries[v] = e;
         else
             e->insert_after(entries[v]);
+        SASSERT(well_formed(m_units[v]));
     }
 
     bool viable::intersect(pvar v, entry* ne) {
@@ -122,6 +128,7 @@ namespace polysat {
             m_alloc.push_back(ne);
             return false;
         }
+        
 
         auto create_entry = [&]() {
             m_trail.push_back({ v, entry_kind::unit_e, ne });
@@ -136,6 +143,13 @@ namespace polysat {
             e->remove_from(m_units[v], e);
         };
 
+        if (ne->interval.is_full()) {
+            while (m_units[v]) 
+                remove_entry(m_units[v]);
+            m_units[v] = create_entry();
+            return true;
+        }
+
         if (!e) 
             m_units[v] = create_entry();
         else {
@@ -257,16 +271,21 @@ namespace polysat {
                     lo = val - lambda_l; 
                     increase_hi(hi);
                 }
-                LOG("forbidden interval " << e->coeff << " * " << e->interval << " [" << lo << ", " << hi << "[");
+                LOG("forbidden interval v" << v << " " << val << " " << e->coeff << " * " << e->interval << " [" << lo << ", " << hi << "[");
                 SASSERT(hi <= mod_value);
-                if (hi == mod_value) hi = 0;
+                bool full = (lo == 0 && hi == mod_value);                    
+                if (hi == mod_value) 
+                    hi = 0;
                 pdd lop = s.var2pdd(v).mk_val(lo);
                 pdd hip = s.var2pdd(v).mk_val(hi);
                 entry* ne = alloc_entry();
                 ne->src = e->src;
                 ne->side_cond = e->side_cond;
                 ne->coeff = 1;
-                ne->interval = eval_interval::proper(lop, lo, hip, hi);
+                if (full)
+                    ne->interval = eval_interval::full();
+                else 
+                    ne->interval = eval_interval::proper(lop, lo, hip, hi);
                 intersect(v, ne);
                 return false;
             }
@@ -371,6 +390,8 @@ namespace polysat {
         entry* first = e;
         entry* last = e->prev();
 
+        if (e->interval.is_full())
+            return false;
         // quick check: last interval doesn't wrap around, so hi_val
         // has not been covered
         if (last->interval.lo_val() < last->interval.hi_val()) 
@@ -413,6 +434,7 @@ namespace polysat {
         return refine_viable(v, val);
     }
 
+
     rational viable::min_viable(pvar v) { 
         refined:
         rational lo(0);
@@ -542,18 +564,15 @@ namespace polysat {
                 auto lhs = hi - next_lo;
                 auto rhs = next_hi - next_lo;
                 signed_constraint c = s.m_constraints.ult(lhs, rhs);
-                core.insert(c);
+                core.propagate(c);
             }
-            for (auto sc : e->side_cond)
-                core.insert(sc);
-            e->src->set_var_dependent(); 
+            for (auto sc : e->side_cond) 
+                core.propagate(sc);
             core.insert(e->src);
             e = n;
         }             
         while (e != first);
 
-        // core.set_bailout();
-        // TODO - review this; c is true under current assignment?
         for (auto c : core) {
             if (c.bvalue(s) == l_false) {
                 core.reset();
@@ -606,7 +625,7 @@ namespace polysat {
 
     std::ostream& viable::display(std::ostream& out) const {
         for (pvar v = 0; v < m_units.size(); ++v)
-            display(out << "v" << v << ": ", v);
+            display(out << "v" << v << ": ", v) << "\n";
         return out;
     }
 
diff --git a/src/test/polysat.cpp b/src/test/polysat.cpp
index c9d737a60..8832833e0 100644
--- a/src/test/polysat.cpp
+++ b/src/test/polysat.cpp
@@ -117,14 +117,14 @@ namespace polysat {
     };
 
 
-class test_polysat {
-public:
+    class test_polysat {
+    public:
 
-    /**
-     * Testing the solver's internal state.
-     */
+        /**
+         * Testing the solver's internal state.
+         */
 
-    /// Creates two separate conflicts (from narrowing) before solving loop is started.
+        /// Creates two separate conflicts (from narrowing) before solving loop is started.
     static void test_add_conflicts() {
         scoped_solver s(__func__);
         auto a = s.var(s.add_var(3));
@@ -137,7 +137,7 @@ public:
         s.expect_unsat();
     }
 
-    /// Has constraints which must be inserted into other watchlist to discover UNSAT
+        /// Has constraints which must be inserted into other watchlist to discover UNSAT
     static void test_wlist() {
         scoped_solver s(__func__);
         auto a = s.var(s.add_var(3));
@@ -153,7 +153,7 @@ public:
         s.expect_unsat();
     }
 
-    /// Has a constraint in cjust[a] where a does not occur.
+        /// Has a constraint in cjust[a] where a does not occur.
     static void test_cjust() {
         scoped_solver s(__func__);
         auto a = s.var(s.add_var(3));
@@ -168,718 +168,739 @@ public:
         s.check();
         s.expect_unsat();
     }
-
-    /**
-     * most basic linear equation solving.
-     * they should be solvable.
-     * they also illustrate some limitations of basic solver even if it solves them.
-     * Example
-     *   the value to a + 1 = 0 is fixed at 3, there should be no search.
-     */
-
-    static void test_l1() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(2));
-        s.add_eq(a + 1);
-        s.check();
-        s.expect_sat({{a, 3}});
-    }
     
-    static void test_l2() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(2));
-        auto b = s.var(s.add_var(2));
-        s.add_eq(2*a + b + 1);
-        s.add_eq(2*b + a);
-        s.check();
-        s.expect_sat({{a, 2}, {b, 3}});
-    }
-
-    static void test_l3() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(2));
-        auto b = s.var(s.add_var(2));
-        s.add_eq(3*b + a + 2);
-        s.check();
-        s.expect_sat();
-    }
-
-    static void test_l4() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(3));
-        s.add_eq(4*a + 2);
-        s.check();
-        s.expect_unsat();
-    }
-
-    static void test_l5() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(3));
-        auto b = s.var(s.add_var(3));
-        s.add_diseq(b);
-        s.add_eq(a + 2*b + 4);
-        s.add_eq(a + 4*b + 4);
-        s.check();
-        s.expect_sat({{a, 4}, {b, 4}});
-    }
-
-    
-    /**
-     * This one is unsat because a*a*(a*a - 1)
-     * is 0 for all values of a.
-     */
-    static void test_p1() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(2));
-        auto p = a*a*(a*a - 1) + 1;
-        s.add_eq(p);
-        s.check();
-        s.expect_unsat();
-    }
-
-    /**
-     * has solutions a = 2 and a = 3
-     */
-    static void test_p2() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(2));
-        auto p = a*(a-1) + 2;
-        s.add_eq(p);
-        s.check();
-        s.expect_sat();
-    }
-
-    /**
-     * unsat
-     * - learns 3*x + 1 == 0 by polynomial resolution
-     * - this forces x == 5, which means the first constraint is unsatisfiable by parity.
-     */
-    static void test_p3() {
-        scoped_solver s(__func__);
-        auto x = s.var(s.add_var(4));
-        auto y = s.var(s.add_var(4));
-        auto z = s.var(s.add_var(4));
-        s.add_eq(x*x*y + 3*y + 7);
-        s.add_eq(2*y + z + 8);
-        s.add_eq(3*x + 4*y*z + 2*z*z + 1);
-        s.check();
-        s.expect_unsat();
-    }
-
-
-    // Unique solution: u = 5
-    static void test_ineq_basic1() {
-        scoped_solver s(__func__);
-        auto u = s.var(s.add_var(4));
-        s.add_ule(u, 5);
-        s.add_ule(5, u);
-        s.check();
-        s.expect_sat({{u, 5}});
-    }
-
-    // Unsatisfiable
-    static void test_ineq_basic2() {
-        scoped_solver s(__func__);
-        auto u = s.var(s.add_var(4));
-        s.add_ult(u, 5);
-        s.add_ule(5, u);
-        s.check();
-        s.expect_unsat();
-    }
-
-    // Solutions with u = v = w
-    static void test_ineq_basic3() {
-        scoped_solver s(__func__);
-        auto u = s.var(s.add_var(4));
-        auto v = s.var(s.add_var(4));
-        auto w = s.var(s.add_var(4));
-        s.add_ule(u, v);
-        s.add_ule(v, w);
-        s.add_ule(w, u);
-        s.check();
-        s.expect_sat();
-        SASSERT_EQ(s.get_value(u.var()), s.get_value(v.var()));
-        SASSERT_EQ(s.get_value(u.var()), s.get_value(w.var()));
-    }
-
-    // Unsatisfiable
-    static void test_ineq_basic4() {
-        scoped_solver s(__func__);
-        auto u = s.var(s.add_var(4));
-        auto v = s.var(s.add_var(4));
-        auto w = s.var(s.add_var(4));
-        s.add_ule(u, v);
-        s.add_ult(v, w);
-        s.add_ule(w, u);
-        s.check();
-        s.expect_unsat();
-    }
-
-    // Satisfiable
-    // Without forbidden intervals, we just try values for u until it works
-    static void test_ineq_basic5() {
-        scoped_solver s(__func__);
-        auto u = s.var(s.add_var(4));
-        auto v = s.var(s.add_var(4));
-        s.add_ule(12, u + v);
-        s.add_ule(v, 2);
-        s.check();
-        s.expect_sat();  // e.g., u = 12, v = 0
-    }
-
-    // Like test_ineq_basic5 but the other forbidden interval will be the longest
-    static void test_ineq_basic6() {
-        scoped_solver s(__func__);
-        auto u = s.var(s.add_var(4));
-        auto v = s.var(s.add_var(4));
-        s.add_ule(14, u + v);
-        s.add_ule(v, 2);
-        s.check();
-        s.expect_sat();
-    }
-
-
-    /**
-     * Check unsat of:
-     * u = v*q + r
-     * r < u
-     * v*q > u
-     */
-    static void test_ineq1() {
-        scoped_solver s(__func__);
-        auto u = s.var(s.add_var(5));
-        auto v = s.var(s.add_var(5));
-        auto q = s.var(s.add_var(5));
-        auto r = s.var(s.add_var(5));
-        s.add_eq(u - (v*q) - r);
-        s.add_ult(r, u);
-        s.add_ult(u, v*q);
-        s.check();
-        s.expect_unsat();
-    }
-
-    /**
-     * Check unsat of:
-     * n*q1 = a - b
-     * n*q2 + r2 = c*a - c*b
-     * n > r2 > 0
-     */
-    static void test_ineq2() {
-        scoped_solver s(__func__);
-        auto n = s.var(s.add_var(5));
-        auto q1 = s.var(s.add_var(5));
-        auto a = s.var(s.add_var(5));
-        auto b = s.var(s.add_var(5));
-        auto c = s.var(s.add_var(5));
-        auto q2 = s.var(s.add_var(5));
-        auto r2 = s.var(s.add_var(5));
-        s.add_eq(n*q1 - a + b);
-        s.add_eq(n*q2 + r2 - c*a + c*b);
-        s.add_ult(r2, n);
-        s.add_diseq(r2);
-        s.check();
-        s.expect_unsat();
-    }
-
-
-    /**
-     * Monotonicity example from certora
-     * 
-     * We do overflow checks by doubling the base bitwidth here.
-     */
-    static void test_monot(unsigned base_bw = 5) {
-        scoped_solver s(std::string{__func__} + "(" + std::to_string(base_bw) + ")");
-
-        auto max_int_const = rational::power_of_two(base_bw) - 1;
-
-        unsigned bw = 2 * base_bw;
-        auto max_int = s.var(s.add_var(bw));
-        s.add_eq(max_int + (-max_int_const));
-
-        auto tb1 = s.var(s.add_var(bw));
-        s.add_ule(tb1, max_int);
-        auto tb2 = s.var(s.add_var(bw));
-        s.add_ule(tb2, max_int);
-        auto a = s.var(s.add_var(bw));
-        s.add_ule(a, max_int);
-        auto v = s.var(s.add_var(bw));
-        s.add_ule(v, max_int);
-        auto base1 = s.var(s.add_var(bw));
-        s.add_ule(base1, max_int);
-        auto base2 = s.var(s.add_var(bw));
-        s.add_ule(base2, max_int);
-        auto elastic1 = s.var(s.add_var(bw));
-        s.add_ule(elastic1, max_int);
-        auto elastic2 = s.var(s.add_var(bw));
-        s.add_ule(elastic2, max_int);
-        auto err = s.var(s.add_var(bw));
-        s.add_ule(err, max_int);
-
-        auto rem1 = s.var(s.add_var(bw));
-        auto quot2 = s.var(s.add_var(bw));
-        s.add_ule(quot2, max_int);
-        auto rem2 = s.var(s.add_var(bw));
-        auto rem3 = s.var(s.add_var(bw));
-        auto quot4 = s.var(s.add_var(bw));
-        s.add_ule(quot4, max_int);
-        auto rem4 = s.var(s.add_var(bw));
-
-        s.add_diseq(elastic1);
-
-        // division: tb1 = (v * base1) / elastic1;
-        s.add_eq((tb1 * elastic1) + rem1 - (v * base1));
-        s.add_ult(rem1, elastic1);
-        s.add_ule((tb1 * elastic1), max_int);
-
-        // division: quot2 = (a * base1) / elastic1
-        s.add_eq((quot2 * elastic1) + rem2 - (a * base1));
-        s.add_ult(rem2, elastic1);
-        s.add_ule((quot2 * elastic1), max_int);
-
-        s.add_eq(base1 + quot2 - base2);
-
-        s.add_eq(elastic1 + a - elastic2);
-
-        // division: tb2 = ((v * base2) / elastic2);
-        s.add_eq((tb2 * elastic2) + rem3 - (v * base2));
-        s.add_ult(rem3, elastic2);
-        s.add_ule((tb2 * elastic2), max_int);
-
-        // division: quot4 = v / elastic2;
-        s.add_eq((quot4 * elastic2) + rem4 - v);
-        s.add_ult(rem4, elastic2);
-        s.add_ule((quot4 * elastic2), max_int);
-
-        s.add_eq(quot4 + 1 - err);
-
-        s.push();
-        s.add_ult(tb1, tb2);
-        s.check();
-        s.expect_unsat();
-        s.pop();
-
-        s.push();
-        s.add_ult(tb2 + err, tb1);
-        s.check();
-        s.expect_unsat();
-        s.pop();
-    }
-
-    /*
-     * Mul-then-div in fixed point arithmetic is (roughly) neutral.
-     *
-     * I.e. we prove "(((a * b) / sf) * sf) / b" to be equal to a, up to some error margin.
-     * 
-     * sf is the scaling factor (we could leave this unconstrained, but non-zero, to make the benchmark a bit harder)
-     * em is the error margin
-     * 
-     * We do overflow checks by doubling the base bitwidth here.
-     */
-    static void test_fixed_point_arith_mul_div_inverse() {
-        scoped_solver s(__func__);
-
-        auto baseBw = 5;
-        auto max_int_const = 31; // (2^5 - 1) -- change this when you change baseBw
-
-        auto bw = 2 * baseBw;
-        auto max_int = s.var(s.add_var(bw));
-        s.add_eq(max_int - max_int_const);
-
-        // "input" variables
-        auto a = s.var(s.add_var(bw));
-        s.add_ule(a, max_int);
-        auto b = s.var(s.add_var(bw));
-        s.add_ule(b, max_int);
-        s.add_ult(0, b); // b > 0
-
-        // scaling factor (setting it, somewhat arbitrarily, to max_int/3)
-        auto sf = s.var(s.add_var(bw));
-        s.add_eq(sf - (max_int_const/3));
-
-        // (a * b) / sf = quot1 <=> quot1 * sf + rem1 - (a * b) = 0
-        auto quot1 = s.var(s.add_var(bw));
-        auto rem1 = s.var(s.add_var(bw));
-        s.add_eq((quot1 * sf) + rem1 - (a * b));
-        s.add_ult(rem1, sf);
-        s.add_ule(quot1 * sf, max_int);
-
-        // (((a * b) / sf) * sf) / b <=> quot2 * b + rem2 - (((a * b) / sf) * sf) = 0
-        auto quot2 = s.var(s.add_var(bw));
-        auto rem2 = s.var(s.add_var(bw));
-        s.add_eq((quot2 * b) + rem2 - (quot1 * sf));
-        s.add_ult(rem2, b);
-        s.add_ule(quot2 * b, max_int);
-
-        // sf / b = quot3 <=> quot3 * b + rem3 = sf
-        auto quot3 = s.var(s.add_var(bw));
-        auto rem3 = s.var(s.add_var(bw));
-        s.add_eq((quot3 * b) + rem3 - sf);
-        s.add_ult(rem3, b);
-        s.add_ule(quot3 * b, max_int);
-
-        // em = sf / b + 1
-        auto em = s.var(s.add_var(bw));
-        s.add_eq(quot3 + 1 - em);
-
-        // we prove quot3 <= a and quot3 + em >= a
-
-        s.push();
-        s.add_ult(a, quot3);
-        s.check();
-        s.expect_unsat();
-        s.pop();
-
-
-        // s.push();
-        // s.add_ult(quot3 + em, a);
-        // s.check();
-        // s.expect_unsat();
-        // s.pop();
-    }
-
-    /*
-     * Div-then-mul in fixed point arithmetic is (roughly) neutral.
-     *
-     * I.e. we prove "(b * ((a * sf) / b)) / sf" to be equal to a, up to some error margin.
-     * 
-     * sf is the scaling factor (we could leave this unconstrained, but non-zero, to make the benchmark a bit harder)
-     * em is the error margin
-     * 
-     * We do overflow checks by doubling the base bitwidth here.
-     */
-    static void test_fixed_point_arith_div_mul_inverse(unsigned base_bw = 5) {
-        scoped_solver s(__func__);
-
-        auto max_int_const = rational::power_of_two(base_bw) - 1;
-
-        auto bw = 2 * base_bw;
-        auto max_int = s.var(s.add_var(bw));
-        s.add_eq(max_int - max_int_const);
-
-        // "input" variables
-        auto a = s.var(s.add_var(bw));
-        s.add_ule(a, max_int);
-        auto b = s.var(s.add_var(bw));
-        s.add_ule(b, max_int);
-        s.add_ult(0, b); // b > 0
-
-        // scaling factor (setting it, somewhat arbitrarily, to max_int/3)
-        auto sf = s.var(s.add_var(bw));
-        s.add_eq(sf - floor(max_int_const/3));
-
-        // (a * sf) / b = quot1 <=> quot1 * b + rem1 - (a * sf) = 0
-        auto quot1 = s.var(s.add_var(bw));
-        auto rem1 = s.var(s.add_var(bw));
-        s.add_eq((quot1 * b) + rem1 - (a * sf));
-        s.add_ult(rem1, b);
-        s.add_ule(quot1 * b, max_int);
-
-        // (b * ((a * sf) / b)) / sf = quot2 <=> quot2 * sf + rem2 - (b * ((a * sf) / b)) = 0
-        auto quot2 = s.var(s.add_var(bw));
-        auto rem2 = s.var(s.add_var(bw));
-        s.add_eq((quot2 * sf) + rem2 - (b * quot1));
-        s.add_ult(rem2, sf);
-        s.add_ule(quot2 * sf, max_int);
-
-        // b / sf = quot3 <=> quot3 * sf + rem3 - b = 0
-        auto quot3 = s.var(s.add_var(bw));
-        auto rem3 = s.var(s.add_var(bw));
-        s.add_eq((quot3 * sf) + rem3 - b);
-        s.add_ult(rem3, sf);
-        s.add_ule(quot3 * sf, max_int);
-
-        // em = b / sf + 1
-        auto em = s.var(s.add_var(bw));
-        s.add_eq(quot3 + 1 - em);
-
-        // we prove quot3 <= a and quot3 + em >= a
-
-        s.push();
-        s.add_ult(quot3 + em, a);
-        s.check();
-        //        s.expect_unsat();
-        s.pop();
-
-        s.push();
-        s.add_ult(a, quot3);
-        s.check();
-//        s.expect_unsat();
-        s.pop();
-
-
-
-        //exit(0);
-    }
-
-    /** Monotonicity under bounds,
-     * puzzle extracted from https://github.com/NikolajBjorner/polysat/blob/main/puzzles/bv.smt2
-     *
-     * x, y, z \in [0..2^64[
-     * x, y, z < 2^32
-     * y <= x
-     * x*z < 2^32
-     * y*z >= 2^32
-     */
-    static void test_monot_bounds(unsigned base_bw = 32) {
-        scoped_solver s(std::string{__func__} + "(" + std::to_string(base_bw) + ")");
-        unsigned bw = 2 * base_bw;
-        auto y = s.var(s.add_var(bw));
-        auto z = s.var(s.add_var(bw));
-        auto x = s.var(s.add_var(bw));
-        auto bound = rational::power_of_two(base_bw);
-#if 1
-        s.add_ult(x, bound);
-        s.add_ult(y, bound);
-        // s.add_ult(z, bound);   // not required
-#else
-        s.add_ule(x, bound - 1);
-        s.add_ule(y, bound - 1);
-        // s.add_ule(z, bound - 1);   // not required
-#endif
-        unsigned a = 13;
-        s.add_ule(z, y);
-        s.add_ult(x*y, a);
-        s.add_ule(a, x*z);
-        s.check();
-        s.expect_unsat();
-    }
-
-    /** Monotonicity under bounds, simplified even more.
-     *
-     * x, y, z \in [0..2^64[
-     * x, y, z < 2^32
-     * z <= y
-     * y*x < z*x
-     */
-    static void test_monot_bounds_simple(unsigned base_bw = 32) {
-        scoped_solver s(__func__);
-        unsigned bw = 2 * base_bw;
-        /*
-        auto z = s.var(s.add_var(bw));
-        auto x = s.var(s.add_var(bw));
-        auto y = s.var(s.add_var(bw));
-        */
-        auto y = s.var(s.add_var(bw));
-        auto x = s.var(s.add_var(bw));
-        auto z = s.var(s.add_var(bw));
-        auto bound = rational::power_of_two(base_bw);
-        s.add_ult(x, bound);
-        s.add_ult(y, bound);
-        s.add_ult(z, bound);
-        s.add_ule(z, y);
-        s.add_ult(y*x, z*x);
-        s.check();
-        s.expect_unsat();
-    }
-
-    /*
-     * Transcribed from https://github.com/NikolajBjorner/polysat/blob/main/puzzles/bv.smt2
-     *
-     * We do overflow checks by doubling the base bitwidth here.
-     */
-    static void test_monot_bounds_full(unsigned base_bw = 5) {
-        scoped_solver s(__func__);
-
-        auto const max_int_const = rational::power_of_two(base_bw) - 1;
-
-        auto const bw = 2 * base_bw;
-        auto const max_int = s.var(s.add_var(bw));
-        s.add_eq(max_int - max_int_const);
-
-        auto const first = s.var(s.add_var(bw));
-        s.add_ule(first, max_int);
-        auto const second = s.var(s.add_var(bw));
-        s.add_ule(second, max_int);
-        auto const idx = s.var(s.add_var(bw));
-        s.add_ule(idx, max_int);
-        auto const q = s.var(s.add_var(bw));
-        s.add_ule(q, max_int);
-        auto const r = s.var(s.add_var(bw));
-        s.add_ule(r, max_int);
-
-        // q = max_int / idx <=> q * idx + r - max_int = 0
-        s.add_eq((q * idx) + r - max_int);
-        s.add_ult(r, idx);
-        s.add_ule(q * idx, max_int);
-
-        /*  last assertion:
-                (not
-                    (=> (bvugt second first) 
-                        (=> 
-                            (=> (not (= idx #x00000000)) 
-                                (bvule (bvsub second first) q)) 
-                            (bvumul_noovfl (bvsub second first) idx))))
-            transforming negated boolean skeleton:
-                (not (=> a (=> (or b c) d))) <=> (and a (not d) (or b c))
+        /**
+         * most basic linear equation solving.
+         * they should be solvable.
+         * they also illustrate some limitations of basic solver even if it solves them.
+         * Example
+         *   the value to a + 1 = 0 is fixed at 3, there should be no search.
          */
 
-        // (bvugt second first)
-        s.add_ult(first, second);
-        // (not (bvumul_noovfl (bvsub second first) idx))
-        s.add_ult(max_int, (second - first) * idx);
-        // s.add_ule((second - first) * idx, max_int);
-
-        // resolving disjunction via push/pop
-
-        // first disjunct: (= idx #x00000000)
-        s.push();
-        s.add_eq(idx);
-        s.check();
-        s.expect_unsat();
-        s.pop();
-
-        // second disjunct: (bvule (bvsub second first) q)
-        s.push();
-        s.add_ule(second - first, q);
-        s.check();
-        s.expect_unsat();
-        s.pop();
-    }
-
-    static void test_var_minimize(unsigned bw = 32) {
-        scoped_solver s(__func__);
-        auto x = s.var(s.add_var(bw));
-        auto y = s.var(s.add_var(bw));
-        auto z = s.var(s.add_var(bw));
-        s.add_eq(x);
-        s.add_eq(4 * y + 8 * z + x + 2); // should only depend on assignment to x
-        s.check();
-        s.expect_unsat();
-    }
-
-
-    /**
-     * x*x <= z
-     * (x+1)*(x+1) <= z
-     * y == x+1
-     * ¬(y*y <= z)
-     *
-     * The original version had signed comparisons but that doesn't matter for the UNSAT result.
-     * UNSAT can be seen easily by substituting the equality.
-     *
-     * Possible ways to solve:
-     * - Integrate AC congruence closure
-     *      See: Deepak Kapur. A Modular Associative Commutative (AC) Congruence Closure Algorithm, FSCD 2021. https://doi.org/10.4230/LIPIcs.FSCD.2021.15
-     * - Propagate equalities as substitutions
-     *      x=t /\ p(x) ==> p(t)
-     *      Ackermann-like reduction
-     *   (index, watch lists over boolean literals)
-     * - Augment explain:
-     *      conflict: y=x+1 /\ y^2 > z
-     *      explain could then derive (x+1)^2 > z
-     */
-    static void test_subst(unsigned bw = 32) {
-        scoped_solver s(__func__);
-        auto x = s.var(s.add_var(bw));
-        auto y = s.var(s.add_var(bw));
-        auto z = s.var(s.add_var(bw));
-        s.add_ule(x * x, z);  // optional
-        s.add_ule((x + 1) * (x + 1), z);
-        s.add_eq(x + 1 - y);
-        s.add_ult(z, y*y);
-        s.check();
-        s.expect_unsat();
-    }
-
-    static void test_subst_signed(unsigned bw = 32) {
-        scoped_solver s(__func__);
-        auto x = s.var(s.add_var(bw));
-        auto y = s.var(s.add_var(bw));
-        auto z = s.var(s.add_var(bw));
-        s.add_sle(x * x, z);  // optional
-        s.add_sle((x + 1) * (x + 1), z);
-        s.add_eq(x + 1 - y);
-        s.add_slt(z, y*y);
-        s.check();
-        s.expect_unsat();
-    }
-
-    // xy < xz and !Omega(x*y) => y < z
-    static void test_ineq_axiom1(unsigned bw = 32, std::optional<unsigned> perm = std::nullopt) {
-        if (perm) {
-            scoped_solver s(std::string(__func__) + " perm=" + std::to_string(*perm));
-            auto const bound = rational::power_of_two(bw/2);
-            auto x = s.var(s.add_var(bw));
-            auto y = s.var(s.add_var(bw));
-            auto z = s.var(s.add_var(bw));
-            permute_args(*perm, x, y, z);
-            s.add_ult(x * y, x * z);
-            s.add_ule(z, y);
-            s.add_ult(x, bound);
-            s.add_ult(y, bound);
-            s.check();
-            s.expect_unsat();
-        } else {
-            for (unsigned i = 0; i < 6; ++i) {
-                test_ineq_axiom1(bw, i);
-            }
-        }
-    }
-
-    static void test_ineq_non_axiom1(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw - 1);
-
-        for (unsigned i = 0; i < 6; ++i) {
+        static void test_l1() {
             scoped_solver s(__func__);
-            auto x = s.var(s.add_var(bw));
-            auto y = s.var(s.add_var(bw));
-            auto z = s.var(s.add_var(bw));
-            permute_args(i, x, y, z);
-            s.add_ult(x * y, x * z);
-            s.add_ule(z, y);
-            //s.add_ult(x, bound);
-            //s.add_ult(y, bound);
+            auto a = s.var(s.add_var(2));
+            s.add_eq(a + 1);
+            s.check();
+            s.expect_sat({{a, 3}});
+        }
+    
+        static void test_l2() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(2));
+            auto b = s.var(s.add_var(2));
+            s.add_eq(2*a + b + 1);
+            s.add_eq(2*b + a);
+            s.check();
+            s.expect_sat({{a, 2}, {b, 3}});
+        }
+
+        static void test_l3() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(2));
+            auto b = s.var(s.add_var(2));
+            s.add_eq(3*b + a + 2);
             s.check();
             s.expect_sat();
         }
-    }
 
-    // xy <= xz & !Omega(x*y) => y <= z or x = 0
-    static void test_ineq_axiom2(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw/2);
-        for (unsigned i = 0; i < 6; ++i) {
+        static void test_l4() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(3));
+            s.add_eq(4*a + 2);
+            s.check();
+            s.expect_unsat();
+        }
+
+        static void test_l5() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(3));
+            auto b = s.var(s.add_var(3));
+            s.add_diseq(b);
+            s.add_eq(a + 2*b + 4);
+            s.add_eq(a + 4*b + 4);
+            s.check();
+            s.expect_sat({{a, 4}, {b, 4}});
+        }
+
+    
+        /**
+         * This one is unsat because a*a*(a*a - 1)
+         * is 0 for all values of a.
+         */
+        static void test_p1() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(2));
+            auto p = a*a*(a*a - 1) + 1;
+            s.add_eq(p);
+            s.check();
+            s.expect_unsat();
+        }
+
+        /**
+         * has solutions a = 2 and a = 3
+         */
+        static void test_p2() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(2));
+            auto p = a*(a-1) + 2;
+            s.add_eq(p);
+            s.check();
+            s.expect_sat();
+        }
+
+        /**
+         * unsat
+         * - learns 3*x + 1 == 0 by polynomial resolution
+         * - this forces x == 5, which means the first constraint is unsatisfiable by parity.
+         */
+        static void test_p3() {
+            scoped_solver s(__func__);
+            auto x = s.var(s.add_var(4));
+            auto y = s.var(s.add_var(4));
+            auto z = s.var(s.add_var(4));
+            s.add_eq(x*x*y + 3*y + 7);
+            s.add_eq(2*y + z + 8);
+            s.add_eq(3*x + 4*y*z + 2*z*z + 1);
+            s.check();
+            s.expect_unsat();
+        }
+
+
+        // Unique solution: u = 5
+        static void test_ineq_basic1() {
+            scoped_solver s(__func__);
+            auto u = s.var(s.add_var(4));
+            s.add_ule(u, 5);
+            s.add_ule(5, u);
+            s.check();
+            s.expect_sat({{u, 5}});
+        }
+
+        // Unsatisfiable
+        static void test_ineq_basic2() {
+            scoped_solver s(__func__);
+            auto u = s.var(s.add_var(4));
+            s.add_ult(u, 5);
+            s.add_ule(5, u);
+            s.check();
+            s.expect_unsat();
+        }
+
+        // Solutions with u = v = w
+        static void test_ineq_basic3() {
+            scoped_solver s(__func__);
+            auto u = s.var(s.add_var(4));
+            auto v = s.var(s.add_var(4));
+            auto w = s.var(s.add_var(4));
+            s.add_ule(u, v);
+            s.add_ule(v, w);
+            s.add_ule(w, u);
+            s.check();
+            s.expect_sat();
+            SASSERT_EQ(s.get_value(u.var()), s.get_value(v.var()));
+            SASSERT_EQ(s.get_value(u.var()), s.get_value(w.var()));
+        }
+
+        // Unsatisfiable
+        static void test_ineq_basic4() {
+            scoped_solver s(__func__);
+            auto u = s.var(s.add_var(4));
+            auto v = s.var(s.add_var(4));
+            auto w = s.var(s.add_var(4));
+            s.add_ule(u, v);
+            s.add_ult(v, w);
+            s.add_ule(w, u);
+            s.check();
+            s.expect_unsat();
+        }
+
+        // Satisfiable
+        // Without forbidden intervals, we just try values for u until it works
+        static void test_ineq_basic5() {
+            scoped_solver s(__func__);
+            auto u = s.var(s.add_var(4));
+            auto v = s.var(s.add_var(4));
+            s.add_ule(12, u + v);
+            s.add_ule(v, 2);
+            s.check();
+            s.expect_sat();  // e.g., u = 12, v = 0
+        }
+
+        // Like test_ineq_basic5 but the other forbidden interval will be the longest
+        static void test_ineq_basic6() {
+            scoped_solver s(__func__);
+            auto u = s.var(s.add_var(4));
+            auto v = s.var(s.add_var(4));
+            s.add_ule(14, u + v);
+            s.add_ule(v, 2);
+            s.check();
+            s.expect_sat();
+        }
+
+
+        /**
+         * Check unsat of:
+         * u = v*q + r
+         * r < u
+         * v*q > u
+         */
+        static void test_ineq1() {
+            scoped_solver s(__func__);
+            auto u = s.var(s.add_var(5));
+            auto v = s.var(s.add_var(5));
+            auto q = s.var(s.add_var(5));
+            auto r = s.var(s.add_var(5));
+            s.add_eq(u - (v*q) - r);
+            s.add_ult(r, u);
+            s.add_ult(u, v*q);
+            s.check();
+            s.expect_unsat();
+        }
+
+        /**
+         * Check unsat of:
+         * n*q1 = a - b
+         * n*q2 + r2 = c*a - c*b
+         * n > r2 > 0
+         */
+        static void test_ineq2() {
+            scoped_solver s(__func__);
+            auto n = s.var(s.add_var(5));
+            auto q1 = s.var(s.add_var(5));
+            auto a = s.var(s.add_var(5));
+            auto b = s.var(s.add_var(5));
+            auto c = s.var(s.add_var(5));
+            auto q2 = s.var(s.add_var(5));
+            auto r2 = s.var(s.add_var(5));
+            s.add_eq(n*q1 - a + b);
+            s.add_eq(n*q2 + r2 - c*a + c*b);
+            s.add_ult(r2, n);
+            s.add_diseq(r2);
+            s.check();
+            s.expect_unsat();
+        }
+
+
+        /**
+         * Monotonicity example from certora
+         * 
+         * We do overflow checks by doubling the base bitwidth here.
+         */
+        static void test_monot(unsigned base_bw = 5) {
+            scoped_solver s(std::string{__func__} + "(" + std::to_string(base_bw) + ")");
+
+            auto max_int_const = rational::power_of_two(base_bw) - 1;
+
+            unsigned bw = 2 * base_bw;
+            auto max_int = s.var(s.add_var(bw));
+            s.add_eq(max_int + (-max_int_const));
+
+            auto tb1 = s.var(s.add_var(bw));
+            s.add_ule(tb1, max_int);
+            auto tb2 = s.var(s.add_var(bw));
+            s.add_ule(tb2, max_int);
+            auto a = s.var(s.add_var(bw));
+            s.add_ule(a, max_int);
+            auto v = s.var(s.add_var(bw));
+            s.add_ule(v, max_int);
+            auto base1 = s.var(s.add_var(bw));
+            s.add_ule(base1, max_int);
+            auto base2 = s.var(s.add_var(bw));
+            s.add_ule(base2, max_int);
+            auto elastic1 = s.var(s.add_var(bw));
+            s.add_ule(elastic1, max_int);
+            auto elastic2 = s.var(s.add_var(bw));
+            s.add_ule(elastic2, max_int);
+            auto err = s.var(s.add_var(bw));
+            s.add_ule(err, max_int);
+
+            auto rem1 = s.var(s.add_var(bw));
+            auto quot2 = s.var(s.add_var(bw));
+            s.add_ule(quot2, max_int);
+            auto rem2 = s.var(s.add_var(bw));
+            auto rem3 = s.var(s.add_var(bw));
+            auto quot4 = s.var(s.add_var(bw));
+            s.add_ule(quot4, max_int);
+            auto rem4 = s.var(s.add_var(bw));
+
+            s.add_diseq(elastic1);
+
+            // division: tb1 = (v * base1) / elastic1;
+            s.add_eq((tb1 * elastic1) + rem1 - (v * base1));
+            s.add_ult(rem1, elastic1);
+            s.add_ule((tb1 * elastic1), max_int);
+
+            // division: quot2 = (a * base1) / elastic1
+            s.add_eq((quot2 * elastic1) + rem2 - (a * base1));
+            s.add_ult(rem2, elastic1);
+            s.add_ule((quot2 * elastic1), max_int);
+
+            s.add_eq(base1 + quot2 - base2);
+
+            s.add_eq(elastic1 + a - elastic2);
+
+            // division: tb2 = ((v * base2) / elastic2);
+            s.add_eq((tb2 * elastic2) + rem3 - (v * base2));
+            s.add_ult(rem3, elastic2);
+            s.add_ule((tb2 * elastic2), max_int);
+
+            // division: quot4 = v / elastic2;
+            s.add_eq((quot4 * elastic2) + rem4 - v);
+            s.add_ult(rem4, elastic2);
+            s.add_ule((quot4 * elastic2), max_int);
+
+            s.add_eq(quot4 + 1 - err);
+
+            s.push();
+            s.add_ult(tb1, tb2);
+            s.check();
+            s.expect_unsat();
+            s.pop();
+
+            s.push();
+            s.add_ult(tb2 + err, tb1);
+            s.check();
+            s.expect_unsat();
+            s.pop();
+        }
+
+        /*
+         * Mul-then-div in fixed point arithmetic is (roughly) neutral.
+         *
+         * I.e. we prove "(((a * b) / sf) * sf) / b" to be equal to a, up to some error margin.
+         * 
+         * sf is the scaling factor (we could leave this unconstrained, but non-zero, to make the benchmark a bit harder)
+         * em is the error margin
+         * 
+         * We do overflow checks by doubling the base bitwidth here.
+         */
+        static void test_fixed_point_arith_mul_div_inverse() {
+            scoped_solver s(__func__);
+        
+            auto baseBw = 5;
+            auto max_int_const = 31; // (2^5 - 1) -- change this when you change baseBw
+
+            auto bw = 2 * baseBw;
+            auto max_int = s.var(s.add_var(bw));
+            s.add_eq(max_int - max_int_const);
+
+            // "input" variables
+            auto a = s.var(s.add_var(bw));
+            s.add_ule(a, max_int);
+            auto b = s.var(s.add_var(bw));
+            s.add_ule(b, max_int);
+            s.add_ult(0, b); // b > 0
+
+            // scaling factor (setting it, somewhat arbitrarily, to max_int/3)
+            auto sf = s.var(s.add_var(bw));
+            s.add_eq(sf - (max_int_const/3));
+
+            // (a * b) / sf = quot1 <=> quot1 * sf + rem1 - (a * b) = 0
+            auto quot1 = s.var(s.add_var(bw));
+            auto rem1 = s.var(s.add_var(bw));
+            s.add_eq((quot1 * sf) + rem1 - (a * b));
+            s.add_ult(rem1, sf);
+            s.add_ule(quot1 * sf, max_int);
+
+            // (((a * b) / sf) * sf) / b <=> quot2 * b + rem2 - (((a * b) / sf) * sf) = 0
+            auto quot2 = s.var(s.add_var(bw));
+            auto rem2 = s.var(s.add_var(bw));
+            s.add_eq((quot2 * b) + rem2 - (quot1 * sf));
+            s.add_ult(rem2, b);
+            s.add_ule(quot2 * b, max_int);
+
+            // sf / b = quot3 <=> quot3 * b + rem3 = sf
+            auto quot3 = s.var(s.add_var(bw));
+            auto rem3 = s.var(s.add_var(bw));
+            s.add_eq((quot3 * b) + rem3 - sf);
+            s.add_ult(rem3, b);
+            s.add_ule(quot3 * b, max_int);
+
+            // em = sf / b + 1
+            auto em = s.var(s.add_var(bw));
+            s.add_eq(quot3 + 1 - em);
+
+            // we prove quot3 <= a and quot3 + em >= a
+
+            s.push();
+            s.add_ult(a, quot3);
+            s.check();
+            s.expect_unsat();
+            s.pop();
+
+
+            // s.push();
+            // s.add_ult(quot3 + em, a);
+            // s.check();
+            // s.expect_unsat();
+            // s.pop();
+        }
+
+        /*
+         * Div-then-mul in fixed point arithmetic is (roughly) neutral.
+         *
+         * I.e. we prove "(b * ((a * sf) / b)) / sf" to be equal to a, up to some error margin.
+         * 
+         * sf is the scaling factor (we could leave this unconstrained, but non-zero, to make the benchmark a bit harder)
+         * em is the error margin
+         * 
+         * We do overflow checks by doubling the base bitwidth here.
+         */
+        static void test_fixed_point_arith_div_mul_inverse(unsigned base_bw = 5) {
+            scoped_solver s(__func__);
+
+            auto max_int_const = rational::power_of_two(base_bw) - 1;
+        
+            auto bw = 2 * base_bw;
+            auto max_int = s.var(s.add_var(bw));
+            s.add_eq(max_int - max_int_const);
+
+            // "input" variables
+            auto a = s.var(s.add_var(bw));
+            s.add_ule(a, max_int);
+            auto b = s.var(s.add_var(bw));
+            s.add_ule(b, max_int);
+            s.add_ult(0, b); // b > 0
+
+            // scaling factor (setting it, somewhat arbitrarily, to max_int/3)
+            auto sf = s.var(s.add_var(bw));
+            s.add_eq(sf - floor(max_int_const/3));
+
+            // (a * sf) / b = quot1 <=> quot1 * b + rem1 - (a * sf) = 0
+            auto quot1 = s.var(s.add_var(bw));
+            auto rem1 = s.var(s.add_var(bw));
+            s.add_eq((quot1 * b) + rem1 - (a * sf));
+            s.add_ult(rem1, b);
+            s.add_ule(quot1 * b, max_int);
+
+            // (b * ((a * sf) / b)) / sf = quot2 <=> quot2 * sf + rem2 - (b * ((a * sf) / b)) = 0
+            auto quot2 = s.var(s.add_var(bw));
+            auto rem2 = s.var(s.add_var(bw));
+            s.add_eq((quot2 * sf) + rem2 - (b * quot1));
+            s.add_ult(rem2, sf);
+            s.add_ule(quot2 * sf, max_int);
+
+            // b / sf = quot3 <=> quot3 * sf + rem3 - b = 0
+            auto quot3 = s.var(s.add_var(bw));
+            auto rem3 = s.var(s.add_var(bw));
+            s.add_eq((quot3 * sf) + rem3 - b);
+            s.add_ult(rem3, sf);
+            s.add_ule(quot3 * sf, max_int);
+
+            // em = b / sf + 1
+            auto em = s.var(s.add_var(bw));
+            s.add_eq(quot3 + 1 - em);
+
+            // we prove quot3 <= a and quot3 + em >= a
+
+            s.push();
+            s.add_ult(quot3 + em, a);
+            s.check();
+            //        s.expect_unsat();
+            s.pop();
+
+            s.push();
+            s.add_ult(a, quot3);
+            s.check();
+            //        s.expect_unsat();
+            s.pop();
+
+
+
+            //exit(0);
+        }
+
+        /** Monotonicity under bounds,
+         * puzzle extracted from https://github.com/NikolajBjorner/polysat/blob/main/puzzles/bv.smt2
+         *
+         * x, y, z \in [0..2^64[
+         * x, y, z < 2^32
+         * y <= x
+         * x*z < 2^32
+         * y*z >= 2^32
+         */
+        static void test_monot_bounds(unsigned base_bw = 32) {
+            scoped_solver s(std::string{__func__} + "(" + std::to_string(base_bw) + ")");
+            unsigned bw = 2 * base_bw;
+            auto y = s.var(s.add_var(bw));
+            auto z = s.var(s.add_var(bw));
+            auto x = s.var(s.add_var(bw));
+            auto bound = rational::power_of_two(base_bw);
+#if 1
+            s.add_ult(x, bound);
+            s.add_ult(y, bound);
+            // s.add_ult(z, bound);   // not required
+#else
+            s.add_ule(x, bound - 1);
+            s.add_ule(y, bound - 1);
+            // s.add_ule(z, bound - 1);   // not required
+#endif
+            unsigned a = 13;
+            s.add_ule(z, y);
+            s.add_ult(x*y, a);
+            s.add_ule(a, x*z);
+            s.check();
+            s.expect_unsat();
+        }
+
+        /** Monotonicity under bounds, simplified even more.
+         *
+         * x, y, z \in [0..2^64[
+         * x, y, z < 2^32
+         * z <= y
+         * y*x < z*x
+         */
+        static void test_monot_bounds_simple(unsigned base_bw = 32) {
+            scoped_solver s(__func__);
+            unsigned bw = 2 * base_bw;
+            /*
+              auto z = s.var(s.add_var(bw));
+              auto x = s.var(s.add_var(bw));
+              auto y = s.var(s.add_var(bw));
+            */
+            auto y = s.var(s.add_var(bw));
+            auto x = s.var(s.add_var(bw));
+            auto z = s.var(s.add_var(bw));
+            auto bound = rational::power_of_two(base_bw);
+            s.add_ult(x, bound);
+            s.add_ult(y, bound);
+            s.add_ult(z, bound);
+            s.add_ule(z, y);
+            s.add_ult(y*x, z*x);
+            s.check();
+            s.expect_unsat();
+        }
+
+        /*
+         * Transcribed from https://github.com/NikolajBjorner/polysat/blob/main/puzzles/bv.smt2
+         *
+         * We do overflow checks by doubling the base bitwidth here.
+         */
+        static void test_monot_bounds_full(unsigned base_bw = 5) {
+            scoped_solver s(__func__);
+        
+            auto const max_int_const = rational::power_of_two(base_bw) - 1;
+        
+            auto const bw = 2 * base_bw;
+            auto const max_int = s.var(s.add_var(bw));
+            s.add_eq(max_int - max_int_const);
+
+            auto const first = s.var(s.add_var(bw));
+            s.add_ule(first, max_int);
+            auto const second = s.var(s.add_var(bw));
+            s.add_ule(second, max_int);
+            auto const idx = s.var(s.add_var(bw));
+            s.add_ule(idx, max_int);
+            auto const q = s.var(s.add_var(bw));
+            s.add_ule(q, max_int);
+            auto const r = s.var(s.add_var(bw));
+            s.add_ule(r, max_int);
+
+            // q = max_int / idx <=> q * idx + r - max_int = 0
+            s.add_eq((q * idx) + r - max_int);
+            s.add_ult(r, idx);
+            s.add_ule(q * idx, max_int);
+
+            /*  last assertion:
+                (not
+                (=> (bvugt second first) 
+                (=> 
+                (=> (not (= idx #x00000000)) 
+                (bvule (bvsub second first) q)) 
+                (bvumul_noovfl (bvsub second first) idx))))
+                transforming negated boolean skeleton:
+                (not (=> a (=> (or b c) d))) <=> (and a (not d) (or b c))
+            */
+
+            // (bvugt second first)
+            s.add_ult(first, second);
+            // (not (bvumul_noovfl (bvsub second first) idx))
+            s.add_ult(max_int, (second - first) * idx);
+            // s.add_ule((second - first) * idx, max_int);
+
+            // resolving disjunction via push/pop
+
+            // first disjunct: (= idx #x00000000)
+            s.push();
+            s.add_eq(idx);
+            s.check();
+            s.expect_unsat();
+            s.pop();
+
+            // second disjunct: (bvule (bvsub second first) q)
+            s.push();
+            s.add_ule(second - first, q);
+            s.check();
+            s.expect_unsat();
+            s.pop();
+        }
+
+        static void test_var_minimize(unsigned bw = 32) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
             auto y = s.var(s.add_var(bw));
             auto z = s.var(s.add_var(bw));
-            permute_args(i, x, y, z);
-            s.add_ult(x * y, x * z);
-            s.add_ult(z, y);
-            s.add_ult(x, bound);
-            s.add_ult(y, bound);
-            s.add_diseq(x);
+            s.add_eq(x);
+            s.add_eq(4 * y + 8 * z + x + 2); // should only depend on assignment to x
             s.check();
             s.expect_unsat();
         }
-    }
 
-    // xy < b & a <= x  & !Omega(x*y) => a*y < b
-    static void test_ineq_axiom3(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw/2);
-        for (unsigned i = 0; i < 24; ++i) {
+
+        /**
+         * x*x <= z
+         * (x+1)*(x+1) <= z
+         * y == x+1
+         * ¬(y*y <= z)
+         *
+         * The original version had signed comparisons but that doesn't matter for the UNSAT result.
+         * UNSAT can be seen easily by substituting the equality.
+         *
+         * Possible ways to solve:
+         * - Integrate AC congruence closure
+         *      See: Deepak Kapur. A Modular Associative Commutative (AC) Congruence Closure Algorithm, FSCD 2021. https://doi.org/10.4230/LIPIcs.FSCD.2021.15
+         * - Propagate equalities as substitutions
+         *      x=t /\ p(x) ==> p(t)
+         *      Ackermann-like reduction
+         *   (index, watch lists over boolean literals)
+         * - Augment explain:
+         *      conflict: y=x+1 /\ y^2 > z
+         *      explain could then derive (x+1)^2 > z
+         */
+        static void test_subst(unsigned bw = 32) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
             auto y = s.var(s.add_var(bw));
-            auto a = s.var(s.add_var(bw));
-            auto b = s.var(s.add_var(bw));
-            permute_args(i, x, y, a, b);
-            s.add_ult(x * y, b);
-            s.add_ule(a, x);
-            s.add_ult(x, bound);
-            s.add_ult(y, bound);
-            s.add_ule(b, a * y);
+            auto z = s.var(s.add_var(bw));
+            s.add_ule(x * x, z);  // optional
+            s.add_ule((x + 1) * (x + 1), z);
+            s.add_eq(x + 1 - y);
+            s.add_ult(z, y*y);
             s.check();
             s.expect_unsat();
         }
-    }
 
-    // x*y <= b & a <= x & !Omega(x*y) => a*y <= b
-    static void test_ineq_axiom4(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw/2);
-        for (unsigned i = 0; i < 24; ++i) {
+        static void test_subst_signed(unsigned bw = 32) {
             scoped_solver s(__func__);
             auto x = s.var(s.add_var(bw));
             auto y = s.var(s.add_var(bw));
+            auto z = s.var(s.add_var(bw));
+            s.add_sle(x * x, z);  // optional
+            s.add_sle((x + 1) * (x + 1), z);
+            s.add_eq(x + 1 - y);
+            s.add_slt(z, y*y);
+            s.check();
+            s.expect_unsat();
+        }
+
+        // xy < xz and !Omega(x*y) => y < z
+        static void test_ineq_axiom1(unsigned bw = 32, std::optional<unsigned> perm = std::nullopt) {
+            if (perm) {
+                scoped_solver s(std::string(__func__) + " perm=" + std::to_string(*perm));
+                auto const bound = rational::power_of_two(bw/2);
+                auto x = s.var(s.add_var(bw));
+                auto y = s.var(s.add_var(bw));
+                auto z = s.var(s.add_var(bw));
+                permute_args(*perm, x, y, z);
+                s.add_ult(x * y, x * z);
+                s.add_ule(z, y);
+                s.add_ult(x, bound);
+                s.add_ult(y, bound);
+                s.check();
+                s.expect_unsat();
+            } 
+            else {
+                for (unsigned i = 0; i < 6; ++i) {
+                    test_ineq_axiom1(bw, i);
+                }
+            }
+        }
+
+        static void test_ineq_non_axiom1(unsigned bw = 32) {
+            auto const bound = rational::power_of_two(bw - 1);
+
+            for (unsigned i = 0; i < 6; ++i) {
+                scoped_solver s(__func__);
+                auto x = s.var(s.add_var(bw));
+                auto y = s.var(s.add_var(bw));
+                auto z = s.var(s.add_var(bw));
+                permute_args(i, x, y, z);
+                s.add_ult(x * y, x * z);
+                s.add_ule(z, y);
+                //s.add_ult(x, bound);
+                //s.add_ult(y, bound);
+                s.check();
+                s.expect_sat();
+            }
+        }
+
+        // xy <= xz & !Omega(x*y) => y <= z or x = 0
+        static void test_ineq_axiom2(unsigned bw = 32) {
+            auto const bound = rational::power_of_two(bw/2);
+            for (unsigned i = 0; i < 6; ++i) {
+                scoped_solver s(__func__);
+                auto x = s.var(s.add_var(bw));
+                auto y = s.var(s.add_var(bw));
+                auto z = s.var(s.add_var(bw));
+                permute_args(i, x, y, z);
+                s.add_ult(x * y, x * z);
+                s.add_ult(z, y);
+                s.add_ult(x, bound);
+                s.add_ult(y, bound);
+                s.add_diseq(x);
+                s.check();
+                s.expect_unsat();
+            }
+        }
+
+        // xy < b & a <= x  & !Omega(x*y) => a*y < b
+        static void test_ineq_axiom3(unsigned bw = 32) {
+            auto const bound = rational::power_of_two(bw/2);
+            for (unsigned i = 0; i < 24; ++i) {
+                scoped_solver s(__func__);
+                auto x = s.var(s.add_var(bw));
+                auto y = s.var(s.add_var(bw));
+                auto a = s.var(s.add_var(bw));
+                auto b = s.var(s.add_var(bw));
+                permute_args(i, x, y, a, b);
+                s.add_ult(x * y, b);
+                s.add_ule(a, x);
+                s.add_ult(x, bound);
+                s.add_ult(y, bound);
+                s.add_ule(b, a * y);
+                s.check();
+                s.expect_unsat();
+            }
+        }
+
+        // x*y <= b & a <= x & !Omega(x*y) => a*y <= b
+        static void test_ineq_axiom4(unsigned bw = 32) {
+            auto const bound = rational::power_of_two(bw/2);
+            for (unsigned i = 0; i < 24; ++i) {
+                scoped_solver s(__func__);
+                auto x = s.var(s.add_var(bw));
+                auto y = s.var(s.add_var(bw));
+                auto a = s.var(s.add_var(bw));
+                auto b = s.var(s.add_var(bw));
+                permute_args(i, x, y, a, b);
+                s.add_ule(x * y, b);
+                s.add_ule(a, x);
+                s.add_ult(x, bound);
+                s.add_ult(y, bound);
+                s.add_ult(b, a * y);
+                s.check();
+                s.expect_unsat();
+            }
+        }
+
+        // x*y <= b & a <= x & !Omega(x*y) => a*y <= b
+        static void test_ineq_non_axiom4(unsigned bw, unsigned i) {
+            auto const bound = rational::power_of_two(bw - 1);
+            scoped_solver s(__func__);
+            LOG("permutation: " << i);
+            auto x = s.var(s.add_var(bw));
+            auto y = s.var(s.add_var(bw));
             auto a = s.var(s.add_var(bw));
             auto b = s.var(s.add_var(bw));
             permute_args(i, x, y, a, b);
@@ -889,310 +910,292 @@ public:
             s.add_ult(y, bound);
             s.add_ult(b, a * y);
             s.check();
-            s.expect_unsat();
+            s.expect_sat();
         }
-    }
-
-    // x*y <= b & a <= x & !Omega(x*y) => a*y <= b
-    static void test_ineq_non_axiom4(unsigned bw, unsigned i) {
-        auto const bound = rational::power_of_two(bw - 1);
-        scoped_solver s(__func__);
-        LOG("permutation: " << i);
-        auto x = s.var(s.add_var(bw));
-        auto y = s.var(s.add_var(bw));
-        auto a = s.var(s.add_var(bw));
-        auto b = s.var(s.add_var(bw));
-        permute_args(i, x, y, a, b);
-        s.add_ule(x * y, b);
-        s.add_ule(a, x);
-        s.add_ult(x, bound);
-        s.add_ult(y, bound);
-        s.add_ult(b, a * y);
-        s.check();
-        s.expect_sat();
-    }
     
-    static void test_ineq_non_axiom4(unsigned bw = 32) {
-        for (unsigned i = 0; i < 24; ++i)
-            test_ineq_non_axiom4(bw, i);
-    }
+        static void test_ineq_non_axiom4(unsigned bw = 32) {
+            for (unsigned i = 0; i < 24; ++i)
+                test_ineq_non_axiom4(bw, i);
+        }
 
-    // a < xy & x <= b  & !Omega(x*y) => a < b*y
-    static void test_ineq_axiom5(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw/2);
-        for (unsigned i = 0; i < 24; ++i) {
+        // a < xy & x <= b  & !Omega(x*y) => a < b*y
+        static void test_ineq_axiom5(unsigned bw = 32) {
+            auto const bound = rational::power_of_two(bw/2);
+            for (unsigned i = 0; i < 24; ++i) {
+                scoped_solver s(__func__);
+                auto x = s.var(s.add_var(bw));
+                auto y = s.var(s.add_var(bw));
+                auto a = s.var(s.add_var(bw));
+                auto b = s.var(s.add_var(bw));
+                permute_args(i, x, y, a, b);
+                s.add_ult(a, x * y);
+                s.add_ule(x, b);
+                s.add_ult(x, bound);
+                s.add_ult(y, bound);
+                s.add_ule(b * y, a);
+                s.check();
+                s.expect_unsat();
+            }
+        }
+
+        // a <= xy & x <= b & !Omega(x*y) => a <= b*y
+        static void test_ineq_axiom6(unsigned bw = 32) {
+            auto const bound = rational::power_of_two(bw/2);
+            for (unsigned i = 0; i < 24; ++i) {
+                scoped_solver s(__func__);
+                auto x = s.var(s.add_var(bw));
+                auto y = s.var(s.add_var(bw));
+                auto a = s.var(s.add_var(bw));
+                auto b = s.var(s.add_var(bw));
+                permute_args(i, x, y, a, b);
+                s.add_ule(a, x * y);
+                s.add_ule(x, b);
+                s.add_ult(x, bound);
+                s.add_ult(y, bound);
+                s.add_ult(b * y, a);
+                s.check();
+                s.expect_unsat();
+            }
+        }
+
+        static void test_quot_rem_incomplete() {
+            unsigned bw = 4;
             scoped_solver s(__func__);
-            auto x = s.var(s.add_var(bw));
-            auto y = s.var(s.add_var(bw));
-            auto a = s.var(s.add_var(bw));
-            auto b = s.var(s.add_var(bw));
-            permute_args(i, x, y, a, b);
-            s.add_ult(a, x * y);
-            s.add_ule(x, b);
-            s.add_ult(x, bound);
-            s.add_ult(y, bound);
-            s.add_ule(b * y, a);
+            s.set_max_conflicts(5);
+            auto quot = s.var(s.add_var(bw));
+            auto rem = s.var(s.add_var(bw));
+            auto a = s.value(rational(2), bw);
+            auto b = s.value(rational(5), bw);
+            // Incomplete axiomatization of quotient/remainder.
+            // quot_rem(2, 5) should have single solution (0, 2),
+            // but with the usual axioms we also get (3, 3).
+            s.add_eq(b * quot + rem - a);
+            s.add_noovfl(b, quot);
+            s.add_ult(rem, b);
+            // To force a solution that's different from the correct one.
+            s.add_diseq(quot - 0);
+            s.check();
+            s.expect_sat({{quot, 3}, {rem, 3}});
+        }
+
+        static void test_quot_rem_fixed() {
+            unsigned bw = 4;
+            scoped_solver s(__func__);
+            s.set_max_conflicts(5);
+            auto a = s.value(rational(2), bw);
+            auto b = s.value(rational(5), bw);
+            auto [quot, rem] = s.quot_rem(a, b);
+            s.add_diseq(quot - 0);  // to force a solution that's different from the correct one
             s.check();
             s.expect_unsat();
         }
-    }
 
-    // a <= xy & x <= b & !Omega(x*y) => a <= b*y
-    static void test_ineq_axiom6(unsigned bw = 32) {
-        auto const bound = rational::power_of_two(bw/2);
-        for (unsigned i = 0; i < 24; ++i) {
+        static void test_quot_rem(unsigned bw = 32) {
             scoped_solver s(__func__);
-            auto x = s.var(s.add_var(bw));
-            auto y = s.var(s.add_var(bw));
+            s.set_max_conflicts(5);
             auto a = s.var(s.add_var(bw));
-            auto b = s.var(s.add_var(bw));
-            permute_args(i, x, y, a, b);
-            s.add_ule(a, x * y);
-            s.add_ule(x, b);
-            s.add_ult(x, bound);
-            s.add_ult(y, bound);
-            s.add_ult(b * y, a);
-            s.check();
-            s.expect_unsat();
-        }
-    }
-
-    static void test_quot_rem_incomplete() {
-        unsigned bw = 4;
-        scoped_solver s(__func__);
-        s.set_max_conflicts(5);
-        auto quot = s.var(s.add_var(bw));
-        auto rem = s.var(s.add_var(bw));
-        auto a = s.value(rational(2), bw);
-        auto b = s.value(rational(5), bw);
-        // Incomplete axiomatization of quotient/remainder.
-        // quot_rem(2, 5) should have single solution (0, 2),
-        // but with the usual axioms we also get (3, 3).
-        s.add_eq(b * quot + rem - a);
-        s.add_noovfl(b, quot);
-        s.add_ult(rem, b);
-        // To force a solution that's different from the correct one.
-        s.add_diseq(quot - 0);
-        s.check();
-        s.expect_sat({{quot, 3}, {rem, 3}});
-    }
-
-    static void test_quot_rem_fixed() {
-        unsigned bw = 4;
-        scoped_solver s(__func__);
-        s.set_max_conflicts(5);
-        auto a = s.value(rational(2), bw);
-        auto b = s.value(rational(5), bw);
-        auto [quot, rem] = s.quot_rem(a, b);
-        s.add_diseq(quot - 0);  // to force a solution that's different from the correct one
-        s.check();
-        s.expect_unsat();
-    }
-
-    static void test_quot_rem(unsigned bw = 32) {
-        scoped_solver s(__func__);
-        s.set_max_conflicts(5);
-        auto a = s.var(s.add_var(bw));
-        auto quot = s.var(s.add_var(bw));
-        auto rem = s.var(s.add_var(bw));
-        auto x = a * 123;
-        auto y = 123;
-        // quot = udiv(a*123, 123)
-        s.add_eq(quot * y + rem - x);
-        s.add_diseq(a - quot);
-        s.add_noovfl(quot, y);     
-        s.add_ult(rem, x);
-        s.check();
-        s.expect_sat();
-    }
-
-    static void test_quot_rem2(unsigned bw = 32) {
-        scoped_solver s(__func__);
-        s.set_max_conflicts(5);
-        auto q = s.var(s.add_var(bw));
-        auto r = s.var(s.add_var(bw));
-        auto idx = s.var(s.add_var(bw));
-        auto second = s.var(s.add_var(bw));
-        auto first = s.var(s.add_var(bw));
-        s.add_eq(q*idx + r, UINT_MAX);
-        s.add_ult(r, idx);
-        s.add_noovfl(q, idx);
-        s.add_ult(first, second);
-        s.add_diseq(idx, 0);
-        s.add_ule(second - first, q);
-        s.add_noovfl(second  - first, idx);
-        s.check();       
-    }
-
-    static void test_band(unsigned bw = 32) {
-        {
-            scoped_solver s(__func__);
-            auto p = s.var(s.add_var(bw));
-            auto q = s.var(s.add_var(bw));
-            s.add_diseq(p - s.band(p, q));
-            s.add_diseq(p - q);
+            auto quot = s.var(s.add_var(bw));
+            auto rem = s.var(s.add_var(bw));
+            auto x = a * 123;
+            auto y = 123;
+            // quot = udiv(a*123, 123)
+            s.add_eq(quot * y + rem - x);
+            s.add_diseq(a - quot);
+            s.add_noovfl(quot, y);     
+            s.add_ult(rem, x);
             s.check();
             s.expect_sat();
         }
-        {
-            scoped_solver s(__func__);
-            auto p = s.var(s.add_var(bw));
-            auto q = s.var(s.add_var(bw));
-            s.add_ult(p, s.band(p, q));
-            s.check();
-            s.expect_unsat();
-        }
-        {
-            scoped_solver s(__func__);
-            auto p = s.var(s.add_var(bw));
-            auto q = s.var(s.add_var(bw));
-            s.add_ult(q, s.band(p, q));
-            s.check();
-            s.expect_unsat();
-        }
-        {
-            scoped_solver s(__func__);
-            auto p = s.var(s.add_var(bw));
-            auto q = s.var(s.add_var(bw));
-            s.add_ule(p, s.band(p, q));
-            s.check();
-            s.expect_sat();
-        }
-        {
-            scoped_solver s(__func__);
-            auto p = s.var(s.add_var(bw));
-            auto q = s.var(s.add_var(bw));
-            s.add_ule(p, s.band(p, q));
-            s.add_diseq(p - s.band(p, q));
-            s.check();
-            s.expect_unsat();
-        }
-    }
 
-    static void test_fi_zero() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(256));
-        auto b = s.var(s.add_var(256));
-        auto c = s.var(s.add_var(256));
-        s.add_eq(a, 0);
-        s.add_eq(c, 0);  // add c to prevent simplification by leading coefficient
-        s.add_eq(4*a - 123456789*b + c);
-        s.check();
-        s.expect_sat({{a, 0}, {b, 0}});
-    }
-
-    static void test_fi_nonzero() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(5));
-        auto b = s.var(s.add_var(5));
-        s.add_ult(b*b*b, 7*a + 6);
-        s.check();
-    }
-
-    static void test_fi_nonmax() {
-        scoped_solver s(__func__);
-        auto a = s.var(s.add_var(5));
-        auto b = s.var(s.add_var(5));
-        s.add_ult(a + 8, b*b*b);
-        s.check();
-    }
-
-    static void test_fi_disequal_mild() {
-        {
-            // small version
+        static void test_quot_rem2(unsigned bw = 32) {
             scoped_solver s(__func__);
-            auto a = s.var(s.add_var(6));
-            auto b = s.var(s.add_var(6));
-            // v > -3*v
-            s.add_eq(a - 3);
-            s.add_ult(-a*b, b);
-            s.check();
+            s.set_max_conflicts(5);
+            auto q = s.var(s.add_var(bw));
+            auto r = s.var(s.add_var(bw));
+            auto idx = s.var(s.add_var(bw));
+            auto second = s.var(s.add_var(bw));
+            auto first = s.var(s.add_var(bw));
+            s.add_eq(q*idx + r, UINT_MAX);
+            s.add_ult(r, idx);
+            s.add_noovfl(q, idx);
+            s.add_ult(first, second);
+            s.add_diseq(idx, 0);
+            s.add_ule(second - first, q);
+            s.add_noovfl(second  - first, idx);
+            s.check();       
         }
-        {
-            // large version
+
+        static void test_band(unsigned bw = 32) {
+            {
+                scoped_solver s(__func__);
+                auto p = s.var(s.add_var(bw));
+                auto q = s.var(s.add_var(bw));
+                s.add_diseq(p - s.band(p, q));
+                s.add_diseq(p - q);
+                s.check();
+                s.expect_sat();
+            }
+            {
+                scoped_solver s(__func__);
+                auto p = s.var(s.add_var(bw));
+                auto q = s.var(s.add_var(bw));
+                s.add_ult(p, s.band(p, q));
+                s.check();
+                s.expect_unsat();
+            }
+            {
+                scoped_solver s(__func__);
+                auto p = s.var(s.add_var(bw));
+                auto q = s.var(s.add_var(bw));
+                s.add_ult(q, s.band(p, q));
+                s.check();
+                s.expect_unsat();
+            }
+            {
+                scoped_solver s(__func__);
+                auto p = s.var(s.add_var(bw));
+                auto q = s.var(s.add_var(bw));
+                s.add_ule(p, s.band(p, q));
+                s.check();
+                s.expect_sat();
+            }
+            {
+                scoped_solver s(__func__);
+                auto p = s.var(s.add_var(bw));
+                auto q = s.var(s.add_var(bw));
+                s.add_ule(p, s.band(p, q));
+                s.add_diseq(p - s.band(p, q));
+                s.check();
+                s.expect_unsat();
+            }
+        }
+
+        static void test_fi_zero() {
             scoped_solver s(__func__);
             auto a = s.var(s.add_var(256));
             auto b = s.var(s.add_var(256));
-            // v > -100*v
-            s.add_eq(a - 100);
-            s.add_ult(-a*b, b);
+            auto c = s.var(s.add_var(256));
+            s.add_eq(a, 0);
+            s.add_eq(c, 0);  // add c to prevent simplification by leading coefficient
+            s.add_eq(4*a - 123456789*b + c);
+            s.check();
+            s.expect_sat({{a, 0}, {b, 0}});
+        }
+
+        static void test_fi_nonzero() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(5));
+            auto b = s.var(s.add_var(5));
+            s.add_ult(b*b*b, 7*a + 6);
             s.check();
         }
-    }
 
-    // Goal: we probably mix up polysat variables and PDD variables at several points; try to uncover such cases
-    // NOTE: actually, add_var seems to keep them in sync, so this is not an issue at the moment (but we should still test it later)
-    // static void test_mixed_vars() {
-    //     scoped_solver s(__func__);
-    //     auto a = s.var(s.add_var(2));
-    //     auto b = s.var(s.add_var(4));
-    //     auto c = s.var(s.add_var(2));
-    //     s.add_eq(a + 2*c + 4);
-    //     s.add_eq(3*b + 4);
-    //     s.check();
-    //     // Expected result:
-    // }
-
-};  // class test_polysat
-
-
-/// Here we deal with linear constraints of the form
-///
-///     a1*x + b1 <= a2*x + b2   (mod m = 2^bw)
-///
-/// and their negation.
-class test_fi {
-
-    static bool is_violated(rational const& a1, rational const& b1, rational const& a2, rational const& b2, rational const& val, bool negated, rational const& m) {
-        rational const lhs = (a1*val + b1) % m;
-        rational const rhs = (a2*val + b2) % m;
-        if (negated)
-            return lhs <= rhs;
-        else
-            return lhs > rhs;
-    }
-
-    // Returns true if the input is valid and the test did useful work
-    static bool check_one(rational const& a1, rational const& b1, rational const& a2, rational const& b2, rational const& val, bool negated, unsigned bw) {
-        rational const m = rational::power_of_two(bw);
-        if (a1.is_zero() && a2.is_zero())
-            return false;
-        if (!is_violated(a1, b1, a2, b2, val, negated, m))
-            return false;
-
-        scoped_solver s(__func__);
-        auto x = s.var(s.add_var(bw));
-        signed_constraint c = s.ule(a1*x + b1, a2*x + b2);
-        if (negated)
-            c.negate();
-        viable& v = s.m_viable;
-        v.intersect(x.var(), c);
-        // Trigger forbidden interval refinement
-        v.is_viable(x.var(), val);
-        auto* e = v.m_units[x.var()];
-        if (!e) {
-            std::cout << "test_fi: no interval for a1=" << a1 << " b1=" << b1 << " a2=" << a2 << " b2=" << b2 << " val=" << val << " neg=" << negated << std::endl;
-            // VERIFY(false);
-            return false;
+        static void test_fi_nonmax() {
+            scoped_solver s(__func__);
+            auto a = s.var(s.add_var(5));
+            auto b = s.var(s.add_var(5));
+            s.add_ult(a + 8, b*b*b);
+            s.check();
         }
-        VERIFY(e);
-        auto* first = e;
-        SASSERT(e->next() == e);  // the result is expected to be a single interval (although for this check it doesn't really matter if there's more...)
-        do {
-            rational const& lo = e->interval.lo_val();
-            rational const& hi = e->interval.hi_val();
-            for (rational x = lo; x != hi; x = (x + 1) % m) {
-                // LOG("lo=" << lo << " hi=" << hi << " x=" << x);
-                if (!is_violated(a1, b1, a2, b2, val, negated, m)) {
-                    std::cout << "test_fi: unsound for a1=" << a1 << " b1=" << b1 << " a2=" << a2 << " b2=" << b2 << " val=" << val << " neg=" << negated << std::endl;
-                    VERIFY(false);
-                }
+
+        static void test_fi_disequal_mild() {
+            {
+                // small version
+                scoped_solver s(__func__);
+                auto a = s.var(s.add_var(6));
+                auto b = s.var(s.add_var(6));
+                // v > -3*v
+                s.add_eq(a - 3);
+                s.add_ult(-a*b, b);
+                s.check();
+            }
+            {
+                // large version
+                scoped_solver s(__func__);
+                auto a = s.var(s.add_var(256));
+                auto b = s.var(s.add_var(256));
+                // v > -100*v
+                s.add_eq(a - 100);
+                s.add_ult(-a*b, b);
+                s.check();
             }
-            e = e->next();
         }
-        while (e != first);
-        return true;
-    }
+    
+        // Goal: we probably mix up polysat variables and PDD variables at several points; try to uncover such cases
+        // NOTE: actually, add_var seems to keep them in sync, so this is not an issue at the moment (but we should still test it later)
+        // static void test_mixed_vars() {
+        //     scoped_solver s(__func__);
+        //     auto a = s.var(s.add_var(2));
+        //     auto b = s.var(s.add_var(4));
+        //     auto c = s.var(s.add_var(2));
+        //     s.add_eq(a + 2*c + 4);
+        //     s.add_eq(3*b + 4);
+        //     s.check();
+        //     // Expected result:
+        // }
+    
+    };  // class test_polysat
+    
+    
+    // Here we deal with linear constraints of the form
+    //
+    //     a1*x + b1 <= a2*x + b2   (mod m = 2^bw)
+    //
+    // and their negation.
+    
+    class test_fi {
+
+        static bool is_violated(rational const& a1, rational const& b1, rational const& a2, rational const& b2, 
+                                rational const& val, bool negated, rational const& m) {
+            rational const lhs = (a1*val + b1) % m;
+            rational const rhs = (a2*val + b2) % m;
+            if (negated)
+                return lhs <= rhs;
+            else
+                return lhs > rhs;
+        }
+
+        // Returns true if the input is valid and the test did useful work
+        static bool check_one(rational const& a1, rational const& b1, rational const& a2, rational const& b2, rational const& val, bool negated, unsigned bw) {
+            rational const m = rational::power_of_two(bw);
+            if (a1.is_zero() && a2.is_zero())
+                return false;
+            if (!is_violated(a1, b1, a2, b2, val, negated, m))
+                return false;
+            
+            scoped_solver s(__func__);
+            auto x = s.var(s.add_var(bw));
+            signed_constraint c = s.ule(a1*x + b1, a2*x + b2);
+            if (negated)
+                c.negate();
+            viable& v = s.m_viable;
+            v.intersect(x.var(), c);
+            // Trigger forbidden interval refinement
+            v.is_viable(x.var(), val);
+            auto* e = v.m_units[x.var()];
+            if (!e) {
+                std::cout << "test_fi: no interval for a1=" << a1 << " b1=" << b1 << " a2=" << a2 << " b2=" << b2 << " val=" << val << " neg=" << negated << std::endl;
+                // VERIFY(false);
+                return false;
+            }
+            VERIFY(e);
+            auto* first = e;
+            SASSERT(e->next() == e);  // the result is expected to be a single interval (although for this check it doesn't really matter if there's more...)
+            do {
+                rational const& lo = e->interval.lo_val();
+                rational const& hi = e->interval.hi_val();
+                for (rational x = lo; x != hi; x = (x + 1) % m) {
+                    // LOG("lo=" << lo << " hi=" << hi << " x=" << x);
+                    if (!is_violated(a1, b1, a2, b2, val, negated, m)) {
+                        std::cout << "test_fi: unsound for a1=" << a1 << " b1=" << b1 << " a2=" << a2 << " b2=" << b2 << " val=" << val << " neg=" << negated << std::endl;
+                        VERIFY(false);
+                    }
+                }
+                e = e->next();
+            }
+            while (e != first);
+            return true;
+        }
 
 public:
     static void exhaustive(unsigned bw = 0) {
@@ -1218,7 +1221,6 @@ public:
                 }
             }
         }
-    }
 
     static void randomized(unsigned num_rounds = 100'000, unsigned bw = 16) {
         std::cout << "test_fi::randomized for bw=" << bw << " (" << num_rounds << " rounds)" << std::endl;
@@ -1242,9 +1244,8 @@ public:
             if (useful)
                 round--;
         }
-    }
-
-};  // class test_fi
+        
+    };  // class test_fi
 
 
     // convert assertions into internal solver state
@@ -1375,20 +1376,20 @@ public:
         for (auto const& [k,v] : expr2pdd)
             dealloc(v);
     }
-
+    
 }  // namespace polysat
 
 
 void tst_polysat() {
     using namespace polysat;
 
+
     test_polysat::test_fi_zero();
     test_polysat::test_fi_nonzero();
     test_polysat::test_fi_nonmax();
     test_polysat::test_fi_disequal_mild();
-    return;
 
-    test_polysat::test_l2();
+
 
 #if 0
     // looks like a fishy conflict lemma?
@@ -1415,10 +1416,6 @@ void tst_polysat() {
 
 #endif
 
-    test_polysat::test_ineq_basic4();
-    //return;
-
-    test_polysat::test_ineq_basic6();
             
     // test_polysat::test_monot_bounds(8);
 
@@ -1442,13 +1439,13 @@ void tst_polysat() {
 
     test_polysat::test_cjust();
     test_polysat::test_subst();
-    test_polysat::test_monot_bounds(2);
 
     test_polysat::test_var_minimize();
 
     test_polysat::test_ineq1();
     test_polysat::test_ineq2();
     test_polysat::test_monot();
+    test_polysat::test_monot_bounds(2);
 
     return;