From a12f4b9686358fdf77bd6ee1583225909408b01d Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sat, 27 Sep 2025 20:33:53 +0300 Subject: [PATCH] prepare for enforcing cheap incremental linearization axioms --- src/math/lp/lar_solver.cpp | 6 +- src/math/lp/nla_stellensatz.cpp | 185 ++++++++++++++++++++++++-------- src/math/lp/nla_stellensatz.h | 7 +- 3 files changed, 153 insertions(+), 45 deletions(-) diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index e65e0a80a..7d1faff28 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -2116,16 +2116,18 @@ namespace lp { mpq lar_solver::adjust_bound_for_int(lpvar j, lconstraint_kind& k, const mpq& bound) { if (!column_is_int(j)) return bound; - if (bound.is_int()) - return bound; switch (k) { case LT: k = LE; + if (bound.is_int()) + return bound - 1; Z3_fallthrough; case LE: return floor(bound); case GT: k = GE; + if (bound.is_int()) + return bound + 1; Z3_fallthrough; case GE: return ceil(bound); diff --git a/src/math/lp/nla_stellensatz.cpp b/src/math/lp/nla_stellensatz.cpp index c43767db6..13ec3e7a0 100644 --- a/src/math/lp/nla_stellensatz.cpp +++ b/src/math/lp/nla_stellensatz.cpp @@ -138,6 +138,103 @@ namespace nla { c().lra_solver().settings().stats().m_nla_stellensatz++; } + // + // TODO - also possible: + // derive units + // x*y = +/-x => y = +/-1 + // + // derive monotonicity: + // x >= 1 => |xy| > |y| + // + // add monotonicity axioms for squares + // + // Use to_refine and model from c().lra_solver() for initial pass. + // Use model from m_solver.lra() for subsequent passes. + // + void stellensatz::saturate_basic_linearize() { + auto &lra = c().lra_solver(); + for (auto j : c().to_refine()) { + auto const &vars = m_mon2vars[j]; + auto val_j = c().val(j); + auto val_vars = m_values[j]; + saturate_signs(j, val_j, vars, val_vars); + saturate_units(j, vars); + } + } + + // + // Sign axioms: + // x = 0 => x*y = 0 + // the equation x*y = 0 is asserted with assumption x = 0 + // x > 0 & y < 0 => x*y < 0 + // + void stellensatz::saturate_signs(lpvar j, rational const& val_j, svector const& vars, + rational const& val_vars) { + if (val_vars == 0 && val_j != 0) { + bound_justifications bounds; + lbool sign = add_bounds(vars, bounds); + SASSERT(sign == l_undef); + auto new_ci = m_solver.lra().add_var_bound(j, lp::lconstraint_kind::EQ, rational(0)); + m_new_mul_constraints.insert(new_ci, bounds); + m_ci2dep.setx(new_ci, nullptr, nullptr); + } + else if (val_j <= 0 && 0 < val_vars) { + bound_justifications bounds; + lbool sign = add_bounds(vars, bounds); + SASSERT(sign == l_false); + auto new_ci = m_solver.lra().add_var_bound(j, lp::lconstraint_kind::GT, rational(0)); + m_new_mul_constraints.insert(new_ci, bounds); + m_ci2dep.setx(new_ci, nullptr, nullptr); + } + else if (val_j >= 0 && 0 > val_vars) { + bound_justifications bounds; + lbool sign = add_bounds(vars, bounds); + SASSERT(sign == l_true); + auto new_ci = m_solver.lra().add_var_bound(j, lp::lconstraint_kind::LT, rational(0)); + m_new_mul_constraints.insert(new_ci, bounds); + m_ci2dep.setx(new_ci, nullptr, nullptr); + } + } + + // + // Unit axioms: + // x = -1 => x*y = -y + // x = 1 => x*y = y + // + void stellensatz::saturate_units(lpvar j, svector const &vars) { + lpvar non_unit = lp::null_lpvar; + bool sign = false; + for (auto v : vars) { + if (m_values[v] == 1) + continue; + if (m_values[v] == -1) { + sign = !sign; + continue; + } + if (non_unit != lp::null_lpvar) + return; + non_unit = v; + } + bound_justifications bounds; + for (auto v : vars) { + if (m_values[v] == 1 || m_values[v] == -1) { + bound b(v, lp::lconstraint_kind::EQ, rational(m_values[v])); + bounds.push_back(b); + } + } + // assert j = +/- non_unit + lp::lar_term new_lhs; + new_lhs.add_monomial(rational(1), j); + new_lhs.add_monomial(rational(sign ? 1 : -1), non_unit); + auto ti = m_solver.lra().add_term(new_lhs.coeffs_as_vector(), m_solver.lra().number_of_vars()); + m_values.push_back(m_values[j] - (sign ? 1 : -1) * m_values[non_unit]); + auto k = lp::lconstraint_kind::EQ; + auto new_ci = m_solver.lra().add_var_bound(ti, k, rational(0)); + SASSERT(m_values.size() - 1 == ti); + m_new_mul_constraints.insert(new_ci, bounds); + m_ci2dep.setx(new_ci, nullptr, nullptr); + } + // record new monomials that are created and recursively down-saturate with respect to these. // this is a simplistic pass void stellensatz::saturate_constraints() { @@ -183,13 +280,11 @@ namespace nla { // multiply by remaining vars void stellensatz::saturate_constraint(lp::constraint_index old_ci, lp::lpvar mi, lpvar x) { lp::lar_base_constraint const& con = m_solver.lra().constraints()[old_ci]; - auto &lra = c().lra_solver(); auto const& lhs = con.coeffs(); auto const& rhs = con.rhs(); auto k = con.kind(); if (k == lp::lconstraint_kind::NE || k == lp::lconstraint_kind::EQ) return; // not supported - auto sign = false, strict = true; svector vars; bool first = true; for (auto v : m_mon2vars[mi]) { @@ -200,45 +295,9 @@ namespace nla { } SASSERT(!first); // v was a member and was removed - vector bounds; + bound_justifications bounds; // compute bounds constraints and sign of vars - - auto add_bound = [&](lpvar v, lp::lconstraint_kind k, int r) { - bound b(v, k, rational(r)); - bounds.push_back(b); - }; - for (auto v : vars) { - // retrieve bounds of v - // if v has positive lower bound add as positive - // if v has negative upper bound add as negative - // otherwise look at the current value of v and add bounds assumption based on current sign. - // todo: detect squares, allow for EQ but skip bounds assumptions. - if (lra.number_of_vars() > v && lra.column_has_lower_bound(v) && lra.get_lower_bound(v).is_pos()) { - bounds.push_back(lra.get_column_lower_bound_witness(v)); - } - else if (lra.number_of_vars() > v && lra.column_has_upper_bound(v) && lra.get_upper_bound(v).is_neg()) { - bounds.push_back(lra.get_column_upper_bound_witness(v)); - sign = !sign; - } - else if (m_values[v].is_neg()) { - if (lra.var_is_int(v)) - add_bound(v, lp::lconstraint_kind::LE, -1); - else - add_bound(v, lp::lconstraint_kind::LT, 0); - sign = !sign; - } - else if (m_values[v].is_pos()) { - if (lra.var_is_int(v)) - add_bound(v, lp::lconstraint_kind::GE, 1); - else - add_bound(v, lp::lconstraint_kind::GT, 0); - } - else { - SASSERT(m_values[v] == 0); - strict = false; - add_bound(v, lp::lconstraint_kind::GE, 0); - } - } + lbool sign = add_bounds(vars, bounds); lp::lar_term new_lhs; rational new_rhs(rhs), term_value(0); for (auto const & [coeff, v] : lhs) { @@ -259,10 +318,10 @@ namespace nla { new_rhs = 0; } - if (sign) + if (sign == l_true) k = swap_side(k); - if (!strict) { + if (sign == l_undef) { switch (k) { case lp::lconstraint_kind::GT: k = lp::lconstraint_kind::GE; @@ -284,6 +343,48 @@ namespace nla { m_ci2dep.setx(new_ci, m_ci2dep.get(old_ci, nullptr), nullptr); } + lbool stellensatz::add_bounds(svector const& vars, bound_justifications& bounds) { + bool sign = false; + auto &lra = c().lra_solver(); + auto add_bound = [&](lpvar v, lp::lconstraint_kind k, int r) { + bound b(v, k, rational(r)); + bounds.push_back(b); + }; + for (auto v : vars) { + if (m_values[v] != 0) + continue; + add_bound(v, lp::lconstraint_kind::EQ, 0); + return l_undef; + } + for (auto v : vars) { + // retrieve bounds of v + // if v has positive lower bound add as positive + // if v has negative upper bound add as negative + // otherwise look at the current value of v and add bounds assumption based on current sign. + // todo: detect squares, allow for EQ but skip bounds assumptions. + if (lra.number_of_vars() > v && lra.column_has_lower_bound(v) && lra.get_lower_bound(v).is_pos()) { + bounds.push_back(lra.get_column_lower_bound_witness(v)); + } else if (lra.number_of_vars() > v && lra.column_has_upper_bound(v) && lra.get_upper_bound(v).is_neg()) { + bounds.push_back(lra.get_column_upper_bound_witness(v)); + sign = !sign; + } else if (m_values[v].is_neg()) { + if (lra.var_is_int(v)) + add_bound(v, lp::lconstraint_kind::LE, -1); + else + add_bound(v, lp::lconstraint_kind::LT, 0); + sign = !sign; + } else if (m_values[v].is_pos()) { + if (lra.var_is_int(v)) + add_bound(v, lp::lconstraint_kind::GE, 1); + else + add_bound(v, lp::lconstraint_kind::GT, 0); + } else { + UNREACHABLE(); + } + } + return sign ? l_true : l_false; + } + // Insert a (new) monomial representing product of vars. // if the length of vars is 1 then the new monomial is vars[0] // if the same monomial was previous defined, return the previously defined monomial. diff --git a/src/math/lp/nla_stellensatz.h b/src/math/lp/nla_stellensatz.h index 263ffa1de..ecf55bbf2 100644 --- a/src/math/lp/nla_stellensatz.h +++ b/src/math/lp/nla_stellensatz.h @@ -34,9 +34,10 @@ namespace nla { rational rhs; }; using bound_justification = std::variant; + using bound_justifications = vector; coi m_coi; - u_map> m_new_mul_constraints; + u_map m_new_mul_constraints; indexed_uint_set m_to_refine; ptr_vector m_ci2dep; vector m_values; @@ -60,8 +61,12 @@ namespace nla { lpvar add_monomial(svector const& vars); bool is_int(svector const& vars) const; 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_basic_linearize(); + void saturate_signs(lpvar j, rational const& val_j, svector const& vars, rational const& val_vars); + void saturate_units(lpvar j, svector const &vars); // lemmas