diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index e00bfd1b0..386ad296d 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -985,6 +985,9 @@ bool core::rm_check(const monic& rm) const {
     return check_monic(m_emons[rm.var()]);
 }
 
+bool core::has_relevant_monomial() const {
+    return any_of(emons(), [&](auto const& m) { return is_relevant(m.var()); });
+}
     
 bool core::find_bfc_to_refine_on_monic(const monic& m, factorization & bf) {
     for (auto f : factorization_factory_imp(m, *this)) {
@@ -1489,6 +1492,11 @@ lbool core::check_power(lpvar r, lpvar x, lpvar y, vector<lemma>& l_vec) {
     return m_powers.check(r, x, y, l_vec);
 }
 
+void core::check_bounded_divisions(vector<lemma>& l_vec) {
+    m_lemma_vec = &l_vec;
+    m_divisions.check_bounded_divisions();
+}
+
 lbool core::check(vector<lemma>& l_vec) {
     lp_settings().stats().m_nla_calls++;
     TRACE("nla_solver", tout << "calls = " << lp_settings().stats().m_nla_calls << "\n";);
@@ -1527,7 +1535,7 @@ lbool core::check(vector<lemma>& l_vec) {
         m_basics.basic_lemma(false);
 
     if (l_vec.empty() && !done())
-        m_divisions.check(l_vec);
+        m_divisions.check();
     
 #if 0
     if (l_vec.empty() && !done() && !run_horner) 
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index eedc671f5..016ff6e9d 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -112,6 +112,9 @@ class core {
     void check_weighted(unsigned sz, std::pair<unsigned, std::function<void(void)>>* checks);
 
 public:    
+    // constructor
+    core(lp::lar_solver& s, reslimit&);
+
     void insert_to_refine(lpvar j);
     void erase_from_to_refine(lpvar j);
     
@@ -120,9 +123,7 @@ public:
 
     void insert_to_active_var_set(unsigned j) const { m_active_var_set.insert(j); }    
 
-    void clear_active_var_set() const {
-        m_active_var_set.clear();
-    }
+    void clear_active_var_set() const { m_active_var_set.clear(); }
 
     void clear_and_resize_active_var_set() const {
         m_active_var_set.clear();
@@ -134,9 +135,9 @@ public:
     reslimit& reslim() { return m_reslim; }  
     emonics& emons() { return m_emons; }
     const emonics& emons() const { return m_emons; }
-    // constructor
-    core(lp::lar_solver& s, reslimit &);
-   
+
+    bool has_relevant_monomial() const;
+
     bool compare_holds(const rational& ls, llc cmp, const rational& rs) const;
     
     rational value(const lp::lar_term& r) const;
@@ -202,8 +203,9 @@ 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) { m_divisions.add_idivision(r, x, y); }
-    void add_rdivision(lpvar r, lpvar x, lpvar y) { m_divisions.add_rdivision(r, x, y); }
+    void add_idivision(lpvar q, lpvar x, lpvar y) { m_divisions.add_idivision(q, x, y); }
+    void add_rdivision(lpvar q, lpvar x, lpvar y) { m_divisions.add_rdivision(q, x, y); }
+    void add_bounded_division(lpvar q, lpvar x, lpvar y) { m_divisions.add_bounded_division(q, 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); }
@@ -381,6 +383,7 @@ public:
     
     lbool check(vector<lemma>& l_vec);
     lbool check_power(lpvar r, lpvar x, lpvar y, vector<lemma>& l_vec);
+    void check_bounded_divisions(vector<lemma>&);
 
     bool  no_lemmas_hold() const;
     
diff --git a/src/math/lp/nla_divisions.cpp b/src/math/lp/nla_divisions.cpp
index 4f5bd4e81..91b674d58 100644
--- a/src/math/lp/nla_divisions.cpp
+++ b/src/math/lp/nla_divisions.cpp
@@ -36,13 +36,22 @@ namespace nla {
         m_core.trail().push(push_back_vector(m_rdivisions));
     }
 
+    void divisions::add_bounded_division(lpvar q, lpvar x, lpvar y) {
+        if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
+            return;
+        if (lp::tv::is_term(x) || lp::tv::is_term(y) || lp::tv::is_term(q))
+            return;
+        m_bounded_divisions.push_back({ q, x, y });
+        m_core.trail().push(push_back_vector(m_bounded_divisions));
+    }
+
     typedef lp::lar_term term;
     
     // y1 >= y2 > 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) {
+    void divisions::check() {
         core& c = m_core;        
         if (c.use_nra_model()) 
             return;
@@ -155,99 +164,45 @@ namespace nla {
     // p <= q * div(r, q) + q - 1 => div(p, q) <= div(r, q)
     // p >= q * div(r, q)         => div(r, q) <= div(p, q)
 
-#if 0
-    bool check_idiv_bounds() {
-        if (m_idiv_terms.empty()) {
-            return true;
-        }
-        bool all_divs_valid = true;
-        unsigned count = 0;
-        unsigned offset = ctx().get_random_value();
-        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));
-            theory_var v1 = internalize_def(to_app(p));
+    void divisions::check_bounded_divisions() {
+        core& c = m_core;
+        unsigned offset = c.random(), sz = m_bounded_divisions.size();        
 
-            if (!is_registered_var(v1))
+        for (unsigned j = 0; j < sz; ++j) {
+            unsigned i = (offset + j) % sz;
+            auto [q, x, y] = m_bounded_divisions[i];
+            if (!c.is_relevant(q))
+                continue;
+            auto xv = c.val(x);
+            auto yv = c.val(y);
+            auto qv = c.val(q);
+            if (xv < 0 || !xv.is_int())
+                continue;
+            if (yv <= 0 || !yv.is_int())
+                continue;
+            if (qv == div(xv, yv))
                 continue;
-            lp::impq q1 = get_ivalue(v1);
-            rational q2;
 
-            if (!q1.x.is_int() || q1.x.is_neg() || !q1.y.is_zero()) {
-                // TBD
-                // q1 = 223/4, q2 = 2, r = 219/8 
-                // take ceil(q1), floor(q1), ceil(q2), floor(q2), for floor(q2) > 0
-                // then 
-                //      p/q <= ceil(q1)/floor(q2) => n <= div(ceil(q1), floor(q2))
-                //      p/q >= floor(q1)/ceil(q2) => n >= div(floor(q1), ceil(q2))
-                continue;
+            rational div_v = div(xv, yv);
+            // y = yv & x <= yv * div(xv, yv) + yv - 1 => div(x, y) <= div(xv, yv)
+            // y = yv & x >= y * div(xv, yv) => div(xv, yv) <= div(x, y)
+            rational mul(1);
+            rational hi = yv * div_v + yv - 1;
+            rational lo = yv * div_v;
+            if (xv > hi) {
+                new_lemma lemma(c, "y = yv & x <= yv * div(xv, yv) + yv - 1 => div(p, y) <= div(xv, yv)");
+                lemma |= ineq(y, llc::NE, yv); 
+                lemma |= ineq(x, llc::GT, hi); 
+                lemma |= ineq(q, llc::LE, div_v);  
+                return;
             }
-
-
-            if (a.is_numeral(q, q2) && q2.is_pos()) {
-                if (!a.is_bounded(n)) {
-                    TRACE("arith", tout << "unbounded " << expr_ref(n, m) << "\n";);
-                    continue;
-                }
-                if (!is_registered_var(v))
-                    continue;
-                lp::impq val_v = get_ivalue(v);
-                if (val_v.y.is_zero() && val_v.x == div(q1.x, q2))
-                    continue;
-
-                TRACE("arith", tout << get_value(v) << " != " << q1 << " div " << q2 << "\n";);
-                rational div_r = div(q1.x, q2);
-                // p <= q * div(q1, q) + q - 1 => div(p, q) <= div(q1, q2)
-                // p >= q * div(q1, q) => div(q1, q) <= div(p, q)
-                rational mul(1);
-                rational hi = q2 * div_r + q2 - 1;
-                rational lo = q2 * 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_q1 = mk_literal(a.mk_le(p, a.mk_numeral(hi, true)));
-                literal p_ge_q1 = 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)));
-                {
-                    scoped_trace_stream _sts(th, ~p_le_q1, n_le_div);
-                    mk_axiom(~p_le_q1, n_le_div);
-                }
-                {
-                    scoped_trace_stream _sts(th, ~p_ge_q1, n_ge_div);
-                    mk_axiom(~p_ge_q1, n_ge_div);
-                }
-
-                all_divs_valid = false;
-                ++count;
-
-
-                TRACE("arith",
-                    tout << q1 << " div " << q2 << "\n";
-                literal_vector lits;
-                lits.push_back(~p_le_q1);
-                lits.push_back(n_le_div);
-                ctx().display_literals_verbose(tout, lits) << "\n\n";
-                lits[0] = ~p_ge_q1;
-                lits[1] = n_ge_div;
-                ctx().display_literals_verbose(tout, lits) << "\n";);
-                continue;
+            if (xv < lo) {
+                new_lemma lemma(c, "y = yv & x >= yv * div(xv, yv) => div(xv, yv) <= div(x, y)");
+                lemma |= ineq(y, llc::NE, yv);
+                lemma |= ineq(x, llc::LT, lo);
+                lemma |= ineq(q, llc::GE, div_v);
+                return;
             }
         }
-
-        return all_divs_valid;
-    }
-#endif
-    
+    }    
 }
diff --git a/src/math/lp/nla_divisions.h b/src/math/lp/nla_divisions.h
index 81fb0574a..80bf5be4e 100644
--- a/src/math/lp/nla_divisions.h
+++ b/src/math/lp/nla_divisions.h
@@ -24,12 +24,14 @@ namespace nla {
         core& m_core;
         vector<std::tuple<lpvar, lpvar, lpvar>> m_idivisions;
         vector<std::tuple<lpvar, lpvar, lpvar>> m_rdivisions;
-        vector<std::tuple<lpvar, lpvar, lpvar, lpvar>> m_bounded_divisions;
+        vector<std::tuple<lpvar, lpvar, lpvar>> m_bounded_divisions;
+        
     public:
         divisions(core& c):m_core(c) {}
         void add_idivision(lpvar q, lpvar x, lpvar y);
         void add_rdivision(lpvar q, lpvar x, lpvar y);
-        void add_bounded_division(lpvar q, lpvar r, lpvar x, lpvar y);
-        void check(vector<lemma>&);
+        void add_bounded_division(lpvar q, lpvar x, lpvar y);
+        void check();
+        void check_bounded_divisions();
     };
 }
diff --git a/src/math/lp/nla_solver.cpp b/src/math/lp/nla_solver.cpp
index 50842f333..bd0f1953c 100644
--- a/src/math/lp/nla_solver.cpp
+++ b/src/math/lp/nla_solver.cpp
@@ -23,12 +23,16 @@ namespace nla {
         m_core->add_monic(v, sz, vs);
     }
 
-    void solver::add_idivision(lpvar r, lpvar x, lpvar y) {
-        m_core->add_idivision(r, x, y);
+    void solver::add_idivision(lpvar q, lpvar x, lpvar y) {
+        m_core->add_idivision(q, x, y);
     }
 
-    void solver::add_rdivision(lpvar r, lpvar x, lpvar y) {
-        m_core->add_rdivision(r, x, y);
+    void solver::add_rdivision(lpvar q, lpvar x, lpvar y) {
+        m_core->add_rdivision(q, x, y);
+    }
+
+    void solver::add_bounded_division(lpvar q, lpvar x, lpvar y) {
+        m_core->add_bounded_division(q, x, y);
     }
 
     void solver::set_relevant(std::function<bool(lpvar)>& is_relevant) {
@@ -39,7 +43,7 @@ namespace nla {
         return m_core->is_monic_var(v);
     }
     
-    bool solver::need_check() { return true; }
+    bool solver::need_check() { return m_core->has_relevant_monomial(); }
     
     lbool solver::check(vector<lemma>& l) {
         return m_core->check(l);
@@ -92,4 +96,8 @@ namespace nla {
         return m_core->check_power(r, x, y, lemmas);
     }
 
+    void solver::check_bounded_divisions(vector<lemma>& lemmas) {
+        m_core->check_bounded_divisions(lemmas);
+    }
+
 }
diff --git a/src/math/lp/nla_solver.h b/src/math/lp/nla_solver.h
index 6667dff31..d04ff8e51 100644
--- a/src/math/lp/nla_solver.h
+++ b/src/math/lp/nla_solver.h
@@ -29,8 +29,10 @@ namespace nla {
         ~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 add_idivision(lpvar q, lpvar x, lpvar y);
+        void add_rdivision(lpvar q, lpvar x, lpvar y);
+        void add_bounded_division(lpvar q, lpvar x, lpvar y);
+        void check_bounded_divisions(vector<lemma>&);
         void set_relevant(std::function<bool(lpvar)>& is_relevant);
         nla_settings& settings();
         void push();
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 0daab80c1..c5b706fde 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -62,7 +62,6 @@ class theory_lra::imp {
 
     struct scope {
         unsigned m_bounds_lim;
-        unsigned m_idiv_lim;
         unsigned m_asserted_qhead;            
         unsigned m_asserted_atoms_lim;
     };
@@ -161,7 +160,6 @@ class theory_lra::imp {
     svector<delayed_atom>  m_asserted_atoms;        
     ptr_vector<expr>       m_not_handled;
     ptr_vector<app>        m_underspecified;
-    ptr_vector<expr>       m_idiv_terms;
     vector<ptr_vector<api_bound> > m_use_list;        // bounds where variables are used.
 
     // attributes for incremental version:
@@ -436,19 +434,23 @@ class theory_lra::imp {
                 }
                 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);
                     ctx().internalize(mod, false);
                     if (ctx().relevancy()) ctx().add_relevancy_dependency(n, mod);
-#if 1
                     if (m_nla && !a.is_numeral(n2)) {
                         // shortcut to create non-linear division axioms.
-                        theory_var r = mk_var(n);
+                        theory_var q = 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
+                        m_nla->add_idivision(register_theory_var_in_lar_solver(q), register_theory_var_in_lar_solver(x), register_theory_var_in_lar_solver(y));
+                    }
+                    if (a.is_numeral(n2) && a.is_bounded(n1)) {
+                        ensure_nla();
+                        theory_var q = mk_var(n);
+                        theory_var x = mk_var(n1);
+                        theory_var y = mk_var(n2);
+                        m_nla->add_bounded_division(register_theory_var_in_lar_solver(q), register_theory_var_in_lar_solver(x), register_theory_var_in_lar_solver(y));
+                    }
                 }
                 else if (a.is_mod(n, n1, n2)) {
                     if (!a.is_numeral(n2, r) || r.is_zero()) found_underspecified(n);
@@ -1077,7 +1079,6 @@ public:
         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_atoms_lim = m_asserted_atoms.size();
         lp().push();
         if (m_nla)
@@ -1092,7 +1093,6 @@ public:
         }
         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_atoms.shrink(m_scopes[old_size].m_asserted_atoms_lim);
         m_asserted_qhead = m_scopes[old_size].m_asserted_qhead;
         m_scopes.resize(old_size);            
@@ -1480,7 +1480,7 @@ public:
     }
     
     void random_update() {
-        if (m_nla)
+        if (m_nla && m_nla->need_check())
             return;
         m_tmp_var_set.clear();
         m_tmp_var_set.resize(th.get_num_vars());
@@ -1799,96 +1799,13 @@ public:
      */
 
     bool check_idiv_bounds() {
-        if (m_idiv_terms.empty()) {
+        if (!m_nla)
             return true;
-        }
-        bool all_divs_valid = true;
-        unsigned count = 0;
-        unsigned offset = ctx().get_random_value(); 
-        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 (false && !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));
-            theory_var v1 = internalize_def(to_app(p));
-
-            if (!is_registered_var(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 (!is_registered_var(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)));
-                {
-                    scoped_trace_stream _sts(th, ~p_le_r1, n_le_div);
-                    mk_axiom(~p_le_r1, n_le_div);
-                }
-                {
-                    scoped_trace_stream _sts(th, ~p_ge_r1, n_ge_div);
-                    mk_axiom(~p_ge_r1, n_ge_div);
-                }
-
-                all_divs_valid = false;
-                ++count;
-
-
-                TRACE("arith",
-                      tout << r1 << " div " << r2 << "\n";
-                      literal_vector lits;
-                      lits.push_back(~p_le_r1);
-                      lits.push_back(n_le_div);
-                      ctx().display_literals_verbose(tout, lits) << "\n\n";
-                      lits[0] = ~p_ge_r1;
-                      lits[1] = n_ge_div;
-                      ctx().display_literals_verbose(tout, lits) << "\n";);                      
-                continue;
-            }
-        }
-        
-        return all_divs_valid;
+        m_nla_lemma_vector.reset();
+        m_nla->check_bounded_divisions(m_nla_lemma_vector);
+        for (auto & lemma : m_nla_lemma_vector)
+            false_case_of_check_nla(lemma);
+        return m_nla_lemma_vector.empty();         
     }
 
     expr_ref var2expr(lpvar v) {
@@ -2105,9 +2022,8 @@ public:
         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);
-            }
+            for (const nla::lemma & l : m_nla_lemma_vector) 
+                false_case_of_check_nla(l);            
             break;
         }
         case l_true:
@@ -2126,11 +2042,11 @@ public:
             TRACE("arith", tout << "canceled\n";);
             return l_undef;
         }
-        if (!m_nla) {
-            TRACE("arith", tout << "no nla\n";);
+        CTRACE("arith",!m_nla, tout << "no nla\n";);
+        if (!m_nla)            
+            return l_true;        
+        if (!m_nla->need_check()) 
             return l_true;
-        }
-        if (!m_nla->need_check()) return l_true;
         return check_nla_continue();
     }