mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-11-25 15:09:32 +00:00
Code Activity

revamp stellensatz to use polynomals

Browse source 
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2025-11-12 20:21:42 -08:00
parent 21b36868b5
commit c7084e9998
4 changed files with 397 additions and 542 deletions
Download patch file Download diff file Expand all files Collapse all files

64
src/math/dd/dd_pdd.cpp
View file

@ -1908,6 +1908,36 @@ namespace dd {
return m->mk_var(var())*h + l;
}
void pdd::merge(unsigned_vector& lo, unsigned_vector& hi, unsigned_vector& common) {
std::sort(lo.begin(), lo.end());
std::sort(hi.begin(), hi.end());
unsigned ir = 0, jr = 0;
for (unsigned i = 0, j = 0; i < lo.size() || j < hi.size(); ) {
if (i == lo.size())
hi[jr++] = hi[j++];
else if (j == hi.size())
lo[ir++] = lo[i++];
else if (lo[i] == hi[j]) {
common.push_back(lo[i]);
++i;
++j;
}
else if (lo[i] > hi[j])
hi[jr++] = hi[j++];
else
lo[ir++] = lo[i++];
}
lo.shrink(ir);
hi.shrink(jr);
}
pdd pdd::mul(unsigned_vector const& vars) const {
pdd r = *this;
for (auto v : vars)
r = m->mk_var(v) * r;
return r;
}
std::pair<unsigned_vector, pdd> pdd::var_factors() const {
if (is_val())
return { unsigned_vector(), *this };
@ -1924,49 +1954,23 @@ namespace dd {
return { unsigned_vector(), *this };
unsigned_vector lo_and_hi;
auto merge = [&](unsigned_vector& lo_vars, unsigned_vector& hi_vars) {
unsigned ir = 0, jr = 0;
for (unsigned i = 0, j = 0; i < lo_vars.size() || j < hi_vars.size(); ) {
if (i == lo_vars.size())
hi_vars[jr++] = hi_vars[j++];
else if (j == hi_vars.size())
lo_vars[ir++] = lo_vars[i++];
else if (lo_vars[i] == hi_vars[j]) {
lo_and_hi.push_back(lo_vars[i]);
++i;
++j;
}
else if (m->m_var2level[lo_vars[i]] > m->m_var2level[hi_vars[j]])
hi_vars[jr++] = hi_vars[j++];
else
lo_vars[ir++] = lo_vars[i++];
}
lo_vars.shrink(ir);
hi_vars.shrink(jr);
};
auto mul = [&](unsigned_vector const& vars, pdd p) {
for (auto v : vars)
p *= m->mk_var(v);
return p;
};
auto [hi_vars, p] = hi().var_factors();
if (lo_vars.back() == v) {
lo_vars.pop_back();
merge(lo_vars, hi_vars);
merge(lo_vars, hi_vars, lo_and_hi);
lo_and_hi.push_back(v);
return { lo_and_hi, mul(lo_vars, q) + mul(hi_vars, p) };
return { lo_and_hi, q.mul(lo_vars) + p.mul(hi_vars) };
}
if (hi_vars.empty())
return { unsigned_vector(), *this };
merge(lo_vars, hi_vars);
merge(lo_vars, hi_vars, lo_and_hi);
hi_vars.push_back(v);
if (lo_and_hi.empty())
return { unsigned_vector(), *this };
else
return { lo_and_hi, mul(lo_vars, q) + mul(hi_vars, p) };
return { lo_and_hi, q.mul(lo_vars) + p.mul(hi_vars) };
}

2
src/math/dd/dd_pdd.h
View file

@ -483,6 +483,8 @@ namespace dd {
* \brief factor out variables
*/
std::pair<unsigned_vector, pdd> var_factors() const;
static void merge(unsigned_vector &lo, unsigned_vector &hi, unsigned_vector& common);
pdd mul(unsigned_vector const &vars) const;
pdd subst_val0(vector<std::pair<unsigned, rational>> const& s) const { return m->subst_val0(*this, s); }
pdd subst_val(pdd const& s) const { VERIFY_EQ(m, s.m); return m->subst_val(*this, s); }

765
src/math/lp/nla_stellensatz.cpp
View file

@ -6,6 +6,7 @@
*/
#include "util/heap.h"
#include "math/dd/pdd_eval.h"
#include "math/lp/nla_core.h"
#include "math/lp/nla_coi.h"
#include "math/lp/nla_stellensatz.h"
@ -29,16 +30,37 @@ namespace nla {
return k1;
}
lp::lconstraint_kind meet(lp::lconstraint_kind k1, lp::lconstraint_kind k2) {
if (k1 == k2)
return k1;
if (k1 == lp::lconstraint_kind::EQ)
return k1;
if (k2 == lp::lconstraint_kind::EQ)
return k2;
if ((k1 == lp::lconstraint_kind::LE && k2 == lp::lconstraint_kind::LT) ||
(k1 == lp::lconstraint_kind::LT && k2 == lp::lconstraint_kind::LE))
return lp::lconstraint_kind::LE;
if ((k1 == lp::lconstraint_kind::GE && k2 == lp::lconstraint_kind::GT) ||
(k1 == lp::lconstraint_kind::GT && k2 == lp::lconstraint_kind::GE))
return lp::lconstraint_kind::GE;
UNREACHABLE();
return k1;
}
stellensatz::stellensatz(core* core) :
common(core),
m_solver(*this),
m_coi(*core)
m_coi(*core),
pddm(m_solver.lra().number_of_vars())
{}
lbool stellensatz::saturate() {
init_solver();
TRACE(arith, display(tout << "stellensatz after saturation\n"));
TRACE(arith, display(tout << "stellensatz before saturation\n"));
eliminate_variables();
// TODO: populate solver
TRACE(arith, display(tout << "stellensatz after saturation\n"));
return l_undef;
lbool r = m_solver.solve();
// IF_VERBOSE(0, verbose_stream() << "stellensatz " << r << "\n");
if (r == l_false)
@ -50,9 +72,7 @@ namespace nla {
m_solver.init();
m_vars2mon.reset();
m_mon2vars.reset();
m_values.reset();
m_justifications.reset();
m_max_monomial_degree = 0;
m_constraints.reset();
m_coi.init();
init_vars();
init_occurs();
@ -60,34 +80,18 @@ namespace nla {
void stellensatz::init_vars() {
auto const& src = c().lra_solver();
auto &dst = m_solver.lra();
auto sz = src.number_of_vars();
for (unsigned v = 0; v < sz; ++v) {
// Declare v into m_solver.lra()
unsigned w;
if (m_coi.mons().contains(v))
init_monomial(v);
else
m_values.push_back(c().val(v));
if (m_coi.terms().contains(v)) {
auto const& t = src.get_term(v);
// justification: variables in coefficients are always declared before term variable.
SASSERT(all_of(t, [&](auto p) { return p.j() < v; }));
w = dst.add_term(t.coeffs_as_vector(), v);
}
else
w = dst.add_var(v, src.var_is_int(v));
// assert bounds on v in the new solver.
VERIFY(w == v);
if (src.column_has_lower_bound(v)) {
auto lo_bound = src.get_lower_bound(v);
SASSERT(lo_bound.y >= 0);
auto k = lo_bound.y > 0 ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE;
auto rhs = lo_bound.x;
auto dep = src.get_column_lower_bound_witness(v);
auto ci = dst.add_var_bound(v, k, rhs);
add_justification(ci, external_justification(dep));
add_var_bound(v, k, rhs, external_justification(dep));
}
if (src.column_has_upper_bound(v)) {
auto hi_bound = src.get_upper_bound(v);
@ -95,13 +99,57 @@ namespace nla {
auto k = hi_bound.y < 0 ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE;
auto rhs = hi_bound.x;
auto dep = src.get_column_upper_bound_witness(v);
auto ci = dst.add_var_bound(v, k, rhs);
add_justification(ci, external_justification(dep));
add_var_bound(v, k, rhs, external_justification(dep));
}
}
m_occurs.reserve(sz);
for (auto const &[v, vars] : m_mon2vars)
m_max_monomial_degree = std::max(m_max_monomial_degree, vars.size());
}
dd::pdd stellensatz::to_pdd(lpvar v) {
if (m_coi.mons().contains(v)) {
auto& mon = c().emons()[v];
dd::pdd prod(pddm);
prod = 1;
for (auto w : mon.vars())
prod *= to_pdd(w);
return prod;
}
if (m_coi.terms().contains(v)) {
auto const& t = c().lra_solver().get_term(v);
dd::pdd sum(pddm);
sum = 0;
for (auto cv : t)
sum *= to_pdd(cv.j())*cv.coeff();
return sum;
}
return pddm.mk_var(v);
}
stellensatz::term_offset stellensatz::to_term_offset(dd::pdd const &p) {
term_offset to;
for (auto const &[coeff, vars] : p) {
if (vars.empty())
to.second += coeff;
else
to.first.add_monomial(coeff, mk_monomial(vars));
}
return to;
}
lpvar stellensatz::mk_monomial(svector<lpvar> const &vars) {
lpvar v;
if (vars.size() == 1)
return vars[0];
svector<lpvar> _vars(vars);
std::sort(_vars.begin(), _vars.end());
if (m_vars2mon.find(_vars, v))
return v;
v = add_var(is_int(vars));
m_vars2mon.insert(_vars, v);
m_mon2vars.insert(v, _vars);
SASSERT(m_values.size() == v);
return v;
}
void stellensatz::init_monomial(unsigned mon_var) {
@ -110,16 +158,37 @@ namespace nla {
std::sort(vars.begin(), vars.end());
m_vars2mon.insert(vars, mon_var);
m_mon2vars.insert(mon_var, vars);
m_values.push_back(cvalue(vars));
}
lp::constraint_index stellensatz::add_var_bound(lp::lpvar v, lp::lconstraint_kind k, rational const& rhs, justification j) {
auto p = to_pdd(v) - rhs;
rational d(1);
for (auto const& [c, is_constant] : p.coefficients())
d = lcm(d, denominator(c));
if (d != 1)
p *= d;
return gcd_normalize(add_constraint(p, k, j));
}
lp::constraint_index stellensatz::add_constraint(dd::pdd p, lp::lconstraint_kind k, justification j) {
if (k == lp::lconstraint_kind::LE) {
k = lp::lconstraint_kind::GE;
p = -p;
}
if (k == lp::lconstraint_kind::LT) {
k = lp::lconstraint_kind::GT;
p = -p;
}
m_constraints.push_back({p, k, j });
return m_constraints.size() - 1;
}
// initialize active set of constraints that evaluate to false in the current model
// loop over active set to eliminate variables.
void stellensatz::eliminate_variables() {
vector<std::pair<lp::constraint_index, uint_set>> active;
for (auto ci : m_solver.lra().constraints().indices()) {
if (constraint_is_true(ci))
continue;
for (unsigned ci = 0; ci < m_constraints.size(); ++ci) {
if (!constraint_is_true(ci))
active.push_back({ci, uint_set()});
}
while (!active.empty()) {
@ -133,121 +202,121 @@ namespace nla {
auto x = select_variable_to_eliminate(ci);
auto f = factor(x, ci);
auto p_value = cvalue(f.p);
if (!f.p.is_empty() && p_value == 0) {
if (!f.p.is_zero() && p_value == 0) {
// add p = 0 as assumption and reduce to q.
auto p_is_0 = assume(f.p, lp::lconstraint_kind::EQ, rational(0));
if (!f.q.is_empty()) {
auto const &con = m_solver.lra().constraints()[ci];
// ci & -p_is_0*f.coeff*x^f.degree => new_ci
if (f.coeff != -1)
p_is_0 = multiply(p_is_0, rational(-f.coeff));
auto p_is_0 = assume(f.p, lp::lconstraint_kind::EQ);
if (!f.q.is_zero()) {
// ci & -p_is_0*x^f.degree => new_ci
p_is_0 = multiply(p_is_0, rational(-1));
for (unsigned i = 0; i < f.degree; ++i)
p_is_0 = multiply_var(p_is_0, x);
auto new_ci = add(ci, p_is_0);
init_occurs(new_ci);
uint_set new_tabu(tabu);
uint_set q_vars;
for (auto cv : f.q)
q_vars.insert(cv.j());
for (auto cv : f.p)
if (!q_vars.contains(cv.j()))
new_tabu.insert(cv.j());
for (auto j : f.q.free_vars())
q_vars.insert(j);
for (auto j : f.p.free_vars())
if (!q_vars.contains(j))
new_tabu.insert(j);
active.push_back({new_ci, new_tabu});
}
continue;
}
for (auto other_ci : m_occurs[x]) {
}
if (other_ci == ci)
continue;
auto const &other_ineq = m_constraints[other_ci];
auto g = factor(x, other_ci);
auto p_value2 = cvalue(g.p);
// check that p_value and p_value2 have different signs
// check that other_ineq.lhs() is disjoint from tabu
// compute overlaps extending x
//
SASSERT(g.degree > 0);
if (g.degree > f.degree) // future: could handle this too by considering tabu to be a map into degrees.
continue;
if (p_value > 0 && p_value2 > 0)
continue;
if (p_value < 0 && p_value2 < 0)
continue;
if (any_of(other_ineq.p.free_vars(), [&](unsigned j) { return tabu.contains(j); }))
continue;
for (unsigned i = 0; i < f.degree; ++i)
f.p *= to_pdd(x);
for (unsigned i = 0; i < g.degree; ++i)
g.p *= to_pdd(x);
auto [m1, f_p] = f.p.var_factors();
auto [m2, g_p] = g.p.var_factors();
unsigned_vector m1m2;
dd::pdd::merge(m1, m2, m1m2);
SASSERT(m1m2.contains(x));
f_p = f_p.mul(m1);
g_p = g_p.mul(m2);
auto sign_f = cvalue(f_p) < 0;
SASSERT(sign_f != cvalue(g_p) < 0);
SASSERT(cvalue(f_p) != 0);
SASSERT(cvalue(g_p) != 0);
// m1m2 * f_p + f.q >= 0
// m1m2 * g_p + g.q >= 0
//
auto ci_a = ci;
auto ci_b = other_ci;
auto aj = assumption_justification();
if (f_p.is_val())
ci_a = multiply(other_ci, f_p.val());
else if (sign_f)
ci_a = multiply(other_ci, add_constraint(f_p, lp::lconstraint_kind::LT, aj));
else
ci_a = multiply(other_ci, add_constraint(f_p, lp::lconstraint_kind::GT, aj));
if (g_p.is_val())
ci_b = multiply(ci, g_p.val());
else if (!sign_f)
ci_b = multiply(ci, add_constraint(g_p, lp::lconstraint_kind::LT, aj));
else
ci_b = multiply(ci, add_constraint(g_p, lp::lconstraint_kind::GT, aj));
auto new_ci = add(ci_a, ci_b);
init_occurs(new_ci);
if (!constraint_is_true(new_ci)) {
auto const &p = m_constraints[ci].p;
auto const &[new_p, new_k, new_j] = m_constraints[new_ci];
uint_set new_tabu(tabu), fv;
for (auto v : new_p.free_vars())
fv.insert(v);
for (auto v : p.free_vars())
if (!fv.contains(v))
new_tabu.insert(v);
active.push_back({new_ci, new_tabu});
}
}
}
}
lp::lpvar stellensatz::select_variable_to_eliminate(lp::constraint_index ci) {
auto &ineq = m_solver.lra().constraints()[ci];
auto const& [p, k, j] = m_constraints[ci];
lpvar best_var = lp::null_lpvar;
for (auto cv : ineq.lhs()) {
auto v = cv.j();
if (best_var == lp::null_lpvar) {
for (auto v : p.free_vars())
if (best_var > v)
best_var = v;
continue;
}
if (is_mon_var(v)) {
for (auto w : m_mon2vars[v]) {
if (best_var > v) {
best_var = v;
}
}
}
else if (best_var > v) {
best_var = v;
}
}
return best_var;
}
unsigned stellensatz::degree_of_var_in_constraint(lpvar var, lp::constraint_index ci) const {
auto &ineq = m_solver.lra().constraints()[ci];
unsigned max_degree = 0;
for (auto cv : ineq.lhs()) {
auto degree = degree_of_var_in_monomial(var, cv.j());
if (degree > max_degree)
max_degree = degree;
}
return max_degree;
}
unsigned stellensatz::degree_of_var_in_monomial(lpvar v, lpvar mi) const {
unsigned degree = 0;
if (is_mon_var(mi)) {
for (auto w : m_mon2vars[mi])
if (v == w)
degree++;
}
else if (mi == v)
++degree;
return degree;
return m_constraints[ci].p.degree(var);
}
stellensatz::factorization stellensatz::factor(lpvar v, lp::constraint_index ci) {
auto &ineq = m_solver.lra().constraints()[ci];
auto const& [p, k, j] = m_constraints[ci];
auto degree = degree_of_var_in_constraint(v, ci);
rational coeff(0);
lp::lar_term p, q;
for (auto cv : ineq.lhs()) {
auto d = degree_of_var_in_monomial(v, cv.j());
if (d < degree)
q.add_monomial(cv.coeff(), cv.j());
else {
auto w = divide(v, degree, cv.j());
if (w == lp::null_lpvar)
coeff += cv.coeff();
else
p.add_monomial(cv.coeff(), w);
}
}
return {coeff, degree, p, q};
}
lp::lpvar stellensatz::divide(lpvar v, unsigned degree, lpvar mi) {
SASSERT(degree_of_var(v, mi) >= degree);
if (is_mon_var(mi)) {
auto const &vars = m_mon2vars[mi];
if (degree == vars.size())
return lp::null_lpvar;
svector<lpvar> _vars;
for (auto w : vars) {
if (w == v && degree > 0)
--degree;
else
_vars.push_back(w);
}
return mk_monomial(_vars);
}
SASSERT(degree == 1 && mi == v);
return lp::null_lpvar;
dd::pdd lc(pddm), rest(pddm);
p.factor(v, degree, lc, rest);
return {degree, lc, rest};
}
//
@ -282,10 +351,7 @@ namespace nla {
if (m_constraints_in_conflict.contains(ci))
return;
m_constraints_in_conflict.insert(ci);
auto just = m_justifications.get(ci);
if (just == nullptr)
return;
auto &b = *just;
auto &[p, k, b] = m_constraints[ci];
if (std::holds_alternative<external_justification>(b)) {
auto dep = std::get<external_justification>(b);
m_solver.lra().push_explanation(dep.dep, ex);
@ -305,8 +371,8 @@ namespace nla {
explain_constraint(new_lemma, m.ci, ex);
}
else if (std::holds_alternative<assumption_justification>(b)) {
auto &con = m_solver.lra().constraints()[ci];
new_lemma |= ineq(con.lhs(), negate(con.kind()), con.rhs());
auto [t, coeff] = to_term_offset(p);
new_lemma |= ineq(t, negate(k), -coeff);
}
else if (std::holds_alternative<gcd_justification>(b)) {
auto& m = std::get<gcd_justification>(b);
@ -316,297 +382,98 @@ namespace nla {
UNREACHABLE();
}
rational stellensatz::cvalue(lp::lar_term const &t) const {
rational r(0);
for (auto cv : t)
r += c().val(cv.j()) * cv.coeff();
return r;
rational stellensatz::cvalue(dd::pdd const& p) const {
dd::pdd_eval eval;
// eval.var2val() = [&](unsigned v) -> void { return c().val(v); };
return eval(p);
}
rational stellensatz::mvalue(lp::lar_term const &t) const {
rational r(0);
for (auto cv : t)
r += m_values[cv.j()] * cv.coeff();
return r;
}
rational stellensatz::cvalue(svector<lpvar> const &prod) const {
rational p(1);
for (auto v : prod)
p *= c().val(v);
return p;
}
rational stellensatz::mvalue(svector<lpvar> const &prod) const {
rational p(1);
for (auto v : prod)
p *= m_values[v];
return p;
}
void stellensatz::gcd_normalize(vector<std::pair<rational, lpvar>> &t, lp::lconstraint_kind k, rational &rhs) {
rational d(1);
for (auto &[c,v] : t)
d = lcm(d, denominator(c));
d = lcm(d, denominator(rhs));
if (d != 1) {
for (auto &[c, v] : t)
c *= d;
rhs *= d;
}
lp::constraint_index stellensatz::gcd_normalize(lp::constraint_index ci) {
auto [p, k, j] = m_constraints[ci];
rational g(0);
for (auto &[c, v] : t)
bool is_int = all_of(p.free_vars(), [&](unsigned v) { return c().lra_solver().var_is_int(v); });
for (auto const& [c, is_constant] : p.coefficients())
if (!is_constant || !is_int)
g = gcd(g, c);
bool is_int = true;
for (auto &[c, v] : t)
is_int &= m_solver.lra().var_is_int(v);
if (!is_int)
g = gcd(g, rhs);
if (g == 1 || g == 0)
return;
if (g == 0 || g == 1)
return ci;
switch (k) {
case lp::lconstraint_kind::GE:
for (auto &[c, v] : t)
c /= g;
rhs /= g;
rhs = ceil(rhs);
break;
case lp::lconstraint_kind::GT:
for (auto &[c, v] : t)
c /= g;
if (is_int)
rhs += 1;
rhs /= g;
rhs = ceil(rhs);
break;
case lp::lconstraint_kind::LE:
for (auto &[c, v] : t)
c /= g;
rhs /= g;
rhs = floor(rhs);
break;
case lp::lconstraint_kind::GE: {
auto q = p * (1/ g);
q += (ceil(q.offset()) - q.offset());
return add_constraint(q, k, gcd_justification(ci));
}
case lp::lconstraint_kind::GT: {
auto q = p;
if (is_int) {
q += rational(1);
k = lp::lconstraint_kind::GE;
}
q *= (1 / g);
q += (ceil(q.offset()) - q.offset());
return add_constraint(q, k, gcd_justification(ci));
}
case lp::lconstraint_kind::LT:
for (auto &[c, v] : t)
c /= g;
if (is_int)
rhs -= 1;
rhs /= g;
rhs = floor(rhs);
break;
case lp::lconstraint_kind::LE:
UNREACHABLE();
case lp::lconstraint_kind::EQ:
case lp::lconstraint_kind::NE:
if (!divides(g, rhs))
break;
for (auto &[c, v] : t)
c /= g;
rhs /= g;
break;
}
}
lp::constraint_index stellensatz::normalize_ineq(lp::constraint_index ci) {
auto const &con = m_solver.lra().constraints()[ci];
auto k = con.kind();
if (k == lp::lconstraint_kind::LE || k == lp::lconstraint_kind::LT) {
auto mjust = multiplication_const_justification({ci, rational(-1)});
k = swap_side(k);
rational value(0);
auto coeffs = con.coeffs();
for (auto &[c, v] : coeffs) {
c = -c;
value += c * m_values[v];
}
auto ti = m_solver.lra().add_term(coeffs, m_solver.lra().number_of_vars());
m_occurs.reserve(m_solver.lra().number_of_vars());
m_values.push_back(value);
ci = m_solver.lra().add_var_bound(ti, k, -con.rhs());
}
if (!divides(g, p.offset()))
return ci;
return add_constraint(p * (1/g), k, gcd_justification(ci));
default:
UNREACHABLE();
return ci;
}
lp::constraint_index stellensatz::add_ineq(
justification const& just, lp::lar_term & t, lp::lconstraint_kind k,
rational const & rhs_) {
auto coeffs = t.coeffs_as_vector();
SASSERT(!coeffs.empty());
auto add = [&](justification const& just, vector<std::pair<rational, lpvar>> const& coeffs, rational const& rhs) {
auto ti = m_solver.lra().add_term(coeffs, m_solver.lra().number_of_vars());
m_occurs.reserve(m_solver.lra().number_of_vars());
m_values.push_back(mvalue(t));
auto new_ci = m_solver.lra().add_var_bound(ti, k, rhs);
TRACE(arith, display(tout, just) << "\n");
SASSERT(m_values.size() - 1 == ti);
add_justification(new_ci, just);
return new_ci;
};
auto new_ci = add(just, coeffs, rhs_);
rational rhs(rhs_);
gcd_normalize(coeffs, k, rhs);
if (rhs != rhs_)
new_ci = add(gcd_justification(new_ci), coeffs, rhs);
init_occurs(new_ci);
return new_ci;
}
lp::constraint_index stellensatz::add_ineq(justification const& just, lpvar j, lp::lconstraint_kind k,
rational const& rhs) {
auto new_ci = m_solver.lra().add_var_bound(j, k, rhs);
add_justification(new_ci, just);
init_occurs(new_ci);
return new_ci;
lp::constraint_index stellensatz::assume(dd::pdd const& p, lp::lconstraint_kind k) {
return add_constraint(p, k, assumption_justification());
}
lp::constraint_index stellensatz::assume(lp::lar_term &t, lp::lconstraint_kind k, rational const &rhs) {
return add_ineq(assumption_justification(), t, k, rhs);
}
// p1 >= lo1, p2 >= lo2 => p1 + p2 >= lo1 + lo2
lp::constraint_index stellensatz::add(lp::constraint_index left, lp::constraint_index right) {
auto const &lcon = m_solver.lra().constraints()[left];
auto const &rcon = m_solver.lra().constraints()[right];
SASSERT(lcon.kind() == lp::lconstraint_kind::GE || lcon.kind() == lp::lconstraint_kind::GT);
SASSERT(rcon.kind() == lp::lconstraint_kind::GE || rcon.kind() == lp::lconstraint_kind::GT);
lp::lar_term t;
for (auto const &[v, c] : lcon.lhs())
t.add_monomial(c, v);
for (auto const &[v, c] : rcon.lhs())
t.add_monomial(c, v);
lp::lconstraint_kind k = join(lcon.kind(), rcon.kind());
rational rhs = lcon.rhs() + rcon.rhs();
auto just = addition_justification{left, right};
return add_ineq(just, t, k, rhs);
auto const &[p, k1, j1] = m_constraints[left];
auto const &[q, k2, j2] = m_constraints[right];
lp::lconstraint_kind k = join(k1, k2);
return gcd_normalize(add_constraint(p + q, k, addition_justification{left, right}));
}
// p >= lo => a * p >= a * lo if a > 0
lp::constraint_index stellensatz::multiply(lp::constraint_index ci, rational const& mul) {
SASSERT(mul != 0);
auto const &con = m_solver.lra().constraints()[ci];
lp::lar_term t;
for (auto const &[v, c] : con.lhs())
t.add_monomial(c * mul, v);
lp::lconstraint_kind k = con.kind();
if (mul < 0)
k = swap_side(k);
rational rhs = con.rhs() * mul;
auto just = multiplication_const_justification{ci, mul};
return add_ineq(just, t, k, rhs);
auto const& [p, k, j] = m_constraints[ci];
return add_constraint(p * mul, mul > 0 ? k : swap_side(k), multiplication_const_justification{ci, mul});
}
lp::constraint_index stellensatz::multiply(lp::constraint_index left, lp::constraint_index right) {
auto const &lcon = m_solver.lra().constraints()[left];
auto const &rcon = m_solver.lra().constraints()[right];
SASSERT(lcon.kind() == lp::lconstraint_kind::GE || lcon.kind() == lp::lconstraint_kind::GT);
SASSERT(rcon.kind() == lp::lconstraint_kind::GE || rcon.kind() == lp::lconstraint_kind::GT);
lp::lar_term t;
for (auto const &[v1, c1] : lcon.lhs()) {
for (auto const &[v2, c2] : rcon.lhs()) {
auto j = mk_monomial(expand_monomial({v1, v2}));
t.add_monomial(c1 * c2, j);
}
}
if (rcon.rhs() != 0)
for (auto const &[v, c] : lcon.lhs())
t.add_monomial(-c * rcon.rhs(), v);
if (lcon.rhs() != 0)
for (auto const &[v, c] : rcon.lhs())
t.add_monomial(-c * lcon.rhs(), v);
lp::lconstraint_kind k = join(lcon.kind(), rcon.kind());
rational rhs{0};
return add_ineq(multiplication_justification{left, right}, t, k, rhs);
auto const &[p, k1, j1] = m_constraints[left];
auto const &[q, k2, j2] = m_constraints[right];
lp::lconstraint_kind k = meet(k1, k2);
return add_constraint(p*q, k, multiplication_justification{left, right});
}
lp::constraint_index stellensatz::multiply_var(lp::constraint_index ci, lpvar x) {
auto const &con = m_solver.lra().constraints()[ci];
lp::lconstraint_kind k = con.kind();
auto const& [p, k, j] = m_constraints[ci];
SASSERT(k == lp::lconstriant_kind::EQ);
SASSERT(con.rhs() == 0);
lp::lar_term t;
for (auto const &[v, c] : con.lhs())
t.add_monomial(c, mk_monomial(expand_monomial({x, v})));
auto just = multiplication_const_justification{ci, rational(1)}; // multiply var
return add_ineq(just, t, k, rational(0));
SASSERT(p.offset() == 0);
return add_constraint(to_pdd(x) * p, k, multiplication_var_justification{ci, x});
}
void stellensatz::init_occurs() {
m_occurs.reset();
m_occurs.reserve(c().lra_solver().number_of_vars());
for (auto ci : m_solver.lra().constraints().indices())
for (lp::constraint_index ci = 0; ci < m_constraints.size(); ++ci)
init_occurs(ci);
}
void stellensatz::init_occurs(lp::constraint_index ci) {
if (ci == lp::null_ci)
return;
ci = normalize_ineq(ci);
auto const &con = m_solver.lra().constraints()[ci];
for (auto cv : con.lhs()) {
auto v = cv.j();
if (is_mon_var(v)) {
for (auto w : m_mon2vars[v])
m_occurs[w].push_back(ci);
}
else
auto const &con = m_constraints[ci];
for (auto v : con.p.free_vars())
m_occurs[v].push_back(ci);
}
}
// 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.
// otherwise create a fresh variable representing the monomial.
// todo: if _vars is a square, then add bound axiom.
lpvar stellensatz::mk_monomial(svector<lpvar> const& vars) {
lpvar v;
if (vars.empty())
return lp::null_lpvar;
if (vars.size() == 1)
return vars[0];
svector<lpvar> _vars(vars);
std::sort(_vars.begin(), _vars.end());
if (m_vars2mon.find(_vars, v))
return v;
v = add_var(is_int(vars));
m_vars2mon.insert(_vars, v);
m_mon2vars.insert(v, _vars);
SASSERT(m_values.size() == v);
m_values.push_back(cvalue(vars));
return v;
}
lpvar stellensatz::mk_monomial(svector<lp::lpvar> const &_vars, lp::lpvar j) {
svector<lpvar> vars(_vars);
if (is_mon_var(j))
vars.append(m_mon2vars[j]);
else
vars.push_back(j);
return mk_monomial(vars);
}
svector<lp::lpvar> stellensatz::expand_monomial(svector<lp::lpvar> const& vars) {
svector<lp::lpvar> result;
for (auto v : vars) {
if (is_mon_var(v))
result.append(m_mon2vars[v]);
else
result.push_back(v);
}
return result;
}
lpvar stellensatz::mk_term(term_offset& t) {
if (t.second != 0) {
auto w = add_var(t.second.is_int()); // or reuse the constant 1 that is already there.
m_values.push_back(t.second);
m_solver.lra().add_var_bound(w, lp::lconstraint_kind::GE, t.second);
m_solver.lra().add_var_bound(w, lp::lconstraint_kind::LE, t.second);
m_justifications.push_back(nullptr);
m_justifications.push_back(nullptr);
t.first.add_monomial(rational(1), w);
t.second = 0;
}
auto ti = m_solver.lra().add_term(t.first.coeffs_as_vector(), m_solver.lra().number_of_vars());
m_values.push_back(cvalue(t.first));
m_occurs.reserve(m_solver.lra().number_of_vars());
return ti;
}
bool stellensatz::is_int(svector<lp::lpvar> const& vars) const {
@ -622,11 +489,9 @@ namespace nla {
}
bool stellensatz::constraint_is_true(lp::constraint_index ci) const {
auto const& con = m_solver.lra().constraints()[ci];
auto lhs = -con.rhs();
for (auto const& cv : con.lhs())
lhs += cv.coeff() * m_values[cv.j()];
switch (con.kind()) {
auto const& [p, k, j] = m_constraints[ci];
auto lhs = cvalue(p);
switch (k) {
case lp::lconstraint_kind::GT: return lhs > 0;
case lp::lconstraint_kind::GE: return lhs >= 0;
case lp::lconstraint_kind::LE: return lhs <= 0;
@ -644,53 +509,18 @@ namespace nla {
display_product(out, vars);
out << "\n";
}
for (auto ci : m_solver.lra().constraints().indices()) {
for (unsigned ci = 0; ci < m_constraints.size(); ++ci) {
out << "(" << ci << ") ";
display_constraint(out, ci);
bool is_true = constraint_is_true(ci);
out << (is_true ? " [true]" : " [false]") << "\n";
auto j = m_justifications.get(ci);
if (!j)
continue;
out << "\t<- ";
display(out, *j);
display(out, m_constraints[ci].j);
out << "\n";
}
return out;
}
std::ostream& stellensatz::display(std::ostream& out, justification const& b) const {
if (std::holds_alternative<external_justification>(b)) {
auto &dep = std::get<external_justification>(b);
unsigned_vector cs;
c().lra_solver().dep_manager().linearize(dep.dep, cs);
for (auto c : cs)
out << "[o " << c << "] "; // constraint index from c().lra_solver.
}
else if (std::holds_alternative<multiplication_justification>(b)) {
auto &m = std::get<multiplication_justification>(b);
out << "mult (" << m.left << ") (" << m.right << ")";
}
else if (std::holds_alternative<addition_justification>(b)) {
auto &m = std::get<addition_justification>(b);
out << "(" << m.left << ") (" << m.right << ")";
}
else if (std::holds_alternative<gcd_justification>(b)) {
auto &m = std::get<gcd_justification>(b);
out << "gcd(" << m.ci << ")";
}
else if (std::holds_alternative<assumption_justification>(b)) {
out << "assumption";
}
else if (std::holds_alternative<multiplication_const_justification>(b)) {
auto &m = std::get<multiplication_const_justification>(b);
out << m.mul << " * (" << m.ci << ")";
}
else
UNREACHABLE();
return out;
}
std::ostream& stellensatz::display_product(std::ostream& out, svector<lpvar> const& vars) const {
bool first = true;
for (auto v : vars) {
@ -716,32 +546,56 @@ namespace nla {
}
std::ostream& stellensatz::display_var(std::ostream& out, lpvar j) const {
if (is_mon_var(j))
display_product(out, m_mon2vars[j]);
else
if (is_mon_var(j)) {
// display_product(out, c().emons()[j]);
}
else {
out << "j" << j;
}
return out;
}
std::ostream& stellensatz::display_constraint(std::ostream& out, lp::constraint_index ci) const {
auto const& con = m_solver.lra().constraints()[ci];
return display_constraint(out, con.lhs(), con.kind(), con.rhs());
return display_constraint(out, m_constraints[ci]);
}
std::ostream& stellensatz::display_constraint(std::ostream& out,
lp::lar_term const& lhs,
lp::lconstraint_kind k, rational const& rhs) const {
bool first = true;
for (auto cv : lhs) {
auto v = cv.j();
c().display_coeff(out, first, cv.coeff());
first = false;
if (is_mon_var(v))
display_product(out, m_mon2vars[v]);
else
out << "j" << v;
std::ostream& stellensatz::display_constraint(std::ostream& out, constraint const &c) const {
auto const &[p, k, j] = c;
p.display(out);
return out << " " << k << " 0";
}
return out << " " << k << " " << rhs;
std::ostream &stellensatz::display(std::ostream &out, justification const &j) const {
if (std::holds_alternative<external_justification>(j)) {
auto dep = std::get<external_justification>(j).dep;
unsigned_vector cs;
c().lra_solver().dep_manager().linearize(dep, cs);
for (auto c : cs)
out << "[o " << c << "] "; // constraint index from c().lra_solver.
}
else if (std::holds_alternative<multiplication_justification>(j)) {
auto m = std::get<multiplication_justification>(j);
out << "(" << m.left << ") * (" << m.right << ")";
}
else if (std::holds_alternative<addition_justification>(j)) {
auto m = std::get<addition_justification>(j);
out << "(" << m.left << ") + (" << m.right << ")";
}
else if (std::holds_alternative<multiplication_const_justification>(j)) {
auto m = std::get<multiplication_const_justification>(j);
out << m.mul << " * (" << m.ci << ")";
}
else if (std::holds_alternative<assumption_justification>(j)) {
out << "assumption";
}
else if (std::holds_alternative<gcd_justification>(j)) {
auto m = std::get<gcd_justification>(j);
out << " gcd (" << m.ci << ")";
}
else
UNREACHABLE();
return out;
}
std::ostream& stellensatz::display_lemma(std::ostream& out, lp::explanation const& ex,
@ -757,32 +611,31 @@ namespace nla {
m_constraints_in_conflict.insert(ci);
out << "(" << ci << ") ";
display_constraint(out, ci) << " ";
auto j = m_justifications.get(ci);
if (!j)
continue;
display(out, *j) << "\n";
if (std::holds_alternative<multiplication_justification>(*j)) {
auto m = std::get<multiplication_justification>(*j);
auto const& j = m_constraints[ci].j;
display(out, j) << "\n";
if (std::holds_alternative<multiplication_justification>(j)) {
auto m = std::get<multiplication_justification>(j);
todo.push_back(m.left);
todo.push_back(m.right);
}
else if (std::holds_alternative<addition_justification>(*j)) {
auto m = std::get<addition_justification>(*j);
else if (std::holds_alternative<addition_justification>(j)) {
auto m = std::get<addition_justification>(j);
todo.push_back(m.left);
todo.push_back(m.right);
}
else if (std::holds_alternative<multiplication_const_justification>(*j)) {
auto m = std::get<multiplication_const_justification>(*j);
else if (std::holds_alternative<multiplication_const_justification>(j)) {
auto m = std::get<multiplication_const_justification>(j);
todo.push_back(m.ci);
}
else if (std::holds_alternative<external_justification>(*j)) {
else if (std::holds_alternative<external_justification>(j)) {
// do nothing
}
else if (std::holds_alternative<assumption_justification>(*j)) {
else if (std::holds_alternative<assumption_justification>(j)) {
// do nothing
}
else if (std::holds_alternative<gcd_justification>(*j)) {
auto m = std::get<gcd_justification>(*j);
else if (std::holds_alternative<gcd_justification>(j)) {
auto m = std::get<gcd_justification>(j);
todo.push_back(m.ci);
}
else
@ -860,18 +713,22 @@ namespace nla {
auto const &value = values[i];
bool is_int = lra_solver->var_is_int(i);
SASSERT(!is_int || value.is_int());
if (s.is_mon_var(i))
s.m_values[i] = s.mvalue(s.m_mon2vars[i]);
else
s.m_values[i] = value;
if (!s.is_mon_var(i))
continue;
rational val(1);
for (auto w : s.c().emons()[i])
val *= values[w];
values[i] = val;
}
auto indices = lra_solver->constraints().indices();
bool satisfies_products = all_of(indices, [&](auto ci) {
NOT_IMPLEMENTED_YET();
// todo: wrong, needs to be at level of lra constraint evaluation
return s.constraint_is_true(ci);});
if (satisfies_products) {
TRACE(arith, tout << "found satisfying model\n");
for (auto v : s.m_coi.vars())
s.c().lra_solver().set_column_value(v, lp::impq(s.m_values[v], rational(0)));
s.c().lra_solver().set_column_value(v, lp::impq(values[v], rational(0)));
}
return satisfies_products;
}

80
src/math/lp/nla_stellensatz.h
View file

@ -5,9 +5,9 @@
#pragma once
#include "util/scoped_ptr_vector.h"
#include "math/dd/dd_pdd.h"
#include "math/lp/nla_coi.h"
#include "math/lp/int_solver.h"
#include "math/polynomial/polynomial.h"
#include <variant>
namespace nla {
@ -26,7 +26,6 @@ namespace nla {
lbool solve_lia();
bool update_values();
vector<std::pair<lpvar, rational>> m_to_refine;
void saturate_basic_linearize();
public:
solver(stellensatz &s) : s(s) {};
void init();
@ -39,12 +38,11 @@ namespace nla {
solver m_solver;
// factor t into coeff*x^degree*p + q, such that degree(x, q) < degree,
// factor t into x^degree*p + q, such that degree(x, q) < degree,
struct factorization {
rational coeff;
unsigned degree;
lp::lar_term p;
lp::lar_term q;
dd::pdd p, q;
};
struct external_justification {
@ -60,6 +58,10 @@ namespace nla {
lp::constraint_index ci;
rational mul;
};
struct multiplication_var_justification {
lp::constraint_index ci;
lpvar v;
};
struct multiplication_justification {
lp::constraint_index left, right;
};
@ -73,16 +75,32 @@ namespace nla {
gcd_justification,
addition_justification,
multiplication_const_justification,
multiplication_var_justification,
multiplication_justification
>;
using term_offset = std::pair<lp::lar_term, rational>;
struct constraint {
dd::pdd p;
lp::lconstraint_kind k;
justification j;
};
coi m_coi;
scoped_ptr_vector<justification> m_justifications;
void add_justification(lp::constraint_index ci, justification const &j) {
VERIFY(ci == m_justifications.size());
m_justifications.push_back(alloc(justification, j));
}
vector<rational> m_values;
dd::pdd_manager pddm;
vector<constraint> m_constraints;
dd::pdd to_pdd(lpvar v);
lpvar mk_monomial(svector<lpvar> const &vars);
void init_monomial(unsigned mon_var);
term_offset to_term_offset(dd::pdd const &p);
lp::constraint_index add_constraint(dd::pdd p, lp::lconstraint_kind k, justification j);
lp::constraint_index add_var_bound(lp::lpvar v, lp::lconstraint_kind k, rational const &rhs, justification j);
struct eq {
bool operator()(unsigned_vector const& a, unsigned_vector const& b) const {
return a == b;
@ -92,55 +110,30 @@ namespace nla {
u_map<unsigned_vector> m_mon2vars;
bool is_mon_var(lpvar v) const { return m_mon2vars.contains(v); }
unsigned m_max_monomial_degree = 0;
vector<svector<lp::constraint_index>> m_occurs; // map from variable to constraints they occur.
bool is_new_inequality(vector<std::pair<rational, lpvar>> lhs, lp::lconstraint_kind k, rational const &rhs);
// initialization
void init_solver();
void init_vars();
void init_monomial(unsigned mon_var);
void init_occurs();
void init_occurs(lp::constraint_index ci);
void eliminate_variables();
lp::lpvar select_variable_to_eliminate(lp::constraint_index ci);
unsigned degree_of_var_in_constraint(lpvar v, lp::constraint_index ci) const;
unsigned degree_of_var_in_monomial(lpvar v, lpvar mi) const;
factorization factor(lpvar v, lp::constraint_index ci);
lp::lpvar divide(lpvar v, unsigned degree, lpvar mi);
bool constraint_is_true(lp::constraint_index ci) const;
// additional variables and monomials and constraints
using term_offset = std::pair<lp::lar_term, rational>; // term and its offset
svector<lp::lpvar> expand_monomial(svector<lp::lpvar> const & vars);
lpvar mk_monomial(svector<lp::lpvar> const& vars);
lpvar mk_monomial(svector<lp::lpvar> const &vars, lp::lpvar j);
lpvar mk_term(term_offset &t);
void gcd_normalize(vector<std::pair<rational, lpvar>> &t, lp::lconstraint_kind k, rational &rhs);
lp::constraint_index add_ineq(justification const& just, lp::lar_term &t, lp::lconstraint_kind k, rational const &rhs);
lp::constraint_index add_ineq(justification const &just, lpvar j, lp::lconstraint_kind k,
rational const &rhs);
lp::constraint_index assume(lp::lar_term &t, lp::lconstraint_kind k, rational const &rhs);
lp::constraint_index normalize_ineq(lp::constraint_index ci);
lp::constraint_index gcd_normalize(lp::constraint_index ci);
lp::constraint_index assume(dd::pdd const& p, lp::lconstraint_kind k);
lp::constraint_index add(lp::constraint_index left, lp::constraint_index right);
lp::constraint_index multiply(lp::constraint_index ci, rational const &mul);
lp::constraint_index multiply(lp::constraint_index left, lp::constraint_index right);
lp::constraint_index multiply_var(lp::constraint_index ci, lpvar x);
bool is_int(svector<lp::lpvar> const& vars) const;
rational cvalue(lp::lar_term const &t) const;
rational mvalue(lp::lar_term const &t) const;
rational cvalue(svector<lpvar> const &prod) const;
rational mvalue(svector<lpvar> const &prod) const;
rational cvalue(dd::pdd const& p) const;
lpvar add_var(bool is_int);
// lemmas
@ -162,12 +155,11 @@ namespace nla {
std::ostream& display(std::ostream& out) const;
std::ostream& display_product(std::ostream& out, svector<lpvar> const& vars) const;
std::ostream& display_constraint(std::ostream& out, lp::constraint_index ci) const;
std::ostream& display_constraint(std::ostream& out, lp::lar_term const& lhs,
lp::lconstraint_kind k, rational const& rhs) const;
std::ostream& display(std::ostream &out, term_offset const &t) const;
std::ostream& display(std::ostream &out, justification const &bounds) const;
std::ostream& display_constraint(std::ostream& out, constraint const& c) const;
std::ostream &display(std::ostream &out, justification const &j) const;
std::ostream &display_var(std::ostream &out, lpvar j) const;
std::ostream &display_lemma(std::ostream &out, lp::explanation const &ex, vector<ineq> const &ineqs);
std::ostream &display(std::ostream &out, term_offset const &t) const;
public:
stellensatz(core* core);