diff --git a/src/math/lp/emonics.cpp b/src/math/lp/emonics.cpp
index 060c2b1b9..21bfe6792 100644
--- a/src/math/lp/emonics.cpp
+++ b/src/math/lp/emonics.cpp
@@ -40,43 +40,39 @@ void emonics::push() {
     TRACE("nla_solver_mons", display(tout << "push\n"););
     SASSERT(invariant());
     m_u_f_stack.push_scope();
-    m_lim.push_back(m_monics.size());
-    m_region.push_scope();
     m_ve.push();
     SASSERT(monics_are_canonized());
     SASSERT(invariant());
 }
 
+void emonics::pop_monic() {
+    m_ve.pop(1);
+    monic& m = m_monics.back();
+    TRACE("nla_solver_mons", display(tout << m << "\n"););
+    remove_cg_mon(m);
+    m_var2index[m.var()] = UINT_MAX;
+    do_canonize(m);
+    // variables in vs are in the same state as they were when add was called
+    lpvar last_var = UINT_MAX;
+    for (lpvar v : m.rvars()) {
+        if (v != last_var) {
+            remove_cell(m_use_lists[v]);
+            last_var = v;
+        }
+    }
+    m_ve.pop(1);
+    m_monics.pop_back();
+}
+
 void emonics::pop(unsigned n) {
     TRACE("nla_solver_mons", tout << "pop: " << n << "\n";);
     SASSERT(invariant());
-    for (unsigned j = 0; j < n; ++j) {
-        unsigned old_sz = m_lim[m_lim.size() - 1];
-        for (unsigned i = m_monics.size(); i-- > old_sz; ) {
-            m_ve.pop(1);
-            monic & m = m_monics[i];
-            TRACE("nla_solver_mons", display(tout << m << "\n"););
-            remove_cg_mon(m);
-            m_var2index[m.var()] = UINT_MAX;
-            do_canonize(m);
-            // variables in vs are in the same state as they were when add was called
-            lpvar last_var = UINT_MAX;
-            for (lpvar v : m.rvars()) {
-                if (v != last_var) {
-                    remove_cell(m_use_lists[v]);
-                    last_var = v;
-                }
-            }
-            m_ve.pop(1);
-        }
-        m_ve.pop(1);
-        m_monics.shrink(old_sz);
-        m_region.pop_scope(1);
-        m_lim.pop_back();
+    for (unsigned i = 0; i < n; ++i) {
         m_u_f_stack.pop_scope(1);
-        SASSERT(invariant());
-        SASSERT(monics_are_canonized());
+        m_ve.pop(1);
     }
+    SASSERT(invariant());
+    SASSERT(monics_are_canonized());
 }
 
 void emonics::remove_cell(head_tail& v) {
@@ -96,7 +92,7 @@ void emonics::remove_cell(head_tail& v) {
 void emonics::insert_cell(head_tail& v, unsigned mIndex) {
     cell*& cur_head = v.m_head;
     cell*& cur_tail = v.m_tail;
-    cell* new_head = new (m_region) cell(mIndex, cur_head);
+    cell* new_head = new (m_u_f_stack.get_region()) cell(mIndex, cur_head);
     cur_head = new_head;
     if (!cur_tail) cur_tail = new_head;
     cur_tail->m_next = new_head;
@@ -331,6 +327,14 @@ void emonics::add(lpvar v, unsigned sz, lpvar const* vs) {
     m_monics.push_back(monic(v, sz, vs, idx));
     do_canonize(m_monics.back());
 
+    class pop_mon : public trail {
+        emonics& p;
+    public:
+        pop_mon(emonics& p) :p(p) {}
+        void undo() override { p.pop_monic(); }
+    };
+    m_u_f_stack.push(pop_mon(*this));
+
     // variables in m_vs are canonical and sorted, 
     // so use last_var to skip duplicates, 
     // while updating use-lists
@@ -351,9 +355,8 @@ void emonics::add(lpvar v, unsigned sz, lpvar const* vs) {
 void emonics::do_canonize(monic & m) const {
     TRACE("nla_solver_mons", tout << m << "\n";);
     m.reset_rfields();
-    for (lpvar v : m.vars()) {
-        m.push_rvar(m_ve.find(v));
-    }
+    for (lpvar v : m.vars()) 
+        m.push_rvar(m_ve.find(v));    
     m.sort_rvars();
     TRACE("nla_solver_mons", tout << m << "\n";);
 }
@@ -365,40 +368,34 @@ bool emonics::is_canonized(const monic & m) const {
 }
 
 void emonics::ensure_canonized() {
-    for (auto & m : m_monics) {
-        do_canonize(m);
-    }
+    for (auto & m : m_monics) 
+        do_canonize(m);    
 }
 
 bool emonics::monics_are_canonized() const {
-    for (auto & m: m_monics) {
-        if (!is_canonized(m)) {
-            return false;
-        }
-    }
+    for (auto & m: m_monics) 
+        if (!is_canonized(m)) 
+            return false;        
     return true;
 }
 
 bool emonics::canonize_divides(monic& m, monic & n) const {
-    if (m.size() > n.size()) return false;
+    if (m.size() > n.size()) 
+        return false;
     unsigned ms = m.size(), ns = n.size();
     unsigned i = 0, j = 0;
     while (true) {
-        if (i == ms) {
-            return true;
-        }
-        else if (j == ns) {
-            return false;
-        }
+        if (i == ms) 
+            return true;        
+        else if (j == ns) 
+            return false;        
         else if (m.rvars()[i] == n.rvars()[j]) {
             ++i; ++j;
         }
-        else if (m.rvars()[i] < n.rvars()[j]) {
-            return false;
-        }
-        else {
-            ++j;
-        }
+        else if (m.rvars()[i] < n.rvars()[j]) 
+            return false;        
+        else 
+            ++j;        
     }
 }
 
diff --git a/src/math/lp/emonics.h b/src/math/lp/emonics.h
index 0d0c20fcb..d26d96ad6 100644
--- a/src/math/lp/emonics.h
+++ b/src/math/lp/emonics.h
@@ -87,15 +87,15 @@ class emonics {
     var_eqs<emonics>&            m_ve;
     mutable vector<monic>        m_monics;     // set of monics
     mutable unsigned_vector      m_var2index;     // var_mIndex -> mIndex
-    unsigned_vector              m_lim;           // backtracking point
     mutable unsigned             m_visited;       // timestamp of visited monics during pf_iterator
-    region                       m_region;        // region for allocating linked lists
     mutable svector<head_tail>   m_use_lists;     // use list of monics where variables occur.
     hash_canonical               m_cg_hash;
     eq_canonical                 m_cg_eq;
     map<lpvar, unsigned_vector, hash_canonical, eq_canonical> m_cg_table; // congruence (canonical) table.
 
 
+    void pop_monic();
+
     void inc_visited() const;
 
     void remove_cell(head_tail& v);
@@ -115,6 +115,8 @@ class emonics {
     std::ostream& display_use(std::ostream& out) const; 
     std::ostream& display_uf(std::ostream& out) const; 
     std::ostream& display(std::ostream& out, cell* c) const;
+
+
 public:
     unsigned number_of_monics() const { return m_monics.size(); }
     /**
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index 0ec7d42c0..e00bfd1b0 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -137,11 +137,6 @@ void core::add_monic(lpvar v, unsigned sz, lpvar const* vs) {
     }
     m_emons.add(v, m_add_buffer);
 }
-
-void core::add_idivision(lpvar r, lpvar x, lpvar y) {
-    m_divisions.add_idivision(r, x, y);
-}
-
     
 void core::push() {
     TRACE("nla_solver_verbose", tout << "\n";);
@@ -164,7 +159,13 @@ rational core::product_value(const monic& m) const {
 }
     
 // return true iff the monic value is equal to the product of the values of the factors
+//                 or if the variable associated with the monomial is not relevant.
 bool core::check_monic(const monic& m) const {    
+#if 0
+    // TODO test this
+    if (!is_relevant(m.var()))
+        return true;
+#endif
     SASSERT((!m_lar_solver.column_is_int(m.var())) || m_lar_solver.get_column_value(m.var()).is_int());
     bool ret = product_value(m) == m_lar_solver.get_column_value(m.var()).x; 
     CTRACE("nla_solver_check_monic", !ret, print_monic(m, tout) << '\n';);
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index c595a0b0f..eedc671f5 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -82,6 +82,7 @@ class core {
 
     lp::lar_solver&          m_lar_solver;
     reslimit&                m_reslim;
+    std::function<bool(lpvar)> m_relevant;
     vector<lemma> *          m_lemma_vec;
     lp::u_set                m_to_refine;
     tangents                 m_tangents;
@@ -201,9 +202,15 @@ public:
     void deregister_monic_from_tables(const monic & m, unsigned i);
 
     void add_monic(lpvar v, unsigned sz, lpvar const* vs);   
-    void add_idivision(lpvar r, lpvar x, lpvar y);
+    void add_idivision(lpvar r, lpvar x, lpvar y) { m_divisions.add_idivision(r, x, y); }
+    void add_rdivision(lpvar r, lpvar x, lpvar y) { m_divisions.add_rdivision(r, x, y); }
+
+    void set_relevant(std::function<bool(lpvar)>& is_relevant) { m_relevant = is_relevant; }
+    bool is_relevant(lpvar v) const { return !m_relevant || m_relevant(v); }
+
     void push();     
     void pop(unsigned n);
+
     trail_stack& trail() { return m_emons.get_trail_stack(); }
 
     rational mon_value_by_vars(unsigned i) const;
diff --git a/src/math/lp/nla_divisions.cpp b/src/math/lp/nla_divisions.cpp
index dd58fcd59..ce875b9c4 100644
--- a/src/math/lp/nla_divisions.cpp
+++ b/src/math/lp/nla_divisions.cpp
@@ -23,43 +23,128 @@ namespace nla {
         m_core.trail().push(push_back_vector(m_idivisions));
     }
 
+    void divisions::add_rdivision(lpvar r, lpvar x, lpvar y) {
+        m_rdivisions.push_back({ r, x, y });
+        m_core.trail().push(push_back_vector(m_rdivisions));
+    }
+
     typedef lp::lar_term term;
     
     // y1 >= y2 > 0 & x1 <= x2 => x1/y1 <= x2/y2
-    // y2 <= y1 < 0 & x1 >= x2 => x1/y1 <= x2/y2
+    // y2 <= y1 < 0 & x1 >= x2 >= 0 => x1/y1 <= x2/y2
+    // y2 <= y1 < 0 & x1 <= x2 <= 0 => x1/y1 >= x2/y2
+
     void divisions::check(vector<lemma>& lemmas) {
         core& c = m_core;        
         if (c.use_nra_model()) 
             return;
 
+        auto monotonicity1 = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& r1, auto& r1val,
+            auto x2, auto& x2val, auto y2, auto& y2val, auto& r2, auto& r2val) {
+                if (y1val >= y2val && y2val > 0 && x1val <= x2val && r1val > r2val) {
+                    new_lemma lemma(c, "y1 >= y2 > 0 & x1 <= x2 => x1/y1 <= x2/y2");
+                    lemma |= ineq(term(y1, rational(-1), y2), llc::LT, 0);
+                    lemma |= ineq(y2, llc::LE, 0);
+                    lemma |= ineq(term(x1, rational(-1), x2), llc::GT, 0);
+                    lemma |= ineq(term(r1, rational(-1), r2), llc::LE, 0);
+                    return true;
+                }
+                return false;
+        };
+
+        auto monotonicity2 = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& r1, auto& r1val,
+            auto x2, auto& x2val, auto y2, auto& y2val, auto& r2, auto& r2val) {
+                if (y2val <= y1val && y1val < 0 && x1val >= x2val && x2val >= 0 && r1val > r2val) {
+                    new_lemma lemma(c, "y2 <= y1 < 0 & x1 >= x2 >= 0 => x1/y1 <= x2/y2");
+                    lemma |= ineq(term(y1, rational(-1), y2), llc::LT, 0);
+                    lemma |= ineq(y1, llc::GE, 0);
+                    lemma |= ineq(term(x1, rational(-1), x2), llc::LT, 0);
+                    lemma |= ineq(x2, llc::LT, 0);
+                    lemma |= ineq(term(r1, rational(-1), r2), llc::LE, 0);
+                    return true;
+                }
+                return false;
+        };
+
+        auto monotonicity3 = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& r1, auto& r1val,
+            auto x2, auto& x2val, auto y2, auto& y2val, auto& r2, auto& r2val) {
+                if (y2val <= y1val && y1val < 0 && x1val <= x2val && x2val <= 0 && r1val < r2val) {
+                    new_lemma lemma(c, "y2 <= y1 < 0 & x1 <= x2 <= 0 => x1/y1 >= x2/y2");
+                    lemma |= ineq(term(y1, rational(-1), y2), llc::LT, 0);
+                    lemma |= ineq(y1, llc::GE, 0);
+                    lemma |= ineq(term(x1, rational(-1), x2), llc::GT, 0);
+                    lemma |= ineq(x2, llc::GT, 0);
+                    lemma |= ineq(term(r1, rational(-1), r2), llc::GE, 0);
+                    return true;
+                }
+                return false;
+        };
+
+        auto monotonicity = [&](auto x1, auto& x1val, auto y1, auto& y1val, auto& r1, auto& r1val,
+            auto x2, auto& x2val, auto y2, auto& y2val, auto& r2, auto& r2val) {
+                if (monotonicity1(x1, x1val, y1, y1val, r1, r1val, x2, x2val, y2, y2val, r2, r2val))
+                    return true;
+                if (monotonicity1(x2, x2val, y2, y2val, r2, r2val, x1, x1val, y1, y1val, r1, r1val))
+                    return true;
+                if (monotonicity2(x1, x1val, y1, y1val, r1, r1val, x2, x2val, y2, y2val, r2, r2val))
+                    return true;
+                if (monotonicity2(x2, x2val, y2, y2val, r2, r2val, x1, x1val, y1, y1val, r1, r1val))
+                    return true;
+                if (monotonicity3(x1, x1val, y1, y1val, r1, r1val, x2, x2val, y2, y2val, r2, r2val))
+                    return true;
+                if (monotonicity3(x2, x2val, y2, y2val, r2, r2val, x1, x1val, y1, y1val, r1, r1val))
+                    return true;
+                return false;
+        };
+
         for (auto const & [r, x, y] : m_idivisions) {
+            if (!c.is_relevant(r))
+                continue;
             auto xval = c.val(x);
             auto yval = c.val(y);
             auto rval = c.val(r);
-            if (!c.var_is_int(x)) 
-                continue;
-            if (yval == 0)
-                continue;
             // idiv semantics
-            if (rval == div(xval, yval))
+            if (!xval.is_int() || !yval.is_int() || yval == 0 || rval == div(xval, yval))
                 continue;
             for (auto const& [r2, x2, y2] : m_idivisions) {
                 if (r2 == r)
                     continue;
+                if (!c.is_relevant(r2))
+                    continue;
                 auto x2val = c.val(x2);
                 auto y2val = c.val(y2);
                 auto r2val = c.val(r2);
-                if (yval >= y2val && y2val > 0 && xval <= x2val && rval > r2val) {
-                    new_lemma lemma(c, "y1 >= y2 > 0 & x1 <= x2 => x1/y1 <= x2/y2");
-                    lemma |= ineq(term(y, rational(-1), y2), llc::LT, rational::zero());
-                    lemma |= ineq(y2, llc::LE, rational::zero());
-                    lemma |= ineq(term(x, rational(-1), x2), llc::GT, rational::zero());
-                    lemma |= ineq(term(r, rational(-1), r2), llc::LE, rational::zero());
+                if (monotonicity(x, xval, y, yval, r, rval, x2, x2val, y2, y2val, r2, r2val))
+                    return;
+            }
+        }
+
+        for (auto const& [r, x, y] : m_rdivisions) {
+            if (!c.is_relevant(r))
+                continue;
+            auto xval = c.val(x);
+            auto yval = c.val(y);
+            auto rval = c.val(r);
+            // / semantics
+            if (yval == 0 || rval == xval / yval)
+                continue;
+            for (auto const& [r2, x2, y2] : m_rdivisions) {
+                if (r2 == r)
+                    continue;
+                if (!c.is_relevant(r2))
+                    continue;
+                auto x2val = c.val(x2);
+                auto y2val = c.val(y2);
+                auto r2val = c.val(r2);
+                if (monotonicity(x, xval, y, yval, r, rval, x2, x2val, y2, y2val, r2, r2val))
                     return;
-                }
             }
         }
         
     }
+
+    // if p is bounded, q a value, r = eval(p):
+    // p <= q * div(r, q) + q - 1 => div(p, q) <= div(r, q)
+    // p >= q * div(r, q)         => div(r, q) <= div(p, q)
     
 }
diff --git a/src/math/lp/nla_divisions.h b/src/math/lp/nla_divisions.h
index 8eb5c676b..56056e1f4 100644
--- a/src/math/lp/nla_divisions.h
+++ b/src/math/lp/nla_divisions.h
@@ -23,9 +23,11 @@ namespace nla {
     class divisions {
         core& m_core;
         vector<std::tuple<lpvar, lpvar, lpvar>> m_idivisions;
+        vector<std::tuple<lpvar, lpvar, lpvar>> m_rdivisions;
     public:
         divisions(core& c):m_core(c) {}
         void add_idivision(lpvar r, lpvar x, lpvar y);
+        void add_rdivision(lpvar r, lpvar x, lpvar y);
         void check(vector<lemma>&);
     };
 }
diff --git a/src/math/lp/nla_solver.cpp b/src/math/lp/nla_solver.cpp
index eb7e5283d..50842f333 100644
--- a/src/math/lp/nla_solver.cpp
+++ b/src/math/lp/nla_solver.cpp
@@ -26,6 +26,14 @@ namespace nla {
     void solver::add_idivision(lpvar r, lpvar x, lpvar y) {
         m_core->add_idivision(r, x, y);
     }
+
+    void solver::add_rdivision(lpvar r, lpvar x, lpvar y) {
+        m_core->add_rdivision(r, x, y);
+    }
+
+    void solver::set_relevant(std::function<bool(lpvar)>& is_relevant) {
+        m_core->set_relevant(is_relevant);
+    }
     
     bool solver::is_monic_var(lpvar v) const {
         return m_core->is_monic_var(v);
diff --git a/src/math/lp/nla_solver.h b/src/math/lp/nla_solver.h
index 9379939db..6667dff31 100644
--- a/src/math/lp/nla_solver.h
+++ b/src/math/lp/nla_solver.h
@@ -24,10 +24,14 @@ namespace nla {
     class solver {
         core* m_core;
     public:
-        void add_monic(lpvar v, unsigned sz, lpvar const* vs);    
-        void add_idivision(lpvar r, lpvar x, lpvar y);
+
         solver(lp::lar_solver& s, reslimit& limit);
         ~solver();
+
+        void add_monic(lpvar v, unsigned sz, lpvar const* vs);
+        void add_idivision(lpvar r, lpvar x, lpvar y);
+        void add_rdivision(lpvar r, lpvar x, lpvar y);
+        void set_relevant(std::function<bool(lpvar)>& is_relevant);
         nla_settings& settings();
         void push();
         void pop(unsigned scopes);
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 880026610..8a9f41e87 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -275,6 +275,11 @@ class theory_lra::imp {
                 (void)_s;
                 m_nla->push();
             }
+            std::function<bool(lpvar)> is_relevant = [&](lpvar v) {
+                theory_var u = lp().local_to_external(v);
+                return ctx().is_relevant(th.get_enode(u));
+            };
+            m_nla->set_relevant(is_relevant);
             smt_params_helper prms(ctx().get_params());
             m_nla->settings().run_order =                   prms.arith_nl_order();
             m_nla->settings().run_tangents =                prms.arith_nl_tangents();
@@ -435,12 +440,14 @@ class theory_lra::imp {
                     app_ref mod(a.mk_mod(n1, n2), m);
                     ctx().internalize(mod, false);
                     if (ctx().relevancy()) ctx().add_relevancy_dependency(n, mod);
-#if 0
-                    // shortcut to create non-linear division axioms.
-                    theory_var r = mk_var(n);
-                    theory_var x = mk_var(n1);
-                    theory_var y = mk_var(n2);
-                    m_nla->add_idivision(get_lpvar(n), get_lpvar(n1), get_lpvar(n2));
+#if 1
+                    if (m_nla && !a.is_numeral(n2)) {
+                        // shortcut to create non-linear division axioms.
+                        theory_var r = mk_var(n);
+                        theory_var x = mk_var(n1);
+                        theory_var y = mk_var(n2);
+                        m_nla->add_idivision(get_lpvar(n), get_lpvar(n1), get_lpvar(n2));
+                }
 #endif
                 }
                 else if (a.is_mod(n, n1, n2)) {
@@ -1801,6 +1808,8 @@ public:
         for (unsigned j = 0; j < m_idiv_terms.size(); ++j) {
             unsigned i =  (offset + j) % m_idiv_terms.size();
             expr* n = m_idiv_terms[i];
+            if (!ctx().is_relevant(n))
+                continue;
             expr* p = nullptr, *q = nullptr;
             VERIFY(a.is_idiv(n, p, q));
             theory_var v = internalize_def(to_app(n));