Introduce new monomials in Horner when shared factors have bounded linear combinations

When solve-eqs eliminates a variable x (= a - b) that appears as a factor
in a nonlinear product x*y, the product splits into a*y - b*y. The NLA
solver then reasons about a*y and b*y independently, losing the tight
bounds that x had. This can cause severe performance degradation (e.g.,
timeout on a QF_UFNIA verification condition that solves in 3s without
solve-eqs).

The Horner module's cross-nested factoring already recovers the factored
form y*(a-b), and interval_from_term (fixed in the previous commit) finds
the LP column for (a-b) with its tight bounds. However, only Horner's
zero-exclusion check used this — the rest of the NLA solver (order lemmas,
tangent planes, bounds propagation) continued reasoning about the split
monomials independently.

This commit adds a new mechanism: when Horner discovers that a linear
sub-expression maps to a bounded LP column, it introduces a new monomial
pairing that column with the shared factor. For example, if y*(a-b) is
discovered and (a-b) maps to LP column j with bounds [L,U], we create a
new monomial m := y*j via add_mul_def and assert the equality
m = a*y - b*y via literals. This allows all NLA modules to generate
lemmas using j's tight bounds.

The feature is gated by smt.arith.nl.horner_max_new_monomials (default 2,
0 to disable). On the motivating benchmark, this changes
simplify+propagate-values+solve-eqs+smt from timeout (30s) to UNSAT in
~15s with no regressions on other configurations.

Files changed:
- horner.cpp: introduce_monomials_from_term_columns() and find_binary_monic()
- horner.h: m_introduced_monomials dedup set
- nla_intervals.cpp/h: m_term_columns to record interval_from_term discoveries
- smt_params_helper.pyg: arith.nl.horner_max_new_monomials parameter

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

src/math/lp/horner.cpp

@@ -99,12 +99,122 @@ bool horner::lemmas_on_row(const T& row) {
 
 }
 
+// Find the binary monomial y*v in emonics, return its variable or null_lpvar.
+static lpvar find_binary_monic(emonics const& emons, lpvar y, lpvar v) {
+    if (!emons.is_used_in_monic(v))
+        return null_lpvar;
+    for (auto const& m : emons.get_use_list(v))
+        if (m.size() == 2 && (m.vars()[0] == y || m.vars()[1] == y))
+            return m.var();
+    return null_lpvar;
+}
+
+// Recover named intermediates destroyed by solve-eqs.
+//
+// When solve-eqs eliminates x = sum(c_i * v_i), product x*y splits into
+// sum(c_i * v_i * y). Horner's cross-nested factoring recovers y*sum(c_i*v_i)
+// and interval_from_term finds the LP term column tc := sum(c_i * v_i)
+// with tight bounds [L, U].
+//
+// We do two things:
+// 1. Create monomial m := y*tc, add LP row m - sum(c_i * mon(y,v_i)) = 0
+// 2. If variable x with monomial x*y exists and val(x) = val(tc),
+//    propagate: tc in [L,U] => x in [L,U]
+void horner::introduce_monomials_from_term_columns() {
+    if (c().params().arith_nl_horner_max_new_monomials() == 0)
+        return;
+    auto const& term_cols = c().m_intervals.term_columns();
+    if (term_cols.empty())
+        return;
+
+    unsigned added = 0;
+    for (lpvar tc : term_cols) {
+        if (!c().lra.column_has_lower_bound(tc) || !c().lra.column_has_upper_bound(tc))
+            continue;
+        if (!c().lra.column_has_term(tc))
+            continue;
+
+        auto const& term = c().lra.get_term(tc);
+
+        for (auto const& ti : term) {
+            lpvar vi = ti.j();
+            if (!c().m_emons.is_used_in_monic(vi))
+                continue;
+
+            for (auto const& m : c().m_emons.get_use_list(vi)) {
+                if (m.size() != 2)
+                    continue;
+                lpvar y = (m.vars()[0] == vi) ? m.vars()[1] : m.vars()[0];
+
+                auto key = std::make_pair(std::min(y, tc), std::max(y, tc));
+                if (m_introduced_monomials.contains(key))
+                    continue;
+
+                // Check that mon(y, v_i) exists for every v_i in tc
+                lp::lar_term eq_term;
+                bool complete = true;
+                for (auto const& tj : term) {
+                    lpvar yv = find_binary_monic(c().m_emons, y, tj.j());
+                    if (yv == null_lpvar) { complete = false; break; }
+                    eq_term.add_monomial(-tj.coeff(), yv);
+                }
+                if (!complete)
+                    continue;
+
+                m_introduced_monomials.push_back(key);
+
+                // (1) m := y * tc, with LP row: m - sum(c_i * mon(y,v_i)) = 0
+                lpvar factors[2] = { y, tc };
+                lpvar new_mon = c().add_mul_def(2, factors);
+                eq_term.add_monomial(rational::one(), new_mon);
+                lp::lpvar eq_col = c().lra.add_term(eq_term.coeffs_as_vector(), UINT_MAX);
+                c().lra.update_column_type_and_bound(eq_col, llc::EQ, rational::zero(), nullptr);
+                c().m_check_feasible = true;
+                TRACE(nla_solver, tout << "introduced monomial j" << new_mon
+                      << " := j" << y << " * j" << tc << "\n";);
+
+                // (2) Propagate tc's bounds to variable x where mon(x, y) exists
+                //     and val(x) = val(tc), i.e., x equals tc in current model.
+                for (auto const& m2 : c().m_emons.get_use_list(y)) {
+                    if (m2.size() != 2)
+                        continue;
+                    lpvar x = (m2.vars()[0] == y) ? m2.vars()[1] : m2.vars()[0];
+                    if (x == tc || c().lra.column_has_term(x))
+                        continue;
+                    if (c().lra.get_column_value(x) != c().lra.get_column_value(tc))
+                        continue;
+                    if (c().lra.column_has_lower_bound(tc)) {
+                        c().lra.update_column_type_and_bound(
+                            x, llc::GE, c().lra.get_lower_bound(tc).x,
+                            c().lra.get_column_lower_bound_witness(tc));
+                        c().m_check_feasible = true;
+                        TRACE(nla_solver, tout << "bound j" << x << " >= "
+                              << c().lra.get_lower_bound(tc).x << " from j" << tc << "\n";);
+                    }
+                    if (c().lra.column_has_upper_bound(tc)) {
+                        c().lra.update_column_type_and_bound(
+                            x, llc::LE, c().lra.get_upper_bound(tc).x,
+                            c().lra.get_column_upper_bound_witness(tc));
+                        c().m_check_feasible = true;
+                        TRACE(nla_solver, tout << "bound j" << x << " <= "
+                              << c().lra.get_upper_bound(tc).x << " from j" << tc << "\n";);
+                    }
+                }
+
+                if (++added >= c().params().arith_nl_horner_max_new_monomials())
+                    return;
+            }
+        }
+    }
+}
+
 bool horner::horner_lemmas() {
     if (!c().params().arith_nl_horner()) {
         TRACE(nla_solver, tout << "not generating horner lemmas\n";);
         return false;
     }
     c().lp_settings().stats().m_horner_calls++;
+    c().m_intervals.clear_term_columns();
     const auto& matrix = c().lra.A_r();
     // choose only rows that depend on m_to_refine variables
     std::set<unsigned> rows_to_check;
@@ -129,6 +239,10 @@ bool horner::horner_lemmas() {
             conflict = true;
         }
     }
+
+    if (!conflict)
+        introduce_monomials_from_term_columns();
+
     return conflict;
 }
 }

src/math/lp/horner.h

@@ -31,7 +31,8 @@ class core;
 
 class horner : common {
     nex_creator::sum_factory  m_row_sum;
-    unsigned         m_row_index;                      
+    unsigned         m_row_index;
+    svector<std::pair<lpvar, lpvar>> m_introduced_monomials; // track what we've already created                      
 public:
     typedef intervals::interval interv;
     horner(core *core);
@@ -49,5 +50,6 @@ public:
     template <typename T> // T has an iterator of (coeff(), var())
     bool row_has_monomial_to_refine(const T&) const;
     bool interval_from_term_with_deps(const nex* e, intervals::interval&) const;
+    void introduce_monomials_from_term_columns();
 }; // end of horner
 }

src/math/lp/nla_intervals.cpp

@@ -292,6 +292,7 @@ bool intervals::interval_from_term(const nex& e, scoped_dep_interval& i) {
     if (j + 1 == 0)
         return false;
 
+    m_term_columns.push_back(j);
     set_var_interval<wd>(j, i);
     interval bi;
     m_dep_intervals.mul<wd>(a, i, bi);

src/math/lp/nla_intervals.h

@@ -21,7 +21,8 @@ class core;
 class intervals {
     mutable dep_intervals     m_dep_intervals;
     core*                     m_core;
-    
+    svector<lpvar>            m_term_columns; // LP columns found by interval_from_term
+
 public:
     typedef dep_intervals::interval interval;
 private:
@@ -88,5 +89,8 @@ public:
     std::ostream& display_separating_interval(std::ostream& out, const nex*n, const scoped_dep_interval& interv_wd, u_dependency* initial_deps);
     bool conflict_u_l(const interval& a, const interval& b) const;
 
+    void clear_term_columns() { m_term_columns.clear(); }
+    const svector<lpvar>& term_columns() const { return m_term_columns; }
+
 }; // end of intervals
 } // end of namespace nla

src/params/smt_params_helper.pyg

@@ -91,6 +91,7 @@ def_module_params(module_name='smt',
                           ('arith.nl.propagate_linear_monomials', BOOL, True, 'propagate linear monomials'),
                           ('arith.nl.optimize_bounds', BOOL, True, 'enable bounds optimization'),
                           ('arith.nl.cross_nested', BOOL, True, 'enable cross-nested consistency checking'),
+                          ('arith.nl.horner_max_new_monomials', UINT, 2, 'maximum number of new monomials introduced per Horner round when a shared factor with a bounded linear combination is discovered (0 to disable)'),
                           ('arith.nl.log', BOOL, False, 'Log lemmas sent to nra solver'),
                           ('arith.propagate_eqs', BOOL, True, 'propagate (cheap) equalities'),
                           ('arith.epsilon', DOUBLE, 1.0, 'initial value of epsilon used for model generation of infinitesimals'),