From 0692c3e01daa0595cbddbc8579c9808bfc284d3b Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Wed, 29 Jul 2026 09:12:56 -0700 Subject: [PATCH] Fstar opt2 (#10261) Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> Copilot-Session: 84b07e32-6458-4ea9-bf14-1cecfb7f1a99 --- src/math/lp/horner.cpp | 7 + src/math/lp/monomial_bounds.cpp | 377 +++++++++++++++++++++++++++++++ src/math/lp/monomial_bounds.h | 38 ++++ src/math/lp/nla_core.cpp | 24 +- src/math/lp/nla_core.h | 6 + src/math/lp/nla_grobner.cpp | 154 +++++++++++++ src/math/lp/nla_grobner.h | 6 + src/params/smt_params_helper.pyg | 1 + 8 files changed, 609 insertions(+), 4 deletions(-) diff --git a/src/math/lp/horner.cpp b/src/math/lp/horner.cpp index ec6412636d..cc30537eff 100644 --- a/src/math/lp/horner.cpp +++ b/src/math/lp/horner.cpp @@ -110,6 +110,13 @@ bool horner::horner_lemmas() { // so the LP maximization only runs when horner is actually scheduled. // optimize_nl_bounds() checks arith.nl.optimize_bounds internally. c().optimize_nl_bounds(); + // optimize_nl_bounds re-calibrated m_to_refine against the model it produced. + // If nothing remains to refine, every monomial is consistent under a feasible + // LP model: the nonlinear goal is satisfied. Declare it and stop. + if (c().to_refine().empty()) { + c().set_nla_satisfied(); + return false; + } 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/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp index 36f498e333..fa0c3251a4 100644 --- a/src/math/lp/monomial_bounds.cpp +++ b/src/math/lp/monomial_bounds.cpp @@ -677,5 +677,382 @@ namespace nla { return new_bound; } + // ================================================================ + // max_min: incremental LP bound optimization. + // + // A direct adaptation of smt::theory_arith::max_min (see + // src/smt/theory_arith_aux.h). We maximize (or minimize) a single + // column 'v' over the current LP tableau by a bounded-effort primal + // simplex walk: repeatedly pick a non-basic variable that improves the + // objective, ratio-test its column to find the tightest blocking basic + // variable, and pivot. The tableau is left at a feasible vertex; the + // implied bound is then read off 'v's tableau row and rounded to respect + // the integrality of integer columns. + // + // Integrality is maintained during the walk (the 'maintain_integrality == + // true' configuration of theory_arith): every move of a column is a multiple + // of the integrality quantum 'min_gain', so integer columns keep integral + // values throughout. The final implied bound is additionally floored/ceiled + // for integer 'v'. + // ================================================================ + + static lp::impq mm_abs(lp::impq const& v) { + return v.is_neg() ? -v : v; + } + + // Round 'val' down to the nearest multiple of the (integral) 'divisor'. + // Mirrors theory_arith::normalize_gain. 'divisor == -1' means "no quantum". + static void mm_round_down(lp::impq& val, rational const& divisor) { + if (divisor.is_one()) + val = lp::impq(lp::floor(val)); + else if (!divisor.is_minus_one()) + val = lp::impq(lp::floor(val / divisor) * divisor); + } + + lpvar monomial_bounds::mm_basic_in_row(unsigned row) const { + return c().lra.get_base_column_in_row(row); + } + + // A gain is safe when the column is unbounded in the chosen direction, or + // the required integral quantum still fits within the maximal feasible move. + // Mirrors theory_arith::safe_gain. + bool monomial_bounds::mm_safe_gain(mm_gain const& g) const { + return g.unbounded || lp::impq(g.min_gain) <= g.max_gain; + } + + // Initialize the gain for moving 'x' in direction 'inc' (increase when inc, + // decrease otherwise) from its own bound. For integer columns the quantum + // 'min_gain' starts at 1. Mirrors theory_arith::init_gains. + monomial_bounds::mm_gain monomial_bounds::mm_init_gains(lpvar x, bool inc) const { + auto& s = c().lra; + mm_gain g; + if (inc && s.column_has_upper_bound(x)) { + g.unbounded = false; + g.max_gain = s.column_upper_bound(x) - s.get_column_value(x); + } + else if (!inc && s.column_has_lower_bound(x)) { + g.unbounded = false; + g.max_gain = s.get_column_value(x) - s.column_lower_bound(x); + } + if (s.column_is_int(x)) + g.min_gain = rational::one(); + return g; + } + + // Tighten 'g' by the room that basic variable 'x_i' (with coefficient 'a_ij' + // on the moving column) has before hitting a bound. When 'x_i' is an integer + // column, the quantum 'min_gain' is raised to the lcm of the denominators of + // the involved coefficients and both gains are rounded down to that quantum, + // so the move keeps 'x_i' integral. Returns true when 'max_gain' was + // strengthened. Mirrors theory_arith::update_gains. + bool monomial_bounds::mm_update_gains(bool inc, lpvar x_i, rational const& a_ij, mm_gain& g) const { + auto& s = c().lra; + SASSERT(!a_ij.is_zero()); + if (!mm_safe_gain(g)) + return false; + + bool decrement_x_i = (inc && a_ij.is_pos()) || (!inc && a_ij.is_neg()); + bool bounded_i = false; + lp::impq max_inc; + if (decrement_x_i && s.column_has_lower_bound(x_i)) { + max_inc = mm_abs((s.get_column_value(x_i) - s.column_lower_bound(x_i)) / a_ij); + bounded_i = true; + } + else if (!decrement_x_i && s.column_has_upper_bound(x_i)) { + max_inc = mm_abs((s.column_upper_bound(x_i) - s.get_column_value(x_i)) / a_ij); + bounded_i = true; + } + + bool xi_int = s.column_is_int(x_i); + rational den_aij(1); + if (xi_int) + den_aij = denominator(a_ij); + SASSERT(den_aij.is_pos() && den_aij.is_int()); + + // Moving 'x_i' by k requires moving the entering column by k/a_ij; to keep + // an integer 'x_i' integral the entering column must step in multiples of + // denominator(a_ij). Accumulate that into the quantum and re-round. + if (xi_int && !den_aij.is_one()) { + if (g.min_gain.is_neg()) + g.min_gain = den_aij; + else + g.min_gain = lcm(g.min_gain, den_aij); + if (!g.unbounded) + mm_round_down(g.max_gain, g.min_gain); + } + if (xi_int && !g.unbounded && !g.max_gain.is_int()) { + g.max_gain = lp::impq(lp::floor(g.max_gain)); + mm_round_down(g.max_gain, g.min_gain); + } + + if (bounded_i) { + if (xi_int) { + max_inc = lp::impq(lp::floor(max_inc)); + mm_round_down(max_inc, g.min_gain); + } + if (g.unbounded) { + g.unbounded = false; + g.max_gain = max_inc; + return true; + } + if (g.max_gain > max_inc) { + g.max_gain = max_inc; + return true; + } + } + return false; + } + + // Ratio test: for entering column 'x_j' moving in direction 'inc', find the + // basic variable 'x_i' that first blocks the move and the maximal gain. + // Returns false (unsafe) when the integrality quantum cannot be satisfied, so + // the caller treats 'x_j' as unusable. Mirrors theory_arith::pick_var_to_leave. + bool monomial_bounds::mm_pick_var_to_leave(lpvar x_j, bool inc, rational& a_ij, mm_gain& g, lpvar& x_i) const { + auto& s = c().lra; + x_i = null_lpvar; + g = mm_init_gains(x_j, inc); + // an integer entering column must sit at an integral value to move in + // integral steps. + if (s.column_is_int(x_j) && !s.get_column_value(x_j).is_int()) + return false; + for (auto const& cell : s.A_r().m_columns[x_j]) { + lpvar si = mm_basic_in_row(cell.var()); + rational const& coeff_ij = s.A_r().get_val(cell); + if (mm_update_gains(inc, si, coeff_ij, g) || + (x_i == null_lpvar && !g.unbounded)) { + x_i = si; + a_ij = coeff_ij; + } + } + return mm_safe_gain(g); + } + + // Apply 'delta' to non-basic column 'j', propagating to dependent basic + // columns (theory_arith::update_value). + void monomial_bounds::mm_update_value(lpvar j, lp::impq const& delta) { + if (delta.is_zero()) + return; + auto& s = c().lra; + lp::impq new_val = s.get_column_value(j) + delta; + s.set_value_for_nbasic_column_report(j, new_val, [](unsigned) {}); + } + + // Move (now non-basic) 'x_i' maximally towards its bound in direction 'inc' + // without violating other columns' bounds, in integral steps when 'x_i' is an + // integer column (theory_arith::move_to_bound). + bool monomial_bounds::mm_move_to_bound(lpvar x_i, bool inc, unsigned& best_efforts) { + auto& s = c().lra; + if (s.column_is_int(x_i) && !s.get_column_value(x_i).is_int()) { + ++best_efforts; + return false; + } + mm_gain g = mm_init_gains(x_i, inc); + for (auto const& cell : s.A_r().m_columns[x_i]) { + lpvar si = mm_basic_in_row(cell.var()); + rational const& coeff = s.A_r().get_val(cell); + mm_update_gains(inc, si, coeff, g); + } + bool result = false; + if (mm_safe_gain(g) && !g.unbounded) { + lp::impq step = g.max_gain; + if (!inc) + step = -step; + mm_update_value(x_i, step); + result = !g.max_gain.is_zero(); + } + if (!result) + ++best_efforts; + return result; + } + + // Primal-simplex walk maximizing/minimizing 'v' (theory_arith::max_min). + void monomial_bounds::mm_optimize(lpvar v, bool maximize) { + auto& s = c().lra; + unsigned best_efforts = 0; + unsigned const max_efforts = 20; + unsigned rounds = 0; + unsigned const max_rounds = 200; + + while (best_efforts < max_efforts && rounds < max_rounds && !c().lp_settings().get_cancel_flag()) { + ++rounds; + lpvar x_j = null_lpvar, x_i = null_lpvar; + rational a_ij(0); + mm_gain best; // gain of the selected move + bool inc = false; + bool has_bound = false; + + // Consider a candidate entering variable 'cand' whose coefficient in + // the objective (v expressed over the non-basic columns) is + // 'obj_coeff'. Returns true to stop scanning (unbounded direction). + auto consider = [&](lpvar cand, rational const& obj_coeff) -> bool { + bool curr_inc = obj_coeff.is_pos() ? maximize : !maximize; + if ((curr_inc && s.column_has_upper_bound(cand)) || + (!curr_inc && s.column_has_lower_bound(cand))) + has_bound = true; + // cannot move a variable already at the relevant bound + if (curr_inc && s.column_has_upper_bound(cand) && + s.get_column_value(cand) == s.column_upper_bound(cand)) + return false; + if (!curr_inc && s.column_has_lower_bound(cand) && + s.get_column_value(cand) == s.column_lower_bound(cand)) + return false; + rational curr_a(0); + mm_gain cur; + lpvar curr_xi = null_lpvar; + bool safe = mm_pick_var_to_leave(cand, curr_inc, curr_a, cur, curr_xi); + if (!safe) { + // the integrality quantum cannot be met on this column + has_bound = true; + ++best_efforts; + return false; + } + if (curr_xi == null_lpvar) { + // limited only by its own bound (or fully unbounded) + x_j = cand; x_i = null_lpvar; inc = curr_inc; best = cur; a_ij = curr_a; + return true; + } + if (cur.max_gain > best.max_gain) { + x_i = curr_xi; x_j = cand; a_ij = curr_a; best = cur; inc = curr_inc; + } + else if (cur.max_gain.is_zero() && (x_i == null_lpvar || curr_xi < x_i)) { + x_i = curr_xi; x_j = cand; a_ij = curr_a; best = cur; inc = curr_inc; + } + return false; + }; + + if (!s.is_base(v)) { + consider(v, rational::one()); + } + else { + unsigned ri = s.r_heading()[v]; + rational a_v(0); + for (auto const& e : s.A_r().m_rows[ri]) + if (e.var() == v) { a_v = e.coeff(); break; } + for (auto const& e : s.A_r().m_rows[ri]) { + if (e.var() == v) + continue; + // v = -(1/a_v) * sum a_e x_e, so d(v)/d(x_e) has the sign of + // -a_e/a_v; only the sign steers the search direction. + rational objc = -e.coeff(); + if (a_v.is_neg()) + objc.neg(); + if (consider(e.var(), objc)) + break; + } + } + + if (!has_bound && x_i == null_lpvar && x_j == null_lpvar) + return; // objective is unbounded in the chosen direction + if (x_j == null_lpvar) + return; // optimized: no improving move remains + + // a non-unit integral quantum means the exact optimum may not be + // reachable in integral steps: count it as best-effort progress. + if (best.min_gain.is_pos() && !best.min_gain.is_one()) + ++best_efforts; + + if (x_i == null_lpvar) { + // move x_j directly to its own bound + if (inc && s.column_has_upper_bound(x_j)) { + if (best.max_gain.is_zero()) + return; + mm_update_value(x_j, best.max_gain); + continue; + } + if (!inc && s.column_has_lower_bound(x_j)) { + if (best.max_gain.is_zero()) + return; + mm_update_value(x_j, -best.max_gain); + continue; + } + return; // unbounded + } + + // x_j can move exactly across to its opposite bound without pivoting + if (s.column_has_lower_bound(x_j) && s.column_has_upper_bound(x_j) && + s.column_lower_bound(x_j) != s.column_upper_bound(x_j) && + (s.column_upper_bound(x_j) - s.column_lower_bound(x_j) == best.max_gain)) { + lp::impq step = best.max_gain; + if (!inc) + step = -step; + mm_update_value(x_j, step); + continue; + } + + // pivot x_j into the basis (x_i leaves); the degenerate pivot keeps + // the current point, then move x_i to its bound to raise v. + s.pivot(x_j, x_i); + bool inc_xi = inc ? a_ij.is_neg() : a_ij.is_pos(); + mm_move_to_bound(x_i, inc_xi, best_efforts); + } + } + + // Read the implied bound on 'v' off its final tableau row and round it to + // respect the integrality of integer columns (theory_arith::mk_bound_from_row + // + normalize_bound). Returns the joined explanation, or nullptr if no bound + // is implied (e.g. a required bound on a row variable is missing). + u_dependency* monomial_bounds::mm_bound_from_row(lpvar v, bool maximize, rational& bound) { + auto& s = c().lra; + if (!s.is_base(v)) + return nullptr; + unsigned ri = s.r_heading()[v]; + auto const& row = s.A_r().m_rows[ri]; + rational a_v(0); + for (auto const& e : row) + if (e.var() == v) { a_v = e.coeff(); break; } + if (a_v.is_zero()) + return nullptr; + lp::impq acc(0); + u_dependency* dep = nullptr; + for (auto const& e : row) { + if (e.var() == v) + continue; + lpvar k = e.var(); + rational ck = -e.coeff() / a_v; // v = sum ck * x_k + if (ck.is_zero()) + continue; + bool use_upper = maximize ? ck.is_pos() : ck.is_neg(); + if (use_upper) { + if (!s.column_has_upper_bound(k)) + return nullptr; + acc += s.column_upper_bound(k) * ck; + dep = s.join_deps(dep, s.get_column_upper_bound_witness(k)); + } + else { + if (!s.column_has_lower_bound(k)) + return nullptr; + acc += s.column_lower_bound(k) * ck; + dep = s.join_deps(dep, s.get_column_lower_bound_witness(k)); + } + } + if (s.column_is_int(v)) + bound = maximize ? lp::floor(acc) : lp::ceil(acc); + else + bound = acc.x; + return dep; + } + + u_dependency* monomial_bounds::improve_bound(lpvar j, bool is_lower, rational& bound) { + auto& s = c().lra; + if (!s.is_feasible()) + return nullptr; + bool maximize = !is_lower; + mm_optimize(j, maximize); + rational b(0); + u_dependency* dep = mm_bound_from_row(j, maximize, b); + if (!dep) + return nullptr; + if (is_lower) { + if (s.column_has_lower_bound(j) && b <= s.column_lower_bound(j).x) + return nullptr; + } + else { + if (s.column_has_upper_bound(j) && b >= s.column_upper_bound(j).x) + return nullptr; + } + bound = b; + return dep; + } + } diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h index 37cdac8490..cc9738eb07 100644 --- a/src/math/lp/monomial_bounds.h +++ b/src/math/lp/monomial_bounds.h @@ -51,10 +51,48 @@ namespace nla { bool propagate_changed_bound(monic & m); bool is_linear(monic const& m, lpvar& w, lpvar & fixed_to_zero); rational fixed_var_product(monic const& m, lpvar w); + + // ---------------------------------------------------------------- + // max_min: incremental LP bound optimization. + // + // Adapted from smt::theory_arith::max_min (theory_arith_aux.h): a + // bounded-variable primal-simplex walk that maximizes/minimizes a + // single column over the current LP tableau by pivoting, leaving the + // tableau at an optimal feasible vertex. Replaces the expensive + // lar_solver::find_improved_bound (which rebuilds the objective and + // re-solves from scratch on every call) and rounds the resulting + // bound to respect the integrality of integer columns. + // ---------------------------------------------------------------- + // + // Gain of a candidate move: how far a non-basic column may travel in the + // chosen direction. 'min_gain' is the integrality quantum -- the + // smallest integral step that keeps every dependent integer column at an + // integral value (-1 means "no quantum", i.e. a rational column); + // 'max_gain' is the largest feasible step (valid only when !unbounded). + // + struct mm_gain { + bool unbounded = true; // no bound blocks the move + rational min_gain = rational::minus_one(); + lp::impq max_gain; + }; + lpvar mm_basic_in_row(unsigned row) const; + bool mm_safe_gain(mm_gain const& g) const; + mm_gain mm_init_gains(lpvar x, bool inc) const; + bool mm_update_gains(bool inc, lpvar x_i, rational const& a_ij, mm_gain& g) const; + bool mm_pick_var_to_leave(lpvar x_j, bool inc, rational& a_ij, mm_gain& g, lpvar& x_i) const; + bool mm_move_to_bound(lpvar x_i, bool inc, unsigned& best_efforts); + void mm_update_value(lpvar j, lp::impq const& delta); + void mm_optimize(lpvar v, bool maximize); + u_dependency* mm_bound_from_row(lpvar v, bool maximize, rational& bound); public: monomial_bounds(core* core); void generate_lemmas(); bool tighten_lp_bounds(); bool propagate_changed_bounds(); + + // Maximize (is_lower == false) or minimize (is_lower == true) column j + // over the LP tableau and, if the resulting bound improves j's current + // bound, return its explanation and set 'bound'; otherwise return nullptr. + u_dependency* improve_bound(lpvar j, bool is_lower, rational& bound); }; } diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index 46f2a7a1cf..8534cfdfce 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -1299,13 +1299,14 @@ lbool core::check(unsigned level) { if (m_to_refine.empty()) return l_true; init_search(); + m_nla_satisfied = false; lbool ret = l_undef; bool run_grobner = need_run_grobner(); bool run_horner = need_run_horner(); bool run_bounds = params().arith_nl_branching(); - auto no_effect = [&]() { return ret == l_undef && !done() && m_lemmas.empty() && m_literals.empty() && !m_check_feasible; }; + auto no_effect = [&]() { return ret == l_undef && !done() && !m_nla_satisfied && m_lemmas.empty() && m_literals.empty() && !m_check_feasible; }; if (no_effect()) m_monomial_bounds.generate_lemmas(); @@ -1329,6 +1330,9 @@ lbool core::check(unsigned level) { return l_undef; if (!m_lemmas.empty() || !m_literals.empty() || m_check_feasible) return l_false; + // bound optimization proved all monomials consistent: goal satisfied. + if (m_nla_satisfied) + return l_true; } if (no_effect() && params().arith_nl_nra_check_assignment() && m_check_assignment_fail_cnt < params().arith_nl_nra_check_assignment_max_fail()) { @@ -1563,8 +1567,11 @@ bool core::optimize_nl_bounds() { if (!lra.is_feasible()) return false; - if (lra.find_feasible_solution() == lp::lp_status::INFEASIBLE) + if (lra.find_feasible_solution() == lp::lp_status::INFEASIBLE) { + // find_feasible_solution moved the model; keep m_to_refine in sync. + init_to_refine(); return false; + } // Gather the candidate columns: every non-fixed leaf variable that // participates in a monomial (mirrors solver=2's max_min_nl_vars). @@ -1605,7 +1612,7 @@ bool core::optimize_nl_bounds() { break; for (bool is_lower : { true, false }) { rational bound; - u_dependency* dep = lra.find_improved_bound(j, is_lower, bound); + u_dependency* dep = m_monomial_bounds.improve_bound(j, is_lower, bound); if (!dep) continue; auto kind = is_lower ? lp::lconstraint_kind::GE : lp::lconstraint_kind::LE; @@ -1613,12 +1620,21 @@ bool core::optimize_nl_bounds() { } } - if (improvements.empty()) + if (improvements.empty()) { + // The exploratory simplex walk in improve_bound/mm_optimize mutated the + // LP model even though no bound was tightened. Restore a clean feasible + // model and re-calibrate m_to_refine so downstream lemma passes (grobner, + // basic_lemma) never see a stale monomial that is now consistent. + lra.find_feasible_solution(); + init_to_refine(); 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(); + // The model changed: re-calibrate m_to_refine against the new assignment. + init_to_refine(); return true; } diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h index e0db94c1ae..81976f07a5 100644 --- a/src/math/lp/nla_core.h +++ b/src/math/lp/nla_core.h @@ -90,6 +90,9 @@ class core { monomial_bounds m_monomial_bounds; unsigned m_conflicts; bool m_check_feasible = false; + // set when bound optimization re-calibrates m_to_refine to empty: every + // monomial is consistent under the optimized model, so the goal is satisfied. + bool m_nla_satisfied = false; horner m_horner; grobner m_grobner; emonics m_emons; @@ -466,6 +469,9 @@ public: return m_to_refine; } + void set_nla_satisfied() { m_nla_satisfied = true; } + bool nla_satisfied() const { return m_nla_satisfied; } + }; // end of core struct pp_mon { diff --git a/src/math/lp/nla_grobner.cpp b/src/math/lp/nla_grobner.cpp index 67e5a60508..609edf96c2 100644 --- a/src/math/lp/nla_grobner.cpp +++ b/src/math/lp/nla_grobner.cpp @@ -719,6 +719,10 @@ namespace nla { TRACE(grobner, m_solver.display(tout << "horner conflict ", e) << "\n"); return true; } + if (c().params().arith_nl_grobner_perfect_squares() && add_perfect_square_conflict(e)) { + TRACE(grobner, m_solver.display(tout << "perfect-square conflict ", e) << "\n"); + return true; + } return false; } evali.get_interval(e.poly(), i_wd); @@ -736,6 +740,156 @@ namespace nla { } } + // Helpers for perfect-square conflict detection (port of theory_arith_nl's + // is_inconsistent2). Variable lists are assumed sorted ascending. + + // Return true and set root=sqrt(coeff) if (coeff, vars) is a perfect square + // a^2 * M^2 (coeff a perfect square and every variable of even multiplicity). + static bool ps_root(rational const& coeff, unsigned_vector const& vars, rational& root) { + if (vars.size() % 2 == 1) + return false; + if (!coeff.is_perfect_square(root)) + return false; + unsigned i = 0, n = vars.size(); + while (i < n) { + unsigned v = vars[i]; + unsigned power = 1; + ++i; + for (; i < n && vars[i] == v; ++i) + ++power; + if (power % 2 == 1) + return false; + } + return true; + } + + // Return true if v12/c12 is the cross term (-2ab)*M1*M2, given that + // (v1,a^2) and (v2,b^2) are perfect squares a^2*M1^2 and b^2*M2^2. + static bool ps_cross(rational const& a, unsigned_vector const& v1, + rational const& b, unsigned_vector const& v2, + rational const& c12, unsigned_vector const& v12) { + if (!c12.is_neg()) + return false; + rational c(-2); + c *= a; + c *= b; + if (c12 != c) + return false; + unsigned n1 = v1.size(), n2 = v2.size(), n12 = v12.size(); + if (n1 + n2 != n12 * 2) + return false; + unsigned i1 = 0, i2 = 0, i12 = 0; + while (true) { + bool e1 = i1 < n1, e2 = i2 < n2, e12 = i12 < n12; + if (!e1 && !e2 && !e12) + return true; + if (!e12) + return false; + unsigned x12 = v12[i12]; + if (e1 && v1[i1] == x12) { i1 += 2; ++i12; } + else if (e2 && v2[i2] == x12) { i2 += 2; ++i12; } + else return false; + } + } + + // r <- coeff * prod(var^power) over the (sorted) variable list. + template + void grobner::monomial_interval(rational const& coeff, unsigned_vector const& vars, scoped_dep_interval& r) { + auto& di = c().m_intervals.get_dep_intervals(); + di.set_value(r, coeff); + unsigned i = 0, n = vars.size(); + while (i < n) { + lpvar v = vars[i]; + unsigned power = 1; + ++i; + for (; i < n && vars[i] == v; ++i) + ++power; + scoped_dep_interval vi(di), vp(di), tmp(di); + c().m_intervals.set_var_interval(v, vi); + di.power(vi, power, vp); + di.mul(r, vp, tmp); + di.set(r, tmp); + } + } + + // Perfect-square conflict detection: within a grobner equation e = 0, find a + // triple of monomials a^2*M1^2, b^2*M2^2, -2ab*M1*M2 forming (a*M1 - b*M2)^2. + // Replace such a triple by [0, oo) when that improves the lower bound, then + // check whether the resulting interval of e separates from zero. + bool grobner::add_perfect_square_conflict(const dd::solver::equation& e) { + vector coeffs; + vector varss; + for (auto const& [coeff, vars] : e.poly()) { + unsigned_vector vs(vars); + std::sort(vs.begin(), vs.end()); + coeffs.push_back(coeff); + varss.push_back(vs); + } + unsigned num = coeffs.size(); + if (num < 3) + return false; + auto& di = c().m_intervals.get_dep_intervals(); + svector deleted; + deleted.resize(num, false); + bool found = false; + for (unsigned i = 0; i < num; ++i) { + if (deleted[i]) + continue; + rational a; + if (!ps_root(coeffs[i], varss[i], a)) + continue; + for (unsigned j = i + 1; j < num && !deleted[i]; ++j) { + if (deleted[j]) + continue; + rational b; + if (!ps_root(coeffs[j], varss[j], b)) + continue; + for (unsigned k = i + 1; k < num; ++k) { + if (k == j || deleted[k]) + continue; + if (!ps_cross(a, varss[i], b, varss[j], coeffs[k], varss[k])) + continue; + // does [0,oo) improve on the summed interval of the triple? + scoped_dep_interval Ii(di), Ij(di), Ik(di), s1(di), s2(di); + monomial_interval(coeffs[i], varss[i], Ii); + monomial_interval(coeffs[j], varss[j], Ij); + monomial_interval(coeffs[k], varss[k], Ik); + di.add(Ii, Ij, s1); + di.add(s1, Ik, s2); + if (di.lower_is_inf(s2) || rational(di.lower(s2)).is_neg()) { + deleted[i] = deleted[j] = deleted[k] = true; + found = true; + break; + } + } + } + } + if (!found) + return false; + // Rebuild the interval of e with the perfect squares replaced by [0, oo) + // (an unconditional lower bound, hence no dependency), tracking deps of + // the remaining monomials, and check whether it separates from zero. + scoped_dep_interval sum(di); + di.reset(sum); + di.set_lower_is_inf(sum, false); + di.set_lower(sum, rational::zero()); + di.set_lower_is_open(sum, false); + di.set_upper_is_inf(sum, true); + for (unsigned i = 0; i < num; ++i) { + if (deleted[i]) + continue; + scoped_dep_interval Ii(di), tmp(di); + monomial_interval(coeffs[i], varss[i], Ii); + di.add(sum, Ii, tmp); + di.set(sum, tmp); + } + std::function f = [this](const lp::explanation& expl) { + lemma_builder lemma(m_core, "pdd-perfect-square"); + lemma &= expl; + }; + return di.check_interval_for_conflict_on_zero(sum, e.dep(), f); + } + bool grobner::propagate_linear_equations() { unsigned changed = 0; m_mon2var.clear(); diff --git a/src/math/lp/nla_grobner.h b/src/math/lp/nla_grobner.h index 4119543431..62877836dc 100644 --- a/src/math/lp/nla_grobner.h +++ b/src/math/lp/nla_grobner.h @@ -75,6 +75,12 @@ namespace nla { bool add_nla_conflict(dd::solver::equation const& eq); void check_missing_propagation(dd::solver::equation const& eq); + // perfect-square conflict detection (arith.nl.grobner_perfect_squares), + // a port of theory_arith_nl's is_inconsistent2. + bool add_perfect_square_conflict(dd::solver::equation const& eq); + template + void monomial_interval(rational const& coeff, unsigned_vector const& vars, scoped_dep_interval& r); + bool equation_is_true(dd::solver::equation const& eq); // adaptive growth (gated by arith.nl.grobner_adaptive) diff --git a/src/params/smt_params_helper.pyg b/src/params/smt_params_helper.pyg index 5dffec3a4a..3d9b54c5a3 100644 --- a/src/params/smt_params_helper.pyg +++ b/src/params/smt_params_helper.pyg @@ -93,6 +93,7 @@ def_module_params(module_name='smt', ('arith.nl.gr_q', UINT, 10, 'grobner\'s quota'), ('arith.nl.grobner_subs_fixed', UINT, 1, '0 - no subs, 1 - substitute, 2 - substitute fixed zeros only'), ('arith.nl.grobner_expand_terms', BOOL, True, 'expand terms before computing grobner basis'), + ('arith.nl.grobner_perfect_squares', BOOL, True, 'expand perfect squares with Grobner'), ('arith.nl.monomial_sandwich', BOOL, False, 'derive bound on a monomial factor by pairing two LP rows that share the other factor'), ('arith.nl.monomial_sandwich.max_fanout', UINT, 0, 'skip monomial sandwich when the conclusion factor appears in more than this many monomials (0 = no limit)'), ('arith.nl.monomial_binomial_sign', BOOL, False, 'derive bound on a binomial-monomial factor anchored on the current LP value of the monomial; replaces order_lemma_on_binomial_sign with a deterministic factor bound conditioned on a one-sided snapshot of the monomial value'),