diff --git a/src/math/polysat/CMakeLists.txt b/src/math/polysat/CMakeLists.txt
index 7b63994b9..2986cb94e 100644
--- a/src/math/polysat/CMakeLists.txt
+++ b/src/math/polysat/CMakeLists.txt
@@ -12,7 +12,6 @@ z3_add_component(polysat
     forbidden_intervals.cpp
     inference_logger.cpp
     justification.cpp
-    linear_solver.cpp
     log.cpp
     naming.cpp
     op_constraint.cpp
@@ -33,7 +32,6 @@ z3_add_component(polysat
     viable.cpp
     viable_fallback.cpp
   COMPONENT_DEPENDENCIES
-    bigfix
     dd
     euf
     interval
diff --git a/src/math/polysat/clause.h b/src/math/polysat/clause.h
index c1c3170bb..b951be24d 100644
--- a/src/math/polysat/clause.h
+++ b/src/math/polysat/clause.h
@@ -24,7 +24,7 @@ namespace polysat {
     // NB code review:
     // right, ref-count is unlikely the right mechanism.
     // In the SAT solver all clauses are managed in one arena (auxiliarary and redundant)
-    // and deleted when they exist the arena.
+    // and deleted when they exit the arena.
     //
     class clause {
     public:
diff --git a/src/math/polysat/fixplex.h b/src/math/polysat/fixplex.h
deleted file mode 100644
index bd68dfb7b..000000000
--- a/src/math/polysat/fixplex.h
+++ /dev/null
@@ -1,483 +0,0 @@
-/*++
-Copyright (c) 2014 Microsoft Corporation
-
-Module Name:
-
-    fixplex.h
-
-Abstract:
-
-    Fixed-precision unsigned integer simplex tableau.
-
-Author:
-
-    Nikolaj Bjorner (nbjorner) 2021-03-19
-    Jakob Rath 2021-04-6
-
---*/
-
-#pragma once
-
-#include <limits>
-#include "math/simplex/sparse_matrix.h"
-#include "fixplex_mod_interval.h"
-#include "util/heap.h"
-#include "util/map.h"
-#include "util/rational.h"
-#include "util/lbool.h"
-#include "util/uint_set.h"
-#include "util/dependency.h"
-#include "util/ref.h"
-#include "util/params.h"
-#include "util/union_find.h"
-
-inline rational to_rational(uint64_t n) { return rational(n, rational::ui64()); }
-inline unsigned trailing_zeros(unsigned short n) { return trailing_zeros((uint32_t)n); }
-inline unsigned trailing_zeros(unsigned char n) { return trailing_zeros((uint32_t)n); }
-inline unsigned numeral2hash(unsigned char const& n) { return n; }
-inline unsigned numeral2hash(unsigned short const& n) { return n; }
-inline unsigned numeral2hash(uint32_t const& n) { return n; }
-inline unsigned numeral2hash(uint64_t const& n) { return static_cast<unsigned>(n ^ (n >> 32)); }
-
-
-namespace polysat {
-
-    typedef unsigned var_t;
-    
-    struct fixplex_base {
-        virtual ~fixplex_base() {}
-        virtual lbool make_feasible() = 0;
-        virtual void add_row(var_t base, unsigned num_vars, var_t const* vars, rational const* coeffs) = 0;
-        virtual void del_row(var_t base_var) = 0;
-        virtual void push() = 0;
-        virtual void pop(unsigned n) = 0;
-        virtual std::ostream& display(std::ostream& out) const = 0;
-        virtual void collect_statistics(::statistics & st) const = 0;
-        virtual void set_bounds(var_t v, rational const& lo, rational const& hi, unsigned dep) = 0;        
-        virtual void set_value(var_t v, rational const& val, unsigned dep) = 0;
-        virtual rational get_value(var_t v) = 0;
-        virtual void restore_bound() = 0;   
-        virtual void add_le(var_t v, var_t w, unsigned dep) = 0;
-        virtual void add_lt(var_t v, var_t w, unsigned dep) = 0;
-        virtual void restore_ineq() = 0;
-        virtual bool inconsistent() const = 0;
-        virtual unsigned_vector const& get_unsat_core() const = 0;
-        virtual void updt_params(params_ref const& p) = 0;
-
-    };
-
-    struct ineq {
-        var_t v = UINT_MAX;
-        var_t w = UINT_MAX;
-        bool strict = false;
-        u_dependency* dep = nullptr;
-        ineq(var_t v, var_t w, u_dependency* dep, bool s) :
-            v(v), w(w), strict(s), dep(dep) {}
-
-        std::ostream& display(std::ostream& out) const {
-            return out << "v" << v << (strict ? " < v" : " <= v") << w;
-        }
-    };
-
-    template<typename Ext>
-    class fixplex : public fixplex_base {
-    public:
-        typedef typename Ext::numeral numeral;
-        typedef typename Ext::scoped_numeral scoped_numeral;
-        typedef typename Ext::manager manager;
-        typedef simplex::sparse_matrix<Ext> matrix;
-        typedef typename matrix::row row;
-        typedef typename matrix::row_iterator row_iterator;
-        typedef typename matrix::col_iterator col_iterator;
-
-        struct var_eq {
-            var_t x, y;
-            u_dependency* dep;
-            var_eq(var_t x, var_t y, u_dependency* dep):
-                x(x), y(y), dep(dep) {}
-        };
-
-    protected:
-        struct var_lt {
-            bool operator()(var_t v1, var_t v2) const { return v1 < v2; }
-        };
-        typedef heap<var_lt> var_heap;
-
-        struct stats {
-            unsigned m_num_pivots;
-            unsigned m_num_infeasible;
-            unsigned m_num_checks;
-            unsigned m_num_approx;
-            stats() { reset(); }
-            void reset() {
-                memset(this, 0, sizeof(*this));
-            }
-        };
-
-        enum pivot_strategy_t {
-            S_BLAND,
-            S_GREATEST_ERROR,
-            S_LEAST_ERROR,
-            S_DEFAULT
-        };
-
-        struct var_info : public mod_interval<numeral> {
-            unsigned      m_base2row:29;
-            unsigned      m_is_base:1;
-            numeral       m_value = 0;
-            u_dependency* m_lo_dep = nullptr;
-            u_dependency* m_hi_dep = nullptr;
-            var_info():
-                m_base2row(0),
-                m_is_base(false)
-            {}
-            ~var_info() override {}
-            var_info& operator&=(mod_interval<numeral> const& range) {
-                mod_interval<numeral>::operator=(range & *this);
-                return *this;
-            }
-            var_info& operator=(mod_interval<numeral> const& range) {
-                mod_interval<numeral>::operator=(range);
-                return *this;                
-            }
-        };
-
-        struct row_info {
-            bool    m_integral { true };
-            var_t   m_base;
-            numeral m_value;
-            numeral m_base_coeff;            
-        };
-
-        struct stashed_bound : var_info {
-            var_t m_var;
-            stashed_bound(var_t v, var_info const& i):
-                var_info(i),
-                m_var(v)
-            {}
-        };
-
-        struct fix_entry {
-            var_t x = null_var;
-            u_dependency* dep = nullptr;
-            fix_entry(var_t x, u_dependency* dep): x(x), dep(dep) {}
-            fix_entry() {}
-        };
-
-        enum trail_i {
-            inc_level_i,
-            set_bound_i,
-            set_inconsistent_i,
-            add_ineq_i,
-            add_row_i,
-            add_eq_i,
-            fixed_val_i
-        };
-
-        static const var_t null_var = UINT_MAX;
-        reslimit&                   m_limit;
-        mutable manager             m;
-        mutable matrix              M;
-        unsigned                    m_max_iterations = UINT_MAX;
-        unsigned                    m_num_non_integral = 0;
-        uint_set                    m_non_integral;
-        var_heap                    m_to_patch;
-        vector<var_info>            m_vars;
-        vector<row_info>            m_rows;
-
-        bool                        m_bland = false ;
-        unsigned                    m_blands_rule_threshold = 1000;
-        unsigned                    m_num_repeated = 0;
-        random_gen                  m_random;
-        uint_set                    m_left_basis;
-        unsigned_vector             m_unsat_core; 
-        bool                        m_inconsistent = false;
-        unsigned_vector             m_base_vars;
-        stats                       m_stats;
-        vector<stashed_bound>       m_stashed_bounds;
-        u_dependency_manager        m_deps;
-        svector<trail_i>            m_trail;
-        svector<var_t>              m_row_trail;
-
-        // euqality propagation
-        union_find_default_ctx      m_union_find_ctx;
-        union_find<>                m_union_find;
-        vector<var_eq>              m_var_eqs;
-        vector<numeral>             m_fixed_vals;
-        map<numeral, fix_entry, typename manager::hash, typename manager::eq> m_value2fixed_var;
-        uint_set                    m_eq_rows;
-
-        // inequalities
-        svector<ineq>               m_ineqs;
-        uint_set                    m_ineqs_to_propagate;
-        uint_set                    m_touched_vars;
-        vector<unsigned_vector>     m_var2ineqs;
-
-        // bound propagation
-        uint_set                    m_bound_rows;
-
-    public:
-        fixplex(params_ref const& p, reslimit& lim):
-            m_limit(lim),
-            M(m),
-            m_to_patch(1024),
-            m_union_find(m_union_find_ctx) {
-            updt_params(p);
-        }
-
-        ~fixplex() override;
-
-        void push() override;
-        void pop(unsigned n) override;
-        bool inconsistent() const override { return m_inconsistent; }
-        void updt_params(params_ref const& p) override;
-
-        lbool make_feasible() override;
-        void add_row(var_t base, unsigned num_vars, var_t const* vars, rational const* coeffs) override;
-        std::ostream& display(std::ostream& out) const override;
-        void collect_statistics(::statistics & st) const override;
-        void del_row(var_t base_var) override;
-        void set_bounds(var_t v, rational const& lo, rational const& hi, unsigned dep) override;
-        void set_value(var_t v, rational const& val, unsigned dep) override;
-        rational get_value(var_t v) override;
-        void restore_bound() override;
-        void add_le(var_t v, var_t w, unsigned dep) override { add_ineq(v, w, dep, false); }
-        void add_lt(var_t v, var_t w, unsigned dep) override { add_ineq(v, w, dep, true); }
-        virtual void restore_ineq() override;
-
-        void set_bounds(var_t v, numeral const& lo, numeral const& hi, unsigned dep);
-        void update_bounds(var_t v, numeral const& l, numeral const& h, u_dependency* dep);
-        void unset_bounds(var_t v) { m_vars[v].set_free(); }
-
-
-        numeral const& lo(var_t var) const { return m_vars[var].lo; }
-        numeral const& hi(var_t var) const { return m_vars[var].hi; }
-        numeral const& value(var_t var) const { return m_vars[var].m_value; }
-        void set_max_iterations(unsigned n) { m_max_iterations = n; }
-        unsigned get_num_vars() const { return m_vars.size(); }
-        void reset();
-
-        vector<var_eq> const& var_eqs() const { return m_var_eqs; }
-
-        void add_row(var_t base, unsigned num_vars, var_t const* vars, numeral const* coeffs);
-        
-        unsigned_vector const& get_unsat_core() const override { return m_unsat_core; }
-
-    private:
-
-        svector<std::pair<unsigned, unsigned>> stack;
-        uint_set on_stack;
-        lbool propagate_ineqs(unsigned idx);
-
-        std::ostream& display_row(std::ostream& out, row const& r, bool values = true) const;
-        var_t get_base_var(row const& r) const { return m_rows[r.id()].m_base; }
-
-        void update_value_core(var_t v, numeral const& delta);
-        void ensure_var(var_t v);
-
-        bool patch();
-        bool propagate_ineqs();
-        bool propagate_row_eqs();
-        bool propagate_row_bounds();
-        bool is_satisfied();
-
-        struct backoff {
-            unsigned m_tries = 0;
-            unsigned m_delay = 1;
-            bool should_propagate() {
-                return ++m_tries >= m_delay;
-            }
-            void update(bool progress) {
-                m_tries = 0;
-                if (progress)
-                    m_delay = 1;
-                else
-                    ++m_delay;
-            }
-        };
-
-        backoff m_propagate_eqs_backoff;
-        backoff m_propagate_bounds_backoff;
-
-
-
-        var_t select_smallest_var() { return m_to_patch.empty()?null_var:m_to_patch.erase_min(); }
-        lbool make_var_feasible(var_t x_i);
-        bool is_infeasible_row(var_t x);
-        bool is_parity_infeasible_row(var_t x);
-        bool is_offset_row(row const& r, numeral& cx, var_t& x, numeral& cy, var_t & y) const;
-        void lookahead_eq(row const& r1, numeral const& cx, var_t x, numeral const& cy, var_t y);
-        void get_offset_eqs(row const& r);
-        void fixed_var_eh(u_dependency* dep, var_t x);
-        var_t find(var_t x) { return m_union_find.find(x); }
-        void merge(var_t x, var_t y) { m_union_find.merge(x, y); }
-        void eq_eh(var_t x, var_t y, u_dependency* dep);
-        bool propagate_row(row const& r);
-        bool propagate_ineq(ineq const& i);
-        bool propagate_strict_bounds(ineq const& i);
-        bool propagate_non_strict_bounds(ineq const& i);
-        bool new_bound(row const& r, var_t x, mod_interval<numeral> const& range);
-        bool new_bound(ineq const& i, var_t x, numeral const& lo, numeral const& hi, u_dependency* a = nullptr, u_dependency* b = nullptr, u_dependency* c = nullptr, u_dependency* d = nullptr);
-        void conflict(ineq const& i, u_dependency* a = nullptr, u_dependency* b = nullptr, u_dependency* c = nullptr, u_dependency* d = nullptr);
-        void conflict(u_dependency* a);
-        void conflict(u_dependency* a, u_dependency* b, u_dependency* c = nullptr, u_dependency* d = nullptr) { conflict(m_deps.mk_join(m_deps.mk_join(a, b), m_deps.mk_join(c, d))); }
-        u_dependency* row2dep(row const& r);
-        void pivot(var_t x_i, var_t x_j, numeral const& b, numeral const& value);
-        numeral value2delta(var_t v, numeral const& new_value) const;
-        numeral value2error(var_t v, numeral const& new_value) const;
-        void update_value(var_t v, numeral const& delta);
-        bool can_pivot(var_t x_i, numeral const& new_value, numeral const& a_ij, var_t x_j);
-        bool has_minimal_trailing_zeros(var_t y, numeral const& b);
-        var_t select_pivot(var_t x_i, numeral const& new_value, numeral& out_b);
-        var_t select_pivot_core(var_t x, numeral const& new_value, numeral& out_b);
-        bool in_bounds(var_t v) const { return in_bounds(v, value(v)); }
-        bool in_bounds(var_t v, numeral const& b) const { return in_bounds(b, m_vars[v]); }
-        bool in_bounds(numeral const& val, mod_interval<numeral> const& range) const { return range.contains(val); }
-        bool is_free(var_t v) const { return lo(v) == hi(v); }
-        bool is_non_free(var_t v) const { return !is_free(v); }
-        bool is_fixed(var_t v) const { return lo(v) + 1 == hi(v); }
-        bool is_valid_variable(var_t v) const { return v < m_vars.size(); }
-        bool is_base(var_t x) const { return m_vars[x].m_is_base; }
-        row base2row(var_t x) const { return row(m_vars[x].m_base2row); }
-        numeral const& row2value(row const& r) const { return m_rows[r.id()].m_value; }
-        numeral const& row2base_coeff(row const& r) const { return m_rows[r.id()].m_base_coeff; }
-        var_t row2base(row const& r) const { return m_rows[r.id()].m_base; }
-        bool row_is_integral(row const& r) const { return m_rows[r.id()].m_integral; }
-        void set_base_value(var_t x); 
-        numeral solve_for(numeral const& row_value, numeral const& coeff);
-        int get_num_non_free_dep_vars(var_t x_j, int best_so_far);
-        void add_patch(var_t v);
-        var_t select_var_to_fix();
-        void check_blands_rule(var_t v);
-        pivot_strategy_t pivot_strategy() { return m_bland ? S_BLAND : S_DEFAULT; }
-        var_t select_error_var(bool least);
-        void set_infeasible_base(var_t v);
-        void set_infeasible_bounds(var_t v);
-
-        u_dependency* mk_leaf(unsigned dep) { return UINT_MAX == dep ? nullptr : m_deps.mk_leaf(dep); }
-
-        // facilities for handling inequalities
-        void add_ineq(var_t v, var_t w, unsigned dep, bool strict);
-        void touch_var(var_t x);
-
-
-        bool is_solved(row const& r) const;
-        bool is_solved(var_t v) const { SASSERT(is_base(v)); return is_solved(base2row(v)); }
-
-        bool well_formed() const;                 
-        bool well_formed_row(row const& r) const;
-
-        void  del_row(row const& r);
-
-        var_t select_pivot_blands(var_t x, numeral const& new_value, numeral& out_b);
-        bool can_improve(var_t x, numeral const& new_value, var_t y, numeral const& b);
-
-        bool pivot_base_vars();
-        bool elim_base(var_t v);
-
-        bool eliminate_var(
-            row const& r_y,
-            col_iterator const& z_col,
-            unsigned tz_b,
-            numeral const& old_value_y);
-    };
-
-    template<typename uint_type>
-    struct generic_uint_ext {
-        typedef uint_type numeral;
-
-        struct manager {
-            typedef uint_type numeral;
-            struct hash {
-                unsigned operator()(numeral const& n) const { 
-                    return numeral2hash(n); 
-                }
-            };
-            struct eq {
-                bool operator()(numeral const& a, numeral const& b) const {
-                    return a == b;
-                }
-            };
-            numeral from_rational(rational const& n) { return static_cast<uint_type>(n.get_uint64()); }
-            rational to_rational(numeral const& n) const { return ::to_rational(n); }
-            void reset() {}
-            void reset(numeral& n) { n = 0; }
-            void del(numeral const& n) {}
-            bool is_zero(numeral const& n) const { return n == 0; }
-            bool is_one(numeral const& n) const { return n == 1; }
-            bool is_even(numeral const& n) const { return (n & 1) == 0; }
-            bool is_minus_one(numeral const& n) const { return n + 1 == 0; }
-            void add(numeral const& a, numeral const& b, numeral& r) { r = a + b; }
-            void sub(numeral const& a, numeral const& b, numeral& r) { r = a - b; }
-            void mul(numeral const& a, numeral const& b, numeral& r) { r = a * b; }
-            void set(numeral& r, numeral const& a) { r = a; }
-            void neg(numeral& a) { a = 0 - a; }
-            numeral inv(numeral const& a) { return 0 - a; }
-            void swap(numeral& a, numeral& b) { std::swap(a, b); }
-            unsigned trailing_zeros(numeral const& a) const { return ::trailing_zeros(a); }
-            numeral mul_inverse(numeral const& x) {
-                if (x == 0)
-                    return 0;
-                numeral t0 = 1, t1 = 0 - 1;
-                numeral r0 = x, r1 = 0 - x;
-                while (r1 != 0) {
-                    numeral q = r0 / r1;
-                    numeral tmp = t1;
-                    t1 = t0 - q * t1;
-                    t0 = tmp;
-                    tmp = r1;
-                    r1 = r0 - q * r1;
-                    r0 = tmp;
-                }
-                return t0;
-            }
-            numeral gcd(numeral x, numeral y) {
-                if (x == 0) 
-                    return y;
-                if (y == 0)
-                    return x;
-                unsigned tz = trailing_zeros(x);
-                numeral shift = std::min(trailing_zeros(y), tz);
-                x >>= tz;
-                if (x == 1) 
-                    return x << shift;
-                if (y == 1) 
-                    return y << shift;
-                if (x == y) 
-                    return x << shift;
-                do {
-                    tz = trailing_zeros(y);
-                    y >>= tz;
-                    if (x > y) 
-                        swap(x, y);
-                    y -= x;
-                }
-                while (y != 0);
-                return x << shift;
-            }
-
-            std::ostream& display(std::ostream& out, numeral const& x) const { 
-                return out << pp(x); 
-            }
-        };
-        typedef _scoped_numeral<manager> scoped_numeral;
-
-    };
-
-    typedef generic_uint_ext<uint64_t> uint64_ext;
-
-
-    template<typename Ext>
-    inline std::ostream& operator<<(std::ostream& out, fixplex<Ext> const& fp) {
-        return fp.display(out);
-    }
-
-
-    inline std::ostream& operator<<(std::ostream& out, ineq const& i) {
-        return i.display(out);
-    }
- 
-
-
-};
-
diff --git a/src/math/polysat/fixplex_def.h b/src/math/polysat/fixplex_def.h
deleted file mode 100644
index db27d8a0d..000000000
--- a/src/math/polysat/fixplex_def.h
+++ /dev/null
@@ -1,1538 +0,0 @@
-/*++
-Copyright (c) 2021 Microsoft Corporation
-
-Module Name:
-
-    fixplex_def.h
-
-Abstract:
-
-    Fixed-precision unsigned integer simplex tableau.
-
-Author:
-
-    Nikolaj Bjorner (nbjorner) 2021-03-19
-    Jakob Rath 2021-04-6
-
-Notes:
-
-Equality pivoting.
-==================
-
-Similar to normal pivoting except base variable must have minimal power of 2
-to ensure that pivoting preserves solutions (Olm-Seidl condition).
-
-Assigning values to base variables could be revised.
-It is desirable to entirely avoid computing values for base variables.
-The requirement is really to establish that there exists a solution within bounds.
-
-
-
-Inequality handling.
-====================
-
-- try patch:
-  x <= y, value(x) > value(y):
-   - x is non-basic: value(x) := value(y); update values of basic.
-   - y is non-basic: value(y) := value(x); update values of basic.
-   - x (y) is basic: pivot, update
-       
-- conflict and bounds:
-  x <= y, lo(x) > hi(y)  
-  x < y, lo(x) >= hi(y)  
-  Conflict detection depends on effectiveness of bounds propagation
-
-  Test case:
-    x <= y, y <= z, z < x should result in a conflict without branching.
-
-- branch (and bound):
-  x <= y, value(x) > value(y):
-     Let delta = (value(x) + value(y)) / 2 (computed as (value(x) - value(y)) / 2 + value(y))
-     case split:
-        x <= delta or x > delta
-     
-     Case x <= delta blocks current solution.
-     Case x > delta incurs bounds propagation on y, y > delta, that also blocks current solution.
-
-- cuts:
-  would be good to understand how to adapt a notion of cuts for modular case.       
-
-
-TODOs:
-- complete inequality propagation using incremental topological sort + patching using assignment to minimal values
-- equality extraction
-- bounds + row propagaiton
-- freedom interval patching
-- cuts
-- branch
-- update synthesis of bounds tightnening
-
---*/
-
-#pragma once
-
-#include "math/polysat/fixplex.h"
-#include "math/simplex/sparse_matrix_def.h"
-#include "fixplex_mod_interval_def.h"
-#include <iostream>
-
-namespace polysat {
-
-    template<typename Ext>
-    fixplex<Ext>::~fixplex() {
-        reset();
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::updt_params(params_ref const& p) {
-        m_max_iterations = p.get_uint("max_iterations", m_max_iterations);
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::push() {
-        m_trail.push_back(trail_i::inc_level_i);
-        m_deps.push_scope();
-        m_union_find_ctx.get_trail_stack().push_scope();
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::pop(unsigned n) {
-        m_deps.pop_scope(n);
-        m_union_find_ctx.get_trail_stack().pop_scope(n);
-        while (n > 0) {
-            switch (m_trail.back()) {
-            case trail_i::inc_level_i:
-                --n;
-                break;
-            case trail_i::set_bound_i:
-                restore_bound();
-                break;
-            case trail_i::add_row_i:
-                del_row(m_row_trail.back());
-                m_row_trail.pop_back();
-                break;
-            case trail_i::add_ineq_i:
-                restore_ineq();
-                break;
-            case trail_i::set_inconsistent_i:
-                SASSERT(m_inconsistent);
-                m_inconsistent = false;
-                m_unsat_core.reset();
-                break;
-            case trail_i::add_eq_i:
-                m_var_eqs.pop_back();
-                break;
-            case trail_i::fixed_val_i:
-                m_value2fixed_var.erase(m_fixed_vals.back());
-                m_fixed_vals.pop_back();
-                break;
-            default:
-                UNREACHABLE();
-            }
-            m_trail.pop_back();
-        }
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::ensure_var(var_t v) {
-        while (v >= m_vars.size()) {
-            M.ensure_var(m_vars.size());
-            m_vars.push_back(var_info());
-            m_union_find.mk_var();
-            m_var2ineqs.push_back(unsigned_vector());
-        }
-        if (m_to_patch.get_bounds() <= v)
-            m_to_patch.set_bounds(2 * v + 1);
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::reset() {
-        M.reset();
-        m_to_patch.reset();
-        m_vars.reset();
-        m_rows.reset();
-        m_left_basis.reset();
-        m_base_vars.reset();
-    }
-
-    // TBD: where does parity test go?
-
-    template<typename Ext>
-    lbool fixplex<Ext>::make_feasible() {
-        ++m_stats.m_num_checks;
-        m_left_basis.reset();
-        m_num_repeated = 0;
-        m_bland = false;
-        unsigned num_iterations = 0;
-        SASSERT(well_formed());
-        while (m_limit.inc() && num_iterations < m_max_iterations) {
-            if (inconsistent())
-                return l_false;
-            if (!propagate_ineqs())
-                return l_false;
-            if (is_satisfied())
-                return l_true;
-            if (!propagate_row_bounds())
-                return l_false;
-            if (!propagate_row_eqs())
-                return l_false;
-            if (!patch())                 
-                return l_undef;
-            
-            ++num_iterations;            
-            SASSERT(well_formed());
-        }
-        return l_undef;
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::patch() {
-        var_t v = select_var_to_fix();
-        if (v == null_var)
-            return false;
-        check_blands_rule(v);
-        switch (make_var_feasible(v)) {
-        case l_true:
-            return true;
-        case l_false:
-            m_to_patch.insert(v);
-            set_infeasible_base(v);
-            return true;
-        case l_undef:
-            m_to_patch.insert(v);
-            return false;
-        }
-        UNREACHABLE();
-        return true;
-    }   
-
-    template<typename Ext>
-    void fixplex<Ext>::add_row(var_t base_var, unsigned num_vars, var_t const* vars, rational const* coeffs) {
-        vector<numeral> _coeffs;
-        for (unsigned i = 0; i < num_vars; ++i)
-            _coeffs.push_back(m.from_rational(coeffs[i]));
-        add_row(base_var, num_vars, vars, _coeffs.data());
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::add_row(var_t base_var, unsigned num_vars, var_t const* vars, numeral const* coeffs) {
-        for (unsigned i = 0; i < num_vars; ++i)
-            ensure_var(vars[i]);
-
-        m_base_vars.reset();
-        row r = M.mk_row();
-        for (unsigned i = 0; i < num_vars; ++i)
-            if (coeffs[i] != 0)
-                M.add_var(r, coeffs[i], vars[i]);
-
-        numeral base_coeff = 0;
-        numeral value = 0;
-        for (auto const& e : M.get_row(r)) {
-            var_t v = e.var();
-            if (v == base_var)
-                base_coeff = e.coeff();
-            else {
-                if (is_base(v))
-                    m_base_vars.push_back(v);
-                value += e.coeff() * m_vars[v].m_value;
-            }
-        }
-        SASSERT(base_coeff != 0);
-        SASSERT(!is_base(base_var));
-        while (m_rows.size() <= r.id())
-            m_rows.push_back(row_info());
-        m_rows[r.id()].m_base = base_var;
-        m_rows[r.id()].m_base_coeff = base_coeff;
-        m_rows[r.id()].m_value = value;
-        m_vars[base_var].m_base2row = r.id();
-        m_vars[base_var].m_is_base = true;
-        set_base_value(base_var);
-        add_patch(base_var);
-        if (!pivot_base_vars())
-            ++m_stats.m_num_approx;
-        SASSERT(well_formed_row(r));
-        SASSERT(well_formed());
-        m_trail.push_back(trail_i::add_row_i);
-        m_row_trail.push_back(base_var);
-        m_eq_rows.insert(r.id());
-        m_bound_rows.insert(r.id());
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::pivot_base_vars() {
-        bool ok = true;
-        for (auto v : m_base_vars)
-            if (!elim_base(v))
-                ok = false;
-        m_base_vars.reset();
-        return ok;
-    }
-
-    /**
-     * Eliminate base variable v from all rows except where v is basic.
-     * Return false if elimination required to multiply a non-basic row with an even number.
-     * This happens when the parity in the non-basic row is smaller than the parity of v in
-     * the basic row. It is expected to be a corner case and we don't try to solve this
-     * inside of fixplex. Instead to deal with the corner case we assume the layer around
-     * fixplex uses a solution from fixplex as a starting point for a complete search (in polysat).
-     */
-    template<typename Ext>
-    bool fixplex<Ext>::elim_base(var_t v) {
-        SASSERT(is_base(v));
-        row r = base2row(v);
-        numeral b = row2base_coeff(r);
-        unsigned tz_b = m.trailing_zeros(b);
-        auto col_it = M.col_begin(v);
-        auto const col_end = M.col_end(v);
-        for (; col_it != col_end; ++col_it) {
-            if (r.id() == col_it.get_row().id())
-                continue;
-            numeral value_v = value(v);
-            if (!eliminate_var(r, col_it, tz_b, value_v))
-                return false;
-        }
-        return true;
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::del_row(row const& r) {
-        var_t var = row2base(r);
-        m_vars[var].m_is_base = false;
-        m_vars[var].set_free();
-        m_rows[r.id()].m_base = null_var;
-        m_non_integral.remove(r.id());
-        M.del(r);
-        m_eq_rows.remove(r.id());
-        m_bound_rows.remove(r.id());
-        SASSERT(M.col_begin(var) == M.col_end(var));
-        SASSERT(well_formed());
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::del_row(var_t var) {
-        TRACE("polysat", tout << var << "\n";);
-        row r;
-        if (is_base(var)) {
-            r = base2row(var);
-        }
-        else {
-            unsigned tz = UINT_MAX;
-            numeral coeff;
-            for (auto [r2, r2_entry] : M.get_col(var)) {
-                unsigned tzc = m.trailing_zeros(r2_entry->coeff());
-                if (tzc < tz) {
-                    r = r2;
-                    tz = tzc;
-                    coeff = r2_entry->coeff();
-                    if (tz == 0)
-                        break;
-                }
-            }
-            if (tz == UINT_MAX)
-                return;
-            var_t old_base = row2base(r);
-            numeral new_value;
-            var_info& vi = m_vars[old_base];
-            if (!vi.contains(value(old_base)))
-                new_value = vi.lo;
-            else
-                new_value = value(old_base);
-            // need to move var such that old_base comes in bound.
-            pivot(old_base, var, coeff, new_value);
-            SASSERT(is_base(var));
-            SASSERT(base2row(var).id() == r.id());
-            SASSERT(vi.contains(value(old_base)));
-        }
-        del_row(r);
-        TRACE("polysat", display(tout););
-        SASSERT(well_formed());
-    }
-
-    /**
-     * increment v by delta
-     */
-    template<typename Ext>
-    void fixplex<Ext>::update_value(var_t v, numeral const& delta) {
-        if (delta == 0)
-            return;
-        m_vars[v].m_value += delta;
-        touch_var(v);
-        SASSERT(!is_base(v));
-        //
-        // v <- v + delta
-        // s*s_coeff + R = 0, where R contains v*v_coeff 
-        // -> 
-        // R.value += delta*v_coeff
-        // s.value = - R.value / s_coeff
-        // 
-        for (auto [r, r_entry] : M.get_col(v)) {
-            row_info& ri = m_rows[r.id()];
-            var_t s = ri.m_base;
-            ri.m_value += delta * r_entry->coeff();
-            set_base_value(s);
-            add_patch(s);
-        }
-    }
-
-
-    /**
-     * Attempt to improve assigment to make x feasible.
-     *
-     * return l_false if x is base variable of infeasible row
-     * return l_true if it is possible to find an assignment that improves
-     * return l_undef if the row could not be used for an improvement
-     *
-     * delta - improvement over previous gap to feasible bound.
-     * new_value - the new value to assign to x within its bounds.
-     */
-
-    template<typename Ext>
-    lbool fixplex<Ext>::make_var_feasible(var_t x) {
-        if (in_bounds(x))
-            return l_true;
-        if (m_vars[x].is_empty())
-            return l_false;
-        numeral new_value = m_vars[x].closest_value(value(x));
-        numeral b;
-        var_t y = select_pivot_core(x, new_value, b);
-        if (y != null_var) 
-            return pivot(x, y, b, new_value), l_true;
-        else if (is_infeasible_row(x))
-            return l_false;
-        else
-            return l_undef;        
-    }
-
-    template<typename Ext>
-    var_t fixplex<Ext>::select_pivot(var_t x, numeral const& new_value, numeral& out_b) {
-        if (m_bland)
-            return select_pivot_blands(x, new_value, out_b);
-        return select_pivot_core(x, new_value, out_b);
-    }
-
-    /**
-       \brief Select a variable y in the row r defining the base var x,
-       s.t. y can be used to patch the error in x_i.  Return null_var
-       if there is no y. Otherwise, return y and store its coefficient
-       in out_b.
-
-       The routine gives up if the coefficients of all free variables do not have the minimal
-       number of trailing zeros.
-    */
-    template<typename Ext>
-    var_t fixplex<Ext>::select_pivot_core(var_t x, numeral const& new_value, numeral& out_b) {
-        SASSERT(is_base(x));
-        var_t max = get_num_vars();
-        var_t result = max;
-        row r = base2row(x);
-        int n = 0;
-        unsigned best_col_sz = UINT_MAX;
-        int best_so_far = INT_MAX;
-        numeral a = row2base_coeff(r);
-        numeral row_value = row2value(r) + a * new_value;
-        numeral delta_y = 0;
-        numeral delta_best = 0;
-        bool best_in_bounds = false;
-
-        for (auto const& r : M.get_row(r)) {
-            var_t y = r.var();
-            numeral const& b = r.coeff();
-            if (x == y)
-                continue;
-            if (!has_minimal_trailing_zeros(y, b))
-                continue;
-            numeral new_y_value = solve_for(row_value - b * value(y), b);
-            bool _in_bounds = in_bounds(y, new_y_value);
-            if (!_in_bounds) {
-                if (lo(y) - new_y_value < new_y_value - hi(y))
-                    delta_y = new_y_value - lo(y);
-                else
-                    delta_y = new_y_value - hi(y) - 1;
-            }
-            int num = get_num_non_free_dep_vars(y, best_so_far);
-            unsigned col_sz = M.column_size(y);
-            bool is_improvement = false, is_plateau = false;
-
-            // improvement criteria would need some scrutiny.
-            if (best_so_far == INT_MAX)
-                is_improvement = true;
-            else if (!best_in_bounds && _in_bounds)
-                is_improvement = true;
-            else if (!best_in_bounds && !_in_bounds && delta_y < delta_best)
-                is_improvement = true;
-            else if (best_in_bounds && _in_bounds && num < best_so_far)
-                is_improvement = true;
-            else if (best_in_bounds && _in_bounds && num == best_so_far && col_sz < best_col_sz)
-                is_improvement = true;
-            else if (!best_in_bounds && !_in_bounds && delta_y == delta_best && best_so_far == num && col_sz == best_col_sz)
-                is_plateau = true;
-            else if (best_in_bounds && _in_bounds && best_so_far == num && col_sz == best_col_sz)
-                is_plateau = true;
-
-            if (is_improvement) {
-                result = y;
-                out_b = b;
-                best_so_far = num;
-                best_col_sz = col_sz;
-                best_in_bounds = _in_bounds;
-                delta_best = delta_y;
-                n = 1;
-            }
-            else if (is_plateau) {
-                n++;
-                if (m_random() % n == 0) {
-                    result = y;
-                    out_b = b;
-                }
-            }
-        }
-        if (result == max)
-            return null_var;
-        if (!best_in_bounds && delta_best >= value2delta(x, new_value))
-            return null_var;
-        return result;
-    }
-
-    template<typename Ext>
-    var_t fixplex<Ext>::select_pivot_blands(var_t x, numeral const& new_value, numeral& out_b) {
-        SASSERT(is_base(x));
-        unsigned max = get_num_vars();
-        var_t result = max;
-        row r = base2row(x);
-        for (auto const& c : M.get_col(r)) {
-            var_t y = c.var();
-            if (x == y || y >= result)
-                continue;
-            numeral const& b = c.coeff();
-            if (can_improve(y, b)) {
-                out_b = b;
-                result = y;
-            }
-        }
-        return result < max ? result : null_var;
-    }
-
-    /**
-     * determine whether setting x := new_value
-     * allows to change the value of y in a direction
-     * that reduces or maintains the overall error.
-     */
-    template<typename Ext>
-    bool fixplex<Ext>::can_improve(var_t x, numeral const& new_x_value, var_t y, numeral const& b) {
-        row r = base2row(x);
-        numeral row_value = row2value(r) + row2base_coeff(r) * new_x_value;
-        numeral new_y_value = solve_for(row_value - b * value(y), b);
-        if (in_bounds(y, new_y_value))
-            return true;
-        return value2error(y, new_y_value) <= value2error(x, value(x));
-    }
-
-    /**
-     * Compute delta to add to the value, such that value + delta is either lo(v), or hi(v) - 1
-     * A pre-condition, when used during pivoting, is that value is not in the interval [lo(v),hi(v)[,
-     * and therefore as a consequence lo(v) != hi(v).
-     */
-    template<typename Ext>
-    typename fixplex<Ext>::numeral
-        fixplex<Ext>::value2delta(var_t v, numeral const& val) const {
-        if (lo(v) - val < val - hi(v))
-            return lo(v) - val;
-        else
-            return hi(v) - val - 1;
-    }
-
-    template<typename Ext>
-    typename fixplex<Ext>::numeral fixplex<Ext>::value2error(var_t v, numeral const& value) const {
-        if (in_bounds(v))
-            return 0;
-        SASSERT(lo(v) != hi(v));
-        if (lo(v) - value < value - hi(v))
-            return lo(v) - value;
-        else
-            return value - hi(v) - 1;
-    }
-
-    /**
-     * The the bounds of variable v.
-     * If the current value of v, value(v), is in bounds, no further updates are made.
-     * If value(v) is outside the the new bounds, then
-     * - the tableau is updated if v is non-basic.
-     * - the variable v is queued to patch if v is basic.
-     */
-    template<typename Ext>
-    void fixplex<Ext>::set_bounds(var_t v, numeral const& l, numeral const& h, unsigned dep) {
-        ensure_var(v);
-        update_bounds(v, l, h, mk_leaf(dep));
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::update_bounds(var_t v, numeral const& l, numeral const& h, u_dependency* dep) {
-        if (inconsistent())
-            return;
-        auto lo0 = m_vars[v].lo;
-        auto hi0 = m_vars[v].hi;
-        m_stashed_bounds.push_back(stashed_bound(v, m_vars[v]));
-        m_trail.push_back(trail_i::set_bound_i);
-       //         std::cout << "new bound " << v << " " << m_vars[v] << " " << mod_interval<numeral>(l, h) << " -> ";        
-        m_vars[v] &= mod_interval(l, h);
-        if (lo0 != m_vars[v].lo)
-            m_vars[v].m_lo_dep = dep;
-        if (hi0 != m_vars[v].hi)
-            m_vars[v].m_hi_dep = dep;
-        // std::cout << m_vars[v] << "\n";
-        if (m_vars[v].is_empty()) {
-            conflict(m_vars[v].m_lo_dep, m_vars[v].m_hi_dep);
-            return;
-        }
-        if (in_bounds(v))
-            return;
-
-        SASSERT(lo(v) != hi(v));
-        if (is_base(v))
-            add_patch(v);
-        else 
-            update_value(v, value2delta(v, value(v)));
-        SASSERT(is_base(v) || in_bounds(v));
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::set_bounds(var_t v, rational const& _lo, rational const& _hi, unsigned dep) {
-        numeral lo = m.from_rational(_lo);
-        numeral hi = m.from_rational(_hi);
-        set_bounds(v, lo, hi, dep);
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::set_value(var_t v, rational const& _val, unsigned dep) {
-        numeral val = m.from_rational(_val);
-        set_bounds(v, val, val + 1, dep);
-    }
-
-    template<typename Ext>
-    rational fixplex<Ext>::get_value(var_t v) {
-        return m.to_rational(m_vars[v].m_value);
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::restore_bound() {
-        auto& b = m_stashed_bounds.back();
-        m_vars[b.m_var].lo = b.lo;
-        m_vars[b.m_var].hi = b.hi;
-        m_vars[b.m_var].m_lo_dep = b.m_lo_dep;
-        m_vars[b.m_var].m_hi_dep = b.m_hi_dep;
-        m_stashed_bounds.pop_back();
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::add_ineq(var_t v, var_t w, unsigned dep, bool strict) {
-        if (inconsistent())
-            return;
-        if (v == w) {
-            if (strict)
-                conflict(mk_leaf(dep));
-            return;
-        }
-        ensure_var(v);
-        ensure_var(w);
-        unsigned idx = m_ineqs.size();        
-        m_var2ineqs[v].push_back(idx);
-        m_var2ineqs[w].push_back(idx);
-        touch_var(v);
-        m_ineqs_to_propagate.insert(idx);
-        m_trail.push_back(trail_i::add_ineq_i);
-        m_ineqs.push_back(ineq(v, w, mk_leaf(dep), strict));
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::restore_ineq() {
-        auto const& ineq = m_ineqs.back();
-        m_var2ineqs[ineq.v].pop_back();
-        m_var2ineqs[ineq.w].pop_back();
-        m_ineqs.pop_back();
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::touch_var(var_t v) {
-        m_touched_vars.insert(v);       
-    }
-
-    /**
-    * Check if the current assignment satisfies the inequalities
-    */
-    template<typename Ext>
-    bool fixplex<Ext>::is_satisfied() {
-        if (!m_to_patch.empty())
-            return false;
-        if (!m_non_integral.empty())
-            return false;
-        for (auto var : m_touched_vars) {
-            for (auto idx : m_var2ineqs[var]) {
-                if (idx >= m_ineqs.size())
-                    continue;
-                auto& ineq = m_ineqs[idx];
-                var_t v = ineq.v;
-                var_t w = ineq.w;
-                if (ineq.strict && value(v) >= value(w))
-                    return false;
-                if (!ineq.strict && value(v) > value(w))
-                    return false;
-            }
-        }
-        m_touched_vars.reset();
-        return true;
-    }
-
-    /**
-    * Propagate bounds and check if the current inequalities are satisfied
-    */
-
-    template<typename Ext>
-    bool fixplex<Ext>::propagate_ineqs() {
-        lbool r = l_true;
-        while (!m_ineqs_to_propagate.empty() && r == l_true) {
-            unsigned idx = *m_ineqs_to_propagate.begin();  
-            if (idx >= m_ineqs.size()) {
-                m_ineqs_to_propagate.remove(idx);
-                continue;
-            }
-            r = propagate_ineqs(idx);
-            if (r == l_undef)
-                return true;
-            m_ineqs_to_propagate.remove(idx);
-        }        
-        return r == l_true;
-    }
-
-    /**
-     * Check if the coefficient b of y has the minimal number of trailing zeros.
-     * In other words, the coefficient b is a multiple of the smallest power of 2.
-     */
-    template<typename Ext>
-    bool fixplex<Ext>::has_minimal_trailing_zeros(var_t y, numeral const& b) {
-        unsigned tz1 = m.trailing_zeros(b);
-        if (tz1 == 0)
-            return true;
-        for (auto [_, r_entry] : M.get_col(y)) {
-            numeral c = r_entry->coeff();
-            unsigned tz2 = m.trailing_zeros(c);
-            if (tz1 > tz2)
-                return false;
-        }
-        return true;
-    }
-
-    /**
-     * Determine if row is linear infeasiable.
-     * A row is linear infeasible if it can be established
-     * that none of the available assignments within current
-     * bounds let the row add up to 0.
-     *
-     * Assume the row is of the form:
-     *   ax + by + cz = 0
-     * with bounds
-     *   x : [lo_x, hi_x[, y : [lo_y, hi_y[, z : [lo_z : hi_z[
-     *
-     * Let range = [lo_x, hi_x[ + [lo_y, hi_y[ + [lo_z : hi_z[
-     * Claim. If range does not contain 0, then the row is infeasible.
-     *
-     */
-    template<typename Ext>
-    bool fixplex<Ext>::is_infeasible_row(var_t x) {
-        SASSERT(is_base(x));
-        auto r = base2row(x);
-        mod_interval<numeral> range(0, 1);
-        for (auto const& e : M.get_row(r)) {
-            var_t v = e.var();
-            numeral const& c = e.coeff();
-            range += m_vars[v] * c;
-            if (range.is_free())
-                return false;
-        }
-        return !range.contains(0);
-    }
-
-    /**
-     * Check if row is infeasible modulo parity constraints.
-     * Let parity be the minimal power of two of coefficients to non-fixed variables.
-     * Let fixed be the sum of fixed variables.
-     * A row is infeasible if parity > the smallest power of two dividing fixed.
-     *
-     * Other parity tests are possible.
-     * The "range" parity test checks if the minimal parities of all but one variable are outside
-     * the range of the value range of a selected variable.
-     */
-    template<typename Ext>
-    bool fixplex<Ext>::is_parity_infeasible_row(var_t x) {
-        SASSERT(is_base(x));
-        auto r = base2row(x);
-        if (row_is_integral(row(r)))
-            return false;
-        numeral fixed = 0;
-        unsigned parity = UINT_MAX;
-        for (auto const& e : M.get_row(row(r))) {
-            var_t v = e.var();
-            auto c = e.coeff();
-            if (is_fixed(v))
-                fixed += value(v) * c;
-            else
-                parity = std::min(m.trailing_zeros(c), parity);
-        }
-
-        if (m.trailing_zeros(fixed) < parity)
-            return true;
-
-        return false;
-    }
-
-
-    /**
-       \brief Given row
-
-         r_x = a*x + b*y + rest = 0
-
-       Assume:
-
-         base(r_x) = x
-         value(r_x) = value(b*y + rest)
-         old_value(y) := value(y)
-
-       Effect:
-
-         base(r_x)  := y
-         value(x)   := new_value
-         value(r_x) := value(r_x) - b*value(y) + a*new_value
-         value(y)   := -value(r_x) / b
-         base_coeff(r_x) := b
-
-       Let r be a row where y has coefficient c != 0.
-       Assume: trailing_zeros(c) >= trailing_zeros(b)
-
-         z = base(r)
-         d = base_coeff(r)
-         b1 = (b >> tz(b))
-         c1 = (c >> tz(b))
-         r <- b1 * r  - c1 * r_x
-         value(r) := b1 * value(r) - b1 * old_value(y) - c1 * value(r_x)
-         value(z) := - value(r) / d
-         base_coeff(r) := b1 * base_coeff(r)
-    */
-    template<typename Ext>
-    void fixplex<Ext>::pivot(var_t x, var_t y, numeral const& b, numeral const& new_value) {
-        ++m_stats.m_num_pivots;
-        SASSERT(is_base(x));
-        SASSERT(!is_base(y));
-        var_info& xI = m_vars[x];
-        var_info& yI = m_vars[y];
-        unsigned rx = xI.m_base2row;
-        auto& row_x = m_rows[rx];
-        row r_x(rx);
-        numeral const& a = row_x.m_base_coeff;
-        numeral old_value_y = yI.m_value;
-        TRACE("fixplex", display_row(tout << "pivot " << x << " " << y << "\n", r_x););
-        row_x.m_base = y;
-        row_x.m_value = row_x.m_value - b * old_value_y + a * new_value;
-        row_x.m_base_coeff = b;
-        yI.m_base2row = rx;
-        yI.m_is_base = true;
-        set_base_value(y);
-        xI.m_is_base = false;
-        xI.m_value = new_value;
-        touch_var(x);
-        add_patch(y);
-        SASSERT(well_formed_row(r_x));
-
-        unsigned tz_b = m.trailing_zeros(b);
-
-        auto col_it = M.col_begin(y);
-        auto const col_end = M.col_end(y);
-        for (; col_it != col_end; ++col_it) {
-            auto r_z = col_it.get_row();
-            unsigned rz = r_z.id();
-            if (rz == rx)
-                continue;
-            TRACE("fixplex", display_row(tout << "eliminate ", r_z, false) << "\n";);
-            VERIFY(eliminate_var(r_x, col_it, tz_b, old_value_y));
-            TRACE("fixplex", display_row(tout << "eliminated ", r_z, false) << "\n";);
-            add_patch(row2base(r_z));
-        }
-        SASSERT(well_formed());
-    }
-
-    /**
-     * r_y         - row where y is base variable
-     * z_col       - column iterator that contains row where z is base variable.
-     * tz_b        - number of trailing zeros to coefficient of y in r_y
-     * old_value_y - the value of y used to compute row2value(r_z)
-     *
-     * Implied variables:
-     * r_z         - row that contains y with z base variable, z != y
-     * c           - coefficient of y in r_z
-     * 
-     * returns true if elimination preserves equivalence (is lossless).
-     * 
-     * TBD: add r_z.id() to m_eq_rows, m_bound_rows with some frequency?
-     * 
-     */
-    template<typename Ext>
-    bool fixplex<Ext>::eliminate_var(
-        row const& r_y,
-        col_iterator const& z_col,
-        unsigned tz_b,
-        numeral const& old_value_y) {
-
-        row const& r_z = z_col.get_row();
-        numeral c = z_col.get_row_entry().coeff();
-        numeral b = row2base_coeff(r_y);
-        auto z = row2base(r_z);        
-        unsigned tz_c = m.trailing_zeros(c);
-        unsigned tz = std::min(tz_b, tz_c);
-        numeral b1 = b >> tz;
-        numeral c1 = 0 - (c >> tz);
-        M.mul(r_z, b1);
-        M.add(r_z, c1, r_y);
-        auto& row_z = m_rows[r_z.id()];
-        row_z.m_value = (b1 * (row2value(r_z) - c * old_value_y)) + c1 * row2value(r_y);
-        row_z.m_base_coeff *= b1;
-        set_base_value(z);
-        m_bound_rows.insert(r_z.id());
-        m_eq_rows.insert(r_z.id());
-        SASSERT(well_formed_row(r_z));
-        return tz_b <= tz_c;
-    }
-
-    /***
-    * Record an infeasible row.
-    */
-    template<typename Ext>
-    void fixplex<Ext>::set_infeasible_base(var_t v) {
-        SASSERT(is_base(v));
-        auto row = base2row(v);
-        ptr_vector<u_dependency> todo;
-        for (auto const& e : M.get_row(row)) {
-            var_t v = e.var();
-            todo.push_back(m_vars[v].m_lo_dep);
-            todo.push_back(m_vars[v].m_hi_dep);
-        }
-        SASSERT(!m_inconsistent);
-        SASSERT(m_unsat_core.empty());
-        m_inconsistent = true;
-        m_trail.push_back(trail_i::set_inconsistent_i);
-        m_deps.linearize(todo, m_unsat_core);
-        ++m_stats.m_num_infeasible;
-    }
-
-    /**
-       \brief Return the number of base variables that are non free and are v dependent.
-       The function adds 1 to the result if v is non free.
-       The function returns with a partial result r if r > best_so_far.
-       This function is used to select the pivot variable.
-    */
-    template<typename Ext>
-    int fixplex<Ext>::get_num_non_free_dep_vars(var_t x_j, int best_so_far) {
-        int result = is_non_free(x_j);
-        for (auto [r, _] : M.get_col(x_j)) {
-            var_t s = row2base(r);
-            result += is_non_free(s);
-            if (result > best_so_far)
-                return result;
-        }
-        return result;
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::add_patch(var_t v) {
-        SASSERT(is_base(v));
-        CTRACE("polysat", !in_bounds(v), tout << "Add patch: v" << v << "\n";);
-        if (!in_bounds(v)) 
-            m_to_patch.insert(v);
-    }
-
-    template<typename Ext>
-    var_t fixplex<Ext>::select_var_to_fix() {
-        switch (pivot_strategy()) {
-        case S_BLAND:
-            return select_smallest_var();
-        case S_GREATEST_ERROR:
-            return select_error_var(false);
-        case S_LEAST_ERROR:
-            return select_error_var(true);
-        default:
-            return select_smallest_var();
-        }
-    }
-
-    template<typename Ext>
-    var_t fixplex<Ext>::select_error_var(bool least) {
-        var_t best = null_var;
-        numeral best_error = 0;
-        for (var_t v : m_to_patch) {
-            numeral curr_error = value2error(v, value(v));
-            if (curr_error == 0)
-                continue;
-            if ((best == null_var) || 
-                (least && curr_error < best_error) ||
-                (!least && curr_error > best_error)) {
-                best = v;
-                best_error = curr_error;
-            }
-        }
-        if (best == null_var)
-            m_to_patch.clear(); // all variables are satisfied
-        else
-            m_to_patch.erase(best);
-        return best;
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::check_blands_rule(var_t v) {
-        if (m_bland) 
-            return;
-        if (!m_left_basis.contains(v)) 
-            m_left_basis.insert(v);
-        else {
-            ++m_num_repeated;
-            m_bland = m_num_repeated > m_blands_rule_threshold;
-            CTRACE("polysat", m_bland, tout << "using blands rule, " << m_num_repeated << "\n";);
-        }
-    }
-
-    /**
-     * Check if row is solved with respect to integrality constraints.
-     * The value of the row is allowed to be off by the base coefficient
-     * representing the case where there is a rational, but not integer solution.
-     */
-    template<typename Ext>                                     
-    bool fixplex<Ext>::is_solved(row const& r) const {
-        return (value(row2base(r)) * row2base_coeff(r)) + row2value(r) == 0;
-    }
-
-    /**
-     * solve for c * x + row_value = 0
-     * Cases 
-     * c = 1: x = -row_value
-     * c = -1: x = row_value
-     * Analytic solutions:
-     * trailing_zeros(c) <= trailing_zeros(row_value):
-     *   x = - inverse(c >> trailing_zeros(c)) * row_value << (trailing_zeros(row_value) - trailing_zeros(c))
-     * trailing_zeros(c) > trailing_zeros(row_value):
-     *   There is no feasible (integral) solution for x
-     *   Possible approximation:
-     *   x = - inverse(c >> trailing_zeros(c)) * row_value >> (trailing_zeros(c) - trailing_zeros(row_value))
-     * 
-     * Approximate approaches:
-     * 0 - c >= c:
-     *   * - row_value / c or (0 - row_value) / c
-     * 0 - c < c
-     *   * row_value / (0 - c) or - (0 - row_value) / (0 - c)
-     * 
-     * Analytic solution requires computing inverse (uses gcd, so multiple divisions).
-     * Approximation can be used to suppress rows that are feasible in a relaxation.
-     * Characteristics of the relaxation(s) requires further analysis.     
-     */
-    template<typename Ext>
-    typename fixplex<Ext>::numeral 
-    fixplex<Ext>::solve_for(numeral const& row_value, numeral const& c) {
-        if (c == 1)
-            return 0 - row_value;
-        if (c + 1 == 0)
-            return row_value;
-        if (0 - c < c)
-            return row_value / (0 - c);
-        return 0 - row_value / c;
-    }
-
-    template<typename Ext>                                     
-    void fixplex<Ext>::set_base_value(var_t x) {
-        SASSERT(is_base(x));
-        row r = base2row(x);
-        m_vars[x].m_value = solve_for(row2value(r), row2base_coeff(r));
-        touch_var(x);
-        m_rows[r.id()].m_integral = is_solved(r);
-        if (row_is_integral(r))
-            m_non_integral.remove(r.id());
-        else 
-            m_non_integral.insert(r.id());            
-    }
-
-    /**
-     * Equality detection.
-     *
-     * Offset equality detection
-     * -------------------------
-     * is_offset_row: determine if row is cx*x + cy*y + k == 0 where k is a constant.
-     * Then walk every row containing x, y, respectively
-     * If there is a row cx*x + cy*z + k' == 0, where y, z are two different variables
-     * but value(y) = value(z), cy is odd
-     * then it follows that k = k' and y = z is implied.
-     * 
-     * Offset equality detection is only applied to integral rows where the current  
-     * evaluation satisfies the row equality.
-     *
-     * Fixed variable equalities
-     * -------------------------
-     * Use persistent hash-table of variables that are fixed at values.
-     * Update table when a variable gets fixed and check for collisions.
-     * 
-     */
-
-    template<typename Ext>
-    bool fixplex<Ext>::propagate_row_eqs() {
-        if (!m_to_patch.empty())
-            return true;
-        if (m_eq_rows.empty())
-            return true;
-        if (!m_propagate_eqs_backoff.should_propagate())
-            return true;
-        unsigned sz = m_var_eqs.size();
-        for (unsigned i : m_eq_rows)
-            get_offset_eqs(row(i));
-        m_propagate_eqs_backoff.update(sz < m_var_eqs.size());
-        m_eq_rows.reset();
-        return !inconsistent();
-    }
-
-
-    template<typename Ext>
-    void fixplex<Ext>::get_offset_eqs(row const& r) {
-        var_t x, y;
-        numeral cx, cy;
-        if (!is_offset_row(r, cx, x, cy, y))
-            return;
-        lookahead_eq(r, cx, x, cy, y);
-        lookahead_eq(r, cy, y, cx, x);
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::is_offset_row(row const& r, numeral& cx, var_t& x, numeral& cy, var_t & y) const {
-        x = null_var;
-        y = null_var;
-        if (!row_is_integral(r))
-            return false;
-        for (auto const& e : M.get_row(r)) {
-            var_t v = e.var();
-            if (is_fixed(v))
-                continue;
-            numeral const& c = e.coeff();
-            if (x == null_var) {
-                cx = c;
-                x = v;
-            }
-            else if (y == null_var) {
-                cy = c;
-                y = v;
-            }
-            else
-                return false;
-        }        
-        return y != null_var;
-    }
-
-
-    template<typename Ext>
-    void fixplex<Ext>::lookahead_eq(row const& r1, numeral const& cx, var_t x, numeral const& cy, var_t y) {   
-        if (m.is_even(cy))
-            return;
-        var_t z, u;
-        numeral cz, cu;
-        for (auto [r2, _] : M.get_col(x)) {
-            if (r1.id() == r2.id())
-                continue;
-            if (!is_offset_row(r2, cz, z, cu, u))
-                continue;
-            if (u == x) {
-                std::swap(z, u);
-                std::swap(cz, cu);
-            }
-            if (z == x && find(u) != find(y) && cx == cz && cu == cy && value(u) == value(y)) 
-                eq_eh(u, y, m_deps.mk_join(row2dep(r1), row2dep(r2)));                
-            if (z == x && find(u) != find(y) && cx + cz == 0 && cu + cy == 0 && value(u) == value(y)) 
-                eq_eh(u, y, m_deps.mk_join(row2dep(r1), row2dep(r2)));
-        }
-    }
-
-    /**
-     * Accumulate equalities between variables fixed to the same values.
-     */
-    template<typename Ext>
-    void fixplex<Ext>::fixed_var_eh(u_dependency* dep, var_t x) {
-        numeral val = value(x);
-        fix_entry e;
-        if (m_value2fixed_var.find(val, e)) {
-            if (find(x) != find(e.x))
-                eq_eh(x, e.x, m_deps.mk_join(e.dep, dep));
-        }
-        else {
-            m_value2fixed_var.insert(val, fix_entry(x, dep));
-            m_fixed_vals.push_back(val);
-            m_trail.push_back(trail_i::fixed_val_i);
-        }
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::eq_eh(var_t x, var_t y, u_dependency* dep) {
-        SASSERT(find(x) != find(y));
-        merge(x, y);
-        m_var_eqs.push_back(var_eq(x, y, dep));
-        m_trail.push_back(trail_i::add_eq_i);
-    }   
-
-    template<typename Ext>
-    bool fixplex<Ext>::propagate_row_bounds() {
-        // return true;
-        // TBD
-        if (inconsistent())
-            return false;
-        if (!m_to_patch.empty())
-            return true;
-        if (m_bound_rows.empty())
-            return true;
-        if (!m_propagate_bounds_backoff.should_propagate())
-            return true;
-        for (unsigned i : m_bound_rows)
-            if (!propagate_row(row(i)))
-                return false;
-        m_bound_rows.reset();
-        m_propagate_bounds_backoff.update(false); // should backoff be adjusted?
-        return true;
-    }
-
-    //
-    // DFS search propagating inequalities
-    // TBDs:
-    // - replace by incremental topological sort (util/topsort.h uses stack and is non-incremental).
-    //    Then patch solutions following a pre-order processing
-    //       value(v) := max({ value(w) | w <= v } u { value(w) + 1 | w < v } u { lo(v) }) unsat if out of bounds.
-    // - combine with row propagation?
-    // - search for shorter cycles? conflicts with literals asserted at lower (glue) levels?
-    // - statistics
-    // use diff-logic propagation instead?
-    // - use heap based on potential
-    // - update gamma without overflow issues
-    //
-
-    template<typename Ext>
-    lbool fixplex<Ext>::propagate_ineqs(unsigned i0_idx) {
-        ineq i0 = m_ineqs[i0_idx];
-        numeral old_lo = m_vars[i0.w].lo;
-        SASSERT(!m_inconsistent);
-        if (!propagate_ineq(i0))
-            return l_false;
-        on_stack.reset();
-        stack.reset();
-        stack.push_back(std::make_pair(0, i0_idx));
-        on_stack.insert(i0.v);
-        while (!stack.empty()) {
-            if (!m_limit.inc())
-                return l_undef;
-            auto [ineq_out, i_idx] = stack.back();        
-            ineq i = m_ineqs[i_idx];
-            auto const& ineqs = m_var2ineqs[i.w];
-            for (; ineq_out < ineqs.size(); ++ineq_out) {
-                auto& i_out = m_ineqs[ineqs[ineq_out]];
-                if (i.w != i_out.v) 
-                    continue;
-                old_lo = m_vars[i_out.w].lo;
-                if (!propagate_ineq(i_out))
-                    return l_false;      
-                
-                bool is_onstack = on_stack.contains(i_out.w);
-                if ((old_lo != m_vars[i_out.w].lo || m_ineqs_to_propagate.contains(ineqs[ineq_out])) && !is_onstack) {
-                    on_stack.insert(i_out.w);
-                    stack.back().first = ineq_out + 1;
-                    stack.push_back(std::make_pair(0, ineqs[ineq_out]));
-                    break;
-                }
-
-                if (!is_onstack) {
-                    on_stack.insert(i_out.w);
-                    continue;
-                }
-
-                bool strict = i_out.strict, found = false, empty = false;
-                auto bound = m_vars[i_out.w];
-                unsigned j = stack.size();
-                for (; !found && j-- > 0; ) {
-                    ineq i2 = m_ineqs[stack[j].second];
-                    strict |= i2.strict;
-                    bound &= m_vars[i2.w];
-                    if (i2.v == i_out.w) 
-                        found = true;
-                    if (bound.is_empty())
-                        empty = true;
-                }
-                if ((empty || strict) && found) {
-                    auto* d = i_out.dep;
-                    for (; j < stack.size(); ++j) 
-                        d = m_deps.mk_join(d, m_ineqs[stack[j].second].dep);
-                    conflict(d);
-                    return l_false;
-                }
-            }
-            if (ineq_out == ineqs.size()) {
-                m_ineqs_to_propagate.remove(i_idx);
-                on_stack.remove(i.w);
-                stack.pop_back();
-            }
-        }
-        return l_true;
-    }
-
-    /**
-     * Bounds propagation
-     * works so far if coefficient to variable is 1 or -1
-     * Generalization is TBD:
-     *  Explore an efficient way to propagate with the following idea:
-     *  For odd c, multiply row by inverse of c and accumulate similar
-     *  propagation.
-     * 
-     * Conflicts spanning multiple rows are TBD:
-     * Idea could be similar to conflicts for inequality propagation.
-     * - create a stack of variables that get tightened.
-     * - walk over every row that contains the top variable on the stack
-     * - perform bounds propagation for the currently examined row
-     * - put newly tightened variables from row on the top of the stack
-     * - if a variable occurs already on the stack determine if the rows on the
-     *   stack resolve into a conflict.
-     * 
-     * TBD: Combination of inequality and row propagation?
-     */
-
-    template<typename Ext>
-    bool fixplex<Ext>::propagate_row(row const& r) {
-        mod_interval<numeral> range(0, 1);
-        numeral free_c = 0;
-        var_t free_v = null_var;
-        for (auto const& e : M.get_row(r)) {
-            var_t v = e.var();
-            numeral const& c = e.coeff();
-            if (is_free(v)) {
-                if (free_v != null_var)
-                    return l_true;
-                free_v = v;
-                free_c = c;
-                continue;
-            }
-            range += m_vars[v] * c;
-            if (range.is_free())
-                return true;
-        }        
-
-        if (free_v != null_var) {
-            range = (-range) * free_c;
-            bool res = new_bound(r, free_v, range);
-            SASSERT(in_bounds(free_v));
-            return res;
-        }
-        for (auto const& e : M.get_row(r)) {
-            var_t v = e.var();
-            SASSERT(!is_free(v));
-            auto range1 = range - m_vars[v] * e.coeff();
-            bool res = new_bound(r, v, range1);
-            if (res != l_true)
-                return res;
-            // SASSERT(in_bounds(v));
-        }
-        return true;
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::propagate_strict_bounds(ineq const& i) {
-        var_t v = i.v, w = i.w;
-        //bool s = i.strict;
-        auto* vlo = m_vars[v].m_lo_dep, *vhi = m_vars[v].m_hi_dep;
-        auto* wlo = m_vars[w].m_lo_dep, *whi = m_vars[w].m_hi_dep;
-
-        if (v == w)
-            return conflict(i, nullptr), false;
-
-        if (m_vars[w].max() == 0)
-            return conflict(i, wlo, whi), false;
-
-        if (m_vars[v].min() + 1 == 0)
-            return conflict(i, vlo, vhi), false;
-
-        if (m_vars[v].min() >= m_vars[w].min() && !new_bound(i, w, m_vars[v].min() + 1, 0, vlo, vhi, wlo, whi))
-            return false;
-
-        if (m_vars[v].max() >= m_vars[w].max() && !new_bound(i, v, 0, m_vars[w].max(), vlo, vhi, wlo, whi))
-            return false;
-
-        if (!is_base(w) && value(v) >= value(w) && value(v) + 1 != 0 &&  m_vars[w].contains(value(v) + 1))
-            update_value(w, value(v) - value(w) + 1);
-        
-        return true;
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::propagate_non_strict_bounds(ineq const& i) {
-        var_t v = i.v, w = i.w;
-        // bool s = i.strict;
-        auto* vlo = m_vars[v].m_lo_dep, *vhi = m_vars[v].m_hi_dep;
-        auto* wlo = m_vars[w].m_lo_dep, *whi = m_vars[w].m_hi_dep;
-
-        if (m_vars[v].min() > m_vars[w].min() && !new_bound(i, w, m_vars[v].min(), 0, vlo, vhi, wlo, whi))
-            return false;
-
-        if (m_vars[v].max() > m_vars[w].max() && !new_bound(i, v, 0, m_vars[w].max() + 1, vlo, vhi, wlo, whi))
-            return false;
-
-        if (value(v) > value(w) && !is_base(w) && m_vars[w].contains(value(v)))
-            update_value(w, value(v) - value(w)); 
-
-        return true;
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::propagate_ineq(ineq const& i) {
-        if (i.strict)
-            return propagate_strict_bounds(i);
-        else
-            return propagate_non_strict_bounds(i);
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::conflict(ineq const& i, u_dependency* a, u_dependency* b, u_dependency* c, u_dependency* d) {
-        conflict(a, m_deps.mk_join(i.dep, m_deps.mk_join(b, m_deps.mk_join(c, d))));
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::conflict(u_dependency* a) {
-        m_unsat_core.reset();
-        m_deps.linearize(a, m_unsat_core);
-        SASSERT(!m_inconsistent);
-        m_inconsistent = true;
-        ++m_stats.m_num_infeasible;
-        m_trail.push_back(trail_i::set_inconsistent_i);
-        TRACE("polysat", tout << "core: " << m_unsat_core << "\n";);
-    }
-
-    template<typename Ext>
-    u_dependency* fixplex<Ext>::row2dep(row const& r) {
-        u_dependency* d = nullptr;
-        for (auto const& e : M.get_row(r)) {
-            var_t v = e.var();
-            d = m_deps.mk_join(m_vars[v].m_lo_dep, d);
-            d = m_deps.mk_join(m_vars[v].m_hi_dep, d);
-        }
-        return d;
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::new_bound(ineq const& i, var_t x, numeral const& l, numeral const& h, u_dependency* a, u_dependency* b, u_dependency* c, u_dependency* d) {
-        bool was_fixed = is_fixed(x);
-        SASSERT(!inconsistent());
-        u_dependency* dep = m_deps.mk_join(i.dep, m_deps.mk_join(a, m_deps.mk_join(b, m_deps.mk_join(c, d))));
-        update_bounds(x, l, h, dep);
-        if (inconsistent())
-            return false;
-        if (!was_fixed && is_fixed(x)) 
-            fixed_var_eh(dep, x);
-        
-        return true;
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::new_bound(row const& r, var_t x, mod_interval<numeral> const& range) {
-        if (range.contains(m_vars[x]))
-            return true;
-        SASSERT(!inconsistent());
-        bool was_fixed = is_fixed(x);
-        u_dependency* dep = row2dep(r);
-        update_bounds(x, range.lo, range.hi, dep);
-        IF_VERBOSE(0, verbose_stream() << "new-bound v" << x << " " << m_vars[x] << "\n");
-        if (inconsistent())
-            return false;
-        if (!was_fixed && is_fixed(x))
-            fixed_var_eh(dep, x);
-        return true;
-    }
-
-    template<typename Ext>
-    std::ostream& fixplex<Ext>::display(std::ostream& out) const {
-        M.display(out);
-        for (unsigned i = 0; i < m_vars.size(); ++i) {
-            var_info const& vi = m_vars[i];
-            out << "v" << i << " " << pp(value(i)) << " " << vi << " ";
-            if (vi.m_is_base) out << "b:" << vi.m_base2row << " " << pp(m_rows[vi.m_base2row].m_value) << " ";
-            out << "\n";
-        }
-        for (ineq const& i : m_ineqs)
-            out << i << "\n";
-        return out;
-    }
-
-    template<typename Ext>
-    std::ostream& fixplex<Ext>::display_row(std::ostream& out, row const& r, bool values) const {
-        out << r.id() << " := " << pp(row2value(r)) << " : ";
-        for (auto const& e : M.get_row(r)) {
-            var_t v = e.var();
-            if (e.coeff() != 1)
-                out << pp(e.coeff()) << " * ";
-            out << "v" << v;
-            if (is_base(v))
-                out << "b";
-            out << " ";
-            if (values)
-                out << pp(value(v)) << " " << m_vars[v] << " ";
-        }
-        return out << "\n";
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::well_formed() const {
-        SASSERT(M.well_formed());
-        for (unsigned i = 0; i < m_rows.size(); ++i) {
-            row r(i);
-            var_t s = row2base(r);
-            if (s == null_var)
-                continue;
-            SASSERT(i == base2row(s).id());
-            
-            VERIFY(well_formed_row(r));
-        }
-        for (unsigned i = 0; i < m_vars.size(); ++i) {
-            SASSERT(is_base(i) || in_bounds(i));
-            if (!is_base(i) && !in_bounds(i))
-                return false;
-        }
-        return true;
-    }
-
-    template<typename Ext>
-    bool fixplex<Ext>::well_formed_row(row const& r) const {
-        var_t s = row2base(r);
-        VERIFY(base2row(s).id() == r.id());
-        VERIFY(m_vars[s].m_is_base);
-        numeral sum = 0;
-        numeral base_coeff = row2base_coeff(r);
-        for (auto const& e : M.get_row(r)) {
-            sum += value(e.var()) * e.coeff();
-            SASSERT(s != e.var() || base_coeff == e.coeff());
-        }
-        if (sum >= base_coeff) {
-            IF_VERBOSE(0, M.display_row(verbose_stream(), r););
-            TRACE("polysat", display(tout << "non-well formed row\n"); M.display_row(tout << "row: ", r););
-            throw default_exception("non-well formed row");
-        }
-
-        if (sum != row2value(r) + base_coeff * value(s)) {
-            std::cout << sum << "\n";
-            std::cout << (row2value(r) + base_coeff * value(s)) << "\n";
-            display_row(std::cout, r, false);
-            display(std::cout);
-        }
-
-        SASSERT(sum == row2value(r) + base_coeff * value(s));
-        return true;
-    }
-
-    template<typename Ext>
-    void fixplex<Ext>::collect_statistics(::statistics & st) const {
-        M.collect_statistics(st);
-        st.update("fixplex num pivots", m_stats.m_num_pivots);
-        st.update("fixplex num infeasible", m_stats.m_num_infeasible);
-        st.update("fixplex num checks", m_stats.m_num_checks);
-        st.update("fixplex num non-integral", m_non_integral.num_elems());
-        st.update("fixplex num approximated row additions", m_stats.m_num_approx);
-    }
-}
diff --git a/src/math/polysat/fixplex_mod_interval.h b/src/math/polysat/fixplex_mod_interval.h
deleted file mode 100644
index 65351c615..000000000
--- a/src/math/polysat/fixplex_mod_interval.h
+++ /dev/null
@@ -1,111 +0,0 @@
-/*++
-Copyright (c) 2014 Microsoft Corporation
-
-Module Name:
-
-    mod_interval.h
-
-Abstract:
-
-    Intervals over fixed precision modular arithmetic
-
-Author:
-
-    Nikolaj Bjorner (nbjorner) 2021-03-19
-    Jakob Rath 2021-04-6
-
---*/
-
-#pragma once
-
-#include "util/rational.h"
-
-
-template<typename Numeral>
-struct pp {
-    Numeral n;
-    pp(Numeral const& n):n(n) {}
-};
-
-
-inline std::ostream& operator<<(std::ostream& out, pp<uint8_t> const& p) {
-    return out << (unsigned)p.n;
-}
-
-template<typename Numeral>
-inline std::ostream& operator<<(std::ostream& out, pp<Numeral> const& p) {
-    if ((Numeral)(0 - p.n) < p.n)
-        return out << "-" << (Numeral)(0 - p.n);
-    return out << p.n;
-}
-
-inline std::ostream& operator<<(std::ostream& out, pp<rational> const& p) {
-    return out << p.n;
-}
-
-template<typename Numeral>
-class mod_interval {
-    bool emp = false;
-public:
-    Numeral lo { 0 };
-    Numeral hi { 0 };
-    mod_interval() {}
-    mod_interval(Numeral const& l, Numeral const& h): lo(l), hi(h) {}
-    virtual ~mod_interval() {}
-    static mod_interval free() { return mod_interval(0, 0); }
-    static mod_interval empty() { mod_interval i(0, 0); i.emp = true; return i; }
-
-    bool is_free() const { return !emp && lo == hi; }
-    bool is_empty() const { return emp; }
-    bool is_singleton() const { return !is_empty() && (lo + 1 == hi || (hi == 0 && is_max(lo))); }
-    bool contains(Numeral const& n) const;
-    bool contains(mod_interval const& other) const;
-    virtual bool is_max(Numeral const& n) const { return (Numeral)(n + 1) == 0; }
-    Numeral max() const;
-    Numeral min() const; 
-
-    void set_free() { lo = hi = 0; emp = false; }
-    void set_bounds(Numeral const& l, Numeral const& h) { lo = l; hi = h; }
-    void set_empty() { emp = true; }
-
-    mod_interval& intersect_ule(Numeral const& h);
-    mod_interval& intersect_uge(Numeral const& l);
-    mod_interval& intersect_ult(Numeral const& h);
-    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;
-    mod_interval operator-(mod_interval const& other) const;
-    mod_interval operator*(mod_interval const& other) const;
-    mod_interval operator-() const;
-    mod_interval operator*(Numeral const& n) const;
-    mod_interval operator+(Numeral const& n) const { return mod_interval(lo + n, hi + n); }
-    mod_interval operator-(Numeral const& n) const { return mod_interval(lo - n, hi - n); }
-    mod_interval& operator+=(mod_interval const& other) { *this = *this + other; return *this; }
-    std::ostream& display(std::ostream& out) const { 
-        if (is_empty()) return out << "empty"; 
-        if (is_free()) return out << "free";
-        return out << "[" << pp(lo) << ", " << pp(hi) << "["; 
-    }
-    Numeral closest_value(Numeral const& n) const;
-    bool operator==(mod_interval const& other) const {
-        if (is_empty())
-            return other.is_empty();
-        if (is_free())
-            return other.is_free();
-        return lo == other.lo && hi == other.hi;
-    }
-    bool operator!=(mod_interval const& other) const {
-        return !(*this == other);
-    }
-};
-
-template<typename Numeral>
-inline std::ostream& operator<<(std::ostream& out, mod_interval<Numeral> const& i) {
-    return i.display(out);
-}
-
diff --git a/src/math/polysat/fixplex_mod_interval_def.h b/src/math/polysat/fixplex_mod_interval_def.h
deleted file mode 100644
index 0bf832432..000000000
--- a/src/math/polysat/fixplex_mod_interval_def.h
+++ /dev/null
@@ -1,333 +0,0 @@
-/*++
-Copyright (c) 2014 Microsoft Corporation
-
-Module Name:
-
-    mod_interval_def.h
-
-Abstract:
-
-    Intervals over fixed precision modular arithmetic
-
-Author:
-
-    Nikolaj Bjorner (nbjorner) 2021-03-19
-    Jakob Rath 2021-04-6
-
---*/
-
-#pragma once
-
-#include "fixplex_mod_interval.h"
-
-template<typename Numeral>
-bool mod_interval<Numeral>::contains(Numeral const& n) const {
-    if (is_empty())
-        return false;
-    if (is_free())
-        return true;
-    if (lo < hi)
-        return lo <= n && n < hi;
-    else
-        return lo <= n || n < hi;
-}
-
-template<typename Numeral>
-bool mod_interval<Numeral>::contains(mod_interval const& other) const {
-    if (is_empty())
-        return other.is_empty();
-    if (is_free())
-        return true;
-    if (hi == 0)
-        return lo <= other.lo && (other.lo < other.hi || other.hi == 0);
-    if (lo < hi)
-        return lo <= other.lo && other.hi <= hi;
-    if (other.lo < other.hi && other.hi <= hi)
-        return true;
-    if (other.lo < other.hi && lo <= other.lo)
-        return true;
-    if (other.hi == 0)
-        return lo <= other.lo;
-    SASSERT(other.hi < other.lo && other.hi != 0);
-    SASSERT(hi < lo && hi != 0);
-    return lo <= other.lo && other.hi <= hi;
-}
-
-template<typename Numeral>
-mod_interval<Numeral> mod_interval<Numeral>::operator+(mod_interval<Numeral> const& other) const {
-    if (is_empty())
-        return *this;
-    if (other.is_empty())
-        return other;
-    if (is_free())
-        return *this;
-    if (other.is_free())
-        return other;
-    Numeral sz = (hi - lo) + (other.hi - other.lo);
-    if (sz < (hi - lo))
-        return mod_interval::free();
-    return mod_interval(lo + other.lo, hi + other.hi);
-}
-
-template<typename Numeral>
-mod_interval<Numeral> mod_interval<Numeral>::operator-(mod_interval<Numeral> const& other) const {
-    return *this + (-other);
-}
-
-template<typename Numeral>
-mod_interval<Numeral> mod_interval<Numeral>::operator-() const {
-    if (is_empty())
-        return *this;
-    if (is_free())
-        return *this;
-    return mod_interval(1 - hi, 1 - lo);
-}
-
-template<typename Numeral>
-mod_interval<Numeral> mod_interval<Numeral>::operator*(Numeral const& n) const {
-    if (is_empty())
-        return *this;
-    if (n == 0)
-        return mod_interval(0, 1);
-    if (n == 1)
-        return *this;
-    if (is_free())
-        return *this;
-    Numeral sz = hi - lo;
-    if (0 - n < n) {
-        Numeral mn = 0 - n;
-        Numeral mz = mn * sz;            
-        if (mz / mn != sz)
-            return mod_interval::free();
-        return mod_interval((hi - 1) * n, n * lo + 1);
-    }
-    else {
-        Numeral mz = n * sz;
-        if (mz / n != sz)
-            return mod_interval::free();
-        return mod_interval(n * lo, n * (hi - 1) + 1);
-    }
-}
-
-template<typename Numeral>
-mod_interval<Numeral> mod_interval<Numeral>::operator&(mod_interval const& other) const {
-    Numeral l, h;
-    if (is_free() || other.is_empty())
-        return other;
-    if (other.is_free() || is_empty())
-        return *this;
-
-    if (lo < hi || hi == 0) {
-        if (other.lo < other.hi || other.hi == 0) {
-            if (hi != 0 && hi <= other.lo)
-                return mod_interval::empty();
-            if (other.hi != 0 && other.hi <= lo)
-                return mod_interval::empty();
-            l = lo >= other.lo ? lo : other.lo;
-            h = hi == 0 ? other.hi : (other.hi == 0 ? hi : (hi <= other.hi ? hi : other.hi));
-            return mod_interval(l, h);
-        }
-        SASSERT(0 < other.hi && other.hi < other.lo);
-        if (other.lo <= lo)
-            return *this;
-        if (other.hi <= lo && lo < hi && hi <= other.lo)
-            return mod_interval::empty();
-        if (lo <= other.hi && other.hi <= hi && hi <= other.lo)
-            return mod_interval(lo, other.hi);
-        if (hi == 0 && lo < other.hi)
-            return *this;
-        if (hi == 0 && other.hi <= lo)
-            return mod_interval(other.lo, hi);
-        if (other.hi <= lo && other.hi <= hi)
-            return mod_interval(other.lo, hi);
-        return *this;
-    }
-    SASSERT(hi < lo);
-    if (other.lo < other.hi || other.hi == 0)
-        return other & *this;
-    SASSERT(other.hi < other.lo);
-    SASSERT(hi != 0);
-    SASSERT(other.hi != 0);
-    if (lo <= other.hi)
-        return *this;
-    if (other.lo <= hi)
-        return other;
-    l = lo <= other.lo ? other.lo : lo;
-    h = hi >= other.hi ? other.hi : hi;
-    return mod_interval(l, h);
-
-}
-
-template<typename Numeral>
-Numeral mod_interval<Numeral>::max() const {
-    if (lo < hi)
-        return hi - 1;
-    else
-        return Numeral(0) - 1;
-}
-
-template<typename Numeral>
-Numeral mod_interval<Numeral>::min() const {
-    if (lo < hi || hi == 0)
-        return lo;
-    else
-        return Numeral(0);
-}
-
-
-template<typename Numeral>
-Numeral mod_interval<Numeral>::closest_value(Numeral const& n) const {
-    if (contains(n))
-        return n;
-    if (is_empty())
-        return n;
-    if ((Numeral)(lo - n) < (Numeral)(n - hi))
-        return lo;
-    return hi - 1;
-}
-
-template<typename Numeral>
-mod_interval<Numeral>& mod_interval<Numeral>::intersect_uge(Numeral const& l) {
-    if (is_empty())
-        return *this;
-    if (lo < hi && hi <= l)
-        set_empty();
-    else if (is_free())
-        lo = l, hi = 0;
-    else if ((lo < hi || hi == 0) && lo < l)
-        lo = l;
-    else if (hi < lo && hi <= l && l <= lo)
-        hi = 0;
-    else if (hi < lo && lo < l)
-        hi = 0, lo = l;
-    else if (0 < l && l < hi && hi < lo)
-        lo = l, hi = 0;
-    return *this;
-}
-
-template<typename Numeral>
-mod_interval<Numeral>& mod_interval<Numeral>::intersect_ugt(Numeral const& l) {
-    if (is_empty())
-        return *this;
-    if (is_max(l))
-        set_empty();
-    else if (is_free())
-        lo = l + 1, hi = 0;
-    else if (lo < hi && hi - 1 <= l)
-        set_empty();
-    else if (lo < hi && l < lo)
-        return *this;
-    else if (lo < hi)
-        lo = l + 1;
-    else if (hi < lo && hi <= l + 1 && l < lo)
-        hi = 0;
-    else if (hi < lo && lo <= l)
-        hi = 0, lo = l + 1;
-    else if (l <= hi && hi < lo)
-        lo = l + 1, hi = 0;
-    return *this;
-}
-
-template<typename Numeral>
-mod_interval<Numeral>& mod_interval<Numeral>::intersect_ule(Numeral const& h) {
-    if (is_empty())
-        return *this;
-    if (is_max(h))
-        return *this;
-    else if (is_free())
-        lo = 0, hi = h + 1;
-    else if (h < lo && (lo < hi || hi == 0))
-        set_empty();
-    else if (hi == 0 && h >= lo)
-        hi = h + 1;
-    else if (lo <= h && h + 1 < hi)
-        hi = h + 1;
-    else if (h < hi && hi < lo)
-        hi = h + 1, lo = 0;
-    else if (hi <= h && h < lo)
-        lo = 0;
-    else if (hi == 0 && hi == h && hi < lo)
-        set_empty();
-    else if (0 < hi && hi == h && hi < lo)
-        lo = 0;
-    else if (0 < hi && hi < h && hi < lo)
-        lo = 0, hi = h + 1;
-    return *this;
-}
-
-template<typename Numeral>
-mod_interval<Numeral>& mod_interval<Numeral>::intersect_ult(Numeral const& h) {
-    if (is_empty())
-        return *this;
-    if (h == 0)
-        set_empty();
-    else if (is_free())
-        lo = 0, hi = h;
-    else if (h <= lo && (lo < hi || hi == 0))
-        set_empty();
-    else if (h > lo && (h < hi || hi == 0))
-        hi = h;
-    else if (hi < lo && h <= hi)
-        hi = h, lo = 0;
-    else if (hi < h && h <= lo)
-        lo = 0;
-    else if (0 < hi && hi < lo && hi + 1 <= h)
-        lo = 0, hi = h;
-    return *this;
-}
-
-
-template<typename Numeral>
-mod_interval<Numeral>& mod_interval<Numeral>::intersect_fixed(Numeral const& a) {
-    if (is_empty())
-        return *this;
-    if (!contains(a))
-        set_empty();
-    else if (is_max(a))
-        lo = a, hi = 0;
-    else
-        lo = a, hi = a + 1;
-    return *this;
-}
-
-template<typename Numeral>
-mod_interval<Numeral>& mod_interval<Numeral>::intersect_diff(Numeral const& a) {
-    if (!contains(a) || is_empty())
-        return *this;
-    if (a == lo && a + 1 == hi)
-        set_empty();
-    else if (a == lo && hi == 0 && is_max(a))
-        set_empty();
-    else if (is_free())
-        lo = a + 1, hi = a;
-    else if (0 < hi && hi < lo && a == lo)        
-        return *this;
-    else if (a == lo && !is_max(a))
-        lo = a + 1;
-    else if (a + 1 == hi)
-        hi = a;
-    else if (hi == 0 && is_max(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/linear_solver.cpp b/src/math/polysat/linear_solver.cpp
deleted file mode 100644
index 2ae28fc04..000000000
--- a/src/math/polysat/linear_solver.cpp
+++ /dev/null
@@ -1,288 +0,0 @@
-/*++
-Copyright (c) 2014 Microsoft Corporation
-
-Module Name:
-
-    linear_solver.cpp
-    
-Author:
-
-    Nikolaj Bjorner (nbjorner) 2021-05-14
-    Jakob Rath 2021-05-14
-
---*/
-
-#include "math/bigfix/u256.h"
-#include "math/polysat/linear_solver.h"
-#include "math/polysat/fixplex_def.h"
-#include "math/polysat/solver.h"
-
-inline unsigned numeral2hash(u256 const& n) { return n.hash(); }
-
-namespace polysat {
-
-    void linear_solver::push() {
-        m_trail.push_back(trail_i::inc_level_i);
-        for (auto* f : m_fix_ptr)
-            f->push();
-    }
-
-    void linear_solver::pop(unsigned n) {
-        for (auto* f : m_fix_ptr)
-            f->pop(n);
-        while (n > 0) {
-            switch (m_trail.back()) {
-            case trail_i::inc_level_i:
-                --n;
-                break;
-            case trail_i::add_var_i: {
-                auto [v, sz] = m_rows.back();
-                --m_sz2num_vars[sz];
-                m_rows.pop_back();
-                break;
-            }
-            case trail_i::add_mono_i: {
-                auto m = m_monomials.back();
-                m_mono2var.erase(m);
-                m_alloc.deallocate(m.num_vars*sizeof(unsigned), m.vars);
-                m_monomials.pop_back();
-                break;
-            }
-            case trail_i::set_active_i: 
-                m_active.pop_back();
-                break;            
-                                    
-            default:
-                UNREACHABLE();
-            }
-            m_trail.pop_back();
-        }
-        m_unsat_f = nullptr;
-    }
-    
-    fixplex_base* linear_solver::sz2fixplex(unsigned sz) { 
-        fixplex_base* b = m_fix.get(sz, nullptr);
-        if (!b) {
-            switch (sz) {
-            case 8:
-                b = alloc(fixplex<generic_uint_ext<unsigned char>>, s.params(), s.m_lim);
-                break;
-            case 16:
-                b = alloc(fixplex<generic_uint_ext<unsigned short>>, s.params(), s.m_lim);
-                break;
-            case 32:
-                b = alloc(fixplex<generic_uint_ext<unsigned>>, s.params(), s.m_lim);
-                break;
-            case 64:
-                b = alloc(fixplex<uint64_ext>, s.params(), s.m_lim);
-                break;
-            case 128:
-                NOT_IMPLEMENTED_YET();
-                break;
-            case 256:
-                b = alloc(fixplex<generic_uint_ext<u256>>, s.params(), s.m_lim);
-                break;
-            default:
-                NOT_IMPLEMENTED_YET();
-                break;
-            }
-            if (b)
-                m_fix_ptr.push_back(b);
-            m_fix.set(sz, b);
-        }
-        return b;
-    }  
-
-
-    var_t linear_solver::internalize_pdd(pdd const& p) {
-        unsigned sz = p.power_of_2();
-        linearize(p);
-        if (m_vars.size() == 1 && m_coeffs.back() == 1) 
-            return m_vars.back();
-        var_t v = fresh_var(sz);
-        m_vars.push_back(v);
-        m_coeffs.push_back(rational::power_of_two(sz) - 1);
-        auto* f = sz2fixplex(sz);
-        if (f)
-            f->add_row(v, m_vars.size(), m_vars.data(), m_coeffs.data());
-        return v;
-    }
-    
-    /**
-     * create the row c.p() - v == 0
-     * when equality is asserted, set range on v as v == 0 or v > 0.
-     */
-        
-
-
-    void linear_solver::new_le(ule_constraint& c) {
-        var_t v = internalize_pdd(c.lhs());
-        var_t w = internalize_pdd(c.rhs());
-        auto pr = std::make_pair(v, w);
-        m_bool_var2row.setx(c.bvar(), pr, pr);
-    }   
-
-    //
-    // v <= w:
-    // static constraints:
-    // - lo(v) <= lo(w)
-    // - hi(v) <= hi(w)
-    //
-    // special case for inequalities with constant bounds
-    // bounds propagation on fp, then bounds strengthening
-    // based on static constraints
-    // internal backtrack search over bounds
-    // inequality graph (with offsets)
-    // 
-
-    void linear_solver::assert_le(bool is_positive, ule_constraint& c) {
-        auto [v, w] = m_bool_var2row[c.bvar()];
-        unsigned sz = c.lhs().power_of_2();
-        auto* fp = sz2fixplex(sz);     
-        if (!fp)
-            return;
-        unsigned c_dep = constraint_idx2dep(m_active.size() - 1);
-        rational z(0);
-        if (c.rhs().is_val()) {
-            bool is_max_value = false; 
-            if (is_positive) 
-                // v <= rhs
-                fp->set_bounds(v, z, c.rhs().val(), c_dep);                 
-            else if (is_max_value) 
-                throw default_exception("conflict not implemented");
-            else 
-                // rhs < v
-                fp->set_bounds(v, c.rhs().val() + 1, z, c_dep);
-            return;
-        }
-
-        if (c.lhs().is_val()) {
-            if (is_positive)
-                // w >= lhs 
-                fp->set_bounds(w, c.lhs().val(), z, c_dep);
-            else if (c.lhs().val() == 0)
-                throw default_exception("conflict not implemented");
-            else 
-                // w < lhs 
-                fp->set_bounds(w, z, c.lhs().val() - 1, c_dep);
-            return;
-        }
-
-        if (is_positive) 
-            fp->add_le(v, w, c_dep);
-        else 
-            fp->add_lt(w, v, c_dep);
-    }
-
-
-    void linear_solver::new_constraint(constraint& c) {
-        if (c.is_ule())
-            new_le(c.to_ule());
-    }
-
-    void linear_solver::activate_constraint(bool is_positive, constraint& c) {
-        if (!c.is_ule())
-            return;
-        m_active.push_back(&c);
-        m_trail.push_back(trail_i::set_active_i);
-        assert_le(is_positive, c.to_ule());
-    }
-
-    void linear_solver::linearize(pdd const& p) {
-        unsigned sz = p.power_of_2();
-        m_vars.reset();
-        m_coeffs.reset();
-        for (auto const& m : p) {
-            m_vars.push_back(mono2var(sz, m.vars));
-            m_coeffs.push_back(m.coeff);
-        }        
-    }
-
-    var_t linear_solver::mono2var(unsigned sz, unsigned_vector const& var) {
-        mono_info m(sz, var.size(), var.data()), m1;
-        if (m_mono2var.find(m, m1))
-            return m1.var;
-        m.vars = static_cast<unsigned*>(m_alloc.allocate(var.size()*sizeof(unsigned)));
-        for (unsigned i = 0; i < var.size(); ++i)
-            m.vars[i] = var[i];
-        m.var = fresh_var(sz);
-        m_mono2var.insert(m);
-        m_monomials.push_back(m);
-        m_trail.push_back(trail_i::add_mono_i);
-        return m.var;
-    }
-
-    var_t linear_solver::pvar2var(unsigned sz, pvar v) {
-        unsigned_vector vars;
-        vars.push_back(v);
-        return mono2var(sz, vars);
-    }
-
-    var_t linear_solver::fresh_var(unsigned sz) {
-        m_sz2num_vars.reserve(sz + 1);
-        m_trail.push_back(trail_i::add_var_i);
-        m_rows.push_back(std::make_pair(0, sz));
-        return m_sz2num_vars[sz]++;
-    }
-
-    void linear_solver::set_value(pvar v, rational const& value, unsigned dep) {
-        unsigned sz = s.size(v);
-        auto* fp = sz2fixplex(sz);
-        if (!fp)
-            return;
-        var_t w = pvar2var(sz, v);
-        fp->set_value(w, value, external_dep2dep(dep));
-    }
-
-    void linear_solver::set_bound(pvar v, rational const& lo, rational const& hi, unsigned dep) {
-        unsigned sz = s.size(v);
-        auto* fp = sz2fixplex(sz);
-        if (!fp)
-            return;
-        var_t w = pvar2var(sz, v);
-        fp->set_bounds(w, lo, hi, external_dep2dep(dep));
-    }
-    
-    /** 
-    * check integer modular feasibility under current bounds.
-    * and inequalities
-    */
-    lbool linear_solver::check() { 
-        lbool res = l_true;
-        m_unsat_f = nullptr;
-        for (auto* fp : m_fix_ptr) {
-            lbool r = fp->make_feasible();
-            if (r == l_false) {
-                m_unsat_f = fp;
-                return r;
-            }
-            if (r == l_undef)
-                res = l_undef;
-        }
-        return res; 
-    }
-
-    void linear_solver::unsat_core(ptr_vector<constraint>& cs, unsigned_vector& deps) {
-        SASSERT(m_unsat_f);
-        deps.reset();
-        cs.reset();
-        for (unsigned dep : m_unsat_f->get_unsat_core()) {
-            if (is_constraint_dep(dep))
-                cs.push_back(m_active[dep2constraint_idx(dep)]);
-            else
-                deps.push_back(dep2external_dep(dep));
-        }
-    }
-
-    // current value assigned to (linear) variable according to tableau.
-    bool linear_solver::value(pvar v, rational& val) { 
-        unsigned sz = s.size(v);
-        auto* fp = sz2fixplex(sz);
-        if (!fp)
-            return false;
-        val = fp->get_value(pvar2var(sz, v));
-        return true;
-    }
-    
-};
-
diff --git a/src/math/polysat/linear_solver.h b/src/math/polysat/linear_solver.h
deleted file mode 100644
index 0e8f93e80..000000000
--- a/src/math/polysat/linear_solver.h
+++ /dev/null
@@ -1,144 +0,0 @@
-/*++
-Copyright (c) 2014 Microsoft Corporation
-
-Module Name:
-
-    linear_solver.h
-
-Abstract:
-
-    Fixed-precision unsigned integer linear solver
-
-    This wraps around fixplex and translates polynomials
-    into linear form.
-    Equalities, Inequalities are converted into rows and slack variables.
-    Slack variables are bounded when constraints are activated.
-    It also handles backtracking state: bounds that are activated inside one 
-    scope are de-activated when exiting the scope.
-    
-Author:
-
-    Nikolaj Bjorner (nbjorner) 2021-05-14
-    Jakob Rath 2021-05-14
-
---*/
-
-#pragma once
-
-#include "util/hashtable.h"
-#include "util/small_object_allocator.h"
-#include "math/polysat/fixplex.h"
-#include "math/polysat/constraint.h"
-#include "math/polysat/ule_constraint.h"
-
-namespace polysat {
-
-    class solver;
-
-    class linear_solver {
-        enum trail_i {
-            inc_level_i,
-            add_var_i,
-            add_mono_i,
-            set_active_i
-        };
-
-        struct mono_info {
-            unsigned  sz { 0 };
-            unsigned  num_vars { 0 }; 
-            unsigned* vars { nullptr };
-            unsigned  var { 0 };            
-            mono_info(unsigned sz, unsigned num_vars, unsigned* vars):
-                sz(sz),
-                num_vars(num_vars),
-                vars(vars)
-            {}
-            mono_info() {}
-            struct hash {
-                unsigned operator()(mono_info const& i) const {
-                    // TODO
-                    return 0;
-                }
-            };
-            struct eq {
-                 bool operator()(mono_info const& a, mono_info const& b) const {
-                     if (a.num_vars != b.num_vars)
-                         return false;
-                     for (unsigned i = 0; i < a.num_vars; ++i)
-                         if (a.vars[i] != b.vars[i])
-                             return false;
-                     return a.sz == b.sz;
-                }
-            };
-        };
-
-        solver&                  s;
-        scoped_ptr_vector<fixplex_base> m_fix;
-        ptr_vector<fixplex_base> m_fix_ptr;
-        svector<trail_i>         m_trail;
-        svector<std::pair<var_t, unsigned>> m_rows; 
-        unsigned_vector          m_var2ext;
-        unsigned_vector          m_ext2var;
-        ptr_vector<constraint>   m_active;
-
-        svector<var_t>           m_vars;
-        vector<rational>         m_coeffs;
-        svector<std::pair<var_t, var_t>> m_bool_var2row;
-
-        hashtable<mono_info, mono_info::hash, mono_info::eq> m_mono2var;
-        unsigned_vector          m_sz2num_vars;
-        small_object_allocator   m_alloc;
-        svector<mono_info>       m_monomials;
-        fixplex_base*            m_unsat_f = nullptr;
-
-        fixplex_base* sz2fixplex(unsigned sz);
-
-        void linearize(pdd const& p);
-        var_t fresh_var(unsigned sz);
-
-        var_t internalize_pdd(pdd const& p);
-        void new_le(ule_constraint& le);
-        void assert_le(bool is_positive, ule_constraint& le);
-
-        // bind monomial to variable.
-        var_t mono2var(unsigned sz, unsigned_vector const& m);
-        var_t pvar2var(unsigned sz, pvar v);
-        // unsigned_vector var2mono(unsigned sz, var_t v) { throw default_exception("nyi"); }
-
-        // distinguish constraint and justification dependencies
-        unsigned external_dep2dep(unsigned dep) const { return dep * 2;  }
-        unsigned constraint_idx2dep(unsigned idx) const { return 1 + (idx * 2);  }
-        bool is_constraint_dep(unsigned dep) const { return 0 != (dep & 0x1); }
-        unsigned dep2constraint_idx(unsigned dep) const { return dep >> 2; }
-        unsigned dep2external_dep(unsigned dep) const { return dep >> 2; }
-
-    public:
-        linear_solver(solver& s):
-            s(s),
-            m_alloc("mononials")
-        {}
-
-        void push();
-        void pop(unsigned n); 
-        void new_constraint(constraint& c);
-        void set_value(pvar v, rational const& value, unsigned dep);
-        void set_bound(pvar v, rational const& lo, rational const& hi, unsigned dep);
-        void activate_constraint(bool is_positive, constraint& c);
-        void activate_constraint(signed_constraint c) { activate_constraint(c.is_positive(), *c.get()); }
-
-        // check integer modular feasibility under current bounds.
-        lbool check();
-        void unsat_core(ptr_vector<constraint>& constraints, unsigned_vector& deps);
-
-        // current value assigned to (linear) variable according to tableau.
-        bool value(pvar v, rational& val);
-
-        std::ostream& display(std::ostream& out) const { return out; }
-        void collect_statistics(::statistics & st) const {}
-
-    };
-
-
-
-};
-
diff --git a/src/math/polysat/pvar_queue.h b/src/math/polysat/pvar_queue.h
index dfb080d13..259ef6036 100644
--- a/src/math/polysat/pvar_queue.h
+++ b/src/math/polysat/pvar_queue.h
@@ -15,6 +15,8 @@ Author:
 
 Revision History:
 
+NSB: I made util/var_queue.h a template. You can feed it a struct comprising of size/activity.
+
 --*/
 #pragma once
 
diff --git a/src/math/polysat/solver.h b/src/math/polysat/solver.h
index a5c1f983b..00669a44c 100644
--- a/src/math/polysat/solver.h
+++ b/src/math/polysat/solver.h
@@ -139,7 +139,6 @@ namespace polysat {
         friend class inference_engine;
         friend class file_inference_logger;
         friend class forbidden_intervals;
-        friend class linear_solver;
         friend class viable;
         friend class viable_fallback;
         friend class search_state;