From 44679d8f5be7b3b09ca32326937f422d7d1d9056 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Fri, 16 Oct 2020 10:49:46 -0700
Subject: [PATCH] arith_solver (#4733)

* porting arithmetic solver

* integrating arithmetic

* lp

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

* na

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

* na

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

* na

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
---
 scripts/mk_project.py             |    2 +-
 src/ast/arith_decl_plugin.cpp     |   61 ++
 src/ast/arith_decl_plugin.h       |    6 +
 src/ast/euf/euf_egraph.cpp        |    6 +-
 src/ast/euf/euf_egraph.h          |    2 +-
 src/ast/static_features.cpp       |   15 +-
 src/ast/static_features.h         |    2 +
 src/math/lp/lar_solver.h          |    2 +-
 src/math/lp/lp_api.h              |  131 +++
 src/math/lp/lp_settings.h         |   21 +
 src/sat/smt/arith_axioms.cpp      |  376 ++++++--
 src/sat/smt/arith_diagnostics.cpp |   52 +-
 src/sat/smt/arith_internalize.cpp |  597 +++++++++++-
 src/sat/smt/arith_solver.cpp      | 1426 ++++++++++++++++++++++++++++-
 src/sat/smt/arith_solver.h        |  400 +++++++-
 src/sat/smt/array_axioms.cpp      |    8 +-
 src/sat/smt/array_internalize.cpp |   23 +-
 src/sat/smt/array_model.cpp       |    5 -
 src/sat/smt/array_solver.cpp      |   10 +
 src/sat/smt/array_solver.h        |    7 +-
 src/sat/smt/bv_internalize.cpp    |    6 -
 src/sat/smt/bv_solver.cpp         |    1 -
 src/sat/smt/bv_solver.h           |    1 -
 src/sat/smt/euf_internalize.cpp   |   27 +-
 src/sat/smt/euf_model.cpp         |    2 +-
 src/sat/smt/euf_relevancy.cpp     |    8 +
 src/sat/smt/euf_solver.cpp        |   15 +-
 src/sat/smt/euf_solver.h          |   11 +-
 src/sat/smt/sat_th.cpp            |   22 +
 src/sat/smt/sat_th.h              |   77 +-
 src/sat/tactic/goal2sat.cpp       |    4 +-
 src/smt/theory_lra.cpp            |  227 +----
 src/util/trail.h                  |   22 +
 33 files changed, 3172 insertions(+), 403 deletions(-)
 create mode 100644 src/math/lp/lp_api.h

diff --git a/scripts/mk_project.py b/scripts/mk_project.py
index c1655d027..0d99dbc39 100644
--- a/scripts/mk_project.py
+++ b/scripts/mk_project.py
@@ -49,7 +49,7 @@ def init_project_def():
     add_lib('core_tactics', ['tactic', 'macros', 'normal_forms', 'rewriter', 'pattern'], 'tactic/core')
     add_lib('arith_tactics', ['core_tactics', 'sat'], 'tactic/arith')
 
-    add_lib('sat_smt', ['sat', 'euf', 'tactic', 'solver', 'smt_params', 'bit_blaster', 'fpa', 'normal_forms'], 'sat/smt')
+    add_lib('sat_smt', ['sat', 'euf', 'tactic', 'solver', 'smt_params', 'bit_blaster', 'fpa', 'normal_forms', 'lp'], 'sat/smt')
     add_lib('sat_tactic', ['tactic', 'sat', 'solver', 'sat_smt'], 'sat/tactic')
     add_lib('nlsat_tactic', ['nlsat', 'sat_tactic', 'arith_tactics'], 'nlsat/tactic')
     add_lib('subpaving_tactic', ['core_tactics', 'subpaving'], 'math/subpaving/tactic')
diff --git a/src/ast/arith_decl_plugin.cpp b/src/ast/arith_decl_plugin.cpp
index 0ac922ecf..bf6b53077 100644
--- a/src/ast/arith_decl_plugin.cpp
+++ b/src/ast/arith_decl_plugin.cpp
@@ -830,6 +830,8 @@ bool arith_util::is_considered_uninterpreted(func_decl* f, unsigned n, expr* con
     return plugin().is_considered_uninterpreted(f);
 }
 
+
+
 func_decl* arith_util::mk_ipower0() {
     sort* s = mk_int();
     sort* rs[2] = { s, s };
@@ -861,3 +863,62 @@ func_decl* arith_util::mk_mod0() {
     sort* rs[2] = { mk_int(), mk_int() };
     return m_manager.mk_func_decl(m_afid, OP_MOD0, 0, nullptr, 2, rs, mk_int());
 }
+
+bool arith_util::is_bounded(expr* n) const {
+    expr* x = nullptr, * y = nullptr;
+    while (true) {
+        if (is_idiv(n, x, y) && is_numeral(y)) {
+            n = x;
+        }
+        else if (is_mod(n, x, y) && is_numeral(y)) {
+            return true;
+        }
+        else if (is_numeral(n)) {
+            return true;
+        }
+        else {
+            return false;
+        }
+    }
+}
+
+bool arith_util::is_extended_numeral(expr* term, rational& r) const {
+    rational mul(1);
+    do {
+        if (is_numeral(term, r)) {
+            r *= mul;
+            return true;
+        }
+        if (is_uminus(term, term)) {
+            mul.neg();
+            continue;
+        }
+        if (is_to_real(term, term)) {
+            continue;
+        }
+        return false;
+    } while (false);
+    return false;
+}
+
+bool arith_util::is_underspecified(expr* e) const {
+    if (!is_app(e))
+        return false;
+    if (to_app(e)->get_family_id() == get_family_id()) {
+        switch (to_app(e)->get_decl_kind()) {
+        case OP_DIV:
+        case OP_IDIV:
+        case OP_REM:
+        case OP_MOD:
+        case OP_DIV0:
+        case OP_IDIV0:
+        case OP_REM0:
+        case OP_MOD0:
+            return true;
+        default:
+            break;
+        }
+    }
+    return false;
+}
+
diff --git a/src/ast/arith_decl_plugin.h b/src/ast/arith_decl_plugin.h
index ada150e27..ff6d21ab1 100644
--- a/src/ast/arith_decl_plugin.h
+++ b/src/ast/arith_decl_plugin.h
@@ -495,6 +495,12 @@ public:
 
     bool is_considered_uninterpreted(func_decl* f, unsigned n, expr* const* args, func_decl_ref& f_out);
 
+    bool is_underspecified(expr* e) const;
+
+    bool is_bounded(expr* e) const;
+
+    bool is_extended_numeral(expr* e, rational& r) const;
+
 };
 
 
diff --git a/src/ast/euf/euf_egraph.cpp b/src/ast/euf/euf_egraph.cpp
index 6c3cd58c1..2fc4629c9 100644
--- a/src/ast/euf/euf_egraph.cpp
+++ b/src/ast/euf/euf_egraph.cpp
@@ -324,6 +324,7 @@ namespace euf {
 
     void egraph::set_value(enode* n, lbool value) {        
         force_push();
+        TRACE("euf", tout << bpp(n) << "\n";);
         SASSERT(n->value() == l_undef);
         n->set_value(value);
         m_updates.push_back(update_record(n, update_record::value_assignment()));
@@ -426,12 +427,11 @@ namespace euf {
             set_conflict(n1, n2, j);
             return;
         }
-        if ((r1->class_size() > r2->class_size() && !r2->interpreted()) || r1->interpreted() || r1->value() != l_undef) {
+        if (!r2->interpreted() && 
+             (r1->class_size() > r2->class_size() || r1->interpreted() || r1->value() != l_undef)) {
             std::swap(r1, r2);
             std::swap(n1, n2);
         }
-        if (r1->value() != l_undef)
-            return;
         if (j.is_congruence() && (m.is_false(r2->get_expr()) || m.is_true(r2->get_expr())))
             add_literal(n1, false);
         if (n1->is_equality() && n1->value() == l_false)
diff --git a/src/ast/euf/euf_egraph.h b/src/ast/euf/euf_egraph.h
index d6cd8e430..48362a672 100644
--- a/src/ast/euf/euf_egraph.h
+++ b/src/ast/euf/euf_egraph.h
@@ -214,7 +214,7 @@ namespace euf {
     public:
         egraph(ast_manager& m);
         ~egraph();
-        enode* find(expr* f) { return m_expr2enode.get(f->get_id(), nullptr); }
+        enode* find(expr* f) const { return m_expr2enode.get(f->get_id(), nullptr); }
         enode* mk(expr* f, unsigned n, enode *const* args);
         enode_vector const& enodes_of(func_decl* f);
         void push() { ++m_num_scopes; }
diff --git a/src/ast/static_features.cpp b/src/ast/static_features.cpp
index 7bb34ed14..defabe07b 100644
--- a/src/ast/static_features.cpp
+++ b/src/ast/static_features.cpp
@@ -280,7 +280,7 @@ void static_features::update_core(expr * e) {
     if (is_app(e) && to_app(e)->get_family_id() == m_srfid) 
         m_has_sr = true;
     if (!m_has_arrays && m_arrayutil.is_array(e)) 
-        m_has_arrays = true;
+        check_array(m.get_sort(e));
     if (!m_has_ext_arrays && m_arrayutil.is_array(e) && 
         !m_arrayutil.is_select(e) && !m_arrayutil.is_store(e)) 
         m_has_ext_arrays = true;
@@ -373,6 +373,16 @@ void static_features::update_core(expr * e) {
     }
 }
 
+void static_features::check_array(sort* s) {
+    if (m_arrayutil.is_array(s)) {
+        m_has_arrays = true;
+        update_core(get_array_range(s));
+        for (unsigned i = get_array_arity(s); i-- > 0; )
+            update_core(get_array_domain(s, i));
+    }
+}
+
+
 void static_features::update_core(sort * s) {
     mark_theory(s->get_family_id());
     if (!m_has_int && m_autil.is_int(s))
@@ -383,8 +393,7 @@ void static_features::update_core(sort * s) {
         m_has_bv = true;
     if (!m_has_fpa && (m_fpautil.is_float(s) || m_fpautil.is_rm(s)))
         m_has_fpa = true;
-    if (!m_has_arrays && m_arrayutil.is_array(s))
-        m_has_arrays = true;
+    check_array(s);
 }
 
 void static_features::process(expr * e, bool form_ctx, bool or_and_ctx, bool ite_ctx, unsigned stack_depth) {
diff --git a/src/ast/static_features.h b/src/ast/static_features.h
index cf584081c..8408c08ef 100644
--- a/src/ast/static_features.h
+++ b/src/ast/static_features.h
@@ -134,6 +134,8 @@ struct static_features {
             m_num_theories++; 
         }
     }
+
+    void check_array(sort* s);
     
     void acc_num(rational const & r) {
         if (r.is_neg())
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index 5af5299c2..4e8eb6933 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -321,7 +321,7 @@ public:
     var_index add_named_var(unsigned ext_j, bool is_integer, const std::string&);
     lp_status maximize_term(unsigned j_or_term, impq &term_max);
     inline 
-    core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) { return
+    core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) const { return
             core_solver_pretty_printer<lp::mpq, lp::impq>(m_mpq_lar_core_solver.m_r_solver, out); }
     void get_infeasibility_explanation(explanation &) const;
     inline void backup_x() { m_backup_x = m_mpq_lar_core_solver.m_r_x; }
diff --git a/src/math/lp/lp_api.h b/src/math/lp/lp_api.h
new file mode 100644
index 000000000..dfccd855b
--- /dev/null
+++ b/src/math/lp/lp_api.h
@@ -0,0 +1,131 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Author:
+
+    Lev Nachmanson (levnach)
+    Nikolaj Bjorner (nbjorner)
+
+--*/
+#pragma once
+
+#include "util/inf_rational.h"
+#include "util/optional.h"
+
+namespace lp_api {
+
+    typedef int bool_var;
+    typedef int theory_var;
+
+    enum bound_kind { lower_t, upper_t };
+
+    inline std::ostream& operator<<(std::ostream& out, bound_kind const& k) {
+        switch (k) {
+        case lower_t: return out << "<=";
+        case upper_t: return out << ">=";
+        }
+        return out;
+    }
+
+    class bound {
+        bool_var    m_bv;
+        theory_var  m_var;
+        lp::lpvar        m_vi;
+        bool             m_is_int;
+        rational         m_value;
+        bound_kind       m_bound_kind;
+        lp::constraint_index m_constraints[2];
+
+    public:
+        bound(bool_var bv, theory_var v, lp::lpvar vi, bool is_int, rational const& val, bound_kind k, lp::constraint_index ct, lp::constraint_index cf) :
+            m_bv(bv),
+            m_var(v),
+            m_vi(vi),
+            m_is_int(is_int),
+            m_value(val),
+            m_bound_kind(k) {
+            m_constraints[0] = cf;
+            m_constraints[1] = ct;
+        }
+
+        virtual ~bound() {}
+
+        theory_var get_var() const { return m_var; }
+
+        lp::tv tv() const { return lp::tv::raw(m_vi); }
+
+        bool_var get_bv() const { return m_bv; }
+
+        bound_kind get_bound_kind() const { return m_bound_kind; }
+
+        bool is_int() const { return m_is_int; }
+
+        rational const& get_value() const { return m_value; }
+
+        lp::constraint_index get_constraint(bool b) const { return m_constraints[b]; }
+
+        inf_rational get_value(bool is_true) const {
+            if (is_true) 
+                return inf_rational(m_value);                         // v >= value or v <= value
+            if (m_is_int) {
+                SASSERT(m_value.is_int());
+                rational const& offset = (m_bound_kind == lower_t) ? rational::minus_one() : rational::one();
+                return inf_rational(m_value + offset); // v <= value - 1 or v >= value + 1
+            }
+            else {
+                return inf_rational(m_value, m_bound_kind != lower_t);  // v <= value - epsilon or v >= value + epsilon                                                                        
+            }
+        }
+
+        virtual std::ostream& display(std::ostream& out) const {
+            return out << m_value << "  " << get_bound_kind() << " v" << get_var();
+        }
+    };
+
+    inline std::ostream& operator<<(std::ostream& out, bound const& b) {
+        return b.display(out);
+    }
+
+
+    typedef optional<inf_rational> opt_inf_rational;
+
+
+    struct stats {
+        unsigned m_assert_lower;
+        unsigned m_assert_upper;
+        unsigned m_bounds_propagations;
+        unsigned m_num_iterations;
+        unsigned m_num_iterations_with_no_progress;
+        unsigned m_need_to_solve_inf;
+        unsigned m_fixed_eqs;
+        unsigned m_conflicts;
+        unsigned m_bound_propagations1;
+        unsigned m_bound_propagations2;
+        unsigned m_assert_diseq;
+        unsigned m_gomory_cuts;
+        unsigned m_assume_eqs;
+        unsigned m_branch;
+        stats() { reset(); }
+        void reset() {
+            memset(this, 0, sizeof(*this));
+        }
+        void collect_statistics(statistics& st) const {
+            st.update("arith-lower", m_assert_lower);
+            st.update("arith-upper", m_assert_upper);
+            st.update("arith-propagations", m_bounds_propagations);
+            st.update("arith-iterations", m_num_iterations);
+            st.update("arith-pivots", m_need_to_solve_inf);
+            st.update("arith-plateau-iterations", m_num_iterations_with_no_progress);
+            st.update("arith-fixed-eqs", m_fixed_eqs);
+            st.update("arith-conflicts", m_conflicts);
+            st.update("arith-bound-propagations-lp", m_bound_propagations1);
+            st.update("arith-bound-propagations-cheap", m_bound_propagations2);
+            st.update("arith-diseq", m_assert_diseq);
+            st.update("arith-gomory-cuts", m_gomory_cuts);
+            st.update("arith-assume-eqs", m_assume_eqs);
+            st.update("arith-branch", m_branch);
+        }
+    };
+
+
+}
diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h
index c4659e3e7..05ce8cdb8 100644
--- a/src/math/lp/lp_settings.h
+++ b/src/math/lp/lp_settings.h
@@ -27,6 +27,7 @@ Revision History:
 #include <cstring>
 #include "math/lp/lp_utils.h"
 #include "util/stopwatch.h"
+#include "util/statistics.h"
 #include "math/lp/lp_types.h"
 
 namespace lp {
@@ -126,6 +127,26 @@ struct statistics {
     unsigned m_cheap_eqs;
     statistics() { reset(); }
     void reset() { memset(this, 0, sizeof(*this)); }
+    void collect_statistics(::statistics& st) const {
+        st.update("arith-factorizations", m_num_factorizations);
+        st.update("arith-make-feasible", m_make_feasible);
+        st.update("arith-max-columns", m_max_cols);
+        st.update("arith-max-rows", m_max_rows);
+        st.update("arith-gcd-calls", m_gcd_calls);
+        st.update("arith-gcd-conflict", m_gcd_conflicts);
+        st.update("arith-cube-calls", m_cube_calls);
+        st.update("arith-cube-success", m_cube_success);
+        st.update("arith-patches", m_patches);
+        st.update("arith-patches-success", m_patches_success);
+        st.update("arith-hnf-calls", m_hnf_cutter_calls);
+        st.update("arith-horner-calls", m_horner_calls);
+        st.update("arith-horner-conflicts", m_horner_conflicts);
+        st.update("arith-horner-cross-nested-forms", m_cross_nested_forms);
+        st.update("arith-grobner-calls", m_grobner_calls);
+        st.update("arith-grobner-conflicts", m_grobner_conflicts);
+        st.update("arith-cheap-eqs", m_cheap_eqs);
+
+    }
 };
 
 struct lp_settings {
diff --git a/src/sat/smt/arith_axioms.cpp b/src/sat/smt/arith_axioms.cpp
index ce4e84225..317f3c1a4 100644
--- a/src/sat/smt/arith_axioms.cpp
+++ b/src/sat/smt/arith_axioms.cpp
@@ -25,15 +25,15 @@ namespace arith {
     void solver::mk_div_axiom(expr* p, expr* q) {
         if (a.is_zero(q)) return;
         literal eqz = eq_internalize(q, a.mk_real(0));
-        literal eq  = eq_internalize(a.mk_mul(q, a.mk_div(p, q)), p);
+        literal eq = eq_internalize(a.mk_mul(q, a.mk_div(p, q)), p);
         add_clause(eqz, eq);
     }
 
     // to_int (to_real x) = x
     // to_real(to_int(x)) <= x < to_real(to_int(x)) + 1
     void solver::mk_to_int_axiom(app* n) {
-        expr* x = nullptr, *y = nullptr;
-        VERIFY (a.is_to_int(n, x));            
+        expr* x = nullptr, * y = nullptr;
+        VERIFY(a.is_to_int(n, x));
         if (a.is_to_real(x, y)) {
             literal eq = eq_internalize(y, n);
             add_clause(eq);
@@ -50,7 +50,7 @@ namespace arith {
     }
 
     // is_int(x) <=> to_real(to_int(x)) = x
-    void solver::mk_is_int_axiom(app* n) {
+    void solver::mk_is_int_axiom(expr* n) {
         expr* x = nullptr;
         VERIFY(a.is_is_int(n, x));
         expr_ref lhs(a.mk_to_real(a.mk_to_int(x)), m);
@@ -59,41 +59,7 @@ namespace arith {
         add_equiv(is_int, eq);
     }
 
-
-#if 0
-
-
-    void mk_ite_axiom(expr* n) {
-        return;
-        expr* c = nullptr, *t = nullptr, *e = nullptr;
-        VERIFY(m.is_ite(n, c, t, e));
-        literal e1 = th.mk_eq(n, t, false);
-        literal e2 = th.mk_eq(n, e, false);
-        scoped_trace_stream sts(th, e1, e2);
-        mk_axiom(e1, e2);
-    }
-
-
-    // create axiom for 
-    //    u = v + r*w,
-    ///   abs(r) > r >= 0
-    void assert_idiv_mod_axioms(theory_var u, theory_var v, theory_var w, rational const& r) {
-        app_ref term(m);
-        term = a.mk_mul(a.mk_numeral(r, true), get_enode(w)->get_owner());
-        term = a.mk_add(get_enode(v)->get_owner(), term);
-        term = a.mk_sub(get_enode(u)->get_owner(), term);
-        theory_var z = internalize_def(term);
-        lpvar zi = register_theory_var_in_lar_solver(z);
-        lpvar vi = register_theory_var_in_lar_solver(v);
-        add_def_constraint_and_equality(zi, lp::GE, rational::zero());
-        add_def_constraint_and_equality(zi, lp::LE, rational::zero());
-        add_def_constraint_and_equality(vi, lp::GE, rational::zero());
-        add_def_constraint_and_equality(vi, lp::LT, abs(r));
-        SASSERT(!is_infeasible());
-        TRACE("arith", tout << term << "\n" << lp().constraints(););
-    }
-
-    void mk_idiv_mod_axioms(expr * p, expr * q) {
+    void solver::mk_idiv_mod_axioms(expr* p, expr* q) {
         if (a.is_zero(q)) {
             return;
         }
@@ -111,30 +77,24 @@ namespace arith {
             literal d_le_0 = mk_literal(a.mk_le(div, zero));
             literal m_ge_0 = mk_literal(a.mk_ge(mod, zero));
             literal m_le_0 = mk_literal(a.mk_le(mod, zero));
-            mk_axiom(q_ge_0, d_ge_0);
-            mk_axiom(q_ge_0, d_le_0);
-            mk_axiom(q_ge_0, m_ge_0);
-            mk_axiom(q_ge_0, m_le_0);
-            mk_axiom(q_le_0, d_ge_0);
-            mk_axiom(q_le_0, d_le_0);
-            mk_axiom(q_le_0, m_ge_0);
-            mk_axiom(q_le_0, m_le_0);
+            add_clause(q_ge_0, d_ge_0);
+            add_clause(q_ge_0, d_le_0);
+            add_clause(q_ge_0, m_ge_0);
+            add_clause(q_ge_0, m_le_0);
+            add_clause(q_le_0, d_ge_0);
+            add_clause(q_le_0, d_le_0);
+            add_clause(q_le_0, m_ge_0);
+            add_clause(q_le_0, m_le_0);
             return;
         }
-#if 0
-        expr_ref eqr(m.mk_eq(a.mk_add(a.mk_mul(q, div), mod), p), m);
-        ctx().get_rewriter()(eqr);
-        literal eq = mk_literal(eqr);
-#else
-        literal eq         = th.mk_eq(a.mk_add(a.mk_mul(q, div), mod), p, false);
-#endif
-        literal mod_ge_0   = mk_literal(a.mk_ge(mod, zero));
+        literal eq = eq_internalize(a.mk_add(a.mk_mul(q, div), mod), p);
+        literal mod_ge_0 = mk_literal(a.mk_ge(mod, zero));
 
         rational k(0);
         expr_ref upper(m);
 
         if (a.is_numeral(q, k)) {
-            if (k.is_pos()) { 
+            if (k.is_pos()) {
                 upper = a.mk_numeral(k - 1, true);
             }
             else if (k.is_neg()) {
@@ -145,19 +105,10 @@ namespace arith {
             k = rational::zero();
         }
 
-        context& c = ctx();
         if (!k.is_zero()) {
-            mk_axiom(eq);
-            mk_axiom(mod_ge_0);
-            mk_axiom(mk_literal(a.mk_le(mod, upper)));
-            {
-                std::function<void(void)> log = [&,this]() {
-                    th.log_axiom_unit(m.mk_implies(m.mk_not(m.mk_eq(q, zero)), c.bool_var2expr(eq.var())));
-                    th.log_axiom_unit(m.mk_implies(m.mk_not(m.mk_eq(q, zero)), c.bool_var2expr(mod_ge_0.var())));
-                    th.log_axiom_unit(m.mk_implies(m.mk_not(m.mk_eq(q, zero)), a.mk_le(mod, upper)));
-                };
-                if_trace_stream _ts(m, log);
-            }
+            add_clause(eq);
+            add_clause(mod_ge_0);
+            add_clause(mk_literal(a.mk_le(mod, upper)));
         }
         else {
 
@@ -175,27 +126,286 @@ namespace arith {
             literal q_ge_0 = mk_literal(a.mk_ge(q, zero));
             literal q_le_0 = mk_literal(a.mk_le(q, zero));
 
-            mk_axiom(q_ge_0, eq);
-            mk_axiom(q_le_0, eq);
-            mk_axiom(q_ge_0, mod_ge_0);
-            mk_axiom(q_le_0, mod_ge_0);
-            mk_axiom(q_le_0, ~mk_literal(a.mk_ge(a.mk_sub(mod, q), zero)));            
-            mk_axiom(q_ge_0, ~mk_literal(a.mk_ge(a.mk_add(mod, q), zero)));        
+            add_clause(q_ge_0, eq);
+            add_clause(q_le_0, eq);
+            add_clause(q_ge_0, mod_ge_0);
+            add_clause(q_le_0, mod_ge_0);
+            add_clause(q_le_0, ~mk_literal(a.mk_ge(a.mk_sub(mod, q), zero)));
+            add_clause(q_ge_0, ~mk_literal(a.mk_ge(a.mk_add(mod, q), zero)));
 
         }
-        if (params().m_arith_enum_const_mod && k.is_pos() && k < rational(8)) {
+        if (get_config().m_arith_enum_const_mod && k.is_pos() && k < rational(8)) {
             unsigned _k = k.get_unsigned();
-            literal_buffer lits;
-            expr_ref_vector exprs(m);
+            literal_vector lits;
             for (unsigned j = 0; j < _k; ++j) {
-                literal mod_j = th.mk_eq(mod, a.mk_int(j), false);
+                literal mod_j = eq_internalize(mod, a.mk_int(j));
                 lits.push_back(mod_j);
-                exprs.push_back(c.bool_var2expr(mod_j.var()));
-                ctx().mark_as_relevant(mod_j);
             }
-            ctx().mk_th_axiom(get_id(), lits.size(), lits.begin());                
-        }            
+            add_clause(lits);
+        }
+    }
+
+    //  n < 0 || rem(a, n) =  mod(a, n)
+    // !n < 0 || rem(a, n) = -mod(a, n)
+    void solver::mk_rem_axiom(expr* dividend, expr* divisor) {
+        expr_ref zero(a.mk_int(0), m);
+        expr_ref rem(a.mk_rem(dividend, divisor), m);
+        expr_ref mod(a.mk_mod(dividend, divisor), m);
+        expr_ref mmod(a.mk_uminus(mod), m);
+        expr_ref degz_expr(a.mk_ge(divisor, zero), m);
+        literal dgez = mk_literal(degz_expr);
+        literal pos = eq_internalize(rem, mod);
+        literal neg = eq_internalize(rem, mmod);
+        add_clause(~dgez, pos);
+        add_clause(dgez, neg);
+    }
+
+    void solver::mk_bound_axioms(lp_api::bound& b) {
+        theory_var v = b.get_var();
+        lp_api::bound_kind kind1 = b.get_bound_kind();
+        rational const& k1 = b.get_value();
+        lp_bounds& bounds = m_bounds[v];
+
+        lp_api::bound* end = nullptr;
+        lp_api::bound* lo_inf = end, * lo_sup = end;
+        lp_api::bound* hi_inf = end, * hi_sup = end;
+
+        for (lp_api::bound* other : bounds) {
+            if (other == &b) continue;
+            if (b.get_bv() == other->get_bv()) continue;
+            lp_api::bound_kind kind2 = other->get_bound_kind();
+            rational const& k2 = other->get_value();
+            if (k1 == k2 && kind1 == kind2) {
+                // the bounds are equivalent.
+                continue;
+            }
+
+            SASSERT(k1 != k2 || kind1 != kind2);
+            if (kind2 == lp_api::lower_t) {
+                if (k2 < k1) {
+                    if (lo_inf == end || k2 > lo_inf->get_value()) {
+                        lo_inf = other;
+                    }
+                }
+                else if (lo_sup == end || k2 < lo_sup->get_value()) {
+                    lo_sup = other;
+                }
+            }
+            else if (k2 < k1) {
+                if (hi_inf == end || k2 > hi_inf->get_value()) {
+                    hi_inf = other;
+                }
+            }
+            else if (hi_sup == end || k2 < hi_sup->get_value()) {
+                hi_sup = other;
+            }
+        }
+        if (lo_inf != end) mk_bound_axiom(b, *lo_inf);
+        if (lo_sup != end) mk_bound_axiom(b, *lo_sup);
+        if (hi_inf != end) mk_bound_axiom(b, *hi_inf);
+        if (hi_sup != end) mk_bound_axiom(b, *hi_sup);
+    }
+
+    void solver::mk_bound_axiom(lp_api::bound& b1, lp_api::bound& b2) {
+        literal   l1(b1.get_bv(), false);
+        literal   l2(b2.get_bv(), false);
+        rational const& k1 = b1.get_value();
+        rational const& k2 = b2.get_value();
+        lp_api::bound_kind kind1 = b1.get_bound_kind();
+        lp_api::bound_kind kind2 = b2.get_bound_kind();
+        bool v_is_int = b1.is_int();
+        SASSERT(b1.get_var() == b2.get_var());
+        if (k1 == k2 && kind1 == kind2) return;
+        SASSERT(k1 != k2 || kind1 != kind2);
+
+        if (kind1 == lp_api::lower_t) {
+            if (kind2 == lp_api::lower_t) {
+                if (k2 <= k1)
+                    add_clause(~l1, l2);
+                else
+                    add_clause(l1, ~l2);
+            }
+            else if (k1 <= k2)
+                // k1 <= k2, k1 <= x or x <= k2
+                add_clause(l1, l2);
+            else {
+                // k1 > hi_inf, k1 <= x => ~(x <= hi_inf)
+                add_clause(~l1, ~l2);
+                if (v_is_int && k1 == k2 + rational(1))
+                    // k1 <= x or x <= k1-1
+                    add_clause(l1, l2);
+            }
+        }
+        else if (kind2 == lp_api::lower_t) {
+            if (k1 >= k2)
+                // k1 >= lo_inf, k1 >= x or lo_inf <= x
+                add_clause(l1, l2);
+            else {
+                // k1 < k2, k2 <= x => ~(x <= k1)
+                add_clause(~l1, ~l2);
+                if (v_is_int && k1 == k2 - rational(1))
+                    // x <= k1 or k1+l <= x
+                    add_clause(l1, l2);
+            }
+        }
+        else {
+            // kind1 == A_UPPER, kind2 == A_UPPER
+            if (k1 >= k2)
+                // k1 >= k2, x <= k2 => x <= k1
+                add_clause(l1, ~l2);
+            else
+                // k1 <= hi_sup , x <= k1 =>  x <= hi_sup
+                add_clause(~l1, l2);
+        }
+    }
+
+    lp_api::bound* solver::mk_var_bound(bool_var bv, theory_var v, lp_api::bound_kind bk, rational const& bound) {
+        scoped_internalize_state st(*this);
+        st.vars().push_back(v);
+        st.coeffs().push_back(rational::one());
+        init_left_side(st);
+        lp::constraint_index cT, cF;
+        bool v_is_int = is_int(v);
+        auto vi = register_theory_var_in_lar_solver(v);
+
+        lp::lconstraint_kind kT = bound2constraint_kind(v_is_int, bk, true);
+        lp::lconstraint_kind kF = bound2constraint_kind(v_is_int, bk, false);
+
+        cT = lp().mk_var_bound(vi, kT, bound);
+        if (v_is_int) {
+            rational boundF = (bk == lp_api::lower_t) ? bound - 1 : bound + 1;
+            cF = lp().mk_var_bound(vi, kF, boundF);
+        }
+        else {
+            cF = lp().mk_var_bound(vi, kF, bound);
+        }
+        add_ineq_constraint(cT, literal(bv, false));
+        add_ineq_constraint(cF, literal(bv, true));
+
+        return alloc(lp_api::bound, bv, v, vi, v_is_int, bound, bk, cT, cF);
+    }
+
+    lp::lconstraint_kind solver::bound2constraint_kind(bool is_int, lp_api::bound_kind bk, bool is_true) {
+        switch (bk) {
+        case lp_api::lower_t:
+            return is_true ? lp::GE : (is_int ? lp::LE : lp::LT);
+        case lp_api::upper_t:
+            return is_true ? lp::LE : (is_int ? lp::GE : lp::GT);
+        }
+        UNREACHABLE();
+        return lp::EQ;
+    }
+
+    void solver::mk_eq_axiom(theory_var v1, theory_var v2) {
+        expr* e1 = var2expr(v1);
+        expr* e2 = var2expr(v2);
+        literal le = b_internalize(a.mk_le(e1, e2));
+        literal ge = b_internalize(a.mk_le(e2, e1));
+        literal eq = eq_internalize(e1, e2);
+        add_clause(~eq, le);
+        add_clause(~eq, ge);
+        add_clause(~le, ~ge, eq);
+    }
+
+
+    // create axiom for 
+    //    u = v + r*w,
+    ///   abs(r) > r >= 0
+    void solver::assert_idiv_mod_axioms(theory_var u, theory_var v, theory_var w, rational const& r) {
+        app_ref term(m);
+        term = a.mk_mul(a.mk_numeral(r, true), var2expr(w));
+        term = a.mk_add(var2expr(v), term);
+        term = a.mk_sub(var2expr(u), term);
+        theory_var z = internalize_def(term);
+        lpvar zi = register_theory_var_in_lar_solver(z);
+        lpvar vi = register_theory_var_in_lar_solver(v);
+        add_def_constraint_and_equality(zi, lp::GE, rational::zero());
+        add_def_constraint_and_equality(zi, lp::LE, rational::zero());
+        add_def_constraint_and_equality(vi, lp::GE, rational::zero());
+        add_def_constraint_and_equality(vi, lp::LT, abs(r));
+        SASSERT(!is_infeasible());
+        TRACE("arith", tout << term << "\n" << lp().constraints(););
+    }
+
+    /**
+      * n = (div p q)
+      *
+      * (div p q) * q + (mod p q) = p, (mod p q) >= 0   
+      *
+      * 0 < q => (p/q <= v(p)/v(q) => n <= floor(v(p)/v(q)))
+      * 0 < q => (v(p)/v(q) <= p/q => v(p)/v(q) - 1 < n)
+      *
+      */
+
+    bool solver::check_idiv_bounds() {
+        if (m_idiv_terms.empty()) {
+            return true;
+        }
+        bool all_divs_valid = true;
+        for (unsigned i = 0; i < m_idiv_terms.size(); ++i) {
+            expr* n = m_idiv_terms[i];
+            expr* p = nullptr, * q = nullptr;
+            VERIFY(a.is_idiv(n, p, q));
+            theory_var v = mk_evar(n);
+            theory_var v1 = mk_evar(p);
+
+            if (!can_get_ivalue(v1))
+                continue;
+            lp::impq r1 = get_ivalue(v1);
+            rational r2;
+
+            if (!r1.x.is_int() || r1.x.is_neg() || !r1.y.is_zero()) {
+                // TBD
+                // r1 = 223/4, r2 = 2, r = 219/8 
+                // take ceil(r1), floor(r1), ceil(r2), floor(r2), for floor(r2) > 0
+                // then 
+                //      p/q <= ceil(r1)/floor(r2) => n <= div(ceil(r1), floor(r2))
+                //      p/q >= floor(r1)/ceil(r2) => n >= div(floor(r1), ceil(r2))
+                continue;
+            }
+
+            if (a.is_numeral(q, r2) && r2.is_pos()) {
+                if (!a.is_bounded(n)) {
+                    TRACE("arith", tout << "unbounded " << expr_ref(n, m) << "\n";);
+                    continue;
+                }
+                if (!can_get_ivalue(v))
+                    continue;
+                lp::impq val_v = get_ivalue(v);
+                if (val_v.y.is_zero() && val_v.x == div(r1.x, r2)) continue;
+
+                TRACE("arith", tout << get_value(v) << " != " << r1 << " div " << r2 << "\n";);
+                rational div_r = div(r1.x, r2);
+                // p <= q * div(r1, q) + q - 1 => div(p, q) <= div(r1, r2)
+                // p >= q * div(r1, q) => div(r1, q) <= div(p, q)
+                rational mul(1);
+                rational hi = r2 * div_r + r2 - 1;
+                rational lo = r2 * div_r;
+
+                // used to normalize inequalities so they 
+                // don't appear as 8*x >= 15, but x >= 2
+                expr* n1 = nullptr, * n2 = nullptr;
+                if (a.is_mul(p, n1, n2) && a.is_extended_numeral(n1, mul) && mul.is_pos()) {
+                    p = n2;
+                    hi = floor(hi / mul);
+                    lo = ceil(lo / mul);
+                }
+                literal p_le_r1 = mk_literal(a.mk_le(p, a.mk_numeral(hi, true)));
+                literal p_ge_r1 = mk_literal(a.mk_ge(p, a.mk_numeral(lo, true)));
+                literal n_le_div = mk_literal(a.mk_le(n, a.mk_numeral(div_r, true)));
+                literal n_ge_div = mk_literal(a.mk_ge(n, a.mk_numeral(div_r, true)));
+
+                add_clause(~p_le_r1, n_le_div);
+                add_clause(~p_ge_r1, n_ge_div);
+
+                all_divs_valid = false;
+
+                TRACE("arith", tout << r1 << " div " << r2 << "\n";);
+                continue;
+            }
+        }
+
+        return all_divs_valid;
     }
-#endif
 
 }
diff --git a/src/sat/smt/arith_diagnostics.cpp b/src/sat/smt/arith_diagnostics.cpp
index 70b91dbbd..8a1f2ca31 100644
--- a/src/sat/smt/arith_diagnostics.cpp
+++ b/src/sat/smt/arith_diagnostics.cpp
@@ -22,8 +22,52 @@ Author:
 namespace arith {
 
     
-    std::ostream& solver::display(std::ostream& out) const { return out; }
-    std::ostream& solver::display_justification(std::ostream& out, sat::ext_justification_idx idx) const { return out;}
-    std::ostream& solver::display_constraint(std::ostream& out, sat::ext_constraint_idx idx) const { return out;}
-    void solver::collect_statistics(statistics& st) const {}
+    std::ostream& solver::display(std::ostream& out) const { 
+            out << lp().constraints();
+            lp().print_terms(out);
+            // the tableau
+            lp().pp(out).print();
+            for (unsigned j = 0; j < lp().number_of_vars(); j++) {
+                lp().print_column_info(j, out);
+            }
+        
+        if (m_nla) {
+            m_nla->display(out);
+        }
+        unsigned nv = get_num_vars();
+        for (unsigned v = 0; v < nv; ++v) {
+            auto t = get_tv(v);
+            auto vi = lp().external_to_column_index(v);
+            out << "v" << v << " ";
+            if (t.is_null()) out << "null"; else out << (t.is_term() ? "t" : "j") << vi;
+            if (m_nla && m_nla->use_nra_model() && can_get_ivalue(v)) {
+                scoped_anum an(m_nla->am());
+                m_nla->am().display(out << " = ", nl_value(v, an));
+            }
+            else if (can_get_value(v)) out << " = " << get_value(v);
+            if (is_int(v)) out << ", int";
+            if (ctx.is_shared(var2enode(v))) out << ", shared";
+            out << " := " << mk_bounded_pp(var2expr(v), m) << "\n";
+        }
+        return out; 
+    }
+
+    std::ostream& solver::display_justification(std::ostream& out, sat::ext_justification_idx idx) const { 
+        auto& jst = euf::th_propagation::from_index(idx);
+        for (auto lit : euf::th_propagation::lits(jst))
+            out << lit << " ";
+        for (auto eq : euf::th_propagation::eqs(jst))
+            out << eq.first->get_expr_id() << " == " << eq.second->get_expr_id() << " ";
+        return out;
+    }
+
+    std::ostream& solver::display_constraint(std::ostream& out, sat::ext_constraint_idx idx) const { 
+        return display_justification(out, idx);
+    }
+
+    void solver::collect_statistics(statistics& st) const {
+        m_stats.collect_statistics(st);
+        lp().settings().stats().collect_statistics(st);
+        if (m_nla) m_nla->collect_statistics(st);
+    }
 }
diff --git a/src/sat/smt/arith_internalize.cpp b/src/sat/smt/arith_internalize.cpp
index 5c39a5683..52bcd5bd3 100644
--- a/src/sat/smt/arith_internalize.cpp
+++ b/src/sat/smt/arith_internalize.cpp
@@ -20,14 +20,591 @@ Author:
 
 namespace arith {
 
-    bool solver::visit(expr* e) { return false; }
-    bool solver::visited(expr* e) { return false; }
-    bool solver::post_visit(expr* e, bool sign, bool root) { return false; }
-    void solver::ensure_var(euf::enode* n) {}
-    sat::literal solver::internalize(expr* e, bool sign, bool root, bool learned) { return sat::null_literal; }
-    void solver::internalize(expr* e, bool redundant) {}
-    euf::theory_var solver::mk_var(euf::enode* n) { return euf::null_theory_var; }
-    bool solver::is_shared(theory_var v) const { return false; }   
-    void solver::pop_core(unsigned n) {}
-        
+    sat::literal solver::internalize(expr* e, bool sign, bool root, bool learned) {
+        flet<bool> _is_learned(m_is_redundant, learned);
+        internalize_atom(e);
+        literal lit = ctx.expr2literal(e);
+        if (sign)
+            lit.neg();
+        return lit;
+    }
+
+    void solver::internalize(expr* e, bool redundant) {
+        flet<bool> _is_learned(m_is_redundant, redundant);
+        internalize_term(e);
+    }
+
+    lpvar solver::get_one(bool is_int) {
+        return add_const(1, is_int ? m_one_var : m_rone_var, is_int);
+    }
+
+    lpvar solver::get_zero(bool is_int) {
+        return add_const(0, is_int ? m_zero_var : m_rzero_var, is_int);
+    }
+
+    void solver::ensure_nla() {
+        if (!m_nla) {
+            m_nla = alloc(nla::solver, *m_solver.get(), m.limit());
+            for (auto const& _s : m_scopes) {
+                (void)_s;
+                m_nla->push();
+            }
+            smt_params_helper prms(s().params());
+            m_nla->settings().run_order() = prms.arith_nl_order();
+            m_nla->settings().run_tangents() = prms.arith_nl_tangents();
+            m_nla->settings().run_horner() = prms.arith_nl_horner();
+            m_nla->settings().horner_subs_fixed() = prms.arith_nl_horner_subs_fixed();
+            m_nla->settings().horner_frequency() = prms.arith_nl_horner_frequency();
+            m_nla->settings().horner_row_length_limit() = prms.arith_nl_horner_row_length_limit();
+            m_nla->settings().run_grobner() = prms.arith_nl_grobner();
+            m_nla->settings().run_nra() = prms.arith_nl_nra();
+            m_nla->settings().grobner_subs_fixed() = prms.arith_nl_grobner_subs_fixed();
+            m_nla->settings().grobner_eqs_growth() = prms.arith_nl_grobner_eqs_growth();
+            m_nla->settings().grobner_expr_size_growth() = prms.arith_nl_grobner_expr_size_growth();
+            m_nla->settings().grobner_expr_degree_growth() = prms.arith_nl_grobner_expr_degree_growth();
+            m_nla->settings().grobner_max_simplified() = prms.arith_nl_grobner_max_simplified();
+            m_nla->settings().grobner_number_of_conflicts_to_report() = prms.arith_nl_grobner_cnfl_to_report();
+            m_nla->settings().grobner_quota() = prms.arith_nl_gr_q();
+            m_nla->settings().grobner_frequency() = prms.arith_nl_grobner_frequency();
+            m_nla->settings().expensive_patching() = prms.arith_nl_expp();
+        }
+    }
+
+    void solver::found_unsupported(expr* n) {
+        ctx.push(value_trail<euf::solver, expr*>(m_not_handled));
+        TRACE("arith", tout << "unsupported " << mk_pp(n, m) << "\n";);
+        m_not_handled = n;
+    }
+
+    void solver::found_underspecified(expr* n) {
+        if (a.is_underspecified(n)) {
+            TRACE("arith", tout << "Unhandled: " << mk_pp(n, m) << "\n";);
+            m_underspecified.push_back(to_app(n));
+        }
+        expr* e = nullptr, * x = nullptr, * y = nullptr;
+        if (a.is_div(n, x, y)) {
+            e = a.mk_div0(x, y);
+        }
+        else if (a.is_idiv(n, x, y)) {
+            e = a.mk_idiv0(x, y);
+        }
+        else if (a.is_rem(n, x, y)) {
+            e = a.mk_rem0(x, y);
+        }
+        else if (a.is_mod(n, x, y)) {
+            e = a.mk_mod0(x, y);
+        }
+        else if (a.is_power(n, x, y)) {
+            e = a.mk_power0(x, y);
+        }
+        if (e) {
+            literal lit = eq_internalize(n, e);
+            add_unit(lit);
+        }
+
+    }
+
+    lpvar solver::add_const(int c, lpvar& var, bool is_int) {
+        if (var != UINT_MAX) {
+            return var;
+        }
+        app_ref cnst(a.mk_numeral(rational(c), is_int), m);
+        mk_enode(cnst);
+        theory_var v = mk_evar(cnst);
+        var = lp().add_var(v, is_int);
+        lp().push();
+        add_def_constraint_and_equality(var, lp::GE, rational(c));
+        add_def_constraint_and_equality(var, lp::LE, rational(c));
+        TRACE("arith", tout << "add " << cnst << ", var = " << var << "\n";);
+        return var;
+    }
+
+    lpvar solver::register_theory_var_in_lar_solver(theory_var v) {
+        lpvar lpv = lp().external_to_local(v);
+        if (lpv != lp::null_lpvar)
+            return lpv;
+        return lp().add_var(v, is_int(v));
+    }
+
+    void solver::linearize_term(expr* term, scoped_internalize_state& st) {
+        st.push(term, rational::one());
+        linearize(st);
+    }
+
+    void solver::linearize_ineq(expr* lhs, expr* rhs, scoped_internalize_state& st) {
+        st.push(lhs, rational::one());
+        st.push(rhs, rational::minus_one());
+        linearize(st);
+    }
+
+    void solver::linearize(scoped_internalize_state& st) {
+        expr_ref_vector& terms = st.terms();
+        svector<theory_var>& vars = st.vars();
+        vector<rational>& coeffs = st.coeffs();
+        rational& offset = st.offset();
+        rational r;
+        expr* n1, * n2;
+        unsigned index = 0;
+        while (index < terms.size()) {
+            SASSERT(index >= vars.size());
+            expr* n = terms[index].get();
+            st.to_ensure_enode().push_back(n);
+            if (a.is_add(n)) {
+                for (expr* arg : *to_app(n)) {
+                    st.push(arg, coeffs[index]);
+                }
+                st.set_back(index);
+            }
+            else if (a.is_sub(n)) {
+                unsigned sz = to_app(n)->get_num_args();
+                terms[index] = to_app(n)->get_arg(0);
+                for (unsigned i = 1; i < sz; ++i) {
+                    st.push(to_app(n)->get_arg(i), -coeffs[index]);
+                }
+            }
+            else if (a.is_mul(n, n1, n2) && a.is_extended_numeral(n1, r)) {
+                coeffs[index] *= r;
+                terms[index] = n2;
+                st.to_ensure_enode().push_back(n1);
+            }
+            else if (a.is_mul(n, n1, n2) && a.is_extended_numeral(n2, r)) {
+                coeffs[index] *= r;
+                terms[index] = n1;
+                st.to_ensure_enode().push_back(n2);
+            }
+            else if (a.is_mul(n)) {
+                theory_var v = internalize_mul(to_app(n));
+                coeffs[vars.size()] = coeffs[index];
+                vars.push_back(v);
+                ++index;
+            }
+            else if (a.is_power(n, n1, n2) && is_app(n1) && a.is_extended_numeral(n2, r) && r.is_unsigned() && r <= 10) {
+                theory_var v = internalize_power(to_app(n), to_app(n1), r.get_unsigned());
+                coeffs[vars.size()] = coeffs[index];
+                vars.push_back(v);
+                ++index;
+            }
+            else if (a.is_numeral(n, r)) {
+                offset += coeffs[index] * r;
+                ++index;
+            }
+            else if (a.is_uminus(n, n1)) {
+                coeffs[index].neg();
+                terms[index] = n1;
+            }
+            else if (a.is_to_real(n, n1)) {
+                terms[index] = n1;
+                if (!ctx.get_enode(n)) {
+                    app* t = to_app(n);
+                    VERIFY(internalize_term(to_app(n1)));
+                    mk_enode(t);
+                    theory_var v = mk_evar(n);
+                    theory_var w = mk_evar(n1);
+                    lpvar vj = register_theory_var_in_lar_solver(v);
+                    lpvar wj = register_theory_var_in_lar_solver(w);
+                    auto lu_constraints = lp().add_equality(vj, wj);
+                    add_def_constraint(lu_constraints.first);
+                    add_def_constraint(lu_constraints.second);
+                }
+            }
+            else if (is_app(n) && a.get_family_id() == to_app(n)->get_family_id()) {
+                bool is_first = nullptr == ctx.get_enode(n);
+                app* t = to_app(n);
+                internalize_args(t);
+                mk_enode(t);
+                theory_var v = mk_evar(n);
+                coeffs[vars.size()] = coeffs[index];
+                vars.push_back(v);
+                ++index;
+                if (!is_first) {
+                    // skip recursive internalization
+                }
+                else if (a.is_to_int(n, n1)) {
+                    mk_to_int_axiom(t);
+                }
+                else if (a.is_idiv(n, n1, n2)) {
+                    if (!a.is_numeral(n2, r) || r.is_zero()) found_underspecified(n);
+                    m_idiv_terms.push_back(n);
+                    app_ref mod(a.mk_mod(n1, n2), m);
+                    internalize(mod, m_is_redundant);
+                    st.to_ensure_var().push_back(n1);
+                    st.to_ensure_var().push_back(n2);
+                }
+                else if (a.is_mod(n, n1, n2)) {
+                    if (!a.is_numeral(n2, r) || r.is_zero()) found_underspecified(n);
+                    mk_idiv_mod_axioms(n1, n2);
+                    st.to_ensure_var().push_back(n1);
+                    st.to_ensure_var().push_back(n2);
+                }
+                else if (a.is_rem(n, n1, n2)) {
+                    if (!a.is_numeral(n2, r) || r.is_zero()) found_underspecified(n);
+                    mk_rem_axiom(n1, n2);
+                    st.to_ensure_var().push_back(n1);
+                    st.to_ensure_var().push_back(n2);
+                }
+                else if (a.is_div(n, n1, n2)) {
+                    if (!a.is_numeral(n2, r) || r.is_zero()) found_underspecified(n);
+                    mk_div_axiom(n1, n2);
+                    st.to_ensure_var().push_back(n1);
+                    st.to_ensure_var().push_back(n2);
+                }
+                else if (!a.is_div0(n) && !a.is_mod0(n) && !a.is_idiv0(n) && !a.is_rem0(n)) {
+                    found_unsupported(n);
+                }
+                else {
+                    // no-op
+                }
+            }
+            else {
+                if (is_app(n)) {
+                    internalize_args(to_app(n));
+                }
+                theory_var v = mk_evar(n);
+                coeffs[vars.size()] = coeffs[index];
+                vars.push_back(v);
+                ++index;
+            }
+        }
+        for (unsigned i = st.to_ensure_enode().size(); i-- > 0; ) {
+            expr* n = st.to_ensure_enode()[i];
+            if (is_app(n)) {
+                mk_enode(to_app(n));
+            }
+        }
+        st.to_ensure_enode().reset();
+        for (unsigned i = st.to_ensure_var().size(); i-- > 0; ) {
+            expr* n = st.to_ensure_var()[i];
+            if (is_app(n)) {
+                internalize_term(to_app(n));
+            }
+        }
+        st.to_ensure_var().reset();
+    }
+
+    bool solver::internalize_atom(expr* atom) {
+        TRACE("arith", tout << mk_pp(atom, m) << "\n";);
+        SASSERT(!ctx.get_enode(atom));
+        expr* n1, * n2;
+        rational r;
+        lp_api::bound_kind k;
+        theory_var v = euf::null_theory_var;
+        bool_var bv = ctx.get_si().add_bool_var(atom);
+        m_bool_var2bound.erase(bv);
+        ctx.attach_lit(literal(bv, false), atom);
+
+        if (a.is_le(atom, n1, n2) && a.is_extended_numeral(n2, r) && is_app(n1)) {
+            v = internalize_def(to_app(n1));
+            k = lp_api::upper_t;
+        }
+        else if (a.is_ge(atom, n1, n2) && a.is_extended_numeral(n2, r) && is_app(n1)) {
+            v = internalize_def(to_app(n1));
+            k = lp_api::lower_t;
+        }
+        else if (a.is_le(atom, n1, n2) && a.is_extended_numeral(n1, r) && is_app(n2)) {
+            v = internalize_def(to_app(n2));
+            k = lp_api::lower_t;
+        }
+        else if (a.is_ge(atom, n1, n2) && a.is_extended_numeral(n1, r) && is_app(n2)) {
+            v = internalize_def(to_app(n2));
+            k = lp_api::upper_t;
+        }
+        else if (a.is_is_int(atom)) {
+            mk_is_int_axiom(atom);
+            return true;
+        }
+        else {
+            TRACE("arith", tout << "Could not internalize " << mk_pp(atom, m) << "\n";);
+            found_unsupported(atom);
+            return true;
+        }
+        enode* n = ctx.get_enode(atom);
+        ctx.attach_th_var(n, this, v);
+
+        if (is_int(v) && !r.is_int()) {
+            r = (k == lp_api::upper_t) ? floor(r) : ceil(r);
+        }
+        lp_api::bound* b = mk_var_bound(bv, v, k, r);
+        m_bounds[v].push_back(b);
+        updt_unassigned_bounds(v, +1);
+        m_bounds_trail.push_back(v);
+        m_bool_var2bound.insert(bv, b);
+        TRACE("arith_verbose", tout << "Internalized " << bv << ": " << mk_pp(atom, m) << "\n";);
+        m_new_bounds.push_back(b);
+        //add_use_lists(b);
+        return true;
+    }
+
+
+    bool solver::internalize_term(expr* term) {
+        if (!has_var(term))
+            internalize_def(term);
+        return true;
+    }
+
+    theory_var solver::internalize_def(expr* term, scoped_internalize_state& st) {
+        TRACE("arith", tout << expr_ref(term, m) << "\n";);
+        if (ctx.get_enode(term))
+            return mk_evar(term);
+
+        linearize_term(term, st);
+        if (is_unit_var(st))
+            return st.vars()[0];
+
+        theory_var v = mk_evar(term);
+        SASSERT(euf::null_theory_var != v);
+        st.coeffs().resize(st.vars().size() + 1);
+        st.coeffs()[st.vars().size()] = rational::minus_one();
+        st.vars().push_back(v);
+        return v;
+    }
+
+    // term - v = 0
+    theory_var solver::internalize_def(expr* term) {
+        scoped_internalize_state st(*this);
+        linearize_term(term, st);
+        return internalize_linearized_def(term, st);
+    }
+
+    void solver::internalize_args(app* t, bool force) {
+        if (!force && !reflect(t))
+            return;
+        for (expr* arg : *t)
+            e_internalize(arg);
+    }
+
+    theory_var solver::internalize_power(app* t, app* n, unsigned p) {
+        internalize_args(t, true);
+        bool _has_var = has_var(t);
+        mk_enode(t);
+        theory_var v = mk_evar(t);
+        if (_has_var)
+            return v;
+        theory_var w = mk_evar(n);
+        svector<lpvar> vars;
+        for (unsigned i = 0; i < p; ++i)
+            vars.push_back(register_theory_var_in_lar_solver(w));
+        ensure_nla();
+        m_solver->register_existing_terms();
+        m_nla->add_monic(register_theory_var_in_lar_solver(v), vars.size(), vars.c_ptr());
+        return v;
+    }
+
+    theory_var solver::internalize_mul(app* t) {
+        SASSERT(a.is_mul(t));
+        internalize_args(t, true);
+        bool _has_var = has_var(t);
+        mk_enode(t);
+        theory_var v = mk_evar(t);
+
+        if (!_has_var) {
+            svector<lpvar> vars;
+            for (expr* n : *t) {
+                if (is_app(n)) VERIFY(internalize_term(to_app(n)));
+                SASSERT(ctx.get_enode(n));
+                theory_var v = mk_evar(n);
+                vars.push_back(register_theory_var_in_lar_solver(v));
+            }
+            TRACE("arith", tout << "v" << v << " := " << mk_pp(t, m) << "\n" << vars << "\n";);
+            m_solver->register_existing_terms();
+            ensure_nla();
+            m_nla->add_monic(register_theory_var_in_lar_solver(v), vars.size(), vars.c_ptr());
+        }
+        return v;
+    }
+
+    theory_var solver::internalize_linearized_def(expr* term, scoped_internalize_state& st) {
+        theory_var v = mk_evar(term);
+        TRACE("arith", tout << mk_bounded_pp(term, m) << " v" << v << "\n";);
+
+        if (is_unit_var(st) && v == st.vars()[0]) {
+            return st.vars()[0];
+        }
+        else if (is_one(st) && a.is_numeral(term)) {
+            return get_one(a.is_int(term));
+        }
+        else if (is_zero(st) && a.is_numeral(term)) {
+            return get_zero(a.is_int(term));
+        }
+        else {
+            init_left_side(st);
+            lpvar vi = get_lpvar(v);
+            if (vi == UINT_MAX) {
+                if (m_left_side.empty()) {
+                    vi = lp().add_var(v, a.is_int(term));
+                    add_def_constraint_and_equality(vi, lp::GE, st.offset());
+                    add_def_constraint_and_equality(vi, lp::LE, st.offset());
+                    return v;
+                }
+                if (!st.offset().is_zero()) {
+                    m_left_side.push_back(std::make_pair(st.offset(), get_one(a.is_int(term))));
+                }
+                if (m_left_side.empty()) {
+                    vi = lp().add_var(v, a.is_int(term));
+                    add_def_constraint_and_equality(vi, lp::GE, rational(0));
+                    add_def_constraint_and_equality(vi, lp::LE, rational(0));
+                }
+                else {
+                    vi = lp().add_term(m_left_side, v);
+                    SASSERT(lp::tv::is_term(vi));
+                    TRACE("arith_verbose",
+                        tout << "v" << v << " := " << mk_pp(term, m)
+                        << " slack: " << vi << " scopes: " << m_scopes.size() << "\n";
+                    lp().print_term(lp().get_term(lp::tv::raw(vi)), tout) << "\n";);
+                }
+            }
+            return v;
+        }
+    }
+
+    bool solver::is_unit_var(scoped_internalize_state& st) {
+        return st.offset().is_zero() && st.vars().size() == 1 && st.coeffs()[0].is_one();
+    }
+
+    bool solver::is_one(scoped_internalize_state& st) {
+        return st.offset().is_one() && st.vars().empty();
+    }
+
+    bool solver::is_zero(scoped_internalize_state& st) {
+        return st.offset().is_zero() && st.vars().empty();
+    }
+
+    void solver::init_left_side(scoped_internalize_state& st) {
+        SASSERT(all_zeros(m_columns));
+        svector<theory_var> const& vars = st.vars();
+        vector<rational> const& coeffs = st.coeffs();
+        for (unsigned i = 0; i < vars.size(); ++i) {
+            theory_var var = vars[i];
+            rational const& coeff = coeffs[i];
+            if (m_columns.size() <= static_cast<unsigned>(var)) {
+                m_columns.setx(var, coeff, rational::zero());
+            }
+            else {
+                m_columns[var] += coeff;
+            }
+        }
+        m_left_side.clear();
+        // reset the coefficients after they have been used.
+        for (unsigned i = 0; i < vars.size(); ++i) {
+            theory_var var = vars[i];
+            rational const& r = m_columns[var];
+            if (!r.is_zero()) {
+                m_left_side.push_back(std::make_pair(r, register_theory_var_in_lar_solver(var)));
+                m_columns[var].reset();
+            }
+        }
+        SASSERT(all_zeros(m_columns));
+    }
+
+
+    enode* solver::mk_enode(expr* e) {
+        TRACE("arith", tout << expr_ref(e, m) << "\n";);
+        enode* n = ctx.get_enode(e);
+        if (n)
+            return n;
+        ptr_buffer<enode> args;
+        if (reflect(e))
+            for (expr* arg : *to_app(e))
+                args.push_back(e_internalize(arg));
+        n = ctx.mk_enode(e, args.size(), args.c_ptr());
+        return n;
+    }
+
+    theory_var solver::mk_evar(expr* n) {
+        enode* e = mk_enode(n);
+        if (e->is_attached_to(get_id()))
+            return e->get_th_var(get_id());
+        theory_var v = mk_var(e);
+        TRACE("arith", tout << "fresh var: v" << v << " " << mk_pp(n, m) << "\n";);
+        SASSERT(m_bounds.size() <= static_cast<unsigned>(v) || m_bounds[v].empty());
+        reserve_bounds(v);
+        ctx.attach_th_var(e, this, v);
+        TRACE("arith", tout << mk_pp(n, m) << " " << v << "\n";);
+        SASSERT(euf::null_theory_var != v);
+        return v;
+    }
+
+    bool solver::has_var(expr* e) {
+        enode* n = ctx.get_enode(e);
+        return n && n->is_attached_to(get_id());
+    }
+
+    void solver::add_eq_constraint(lp::constraint_index index, enode* n1, enode* n2) {
+        m_constraint_sources.setx(index, equality_source, null_source);
+        m_equalities.setx(index, enode_pair(n1, n2), enode_pair(nullptr, nullptr));
+    }
+
+    void solver::add_ineq_constraint(lp::constraint_index index, literal lit) {
+        m_constraint_sources.setx(index, inequality_source, null_source);
+        m_inequalities.setx(index, lit, sat::null_literal);
+    }
+
+    void solver::add_def_constraint(lp::constraint_index index) {
+        m_constraint_sources.setx(index, definition_source, null_source);
+        m_definitions.setx(index, euf::null_theory_var, euf::null_theory_var);
+    }
+
+    void solver::add_def_constraint(lp::constraint_index index, theory_var v) {
+        m_constraint_sources.setx(index, definition_source, null_source);
+        m_definitions.setx(index, v, euf::null_theory_var);
+    }
+
+    void solver::add_def_constraint_and_equality(lpvar vi, lp::lconstraint_kind kind,
+        const rational& bound) {
+        lpvar vi_equal;
+        lp::constraint_index ci = lp().add_var_bound_check_on_equal(vi, kind, bound, vi_equal);
+        add_def_constraint(ci);
+        if (vi_equal != lp::null_lpvar)
+            report_equality_of_fixed_vars(vi, vi_equal);
+    }
+
+    bool solver::reflect(expr* n) const {
+        return get_config().m_arith_reflect || a.is_underspecified(n);
+    }
+
+    lpvar solver::get_lpvar(theory_var v) const {
+        return lp().external_to_local(v);
+    }
+
+    lp::tv solver::get_tv(theory_var v) const {
+        return lp::tv::raw(get_lpvar(v));
+    }
+
+    /**
+   \brief We must redefine this method, because theory of arithmetic contains
+   underspecified operators such as division by 0.
+   (/ a b) is essentially an uninterpreted function when b = 0.
+   Thus, 'a' must be considered a shared var if it is the child of an underspecified operator.
+
+   if merge(a / b, x + y) and a / b is root, then x + y become shared and all z + u in equivalence class of x + y.
+
+
+   TBD: when the set of underspecified subterms is small, compute the shared variables below it.
+   Recompute the set if there are merges that invalidate it.
+   Use the set to determine if a variable is shared.
+*/
+    bool solver::is_shared(theory_var v) const {
+        if (m_underspecified.empty()) {
+            return false;
+        }
+        enode* n = var2enode(v);
+        enode* r = n->get_root();
+        unsigned usz = m_underspecified.size();
+        if (r->num_parents() > 2 * usz) {
+            for (unsigned i = 0; i < usz; ++i) {
+                app* u = m_underspecified[i];
+                unsigned sz = u->get_num_args();
+                for (unsigned j = 0; j < sz; ++j)
+                    if (expr2enode(u->get_arg(j))->get_root() == r)
+                        return true;
+            }
+        }
+        else {
+            for (enode* parent : euf::enode_parents(r))
+                if (a.is_underspecified(parent->get_expr()))
+                    return true;
+        }
+        return false;
+    }
+
 }
+
diff --git a/src/sat/smt/arith_solver.cpp b/src/sat/smt/arith_solver.cpp
index 07df6f73e..d72180460 100644
--- a/src/sat/smt/arith_solver.cpp
+++ b/src/sat/smt/arith_solver.cpp
@@ -21,19 +21,1423 @@ Author:
 
 namespace arith {
 
-    solver::solver(euf::solver& ctx, theory_id id):
+    solver::solver(euf::solver& ctx, theory_id id) :
         th_euf_solver(ctx, symbol("arith"), id),
-        a(m)
-    {}
+        m_model_eqs(DEFAULT_HASHTABLE_INITIAL_CAPACITY, var_value_hash(*this), var_value_eq(*this)),
+        m_resource_limit(*this),
+        m_bp(*this),
+        a(m),
+        m_bound_terms(m),
+        m_bound_predicate(m)
+    {
+        reset_variable_values();
+        m_solver = alloc(lp::lar_solver);
+        // initialize 0, 1 variables:
+        get_one(true);
+        get_one(false);
+        get_zero(true);
+        get_zero(false);
+
+        smt_params_helper lpar(ctx.s().params());
+        lp().settings().set_resource_limit(m_resource_limit);
+        lp().settings().simplex_strategy() = static_cast<lp::simplex_strategy_enum>(lpar.arith_simplex_strategy());
+        lp().settings().bound_propagation() = bound_prop_mode::BP_NONE != propagation_mode();
+        lp().settings().enable_hnf() = lpar.arith_enable_hnf();
+        lp().settings().print_external_var_name() = lpar.arith_print_ext_var_names();
+        lp().set_track_pivoted_rows(lpar.arith_bprop_on_pivoted_rows());
+        lp().settings().report_frequency = lpar.arith_rep_freq();
+        lp().settings().print_statistics = lpar.arith_print_stats();
+        lp().settings().cheap_eqs() = lpar.arith_cheap_eqs();
+        lp().set_cut_strategy(get_config().m_arith_branch_cut_ratio);
+        lp().settings().int_run_gcd_test() = get_config().m_arith_gcd_test;
+        lp().settings().set_random_seed(get_config().m_random_seed);
+
+        m_lia = alloc(lp::int_solver, *m_solver.get());
+    }
 
     solver::~solver() {}
 
-    void solver::asserted(literal l) {}
-    sat::check_result solver::check() { return sat::check_result::CR_DONE; }
-    
-    euf::th_solver* solver::clone(sat::solver* s, euf::solver& ctx) { return nullptr;}
-    void solver::new_eq_eh(euf::th_eq const& eq) {}
-    bool solver::unit_propagate() { return false; }
-    void solver::add_value(euf::enode* n, model& mdl, expr_ref_vector& values) {}
-    void solver::add_dep(euf::enode* n, top_sort<euf::enode>& dep) {}
+    void solver::asserted(literal l) {
+        m_asserted.push_back(l);
+    }
+
+    euf::th_solver* solver::clone(sat::solver* s, euf::solver& ctx) {
+        arith::solver* result = alloc(arith::solver, ctx, get_id());
+        ast_translation tr(m, ctx.get_manager());
+        for (unsigned i = result->get_num_vars(); i < get_num_vars(); ++i) {
+            expr* e1 = var2expr(i);
+            expr* e2 = tr(e1);
+            result->mk_evar(e2);
+        }
+
+        unsigned v = 0;
+        result->m_bounds.resize(m_bounds.size());
+        for (auto const& bounds : m_bounds) {            
+            for (auto* b : bounds) {
+                auto* b2 = result->mk_var_bound(b->get_bv(), v, b->get_bound_kind(), b->get_value());
+                result->m_bounds[v].push_back(b2);
+                result->m_bounds_trail.push_back(v);
+                result->updt_unassigned_bounds(v, +1);
+                result->m_bool_var2bound.insert(b->get_bv(), b2);
+                result->m_new_bounds.push_back(b2);
+            }
+            ++v;
+        }
+
+        // clone rows into m_solver, m_nla, m_lia
+        NOT_IMPLEMENTED_YET();
+
+        return result;        
+    }
+
+    bool solver::unit_propagate() {
+        if (m_new_bounds.empty() && m_asserted_qhead == m_asserted.size())
+            return false;
+
+        flush_bound_axioms();
+
+        unsigned qhead = m_asserted_qhead;
+        while (m_asserted_qhead < m_asserted.size() && !s().inconsistent() && m.inc()) {
+            literal lit = m_asserted[m_asserted_qhead];
+            lp_api::bound* b = nullptr;
+            TRACE("arith", tout << "propagate: " << lit << "\n";);
+            CTRACE("arith", !m_bool_var2bound.contains(lit.var()), tout << "not found " << lit << "\n";);
+            if (m_bool_var2bound.find(lit.var(), b)) 
+                assert_bound(!lit.sign(), *b);
+            ++m_asserted_qhead;
+        }
+        if (s().inconsistent())
+            return true;
+
+        lbool lbl = make_feasible();
+        if (!m.inc())
+            return false;
+
+        switch (lbl) {
+        case l_false:
+            TRACE("arith", tout << "propagation conflict\n";);
+            get_infeasibility_explanation_and_set_conflict();
+            break;
+        case l_true:
+            propagate_basic_bounds(qhead);
+            propagate_bounds_with_lp_solver();
+            break;
+        case l_undef:
+            break;
+        }
+        return true;
+    }
+
+    void solver::propagate_basic_bounds(unsigned qhead) {
+        lp_api::bound* b = nullptr;
+        for (; qhead < m_asserted_qhead && !s().inconsistent(); ++qhead) {
+            literal lit = m_asserted[qhead];
+            if (m_bool_var2bound.find(lit.var(), b)) 
+                propagate_bound(lit, *b);            
+        }
+    }
+
+    // for glb lo': lo' < lo:
+    //   lo <= x -> lo' <= x 
+    //   lo <= x -> ~(x <= lo')
+    // for lub hi': hi' > hi
+    //   x <= hi -> x <= hi'
+    //   x <= hi -> ~(x >= hi')
+
+    void solver::propagate_bound(literal lit1, lp_api::bound& b) {
+        if (bound_prop_mode::BP_NONE == propagation_mode())
+            return;
+
+        bool is_true = !lit1.sign();
+        lp_api::bound_kind k = b.get_bound_kind();
+        theory_var v = b.get_var();
+        inf_rational val = b.get_value(is_true);
+        lp_bounds const& bounds = m_bounds[v];
+        SASSERT(!bounds.empty());
+        if (bounds.size() == 1) return;
+        if (m_unassigned_bounds[v] == 0) return;
+        bool v_is_int = b.is_int();
+        literal lit2 = sat::null_literal;
+        bool find_glb = (is_true == (k == lp_api::lower_t));
+        TRACE("arith", tout << "v" << v << " find_glb: " << find_glb << " is_true: " << is_true << " k: " << k << " is_lower: " << (k == lp_api::lower_t) << "\n";);
+        if (find_glb) {
+            rational glb;
+            lp_api::bound* lb = nullptr;
+            for (lp_api::bound* b2 : bounds) {
+                if (b2 == &b) continue;
+                rational const& val2 = b2->get_value();
+                if (((is_true || v_is_int) ? val2 < val : val2 <= val) && (!lb || glb < val2)) {
+                    lb = b2;
+                    glb = val2;
+                }
+            }
+            if (!lb) return;
+            bool sign = lb->get_bound_kind() != lp_api::lower_t;
+            lit2 = literal(lb->get_bv(), sign);
+        }
+        else {
+            rational lub;
+            lp_api::bound* ub = nullptr;
+            for (lp_api::bound* b2 : bounds) {
+                if (b2 == &b) continue;
+                rational const& val2 = b2->get_value();
+                if (((is_true || v_is_int) ? val < val2 : val <= val2) && (!ub || val2 < lub)) {
+                    ub = b2;
+                    lub = val2;
+                }
+            }
+            if (!ub) return;
+            bool sign = ub->get_bound_kind() != lp_api::upper_t;
+            lit2 = literal(ub->get_bv(), sign);
+        }
+        updt_unassigned_bounds(v, -1);
+        ++m_stats.m_bound_propagations2;
+        reset_evidence();
+        m_core.push_back(lit1);
+        m_params.push_back(parameter(m_farkas));
+        m_params.push_back(parameter(rational(1)));
+        m_params.push_back(parameter(rational(1)));
+        assign(lit2, m_core, m_eqs, m_params);
+        ++m_stats.m_bounds_propagations;
+    }
+
+    void solver::propagate_bounds_with_lp_solver() {
+        if (!should_propagate())
+            return;
+
+        m_bp.init();
+        lp().propagate_bounds_for_touched_rows(m_bp);
+
+        if (!m.inc())
+            return;
+
+        if (is_infeasible())
+            get_infeasibility_explanation_and_set_conflict();
+        else
+            for (auto& ib : m_bp.ibounds())
+                if (m.inc() && !s().inconsistent())
+                    propagate_lp_solver_bound(ib);
+    }
+
+    void solver::propagate_lp_solver_bound(const lp::implied_bound& be) {
+        lpvar vi = be.m_j;
+        theory_var v = lp().local_to_external(vi);
+
+        if (v == euf::null_theory_var)
+            return;
+
+        TRACE("arith", tout << "v" << v << " " << be.kind() << " " << be.m_bound << "\n";);
+
+        reserve_bounds(v);
+
+        if (m_unassigned_bounds[v] == 0 && !should_refine_bounds()) {
+            TRACE("arith", tout << "return\n";);
+            return;
+        }
+        lp_bounds const& bounds = m_bounds[v];
+        bool first = true;
+        for (unsigned i = 0; i < bounds.size(); ++i) {
+            lp_api::bound* b = bounds[i];
+            if (s().value(b->get_bv()) != l_undef) {
+                continue;
+            }
+            literal lit = is_bound_implied(be.kind(), be.m_bound, *b);
+            if (lit == sat::null_literal) {
+                continue;
+            }
+            TRACE("arith", tout << lit << " bound: " << *b << " first: " << first << "\n";);
+
+            lp().settings().stats().m_num_of_implied_bounds++;
+            if (first) {
+                first = false;
+                reset_evidence();
+                m_explanation.clear();
+                lp().explain_implied_bound(be, m_bp);
+            }
+            CTRACE("arith", m_unassigned_bounds[v] == 0, tout << "missed bound\n";);
+            updt_unassigned_bounds(v, -1);
+            DEBUG_CODE(for (auto& lit : m_core) { VERIFY(s().value(lit) == l_true); });
+            ++m_stats.m_bound_propagations1;
+            assign(lit, m_core, m_eqs, m_params);
+        }
+
+        if (should_refine_bounds() && first)
+            refine_bound(v, be);
+    }
+
+    literal solver::is_bound_implied(lp::lconstraint_kind k, rational const& value, lp_api::bound const& b) const {
+        if ((k == lp::LE || k == lp::LT) && b.get_bound_kind() == lp_api::upper_t && value <= b.get_value()) {
+            TRACE("arith", tout << "v <= value <= b.get_value() => v <= b.get_value() \n";);
+            return literal(b.get_bv(), false);
+        }
+        if ((k == lp::GE || k == lp::GT) && b.get_bound_kind() == lp_api::lower_t && b.get_value() <= value) {
+            TRACE("arith", tout << "b.get_value() <= value <= v => b.get_value() <= v \n";);
+            return literal(b.get_bv(), false);
+        }
+        if (k == lp::LE && b.get_bound_kind() == lp_api::lower_t && value < b.get_value()) {
+            TRACE("arith", tout << "v <= value < b.get_value() => v < b.get_value()\n";);
+            return literal(b.get_bv(), true);
+        }
+        if (k == lp::LT && b.get_bound_kind() == lp_api::lower_t && value <= b.get_value()) {
+            TRACE("arith", tout << "v < value <= b.get_value() => v < b.get_value()\n";);
+            return literal(b.get_bv(), true);
+        }
+        if (k == lp::GE && b.get_bound_kind() == lp_api::upper_t && b.get_value() < value) {
+            TRACE("arith", tout << "b.get_value() < value <= v => b.get_value() < v\n";);
+            return literal(b.get_bv(), true);
+        }
+        if (k == lp::GT && b.get_bound_kind() == lp_api::upper_t && b.get_value() <= value) {
+            TRACE("arith", tout << "b.get_value() <= value < v => b.get_value() < v\n";);
+            return literal(b.get_bv(), true);
+        }
+        return sat::null_literal;
+    }
+
+
+    void solver::consume(rational const& v, lp::constraint_index j) {
+        set_evidence(j, m_core, m_eqs);
+        m_explanation.add_pair(j, v);
+    }
+
+    void solver::add_eq(lpvar u, lpvar v, lp::explanation const& e) {
+        if (s().inconsistent())
+            return;
+        theory_var uv = lp().local_to_external(u); // variables that are returned should have external representations
+        theory_var vv = lp().local_to_external(v); // so maybe better to have them already transformed to external form
+        if (is_equal(uv, vv))
+            return;
+        enode* n1 = var2enode(uv);
+        enode* n2 = var2enode(vv);
+        if (m.get_sort(n1->get_expr()) != m.get_sort(n2->get_expr()))
+            return;
+        reset_evidence();
+        for (auto const& ev : e)
+            set_evidence(ev.ci(), m_core, m_eqs);
+        auto* jst = euf::th_propagation::mk(*this, m_core, m_eqs);
+        ctx.propagate(n1, n2, jst->to_index());
+    }
+
+    bool solver::bound_is_interesting(unsigned vi, lp::lconstraint_kind kind, const rational& bval) const {
+        theory_var v = lp().local_to_external(vi);
+        if (v == euf::null_theory_var)
+            return false;
+
+        if (should_refine_bounds())
+            return true;
+
+        if (m_bounds.size() <= static_cast<unsigned>(v) || m_unassigned_bounds[v] == 0)
+            return false;
+
+        for (lp_api::bound* b : m_bounds[v])
+            if (s().value(b->get_bv()) == l_undef && sat::null_literal != is_bound_implied(kind, bval, *b))
+                return true;
+
+        return false;
+    }
+
+    void solver::refine_bound(theory_var v, const lp::implied_bound& be) {
+        lpvar vi = be.m_j;
+        if (lp::tv::is_term(vi))
+            return;
+        expr_ref w(var2expr(v), m);
+        if (a.is_add(w) || a.is_numeral(w) || m.is_ite(w))
+            return;
+        literal bound = sat::null_literal;
+        switch (be.kind()) {
+        case lp::LE:
+            if (is_int(v) && (lp().column_has_lower_bound(vi) || !lp().column_has_upper_bound(vi)))
+                bound = mk_literal(a.mk_le(w, a.mk_numeral(floor(be.m_bound), a.is_int(w))));
+            if (is_real(v) && !lp().column_has_upper_bound(vi))
+                bound = mk_literal(a.mk_le(w, a.mk_numeral(be.m_bound, a.is_int(w))));
+            break;
+        case lp::GE:
+            if (is_int(v) && (lp().column_has_upper_bound(vi) || !lp().column_has_lower_bound(vi)))
+                bound = mk_literal(a.mk_ge(w, a.mk_numeral(ceil(be.m_bound), a.is_int(w))));
+            if (is_real(v) && !lp().column_has_lower_bound(vi))
+                bound = mk_literal(a.mk_ge(w, a.mk_numeral(be.m_bound, a.is_int(w))));
+            break;
+        default:
+            break;
+        }
+        if (bound == sat::null_literal)
+            return;
+        if (s().value(bound) == l_true)
+            return;
+
+        ++m_stats.m_bound_propagations1;
+        reset_evidence();
+        m_explanation.clear();
+        lp().explain_implied_bound(be, m_bp);
+        assign(bound, m_core, m_eqs, m_params);
+    }
+
+
+    void solver::assert_bound(bool is_true, lp_api::bound& b) {
+        TRACE("arith", tout << b << "\n";);
+        lp::constraint_index ci = b.get_constraint(is_true);
+        lp().activate(ci);
+        if (is_infeasible()) {
+            return;
+        }
+        lp::lconstraint_kind k = bound2constraint_kind(b.is_int(), b.get_bound_kind(), is_true);
+        if (k == lp::LT || k == lp::LE) {
+            ++m_stats.m_assert_lower;
+        }
+        else {
+            ++m_stats.m_assert_upper;
+        }
+        inf_rational value = b.get_value(is_true);
+        if (propagate_eqs() && value.is_rational()) {
+            propagate_eqs(b.tv(), ci, k, b, value.get_rational());
+        }
+#if 0
+        if (propagation_mode() != BP_NONE)
+            lp().mark_rows_for_bound_prop(b.tv().id());
+#endif
+
+    }
+
+    void solver::propagate_eqs(lp::tv t, lp::constraint_index ci, lp::lconstraint_kind k, lp_api::bound& b, rational const& value) {
+        if (k == lp::GE && set_lower_bound(t, ci, value) && has_upper_bound(t.index(), ci, value)) {
+            fixed_var_eh(b.get_var(), value);
+        }
+        else if (k == lp::LE && set_upper_bound(t, ci, value) && has_lower_bound(t.index(), ci, value)) {
+            fixed_var_eh(b.get_var(), value);
+        }
+    }
+
+
+    bool solver::set_bound(lp::tv tv, lp::constraint_index ci, rational const& v, bool is_lower) {
+        if (tv.is_term()) {
+            lpvar ti = tv.id();
+            auto& vec = is_lower ? m_lower_terms : m_upper_terms;
+            if (vec.size() <= ti) {
+                vec.resize(ti + 1, constraint_bound(UINT_MAX, rational()));
+            }
+            constraint_bound& b = vec[ti];
+            if (b.first == UINT_MAX || (is_lower ? b.second < v : b.second > v)) {
+                TRACE("arith", tout << "tighter bound " << tv.to_string() << "\n";);
+                m_history.push_back(vec[ti]);
+                ctx.push(history_trail<euf::solver, constraint_bound>(vec, ti, m_history));
+                b.first = ci;
+                b.second = v;
+            }
+            return true;
+        }
+        else {
+            TRACE("arith", tout << "not a term " << tv.to_string() << "\n";);
+            // m_solver already tracks bounds on proper variables, but not on terms.
+            bool is_strict = false;
+            rational b;
+            if (is_lower) {
+                return lp().has_lower_bound(tv.id(), ci, b, is_strict) && !is_strict && b == v;
+            }
+            else {
+                return lp().has_upper_bound(tv.id(), ci, b, is_strict) && !is_strict && b == v;
+            }
+        }
+    }
+
+    void solver::flush_bound_axioms() {
+        CTRACE("arith", !m_new_bounds.empty(), tout << "flush bound axioms\n";);
+
+        while (!m_new_bounds.empty()) {
+            lp_bounds atoms;
+            atoms.push_back(m_new_bounds.back());
+            m_new_bounds.pop_back();
+            theory_var v = atoms.back()->get_var();
+            for (unsigned i = 0; i < m_new_bounds.size(); ++i) {
+                if (m_new_bounds[i]->get_var() == v) {
+                    atoms.push_back(m_new_bounds[i]);
+                    m_new_bounds[i] = m_new_bounds.back();
+                    m_new_bounds.pop_back();
+                    --i;
+                }
+            }
+            CTRACE("arith_verbose", !atoms.empty(),
+                for (unsigned i = 0; i < atoms.size(); ++i) {
+                    atoms[i]->display(tout); tout << "\n";
+                });
+            lp_bounds occs(m_bounds[v]);
+
+            std::sort(atoms.begin(), atoms.end(), compare_bounds());
+            std::sort(occs.begin(), occs.end(), compare_bounds());
+
+            iterator begin1 = occs.begin();
+            iterator begin2 = occs.begin();
+            iterator end = occs.end();
+            begin1 = first(lp_api::lower_t, begin1, end);
+            begin2 = first(lp_api::upper_t, begin2, end);
+
+            iterator lo_inf = begin1, lo_sup = begin1;
+            iterator hi_inf = begin2, hi_sup = begin2;
+            bool flo_inf, fhi_inf, flo_sup, fhi_sup;
+            ptr_addr_hashtable<lp_api::bound> visited;
+            for (unsigned i = 0; i < atoms.size(); ++i) {
+                lp_api::bound* a1 = atoms[i];
+                iterator lo_inf1 = next_inf(a1, lp_api::lower_t, lo_inf, end, flo_inf);
+                iterator hi_inf1 = next_inf(a1, lp_api::upper_t, hi_inf, end, fhi_inf);
+                iterator lo_sup1 = next_sup(a1, lp_api::lower_t, lo_sup, end, flo_sup);
+                iterator hi_sup1 = next_sup(a1, lp_api::upper_t, hi_sup, end, fhi_sup);
+                if (lo_inf1 != end) lo_inf = lo_inf1;
+                if (lo_sup1 != end) lo_sup = lo_sup1;
+                if (hi_inf1 != end) hi_inf = hi_inf1;
+                if (hi_sup1 != end) hi_sup = hi_sup1;
+                if (!flo_inf) lo_inf = end;
+                if (!fhi_inf) hi_inf = end;
+                if (!flo_sup) lo_sup = end;
+                if (!fhi_sup) hi_sup = end;
+                visited.insert(a1);
+                if (lo_inf1 != end && lo_inf != end && !visited.contains(*lo_inf)) mk_bound_axiom(*a1, **lo_inf);
+                if (lo_sup1 != end && lo_sup != end && !visited.contains(*lo_sup)) mk_bound_axiom(*a1, **lo_sup);
+                if (hi_inf1 != end && hi_inf != end && !visited.contains(*hi_inf)) mk_bound_axiom(*a1, **hi_inf);
+                if (hi_sup1 != end && hi_sup != end && !visited.contains(*hi_sup)) mk_bound_axiom(*a1, **hi_sup);
+            }
+        }
+    }
+
+    lp_bounds::iterator solver::first(
+        lp_api::bound_kind kind,
+        iterator it,
+        iterator end) {
+        for (; it != end; ++it) {
+            lp_api::bound* a = *it;
+            if (a->get_bound_kind() == kind) return it;
+        }
+        return end;
+    }
+
+    lp_bounds::iterator solver::next_inf(
+        lp_api::bound* a1,
+        lp_api::bound_kind kind,
+        iterator it,
+        iterator end,
+        bool& found_compatible) {
+        rational const& k1(a1->get_value());
+        iterator result = end;
+        found_compatible = false;
+        for (; it != end; ++it) {
+            lp_api::bound* a2 = *it;
+            if (a1 == a2) continue;
+            if (a2->get_bound_kind() != kind) continue;
+            rational const& k2(a2->get_value());
+            found_compatible = true;
+            if (k2 <= k1) {
+                result = it;
+            }
+            else {
+                break;
+            }
+        }
+        return result;
+    }
+
+    lp_bounds::iterator solver::next_sup(
+        lp_api::bound* a1,
+        lp_api::bound_kind kind,
+        iterator it,
+        iterator end,
+        bool& found_compatible) {
+        rational const& k1(a1->get_value());
+        found_compatible = false;
+        for (; it != end; ++it) {
+            lp_api::bound* a2 = *it;
+            if (a1 == a2) continue;
+            if (a2->get_bound_kind() != kind) continue;
+            rational const& k2(a2->get_value());
+            found_compatible = true;
+            if (k1 < k2) {
+                return it;
+            }
+        }
+        return end;
+    }
+
+    void solver::add_value(euf::enode* n, model& mdl, expr_ref_vector& values) {
+        theory_var v = n->get_th_var(get_id());
+        expr* o = n->get_expr();
+        expr_ref value(m);
+        if (use_nra_model() && lp().external_to_local(v) != lp::null_lpvar) {
+            anum const& an = nl_value(v, *m_a1);
+            if (a.is_int(o) && !m_nla->am().is_int(an))
+                value = a.mk_numeral(rational::zero(), a.is_int(o));
+            else
+                value = a.mk_numeral(m_nla->am(), nl_value(v, *m_a1), a.is_int(o));
+        }
+        else {
+            rational r = get_value(v);
+            TRACE("arith", tout << mk_pp(o, m) << " v" << v << " := " << r << "\n";);
+            SASSERT("integer variables should have integer values: " && (!a.is_int(o) || r.is_int() || m.limit().is_canceled()));
+            if (a.is_int(o) && !r.is_int()) 
+                r = floor(r);
+            value = a.mk_numeral(r, m.get_sort(o));
+        }
+        values.set(n->get_root_id(), value);
+    }
+
+    void solver::push_core() {
+        TRACE("arith", tout << "push\n";);
+        m_scopes.push_back(scope());
+        scope& sc = m_scopes.back();
+        sc.m_bounds_lim = m_bounds_trail.size();
+        sc.m_asserted_qhead = m_asserted_qhead;
+        sc.m_idiv_lim = m_idiv_terms.size();
+        sc.m_asserted_lim = m_asserted.size();
+        sc.m_not_handled = m_not_handled;
+        sc.m_underspecified_lim = m_underspecified.size();
+        lp().push();
+        if (m_nla)
+            m_nla->push();
+
+    }
+
+    void solver::pop_core(unsigned num_scopes) {
+        TRACE("arith", tout << "pop " << num_scopes << "\n";);
+        unsigned old_size = m_scopes.size() - num_scopes;
+        del_bounds(m_scopes[old_size].m_bounds_lim);
+        m_idiv_terms.shrink(m_scopes[old_size].m_idiv_lim);
+        m_asserted.shrink(m_scopes[old_size].m_asserted_lim);
+        m_asserted_qhead = m_scopes[old_size].m_asserted_qhead;
+        m_underspecified.shrink(m_scopes[old_size].m_underspecified_lim);
+        m_not_handled = m_scopes[old_size].m_not_handled;
+        m_scopes.resize(old_size);
+        lp().pop(num_scopes);
+        m_new_bounds.reset();
+        if (m_nla)
+            m_nla->pop(num_scopes);
+        TRACE("arith", tout << "num scopes: " << num_scopes << " new scope level: " << m_scopes.size() << "\n";);
+    }
+
+    void solver::del_bounds(unsigned old_size) {
+        for (unsigned i = m_bounds_trail.size(); i-- > old_size; ) {
+            unsigned v = m_bounds_trail[i];
+            lp_api::bound* b = m_bounds[v].back();
+            // del_use_lists(b);
+            dealloc(b);
+            m_bounds[v].pop_back();
+        }
+        m_bounds_trail.shrink(old_size);
+    }
+
+    void solver::report_equality_of_fixed_vars(unsigned vi1, unsigned vi2) {
+        rational bound;
+        lp::constraint_index ci1, ci2, ci3, ci4;
+        theory_var v1 = lp().local_to_external(vi1);
+        theory_var v2 = lp().local_to_external(vi2);
+        TRACE("arith", tout << "fixed: " << mk_pp(var2expr(v1), m) << " " << mk_pp(var2expr(v2), m) << "\n";);
+        // we expect lp() to ensure that none of these returns happen.
+
+        if (is_equal(v1, v2))
+            return;
+        if (is_int(v1) != is_int(v2))
+            return;
+        if (!has_lower_bound(vi1, ci1, bound))
+            return;
+        if (!has_upper_bound(vi1, ci2, bound))
+            return;
+        if (!has_lower_bound(vi2, ci3, bound))
+            return;
+        if (!has_upper_bound(vi2, ci4, bound))
+            return;
+
+        ++m_stats.m_fixed_eqs;
+        reset_evidence();
+        set_evidence(ci1, m_core, m_eqs);
+        set_evidence(ci2, m_core, m_eqs);
+        set_evidence(ci3, m_core, m_eqs);
+        set_evidence(ci4, m_core, m_eqs);
+        enode* x = var2enode(v1);
+        enode* y = var2enode(v2);
+        auto* jst = euf::th_propagation::mk(*this, m_core, m_eqs);
+        ctx.propagate(x, y, jst->to_index());
+    }
+
+    bool solver::is_equal(theory_var x, theory_var y) const {
+        return x == y || var2enode(x)->get_root() == var2enode(y)->get_root();
+    }
+
+    bool solver::has_upper_bound(lpvar vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, false); }
+
+    bool solver::has_lower_bound(lpvar vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, true); }
+
+    bool solver::has_bound(lpvar vi, lp::constraint_index& ci, rational const& bound, bool is_lower) {
+        if (lp::tv::is_term(vi)) {
+            theory_var v = lp().local_to_external(vi);
+            rational val;
+            TRACE("arith", tout << lp().get_variable_name(vi) << " " << v << "\n";);
+            if (v != euf::null_theory_var && a.is_numeral(var2expr(v), val) && bound == val) {
+                ci = UINT_MAX;
+                return bound == val;
+            }
+
+            auto& vec = is_lower ? m_lower_terms : m_upper_terms;
+            lpvar ti = lp::tv::unmask_term(vi);
+            if (vec.size() > ti) {
+                constraint_bound& b = vec[ti];
+                ci = b.first;
+                return ci != UINT_MAX && bound == b.second;
+            }
+            else {
+                return false;
+            }
+        }
+        else {
+            bool is_strict = false;
+            rational b;
+            if (is_lower) {
+                return lp().has_lower_bound(vi, ci, b, is_strict) && b == bound && !is_strict;
+            }
+            else {
+                return lp().has_upper_bound(vi, ci, b, is_strict) && b == bound && !is_strict;
+            }
+        }
+    }
+
+    void solver::updt_unassigned_bounds(theory_var v, int inc) {
+        TRACE("arith", tout << "v" << v << " " << m_unassigned_bounds[v] << " += " << inc << "\n";);
+        ctx.push(vector_value_trail<euf::solver, unsigned, false>(m_unassigned_bounds, v));
+        m_unassigned_bounds[v] += inc;
+    }
+
+    void solver::reserve_bounds(theory_var v) {
+        while (m_bounds.size() <= static_cast<unsigned>(v)) {
+            m_bounds.push_back(lp_bounds());
+            m_unassigned_bounds.push_back(0);
+        }
+    }
+
+    bool solver::all_zeros(vector<rational> const& v) const {
+        for (rational const& r : v)
+            if (!r.is_zero())
+                return false;
+        return true;
+    }
+
+    bound_prop_mode solver::propagation_mode() const {
+        return m_num_conflicts < get_config().m_arith_propagation_threshold ?
+            get_config().m_arith_bound_prop :
+            bound_prop_mode::BP_NONE;
+    }
+
+    void solver::init_variable_values() {
+        reset_variable_values();
+        if (m.inc() && m_solver.get() && get_num_vars() > 0) {
+            TRACE("arith", display(tout << "update variable values\n"););
+            lp().get_model(m_variable_values);
+        }
+    }
+
+    void solver::reset_variable_values() {
+        m_variable_values.clear();
+    }
+
+    lbool solver::get_phase(bool_var v) {
+        lp_api::bound* b;
+        if (!m_bool_var2bound.find(v, b)) {
+            return l_undef;
+        }
+        lp::lconstraint_kind k = lp::EQ;
+        switch (b->get_bound_kind()) {
+        case lp_api::lower_t:
+            k = lp::GE;
+            break;
+        case lp_api::upper_t:
+            k = lp::LE;
+            break;
+        default:
+            break;
+        }
+        auto vi = register_theory_var_in_lar_solver(b->get_var());
+        if (vi == lp::null_lpvar) {
+            return l_undef;
+        }
+        return lp().compare_values(vi, k, b->get_value()) ? l_true : l_false;
+    }
+
+    bool solver::can_get_value(theory_var v) const {
+        return can_get_bound(v) && !m_variable_values.empty();
+    }
+
+    bool solver::can_get_bound(theory_var v) const {
+        return v != euf::null_theory_var && lp().external_is_used(v);
+    }
+
+    bool solver::can_get_ivalue(theory_var v) const {
+        return can_get_bound(v);
+    }
+
+    void solver::ensure_column(theory_var v) {
+        if (!lp().external_is_used(v))
+            register_theory_var_in_lar_solver(v);
+    }
+
+    lp::impq solver::get_ivalue(theory_var v) const {
+        SASSERT(can_get_ivalue(v));
+        auto t = get_tv(v);
+        if (!t.is_term())
+            return lp().get_column_value(t.id());
+        m_todo_terms.push_back(std::make_pair(t, rational::one()));
+        lp::impq result(0);
+        while (!m_todo_terms.empty()) {
+            t = m_todo_terms.back().first;
+            rational coeff = m_todo_terms.back().second;
+            m_todo_terms.pop_back();
+            if (t.is_term()) {
+                const lp::lar_term& term = lp().get_term(t);
+                for (const auto& i : term) {
+                    m_todo_terms.push_back(std::make_pair(lp().column2tv(i.column()), coeff * i.coeff()));
+                }
+            }
+            else {
+                result += lp().get_column_value(t.id()) * coeff;
+            }
+        }
+        return result;
+    }
+
+    rational solver::get_value(theory_var v) const {
+        if (v == euf::null_theory_var || !lp().external_is_used(v)) {
+            return rational::zero();
+        }
+
+        auto t = get_tv(v);
+        if (m_variable_values.count(t.index()) > 0)
+            return m_variable_values[t.index()];
+
+        if (!t.is_term() && lp().is_fixed(t.id()))
+            return lp().column_lower_bound(t.id()).x;
+
+        if (!t.is_term())
+            return rational::zero();
+
+        m_todo_terms.push_back(std::make_pair(t, rational::one()));
+        rational result(0);
+        while (!m_todo_terms.empty()) {
+            auto t2 = m_todo_terms.back().first;
+            rational coeff = m_todo_terms.back().second;
+            m_todo_terms.pop_back();
+            if (t2.is_term()) {
+                const lp::lar_term& term = lp().get_term(t2);
+                for (const auto& i : term) {
+                    auto tv = lp().column2tv(i.column());
+                    if (m_variable_values.count(tv.index()) > 0) {
+                        result += m_variable_values[tv.index()] * coeff * i.coeff();
+                    }
+                    else {
+                        m_todo_terms.push_back(std::make_pair(tv, coeff * i.coeff()));
+                    }
+                }
+            }
+            else {
+                result += m_variable_values[t2.index()] * coeff;
+            }
+        }
+        m_variable_values[t.index()] = result;
+        return result;
+    }
+
+    void solver::random_update() {
+        if (m_nla)
+            return;
+        m_tmp_var_set.clear();
+        m_tmp_var_set.resize(get_num_vars());
+        m_model_eqs.reset();
+        svector<lpvar> vars;
+        theory_var sz = static_cast<theory_var>(get_num_vars());
+        for (theory_var v = 0; v < sz; ++v) {
+            ensure_column(v);
+            lp::column_index vj = lp().to_column_index(v);
+            SASSERT(!vj.is_null());
+            theory_var other = m_model_eqs.insert_if_not_there(v);
+            if (is_equal(v, other))
+                continue;
+            if (!lp().is_fixed(vj))
+                vars.push_back(vj.index());
+            else if (!m_tmp_var_set.contains(other)) {
+                lp::column_index other_j = lp().to_column_index(other);
+                if (!lp().is_fixed(other_j)) {
+                    m_tmp_var_set.insert(other);
+                    vars.push_back(other_j.index());
+                }
+            }
+        }
+        if (!vars.empty())
+            lp().random_update(vars.size(), vars.c_ptr());
+    }
+
+    bool solver::assume_eqs() {
+        TRACE("arith", display(tout););
+        random_update();
+        m_model_eqs.reset();
+        theory_var sz = static_cast<theory_var>(get_num_vars());
+        unsigned old_sz = m_assume_eq_candidates.size();
+        int start = s().rand()();
+        for (theory_var i = 0; i < sz; ++i) {
+            theory_var v = (i + start) % sz;
+            enode* n1 = var2enode(v);
+            ensure_column(v);
+            if (!can_get_ivalue(v))
+                continue;
+            theory_var other = m_model_eqs.insert_if_not_there(v);
+            TRACE("arith", tout << "insert: v" << v << " := " << get_value(v) << " found: v" << other << "\n";);
+            if (!is_equal(other, v))
+                m_assume_eq_candidates.push_back(std::make_pair(v, other));
+        }
+
+        if (m_assume_eq_candidates.size() > old_sz)
+            ctx.push(restore_size_trail<euf::solver, std::pair<theory_var, theory_var>, false>(m_assume_eq_candidates, old_sz));
+
+        return delayed_assume_eqs();
+    }
+
+    bool solver::delayed_assume_eqs() {
+        if (m_assume_eq_head == m_assume_eq_candidates.size())
+            return false;
+
+        ctx.push(value_trail<euf::solver, unsigned>(m_assume_eq_head));
+        while (m_assume_eq_head < m_assume_eq_candidates.size()) {
+            std::pair<theory_var, theory_var> const& p = m_assume_eq_candidates[m_assume_eq_head];
+            theory_var v1 = p.first;
+            theory_var v2 = p.second;
+            enode* n1 = var2enode(v1);
+            enode* n2 = var2enode(v2);
+            m_assume_eq_head++;
+            CTRACE("arith",
+                is_eq(v1, v2) && n1->get_root() != n2->get_root(),
+                tout << "assuming eq: v" << v1 << " = v" << v2 << "\n";);
+            if (!is_eq(v1, v2))
+                continue;
+            if (n1->get_root() == n2->get_root())
+                continue;
+            literal eq = eq_internalize(n1->get_expr(), n2->get_expr());
+            if (s().value(eq) != l_true)
+                return true;
+        }
+        return false;
+    }
+
+    bool solver::use_nra_model() {
+        if (m_nla && m_nla->use_nra_model()) {
+            if (!m_a1) {
+                m_a1 = alloc(scoped_anum, m_nla->am());
+                m_a2 = alloc(scoped_anum, m_nla->am());
+            }
+            return true;
+        }
+        return false;
+    }
+
+    bool solver::is_eq(theory_var v1, theory_var v2) {
+        if (use_nra_model()) {
+            return m_nla->am().eq(nl_value(v1, *m_a1), nl_value(v2, *m_a2));
+        }
+        else {
+            return get_ivalue(v1) == get_ivalue(v2);
+        }
+    }
+
+    sat::check_result solver::check() {
+        flet<bool> _is_learned(m_is_redundant, true);
+        reset_variable_values();
+        IF_VERBOSE(12, verbose_stream() << "final-check " << lp().get_status() << "\n");
+        SASSERT(lp().ax_is_correct());
+
+        if (lp().get_status() != lp::lp_status::OPTIMAL) {
+            switch (make_feasible()) {
+            case l_false:
+                get_infeasibility_explanation_and_set_conflict();
+                return sat::check_result::CR_CONTINUE;
+            case l_undef:
+                TRACE("arith", tout << "check feasiable is undef\n";);
+                return m.inc() ? sat::check_result::CR_CONTINUE : sat::check_result::CR_GIVEUP;
+            case l_true:
+                break;
+            default:
+                UNREACHABLE();
+            }
+        }
+
+        auto st = sat::check_result::CR_DONE;
+
+        TRACE("arith", ctx.display(tout););
+
+        switch (check_lia()) {
+        case l_true:
+            break;
+        case l_false:
+            return sat::check_result::CR_CONTINUE;
+        case l_undef:
+            TRACE("arith", tout << "check-lia giveup\n";);
+            st = sat::check_result::CR_CONTINUE;
+            break;
+        }
+
+        switch (check_nla()) {
+        case l_true:
+            break;
+        case l_false:
+            return sat::check_result::CR_CONTINUE;
+        case l_undef:
+            TRACE("arith", tout << "check-nra giveup\n";);
+            st = sat::check_result::CR_GIVEUP;
+            break;
+        }
+
+        if (delayed_assume_eqs()) {
+            ++m_stats.m_assume_eqs;
+            return sat::check_result::CR_CONTINUE;
+        }
+        if (assume_eqs()) {
+            ++m_stats.m_assume_eqs;
+            return sat::check_result::CR_CONTINUE;
+        }
+        if (m_not_handled != nullptr) {
+            TRACE("arith", tout << "unhandled operator " << mk_pp(m_not_handled, m) << "\n";);
+            st = sat::check_result::CR_GIVEUP;
+        }
+        return st;
+    }
+
+    nlsat::anum const& solver::nl_value(theory_var v, scoped_anum& r) const {
+        SASSERT(m_nla);
+        SASSERT(m_nla->use_nra_model());
+        auto t = get_tv(v);
+        if (t.is_term()) {
+
+            m_todo_terms.push_back(std::make_pair(t, rational::one()));
+            TRACE("nl_value", tout << "v" << v << " " << t.to_string() << "\n";);
+            TRACE("nl_value", tout << "v" << v << " := w" << t.to_string() << "\n";
+            lp().print_term(lp().get_term(t), tout) << "\n";);
+
+            m_nla->am().set(r, 0);
+            while (!m_todo_terms.empty()) {
+                rational wcoeff = m_todo_terms.back().second;
+                t = m_todo_terms.back().first;
+                m_todo_terms.pop_back();
+                lp::lar_term const& term = lp().get_term(t);
+                TRACE("nl_value", lp().print_term(term, tout) << "\n";);
+                scoped_anum r1(m_nla->am());
+                rational c1(0);
+                m_nla->am().set(r1, c1.to_mpq());
+                m_nla->am().add(r, r1, r);
+                for (auto const& arg : term) {
+                    auto wi = lp().column2tv(arg.column());
+                    c1 = arg.coeff() * wcoeff;
+                    if (wi.is_term()) {
+                        m_todo_terms.push_back(std::make_pair(wi, c1));
+                    }
+                    else {
+                        m_nla->am().set(r1, c1.to_mpq());
+                        m_nla->am().mul(m_nla->am_value(wi.id()), r1, r1);
+                        m_nla->am().add(r1, r, r);
+                    }
+                }
+            }
+            return r;
+        }
+        else {
+            return m_nla->am_value(t.id());
+        }
+    }
+
+    lbool solver::make_feasible() {
+        TRACE("pcs", tout << lp().constraints(););
+        auto status = lp().find_feasible_solution();
+        TRACE("arith_verbose", display(tout););
+        switch (status) {
+        case lp::lp_status::INFEASIBLE:
+            return l_false;
+        case lp::lp_status::FEASIBLE:
+        case lp::lp_status::OPTIMAL:
+            return l_true;
+        case lp::lp_status::TIME_EXHAUSTED:
+        default:
+            TRACE("arith", tout << "status treated as inconclusive: " << status << "\n";);
+            return l_undef;
+        }
+    }
+
+    lbool solver::check_lia() {
+        TRACE("arith", );
+        if (!m.inc())
+            return l_undef;
+        lbool lia_check = l_undef;
+        if (!check_idiv_bounds())
+            return l_false;
+
+        switch (m_lia->check(&m_explanation)) {
+        case lp::lia_move::sat:
+            lia_check = l_true;
+            break;
+
+        case lp::lia_move::branch: {
+            TRACE("arith", tout << "branch\n";);
+            app_ref b(m);
+            bool u = m_lia->is_upper();
+            auto const& k = m_lia->get_offset();
+            rational offset;
+            expr_ref t(m);
+            b = mk_bound(m_lia->get_term(), k, !u, offset, t);
+            IF_VERBOSE(4, verbose_stream() << "branch " << b << "\n";);
+            // branch on term >= k + 1
+            // branch on term <= k
+            // TBD: ctx().force_phase(ctx().get_literal(b));
+            // at this point we have a new unassigned atom that the 
+            // SAT core assigns a value to
+            lia_check = l_false;
+            ++m_stats.m_branch;
+            add_variable_bound(t, offset);
+            break;
+        }
+        case lp::lia_move::cut: {
+            TRACE("arith", tout << "cut\n";);
+            ++m_stats.m_gomory_cuts;
+            // m_explanation implies term <= k
+            reset_evidence();
+            for (auto const& ev : m_explanation)
+                set_evidence(ev.ci(), m_core, m_eqs);
+            // The call mk_bound() can set the m_infeasible_column in lar_solver
+            // so the explanation is safer to take before this call.
+            app_ref b = mk_bound(m_lia->get_term(), m_lia->get_offset(), !m_lia->is_upper());
+            IF_VERBOSE(4, verbose_stream() << "cut " << b << "\n");
+            literal lit = expr2literal(b);
+            assign(lit, m_core, m_eqs, m_params);
+            lia_check = l_false;
+            break;
+        }
+        case lp::lia_move::conflict:
+            TRACE("arith", tout << "conflict\n";);
+            // ex contains unsat core
+            set_conflict();
+            return l_false;
+        case lp::lia_move::undef:
+            TRACE("arith", tout << "lia undef\n";);
+            lia_check = l_undef;
+            break;
+        case lp::lia_move::continue_with_check:
+            lia_check = l_undef;
+            break;
+        default:
+            UNREACHABLE();
+        }
+        return lia_check;
+    }
+
+    void solver::assign(literal lit, literal_vector const& core, svector<enode_pair> const& eqs, vector<parameter> const& params) {
+        if (core.size() < small_lemma_size() && eqs.empty()) {
+            m_core2.reset();
+            for (auto const& c : core)
+                m_core2.push_back(~c);
+            m_core2.push_back(lit);
+            add_clause(m_core2);
+        }
+        else {
+            auto* jst = euf::th_propagation::mk(*this, core, eqs);
+            ctx.propagate(lit, jst->to_index());
+        }
+    }
+
+    void solver::get_infeasibility_explanation_and_set_conflict() {
+        m_explanation.clear();
+        lp().get_infeasibility_explanation(m_explanation);
+        set_conflict();
+    }
+
+    void solver::set_conflict() {
+        literal_vector core;
+        set_conflict_or_lemma(core, true);
+    }
+
+    void solver::set_conflict_or_lemma(literal_vector const& core, bool is_conflict) {
+        reset_evidence();
+        m_core.append(core);
+        ++m_num_conflicts;
+        ++m_stats.m_conflicts;
+        TRACE("arith", display(tout << "is-conflict: " << is_conflict << "\n"););
+        for (auto const& ev : m_explanation)
+            set_evidence(ev.ci(), m_core, m_eqs);
+        for (auto const& eq : m_eqs)
+            m_core.push_back(eq_internalize(eq.first->get_expr(), eq.second->get_expr()));
+        for (literal& c : m_core)
+            c.neg();
+        add_clause(m_core);
+    }
+
+    bool solver::is_infeasible() const {
+        return lp().get_status() == lp::lp_status::INFEASIBLE;
+    }
+
+    void solver::set_evidence(lp::constraint_index idx, literal_vector& core, svector<enode_pair>& eqs) {
+        if (idx == UINT_MAX) {
+            return;
+        }
+        switch (m_constraint_sources[idx]) {
+        case inequality_source: {
+            literal lit = m_inequalities[idx];
+            SASSERT(lit != sat::null_literal);
+            core.push_back(lit);
+            break;
+        }
+        case equality_source:
+            SASSERT(m_equalities[idx].first != nullptr);
+            SASSERT(m_equalities[idx].second != nullptr);
+            m_eqs.push_back(m_equalities[idx]);
+            break;
+        case definition_source:
+            // skip definitions (these are treated as hard constraints)
+            break;
+        default:
+            UNREACHABLE();
+            break;
+        }
+    }
+
+    void solver::add_variable_bound(expr* t, rational const& offset) {
+        if (!use_bounded_expansion())
+            return;
+        if (m_bounded_range_lit == sat::null_literal)
+            return;
+        // if term is not already bounded, add a range and assert max_bound => range
+        bound_info bi(offset, init_range());
+        if (m_term2bound_info.find(t, bi))
+            return;
+        expr_ref hi(a.mk_le(t, a.mk_int(offset + bi.m_range)), m);
+        expr_ref lo(a.mk_ge(t, a.mk_int(offset - bi.m_range)), m);
+        add_clause(~m_bounded_range_lit, mk_literal(hi));
+        add_clause(~m_bounded_range_lit, mk_literal(lo));
+        m_bound_terms.push_back(lo);
+        m_bound_terms.push_back(hi);
+        m_bound_terms.push_back(t);
+        m_predicate2term.insert(lo, t);
+        m_predicate2term.insert(hi, t);
+        m_term2bound_info.insert(t, bi);
+    }
+
+    void solver::reset_evidence() {
+        m_core.reset();
+        m_eqs.reset();
+        m_params.reset();
+    }
+
+    // create an eq atom representing "term = offset"
+    literal solver::mk_eq(lp::lar_term const& term, rational const& offset) {
+        u_map<rational> coeffs;
+        term2coeffs(term, coeffs);
+        bool isint = offset.is_int();
+        for (auto const& kv : coeffs) isint &= is_int(kv.m_key) && kv.m_value.is_int();
+        app_ref t = coeffs2app(coeffs, rational::zero(), isint);
+        app_ref s(a.mk_numeral(offset, isint), m);
+        return eq_internalize(t, s);
+    }
+    // create a bound atom representing term >= k is lower_bound is true, and term <= k if it is false
+    app_ref solver::mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound) {
+        rational offset;
+        expr_ref t(m);
+        return mk_bound(term, k, lower_bound, offset, t);
+    }
+
+    app_ref solver::mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound, rational& offset, expr_ref& t) {
+        offset = k;
+        u_map<rational> coeffs;
+        term2coeffs(term, coeffs);
+        bool is_int = true;
+        rational lc = denominator(k);
+        for (auto const& kv : coeffs) {
+            theory_var w = kv.m_key;
+            expr* o = var2expr(w);
+            is_int = a.is_int(o);
+            if (!is_int) break;
+            lc = lcm(lc, denominator(kv.m_value));
+        }
+
+        // ensure that coefficients are integers when all variables are integers as well.
+        if (is_int && !lc.is_one()) {
+            SASSERT(lc.is_pos());
+            offset *= lc;
+            for (auto& kv : coeffs) kv.m_value *= lc;
+        }
+
+        if (is_int) {
+            // 3x + 6y >= 5 -> x + 3y >= 5/3, then x + 3y >= 2
+            // 3x + 6y <= 5 -> x + 3y <= 1
+            rational g = gcd_reduce(coeffs);
+            if (!g.is_one()) {
+                if (lower_bound) {
+                    TRACE("arith", tout << "lower: " << offset << " / " << g << " = " << offset / g << " >= " << ceil(offset / g) << "\n";);
+                    offset = ceil(offset / g);
+                }
+                else {
+                    TRACE("arith", tout << "upper: " << offset << " / " << g << " = " << offset / g << " <= " << floor(offset / g) << "\n";);
+                    offset = floor(offset / g);
+                }
+            }
+        }
+        if (!coeffs.empty() && coeffs.begin()->m_value.is_neg()) {
+            offset.neg();
+            lower_bound = !lower_bound;
+            for (auto& kv : coeffs) kv.m_value.neg();
+        }
+
+        // CTRACE("arith", is_int,
+        //        lp().print_term(term, tout << "term: ") << "\n";
+        //        tout << "offset: " << offset << " gcd: " << g << "\n";);
+
+        app_ref atom(m);
+        t = coeffs2app(coeffs, rational::zero(), is_int);
+        if (lower_bound) {
+            atom = a.mk_ge(t, a.mk_numeral(offset, is_int));
+        }
+        else {
+            atom = a.mk_le(t, a.mk_numeral(offset, is_int));
+        }
+
+        TRACE("arith", tout << t << ": " << atom << "\n";
+        lp().print_term(term, tout << "bound atom: ") << (lower_bound ? " >= " : " <= ") << k << "\n";);
+        b_internalize(atom);
+        return atom;
+    }
+
+    void solver::term2coeffs(lp::lar_term const& term, u_map<rational>& coeffs) {
+        term2coeffs(term, coeffs, rational::one());
+    }
+
+    void solver::term2coeffs(lp::lar_term const& term, u_map<rational>& coeffs, rational const& coeff) {
+        TRACE("arith", lp().print_term(term, tout) << "\n";);
+        for (const auto& ti : term) {
+            theory_var w;
+            auto tv = lp().column2tv(ti.column());
+            if (tv.is_term()) {
+                lp::lar_term const& term1 = lp().get_term(tv);
+                rational coeff2 = coeff * ti.coeff();
+                term2coeffs(term1, coeffs, coeff2);
+                continue;
+            }
+            else {
+                w = lp().local_to_external(tv.id());
+                SASSERT(w >= 0);
+                TRACE("arith", tout << (tv.id()) << ": " << w << "\n";);
+            }
+            rational c0(0);
+            coeffs.find(w, c0);
+            coeffs.insert(w, c0 + ti.coeff() * coeff);
+        }
+    }
+
+    app_ref solver::coeffs2app(u_map<rational> const& coeffs, rational const& offset, bool is_int) {
+        expr_ref_vector args(m);
+        for (auto const& kv : coeffs) {
+            theory_var w = kv.m_key;
+            expr* o = var2expr(w);
+            if (kv.m_value.is_zero())
+                continue;
+            else if (kv.m_value.is_one())
+                args.push_back(o);
+            else
+                args.push_back(a.mk_mul(a.mk_numeral(kv.m_value, is_int), o));
+        }
+
+        if (!offset.is_zero() || args.empty()) 
+            args.push_back(a.mk_numeral(offset, is_int));
+
+        if (args.size() == 1)
+            return app_ref(to_app(args[0].get()), m);
+
+        return app_ref(a.mk_add(args.size(), args.c_ptr()), m);        
+    }
+
+    app_ref solver::mk_term(lp::lar_term const& term, bool is_int) {
+        u_map<rational> coeffs;
+        term2coeffs(term, coeffs);
+        return coeffs2app(coeffs, rational::zero(), is_int);
+    }
+
+    rational solver::gcd_reduce(u_map<rational>& coeffs) {
+        rational g(0);
+        for (auto const& kv : coeffs) 
+            g = gcd(g, kv.m_value);        
+        if (g.is_zero())
+            return rational::one();
+        if (!g.is_one()) 
+            for (auto& kv : coeffs) 
+                kv.m_value /= g;                    
+        return g;
+    }
+
+    void solver::false_case_of_check_nla(const nla::lemma& l) {
+        m_lemma = l; //todo avoid the copy
+        m_explanation = l.expl();
+        literal_vector core;
+        for (auto const& ineq : m_lemma.ineqs()) {
+            bool is_lower = true, pos = true, is_eq = false;
+            switch (ineq.cmp()) {
+            case lp::LE: is_lower = false; pos = false;  break;
+            case lp::LT: is_lower = true;  pos = true; break;
+            case lp::GE: is_lower = true;  pos = false;  break;
+            case lp::GT: is_lower = false; pos = true; break;
+            case lp::EQ: is_eq = true; pos = false; break;
+            case lp::NE: is_eq = true; pos = true; break;
+            default: UNREACHABLE();
+            }
+            TRACE("arith", tout << "is_lower: " << is_lower << " pos " << pos << "\n";);
+            app_ref atom(m);
+            // TBD utility: lp::lar_term term = mk_term(ineq.m_poly);
+            // then term is used instead of ineq.m_term
+            if (is_eq) {
+                core.push_back(~mk_eq(ineq.term(), ineq.rs()));
+            }
+            else {
+                // create term >= 0 (or term <= 0)
+                core.push_back(~ctx.expr2literal(mk_bound(ineq.term(), ineq.rs(), is_lower)));
+            }
+        }
+        set_conflict_or_lemma(core, false);
+    }
+
+    lbool solver::check_nla() {
+        if (!m.inc()) {
+            TRACE("arith", tout << "canceled\n";);
+            return l_undef;
+        }
+        if (!m_nla) {
+            TRACE("arith", tout << "no nla\n";);
+            return l_true;
+        }
+        if (!m_nla->need_check())
+            return l_true;
+
+        m_a1 = nullptr; m_a2 = nullptr;
+        lbool r = m_nla->check(m_nla_lemma_vector);
+        switch (r) {
+        case l_false: {
+            for (const nla::lemma& l : m_nla_lemma_vector)
+                false_case_of_check_nla(l);
+            break;
+        }
+        case l_true:
+            if (assume_eqs())
+                return l_false;
+            break;
+        case l_undef:
+            break;
+        }
+        return r;
+    }
+
+    void solver::get_antecedents(literal l, sat::ext_justification_idx idx, literal_vector& r, bool probing) {
+        auto& jst = euf::th_propagation::from_index(idx);
+        for (auto lit : euf::th_propagation::lits(jst))
+            r.push_back(lit);
+        for (auto eq : euf::th_propagation::eqs(jst))
+            ctx.add_antecedent(eq.first, eq.second);
+
+        if (!probing && ctx.use_drat()) {
+            literal_vector lits;
+            for (auto lit : euf::th_propagation::lits(jst))
+                lits.push_back(~lit);
+            lits.push_back(l);
+            ctx.get_drat().add(lits, status());
+            for (auto eq : euf::th_propagation::eqs(jst))
+                IF_VERBOSE(0, verbose_stream() << "drat-log with equalities is TBD " << eq.first->get_expr_id() << "\n");
+        }
+    }
 }
diff --git a/src/sat/smt/arith_solver.h b/src/sat/smt/arith_solver.h
index 4e24119a9..cb0380766 100644
--- a/src/sat/smt/arith_solver.h
+++ b/src/sat/smt/arith_solver.h
@@ -19,6 +19,16 @@ Author:
 #include "ast/ast_trail.h"
 #include "sat/smt/sat_th.h"
 #include "ast/arith_decl_plugin.h"
+#include "math/lp/lp_solver.h"
+#include "math/lp/lp_primal_simplex.h"
+#include "math/lp/lp_dual_simplex.h"
+#include "math/lp/indexed_value.h"
+#include "math/lp/lar_solver.h"
+#include "math/lp/nla_solver.h"
+#include "math/lp/lp_types.h"
+#include "math/lp/lp_api.h"
+#include "math/polynomial/algebraic_numbers.h"
+#include "math/polynomial/polynomial.h"
 
 namespace euf {
     class solver;
@@ -26,43 +36,382 @@ namespace euf {
 
 namespace arith {
 
+    typedef ptr_vector<lp_api::bound> lp_bounds;
+    typedef lp::var_index lpvar;
+    typedef euf::theory_var theory_var;
+    typedef euf::theory_id theory_id;
+    typedef euf::enode enode;
+    typedef euf::enode_pair enode_pair;
+    typedef sat::literal literal;
+    typedef sat::bool_var bool_var;
+    typedef sat::literal_vector literal_vector;
+
     class solver : public euf::th_euf_solver {
-        typedef euf::theory_var theory_var;
-        typedef euf::theory_id theory_id;
-        typedef sat::literal literal;
-        typedef sat::bool_var bool_var;
-        typedef sat::literal_vector literal_vector;
-        typedef union_find<solver, euf::solver>  array_union_find;
 
-
-        struct stats {
-            void reset() { memset(this, 0, sizeof(*this)); }
-            stats() { reset(); }
+        struct scope {
+            unsigned m_bounds_lim;
+            unsigned m_idiv_lim;
+            unsigned m_asserted_qhead;
+            unsigned m_asserted_lim;
+            unsigned m_underspecified_lim;
+            expr* m_not_handled;
         };
 
+        class resource_limit : public lp::lp_resource_limit {
+            solver& m_imp;
+        public:
+            resource_limit(solver& i) : m_imp(i) { }
+            bool get_cancel_flag() override { return !m_imp.m.inc(); }
+        };
+
+        struct var_value_eq {
+            solver& m_th;
+            var_value_eq(solver& th) :m_th(th) {}
+            bool operator()(theory_var v1, theory_var v2) const {
+                if (m_th.is_int(v1) != m_th.is_int(v2)) {
+                    return false;
+                }
+                return m_th.is_eq(v1, v2);
+            }
+        };
+
+        struct var_value_hash {
+            solver& m_th;
+            var_value_hash(solver& th) :m_th(th) {}
+            unsigned operator()(theory_var v) const {
+                if (m_th.use_nra_model()) {
+                    return m_th.is_int(v);
+                }
+                else {
+                    return (unsigned)std::hash<lp::impq>()(m_th.get_ivalue(v));
+                }
+            }
+        };
+        int_hashtable<var_value_hash, var_value_eq>   m_model_eqs;
+
+
+
+
+        // temporary values kept during internalization
+        struct internalize_state {
+            expr_ref_vector     m_terms;
+            vector<rational>    m_coeffs;
+            svector<theory_var> m_vars;
+            rational            m_offset;
+            ptr_vector<expr>    m_to_ensure_enode, m_to_ensure_var;
+            internalize_state(ast_manager& m) : m_terms(m) {}
+            void reset() {
+                m_terms.reset();
+                m_coeffs.reset();
+                m_offset.reset();
+                m_vars.reset();
+                m_to_ensure_enode.reset();
+                m_to_ensure_var.reset();
+            }
+        };
+        ptr_vector<internalize_state> m_internalize_states;
+        unsigned                      m_internalize_head{ 0 };
+
+        class scoped_internalize_state {
+            solver& m_imp;
+            internalize_state& m_st;
+
+            internalize_state& push_internalize(solver& i) {
+                if (i.m_internalize_head == i.m_internalize_states.size()) {
+                    i.m_internalize_states.push_back(alloc(internalize_state, i.m));
+                }
+                internalize_state& st = *i.m_internalize_states[i.m_internalize_head++];
+                st.reset();
+                return st;
+            }
+        public:
+            scoped_internalize_state(solver& i) : m_imp(i), m_st(push_internalize(i)) {}
+            ~scoped_internalize_state() { --m_imp.m_internalize_head; }
+            expr_ref_vector& terms() { return m_st.m_terms; }
+            vector<rational>& coeffs() { return m_st.m_coeffs; }
+            svector<theory_var>& vars() { return m_st.m_vars; }
+            rational& offset() { return m_st.m_offset; }
+            ptr_vector<expr>& to_ensure_enode() { return m_st.m_to_ensure_enode; }
+            ptr_vector<expr>& to_ensure_var() { return m_st.m_to_ensure_var; }
+            void push(expr* e, rational c) { m_st.m_terms.push_back(e); m_st.m_coeffs.push_back(c); }
+            void set_back(unsigned i) {
+                if (terms().size() == i + 1) return;
+                terms()[i] = terms().back();
+                coeffs()[i] = coeffs().back();
+                terms().pop_back();
+                coeffs().pop_back();
+            }
+        };
+
+        typedef vector<std::pair<rational, lpvar>> var_coeffs;
+        vector<rational>         m_columns;
+        var_coeffs               m_left_side;              // constraint left side
+
+        mutable std::unordered_map<lpvar, rational> m_variable_values; // current model
+        lpvar m_one_var   { UINT_MAX };
+        lpvar m_zero_var  { UINT_MAX };
+        lpvar m_rone_var  { UINT_MAX };
+        lpvar m_rzero_var { UINT_MAX };
+
+        enum constraint_source {
+            inequality_source,
+            equality_source,
+            definition_source,
+            null_source
+        };
+        svector<constraint_source>                    m_constraint_sources;
+        svector<literal>                              m_inequalities;    // asserted rows corresponding to inequality literals.
+        svector<euf::enode_pair>                      m_equalities;      // asserted rows corresponding to equalities.
+        svector<theory_var>                           m_definitions;     // asserted rows corresponding to definitions
+
+        literal_vector  m_asserted;
+        expr* m_not_handled{ nullptr };
+        ptr_vector<app>        m_underspecified;
+        ptr_vector<expr>       m_idiv_terms;
+        vector<ptr_vector<lp_api::bound> > m_use_list;        // bounds where variables are used.
+
+        // attributes for incremental version:
+        u_map<lp_api::bound*>      m_bool_var2bound;
+        vector<lp_bounds>      m_bounds;
+        unsigned_vector        m_unassigned_bounds;
+        unsigned_vector        m_bounds_trail;
+        unsigned               m_asserted_qhead{ 0 };
+
+        svector<std::pair<theory_var, theory_var> >       m_assume_eq_candidates;
+        unsigned                                          m_assume_eq_head{ 0 };
+        lp::u_set                                         m_tmp_var_set;
+
+        unsigned                                          m_num_conflicts{ 0 };
+        lp_api::stats                                     m_stats;
+        svector<scope>                                    m_scopes;
+
+        // non-linear arithmetic
+        scoped_ptr<nla::solver>  m_nla;
+        scoped_ptr<scoped_anum>  m_a1, m_a2;
+
+        // integer arithmetic
+        scoped_ptr<lp::int_solver>   m_lia;
+
+
+        scoped_ptr<lp::lar_solver>   m_solver;
+        resource_limit               m_resource_limit;
+        lp_bounds                    m_new_bounds;
+        symbol                       m_farkas;
+        lp::lp_bound_propagator<solver> m_bp;
+        mutable vector<std::pair<lp::tv, rational>> m_todo_terms;
+
+        // lemmas
+        lp::explanation     m_explanation;
+        vector<nla::lemma>  m_nla_lemma_vector;
+        literal_vector      m_core, m_core2;
+        svector<enode_pair> m_eqs;
+        vector<parameter>   m_params;
+        nla::lemma          m_lemma;
+
 
         arith_util        a;
-        stats             m_stats;
+
+        bool is_int(theory_var v) const { return is_int(var2enode(v)); }
+        bool is_int(euf::enode* n) const { return a.is_int(n->get_expr()); }
+        bool is_real(theory_var v) const { return is_real(var2enode(v)); }
+        bool is_real(euf::enode* n) const { return a.is_real(n->get_expr()); }
+        
 
         // internalize
-        bool visit(expr* e) override;
-        bool visited(expr* e) override;
-        bool post_visit(expr* e, bool sign, bool root) override;
-        void ensure_var(euf::enode* n);
+        lpvar add_const(int c, lpvar& var, bool is_int);
+        lpvar get_one(bool is_int);
+        lpvar get_zero(bool is_int);
+        void ensure_nla();
+        void found_unsupported(expr* n);
+        void found_underspecified(expr* n);
+
+        void linearize_term(expr* term, scoped_internalize_state& st);
+        void linearize_ineq(expr* lhs, expr* rhs, scoped_internalize_state& st);
+        void linearize(scoped_internalize_state& st);
+
+        void add_eq_constraint(lp::constraint_index index, enode* n1, enode* n2);
+        void add_ineq_constraint(lp::constraint_index index, literal lit);
+        void add_def_constraint(lp::constraint_index index);
+        void add_def_constraint(lp::constraint_index index, theory_var v);
+        void add_def_constraint_and_equality(lpvar vi, lp::lconstraint_kind kind, const rational& bound);
+        void internalize_args(app* t, bool force = false);
+        theory_var internalize_power(app* t, app* n, unsigned p);
+        theory_var internalize_mul(app* t);
+        theory_var internalize_def(expr* term);
+        theory_var internalize_def(expr* term, scoped_internalize_state& st);
+        theory_var internalize_linearized_def(expr* term, scoped_internalize_state& st);
+        void init_left_side(scoped_internalize_state& st);
+        bool internalize_term(expr* term);
+        bool internalize_atom(expr* atom);
+        bool is_unit_var(scoped_internalize_state& st);
+        bool is_one(scoped_internalize_state& st);
+        bool is_zero(scoped_internalize_state& st);
+        enode* mk_enode(expr* e);
+
+        lpvar register_theory_var_in_lar_solver(theory_var v);
+        theory_var mk_evar(expr* e);
+        bool has_var(expr* e);
+        bool reflect(expr* n) const;
+
+        lpvar get_lpvar(theory_var v) const;
+        lp::tv get_tv(theory_var v) const;
 
         // axioms
         void mk_div_axiom(expr* p, expr* q);
         void mk_to_int_axiom(app* n);
-        void mk_is_int_axiom(app* n);
+        void mk_is_int_axiom(expr* n);
+        void mk_idiv_mod_axioms(expr* p, expr* q);
+        void mk_rem_axiom(expr* dividend, expr* divisor);
+        void mk_bound_axioms(lp_api::bound& b);
+        void mk_bound_axiom(lp_api::bound& b1, lp_api::bound& b2);
+        void flush_bound_axioms();
+
+        // bounds
+        struct compare_bounds {
+            bool operator()(lp_api::bound* a1, lp_api::bound* a2) const { return a1->get_value() < a2->get_value(); }
+        };
+
+        typedef lp_bounds::iterator iterator;
+
+        lp_bounds::iterator first(
+            lp_api::bound_kind kind,
+            iterator it,
+            iterator end);
+
+        lp_bounds::iterator next_inf(
+            lp_api::bound* a1,
+            lp_api::bound_kind kind,
+            iterator it,
+            iterator end,
+            bool& found_compatible);
+
+        lp_bounds::iterator next_sup(
+            lp_api::bound* a1,
+            lp_api::bound_kind kind,
+            iterator it,
+            iterator end,
+            bool& found_compatible);
+
+        void propagate_eqs(lp::tv t, lp::constraint_index ci, lp::lconstraint_kind k, lp_api::bound& b, rational const& value);
+        void propagate_basic_bounds(unsigned qhead);
+        void propagate_bounds_with_lp_solver();
+        void propagate_bound(literal lit, lp_api::bound& b);
+        void propagate_lp_solver_bound(const lp::implied_bound& be);
+        void refine_bound(theory_var v, const lp::implied_bound& be);
+        literal is_bound_implied(lp::lconstraint_kind k, rational const& value, lp_api::bound const& b) const;
+        void assert_bound(bool is_true, lp_api::bound& b);
+        void mk_eq_axiom(theory_var v1, theory_var v2);
+        void assert_idiv_mod_axioms(theory_var u, theory_var v, theory_var w, rational const& r);
+        lp_api::bound* mk_var_bound(bool_var bv, theory_var v, lp_api::bound_kind bk, rational const& bound);
+        lp::lconstraint_kind bound2constraint_kind(bool is_int, lp_api::bound_kind bk, bool is_true);
+        void fixed_var_eh(theory_var v1, rational const& bound) {}
+        bool set_upper_bound(lp::tv t, lp::constraint_index ci, rational const& v) { return set_bound(t, ci, v, false); }
+        bool set_lower_bound(lp::tv t, lp::constraint_index ci, rational const& v) { return set_bound(t, ci, v, true); }
+        bool set_bound(lp::tv tv, lp::constraint_index ci, rational const& v, bool is_lower);
+
+        typedef std::pair<lp::constraint_index, rational> constraint_bound;
+        vector<constraint_bound>        m_lower_terms;
+        vector<constraint_bound>        m_upper_terms;
+        vector<constraint_bound>        m_history;
+
+
+        // solving
+        void report_equality_of_fixed_vars(unsigned vi1, unsigned vi2);
+        void reserve_bounds(theory_var v);
+
+        void updt_unassigned_bounds(theory_var v, int inc);
 
         void pop_core(unsigned n) override;
-        
+        void push_core() override;
+        void del_bounds(unsigned old_size);
+
+        bool all_zeros(vector<rational> const& v) const;
+
+        bound_prop_mode propagation_mode() const;
+        void init_variable_values();
+        void reset_variable_values();
+
+        bool can_get_value(theory_var v) const;
+        bool can_get_bound(theory_var v) const;
+        bool can_get_ivalue(theory_var v) const;
+        void ensure_column(theory_var v);
+        lp::impq get_ivalue(theory_var v) const;
+        rational get_value(theory_var v) const;
+
+        void random_update();
+        bool assume_eqs();
+        bool delayed_assume_eqs();
+        bool is_eq(theory_var v1, theory_var v2);
+        bool use_nra_model();
+
+        lbool make_feasible();
+        lbool check_lia();
+        lbool check_nla();
+        void add_variable_bound(expr* t, rational const& offset);
+        bool is_infeasible() const;
+
+        nlsat::anum const& nl_value(theory_var v, scoped_anum& r) const;
+
+
+        bool has_bound(lpvar vi, lp::constraint_index& ci, rational const& bound, bool is_lower);
+        bool has_lower_bound(lpvar vi, lp::constraint_index& ci, rational const& bound);
+        bool has_upper_bound(lpvar vi, lp::constraint_index& ci, rational const& bound);
+
+        /*
+         * Facility to put a small box around integer variables used in branch and bounds.
+         */
+
+        struct bound_info {
+            rational m_offset;
+            unsigned m_range;
+            bound_info() {}
+            bound_info(rational const& o, unsigned r) :m_offset(o), m_range(r) {}
+        };
+        unsigned                  m_bounded_range_idx;  // current size of bounded range.
+        literal                   m_bounded_range_lit;  // current bounded range literal
+        expr_ref_vector           m_bound_terms; // predicates used for bounds
+        expr_ref                  m_bound_predicate;
+
+        obj_map<expr, expr*>      m_predicate2term;
+        obj_map<expr, bound_info> m_term2bound_info;
+
+        bool use_bounded_expansion() const { return get_config().m_arith_bounded_expansion; }
+        unsigned small_lemma_size() const { return get_config().m_arith_small_lemma_size; }
+        bool propagate_eqs() const { return get_config().m_arith_propagate_eqs && m_num_conflicts < get_config().m_arith_propagation_threshold; }
+        bool should_propagate() const { return bound_prop_mode::BP_NONE != propagation_mode(); }
+        bool should_refine_bounds() const { return bound_prop_mode::BP_REFINE == propagation_mode() && s().at_search_lvl(); }
+
+        unsigned init_range() const { return 5; }
+        unsigned max_range() const { return 20; }
+
+        void reset_evidence();
+        bool check_idiv_bounds();
+        bool is_bounded(expr* n);
+
+        app_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound);
+        app_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound, rational& offset, expr_ref& t);
+        literal mk_eq(lp::lar_term const& term, rational const& offset);
+
+        rational gcd_reduce(u_map<rational>& coeffs);
+        app_ref mk_term(lp::lar_term const& term, bool is_int);
+        app_ref coeffs2app(u_map<rational> const& coeffs, rational const& offset, bool is_int);
+        void term2coeffs(lp::lar_term const& term, u_map<rational>& coeffs, rational const& coeff);
+        void term2coeffs(lp::lar_term const& term, u_map<rational>& coeffs);
+
+        void get_infeasibility_explanation_and_set_conflict();
+        void set_conflict();
+        void set_conflict_or_lemma(literal_vector const& core, bool is_conflict);
+        void set_evidence(lp::constraint_index idx, literal_vector& core, svector<enode_pair>& eqs);
+        void assign(literal lit, literal_vector const& core, svector<enode_pair> const& eqs, vector<parameter> const& params);
+
+        void false_case_of_check_nla(const nla::lemma& l);        
+
     public:
         solver(euf::solver& ctx, theory_id id);
         ~solver() override;
         bool is_external(bool_var v) override { return false; }
         bool propagate(literal l, sat::ext_constraint_idx idx) override { UNREACHABLE(); return false; }
-        void get_antecedents(literal l, sat::ext_justification_idx idx, literal_vector& r, bool probing) override {}
+        void get_antecedents(literal l, sat::ext_justification_idx idx, literal_vector& r, bool probing) override;
         void asserted(literal l) override;
         sat::check_result check() override;
 
@@ -71,14 +420,23 @@ namespace arith {
         std::ostream& display_constraint(std::ostream& out, sat::ext_constraint_idx idx) const override;
         void collect_statistics(statistics& st) const override;
         euf::th_solver* clone(sat::solver* s, euf::solver& ctx) override;
-        void new_eq_eh(euf::th_eq const& eq) override;
+        bool use_diseqs() const override { return true; }
+        void new_eq_eh(euf::th_eq const& eq) override { mk_eq_axiom(eq.v1(), eq.v2()); }
+        void new_diseq_eh(euf::th_eq const& de) override { mk_eq_axiom(de.v1(), de.v2()); }
         bool unit_propagate() override;
         void add_value(euf::enode* n, model& mdl, expr_ref_vector& values) override;
-        void add_dep(euf::enode* n, top_sort<euf::enode>& dep) override;
         sat::literal internalize(expr* e, bool sign, bool root, bool learned) override;
         void internalize(expr* e, bool redundant) override;
-        euf::theory_var mk_var(euf::enode* n) override;
         void apply_sort_cnstr(euf::enode* n, sort* s) override {}
         bool is_shared(theory_var v) const override;
+        lbool get_phase(bool_var v) override;
+
+        // bounds and equality propagation callbacks
+        lp::lar_solver& lp() { return *m_solver; }
+        lp::lar_solver const& lp() const { return *m_solver; }
+        bool is_equal(theory_var x, theory_var y) const;
+        void add_eq(lpvar u, lpvar v, lp::explanation const& e);
+        void consume(rational const& v, lp::constraint_index j);
+        bool bound_is_interesting(unsigned vi, lp::lconstraint_kind kind, const rational& bval) const;
     };
 }
diff --git a/src/sat/smt/array_axioms.cpp b/src/sat/smt/array_axioms.cpp
index 7b7419efd..157a3b79a 100644
--- a/src/sat/smt/array_axioms.cpp
+++ b/src/sat/smt/array_axioms.cpp
@@ -49,7 +49,7 @@ namespace array {
         case axiom_record::kind_t::is_default:
             return assert_default(r);
         case axiom_record::kind_t::is_extensionality:
-            return assert_extensionality(r.n->get_arg(0)->get_expr(), r.n->get_arg(1)->get_expr());
+            return assert_extensionality(r.n->get_expr(), r.select->get_expr());
         case axiom_record::kind_t::is_congruence:
             return assert_congruent_axiom(r.n->get_expr(), r.select->get_expr());
         default:
@@ -214,12 +214,11 @@ namespace array {
     bool solver::assert_extensionality(expr* e1, expr* e2) {
         TRACE("array", tout << "extensionality-axiom: " << mk_bounded_pp(e1, m) << " == " << mk_bounded_pp(e2, m) << "\n";);
         ++m_stats.m_num_extensionality_axiom;
-        func_decl_ref_vector* funcs = nullptr;
-        VERIFY(m_sort2diff.find(m.get_sort(e1), funcs));
+        func_decl_ref_vector const& funcs = sort2diff(m.get_sort(e1));
         expr_ref_vector args1(m), args2(m);
         args1.push_back(e1);
         args2.push_back(e2);
-        for (func_decl* f : *funcs) {
+        for (func_decl* f : funcs) {
             expr* k = m.mk_app(f, e1, e2);
             args1.push_back(k);
             args2.push_back(k);
@@ -527,6 +526,7 @@ namespace array {
         unsigned num_vars = get_num_vars();
         for (unsigned i = 0; i < num_vars; i++) {
             euf::enode * n = var2enode(i);
+            
             if (!a.is_array(n->get_expr())) {
                 continue;
             }
diff --git a/src/sat/smt/array_internalize.cpp b/src/sat/smt/array_internalize.cpp
index f612be3c5..44e506d93 100644
--- a/src/sat/smt/array_internalize.cpp
+++ b/src/sat/smt/array_internalize.cpp
@@ -22,8 +22,10 @@ namespace array {
 
     sat::literal solver::internalize(expr* e, bool sign, bool root, bool redundant) { 
         SASSERT(m.is_bool(e));
-        if (!visit_rec(m, e, sign, root, redundant))
+        if (!visit_rec(m, e, sign, root, redundant)) {
+            TRACE("array", tout << mk_pp(e, m) << "\n";);
             return sat::null_literal;
+        }
         return expr2literal(e);
     }
 
@@ -81,7 +83,7 @@ namespace array {
     }
 
     void solver::internalize_ext(euf::enode* n) {
-        push_axiom(extensionality_axiom(n));
+        push_axiom(extensionality_axiom(n->get_arg(0), n->get_arg(1)));
     }
 
     void solver::internalize_default(euf::enode* n) {
@@ -95,6 +97,8 @@ namespace array {
     }
 
     bool solver::visit(expr* e) {
+        if (visited(e))
+            return true;
         if (!is_app(e) || to_app(e)->get_family_id() != get_id()) {
             ctx.internalize(e, m_is_redundant);
             euf::enode* n = expr2enode(e);
@@ -192,5 +196,20 @@ namespace array {
         return false;
     }
 
+    func_decl_ref_vector const& solver::sort2diff(sort* s) {
+        func_decl_ref_vector* result = nullptr;
+        if (m_sort2diff.find(s, result))
+            return *result;
+        
+        unsigned dimension = get_array_arity(s);
+        result = alloc(func_decl_ref_vector, m);
+        for (unsigned i = 0; i < dimension; ++i) 
+            result->push_back(a.mk_array_ext(s, i));
+        m_sort2diff.insert(s, result);
+        ctx.push(insert_map<euf::solver, obj_map<sort, func_decl_ref_vector*>, sort*>(m_sort2diff, s));
+        ctx.push(new_obj_trail<euf::solver,func_decl_ref_vector>(result));
+        return *result;
+    }
+
 }
 
diff --git a/src/sat/smt/array_model.cpp b/src/sat/smt/array_model.cpp
index 25ebb13ef..1ea6eb4de 100644
--- a/src/sat/smt/array_model.cpp
+++ b/src/sat/smt/array_model.cpp
@@ -92,11 +92,6 @@ namespace array {
                 if (!value || value == fi->get_else())
                     continue;
                 args.reset();
-                bool relevant = true;
-                for (unsigned i = 1; relevant && i < p->num_args(); ++i)
-                    relevant = ctx.is_relevant(p->get_arg(i)->get_root());
-                if (!relevant)
-                    continue;
                 for (unsigned i = 1; i < p->num_args(); ++i) 
                     args.push_back(values.get(p->get_arg(i)->get_root_id()));    
                 fi->insert_entry(args.c_ptr(), value);
diff --git a/src/sat/smt/array_solver.cpp b/src/sat/smt/array_solver.cpp
index 95cb52ff2..4df6df2a5 100644
--- a/src/sat/smt/array_solver.cpp
+++ b/src/sat/smt/array_solver.cpp
@@ -88,6 +88,8 @@ namespace array {
         m_constraint->initialize(m_constraint.get(), this);
     }
 
+    solver::~solver() {}
+
     sat::check_result solver::check() {
         force_push();
         // flet<bool> _is_redundant(m_is_redundant, true);
@@ -108,6 +110,8 @@ namespace array {
     }
 
     std::ostream& solver::display(std::ostream& out) const {
+        if (get_num_vars() > 0)
+            out << "array\n";
         for (unsigned i = 0; i < get_num_vars(); ++i) {
             auto& d = get_var_data(i);
             out << var2enode(i)->get_expr_id() << " " << mk_bounded_pp(var2expr(i), m, 2) << "\n";
@@ -117,6 +121,7 @@ namespace array {
         }
         return out;
     }
+
     std::ostream& solver::display_info(std::ostream& out, char const* id, euf::enode_vector const& v) const {
         if (v.empty())
             return out;
@@ -163,6 +168,11 @@ namespace array {
         m_find.merge(eq.v1(), eq.v2());
     }
 
+    void solver::new_diseq_eh(euf::th_eq const& eq) {
+        force_push();
+        push_axiom(extensionality_axiom(var2enode(eq.v1()), var2enode(eq.v2())));
+    }
+
     bool solver::unit_propagate() {
         if (m_qhead == m_axiom_trail.size())
             return false;
diff --git a/src/sat/smt/array_solver.h b/src/sat/smt/array_solver.h
index 0d7747ee0..b46dfd761 100644
--- a/src/sat/smt/array_solver.h
+++ b/src/sat/smt/array_solver.h
@@ -67,6 +67,7 @@ namespace array {
         array_union_find                     m_find;
 
         theory_var find(theory_var v) { return m_find.find(v); }
+        func_decl_ref_vector const& sort2diff(sort* s);
 
         // internalize
         bool visit(expr* e) override;
@@ -131,7 +132,7 @@ namespace array {
         axiom_record select_axiom(euf::enode* s, euf::enode* n) { return axiom_record(axiom_record::kind_t::is_select, n, s); }
         axiom_record default_axiom(euf::enode* n) { return axiom_record(axiom_record::kind_t::is_default, n); }
         axiom_record store_axiom(euf::enode* n) { return axiom_record(axiom_record::kind_t::is_store, n); }
-        axiom_record extensionality_axiom(euf::enode* n) { return axiom_record(axiom_record::kind_t::is_extensionality, n); }
+        axiom_record extensionality_axiom(euf::enode* x, euf::enode* y) { return axiom_record(axiom_record::kind_t::is_extensionality, x, y); }
         axiom_record congruence_axiom(euf::enode* a, euf::enode* b) { return axiom_record(axiom_record::kind_t::is_congruence, a, b); }
 
         scoped_ptr<sat::constraint_base> m_constraint;
@@ -189,7 +190,7 @@ namespace array {
         std::ostream& display_info(std::ostream& out, char const* id, euf::enode_vector const& v) const;
     public:
         solver(euf::solver& ctx, theory_id id);
-        ~solver() override {}
+        ~solver() override;
         bool is_external(bool_var v) override { return false; }
         bool propagate(literal l, sat::ext_constraint_idx idx) override { UNREACHABLE(); return false; }
         void get_antecedents(literal l, sat::ext_justification_idx idx, literal_vector& r, bool probing) override {}
@@ -202,6 +203,8 @@ namespace array {
         void collect_statistics(statistics& st) const override;
         euf::th_solver* clone(sat::solver* s, euf::solver& ctx) override;
         void new_eq_eh(euf::th_eq const& eq) override;
+        bool use_diseqs() const override { return true; }
+        void new_diseq_eh(euf::th_eq const& eq) override;
         bool unit_propagate() override;
         void add_value(euf::enode* n, model& mdl, expr_ref_vector& values) override;
         void add_dep(euf::enode* n, top_sort<euf::enode>& dep) override;
diff --git a/src/sat/smt/bv_internalize.cpp b/src/sat/smt/bv_internalize.cpp
index 0766510a5..0690e6f06 100644
--- a/src/sat/smt/bv_internalize.cpp
+++ b/src/sat/smt/bv_internalize.cpp
@@ -257,12 +257,6 @@ namespace bv {
     }
 
     void solver::get_bits(theory_var v, expr_ref_vector& r) {
-        for (literal lit : m_bits[v]) {
-            if (!literal2expr(lit))
-                ctx.display(std::cout << "Missing expression: " << lit << "\n");
-            SASSERT(literal2expr(lit));
-        }
-
         for (literal lit : m_bits[v]) 
             r.push_back(literal2expr(lit));        
     }
diff --git a/src/sat/smt/bv_solver.cpp b/src/sat/smt/bv_solver.cpp
index b1c9fcf2b..8165b2743 100644
--- a/src/sat/smt/bv_solver.cpp
+++ b/src/sat/smt/bv_solver.cpp
@@ -56,7 +56,6 @@ namespace bv {
         m_ackerman(*this),
         m_bb(m, get_config()),
         m_find(*this) {
-        ctx.get_egraph().set_th_propagates_diseqs(id);
     }
 
     void solver::fixed_var_eh(theory_var v1) {
diff --git a/src/sat/smt/bv_solver.h b/src/sat/smt/bv_solver.h
index d3deea85d..b9894a34f 100644
--- a/src/sat/smt/bv_solver.h
+++ b/src/sat/smt/bv_solver.h
@@ -221,7 +221,6 @@ namespace bv {
         void get_arg_bits(app* n, unsigned idx, expr_ref_vector& r);
         void fixed_var_eh(theory_var v);
         bool is_bv(theory_var v) const { return bv.is_bv(var2expr(v)); }
-        sat::status status() const { return sat::status::th(m_is_redundant, get_id());  }
         void register_true_false_bit(theory_var v, unsigned i);
         void add_bit(theory_var v, sat::literal lit);
         atom* mk_atom(sat::bool_var b);
diff --git a/src/sat/smt/euf_internalize.cpp b/src/sat/smt/euf_internalize.cpp
index 263910669..9f61d6f48 100644
--- a/src/sat/smt/euf_internalize.cpp
+++ b/src/sat/smt/euf_internalize.cpp
@@ -21,6 +21,9 @@ Author:
 namespace euf {
 
     void solver::internalize(expr* e, bool redundant) {
+        SASSERT(!get_enode(e) || get_enode(e)->bool_var() < UINT_MAX);
+        if (get_enode(e))
+            return;
         if (si.is_bool_op(e))
             attach_lit(si.internalize(e, redundant), e);
         else if (auto* ext = expr2solver(e))
@@ -31,23 +34,28 @@ namespace euf {
     }
 
     sat::literal solver::internalize(expr* e, bool sign, bool root, bool redundant) {
-        euf::enode* n = m_egraph.find(e);
+        euf::enode* n = get_enode(e);
         if (n) {
             if (m.is_bool(e)) {
                 VERIFY(!s().was_eliminated(n->bool_var()));
+                SASSERT(n->bool_var() != UINT_MAX);
                 return literal(n->bool_var(), sign);
             }
+            TRACE("euf", tout << "non-bool\n";);
             return sat::null_literal;
         }
         if (si.is_bool_op(e))
             return attach_lit(si.internalize(e, redundant), e);
         if (auto* ext = expr2solver(e))
             return ext->internalize(e, sign, root, redundant);
-        if (!visit_rec(m, e, sign, root, redundant))
+        if (!visit_rec(m, e, sign, root, redundant)) {
+            TRACE("euf", tout << "visit-rec\n";);          
             return sat::null_literal;
-        SASSERT(m_egraph.find(e));
+        }
+        SASSERT(get_enode(e));
         if (m.is_bool(e))
             return literal(si.to_bool_var(e), sign);
+        std::cout << "internalize-non-bool\n";
         return sat::null_literal;
     }
 
@@ -138,8 +146,9 @@ namespace euf {
         enode* n = m_egraph.find(e);
         if (!n) 
             n = m_egraph.mk(e, 0, nullptr); 
+        SASSERT(n->bool_var() == UINT_MAX || n->bool_var() == v);
         m_egraph.set_bool_var(n, v);
-        if (!m.is_true(e) && !m.is_false(e))
+        if (m.is_eq(e) || m.is_or(e) || m.is_and(e) || m.is_not(e))
             m_egraph.set_merge_enabled(n, false);
         return lit;
     }
@@ -174,7 +183,9 @@ namespace euf {
                 }
             }
             s().mk_clause(lits, st);
-        }
+            if (relevancy_enabled())
+                add_root(lits.size(), lits.c_ptr());
+    }
         else {
             // g(f(x_i)) = x_i
             // f(x_1) = a + .... + f(x_n) = a >= 2
@@ -198,6 +209,8 @@ namespace euf {
             expr_ref at_least2(pb.mk_at_least_k(eqs.size(), eqs.c_ptr(), 2), m);
             sat::literal lit = si.internalize(at_least2, m_is_redundant);
             s().mk_clause(1, &lit, st);
+            if (relevancy_enabled())
+                add_root(1, &lit);
         }
     }
 
@@ -216,6 +229,8 @@ namespace euf {
                     expr_ref eq = mk_eq(args[i]->get_expr(), args[j]->get_expr());
                     sat::literal lit = internalize(eq, true, false, m_is_redundant);
                     s().add_clause(1, &lit, st);
+                    if (relevancy_enabled())
+                        add_root(1, &lit);
                 }
             }
         }
@@ -233,6 +248,8 @@ namespace euf {
                 expr_ref eq = mk_eq(fapp, fresh);
                 sat::literal lit = internalize(eq, false, false, m_is_redundant);
                 s().add_clause(1, &lit, st);
+                if (relevancy_enabled())
+                    add_root(1, &lit);
             }
         }
     }
diff --git a/src/sat/smt/euf_model.cpp b/src/sat/smt/euf_model.cpp
index 411087d1e..0be72a4df 100644
--- a/src/sat/smt/euf_model.cpp
+++ b/src/sat/smt/euf_model.cpp
@@ -139,7 +139,7 @@ namespace euf {
                 mbS->add_value(n, *mdl, m_values);
             else if (auto* mbE = expr2solver(e))
                 mbE->add_value(n, *mdl, m_values);
-            else  {
+            else {
                 IF_VERBOSE(1, verbose_stream() << "no model values created for " << mk_pp(e, m) << "\n");
             }                
         }
diff --git a/src/sat/smt/euf_relevancy.cpp b/src/sat/smt/euf_relevancy.cpp
index 5f0d9010c..8896954c0 100644
--- a/src/sat/smt/euf_relevancy.cpp
+++ b/src/sat/smt/euf_relevancy.cpp
@@ -22,6 +22,14 @@ Author:
 
 namespace euf {
 
+    bool solver::is_relevant(expr* e) const { 
+        return m_relevant_expr_ids.get(e->get_id(), true); 
+    }
+
+    bool solver::is_relevant(enode* n) const { 
+        return m_relevant_expr_ids.get(n->get_expr_id(), true); 
+    }
+
     void solver::ensure_dual_solver() {
         if (!m_dual_solver)
             m_dual_solver = alloc(sat::dual_solver, s().rlimit());
diff --git a/src/sat/smt/euf_solver.cpp b/src/sat/smt/euf_solver.cpp
index 0be711400..2b1da2240 100644
--- a/src/sat/smt/euf_solver.cpp
+++ b/src/sat/smt/euf_solver.cpp
@@ -23,6 +23,7 @@ Author:
 #include "sat/smt/bv_solver.h"
 #include "sat/smt/euf_solver.h"
 #include "sat/smt/array_solver.h"
+#include "sat/smt/arith_solver.h"
 #include "sat/smt/q_solver.h"
 #include "sat/smt/fpa_solver.h"
 
@@ -99,6 +100,7 @@ namespace euf {
         bv_util bvu(m);
         array_util au(m);
         fpa_util fpa(m);
+        arith_util arith(m);
         if (pb.get_family_id() == fid) 
             ext = alloc(sat::ba_solver, *this, fid);
         else if (bvu.get_family_id() == fid) 
@@ -107,13 +109,17 @@ namespace euf {
             ext = alloc(array::solver, *this, fid);
         else if (fpa.get_family_id() == fid) 
             ext = alloc(fpa::solver, *this);
-
+        else if (arith.get_family_id() == fid)
+            ext = alloc(arith::solver, *this, fid);
+        
         if (ext) {
             if (use_drat())
                 s().get_drat().add_theory(fid, ext->name());
             ext->set_solver(m_solver);
             ext->push_scopes(s().num_scopes());
             add_solver(fid, ext);
+            if (ext->use_diseqs())
+                m_egraph.set_th_propagates_diseqs(fid);
         }
         else if (f) 
             unhandled_function(f);
@@ -260,7 +266,7 @@ namespace euf {
             euf::enode* nb = sign ? mk_false() : mk_true();
             m_egraph.merge(n, nb, c);
         }
-        else if (sign && n->is_equality())
+        else if (sign && n->is_equality()) 
             m_egraph.new_diseq(n);        
     }
 
@@ -419,7 +425,7 @@ namespace euf {
         m_var_trail.shrink(s.m_var_lim);        
         m_scopes.shrink(m_scopes.size() - n);
         SASSERT(m_egraph.num_scopes() == m_scopes.size());
-        TRACE("euf", tout << "pop to: " << m_scopes.size() << "\n";);
+        TRACE("euf", display(tout << "pop to: " << m_scopes.size() << "\n"););
     }
 
     void solver::start_reinit(unsigned n) {
@@ -472,7 +478,7 @@ namespace euf {
             VERIFY(lit.var() == kv.m_value);
             attach_lit(lit, kv.m_key);            
         }
-        TRACE("euf", tout << "replay done\n";);
+        TRACE("euf", display(tout << "replay done\n"););
     }
 
     void solver::pre_simplify() {
@@ -493,7 +499,6 @@ namespace euf {
     }
 
     lbool solver::get_phase(bool_var v) { 
-        TRACE("euf", tout << "phase: " << v << "\n";);            
         auto* ext = bool_var2solver(v);
         if (ext)
             return ext->get_phase(v);
diff --git a/src/sat/smt/euf_solver.h b/src/sat/smt/euf_solver.h
index 44b3dae9c..8cee826ca 100644
--- a/src/sat/smt/euf_solver.h
+++ b/src/sat/smt/euf_solver.h
@@ -192,6 +192,7 @@ namespace euf {
             if (m_conflict) dealloc(sat::constraint_base::mem2base_ptr(m_conflict));
             if (m_eq) dealloc(sat::constraint_base::mem2base_ptr(m_eq));
             if (m_lit) dealloc(sat::constraint_base::mem2base_ptr(m_lit));
+            m_trail.reset();
         }
 
         struct scoped_set_translate {
@@ -211,9 +212,9 @@ namespace euf {
         
         sat::sat_internalizer& get_si() { return si; }
         ast_manager& get_manager() { return m; }
-        enode* get_enode(expr* e) { return m_egraph.find(e); }
-        sat::literal expr2literal(expr* e) const { return literal(si.to_bool_var(e), false); }
-        sat::literal enode2literal(enode* e) const { return expr2literal(e->get_expr()); }
+        enode* get_enode(expr* e) const { return m_egraph.find(e); }
+        sat::literal expr2literal(expr* e) const { return enode2literal(get_enode(e)); }
+        sat::literal enode2literal(enode* e) const { return sat::literal(e->bool_var(), false); }
         smt_params const& get_config() const { return m_config; }
         region& get_region() { return m_trail.get_region(); }
         egraph& get_egraph() { return m_egraph; }
@@ -299,8 +300,8 @@ namespace euf {
         void add_root(unsigned n, sat::literal const* lits);
         void add_aux(unsigned n, sat::literal const* lits);
         void track_relevancy(sat::bool_var v);
-        bool is_relevant(expr* e) const { return m_relevant_expr_ids.get(e->get_id(), true); }
-        bool is_relevant(enode* n) const { return m_relevant_expr_ids.get(n->get_expr_id(), true); }
+        bool is_relevant(expr* e) const;
+        bool is_relevant(enode* n) const;
 
         // model construction
         void update_model(model_ref& mdl);
diff --git a/src/sat/smt/sat_th.cpp b/src/sat/smt/sat_th.cpp
index 31fa24078..137c7fb31 100644
--- a/src/sat/smt/sat_th.cpp
+++ b/src/sat/smt/sat_th.cpp
@@ -204,5 +204,27 @@ namespace euf {
         return b_internalize(ctx.mk_eq(a, b)); 
     }
 
+    unsigned th_propagation::get_obj_size(unsigned num_lits, unsigned num_eqs) {
+        return sizeof(th_propagation) + sizeof(sat::literal) * num_lits + sizeof(enode_pair) * num_eqs;
+    }
 
+    th_propagation::th_propagation(sat::literal_vector const& lits, enode_pair_vector const& eqs) {
+        m_num_literals = lits.size();
+        m_num_eqs = eqs.size();
+        m_literals = reinterpret_cast<literal*>(reinterpret_cast<char*>(this) + sizeof(th_propagation));
+        unsigned i = 0;
+        for (sat::literal lit : lits)
+            m_literals[i++] = lit;
+        m_eqs = reinterpret_cast<enode_pair*>(reinterpret_cast<char*>(this) + sizeof(th_propagation) + sizeof(literal) * m_num_literals);
+        i = 0;
+        for (auto eq : eqs)
+            m_eqs[i++] = eq;
+    }
+
+    th_propagation* th_propagation::mk(th_euf_solver& th, sat::literal_vector const& lits, enode_pair_vector const& eqs) {
+        region& r = th.ctx.get_region();
+        void* mem = r.allocate(get_obj_size(lits.size(), eqs.size()));
+        sat::constraint_base::initialize(mem, &th);
+        return new (sat::constraint_base::ptr2mem(mem)) th_propagation(lits, eqs);
+    }
 }
diff --git a/src/sat/smt/sat_th.h b/src/sat/smt/sat_th.h
index 2a1f04912..e6bb2be0a 100644
--- a/src/sat/smt/sat_th.h
+++ b/src/sat/smt/sat_th.h
@@ -25,12 +25,12 @@ Author:
 namespace euf {
 
     class solver;
-    
+
     class th_internalizer {
     protected:
         euf::enode_vector     m_args;
         svector<sat::eframe>  m_stack;
-        bool                  m_is_redundant { false };
+        bool                  m_is_redundant{ false };
 
         bool visit_rec(ast_manager& m, expr* e, bool sign, bool root, bool redundant);
 
@@ -47,10 +47,12 @@ namespace euf {
 
         sat::literal b_internalize(expr* e) { return internalize(e, false, false, m_is_redundant); }
 
+        sat::literal mk_literal(expr* e) { return b_internalize(e); }
+
         /**
            \brief Apply (interpreted) sort constraints on the given enode.
         */
-        virtual void apply_sort_cnstr(enode * n, sort * s) {}
+        virtual void apply_sort_cnstr(enode* n, sort* s) {}
 
         /**
            \record that an equality has been internalized.
@@ -72,7 +74,7 @@ namespace euf {
         virtual ~th_model_builder() {}
 
         /**
-           \brief compute the value for enode \c n and store the value in \c values 
+           \brief compute the value for enode \c n and store the value in \c values
            for the root of the class of \c n.
          */
         virtual void add_value(euf::enode* n, model& mdl, expr_ref_vector& values) {}
@@ -80,7 +82,7 @@ namespace euf {
         /**
            \brief compute dependencies for node n
          */
-        virtual void add_dep(euf::enode* n, top_sort<euf::enode>& dep) { dep.insert(n, nullptr);  }
+        virtual void add_dep(euf::enode* n, top_sort<euf::enode>& dep) { dep.insert(n, nullptr); }
 
         /**
            \brief should function be included in model.
@@ -95,11 +97,11 @@ namespace euf {
 
     class th_solver : public sat::extension, public th_model_builder, public th_decompile, public th_internalizer {
     protected:
-        ast_manager &       m;
+        ast_manager& m;
     public:
-        th_solver(ast_manager& m, symbol const& name, euf::theory_id id): extension(name, id), m(m) {}
+        th_solver(ast_manager& m, symbol const& name, euf::theory_id id) : extension(name, id), m(m) {}
 
-        virtual th_solver* clone(sat::solver* s, euf::solver& ctx) = 0;  
+        virtual th_solver* clone(sat::solver* s, euf::solver& ctx) = 0;
 
         virtual void new_eq_eh(euf::th_eq const& eq) {}
 
@@ -112,11 +114,13 @@ namespace euf {
         */
         virtual bool is_shared(theory_var v) const { return false; }
 
+        sat::status status() const { return sat::status::th(m_is_redundant, get_id()); }
+
     };
 
     class th_euf_solver : public th_solver {
     protected:
-        solver &            ctx;
+        solver& ctx;
         euf::enode_vector   m_var2enode;
         unsigned_vector     m_var2enode_lim;
         unsigned            m_num_scopes{ 0 };
@@ -124,7 +128,7 @@ namespace euf {
         smt_params const& get_config() const;
         sat::literal expr2literal(expr* e) const;
         region& get_region();
-        
+
 
         sat::status mk_status();
         bool add_unit(sat::literal lit);
@@ -136,7 +140,7 @@ namespace euf {
         bool add_clause(sat::literal_vector const& lits);
         void add_equiv(sat::literal a, sat::literal b);
         void add_equiv_and(sat::literal a, sat::literal_vector const& bs);
-       
+
 
         bool is_true(sat::literal lit);
         bool is_true(sat::literal a, sat::literal b) { return is_true(a) || is_true(b); }
@@ -154,17 +158,19 @@ namespace euf {
 
         virtual void push_core();
         virtual void pop_core(unsigned n);
-        void force_push() { 
+        void force_push() {
             CTRACE("euf", m_num_scopes > 0, tout << "push-core " << m_num_scopes << "\n";);
-            for (; m_num_scopes > 0; --m_num_scopes) push_core(); 
+            for (; m_num_scopes > 0; --m_num_scopes) push_core();
         }
 
+        friend class th_propagation;
+
     public:
         th_euf_solver(euf::solver& ctx, symbol const& name, euf::theory_id id);
         virtual ~th_euf_solver() {}
-        virtual theory_var mk_var(enode * n);
-        unsigned get_num_vars() const { return m_var2enode.size();}
-        enode* expr2enode(expr* e) const; 
+        virtual theory_var mk_var(enode* n);
+        unsigned get_num_vars() const { return m_var2enode.size(); }
+        enode* expr2enode(expr* e) const;
         enode* var2enode(theory_var v) const { return m_var2enode[v]; }
         expr* var2expr(theory_var v) const { return var2enode(v)->get_expr(); }
         expr* bool_var2expr(sat::bool_var v) const;
@@ -172,11 +178,46 @@ namespace euf {
         expr_ref literal2expr(sat::literal lit) const { expr* e = bool_var2expr(lit.var()); return lit.sign() ? expr_ref(m.mk_not(e), m) : expr_ref(e, m); }
         theory_var get_th_var(enode* n) const { return n->get_th_var(get_id()); }
         theory_var get_th_var(expr* e) const;
-        trail_stack<euf::solver> & get_trail_stack();
+        trail_stack<euf::solver>& get_trail_stack();
         bool is_attached_to_var(enode* n) const;
         bool is_root(theory_var v) const { return var2enode(v)->is_root(); }
         void push() override { m_num_scopes++; }
-        void pop(unsigned n) override;       
+        void pop(unsigned n) override;
+    };
+
+
+    class th_propagation {
+        unsigned    m_num_literals;
+        unsigned    m_num_eqs;
+        sat::literal*    m_literals;
+        enode_pair* m_eqs;
+        static unsigned get_obj_size(unsigned num_lits, unsigned num_eqs);
+        th_propagation(sat::literal_vector const& lits, enode_pair_vector const& eqs);
+    public:
+        static th_propagation* mk(th_euf_solver& th, sat::literal_vector const& lits, enode_pair_vector const& eqs);
+
+        sat::ext_constraint_idx to_index() const {
+            return sat::constraint_base::mem2base(this);
+        }
+        static th_propagation& from_index(size_t idx) {
+            return *reinterpret_cast<th_propagation*>(sat::constraint_base::from_index(idx)->mem());
+        }
+
+        class lits {
+            th_propagation& th;
+        public:
+            lits(th_propagation& th) : th(th) {}
+            sat::literal const* begin() const { return th.m_literals; }
+            sat::literal const* end() const { return th.m_literals + th.m_num_literals; }
+        };
+
+        class eqs {
+            th_propagation& th;
+        public:
+            eqs(th_propagation& th) : th(th) {}
+            enode_pair const* begin() const { return th.m_eqs; }
+            enode_pair const* end() const { return th.m_eqs + th.m_num_eqs; }
+        };
 
     };
 
diff --git a/src/sat/tactic/goal2sat.cpp b/src/sat/tactic/goal2sat.cpp
index 1cad3c0ae..7cabedf25 100644
--- a/src/sat/tactic/goal2sat.cpp
+++ b/src/sat/tactic/goal2sat.cpp
@@ -285,7 +285,7 @@ struct goal2sat::imp : public sat::sat_internalizer {
             else {                
                 if (!is_uninterp_const(t)) {
                     if (m_euf) {
-                        convert_euf(t, root, sign);                        
+                        convert_euf(t, root, sign);  
                         return;
                     }
                     else
@@ -632,7 +632,7 @@ struct goal2sat::imp : public sat::sat_internalizer {
             flet<bool> _top(m_top_level, false);
             lit = euf->internalize(e, sign, root, m_is_redundant);           
         }
-        if (lit == sat::null_literal)
+        if (lit == sat::null_literal) 
             return;
         if (top_level_relevant())
             euf->track_relevancy(lit.var());
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 4e158534d..5f35dd5d3 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -26,6 +26,7 @@
 #include "math/lp/lar_solver.h"
 #include "math/lp/nla_solver.h"
 #include "math/lp/lp_types.h"
+#include "math/lp/lp_api.h"
 #include "math/polynomial/algebraic_numbers.h"
 #include "math/polynomial/polynomial.h"
 #include "util/nat_set.h"
@@ -50,91 +51,6 @@
 
 typedef lp::var_index lpvar;
 
-namespace lp_api {
-enum bound_kind { lower_t, upper_t };
-
-std::ostream& operator<<(std::ostream& out, bound_kind const& k) {
-    switch (k) {
-    case lower_t: return out << "<=";
-    case upper_t: return out << ">=";
-    }
-    return out;
-}
-
-class bound { 
-    smt::bool_var    m_bv;
-    smt::theory_var  m_var;
-    lpvar            m_vi;
-    bool             m_is_int;
-    rational         m_value;
-    bound_kind       m_bound_kind;
-    lp::constraint_index m_constraints[2];
-
-public:
-    bound(smt::bool_var bv, smt::theory_var v, lpvar vi, bool is_int, rational const & val, bound_kind k, lp::constraint_index ct, lp::constraint_index cf):
-        m_bv(bv),
-        m_var(v),
-        m_vi(vi),
-        m_is_int(is_int),
-        m_value(val),
-        m_bound_kind(k) {
-        m_constraints[0] = cf;
-        m_constraints[1] = ct;
-    }
-    virtual ~bound() {}
-    smt::theory_var get_var() const { return m_var; }
-    lp::tv tv() const { return lp::tv::raw(m_vi); }
-    smt::bool_var get_bv() const { return m_bv; }
-    bound_kind get_bound_kind() const { return m_bound_kind; }
-    bool is_int() const { return m_is_int; }
-    rational const& get_value() const { return m_value; }
-    lp::constraint_index get_constraint(bool b) const { return m_constraints[b]; }
-    inf_rational get_value(bool is_true) const { 
-        if (is_true) return inf_rational(m_value);                         // v >= value or v <= value
-        if (m_is_int) {
-            SASSERT(m_value.is_int());
-            if (m_bound_kind == lower_t) return inf_rational(m_value - rational::one()); // v <= value - 1
-            return inf_rational(m_value + rational::one());                              // v >= value + 1
-        }
-        else {
-            if (m_bound_kind == lower_t) return inf_rational(m_value, false);  // v <= value - epsilon
-            return inf_rational(m_value, true);                                // v >= value + epsilon
-        }
-    } 
-    virtual std::ostream& display(std::ostream& out) const {
-        return out << m_value << "  " << get_bound_kind() << " v" << get_var();
-    }
-};
-
-std::ostream& operator<<(std::ostream& out, bound const& b) {
-    return b.display(out);
-}
-
-struct stats {
-    unsigned m_assert_lower;
-    unsigned m_assert_upper;
-    unsigned m_bounds_propagations;
-    unsigned m_num_iterations;
-    unsigned m_num_iterations_with_no_progress;
-    unsigned m_need_to_solve_inf;
-    unsigned m_fixed_eqs;
-    unsigned m_conflicts;
-    unsigned m_bound_propagations1;
-    unsigned m_bound_propagations2;
-    unsigned m_assert_diseq;
-    unsigned m_gomory_cuts;
-    unsigned m_assume_eqs;
-    unsigned m_branch;
-    stats() { reset(); }
-    void reset() {
-        memset(this, 0, sizeof(*this));
-    }
-};
-
-typedef optional<inf_rational> opt_inf_rational;
-
-
-}
 
 namespace smt {
 
@@ -387,7 +303,7 @@ class theory_lra::imp {
     }
 
     void found_underspecified(expr* n) {
-        if (is_app(n) && is_underspecified(to_app(n))) {
+        if (a.is_underspecified(n)) {
             TRACE("arith", tout << "Unhandled: " << mk_pp(n, m) << "\n";);
             m_underspecified.push_back(to_app(n));
         }
@@ -415,26 +331,6 @@ class theory_lra::imp {
 
     }
 
-    bool is_numeral(expr* term, rational& r) {
-        rational mul(1);
-        do {
-            if (a.is_numeral(term, r)) {
-                r *= mul;
-                return true;
-            }
-            if (a.is_uminus(term, term)) {
-                mul.neg();
-                continue;
-            }
-            if (a.is_to_real(term, term)) {
-                continue;
-            }                
-            return false;
-        }
-        while (false);
-        return false;
-    }
-
     void linearize_term(expr* term, scoped_internalize_state& st) {
         st.push(term, rational::one());
         linearize(st);
@@ -471,12 +367,12 @@ class theory_lra::imp {
                     st.push(to_app(n)->get_arg(i), -coeffs[index]);
                 }
             }
-            else if (a.is_mul(n, n1, n2) && is_numeral(n1, r)) {
+            else if (a.is_mul(n, n1, n2) && a.is_extended_numeral(n1, r)) {
                 coeffs[index] *= r;
                 terms[index] = n2;
                 st.to_ensure_enode().push_back(n1);
             }
-            else if (a.is_mul(n, n1, n2) && is_numeral(n2, r)) {
+            else if (a.is_mul(n, n1, n2) && a.is_extended_numeral(n2, r)) {
                 coeffs[index] *= r;
                 terms[index] = n1;
                 st.to_ensure_enode().push_back(n2);
@@ -487,7 +383,7 @@ class theory_lra::imp {
                 vars.push_back(v);
                 ++index;
             }
-            else if(a.is_power(n, n1, n2) && is_app(n1) && is_numeral(n2, r) && r.is_unsigned() && r <= 10) {
+            else if(a.is_power(n, n1, n2) && is_app(n1) && a.is_extended_numeral(n2, r) && r.is_unsigned() && r <= 10) {
                 theory_var v = internalize_power(to_app(n), to_app(n1), r.get_unsigned());
                 coeffs[vars.size()] = coeffs[index];
                 vars.push_back(v);
@@ -672,27 +568,9 @@ class theory_lra::imp {
         ctx().mk_th_axiom(get_id(), l1, l2, l3, num_params, params);
     }
 
-    bool is_underspecified(app* n) const {
-        if (n->get_family_id() == get_id()) {
-            switch (n->get_decl_kind()) {
-            case OP_DIV:
-            case OP_IDIV:
-            case OP_REM:
-            case OP_MOD:
-            case OP_DIV0:
-            case OP_IDIV0:
-            case OP_REM0:
-            case OP_MOD0:
-                return true;
-            default:
-                break;
-            }
-        }
-        return false;
-    }
 
     bool reflect(app* n) const {
-        return params().m_arith_reflect || is_underspecified(n);          
+        return params().m_arith_reflect || a.is_underspecified(n);          
     }
 
     bool has_var(expr* n) {
@@ -703,7 +581,7 @@ class theory_lra::imp {
         return th.is_attached_to_var(e);
     }
 
-    void ensure_bounds(theory_var v) {
+    void reserve_bounds(theory_var v) {
         while (m_bounds.size() <= static_cast<unsigned>(v)) {
             m_bounds.push_back(lp_bounds());
             m_unassigned_bounds.push_back(0);
@@ -720,7 +598,7 @@ class theory_lra::imp {
             v = th.mk_var(e);
             TRACE("arith", tout << "fresh var: v" << v << " " << mk_pp(n, m) << "\n";);
             SASSERT(m_bounds.size() <= static_cast<unsigned>(v) || m_bounds[v].empty());
-            ensure_bounds(v);
+            reserve_bounds(v);
             ctx().attach_th_var(e, &th, v);
         }
         else {
@@ -1032,19 +910,19 @@ public:
         bool_var bv = ctx().mk_bool_var(atom);
         m_bool_var2bound.erase(bv);
         ctx().set_var_theory(bv, get_id());
-        if (a.is_le(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) {
+        if (a.is_le(atom, n1, n2) && a.is_extended_numeral(n2, r) && is_app(n1)) {
             v = internalize_def(to_app(n1));
             k = lp_api::upper_t;
         }
-        else if (a.is_ge(atom, n1, n2) && is_numeral(n2, r) && is_app(n1)) {
+        else if (a.is_ge(atom, n1, n2) && a.is_extended_numeral(n2, r) && is_app(n1)) {
             v = internalize_def(to_app(n1));
             k = lp_api::lower_t;
         }  
-        else if (a.is_le(atom, n1, n2) && is_numeral(n1, r) && is_app(n2)) {
+        else if (a.is_le(atom, n1, n2) && a.is_extended_numeral(n1, r) && is_app(n2)) {
             v = internalize_def(to_app(n2));
             k = lp_api::lower_t;
         }
-        else if (a.is_ge(atom, n1, n2) && is_numeral(n1, r) && is_app(n2)) {
+        else if (a.is_ge(atom, n1, n2) && a.is_extended_numeral(n1, r) && is_app(n2)) {
             v = internalize_def(to_app(n2));
             k = lp_api::upper_t;
         }
@@ -1891,24 +1769,6 @@ public:
      * 
      */
 
-    bool is_bounded(expr* n) {
-        expr* x = nullptr, *y = nullptr;
-        while (true) {
-            if (a.is_idiv(n, x, y) && a.is_numeral(y)) {
-                n = x;
-            }
-            else if (a.is_mod(n, x, y) && a.is_numeral(y)) {
-                return true;
-            }
-            else if (a.is_numeral(n)) {
-                return true;
-            }
-            else {
-                return false;
-            }
-        }
-    }
-
     bool check_idiv_bounds() {
         if (m_idiv_terms.empty()) {
             return true;
@@ -1918,7 +1778,7 @@ public:
             expr* n = m_idiv_terms[i];
             expr* p = nullptr, *q = nullptr;
             VERIFY(a.is_idiv(n, p, q));
-            theory_var v  = mk_var(n);
+            theory_var v = mk_var(n);
             theory_var v1 = mk_var(p);
 
             if (!can_get_ivalue(v1))
@@ -1937,7 +1797,7 @@ public:
             }
 
             if (a.is_numeral(q, r2) && r2.is_pos()) {
-                if (!is_bounded(n)) {
+                if (!a.is_bounded(n)) {
                     TRACE("arith", tout << "unbounded " << expr_ref(n, m) << "\n";);
                     continue;
                 }
@@ -1957,7 +1817,7 @@ public:
                 // used to normalize inequalities so they 
                 // don't appear as 8*x >= 15, but x >= 2
                 expr *n1 = nullptr, *n2 = nullptr;
-                if (a.is_mul(p, n1, n2) && is_numeral(n1, mul) && mul.is_pos()) {
+                if (a.is_mul(p, n1, n2) && a.is_extended_numeral(n1, mul) && mul.is_pos()) {
                     p = n2;
                     hi = floor(hi/mul);
                     lo = ceil(lo/mul);
@@ -2271,7 +2131,7 @@ public:
         }
         else {
             for (enode * parent : r->get_const_parents()) {
-                if (is_underspecified(parent->get_owner())) {
+                if (a.is_underspecified(parent->get_owner())) {
                     return true;
                 }
             }
@@ -2391,7 +2251,7 @@ public:
 
         TRACE("arith", tout << "v" << v << " " << be.kind() << " " << be.m_bound << "\n";);
 
-        ensure_bounds(v);
+        reserve_bounds(v);
 
             
         if (m_unassigned_bounds[v] == 0 && !should_refine_bounds()) {
@@ -3132,26 +2992,6 @@ public:
     bool set_lower_bound(lp::tv t, lp::constraint_index ci, rational const& v) { return set_bound(t, ci, v, true);   }
 
     vector<constraint_bound> m_history;
-    template<typename Ctx, typename T, bool CallDestructors=true>
-    class history_trail : public trail<Ctx> {
-        vector<T, CallDestructors> & m_dst;
-        unsigned                     m_idx;
-        vector<T, CallDestructors> & m_hist;
-    public:
-        history_trail(vector<T, CallDestructors> & v, unsigned idx, vector<T, CallDestructors> & hist):
-            m_dst(v),
-            m_idx(idx),
-            m_hist(hist) {}
-
-        ~history_trail() override {
-        }
-
-        void undo(Ctx & ctx) override {
-            m_dst[m_idx] = m_hist.back();
-            m_hist.pop_back();
-        }
-    };
-
 
     bool set_bound(lp::tv tv, lp::constraint_index ci, rational const& v, bool is_lower) {
         if (tv.is_term()) {
@@ -3923,38 +3763,9 @@ public:
 
     void collect_statistics(::statistics & st) const {
         m_arith_eq_adapter.collect_statistics(st);
-        st.update("arith-lower", m_stats.m_assert_lower);
-        st.update("arith-upper", m_stats.m_assert_upper);
-        st.update("arith-propagations", m_stats.m_bounds_propagations);
-        st.update("arith-iterations", m_stats.m_num_iterations);
-        st.update("arith-factorizations", lp().settings().stats().m_num_factorizations);
-        st.update("arith-pivots", m_stats.m_need_to_solve_inf);
-        st.update("arith-plateau-iterations", m_stats.m_num_iterations_with_no_progress);
-        st.update("arith-fixed-eqs", m_stats.m_fixed_eqs);
-        st.update("arith-conflicts", m_stats.m_conflicts);
-        st.update("arith-bound-propagations-lp", m_stats.m_bound_propagations1);
-        st.update("arith-bound-propagations-cheap", m_stats.m_bound_propagations2);
-        st.update("arith-diseq", m_stats.m_assert_diseq);
-        st.update("arith-make-feasible", lp().settings().stats().m_make_feasible);
-        st.update("arith-max-columns", lp().settings().stats().m_max_cols);
-        st.update("arith-max-rows", lp().settings().stats().m_max_rows);
-        st.update("arith-gcd-calls", lp().settings().stats().m_gcd_calls);
-        st.update("arith-gcd-conflict", lp().settings().stats().m_gcd_conflicts);
-        st.update("arith-cube-calls", lp().settings().stats().m_cube_calls);
-        st.update("arith-cube-success", lp().settings().stats().m_cube_success);
-        st.update("arith-patches", lp().settings().stats().m_patches);
-        st.update("arith-patches-success", lp().settings().stats().m_patches_success);
-        st.update("arith-hnf-calls", lp().settings().stats().m_hnf_cutter_calls);
-        st.update("arith-horner-calls", lp().settings().stats().m_horner_calls);
-        st.update("arith-horner-conflicts", lp().settings().stats().m_horner_conflicts);
-        st.update("arith-horner-cross-nested-forms", lp().settings().stats().m_cross_nested_forms);
-        st.update("arith-grobner-calls", lp().settings().stats().m_grobner_calls);
-        st.update("arith-grobner-conflicts", lp().settings().stats().m_grobner_conflicts);
+        m_stats.collect_statistics(st);
+        lp().settings().stats().collect_statistics(st);
         if (m_nla) m_nla->collect_statistics(st);
-        st.update("arith-gomory-cuts", m_stats.m_gomory_cuts);
-        st.update("arith-assume-eqs", m_stats.m_assume_eqs);
-        st.update("arith-branch", m_stats.m_branch);
-        st.update("arith-cheap-eqs", lp().settings().stats().m_cheap_eqs);
     }        
 
     /*
diff --git a/src/util/trail.h b/src/util/trail.h
index a2ce7a92a..1cdd07ef9 100644
--- a/src/util/trail.h
+++ b/src/util/trail.h
@@ -317,6 +317,28 @@ public:
     }
 };
 
+
+template<typename Ctx, typename T, bool CallDestructors = true>
+class history_trail : public trail<Ctx> {
+    vector<T, CallDestructors>& m_dst;
+    unsigned                     m_idx;
+    vector<T, CallDestructors>& m_hist;
+public:
+    history_trail(vector<T, CallDestructors>& v, unsigned idx, vector<T, CallDestructors>& hist) :
+        m_dst(v),
+        m_idx(idx),
+        m_hist(hist) {}
+    
+    ~history_trail() override {
+    }
+    
+    void undo(Ctx& ctx) override {
+        m_dst[m_idx] = m_hist.back();
+        m_hist.pop_back();
+    }
+};
+
+
 template<typename Ctx, typename T>
 class new_obj_trail : public trail<Ctx> {
     T * m_obj;