diff --git a/src/math/polysat/conflict.cpp b/src/math/polysat/conflict.cpp
index 41ba6614e..dd190915c 100644
--- a/src/math/polysat/conflict.cpp
+++ b/src/math/polysat/conflict.cpp
@@ -461,15 +461,9 @@ namespace polysat {
         s.inc_activity(v);
 
         m_vars.remove(v);
-        for (signed_constraint const c : s.m_viable.get_constraints(v))
-            if (insert_or_replace(c))
+        for (sat::literal const lit : s.m_viable.get_propagation_reason(v))
+            if (insert_or_replace(s.lit2cnstr(lit)))
                 return;
-        for (auto const& i : s.m_viable.units(v)) {
-            if (insert_or_replace(s.eq(i.lo(), i.lo_val())))
-                return;
-            if (insert_or_replace(s.eq(i.hi(), i.hi_val())))
-                return;
-        }
         logger().log(inf_resolve_value(s, v));
 
         revert_pvar(v);
@@ -627,12 +621,8 @@ namespace polysat {
 
                 auto const& j = s.m_justification[v];
                 if (j.is_propagation_by_viable()) {
-                    for (signed_constraint c : s.m_viable.get_constraints(v))
-                        enqueue_constraint(c);
-                    for (auto const& i : s.m_viable.units(v)) {
-                        enqueue_constraint(s.eq(i.lo(), i.lo_val()));
-                        enqueue_constraint(s.eq(i.hi(), i.hi_val()));
-                    }
+                    for (sat::literal const lit : s.m_viable.get_propagation_reason(v))
+                        enqueue_lit(lit);
                 }
                 else if (j.is_propagation_by_slicing()) {
                     s.m_slicing.explain_value(v, enqueue_lit, enqueue_var);
diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp
index c10e71a00..b8eddb457 100644
--- a/src/math/polysat/solver.cpp
+++ b/src/math/polysat/solver.cpp
@@ -740,6 +740,10 @@ namespace polysat {
                 m_viable_fallback.pop_constraint();
                 break;
             }
+            case trail_instr_t::viable_propagation_i: {
+                m_viable.pop_propagation_reason();
+                break;
+            }
             case trail_instr_t::assign_i: {
                 auto v = m_search.back().var();
                 LOG_V(20, "Undo assign_i: v" << v);
@@ -866,7 +870,9 @@ namespace polysat {
             // The fallback solver currently does not detect propagations, because we would need to handle justifications differently.
             // However, this case may still occur if during viable::intersect, we run into the refinement budget,
             // but here, we continue refinement and actually succeed until propagation.
-            assign_propagate_by_viable(v, val);
+            // ???
+            UNREACHABLE();
+            SASSERT(m_justification[v].is_propagation_by_viable());
             return;
         case find_t::multiple:
             j = justification::decision(m_level + 1);
@@ -886,19 +892,10 @@ namespace polysat {
     }
 
     void solver::assign_propagate_by_viable(pvar v, rational const& val) {
-        // NOTE:
-        // The propagation v := val might depend on a lower level than the current level (m_level).
-        // This can happen if the constraints that cause the propagation have been bool-propagated at an earlier level,
-        // but appear later in the stack (cf. replay).
-        // The level of v should then also be the earlier level instead of m_level.
         unsigned lvl = 0;
-        for (signed_constraint c : m_viable.get_constraints(v)) {
-            LOG("due to: " << lit_pp(*this, c));
-            VERIFY(m_bvars.is_assigned(c.blit()));
-            lvl = std::max(lvl, m_bvars.level(c.blit()));
-            for (pvar w : c->vars())
-                if (is_assigned(w))  // TODO: question of which variables are relevant. e.g., v1 > v0 implies v1 > 0 without dependency on v0. maybe add a lemma v1 > v0 --> v1 > 0 on the top level to reduce false variable dependencies? instead of handling such special cases separately everywhere.
-                    lvl = std::max(lvl, get_level(w));
+        for (sat::literal const lit : m_viable.get_propagation_reason(v)) {
+            VERIFY(m_bvars.is_assigned(lit));
+            lvl = std::max(lvl, m_bvars.level(lit));
         }
         assign_propagate(v, val, justification::propagation_by_viable(lvl));
     }
@@ -1597,8 +1594,8 @@ namespace polysat {
                 auto const& j = m_justification[v];
                 out << "\t" << assignment_pp(*this, v, get_value(v)) << " @" << j.level() << " ";
                 if (j.is_propagation_by_viable())
-                    for (auto const& c : m_viable.get_constraints(v))
-                        out << c << " ";
+                    for (sat::literal const lit : m_viable.get_propagation_reason(v))
+                        out << lit2cnstr(lit) << " ";
                 if (j.is_propagation_by_slicing()) {
                     out << "by slicing: ";
                     const_cast<slicing&>(m_slicing).explain_value(v,
diff --git a/src/math/polysat/superposition.cpp b/src/math/polysat/superposition.cpp
index f5590bdda..9954872c7 100644
--- a/src/math/polysat/superposition.cpp
+++ b/src/math/polysat/superposition.cpp
@@ -64,8 +64,9 @@ namespace polysat {
     // c2 ... constraint that is currently false
     // Try to replace it with a new false constraint (derived from superposition with a true constraint)
     lbool ex_polynomial_superposition::find_replacement(signed_constraint c2, pvar v, conflict& core) {
-        for (auto c1 : s.m_viable.get_constraints(v)) {
-            if (!c1->contains_var(v))  // side conditions do not contain v; skip them here
+        for (sat::literal lit1 : s.m_viable.get_propagation_reason(v)) {
+            signed_constraint c1 = s.lit2cnstr(lit1);
+            if (!c1->contains_var(v))  // side conditions and interval endpoint constraints do not contain v; skip them here
                 continue;
             if (!c1.is_pos_eq())
                 continue;
diff --git a/src/math/polysat/trail.h b/src/math/polysat/trail.h
index 4aede79a2..ac0907978 100644
--- a/src/math/polysat/trail.h
+++ b/src/math/polysat/trail.h
@@ -16,7 +16,7 @@ Author:
 
 namespace polysat {
 
-    enum class trail_instr_t {
+    enum class trail_instr_t : std::uint8_t {
         qhead_i,
         lemma_qhead_i,
         add_lemma_i,
@@ -26,6 +26,7 @@ namespace polysat {
         viable_add_i,
         viable_rem_i,
         viable_constraint_i,
+        viable_propagation_i,
         assign_i,
         assign_bool_i
     };
diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp
index 00ed7d394..9d93cb0c0 100644
--- a/src/math/polysat/viable.cpp
+++ b/src/math/polysat/viable.cpp
@@ -25,6 +25,8 @@ TODO: plan to fix the FI "pumping":
         - for inequalities, a coefficient 2^k*a means that intervals are periodic because the upper k bits of x are irrelevant;
           storing the interval for x[K-k:0] would take care of this.
 
+TODO: we can now support propagations by fallback solver
+
 --*/
 
 
@@ -59,12 +61,14 @@ namespace polysat {
         m_units.push_back({});
         m_equal_lin.push_back(nullptr);
         m_diseq_lin.push_back(nullptr);
+        m_propagation_reasons.push_back({});
     }
 
     void viable::pop_var() {
         m_units.pop_back();
         m_equal_lin.pop_back();
         m_diseq_lin.pop_back();
+        m_propagation_reasons.pop_back();
     }
 
     viable::entry* viable::alloc_entry(pvar var) {
@@ -176,7 +180,7 @@ namespace polysat {
             rational val;
             switch (find_viable(v, val)) {
             case find_t::singleton:
-                propagate(v, val);
+                SASSERT(s.m_justification[v].is_propagation_by_viable());
                 prop = true;
                 break;
             case find_t::empty:
@@ -194,7 +198,16 @@ namespace polysat {
         return prop;
     }
 
-    void viable::propagate(pvar v, rational const& val) {
+    void viable::propagate(pvar v, rational const& val, ptr_vector<entry> const& reason) {
+        SASSERT(m_propagation_reasons[v].begin == 0);
+        SASSERT(m_propagation_reasons[v].len == 0);
+        unsigned const begin = m_propagation_reason_storage.size();
+
+        auto add_reason = [this](signed_constraint c) {
+            s.try_assign_eval(c);
+            m_propagation_reason_storage.push_back(c.blit());
+        };
+
         // NOTE: all propagations must be justified by a prefix of \Gamma,
         //       otherwise dependencies may be missed during conflict resolution.
         //       The propagation reason for v := val consists of the following constraints:
@@ -202,31 +215,36 @@ namespace polysat {
         //       - side conditions
         //       - i.lo() == i.lo_val() for each unit interval i
         //       - i.hi() == i.hi_val() for each unit interval i
-
-        // NSB review:
-        // the bounds added by x < p and p < x in forbidden_intervals
-        // match_non_max, match_non_zero
-        // use values that are approximations. Then the propagations in
-        // try_assign_eval are incorrect.
-        // For example, x > p means x has forbidden interval [0, p + 1[,
-        // the numeric interval is [0, 1[, but p + 1 == 1 is not ensured
-        // even p may have free variables.
-        // the proper side condition on p + 1 is -1 > p or -2 >= p or p + 1 != 0
-        // I am disabling match_non_max and match_non_zero from forbidden_interval
-        // The narrowing rules in ule_constraint already handle the bounds propagaitons
-        // as it propagates p != -1 and 0 != q (p < -1, q > 0),
-        //
-
-        for (auto const& c : get_constraints(v)) {
-            s.try_assign_eval(c);
-        }
-        for (auto const& i : units(v)) {
-            s.try_assign_eval(s.eq(i.lo(), i.lo_val()));
-            s.try_assign_eval(s.eq(i.hi(), i.hi_val()));
+        for (entry* e : reason) {
+            for (signed_constraint c : e->side_cond)
+                add_reason(c);
+            for (signed_constraint c : e->src)
+                add_reason(c);
+            auto const& i = e->interval;
+            add_reason(s.eq(i.lo(), i.lo_val()));
+            add_reason(s.eq(i.hi(), i.hi_val()));
         }
+
+        unsigned const len = m_propagation_reason_storage.size() - begin;
+        SASSERT(begin + len == m_propagation_reason_storage.size());
+        m_propagation_reasons[v].begin = begin;
+        m_propagation_reasons[v].len = len;
+        m_propagation_trail.push_back(v);
+        s.m_trail.push_back(trail_instr_t::viable_propagation_i);
+
         s.assign_propagate_by_viable(v, val);
     }
 
+    void viable::pop_propagation_reason() {
+        pvar const v = m_propagation_trail.back();
+        m_propagation_trail.pop_back();
+        unsigned begin = m_propagation_reasons[v].begin;
+        unsigned len = m_propagation_reasons[v].len;
+        SASSERT(len > 0);
+        SASSERT(begin + len == m_propagation_reason_storage.size());
+        m_propagation_reason_storage.shrink(begin);
+    }
+
     bool viable::intersect(pvar v, signed_constraint const& c) {
         LOG("intersect v" << v << " in " << lit_pp(s, c));
         if (s.is_assigned(v)) {
@@ -1161,14 +1179,19 @@ namespace polysat {
 
     find_t viable::find_viable(pvar v, rational& lo) {
         rational hi;
-        switch (find_viable2(v, lo, hi)) {
+        ptr_vector<entry> relevant_entries;
+        switch (find_viable2(v, lo, hi, relevant_entries)) {
         case l_true:
             if (hi < 0) {
                 // fallback solver, treat propagations as decisions for now
                 // (this is because the propagation justification currently always uses intervals, which is unsound in this case)
                 return find_t::multiple;
             }
-            return (lo == hi) ? find_t::singleton : find_t::multiple;
+            if (lo == hi) {
+                propagate(v, lo, relevant_entries);
+                return find_t::singleton;
+            }
+            return find_t::multiple;
         case l_false:
             return find_t::empty;
         default:
@@ -1177,6 +1200,7 @@ namespace polysat {
     }
 
     bool viable::has_upper_bound(pvar v, rational& out_hi, vector<signed_constraint>& out_c) {
+        NOT_IMPLEMENTED_YET();
         entry const* first = m_units[v].get_entries(s.size(v));  // TODO: take other sizes into account
         entry const* e = first;
         bool found = false;
@@ -1213,6 +1237,7 @@ namespace polysat {
     }
 
     bool viable::has_lower_bound(pvar v, rational& out_lo, vector<signed_constraint>& out_c) {
+        NOT_IMPLEMENTED_YET();
         entry const* first = m_units[v].get_entries(s.size(v));  // TODO: take other sizes into account
         entry const* e = first;
         bool found = false;
@@ -1249,6 +1274,7 @@ namespace polysat {
     }
 
     bool viable::has_max_forbidden(pvar v, signed_constraint const& c, rational& out_lo, rational& out_hi, vector<signed_constraint>& out_c) {
+        NOT_IMPLEMENTED_YET();
         // TODO:
         // - skip intervals adjacent to c's interval if they contain side conditions on y?
         //      constraints over y are allowed if level(c) < level(y) (e.g., boolean propagated)
@@ -1347,8 +1373,6 @@ namespace polysat {
         return true;
     }
 
-
-
     // When iterating over intervals:
     // - instead of only intervals of v, go over intervals of each entry of overlaps
     // - need a function to map interval from overlap into an interval over v
@@ -1438,7 +1462,7 @@ namespace polysat {
     // - how to integrate fallback solver?
     //   when lowest level fails, we can try more refinement there.
     //   in case of refinement loop, try fallback solver with constraints only from lower level.
-    lbool viable::find_viable2(pvar v, rational& lo, rational& hi) {
+    lbool viable::find_viable2(pvar v, rational& lo, rational& hi, ptr_vector<entry>& relevant_entries) {
         fixed_bits_info fbi;
 
         if (!collect_bit_information(v, true, fbi))
@@ -1472,7 +1496,7 @@ namespace polysat {
 
         rational const& max_value = s.var2pdd(v).max_value();
 
-        lbool result_lo = find_on_layers(v, widths, overlaps, fbi, rational::zero(), max_value, lo);
+        lbool result_lo = find_on_layers(v, widths, overlaps, fbi, rational::zero(), max_value, lo, relevant_entries);
         if (result_lo == l_false)
             return l_false;  // conflict
         if (result_lo == l_undef)
@@ -1483,7 +1507,7 @@ namespace polysat {
             return l_true;
         }
 
-        lbool result_hi = find_on_layers(v, widths, overlaps, fbi, lo + 1, max_value, hi);
+        lbool result_hi = find_on_layers(v, widths, overlaps, fbi, lo + 1, max_value, hi, relevant_entries);
         if (result_hi == l_false)
             hi = lo;  // no other viable value
         if (result_hi == l_undef)
@@ -1526,10 +1550,10 @@ namespace polysat {
         fixed_bits_info const& fbi,
         rational const& to_cover_lo,
         rational const& to_cover_hi,
-        rational& val
+        rational& val,
+        ptr_vector<entry>& relevant_entries
     ) {
         ptr_vector<entry> refine_todo;
-        ptr_vector<entry> relevant_entries;
 
         // max number of interval refinements before falling back to the univariate solver
         unsigned const refinement_budget = 100;
@@ -1938,7 +1962,12 @@ namespace polysat {
             entry** intervals_begin = intervals.data() + first_interval;
             unsigned num_intervals = intervals.size() - first_interval - 1;
 
-            if (!set_conflict_by_interval_rec(v, w, intervals_begin, num_intervals, core, create_lemma, lemma, vars_to_explain))
+            bool const ok = set_conflict_by_interval_rec(v, w, intervals_begin, num_intervals, core, create_lemma, lemma, vars_to_explain);
+
+            SASSERT(!intervals.back());
+            // intervals.pop_back();  // get rid of the dummy space (we don't use the intervals anymore after setting the conflict, so this is not necessary)
+
+            if (!ok)
                 return false;
         }
 
diff --git a/src/math/polysat/viable.h b/src/math/polysat/viable.h
index 753cf3d4a..1a48ae128 100644
--- a/src/math/polysat/viable.h
+++ b/src/math/polysat/viable.h
@@ -93,6 +93,14 @@ namespace polysat {
         ptr_vector<entry>       m_diseq_lin;    // entries that have distinct non-zero multipliers
         svector<std::tuple<pvar, entry_kind, entry*>>     m_trail; // undo stack
 
+        struct reason_ptr {
+            unsigned begin = 0;
+            unsigned len = 0;
+        };
+        sat::literal_vector     m_propagation_reason_storage;
+        pvar_vector             m_propagation_trail;
+        svector<reason_ptr>     m_propagation_reasons;
+
         unsigned size(pvar v) const;
 
         bool well_formed(entry* e);
@@ -185,7 +193,7 @@ namespace polysat {
 
         void insert(entry* e, pvar v, ptr_vector<entry>& entries, entry_kind k);
 
-        void propagate(pvar v, rational const& val);
+        void propagate(pvar v, rational const& val, ptr_vector<entry> const& reason);
 
         /**
          * Enter conflict state when intervals cover the full domain.
@@ -214,7 +222,7 @@ namespace polysat {
          * NOTE: out_hi is set to -1 by the fallback solver.
          * \return l_true on success, l_false on conflict, l_undef on resource limit
          */
-        lbool find_viable2(pvar v, rational& out_lo, rational& out_hi);
+        lbool find_viable2(pvar v, rational& out_lo, rational& out_hi, ptr_vector<entry>& out_relevant_entries);
 
         lbool find_on_layers(
             pvar v,
@@ -223,7 +231,8 @@ namespace polysat {
             fixed_bits_info const& fbi,
             rational const& to_cover_lo,
             rational const& to_cover_hi,
-            rational& out_val);
+            rational& out_val,
+            ptr_vector<entry>& out_relevant_entries);
 
         lbool find_on_layer(
             pvar v,
@@ -234,8 +243,8 @@ namespace polysat {
             rational const& to_cover_lo,
             rational const& to_cover_hi,
             rational& out_val,
-            ptr_vector<entry>& refine_todo,
-            ptr_vector<entry>& relevant_entries);
+            ptr_vector<entry>& out_refine_todo,
+            ptr_vector<entry>& out_relevant_entries);
 
         std::pair<entry*, bool> find_value(rational const& val, entry* entries);
 
@@ -252,6 +261,8 @@ namespace polysat {
         void pop_viable();
         void push_viable();
 
+        void pop_propagation_reason();
+
         /**
          * update state of viable for pvar v
          * based on affine constraints.
@@ -297,6 +308,7 @@ namespace polysat {
 
         /**
          * Find a next viable value for variable.
+         * Side effect: propagate v in the solver if only a single value is viable.
          */
         find_t find_viable(pvar v, rational& out_val);
 
@@ -310,90 +322,24 @@ namespace polysat {
 
         std::ostream& display(std::ostream& out, char const* delimiter = "") const;
 
-        class iterator {
-            entry* curr = nullptr;
-            bool   visited = false;
-            unsigned idx = 0;
+        class propagation_reason {
+            sat::literal const* m_begin;
+            sat::literal const* m_end;
+            propagation_reason(sat::literal const* begin, sat::literal const* end): m_begin(begin), m_end(end) {}
+            friend class viable;
         public:
-            iterator(entry* curr, bool visited) :
-                curr(curr), visited(visited || !curr) {}
-
-            iterator& operator++() {
-                if (idx < curr->side_cond.size() + curr->src.size() - 1)
-                    ++idx;
-                else {
-                    idx = 0;
-                    visited = true;
-                    curr = curr->next();
-                }
-                return *this;
-            }
-
-            signed_constraint& operator*() {
-                return idx < curr->side_cond.size() ? curr->side_cond[idx] : curr->src[idx - curr->side_cond.size()];
-            }
-
-            bool operator==(iterator const& other) const {
-                return visited == other.visited && curr == other.curr;
-            }
-
-            bool operator!=(iterator const& other) const {
-                return !(*this == other);
-            }
+            sat::literal const* begin() const { return m_begin; }
+            sat::literal const* end() const { return m_end; }
         };
 
-        class constraints {
-            viable const& v;
-            pvar var;
-        public:
-            constraints(viable const& v, pvar var) : v(v), var(var) {}
-            // TODO: take other widths into account!
-            iterator begin() const { return iterator(v.m_units[var].get_entries(v.size(var)), false); }
-            iterator end() const { return iterator(v.m_units[var].get_entries(v.size(var)), true); }
-        };
-
-        constraints get_constraints(pvar v) const {
-            return constraints(*this, v);
+        propagation_reason get_propagation_reason(pvar v) const {
+            auto const& r = m_propagation_reasons[v];
+            SASSERT(r.len > 0);
+            sat::literal const* begin = &m_propagation_reason_storage[r.begin];
+            sat::literal const* end = begin + r.len;
+            return { begin, end };
         }
 
-        class int_iterator {
-            entry* curr = nullptr;
-            bool visited = false;
-        public:
-            int_iterator(entry* curr, bool visited) :
-                curr(curr), visited(visited || !curr) {}
-            int_iterator& operator++() {
-                visited = true;
-                curr = curr->next();
-                return *this;
-            }
-
-            eval_interval const& operator*() {
-                return curr->interval;
-            }
-
-            bool operator==(int_iterator const& other) const {
-                return visited == other.visited && curr == other.curr;
-            }
-
-            bool operator!=(int_iterator const& other) const {
-                return !(*this == other);
-            }
-
-        };
-
-        class intervals {
-            viable const& v;
-            pvar var;
-        public:
-            intervals(viable const& v, pvar var): v(v), var(var) {}
-            // TODO: take other widths into account!
-            int_iterator begin() const { return int_iterator(v.m_units[var].get_entries(v.size(var)), false); }
-            int_iterator end() const { return int_iterator(v.m_units[var].get_entries(v.size(var)), true); }
-        };
-
-        intervals units(pvar v) { return intervals(*this, v); }
-
         std::ostream& display(std::ostream& out, pvar v, char const* delimiter = "") const;
 
         struct var_pp {