From aba41d026f61e4fb21c70bc57792e8586f30b3cf Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Thu, 23 Jul 2026 09:00:46 -0700 Subject: [PATCH] =?UTF-8?q?nla:=20add=20LP-based=20nonlinear=20bound=20opt?= =?UTF-8?q?imization=20for=20cross-nested=20confl=E2=80=A6=20(#10180)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit …icts Add core::optimize_nl_bounds() (gated by arith.nl.optimize_bounds) which runs LP max/min over monomial leaf variables inside core::propagate(), analogous to solver=2's max_min_nl_vars, so nla (arith.solver=6) can detect cross-nested conflicts previously missed. Collect improved bounds first, then apply them and re-establish feasibility once; reconcile the core solver via find_feasible_solution before the raw maximize solves to preserve inf_heap_is_correct(). Skip null witnesses in get_dependencies_of_maximum for implied/unconditional bounds. On FStar-UInt128-divergence solver=6 this yields unsat in 2 final-checks, seed-insensitive (seeds 1-10). Copilot-Session: ac36bb84-de91-4e6c-86df-6008c7396ceb --------- Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> Copilot-Session: ac36bb84-de91-4e6c-86df-6008c7396ceb Copilot-Session: 96a14756-2ffe-4cc3-87e7-49fda1b6113a --- src/math/interval/dep_intervals.h | 12 + src/math/lp/horner.cpp | 6 + src/math/lp/lar_solver.cpp | 1 - src/math/lp/monomial_bounds.cpp | 479 ++++++++++++++---------------- src/math/lp/monomial_bounds.h | 42 +-- src/math/lp/nla_core.cpp | 91 +++++- src/math/lp/nla_core.h | 5 +- src/math/lp/nla_solver.cpp | 4 +- src/math/lp/nla_solver.h | 2 +- src/params/smt_params_helper.pyg | 1 + src/smt/theory_lra.cpp | 11 +- src/util/mpz.cpp | 6 +- src/util/mpz.h | 2 +- 13 files changed, 370 insertions(+), 292 deletions(-) diff --git a/src/math/interval/dep_intervals.h b/src/math/interval/dep_intervals.h index eab62821af..aebdad4ea3 100644 --- a/src/math/interval/dep_intervals.h +++ b/src/math/interval/dep_intervals.h @@ -92,6 +92,12 @@ private: unsynch_mpq_manager::is_zero(a.m_lower) && unsynch_mpq_manager::is_zero(a.m_upper); } + bool is_ge_0(interval const& a) const { + return !lower_is_inf(a) && !lower_is_open(a) && unsynch_mpq_manager::is_zero(a.m_lower); + } + bool is_le_0(interval const &a) const { + return !upper_is_inf(a) && !upper_is_open(a) && unsynch_mpq_manager::is_zero(a.m_upper); + } // Setters void set_lower(interval& a, mpq const& n) const { m_manager.set(a.m_lower, n); } @@ -191,6 +197,12 @@ public: bool lower_is_inf(const interval& a) const { return m_config.lower_is_inf(a); } bool lower_is_open(const interval& a) const { return m_config.lower_is_open(a); } bool upper_is_open(const interval& a) const { return m_config.upper_is_open(a); } + bool is_ge_0(const interval &a) const { return m_config.is_ge_0(a); } + bool is_le_0(const interval &a) const { + return m_config.is_le_0(a); + } + + template void get_upper_dep(const interval& a, T& expl) { linearize(a.m_upper_dep, expl); } template diff --git a/src/math/lp/horner.cpp b/src/math/lp/horner.cpp index fffb4357c7..ec6412636d 100644 --- a/src/math/lp/horner.cpp +++ b/src/math/lp/horner.cpp @@ -104,6 +104,12 @@ bool horner::horner_lemmas() { TRACE(nla_solver, tout << "not generating horner lemmas\n";); return false; } + // Tighten the bounds of nonlinear variables over the LP tableau before the + // cross-nested/horner interval evaluation, so the tighter implied bounds can + // exclude zero and expose conflicts. Done here (instead of core::propagate) + // so the LP maximization only runs when horner is actually scheduled. + // optimize_nl_bounds() checks arith.nl.optimize_bounds internally. + c().optimize_nl_bounds(); c().lp_settings().stats().m_horner_calls++; const auto& matrix = c().lra.A_r(); // choose only rows that depend on m_to_refine variables diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 2f7dc11d25..dcf5e2ef49 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -644,7 +644,6 @@ namespace lp { for (auto c : cs) m_imp->m_constraints.display(tout, c) << "\n"; }); - SASSERT(bound_dep != nullptr); dep = dep_manager().mk_join(dep, bound_dep); } return dep; diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp index ff31f7ef53..2288cfbee7 100644 --- a/src/math/lp/monomial_bounds.cpp +++ b/src/math/lp/monomial_bounds.cpp @@ -14,247 +14,61 @@ namespace nla { - monomial_bounds::monomial_bounds(core* c): - common(c), - dep(c->m_intervals.get_dep_intervals()) { - - std::function fixed_eh = [c, this](lpvar v) { - c->trail().push(push_back_vector(m_fixed_var_trail)); - m_fixed_var_trail.push_back(v); - }; -// uncomment to enable: -// c->lra.m_fixed_var_eh = fixed_eh; - } + monomial_bounds::monomial_bounds(core *c) : common(c), dep(c->m_intervals.get_dep_intervals()) {} - void monomial_bounds::propagate() { - for (lpvar v : c().m_to_refine) { - propagate(c().emon(v)); - if (add_lemma()) + void monomial_bounds::generate_lemmas() { + for (auto v : c().m_to_refine) { + generate_lemma(c().emon(v)); + if (add_lemma()) break; } - propagate_fixed_vars(); } - void monomial_bounds::propagate_fixed_vars() { - if (m_fixed_var_qhead == m_fixed_var_trail.size()) - return; - c().trail().push(value_trail(m_fixed_var_qhead)); - while (m_fixed_var_qhead < m_fixed_var_trail.size()) - propagate_fixed_var(m_fixed_var_trail[m_fixed_var_qhead++]); - } - - void monomial_bounds::propagate_fixed_var(lpvar v) { - SASSERT(c().var_is_fixed(v)); - TRACE(nla_solver, tout << "propagate fixed var: " << c().var_str(v) << "\n";); - for (auto const& m : c().emons().get_use_list(v)) - propagate_fixed_var(m, v); - } - - void monomial_bounds::propagate_fixed_var(monic const& m, lpvar v) { - unsigned num_free = 0; - lpvar free_var = null_lpvar; - for (auto w : m) - if (!c().var_is_fixed(w)) - ++num_free, free_var = w; - if (num_free != 1) - return; - u_dependency* d = nullptr; - auto& lra = c().lra; - lp::mpq coeff(1); - for (auto w : m) { - if (c().var_is_fixed(w)) { - d = lra.join_deps(d, lra.get_bound_constraint_witnesses_for_column(w)); - coeff *= lra.get_lower_bound(w).x; - } - } - vector> coeffs; - coeffs.push_back({coeff, free_var}); - coeffs.push_back({mpq(-1), m.var()}); - lpvar j = lra.add_term(coeffs, UINT_MAX); - lra.update_column_type_and_bound(j, llc::EQ, mpq(0), d); - } - - bool monomial_bounds::is_too_big(mpq const& q) const { + bool monomial_bounds::is_too_big(mpq const &q) const { return rational(q).bitsize() > 256; } /** * Accumulate product of variables in monomial starting at position 'start' */ - void monomial_bounds::compute_product(unsigned start, monic const& m, scoped_dep_interval& product) { + void monomial_bounds::compute_product(unsigned start, monic const &m, scoped_dep_interval &product) { scoped_dep_interval vi(dep); unsigned power = 1; - for (unsigned i = start; i < m.size(); ) { + for (unsigned i = start; i < m.size();) { lpvar v = m.vars()[i]; var2interval(v, vi); ++i; - for (power = 1; i < m.size() && m.vars()[i] == v; ++i, ++power); - dep.power(vi, power, vi); + for (power = 1; i < m.size() && m.vars()[i] == v; ++i, ++power) + ; + dep.power(vi, power, vi); dep.mul(product, vi, product); } } - /** - * Monomial definition implies that a variable v is within 'range' - * If the current value of v is outside of the range, we add - * a bounds axiom. - */ - bool monomial_bounds::propagate_value(dep_interval& range, lpvar v) { - - bool propagated = false; - if (should_propagate_upper(range, v, 1)) { - auto const& upper = dep.upper(range); - auto cmp = dep.upper_is_open(range) ? llc::LT : llc::LE; - ++c().lra.settings().stats().m_nla_propagate_bounds; - lp::explanation ex; - dep.get_upper_dep(range, ex); - if (is_too_big(upper)) - return false; - lemma_builder lemma(c(), "propagate value - upper bound of range is below value"); - lemma &= ex; - lemma |= ineq(v, cmp, upper); - TRACE(nla_solver, dep.display(tout << c().val(v) << " > ", range) << "\n" << lemma << "\n";); - propagated = true; - } - if (should_propagate_lower(range, v, 1)) { - auto const& lower = dep.lower(range); - auto cmp = dep.lower_is_open(range) ? llc::GT : llc::GE; - ++c().lra.settings().stats().m_nla_propagate_bounds; - lp::explanation ex; - dep.get_lower_dep(range, ex); - if (is_too_big(lower)) - return false; - lemma_builder lemma(c(), "propagate value - lower bound of range is above value"); - lemma &= ex; - lemma |= ineq(v, cmp, lower); - TRACE(nla_solver, dep.display(tout << c().val(v) << " < ", range) << "\n" << lemma << "\n";); - propagated = true; - } - return propagated; - } - bool monomial_bounds::should_propagate_lower(dep_interval const& range, lpvar v, unsigned p) { + + bool monomial_bounds::should_propagate_lower(dep_interval const &range, lpvar v, unsigned p) { if (dep.lower_is_inf(range)) return false; auto bound = c().val(v); - auto const& lower = dep.lower(range); + auto const &lower = dep.lower(range); if (p > 1) bound = power(bound, p); return bound < lower; } - bool monomial_bounds::should_propagate_upper(dep_interval const& range, lpvar v, unsigned p) { + bool monomial_bounds::should_propagate_upper(dep_interval const &range, lpvar v, unsigned p) { if (dep.upper_is_inf(range)) return false; auto bound = c().val(v); - auto const& upper = dep.upper(range); + auto const &upper = dep.upper(range); if (p > 1) bound = power(bound, p); return bound > upper; } - /** - * Ensure that bounds are integral when the variable is integer. - */ - void monomial_bounds::propagate_bound(lpvar v, lp::lconstraint_kind cmp, rational const& q, u_dependency* d) { - SASSERT(cmp != llc::EQ && cmp != llc::NE); - if (!c().var_is_int(v)) - c().lra.update_column_type_and_bound(v, cmp, q, d); - else if (q.is_int()) { - if (cmp == llc::GT) - c().lra.update_column_type_and_bound(v, llc::GE, q + 1, d); - else if(cmp == llc::LT) - c().lra.update_column_type_and_bound(v, llc::LE, q - 1, d); - else - c().lra.update_column_type_and_bound(v, cmp, q, d); - } - else if (cmp == llc::GE || cmp == llc::GT) - c().lra.update_column_type_and_bound(v, llc::GE, ceil(q), d); - else - c().lra.update_column_type_and_bound(v, llc::LE, floor(q), d); - } - - - /** - * val(v)^p should be in range. - * if val(v)^p > upper(range) add - * v <= root(p, upper(range)) and v >= -root(p, upper(range)) if p is even - * v <= root(p, upper(range)) if p is odd - * if val(v)^p < lower(range) add - * v >= root(p, lower(range)) or v <= -root(p, lower(range)) if p is even - * v >= root(p, lower(range)) if p is odd - */ - - bool monomial_bounds::propagate_value(dep_interval& range, lpvar v, unsigned p) { - SASSERT(p > 0); - if (p == 1) - return propagate_value(range, v); - rational r; - if (should_propagate_upper(range, v, p)) { // v.upper^p > range.upper - lp::explanation ex; - dep.get_upper_dep(range, ex); - // p even, range.upper < 0, v^p >= 0 -> infeasible - if (p % 2 == 0 && rational(dep.upper(range)).is_neg()) { - ++c().lra.settings().stats().m_nla_propagate_bounds; - lemma_builder lemma(c(), "range requires a non-negative upper bound"); - lemma &= ex; - return true; - } - - if (rational(dep.upper(range)).root(p, r)) { - // v = -2, [-4,-3]^3 < v^3 -> add bound v <= -3 - // v = -2, [-1,+1]^2 < v^2 -> add bound v >= -1 - - if ((p % 2 == 1) || c().val(v).is_pos()) { - ++c().lra.settings().stats().m_nla_propagate_bounds; - auto le = dep.upper_is_open(range) ? llc::LT : llc::LE; - lemma_builder lemma(c(), "propagate value - root case - upper bound of range is below value"); - lemma &= ex; - lemma |= ineq(v, le, r); - return true; - } - - if (p % 2 == 0 && c().val(v).is_neg()) { - ++c().lra.settings().stats().m_nla_propagate_bounds; - SASSERT(!r.is_neg()); - auto ge = dep.upper_is_open(range) ? llc::GT : llc::GE; - lemma_builder lemma(c(), "propagate value - root case - upper bound of range is below negative value"); - lemma &= ex; - lemma |= ineq(v, ge, -r); - return true; - } - } - } - - if (should_propagate_lower(range, v, p)) { // v.lower^p < range.lower - // - // range.lower < 0 -> v.lower >= root(p, range.lower) - // range.lower >= 0, p odd -> v.lower >= root(p, range.lower) - // range.lower >= 0, p even, v.lower >= 0 -> v.lower >= root(p, range.lower) - // default: - // v.lower >= root(p, range.lower) || (p even & v.upper <= -root(p, range.lower)) - // - // pre-condition: p even -> range.lower >= 0 - // - if (rational(dep.lower(range)).root(p, r)) { - ++c().lra.settings().stats().m_nla_propagate_bounds; - auto ge = dep.lower_is_open(range) ? llc::GT : llc::GE; - auto le = dep.lower_is_open(range) ? llc::LT : llc::LE; - lp::explanation ex; - dep.get_lower_dep(range, ex); - lemma_builder lemma(c(), "propagate value - root case - lower bound of range is above value"); - lemma &= ex; - lemma |= ineq(v, ge, r); - if (p % 2 == 0) - lemma |= ineq(v, le, -r); - return true; - } - } - return false; - } - - void monomial_bounds::var2interval(lpvar v, scoped_dep_interval& i) { - u_dependency* d = nullptr; + void monomial_bounds::var2interval(lpvar v, scoped_dep_interval &i) { + u_dependency *d = nullptr; rational bound; bool is_strict; if (c().has_lower_bound(v, d, bound, is_strict)) { @@ -269,7 +83,7 @@ namespace nla { if (c().has_upper_bound(v, d, bound, is_strict)) { dep.set_upper_is_open(i, is_strict); dep.set_upper(i, bound); - dep.set_upper_dep(i, d); + dep.set_upper_dep(i, d); dep.set_upper_is_inf(i, false); } else { @@ -278,55 +92,72 @@ namespace nla { } /** - * Propagate bounds for monomial 'm'. - * For each variable v in m, compute the intervals of the remaining variables in m. - * Compute also the interval for m.var() as mi - * If the value of v is outside of mi / product_of_other, add a bounds lemma. - * If the value of m.var() is outside of product_of_all_vars, add a bounds lemma. + * Interval-based lemma generation for monomial 'm'. + * Runs the shared-factor (sandwich) and binomial-sign propagators. + * These emit lemmas; they do not tighten LP bounds. */ - bool monomial_bounds::propagate(monic const& m) { + bool monomial_bounds::generate_lemma(monic const &m) { unsigned num_free, power; lpvar free_var; analyze_monomial(m, num_free, free_var, power); - bool do_propagate_up = num_free == 0; + bool do_propagate_down = !is_free(m.var()) && num_free <= 1; + if (do_propagate_down && c().params().arith_nl_monomial_sandwich() && propagate_shared_factor(m)) + return true; + if (c().params().arith_nl_monomial_binomial_sign() && propagate_binomial_sign(m)) + return true; + return false; + } + + /** + * LP-bound tightening for monomial 'm'. + * For each variable v in m, divide the interval of m.var() by the product of + * the other variables and strengthen v's LP bounds (down-propagation). + * Finally strengthen the LP bounds of m.var() from the product interval. + * Unlike generate_lemma(), this emits no lemmas -- it only tightens LP bounds. + */ + bool monomial_bounds::tighten_lp(monic const &m) { + unsigned num_free, power; + lpvar free_var; + analyze_monomial(m, num_free, free_var, power); + bool do_propagate_up = num_free == 0; bool do_propagate_down = !is_free(m.var()) && num_free <= 1; if (!do_propagate_up && !do_propagate_down) return false; + scoped_dep_interval product(dep); scoped_dep_interval vi(dep), mi(dep); scoped_dep_interval other_product(dep); var2interval(m.var(), mi); dep.set_value(product, rational::one()); - for (unsigned i = 0; i < m.size(); ) { + bool tightened = false; + for (unsigned i = 0; i < m.size();) { lpvar v = m.vars()[i]; ++i; - for (power = 1; i < m.size() && v == m.vars()[i]; ++i, ++power); + for (power = 1; i < m.size() && v == m.vars()[i]; ++i, ++power) + ; var2interval(v, vi); - dep.power(vi, power, vi); + dep.power(vi, power, vi); if (do_propagate_down && (num_free == 0 || free_var == v)) { dep.set(other_product, product); compute_product(i, m, other_product); - if (propagate_down(m, mi, v, power, other_product)) - return true; + if (tighten_lp_bound(mi, v, power, other_product)) + tightened = true; } dep.mul(product, vi, product); } - if (do_propagate_down && c().params().arith_nl_monomial_sandwich() && - propagate_shared_factor(m)) - return true; - if (c().params().arith_nl_monomial_binomial_sign() && - propagate_binomial_sign(m)) - return true; - return do_propagate_up && propagate_value(product, m.var()); + if (!do_propagate_up) + return tightened; + return tighten_lp_bound(product, m.var(), 1) || tightened; } - bool monomial_bounds::propagate_down(monic const& m, dep_interval& mi, lpvar v, unsigned power, dep_interval& product) { - if (!dep.separated_from_zero(product)) + bool monomial_bounds::tighten_lp_bound(dep_interval &mi, lpvar v, unsigned power, + dep_interval &product) { + if (!dep.separated_from_zero(product)) return false; scoped_dep_interval range(dep); dep.div(mi, product, range); - return propagate_value(range, v, power); + return tighten_lp_bound(range, v, power); } bool monomial_bounds::is_free(lpvar v) const { @@ -368,17 +199,18 @@ namespace nla { } } - void monomial_bounds::unit_propagate() { + bool monomial_bounds::propagate_linear_monomials() { + bool propagated = false; for (lpvar v : c().m_monics_with_changed_bounds) { if (!c().is_monic_var(v)) continue; monic& m = c().emon(v); - unit_propagate(m); - if (add_lemma()) - break; - if (c().m_conflicts > 0) + if (propagate_linear_monomial(m)) + propagated = true; + if (c().lra.get_status() == lp::lp_status::INFEASIBLE) break; } + return propagated; } bool monomial_bounds::add_lemma() { @@ -391,27 +223,30 @@ namespace nla { return true; } - void monomial_bounds::unit_propagate(monic & m) { + bool monomial_bounds::propagate_linear_monomial(monic & m) { if (m.is_propagated()) - return; + return false; lpvar w, fixed_to_zero; if (!is_linear(m, w, fixed_to_zero)) - return; + return false; c().emons().set_propagated(m); + bool propagated = false; if (fixed_to_zero != null_lpvar) { - propagate_fixed_to_zero(m, fixed_to_zero); + propagated = propagate_fixed_to_zero(m, fixed_to_zero); } else { rational k = fixed_var_product(m, w); if (w == null_lpvar) - propagate_fixed(m, k); + propagated = propagate_fixed(m, k); else - propagate_nonfixed(m, k, w); + propagated = propagate_nonfixed(m, k, w); } - ++c().lra.settings().stats().m_nla_propagate_eq; + if (propagated) + ++c().lra.settings().stats().m_nla_propagate_eq; + return propagated; } lp::explanation monomial_bounds::get_explanation(u_dependency* dep) { @@ -423,25 +258,31 @@ namespace nla { return exp; } - void monomial_bounds::propagate_fixed_to_zero(monic const& m, lpvar fixed_to_zero) { + bool monomial_bounds::propagate_fixed_to_zero(monic const& m, lpvar fixed_to_zero) { + if (c().var_is_fixed_to_zero(m.var())) + return false; auto* dep = c().lra.get_bound_constraint_witnesses_for_column(fixed_to_zero); TRACE(nla_solver, tout << "propagate fixed " << m << " = 0, fixed_to_zero = " << fixed_to_zero << "\n";); c().lra.update_column_type_and_bound(m.var(), lp::lconstraint_kind::EQ, rational(0), dep); // propagate fixed equality c().add_fixed_equality(m.var(), rational(0), get_explanation(dep)); + return true; } - void monomial_bounds::propagate_fixed(monic const& m, rational const& k) { + bool monomial_bounds::propagate_fixed(monic const& m, rational const& k) { + if (c().var_is_fixed(m.var()) && c().get_lower_bound(m.var()) == k) + return false; auto* dep = explain_fixed(m, k); TRACE(nla_solver, tout << "propagate fixed " << m << " = " << k << "\n";); c().lra.update_column_type_and_bound(m.var(), lp::lconstraint_kind::EQ, k, dep); // propagate fixed equality c().add_fixed_equality(m.var(), k, get_explanation(dep)); + return true; } - void monomial_bounds::propagate_nonfixed(monic const& m, rational const& k, lpvar w) { + bool monomial_bounds::propagate_nonfixed(monic const& m, rational const& k, lpvar w) { vector> coeffs; coeffs.push_back({-k, w}); coeffs.push_back({rational::one(), m.var()}); @@ -453,6 +294,7 @@ namespace nla { if (k == 1) { c().add_equality(m.var(), w, get_explanation(dep)); } + return true; } u_dependency* monomial_bounds::explain_fixed(monic const& m, rational const& k) { @@ -506,13 +348,6 @@ namespace nla { return r; } - lpvar monomial_bounds::non_fixed_var(monic const& m) { - for (lpvar v : m) - if (!c().var_is_fixed(v)) - return v; - return null_lpvar; - } - /** * Dual-row shared-factor sandwich. For a binary monomial m = u*v, find LP * term columns whose term has shape a_m * m + a_v * v (exactly two @@ -592,7 +427,7 @@ namespace nla { << " m=" << m.var() << " v=" << v << " u=" << u << " a_m=" << a_m << " a_v=" << a_v << "\n";); - if (propagate_value(u_int, u)) + if (tighten_lp_bound(u_int, u, 1)) return true; // one lemma per call to keep the channel quiet } return false; @@ -698,5 +533,149 @@ namespace nla { return try_anchor(f1, f0) || try_anchor(f0, f1); } + /** + * range is an interval that v^p is guaranteed to lie in. + * Strengthen the *upper* bound of v from range, analogously to the upper + * branch of propagate_value(range, v, p), but only when a single bound on v + * follows (no lemmas). We use the existing bounds of v -- not its value -- + * to resolve the sign for even powers. + * + * An upper bound on v is implied by: + * range.upper = U: + * p odd -> v <= root(p, U) + * p even, U >= 0 -> v <= root(p, U) (|v| <= root(p, U)) + * range.lower = L, p even, v known <= 0: + * v <= -root(p, L) (resolves the disjunction) + * Only exact rational roots are used, so bounds that are not obtained from + * propagation are out of scope. + */ + bool monomial_bounds::tighten_lp_upper_bound(dep_interval const &range, lpvar v, unsigned p) { + SASSERT(p > 0); + auto improves_upper = [&](rational const& cand) { + return !c().has_upper_bound(v) || cand < c().get_upper_bound(v); + }; + bool tightened = false; + rational r; + // From range.upper: v <= root(p, U). + if (!dep.upper_is_inf(range)) { + rational U(dep.upper(range)); + if (U.root(p, r) && improves_upper(r)) { + auto cmp = dep.upper_is_open(range) ? llc::LT : llc::LE; + propagate_lp_bound(v, cmp, r, dep.get_upper_dep(range)); + tightened = true; + } + } + // Even power, v known non-positive: range.lower gives v <= -root(p, L). + if ((p & 1) == 0 && !dep.lower_is_inf(range) && + c().has_upper_bound(v) && !c().get_upper_bound(v).is_pos()) { + rational L(dep.lower(range)); + if (!L.is_neg() && L.root(p, r) && improves_upper(-r)) { + auto cmp = dep.lower_is_open(range) ? llc::LT : llc::LE; + u_dependency* d = c().lra.join_deps(dep.get_lower_dep(range), + c().lra.get_column_upper_bound_witness(v)); + propagate_lp_bound(v, cmp, -r, d); + tightened = true; + } + } + return tightened; + } + + /** + * range is an interval that v^p is guaranteed to lie in. + * Strengthen the *lower* bound of v from range (mirror of the above). + * + * A lower bound on v is implied by: + * range.lower = L: + * p odd -> v >= root(p, L) + * range.upper = U, p even, U >= 0: + * v >= -root(p, U) (|v| <= root(p, U)) + * range.lower = L, p even, v known >= 0: + * v >= root(p, L) (resolves the disjunction) + */ + bool monomial_bounds::tighten_lp_lower_bound(dep_interval const &range, lpvar v, unsigned p) { + SASSERT(p > 0); + auto improves_lower = [&](rational const& cand) { + return !c().has_lower_bound(v) || cand > c().get_lower_bound(v); + }; + bool tightened = false; + rational r; + if ((p & 1) == 1) { + // From range.lower: v >= root(p, L). + if (!dep.lower_is_inf(range)) { + rational L(dep.lower(range)); + if (L.root(p, r) && improves_lower(r)) { + auto cmp = dep.lower_is_open(range) ? llc::GT : llc::GE; + propagate_lp_bound(v, cmp, r, dep.get_lower_dep(range)); + tightened = true; + } + } + return tightened; + } + // Even power. From range.upper: v >= -root(p, U). + if (!dep.upper_is_inf(range)) { + rational U(dep.upper(range)); + if (!U.is_neg() && U.root(p, r) && improves_lower(-r)) { + auto cmp = dep.upper_is_open(range) ? llc::GT : llc::GE; + propagate_lp_bound(v, cmp, -r, dep.get_upper_dep(range)); + tightened = true; + } + } + // Even power, v known non-negative: range.lower gives v >= root(p, L). + if (!dep.lower_is_inf(range) && + c().has_lower_bound(v) && !c().get_lower_bound(v).is_neg()) { + rational L(dep.lower(range)); + if (!L.is_neg() && L.root(p, r) && improves_lower(r)) { + auto cmp = dep.lower_is_open(range) ? llc::GT : llc::GE; + u_dependency* d = c().lra.join_deps(dep.get_lower_dep(range), + c().lra.get_column_lower_bound_witness(v)); + propagate_lp_bound(v, cmp, r, d); + tightened = true; + } + } + return tightened; + } + + /** + * Ensure that bounds are integral when the variable is integer. + */ + void monomial_bounds::propagate_lp_bound(lpvar v, lp::lconstraint_kind cmp, rational const &q, u_dependency *d) { + SASSERT(cmp != llc::EQ && cmp != llc::NE); + IF_VERBOSE(1, verbose_stream() << "propagate_lp_bound: v=" << v << " cmp=" << cmp << " q=" << q << "\n";); + if (!c().var_is_int(v)) + c().lra.update_column_type_and_bound(v, cmp, q, d); + else if (q.is_int()) { + if (cmp == llc::GT) + c().lra.update_column_type_and_bound(v, llc::GE, q + 1, d); + else if (cmp == llc::LT) + c().lra.update_column_type_and_bound(v, llc::LE, q - 1, d); + else + c().lra.update_column_type_and_bound(v, cmp, q, d); + } + else if (cmp == llc::GE || cmp == llc::GT) + c().lra.update_column_type_and_bound(v, llc::GE, ceil(q), d); + else + c().lra.update_column_type_and_bound(v, llc::LE, floor(q), d); + } + + bool monomial_bounds::tighten_lp_bound(dep_interval const &range, lpvar v, unsigned power) { + bool propagated = false; + if (tighten_lp_upper_bound(range, v, power)) + propagated = true; + if (tighten_lp_lower_bound(range, v, power)) + propagated = true; + return propagated; + } + + bool monomial_bounds::tighten_lp_bounds() { + bool new_bound = false; + for (auto &m : c().emons()) + if (tighten_lp(m)) + new_bound = true; + if (propagate_linear_monomials()) + new_bound = true; + IF_VERBOSE(1, verbose_stream() << "tighten_lp_bounds: new_bound=" << new_bound << "\n";); + return new_bound; + } + } diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h index 564fda6987..a018af4539 100644 --- a/src/math/lp/monomial_bounds.h +++ b/src/math/lp/monomial_bounds.h @@ -17,22 +17,27 @@ namespace nla { class monomial_bounds : common { dep_intervals& dep; + bool tighten_lp_bound(dep_interval const &range, lpvar v, unsigned p); + bool tighten_lp_upper_bound(dep_interval const& range, lpvar v, unsigned p); + bool tighten_lp_lower_bound(dep_interval const& range, lpvar v, unsigned p); + bool tighten_lp_bound(dep_interval &mi, lpvar v, unsigned power, dep_interval &product); + + void propagate_lp_bound(lpvar v, lp::lconstraint_kind cmp, rational const &q, u_dependency *d); + + bool should_propagate_lower(dep_interval const& range, lpvar v, unsigned p); bool should_propagate_upper(dep_interval const& range, lpvar v, unsigned p); - void propagate_bound(lpvar v, lp::lconstraint_kind cmp, rational const& q, u_dependency* d); + void var2interval(lpvar v, scoped_dep_interval& i); bool is_too_big(mpq const& q) const; - bool propagate_down(monic const& m, lpvar u); - bool propagate_value(dep_interval& range, lpvar v); - bool propagate_value(dep_interval& range, lpvar v, unsigned power); void compute_product(unsigned start, monic const& m, scoped_dep_interval& i); - bool propagate(monic const& m); - void propagate_fixed_to_zero(monic const& m, lpvar fixed_to_zero); - void propagate_fixed(monic const& m, rational const& k); - void propagate_nonfixed(monic const& m, rational const& k, lpvar w); + bool generate_lemma(monic const& m); + bool tighten_lp(monic const& m); + bool propagate_fixed_to_zero(monic const& m, lpvar fixed_to_zero); + bool propagate_fixed(monic const& m, rational const& k); + bool propagate_nonfixed(monic const& m, rational const& k, lpvar w); u_dependency* explain_fixed(monic const& m, rational const& k); lp::explanation get_explanation(u_dependency* dep); - bool propagate_down(monic const& m, dep_interval& mi, lpvar v, unsigned power, dep_interval& product); bool propagate_shared_factor(monic const& m); bool propagate_binomial_sign(monic const& m); void analyze_monomial(monic const& m, unsigned& num_free, lpvar& free_v, unsigned& power) const; @@ -40,21 +45,16 @@ namespace nla { bool is_zero(lpvar v) const; bool add_lemma(); - // monomial propagation - void unit_propagate(monic & m); + // linear-monomial equality propagation: + // when all but one variable of a monomial are fixed, the monomial is + // linear and its value/equality can be propagated into the LP solver. + bool propagate_linear_monomial(monic & m); + bool propagate_linear_monomials(); bool is_linear(monic const& m, lpvar& w, lpvar & fixed_to_zero); rational fixed_var_product(monic const& m, lpvar w); - lpvar non_fixed_var(monic const& m); - - // fixed variable propagation - unsigned m_fixed_var_qhead = 0; - unsigned_vector m_fixed_var_trail; - void propagate_fixed_vars(); - void propagate_fixed_var(lpvar v); - void propagate_fixed_var(monic const& m, lpvar v); public: monomial_bounds(core* core); - void propagate(); - void unit_propagate(); + void generate_lemmas(); + bool tighten_lp_bounds(); }; } diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index a78dfc451c..2cfe1d2efa 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -1308,7 +1308,7 @@ lbool core::check(unsigned level) { auto no_effect = [&]() { return ret == l_undef && !done() && m_lemmas.empty() && m_literals.empty() && !m_check_feasible; }; if (no_effect()) - m_monomial_bounds.propagate(); + m_monomial_bounds.generate_lemmas(); if (no_effect() && refine_pseudo_linear()) return l_false; @@ -1522,19 +1522,94 @@ void core::set_use_nra_model(bool m) { m_use_nra_model = m; } } + -void core::propagate() { -#if Z3DEBUG - flet f(lra.validate_blocker(), true); -#endif - clear(); - m_monomial_bounds.unit_propagate(); +bool core::propagate() { + clear(); + bool propagated = m_monomial_bounds.tighten_lp_bounds(); m_monics_with_changed_bounds.reset(); + return propagated; } +/** + \brief Tighten the bounds of variables occurring in nonlinear monomials by + maximizing/minimizing them over the LP tableau (analogous to theory_arith's + max_min_nl_vars). The tighter implied bounds, each carrying an LP explanation, + let the subsequent horner/cross-nested interval evaluation exclude zero and + detect a conflict that would otherwise be missed with only the propagated + bounds. +*/ +bool core::optimize_nl_bounds() { + if (!params().arith_nl_optimize_bounds() || !m_bounds_optimization_enabled) + return false; + + trail().push(value_trail(m_bounds_optimization_enabled)); + m_bounds_optimization_enabled = false; + + if (!lra.is_feasible()) + return false; + if (lra.find_feasible_solution() == lp::lp_status::INFEASIBLE) + return false; + + // Gather the candidate columns: every non-fixed leaf variable that + // participates in a monomial (mirrors solver=2's max_min_nl_vars). + svector cands; + auto add = [&](lpvar j) { + if (active_var_set_contains(j)) + return; + insert_to_active_var_set(j); + if (lra.column_is_fixed(j)) + return; + cands.push_back(j); + }; + clear_active_var_set(); + for (auto const& m : m_emons) { + add(m.var()); + for (lpvar k : m.vars()) + add(k); + } + + // Throttle: the LP maximize/minimize cost scales with the number of + // candidate variables (two LP optimizations each). On large nonlinear + // problems this pass is expensive and rarely productive, so skip it when the + // candidate set exceeds the threshold (0 = unlimited). + unsigned const max_vars = params().arith_nl_optimize_bounds_lp_max_vars(); + if (max_vars != 0 && cands.size() > max_vars) + return false; + + // Collect improved bounds first (each find_improved_bound maximizes a term + // over the *unchanged* constraint set, so all improvements are valid implied + // bounds), then apply them together and re-establish feasibility once. + // Interleaving update_column_type_and_bound between the maximize calls + // corrupts the core solver's x/inf_heap (maximize_term_on_tableau issues a + // raw solve() that does not reconcile pending bound changes). + struct improved_bound { lpvar j; lp::lconstraint_kind kind; rational bound; u_dependency* dep; }; + vector improvements; + for (lpvar j : cands) { + if (!lra.is_feasible()) + break; + for (bool is_lower : { true, false }) { + rational bound; + u_dependency* dep = lra.find_improved_bound(j, is_lower, bound); + if (!dep) + continue; + auto kind = is_lower ? lp::lconstraint_kind::GE : lp::lconstraint_kind::LE; + improvements.push_back({ j, kind, bound, dep }); + } + } + + if (improvements.empty()) + return false; + + for (auto const& ib : improvements) + lra.update_column_type_and_bound(ib.j, ib.kind, ib.bound, ib.dep); + lra.find_feasible_solution(); + return true; +} + + void core::simplify() { // in-processing simplifiation can go here, such as bounds improvements. - } bool core::is_pseudo_linear(monic const& m) const { diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h index c4773ffa78..3d6eaafabb 100644 --- a/src/math/lp/nla_core.h +++ b/src/math/lp/nla_core.h @@ -108,6 +108,7 @@ class core { nla_throttle m_throttle; bool m_throttle_enabled = true; + bool m_bounds_optimization_enabled = true; @@ -118,6 +119,8 @@ class core { bool is_pseudo_linear(monic const& m) const; void refine_pseudo_linear(monic const& m); + bool optimize_nl_bounds(); + std::ostream& display_constraint_smt(std::ostream& out, unsigned id, lp::lar_base_constraint const& c) const; std::ostream& display_declarations_smt(std::ostream& out) const; @@ -406,7 +409,7 @@ public: bool no_lemmas_hold() const; - void propagate(); + bool propagate(); void simplify(); diff --git a/src/math/lp/nla_solver.cpp b/src/math/lp/nla_solver.cpp index da1e1b3a95..08f8c23433 100644 --- a/src/math/lp/nla_solver.cpp +++ b/src/math/lp/nla_solver.cpp @@ -54,8 +54,8 @@ namespace nla { return m_core->check(level); } - void solver::propagate() { - m_core->propagate(); + bool solver::propagate() { + return m_core->propagate(); } void solver::push(){ diff --git a/src/math/lp/nla_solver.h b/src/math/lp/nla_solver.h index 2e85bb18f4..79556c50cd 100644 --- a/src/math/lp/nla_solver.h +++ b/src/math/lp/nla_solver.h @@ -39,7 +39,7 @@ namespace nla { void pop(unsigned scopes); bool need_check(); lbool check(unsigned level); - void propagate(); + bool propagate(); void simplify() { m_core->simplify(); } lbool check_power(lpvar r, lpvar x, lpvar y); bool is_monic_var(lpvar) const; diff --git a/src/params/smt_params_helper.pyg b/src/params/smt_params_helper.pyg index 67c19ae9f1..5dffec3a4a 100644 --- a/src/params/smt_params_helper.pyg +++ b/src/params/smt_params_helper.pyg @@ -100,6 +100,7 @@ def_module_params(module_name='smt', ('arith.nl.delay', UINT, 10, 'number of calls to final check before invoking bounded nlsat check'), ('arith.nl.propagate_linear_monomials', BOOL, True, 'propagate linear monomials'), ('arith.nl.optimize_bounds', BOOL, True, 'enable bounds optimization'), + ('arith.nl.optimize_bounds_lp_max_vars', UINT, 120, 'skip LP-based nonlinear bounds optimization when the number of candidate monomial variables exceeds this threshold (0 = unlimited)'), ('arith.nl.cross_nested', BOOL, True, 'enable cross-nested consistency checking'), ('arith.nl.log', BOOL, False, 'Log lemmas sent to nra solver'), ('arith.propagate_eqs', BOOL, True, 'propagate (cheap) equalities'), diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 0001428726..0248fe2dd0 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -1349,7 +1349,6 @@ public: literal eqz = mk_literal(m.mk_eq(q, zero)); literal mod_ge_0 = mk_literal(a.mk_ge(mod, zero)); - // q = 0 or p = (p mod q) + q * (p div q) // q = 0 or (p mod q) >= 0 // q >= 0 or (p mod q) + q <= -1 @@ -1746,6 +1745,7 @@ public: IF_VERBOSE(12, verbose_stream() << "final-check " << lp().get_status() << "\n"); lbool is_sat = l_true; SASSERT(lp().ax_is_correct()); + propagate_nla(); if (!lp().is_feasible() || lp().has_changed_columns()) is_sat = make_feasible(); final_check_status st = FC_DONE; @@ -2126,7 +2126,6 @@ public: default: UNREACHABLE(); } - TRACE(arith, tout << "is_lower: " << is_lower << " pos " << pos << "\n";); expr_ref atom(m); // TBD utility: lp::lar_term term = mk_term(ineq.m_poly); // then term is used instead of ineq.m_term @@ -2135,6 +2134,7 @@ public: else // create term >= 0 (or term <= 0) atom = mk_bound(ineq.term(), ineq.rs(), is_lower); + TRACE(arith, tout << "is_lower: " << is_lower << " pos " << pos << " " << atom << "\n";); return literal(ctx().get_bool_var(atom), pos); } @@ -2265,7 +2265,6 @@ public: bool propagate_core() { m_model_is_initialized = false; flush_bound_axioms(); - propagate_nla(); if (ctx().inconsistent()) return true; if (!can_propagate_core()) @@ -2329,12 +2328,14 @@ public: return true; } - void propagate_nla() { + bool propagate_nla() { + bool propagated = false; if (m_nla) { - m_nla->propagate(); + propagated = m_nla->propagate() || propagated; add_lemmas(); lp().collect_more_rows_for_lp_propagation(); } + return propagated; } void add_equality(lpvar j, rational const& k, lp::explanation const& exp) { diff --git a/src/util/mpz.cpp b/src/util/mpz.cpp index 7fd91404cb..465b5247c9 100644 --- a/src/util/mpz.cpp +++ b/src/util/mpz.cpp @@ -2459,8 +2459,10 @@ static unsigned div_u(unsigned k, unsigned n) { template bool mpz_manager::root(mpz & a, unsigned n) { - SASSERT(n % 2 != 0 || is_nonneg(a)); - if (is_zero(a)) { + SASSERT(n != 0); + if (n % 2 == 0 && is_neg(a)) + return false; + if (is_zero(a) || n == 1) { return true; // precise } diff --git a/src/util/mpz.h b/src/util/mpz.h index 7b714b9bed..d00425349c 100644 --- a/src/util/mpz.h +++ b/src/util/mpz.h @@ -703,7 +703,7 @@ public: Otherwise return false, and update a with the smallest integer r such that r*r > n. - \remark This method assumes that if n is even, then a is nonegative + \remark This method returns false if a is negative and n is even. */ bool root(mpz & a, unsigned n); bool root(mpz const & a, unsigned n, mpz & r) { set(r, a); return root(r, n); }