From c621f5974088f455f5dc3b5ba2cdf45cc9c7f444 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sun, 28 Sep 2025 14:36:53 +0300 Subject: [PATCH] fix bug with saturation of monotonicity, and add more general case for downward saturation Signed-off-by: Nikolaj Bjorner --- src/math/lp/nla_stellensatz.cpp | 78 ++++++++++++++++++++++++--------- src/math/lp/nla_stellensatz.h | 8 ++-- 2 files changed, 62 insertions(+), 24 deletions(-) diff --git a/src/math/lp/nla_stellensatz.cpp b/src/math/lp/nla_stellensatz.cpp index f399afe2b..151b46eff 100644 --- a/src/math/lp/nla_stellensatz.cpp +++ b/src/math/lp/nla_stellensatz.cpp @@ -46,6 +46,7 @@ namespace nla { init_solver(); saturate_constraints(); saturate_basic_linearize(); + TRACE(arith, display(tout << "stellensatz after saturation\n")); lbool r = m_solver.solve(ex); if (r == l_false) add_lemma(ex); @@ -137,7 +138,7 @@ namespace nla { } // - // a constraint can be explained by a set of bounds used as assumptions for the constraint + // a constraint can be explained by a set of bounds used as assumptions // and by an original constraint. // void stellensatz::explain_constraint(lemma_builder& new_lemma, lp::constraint_index ci, lp::explanation& ex) { @@ -156,7 +157,13 @@ namespace nla { if (std::holds_alternative(b)) { auto dep = *std::get_if(&b); m_solver.lra().push_explanation(dep, ex); - } else { + } + else { + // + // the inequality was posted as an assumption + // negate it to add to the lemma + // recall that lemmas are represented in the form: /\ Assumptions => \/ C + // auto [v, k, rhs] = *std::get_if(&b); k = negate(k); if (m_solver.lra().var_is_int(v)) { @@ -296,13 +303,13 @@ namespace nla { auto valx = m_values[x]; auto valy = m_values[y]; if (valx > 1 && valy > 0 && val_j <= valy) - saturate_monotonicity(j, val_j, x, false, y, false); + saturate_monotonicity(j, val_j, x, 1, y, 1); else if (valx > 1 && valy < 0 && -val_j <= -valy) - saturate_monotonicity(j, val_j, x, false, y, true); + saturate_monotonicity(j, val_j, x, 1, y, -1); else if (valx < -1 && valy > 0 && -val_j <= valy) - saturate_monotonicity(j, val_j, x, true, y, false); + saturate_monotonicity(j, val_j, x, -1, y, 1); else if (valx < -1 && valy < 0 && val_j <= -valy) - saturate_monotonicity(j, val_j, x, true, y, true); + saturate_monotonicity(j, val_j, x, -1, y, -1); } // x > 1, y > 0 => xy > y @@ -349,6 +356,7 @@ namespace nla { lp::constraint_index stellensatz::add_ineq(bound_justifications const& bounds, lp::lar_term const& t, lp::lconstraint_kind k, rational const& rhs) { + SASSERT(!t.coeffs_as_vector().empty()); auto ti = m_solver.lra().add_term(t.coeffs_as_vector(), m_solver.lra().number_of_vars()); m_values.push_back(value(t)); auto new_ci = m_solver.lra().add_var_bound(ti, k, rhs); @@ -369,7 +377,8 @@ namespace nla { // record new monomials that are created and recursively down-saturate with respect to these. // this is a simplistic pass void stellensatz::saturate_constraints() { - vector> var2cs; + vector> var2cs; // cs contain a term using v + vector> vars2cs; // cs contain a term with a monomial using v // current approach: only resolve against var2cs, which is initialized // with monomials in the input. @@ -380,25 +389,55 @@ namespace nla { if (v >= var2cs.size()) var2cs.resize(v + 1); var2cs[v].push_back(ci); + if (m_mon2vars.contains(v)) { + for (auto w : m_mon2vars[v]) { + if (w >= vars2cs.size()) + vars2cs.resize(w + 1); + vars2cs[w].push_back(ci); + } + } } // insert monomials to be refined insert_monomials_from_constraint(ci); } + auto is_subset = [&](svector const &a, svector const& b) { + if (a.size() >= b.size()) + return false; + // check if a is a subset of b, counting multiplicies, assume a, b are sorted + unsigned i = 0, j = 0; + while (i < a.size() && j < b.size()) { + if (a[i] == b[j]) + ++i; + ++j; + } + return i == a.size(); + }; + for (unsigned it = 0; it < m_to_refine.size(); ++it) { auto j = m_to_refine[it]; auto vars = m_mon2vars[j]; for (auto v : vars) { if (v >= var2cs.size()) continue; - auto cs = var2cs[v]; - for (auto ci : cs) { + svector _vars; + _vars.push_back(v); + for (auto ci : var2cs[v]) { for (auto [coeff, u] : m_solver.lra().constraints()[ci].coeffs()) { if (u == v) - saturate_constraint(ci, j, v); + saturate_constraint(ci, j, _vars); } } } + for (auto v : vars) { + if (v >= vars2cs.size()) + continue; + for (auto ci : vars2cs[v]) { + for (auto [coeff, u] : m_solver.lra().constraints()[ci].coeffs()) + if (m_mon2vars.contains(u) && is_subset(m_mon2vars[u], vars)) + saturate_constraint(ci, j, m_mon2vars[u]); + } + } } IF_VERBOSE(5, c().lra_solver().display(verbose_stream() << "original\n"); @@ -407,8 +446,8 @@ namespace nla { } // multiply by remaining vars - void stellensatz::saturate_constraint(lp::constraint_index old_ci, lpvar mi, lpvar x) { - resolvent r = {old_ci, mi, x}; + void stellensatz::saturate_constraint(lp::constraint_index old_ci, lpvar mi, svector const& xs) { + resolvent r = {old_ci, mi, xs}; if (m_resolvents.contains(r)) return; m_resolvents.insert(r); @@ -419,15 +458,14 @@ namespace nla { if (k == lp::lconstraint_kind::NE || k == lp::lconstraint_kind::EQ) return; // not supported - svector vars; - bool first = true; - for (auto v : m_mon2vars[mi]) { - if (v != x || !first) - vars.push_back(v); - else - first = false; + + // xs is a proper subset of vars in mi + svector vars(m_mon2vars[mi]); + for (auto x : xs) { + SASSERT(vars.contains(x)); + vars.erase(x); } - SASSERT(!first); // v was a member and was removed + SASSERT(!vars.empty()); bound_justifications bounds; // compute bounds constraints and sign of vars diff --git a/src/math/lp/nla_stellensatz.h b/src/math/lp/nla_stellensatz.h index 5ee07eb48..6b4683a53 100644 --- a/src/math/lp/nla_stellensatz.h +++ b/src/math/lp/nla_stellensatz.h @@ -58,15 +58,15 @@ namespace nla { struct resolvent { lp::constraint_index old_ci; lpvar mi; - lpvar x; + svector xs; struct eq { bool operator()(resolvent const& a, resolvent const& b) const { - return a.old_ci == b.old_ci && a.mi == b.mi && a.x == b.x; + return a.old_ci == b.old_ci && a.mi == b.mi && a.xs == b.xs; } }; struct hash { unsigned operator()(resolvent const& a) const { - return hash_u_u(a.old_ci, hash_u_u(a.mi, a.x)); + return hash_u_u(a.old_ci, hash_u_u(a.mi, svector_hash()(a.xs))); } }; }; @@ -92,7 +92,7 @@ namespace nla { lpvar add_var(bool is_int); lbool add_bounds(svector const &vars, bound_justifications &bounds); void saturate_constraints(); - void saturate_constraint(lp::constraint_index con_id, lp::lpvar mi, lpvar x); + void saturate_constraint(lp::constraint_index con_id, lp::lpvar mi, svector const & xs); void saturate_basic_linearize(); void saturate_basic_linearize(lpvar j, rational const &val_j, svector const &vars, rational const &val_vars);