Add mod-factor-propagation lemma to NLA divisions solver

When a monic x*y has a factor x with mod(x, p) = 0 (fixed), propagate
mod(x*y, p) = 0. This enables Z3 to prove divisibility properties like
x mod p = 0 => (x*y) mod p = 0, which previously timed out even for
p = 2. The lemma fires in the NLA divisions check and allows Gröbner
basis and LIA to subsequently derive distributivity of div over addition.

Extends division tuples from (q, x, y) to (q, x, y, r) to track the
mod lpvar. Also registers bounded divisions from the mod internalization
path in theory_lra, not just the idiv path.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

src/math/lp/nla_core.h

@@ -215,9 +215,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 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 add_idivision(lpvar q, lpvar x, lpvar y, lpvar r) { m_divisions.add_idivision(q, x, y, r); }
+    void add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r) { m_divisions.add_rdivision(q, x, y, r); }
+    void add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r) { m_divisions.add_bounded_division(q, x, y, r); }
     void set_add_mul_def_hook(std::function<lpvar(unsigned, lpvar const*)> const& f) { m_add_mul_def_hook = f; }
     lpvar add_mul_def(unsigned sz, lpvar const* vs) { SASSERT(m_add_mul_def_hook); lpvar v = m_add_mul_def_hook(sz, vs); add_monic(v, sz, vs); return v; }
 

src/math/lp/nla_divisions.cpp

@@ -18,26 +18,26 @@ Description:
 
 namespace nla {
 
-    void divisions::add_idivision(lpvar q, lpvar x, lpvar y) {
-        if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
+    void divisions::add_idivision(lpvar q, lpvar x, lpvar y, lpvar r) {
+        if (x == null_lpvar || y == null_lpvar || q == null_lpvar || r == null_lpvar)
             return;
-        m_idivisions.push_back({q, x, y});
+        m_idivisions.push_back({q, x, y, r});
         m_core.trail().push(push_back_vector(m_idivisions));
     }
 
-    void divisions::add_rdivision(lpvar q, lpvar x, lpvar y) {
-        if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
+    void divisions::add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r) {
+        if (x == null_lpvar || y == null_lpvar || q == null_lpvar || r == null_lpvar)
             return;
-        m_rdivisions.push_back({ q, x, y });
+        m_rdivisions.push_back({ q, x, y, r });
         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)
+    void divisions::add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r) {
+        if (x == null_lpvar || y == null_lpvar || q == null_lpvar || r == null_lpvar)
             return;
         if (m_core.lra.column_has_term(x) || m_core.lra.column_has_term(y) ||  m_core.lra.column_has_term(q))
             return;
-        m_bounded_divisions.push_back({ q, x, y });
+        m_bounded_divisions.push_back({ q, x, y, r });
         m_core.trail().push(push_back_vector(m_bounded_divisions));
     }
 
@@ -111,7 +111,7 @@ namespace nla {
                 return false;
         };
 
-        for (auto const & [r, x, y] : m_idivisions) {
+        for (auto const & [r, x, y, md] : m_idivisions) {
             if (!c.is_relevant(r))
                 continue;
             auto xval = c.val(x);
@@ -120,7 +120,7 @@ namespace nla {
             // idiv semantics
             if (!xval.is_int() || !yval.is_int() || yval == 0 || rval == div(xval, yval))
                 continue;
-            for (auto const& [q2, x2, y2] : m_idivisions) {
+            for (auto const& [q2, x2, y2, md2] : m_idivisions) {
                 if (q2 == r)
                     continue;
                 if (!c.is_relevant(q2))
@@ -133,7 +133,7 @@ namespace nla {
             }
         }
 
-        for (auto const& [r, x, y] : m_rdivisions) {
+        for (auto const& [r, x, y, md] : m_rdivisions) {
             if (!c.is_relevant(r))
                 continue;
             auto xval = c.val(x);
@@ -142,7 +142,7 @@ namespace nla {
             // / semantics
             if (yval == 0 || rval == xval / yval)
                 continue;
-            for (auto const& [q2, x2, y2] : m_rdivisions) {
+            for (auto const& [q2, x2, y2, md2] : m_rdivisions) {
                 if (q2 == r)
                     continue;
                 if (!c.is_relevant(q2))
@@ -154,7 +154,8 @@ namespace nla {
                     return;
             }
         }
-        
+
+        check_mod_mult();
     }
 
     // if p is bounded, q a value, r = eval(p):
@@ -163,11 +164,11 @@ namespace nla {
 
     void divisions::check_bounded_divisions() {
         core& c = m_core;
-        unsigned offset = c.random(), sz = m_bounded_divisions.size();        
+        unsigned offset = c.random(), sz = m_bounded_divisions.size();
 
         for (unsigned j = 0; j < sz; ++j) {
             unsigned i = (offset + j) % sz;
-            auto [q, x, y] = m_bounded_divisions[i];
+            auto [q, x, y, r] = m_bounded_divisions[i];
             if (!c.is_relevant(q))
                 continue;
             auto xv = c.val(x);
@@ -188,9 +189,9 @@ namespace nla {
             rational lo = yv * div_v;
             if (xv > hi) {
                 lemma_builder 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);  
+                lemma |= ineq(y, llc::NE, yv);
+                lemma |= ineq(x, llc::GT, hi);
+                lemma |= ineq(q, llc::LE, div_v);
                 return;
             }
             if (xv < lo) {
@@ -201,5 +202,45 @@ namespace nla {
                 return;
             }
         }
-    }    
+    }
+
+    // mod(factor, p) = 0 => mod(factor * k, p) = 0
+    // For each division (q, x, y, r) where x is a monic m = f1 * f2 * ... * fk,
+    // if some factor fi has mod(fi, p) = 0 (fixed), then mod(x, p) = 0.
+    void divisions::check_mod_mult() {
+        core& c = m_core;
+        unsigned offset = c.random(), sz = m_bounded_divisions.size();
+
+        for (unsigned j = 0; j < sz; ++j) {
+            unsigned i = (offset + j) % sz;
+            auto [q, x, y, r] = m_bounded_divisions[i];
+            if (!c.is_relevant(q))
+                continue;
+            if (c.var_is_fixed_to_zero(r))
+                continue;
+            if (c.val(r).is_zero())
+                continue;
+            if (!c.is_monic_var(x))
+                continue;
+            auto yv = c.val(y);
+            if (yv <= 0 || !yv.is_int())
+                continue;
+            auto const& m = c.emons()[x];
+            for (lpvar f : m.vars()) {
+                for (auto const& [q2, x2, y2, r2] : m_bounded_divisions) {
+                    if (x2 != f)
+                        continue;
+                    if (c.val(y2) != yv)
+                        continue;
+                    if (!c.var_is_fixed_to_zero(r2))
+                        continue;
+                    // mod(factor, p) = 0 => mod(product, p) = 0
+                    lemma_builder lemma(c, "mod(factor, p) = 0 => mod(factor * k, p) = 0");
+                    lemma |= ineq(r2, llc::NE, 0);
+                    lemma |= ineq(r, llc::EQ, 0);
+                    return;
+                }
+            }
+        }
+    }
 }

src/math/lp/nla_divisions.h

@@ -22,16 +22,17 @@ namespace nla {
     
     class divisions {
         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>> m_bounded_divisions;
-        
+        vector<std::tuple<lpvar, lpvar, lpvar, lpvar>> m_idivisions;
+        vector<std::tuple<lpvar, lpvar, lpvar, lpvar>> m_rdivisions;
+        vector<std::tuple<lpvar, 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 x, lpvar y);
+        void add_idivision(lpvar q, lpvar x, lpvar y, lpvar r);
+        void add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r);
+        void add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r);
         void check();
         void check_bounded_divisions();
+        void check_mod_mult();
     };
 }

src/math/lp/nla_solver.cpp

@@ -20,16 +20,16 @@ namespace nla {
         m_core->add_monic(v, sz, vs);
     }
 
-    void solver::add_idivision(lpvar q, lpvar x, lpvar y) {
-        m_core->add_idivision(q, x, y);
+    void solver::add_idivision(lpvar q, lpvar x, lpvar y, lpvar r) {
+        m_core->add_idivision(q, x, y, r);
     }
 
-    void solver::add_rdivision(lpvar q, lpvar x, lpvar y) {
-        m_core->add_rdivision(q, x, y);
+    void solver::add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r) {
+        m_core->add_rdivision(q, x, y, r);
     }
 
-    void solver::add_bounded_division(lpvar q, lpvar x, lpvar y) {
-        m_core->add_bounded_division(q, x, y);
+    void solver::add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r) {
+        m_core->add_bounded_division(q, x, y, r);
     }
 
     void solver::set_relevant(std::function<bool(lpvar)>& is_relevant) {

src/math/lp/nla_solver.h

@@ -28,9 +28,9 @@ namespace nla {
         ~solver();
         const auto& monics_with_changed_bounds() const { return m_core->monics_with_changed_bounds(); }
         void add_monic(lpvar v, unsigned sz, lpvar const* vs);
-        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 add_idivision(lpvar q, lpvar x, lpvar y, lpvar r);
+        void add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r);
+        void add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r);
         void check_bounded_divisions();
         void set_relevant(std::function<bool(lpvar)>& is_relevant);
         void updt_params(params_ref const& p);

src/smt/theory_lra.cpp

@@ -449,25 +449,43 @@ class theory_lra::imp {
                         internalize_term(to_app(n));
                         internalize_term(to_app(n1));
                         internalize_term(to_app(n2));
+                        internalize_term(to_app(mod));
                         theory_var q = mk_var(n);
                         theory_var x = mk_var(n1);
                         theory_var y = mk_var(n2);
-                        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));
+                        theory_var rv = mk_var(mod);
+                        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), register_theory_var_in_lar_solver(rv));
                     }
                     if (a.is_numeral(n2) && a.is_bounded(n1)) {
                         ensure_nla();
                         internalize_term(to_app(n));
                         internalize_term(to_app(n1));
                         internalize_term(to_app(n2));
+                        internalize_term(to_app(mod));
                         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));
+                        theory_var rv = mk_var(mod);
+                        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), register_theory_var_in_lar_solver(rv));
                     }
                 }
                 else if (a.is_mod(n, n1, n2)) {
                     if (!a.is_numeral(n2, r) || r.is_zero()) found_underspecified(n);
-                    if (!ctx().relevancy()) mk_idiv_mod_axioms(n1, n2);    
+                    if (!ctx().relevancy()) mk_idiv_mod_axioms(n1, n2);
+                    if (a.is_numeral(n2) && !r.is_zero()) {
+                        ensure_nla();
+                        app_ref div(a.mk_idiv(n1, n2), m);
+                        ctx().internalize(div, false);
+                        internalize_term(to_app(div));
+                        internalize_term(to_app(n1));
+                        internalize_term(to_app(n2));
+                        internalize_term(t);
+                        theory_var q = mk_var(div);
+                        theory_var x = mk_var(n1);
+                        theory_var y = mk_var(n2);
+                        theory_var rv = mk_var(n);
+                        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), register_theory_var_in_lar_solver(rv));
+                    }
                 }
                 else if (a.is_rem(n, n1, n2)) {
                     if (!a.is_numeral(n2, r) || r.is_zero()) found_underspecified(n);

src/test/CMakeLists.txt

@@ -87,6 +87,7 @@ add_executable(test-z3
   memory.cpp
   model2expr.cpp
   model_based_opt.cpp
+  mod_factor.cpp
   model_evaluator.cpp
   model_retrieval.cpp
   monomial_bounds.cpp

src/test/main.cpp

@@ -154,6 +154,7 @@
     X(hilbert_basis) \
     X(heap_trie) \
     X(karr) \
+    X(mod_factor) \
     X(no_overflow) \
     X(datalog_parser) \
     X_ARGV(datalog_parser_file) \

src/test/mod_factor.cpp

@@ -0,0 +1,65 @@
+/*++
+Copyright (c) 2025 Microsoft Corporation
+--*/
+
+#include "api/z3.h"
+#include "util/util.h"
+
+// x mod 7 = 0 & (x*y) mod 7 != 0 should be unsat
+// Exercises: mod internalization path (is_mod with numeric divisor)
+static void test_mod_factor_mod_path() {
+    Z3_config cfg = Z3_mk_config();
+    Z3_context ctx = Z3_mk_context(cfg);
+    Z3_solver s = Z3_mk_solver(ctx);
+    Z3_solver_inc_ref(ctx, s);
+    Z3_sort int_sort = Z3_mk_int_sort(ctx);
+    Z3_ast x = Z3_mk_const(ctx, Z3_mk_string_symbol(ctx, "x"), int_sort);
+    Z3_ast y = Z3_mk_const(ctx, Z3_mk_string_symbol(ctx, "y"), int_sort);
+    Z3_ast seven = Z3_mk_int(ctx, 7, int_sort);
+    Z3_ast zero = Z3_mk_int(ctx, 0, int_sort);
+    Z3_solver_assert(ctx, s, Z3_mk_eq(ctx, Z3_mk_mod(ctx, x, seven), zero));
+    Z3_ast xy_args[] = {x, y};
+    Z3_ast xy = Z3_mk_mul(ctx, 2, xy_args);
+    Z3_solver_assert(ctx, s, Z3_mk_not(ctx, Z3_mk_eq(ctx, Z3_mk_mod(ctx, xy, seven), zero)));
+    ENSURE(Z3_solver_check(ctx, s) == Z3_L_FALSE);
+    Z3_solver_dec_ref(ctx, s);
+    Z3_del_config(cfg);
+    Z3_del_context(ctx);
+}
+
+// (x mod 100) mod 7 = 0 => ((x mod 100) * y) mod 7 = 0
+// Exercises: idiv internalization path (is_idiv + numeric divisor + bounded dividend)
+// because (x mod 100) is recognized as bounded by is_bounded()
+static void test_mod_factor_idiv_path() {
+    Z3_config cfg = Z3_mk_config();
+    Z3_context ctx = Z3_mk_context(cfg);
+    Z3_solver s = Z3_mk_solver(ctx);
+    Z3_solver_inc_ref(ctx, s);
+    Z3_sort int_sort = Z3_mk_int_sort(ctx);
+    Z3_ast x = Z3_mk_const(ctx, Z3_mk_string_symbol(ctx, "x"), int_sort);
+    Z3_ast y = Z3_mk_const(ctx, Z3_mk_string_symbol(ctx, "y"), int_sort);
+    Z3_ast seven = Z3_mk_int(ctx, 7, int_sort);
+    Z3_ast zero = Z3_mk_int(ctx, 0, int_sort);
+    Z3_ast hundred = Z3_mk_int(ctx, 100, int_sort);
+    // xm = x mod 100 (bounded by is_bounded)
+    Z3_ast xm = Z3_mk_mod(ctx, x, hundred);
+    // xm mod 7 = 0
+    Z3_solver_assert(ctx, s, Z3_mk_eq(ctx, Z3_mk_mod(ctx, xm, seven), zero));
+    // (xm * y) div 7 — triggers idiv internalization with bounded dividend
+    Z3_ast xm_y_args[] = {xm, y};
+    Z3_ast xm_y = Z3_mk_mul(ctx, 2, xm_y_args);
+    Z3_ast xm_y_div = Z3_mk_div(ctx, xm_y, seven);
+    // assert (xm * y) mod 7 != 0
+    Z3_solver_assert(ctx, s, Z3_mk_not(ctx, Z3_mk_eq(ctx, Z3_mk_mod(ctx, xm_y, seven), zero)));
+    // use div to keep it alive
+    Z3_solver_assert(ctx, s, Z3_mk_ge(ctx, xm_y_div, zero));
+    ENSURE(Z3_solver_check(ctx, s) == Z3_L_FALSE);
+    Z3_solver_dec_ref(ctx, s);
+    Z3_del_config(cfg);
+    Z3_del_context(ctx);
+}
+
+void tst_mod_factor() {
+    test_mod_factor_mod_path();
+    test_mod_factor_idiv_path();
+}