diff --git a/src/math/bigfix/u256.cpp b/src/math/bigfix/u256.cpp
index 35182b1d2..cd3dba633 100644
--- a/src/math/bigfix/u256.cpp
+++ b/src/math/bigfix/u256.cpp
@@ -1,3 +1,4 @@
+#include "util/mpn.h"
 #include "math/bigfix/u256.h"
 #include "math/bigfix/Hacl_Bignum256.h"
 #include <memory>
@@ -13,6 +14,13 @@ u256::u256(uint64_t n) {
 }
 
 u256::u256(rational const& n) {
+#if 1
+    for (unsigned i = 0; i < 4; ++i) {
+        m_num[i] = 0;
+        for (unsigned j = 0; j < 64; ++j)
+            m_num[i] |= n.get_bit(i * 64 + j) << j;
+    }
+#else
     uint8_t bytes[32];
     for (unsigned i = 0; i < 32; ++i)
         bytes[i] = 0;
@@ -21,6 +29,7 @@ u256::u256(rational const& n) {
     auto* v = Hacl_Bignum256_new_bn_from_bytes_be(32, bytes);
     std::uninitialized_copy(v, v + 4, m_num);    
     free(v);
+#endif
 }
 
 
@@ -29,6 +38,7 @@ u256::u256(uint64_t const* v) {
 }
 
 u256 u256::operator*(u256 const& other) const {
+    // TBD: maybe just use mpn_manager::mul?
     uint64_t result[8];
     Hacl_Bignum256_mul(const_cast<uint64_t*>(m_num), const_cast<uint64_t*>(other.m_num), result);
     return u256(result);
@@ -63,7 +73,7 @@ u256 u256::operator<<(uint64_t sh) const {
 u256 u256::operator>>(uint64_t sh) const {
     u256 r;
     if (0 == sh || sh >= 256) 
-        ;
+	;
     else if (sh >= 176) 
         r.m_num[0] = m_num[3] >> (sh - 176);
     else if (sh >= 128) {
@@ -132,6 +142,10 @@ u256 u256::mod(u256 const& other) const {
         VERIFY(Hacl_Bignum256_mod(const_cast<uint64_t*>(other.m_num), a, r.m_num));
         return r;
     }
+    
+    // claim:
+    // a mod 2^k*b = ((a >> k) mod b) << k | (a & ((1 << k) - 1))
+    
     unsigned tz = other.trailing_zeros();
     u256 thz = *this >> tz;
     u256 n = other >> tz;
@@ -156,14 +170,8 @@ u256 u256::mul_inverse() const {
     u256 r0(*this);
     u256 r1(-r0); 
     while (!r1.is_zero()) {
-        auto tmp = t1;
-    }
-#if 0
-    numeral t0 = 1, t1 = 0 - 1;
-    numeral r0 = x, r1 = 0 - x;
-    while (r1 != 0) {
-        numeral q = r0 / r1;
-        numeral tmp = t1;
+        u256 q = r0 / r1;
+        u256 tmp = t1;
         t1 = t0 - q * t1;
         t0 = tmp;
         tmp = r1;
@@ -171,10 +179,6 @@ u256 u256::mul_inverse() const {
         r0 = tmp;
     }
     return t0;
-#endif
-    
-    NOT_IMPLEMENTED_YET();
-    return *this;
 }
 
 unsigned u256::trailing_zeros() const {
@@ -237,11 +241,41 @@ bool u256::operator>(uint64_t other) const {
     return 0 != Hacl_Bignum256_lt_mask(_other, const_cast<uint64_t*>(m_num));
 }
 
-
-std::ostream& u256::display(std::ostream& out) const {
+rational u256::to_rational() const {
     rational n;
     for (unsigned i = 0; i < 4; ++i) 
         if (m_num[i] != 0) 
             n += rational(m_num[i], rational::ui64()) * rational::power_of_two(i * 64);
-    return out << n;
+    return n;
+}
+
+std::ostream& u256::display(std::ostream& out) const {
+    return out << to_rational();
+}
+
+// mpn implements the main functionality needed for unsigned fixed-point arithmetic
+// we could use mpn for add/sub/mul as well and maybe punt on Hacl dependency.
+
+u256 u256::operator/(u256 const& other) const {
+    u256 result;
+    mpn_manager m;
+    mpn_digit rem[8];
+    unsigned n1 = 0, n2 = 0;
+    for (unsigned i = 4; i-- > 0; ) {
+        if (m_num[i]) {
+            n1 = 2 * (i + 1);
+            break;
+        }
+    }
+    for (unsigned i = 4; i-- > 0; ) {
+        if (other.m_num[i]) {
+            n2 = 2 * (i + 1);
+            break;
+        }
+    }
+    VERIFY(m.div(reinterpret_cast<mpn_digit const*>(m_num), n1,
+        reinterpret_cast<mpn_digit const*>(other.m_num), n2,
+        reinterpret_cast<mpn_digit*>(result.m_num),
+        rem));
+    return result;
 }
diff --git a/src/math/bigfix/u256.h b/src/math/bigfix/u256.h
index 0f38e362c..be60e6050 100644
--- a/src/math/bigfix/u256.h
+++ b/src/math/bigfix/u256.h
@@ -10,9 +10,11 @@ public:
     u256();
     u256(uint64_t n);
     u256(rational const& n);
+    rational to_rational() const;
     u256 operator*(u256 const& other) const;
     u256 operator+(u256 const& other) const { u256 r = *this; return r += other; }
     u256 operator-(u256 const& other) const { u256 r = *this; return r -= other; }
+    u256 operator/(u256 const& other) const;
     u256 operator-() const { u256 r = *this; return r.uminus(); }
     u256 operator<<(uint64_t sh) const;
     u256 operator>>(uint64_t sh) const;
@@ -28,6 +30,9 @@ public:
     u256& operator+=(u256 const& other);
     u256& operator*=(u256 const& other);
     u256& operator-=(u256 const& other);
+    u256& operator/=(u256 const& other) { *this = *this / other; return *this; }
+    u256& operator>>=(uint64_t sh) { *this = *this >> sh; return *this; }
+    u256& operator<<=(uint64_t sh) { *this = *this << sh; return *this; }
     u256& uminus();  /* unary minus */
 
     // comparisons
diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h
index a9a45fb46..f412d325a 100644
--- a/src/math/dd/dd_pdd.h
+++ b/src/math/dd/dd_pdd.h
@@ -422,6 +422,7 @@ namespace dd {
         std::ostream& display(std::ostream& out) const { return m.display(out, *this); }
         bool operator==(pdd const& other) const { return root == other.root; }
         bool operator!=(pdd const& other) const { return root != other.root; }
+        unsigned hash() const { return root; }
 
         unsigned power_of_2() const { return m.power_of_2(); }
 
diff --git a/src/math/interval/mod_interval.h b/src/math/interval/mod_interval.h
index 56929acc4..899395fd0 100644
--- a/src/math/interval/mod_interval.h
+++ b/src/math/interval/mod_interval.h
@@ -66,6 +66,8 @@ public:
     mod_interval& intersect_ugt(Numeral const& l);
     mod_interval& intersect_fixed(Numeral const& n);
     mod_interval& intersect_diff(Numeral const& n);
+    mod_interval& update_lo(Numeral const& new_lo);
+    mod_interval& update_hi(Numeral const& new_hi);
 
     mod_interval operator&(mod_interval const& other) const;
     mod_interval operator+(mod_interval const& other) const;
diff --git a/src/math/interval/mod_interval_def.h b/src/math/interval/mod_interval_def.h
index da3e1b843..4a02e41b2 100644
--- a/src/math/interval/mod_interval_def.h
+++ b/src/math/interval/mod_interval_def.h
@@ -229,3 +229,23 @@ mod_interval<Numeral>& mod_interval<Numeral>::intersect_diff(Numeral const& a) {
         hi = a;
     return *this;
 }
+
+template<typename Numeral>
+mod_interval<Numeral>& mod_interval<Numeral>::update_lo(Numeral const& new_lo) {
+    SASSERT(lo <= new_lo);
+    if (lo < hi && hi <= new_lo)
+	set_empty();
+    else
+	lo = new_lo;
+    return *this;
+}
+
+template<typename Numeral>
+mod_interval<Numeral>& mod_interval<Numeral>::update_hi(Numeral const& new_hi) {
+    SASSERT(new_hi <= hi);
+    if (new_hi <= lo && lo < hi)
+	set_empty();
+    else
+	hi = new_hi;
+    return *this;
+}
diff --git a/src/math/polysat/boolean.cpp b/src/math/polysat/boolean.cpp
index eaf99995a..1185c881a 100644
--- a/src/math/polysat/boolean.cpp
+++ b/src/math/polysat/boolean.cpp
@@ -7,7 +7,6 @@ Module Name:
 
 Author:
 
-    Nikolaj Bjorner (nbjorner) 2021-03-19
     Jakob Rath 2021-04-6
 
 --*/
@@ -25,10 +24,15 @@ namespace polysat {
             m_reason.push_back(nullptr);
             m_lemma.push_back(nullptr);
             return var;
-        } else {
+        }
+	else {
             sat::bool_var var = m_unused.back();
             m_unused.pop_back();
             SASSERT_EQ(m_level[var], UINT_MAX);
+	    SASSERT_EQ(m_value[2*var], l_undef);
+	    SASSERT_EQ(m_value[2*var+1], l_undef);
+	    SASSERT_EQ(m_reason[var], nullptr);
+	    SASSERT_EQ(m_lemma[var], nullptr);
             return var;
         }
     }
diff --git a/src/math/polysat/boolean.h b/src/math/polysat/boolean.h
index c7e40bd98..5de29693e 100644
--- a/src/math/polysat/boolean.h
+++ b/src/math/polysat/boolean.h
@@ -7,7 +7,6 @@ Module Name:
 
 Author:
 
-    Nikolaj Bjorner (nbjorner) 2021-03-19
     Jakob Rath 2021-04-6
 
 --*/
diff --git a/src/math/polysat/clause_builder.cpp b/src/math/polysat/clause_builder.cpp
index 5bfdfa805..792da5ce3 100644
--- a/src/math/polysat/clause_builder.cpp
+++ b/src/math/polysat/clause_builder.cpp
@@ -7,7 +7,6 @@ Module Name:
 
 Author:
 
-    Nikolaj Bjorner (nbjorner) 2021-03-19
     Jakob Rath 2021-04-6
 
 --*/
diff --git a/src/math/polysat/clause_builder.h b/src/math/polysat/clause_builder.h
index 26be8ffc2..6942e3498 100644
--- a/src/math/polysat/clause_builder.h
+++ b/src/math/polysat/clause_builder.h
@@ -7,9 +7,15 @@ Module Name:
 
 Author:
 
-    Nikolaj Bjorner (nbjorner) 2021-03-19
     Jakob Rath 2021-04-6
 
+Notes:
+
+ Builds a clause from literals and constraints.
+ Takes care to
+ - resolve with unit clauses and accumulate their dependencies,
+ - skip trivial new constraints such as "4 <= 1".
+
 --*/
 #pragma once
 #include "math/polysat/constraint.h"
@@ -17,16 +23,12 @@ Author:
 
 namespace polysat {
 
-    /// Builds a clause from literals and constraints.
-    /// Takes care to
-    /// - resolve with unit clauses and accumulate their dependencies,
-    /// - skip trivial new constraints such as "4 <= 1".
     class clause_builder {
-        solver& m_solver;
-        sat::literal_vector m_literals;
+        solver&               m_solver;
+        sat::literal_vector   m_literals;
         constraint_ref_vector m_new_constraints;
-        p_dependency_ref m_dep;
-        unsigned m_level = 0;
+        p_dependency_ref      m_dep;
+        unsigned              m_level = 0;
 
     public:
         clause_builder(solver& s);
@@ -42,6 +44,7 @@ namespace polysat {
         /// Add a literal to the clause.
         /// Intended to be used for literals representing a constraint that already exists.
         void push_literal(sat::literal lit);
+      
         /// Add a constraint to the clause that does not yet exist in the solver so far.
         void push_new_constraint(constraint_literal c);
     };
diff --git a/src/math/polysat/constraint.h b/src/math/polysat/constraint.h
index 9b4060c66..6ae437671 100644
--- a/src/math/polysat/constraint.h
+++ b/src/math/polysat/constraint.h
@@ -42,7 +42,6 @@ namespace polysat {
         friend class constraint;
 
         bool_var_manager& m_bvars;
-        // poly_dep_manager& m_dm;
 
         // Association to boolean variables
         ptr_vector<constraint>   m_bv2constraint;
@@ -135,6 +134,9 @@ namespace polysat {
             m_manager->m_bvars.del_var(m_bvar);
         }
 
+	    virtual unsigned hash() const = 0;
+	    virtual bool operator==(constraint const& other) const = 0;
+
         bool is_eq() const { return m_kind == ckind_t::eq_t; }
         bool is_ule() const { return m_kind == ckind_t::ule_t; }
         ckind_t kind() const { return m_kind; }
@@ -360,6 +362,7 @@ namespace polysat {
             else
                 SASSERT_EQ(c->blit(), lit);
         }
+      // NSB review: assumes life-time of c extends use in tmp_assign.
         tmp_assign(constraint_ref const& c, sat::literal lit): tmp_assign(c.get(), lit) {}
         void revert() {
             if (m_should_unassign) {
diff --git a/src/math/polysat/eq_constraint.cpp b/src/math/polysat/eq_constraint.cpp
index daf4b0193..d7cf64bd5 100644
--- a/src/math/polysat/eq_constraint.cpp
+++ b/src/math/polysat/eq_constraint.cpp
@@ -202,4 +202,12 @@ namespace polysat {
             return inequality(zero, p(), true, this);
     }
 
+    unsigned eq_constraint::hash() const {
+    	return p().hash();
+    }
+    
+    bool eq_constraint::operator==(constraint const& other) const {
+    	return other.is_eq() && p() == other.to_eq().p();
+    }
+ 
 }
diff --git a/src/math/polysat/eq_constraint.h b/src/math/polysat/eq_constraint.h
index b1790303b..4974c730e 100644
--- a/src/math/polysat/eq_constraint.h
+++ b/src/math/polysat/eq_constraint.h
@@ -33,6 +33,8 @@ namespace polysat {
         void narrow(solver& s) override;
         bool forbidden_interval(solver& s, pvar v, eval_interval& out_interval, constraint_literal& out_neg_cond) override;
         inequality as_inequality() const override;
+        unsigned hash() const override;
+        bool operator==(constraint const& other) const override;
     };
 
 }
diff --git a/src/math/polysat/ule_constraint.cpp b/src/math/polysat/ule_constraint.cpp
index e53bb26f4..a3b0c7930 100644
--- a/src/math/polysat/ule_constraint.cpp
+++ b/src/math/polysat/ule_constraint.cpp
@@ -293,4 +293,13 @@ namespace polysat {
         else 
             return inequality(rhs(), lhs(), true, this);
     }
+
+    unsigned ule_constraint::hash() const {
+    	return mk_mix(lhs().hash(), rhs().hash(), 23);
+    }
+    
+    bool ule_constraint::operator==(constraint const& other) const {
+        return other.is_ule() && lhs() == other.to_ule().lhs() && rhs() == other.to_ule().rhs();
+    }
+
 }
diff --git a/src/math/polysat/ule_constraint.h b/src/math/polysat/ule_constraint.h
index f0e94e5bd..cc148a811 100644
--- a/src/math/polysat/ule_constraint.h
+++ b/src/math/polysat/ule_constraint.h
@@ -39,6 +39,8 @@ namespace polysat {
         void narrow(solver& s) override;
         bool forbidden_interval(solver& s, pvar v, eval_interval& out_interval, constraint_literal& out_neg_cond) override;
         inequality as_inequality() const override;
+        unsigned hash() const override;
+        bool operator==(constraint const& other) const override;
     };
 
 }
diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp
index 40c79975d..8201f7b87 100644
--- a/src/math/polysat/viable.cpp
+++ b/src/math/polysat/viable.cpp
@@ -19,104 +19,16 @@ and narrow the range using the BDDs that are cached.
 
 --*/
 
+
 #include "math/polysat/viable.h"
 #include "math/polysat/solver.h"
-#include "math/interval/mod_interval_def.h"
+#if NEW_VIABLE
+#include "math/polysat/viable_set_def.h"
+#endif
 
 
 namespace polysat {
 
-#if NEW_VIABLE
-
-    dd::find_t viable_set::find_hint(rational const& d, rational& val) const {
-        if (is_empty())
-            return dd::find_t::empty;         
-        if (contains(d))
-            val = d;
-        else 
-            val = lo;
-        if (lo + 1 == hi || hi == 0 && is_max(lo))
-            return dd::find_t::singleton;            
-        return dd::find_t::multiple;
-    }
-
-    bool viable_set::is_max(rational const& a) const {
-        return a + 1 == rational::power_of_two(m_num_bits);
-    }
-    
-    void viable_set::intersect_eq(rational const& a, bool is_positive) {
-        if (is_positive) 
-            intersect_fixed(a);        
-        else 
-            intersect_diff(a);        
-    }
-
-    bool viable_set::intersect_eq(rational const& a, rational const& b, bool is_positive) {
-        if (!a.is_odd()) {
-            std::function<bool(rational const&)> eval = [&](rational const& x) {
-                return is_positive == (mod(a * x + b, p2()) == 0);
-            };       
-            return narrow(eval);            
-        }
-        if (b == 0)
-            intersect_eq(b, is_positive);
-        else {
-            rational a_inv;
-            VERIFY(a.mult_inverse(m_num_bits, a_inv));
-            intersect_eq(mod(a_inv * -b, p2()), is_positive);
-        }
-        return true;
-    }
-
-    bool viable_set::intersect_le(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) {
-        // x <= 0
-        if (a.is_odd() && b == 0 && c == 0 && d == 0)
-            intersect_eq(b, is_positive);
-        else if (a == 1 && b == 0 && c == 0) {
-            // x <= d or x > d
-            if (is_positive)
-                intersect_ule(d);
-            else 
-                intersect_ugt(d);
-        }
-        else if (a == 0 && c == 1 && d == 0) {
-            // x >= b or x < b
-            if (is_positive)
-                intersect_uge(b);
-            else
-                intersect_ult(b);
-        }
-        // TBD: can also handle wrap-around semantics (for signed comparison)
-        else {
-            std::function<bool(rational const&)> eval = [&](rational const& x) {
-                return is_positive == mod(a * x + b, p2()) <= mod(c * x + d, p2());
-            };       
-            return narrow(eval);
-        }
-        return true;
-    }
-
-    rational viable_set::prev(rational const& p) const {
-        if (p > 0)
-            return p - 1;
-        else
-            return rational::power_of_two(m_num_bits) - 1;
-    }
-
-    bool viable_set::narrow(std::function<bool(rational const&)>& eval) {
-        unsigned budget = 10;
-        while (budget > 0 && !is_empty() && !eval(lo)) {
-            --budget;
-            intersect_diff(lo);
-        }
-        while (budget > 0 && !is_empty() && !eval(prev(hi))) {
-            --budget;
-            intersect_diff(prev(hi));
-        }
-        return 0 < budget;
-    }
-
-#endif
 
     viable::viable(solver& s):
         s(s),
@@ -134,15 +46,19 @@ namespace polysat {
 #endif
     }
 
-
     void viable::push_viable(pvar v) {
         s.m_trail.push_back(trail_instr_t::viable_i);
+#if NEW_VIALBLE
+        m_viable_trail.push_back(std::make_pair(v, alloc(viable_set, *m_viable[v])));
+#else
         m_viable_trail.push_back(std::make_pair(v, m_viable[v]));
+#endif
+
     }
 
     void viable::pop_viable() {
         auto p = m_viable_trail.back();
-        m_viable[p.first] = p.second;
+        m_viable.set(p.first, p.second);
         m_viable_trail.pop_back();
     }
 
@@ -151,9 +67,9 @@ namespace polysat {
     void viable::intersect_eq(rational const& a, pvar v, rational const& b, bool is_positive) {
 #if NEW_VIABLE
         push_viable(v);
-        if (!m_viable[v].intersect_eq(a, b, is_positive)) 
+        if (!m_viable[v]->intersect_eq(a, b, is_positive)) 
             intersect_eq_bdd(v, a, b, is_positive);
-        if (m_viable[v].is_empty())
+        if (m_viable[v]->is_empty())
             s.set_conflict(v);
 #else
         
@@ -184,9 +100,9 @@ namespace polysat {
     void viable::intersect_ule(pvar v, rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) {
 #if NEW_VIABLE
         push_viable(v);
-        if (!m_viable[v].intersect_le(a, b, c, d, is_positive)) 
+        if (!m_viable[v]->intersect_le(a, b, c, d, is_positive)) 
             intersect_ule_bdd(v, a, b, c, d, is_positive);
-        if (m_viable[v].is_empty())
+        if (m_viable[v]->is_empty())
             s.set_conflict(v);
 #else
         bddv const& x = var2bits(v).var();
@@ -230,6 +146,8 @@ namespace polysat {
             for (auto* e : m_constraint_cache) 
                 entries.push_back(e);
             std::stable_sort(entries.begin(), entries.end(), [&](cached_constraint* a, cached_constraint* b) { return a->m_activity < b->m_activity; });
+	    for (auto* e : entries)
+		e->m_activity /= 2;
             for (unsigned i = 0; i < max_entries/2; ++i) {
                 m_constraint_cache.remove(entries[i]);
                 dealloc(entries[i]);
@@ -238,12 +156,12 @@ namespace polysat {
     }
 
     void viable::narrow(pvar v, bdd const& is_false) {        
-        rational bound = m_viable[v].lo;
+        rational bound = m_viable[v]->lo;
         if (var2bits(v).sup(is_false, bound)) 
-            m_viable[v].intersect_ugt(bound);
-        bound = m_viable[v].prev(m_viable[v].hi);
+            m_viable[v]->update_lo(m_viable[v]->next(bound));
+        bound = m_viable[v]->prev(m_viable[v]->hi);
         if (var2bits(v).inf(is_false, bound)) 
-            m_viable[v].intersect_ult(bound);
+            m_viable[v]->update_hi(m_viable[v]->prev(bound));	
     }
 
     void viable::intersect_ule_bdd(pvar v, rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) {
@@ -277,7 +195,7 @@ namespace polysat {
 
     bool viable::has_viable(pvar v) {
 #if NEW_VIABLE
-        return !m_viable[v].is_empty();
+        return !m_viable[v]->is_empty();
 #else
         return !m_viable[v].is_false();
 #endif
@@ -285,7 +203,7 @@ namespace polysat {
 
     bool viable::is_viable(pvar v, rational const& val) {
 #if NEW_VIABLE
-        return m_viable[v].contains(val);
+        return m_viable[v]->contains(val);
 #else
         return var2bits(v).contains(m_viable[v], val);
 #endif
@@ -295,8 +213,8 @@ namespace polysat {
 #if NEW_VIABLE
         push_viable(v);
         IF_VERBOSE(10, verbose_stream() << " v" << v << " != " << val << "\n");
-        m_viable[v].intersect_diff(val);
-        if (m_viable[v].is_empty())
+        m_viable[v]->intersect_diff(val);
+        if (m_viable[v]->is_empty())
             s.set_conflict(v);
 #else
         LOG("pvar " << v << " /= " << val);
@@ -317,7 +235,7 @@ namespace polysat {
 
     dd::find_t viable::find_viable(pvar v, rational & val) {
 #if NEW_VIABLE
-        return m_viable[v].find_hint(s.m_value[v], val);
+        return m_viable[v]->find_hint(s.m_value[v], val);
 #else
         return var2bits(v).find_hint(m_viable[v], s.m_value[v], val);
 #endif
diff --git a/src/math/polysat/viable.h b/src/math/polysat/viable.h
index 479379de5..df57892cb 100644
--- a/src/math/polysat/viable.h
+++ b/src/math/polysat/viable.h
@@ -14,7 +14,51 @@ Author:
 Notes:
 
     NEW_VIABLE uses cheaper book-keeping, but is partial.
-    
+   
+
+Alternative to using rational, instead use fixed-width numerals.
+
+
+map from num_bits to template set
+
+class viable_trail_base {
+public:
+   virtual pop(pvar v) = 0;
+   virtual push(pvar v) = 0;
+   static viable_trail_base* mk_trail(unsigned num_bits);
+};
+
+class viable_trail<Numeral> : public viable_trail_base {
+    vector<viable_set<Numeral>> m_viable;
+    vector<viable_set<Numeral>> m_trail;
+public:
+    void pop(pvar v) override {
+         m_viable[v] = m_trail.back();
+         m_trail.pop_back();
+    }
+    void push(pvar v) override {
+         m_trail.push_back(m_viable[v]);
+    }
+};
+
+// from num-bits to viable_trail_base*
+scoped_ptr_vector<viable_trail_base> m_viable_trails;
+
+viable_set_base& to_viable(pvar v) {
+  return (*m_viable_trails[num_bits(v)])[v];
+}
+
+viable_set_base is required to expose functions:
+lo, hi, 
+prev, next           alternative as bit-vectors
+update_lo (a) 
+update_hi (a) 
+intersect_le  (a, b, c, d)
+intersect_diff (a, b)
+intersect_eq (a, b)
+is_empty
+contains
+
 --*/
 #pragma once
 
@@ -24,37 +68,14 @@ Notes:
 
 #include "math/dd/dd_bdd.h"
 #include "math/polysat/types.h"
-#include "math/interval/mod_interval.h"
-
+#if NEW_VIABLE
+#include "math/polysat/viable_set.h"
+#endif
 
 namespace polysat {
 
     class solver;
 
-#if NEW_VIABLE
-    //
-    // replace BDDs by viable sets that emulate affine relations.
-    // viable_set has an interval of feasible values.
-    // it also can use ternary bit-vectors.
-    // or we could also just use a vector of lbool instead of ternary bit-vectors
-    // updating them at individual positions is relatively cheap instead of copying the
-    // vectors every time a range is narrowed.
-    //
-    class viable_set : public mod_interval<rational> {
-        unsigned m_num_bits;
-        rational p2() const { return rational::power_of_two(m_num_bits); }
-        bool is_max(rational const& a) const override;
-        void intersect_eq(rational const& a, bool is_positive);
-        bool narrow(std::function<bool(rational const&)>& eval);
-    public:
-        viable_set(unsigned num_bits): m_num_bits(num_bits) {}
-        ~viable_set() override {}
-        dd::find_t find_hint(rational const& c, rational& val) const;
-        bool intersect_eq(rational const& a, rational const& b, bool is_positive);
-        bool intersect_le(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive);
-        rational prev(rational const& p) const;
-    };
-#endif
 
     class viable {
         typedef dd::bdd bdd;
@@ -86,8 +107,8 @@ namespace polysat {
                 }
             };
         };
-        vector<viable_set>                                       m_viable; 
-        vector<std::pair<pvar, viable_set>>                      m_viable_trail;
+        scoped_ptr_vector<viable_set>                             m_viable; 
+        vector<std::pair<pvar, viable_set*>>                      m_viable_trail;
         hashtable<cached_constraint*, cached_constraint::hash, cached_constraint::eq> m_constraint_cache;
 
         void intersect_ule_bdd(pvar v, rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive);
@@ -119,7 +140,7 @@ namespace polysat {
 
         void push(unsigned num_bits) { 
 #if NEW_VIABLE
-            m_viable.push_back(viable_set(num_bits)); 
+            m_viable.push_back(alloc(viable_set, num_bits)); 
 #else
             m_viable.push_back(m_bdd.mk_true()); 
 #endif
diff --git a/src/math/polysat/viable_set.h b/src/math/polysat/viable_set.h
new file mode 100644
index 000000000..69f513ff3
--- /dev/null
+++ b/src/math/polysat/viable_set.h
@@ -0,0 +1,56 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    set of viable values as wrap-around interval
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2021-03-19
+    Jakob Rath 2021-04-6
+
+Notes:
+
+    
+  replace BDDs by viable sets that emulate affine relations.
+  viable_set has an interval of feasible values.
+  it also can use ternary bit-vectors.
+  or we could also just use a vector of lbool instead of ternary bit-vectors
+  updating them at individual positions is relatively cheap instead of copying the
+  vectors every time a range is narrowed.
+    
+    
+--*/
+#pragma once
+
+
+#include <limits>
+
+#include "math/dd/dd_bdd.h"
+#include "math/polysat/types.h"
+#include "math/interval/mod_interval.h"
+
+
+namespace polysat {
+
+  
+    class viable_set : public mod_interval<rational> {
+        unsigned m_num_bits;
+        rational p2() const { return rational::power_of_two(m_num_bits); }
+        bool is_max(rational const& a) const override;
+        void intersect_eq(rational const& a, bool is_positive);
+        bool narrow(std::function<bool(rational const&)>& eval);
+    public:
+        viable_set(unsigned num_bits): m_num_bits(num_bits) {}
+        ~viable_set() override {}
+        dd::find_t find_hint(rational const& c, rational& val) const;
+        bool intersect_eq(rational const& a, rational const& b, bool is_positive);
+        bool intersect_le(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive);
+        rational prev(rational const& p) const;
+	    rational next(rational const& p) const;
+    };
+
+}
+
+
diff --git a/src/math/polysat/viable_set_def.h b/src/math/polysat/viable_set_def.h
new file mode 100644
index 000000000..6695997ab
--- /dev/null
+++ b/src/math/polysat/viable_set_def.h
@@ -0,0 +1,121 @@
+/*++
+Copyright (c) 2021 Microsoft Corporation
+
+Module Name:
+
+    set of viable values as wrap-around interval
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2021-03-19
+    Jakob Rath 2021-04-6
+
+    
+--*/
+#pragma once
+
+
+#include "math/polysat/viable_set.h"
+#include "math/interval/mod_interval_def.h"
+
+namespace polysat {
+
+      dd::find_t viable_set::find_hint(rational const& d, rational& val) const {
+        if (is_empty())
+            return dd::find_t::empty;         
+        if (contains(d))
+            val = d;
+        else 
+            val = lo;
+        if (lo + 1 == hi || hi == 0 && is_max(lo))
+            return dd::find_t::singleton;            
+        return dd::find_t::multiple;
+    }
+
+    bool viable_set::is_max(rational const& a) const {
+        return a + 1 == rational::power_of_two(m_num_bits);
+    }
+    
+    void viable_set::intersect_eq(rational const& a, bool is_positive) {
+        if (is_positive) 
+            intersect_fixed(a);        
+        else 
+            intersect_diff(a);        
+    }
+
+    bool viable_set::intersect_eq(rational const& a, rational const& b, bool is_positive) {
+        if (!a.is_odd()) {
+            std::function<bool(rational const&)> eval = [&](rational const& x) {
+                return is_positive == (mod(a * x + b, p2()) == 0);
+            };       
+            return narrow(eval);            
+        }
+        if (b == 0)
+            intersect_eq(b, is_positive);
+        else {
+            rational a_inv;
+            VERIFY(a.mult_inverse(m_num_bits, a_inv));
+            intersect_eq(mod(a_inv * -b, p2()), is_positive);
+        }
+        return true;
+    }
+
+    bool viable_set::intersect_le(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) {
+        // x <= 0
+        if (a.is_odd() && b == 0 && c == 0 && d == 0)
+            intersect_eq(b, is_positive);
+        else if (a == 1 && b == 0 && c == 0) {
+            // x <= d or x > d
+            if (is_positive)
+                intersect_ule(d);
+            else 
+                intersect_ugt(d);
+        }
+        else if (a == 0 && c == 1 && d == 0) {
+            // x >= b or x < b
+            if (is_positive)
+                intersect_uge(b);
+            else
+                intersect_ult(b);
+        }
+        // TBD: can also handle wrap-around semantics (for signed comparison)
+        else {
+            std::function<bool(rational const&)> eval = [&](rational const& x) {
+                return is_positive == mod(a * x + b, p2()) <= mod(c * x + d, p2());
+            };       
+            return narrow(eval);
+        }
+        return true;
+    }
+
+    rational viable_set::prev(rational const& p) const {
+        if (p > 0)
+            return p - 1;
+        else
+            return rational::power_of_two(m_num_bits) - 1;
+    }
+
+    rational viable_set::next(rational const& p) const {
+        if (is_max(p))
+	    return rational(0);
+        else
+            return p + 1;
+    }
+
+    bool viable_set::narrow(std::function<bool(rational const&)>& eval) {
+        unsigned budget = 10;
+        while (budget > 0 && !is_empty() && !eval(lo)) {
+            --budget;
+            intersect_diff(lo);
+        }
+        while (budget > 0 && !is_empty() && !eval(prev(hi))) {
+            --budget;
+            intersect_diff(prev(hi));
+        }
+        return 0 < budget;
+    }
+
+
+}
+
+