diff --git a/src/math/lp/horner.cpp b/src/math/lp/horner.cpp index fffb4357c..33fa4e349 100644 --- a/src/math/lp/horner.cpp +++ b/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 rows_to_check; @@ -129,6 +239,10 @@ bool horner::horner_lemmas() { conflict = true; } } + + if (!conflict) + introduce_monomials_from_term_columns(); + return conflict; } } diff --git a/src/math/lp/horner.h b/src/math/lp/horner.h index 9d530acee..40f5b2153 100644 --- a/src/math/lp/horner.h +++ b/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> m_introduced_monomials; // track what we've already created public: typedef intervals::interval interv; horner(core *core); @@ -49,5 +50,6 @@ public: template // 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 } diff --git a/src/math/lp/nla_intervals.cpp b/src/math/lp/nla_intervals.cpp index f6c15bd8b..79740fd13 100644 --- a/src/math/lp/nla_intervals.cpp +++ b/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(j, i); interval bi; m_dep_intervals.mul(a, i, bi); diff --git a/src/math/lp/nla_intervals.h b/src/math/lp/nla_intervals.h index 514049aca..48b005adf 100644 --- a/src/math/lp/nla_intervals.h +++ b/src/math/lp/nla_intervals.h @@ -21,7 +21,8 @@ class core; class intervals { mutable dep_intervals m_dep_intervals; core* m_core; - + svector 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& term_columns() const { return m_term_columns; } + }; // end of intervals } // end of namespace nla diff --git a/src/params/smt_params_helper.pyg b/src/params/smt_params_helper.pyg index 708d88bb6..23be55350 100644 --- a/src/params/smt_params_helper.pyg +++ b/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'),