mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-11-22 05:36:41 +00:00
Code Activity

v2 running

Browse source 
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2025-11-13 13:59:12 -08:00
parent c7084e9998
commit 5ea25dcf60
2 changed files with 282 additions and 189 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -47,32 +47,70 @@ namespace nla {
return k1;
}
lpvar stellensatz::monomial_factory::mk_monomial(lp::lar_solver &lra, svector<lpvar> const &vars) {
lpvar v = lp::null_lpvar;
if (vars.empty())
return 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;
auto is_int = all_of(vars, [&](lpvar v) { return lra.var_is_int(v); });
auto nv = lra.number_of_vars();
v = lra.add_var(nv, is_int);
m_vars2mon.insert(_vars, v);
m_mon2vars.insert(v, _vars);
return v;
}
lpvar stellensatz::monomial_factory::get_monomial(svector<lpvar> const &vars) {
lpvar v = lp::null_lpvar;
if (vars.empty())
return 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;
NOT_IMPLEMENTED_YET();
return lp::null_lpvar;
}
void stellensatz::monomial_factory::init(lpvar v, svector<lpvar> const &_vars) {
svector<lpvar> vars(_vars);
std::sort(vars.begin(), vars.end());
m_vars2mon.insert(vars, v);
m_mon2vars.insert(v, vars);
}
stellensatz::stellensatz(core* core) :
common(core),
m_solver(*this),
m_coi(*core),
pddm(m_solver.lra().number_of_vars())
m_coi(*core),
pddm(core->lra_solver().number_of_vars())
{}
lbool stellensatz::saturate() {
init_solver();
TRACE(arith, display(tout << "stellensatz before saturation\n"));
eliminate_variables();
auto r = eliminate_variables();
if (r == l_false)
return r;
// TODO: populate solver
TRACE(arith, display(tout << "stellensatz after saturation\n"));
return l_undef;
lbool r = m_solver.solve();
r = m_solver.solve();
// IF_VERBOSE(0, verbose_stream() << "stellensatz " << r << "\n");
if (r == l_false)
add_lemma();
return r;
}
void stellensatz::init_solver() {
m_solver.init();
m_vars2mon.reset();
m_mon2vars.reset();
m_constraints.reset();
m_monomial_factory.reset();
m_coi.init();
init_vars();
init_occurs();
@ -82,7 +120,9 @@ namespace nla {
auto const& src = c().lra_solver();
auto sz = src.number_of_vars();
for (unsigned v = 0; v < sz; ++v) {
if (m_coi.mons().contains(v))
if (!m_coi.vars().contains(v))
continue;
if (c().is_monic_var(v))
init_monomial(v);
// assert bounds on v in the new solver.
if (src.column_has_lower_bound(v)) {
@ -119,7 +159,7 @@ namespace nla {
dd::pdd sum(pddm);
sum = 0;
for (auto cv : t)
sum *= to_pdd(cv.j())*cv.coeff();
sum += to_pdd(cv.j())*cv.coeff();
return sum;
}
return pddm.mk_var(v);
@ -131,33 +171,14 @@ namespace nla {
if (vars.empty())
to.second += coeff;
else
to.first.add_monomial(coeff, mk_monomial(vars));
to.first.add_monomial(coeff, m_monomial_factory.get_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) {
auto &mon = c().emons()[mon_var];
svector<lpvar> vars(mon.vars());
std::sort(vars.begin(), vars.end());
m_vars2mon.insert(vars, mon_var);
m_mon2vars.insert(mon_var, vars);
m_monomial_factory.init(mon_var, mon.vars());
}
lp::constraint_index stellensatz::add_var_bound(lp::lpvar v, lp::lconstraint_kind k, rational const& rhs, justification j) {
@ -179,13 +200,17 @@ namespace nla {
k = lp::lconstraint_kind::GT;
p = -p;
}
if (k == lp::lconstraint_kind::GT && is_int(p)) {
k = lp::lconstraint_kind::GE;
p += rational(1);
}
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() {
lbool stellensatz::eliminate_variables() {
vector<std::pair<lp::constraint_index, uint_set>> active;
for (unsigned ci = 0; ci < m_constraints.size(); ++ci) {
if (!constraint_is_true(ci))
@ -200,17 +225,24 @@ namespace nla {
auto [ci, tabu] = active.back();
active.pop_back();
auto x = select_variable_to_eliminate(ci);
if (x == lp::null_lpvar)
continue;
auto f = factor(x, ci);
TRACE(arith, tout << "factor " << m_constraints[ci].p << " -> j" << x << "^" << f.degree << " * " << f.p << " + "
<< f.q << "\n");
auto p_value = cvalue(f.p);
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);
if (!f.q.is_zero()) {
// ci & -p_is_0*x^f.degree => new_ci
p_is_0 = multiply(p_is_0, rational(-1));
dd::pdd r = pddm.mk_val(rational(-1));
for (unsigned i = 0; i < f.degree; ++i)
p_is_0 = multiply_var(p_is_0, x);
r = r * pddm.mk_var(x);
p_is_0 = multiply(p_is_0, r);
auto new_ci = add(ci, p_is_0);
TRACE(arith,
display_constraint(tout << "reduced", new_ci) << "\n");
init_occurs(new_ci);
uint_set new_tabu(tabu);
uint_set q_vars;
@ -223,7 +255,13 @@ namespace nla {
}
continue;
}
for (auto other_ci : m_occurs[x]) {
unsigned num_x = m_occurs[x].size();
for (unsigned i = 0; i < f.degree; ++i)
f.p *= to_pdd(x);
auto [_m1, _f_p] = f.p.var_factors();
for (unsigned cx = 0; cx < num_x; ++cx) {
auto other_ci = m_occurs[x][cx];
if (other_ci == ci)
continue;
auto const &other_ineq = m_constraints[other_ci];
@ -236,6 +274,8 @@ namespace nla {
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_value2 == 0)
continue;
if (p_value > 0 && p_value2 > 0)
continue;
if (p_value < 0 && p_value2 < 0)
@ -243,18 +283,24 @@ namespace nla {
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);
TRACE(arith, tout << "factor " << other_ineq.p << " -> j" << x << "^" << g.degree << " * " << g.p
<< " + " << g.q << "\n");
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;
auto m1(_m1);
auto f_p(_f_p);
dd::pdd::merge(m1, m2, m1m2);
SASSERT(m1m2.contains(x));
f_p = f_p.mul(m1);
g_p = g_p.mul(m2);
TRACE(arith, tout << "m1 " << m1 << " m2 " << m2 << " m1m2: " << m1m2 << "\n");
auto sign_f = cvalue(f_p) < 0;
SASSERT(sign_f != cvalue(g_p) < 0);
@ -264,38 +310,58 @@ namespace nla {
// m1m2 * f_p + f.q >= 0
// m1m2 * g_p + g.q >= 0
//
auto ci_a = ci;
auto ci_b = other_ci;
lp::constraint_index ci_a, ci_b;
auto aj = assumption_justification();
if (f_p.is_val())
ci_a = multiply(other_ci, f_p.val());
ci_b = multiply(other_ci, pddm.mk_val(f_p.val()));
else if (sign_f)
ci_a = multiply(other_ci, add_constraint(f_p, lp::lconstraint_kind::LT, aj));
ci_b = 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));
ci_b = multiply(other_ci, add_constraint(f_p, lp::lconstraint_kind::GT, aj));
if (g_p.is_val())
ci_b = multiply(ci, g_p.val());
ci_a = multiply(ci, pddm.mk_val(g_p.val()));
else if (!sign_f)
ci_b = multiply(ci, add_constraint(g_p, lp::lconstraint_kind::LT, aj));
ci_a = 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));
ci_a = multiply(ci, add_constraint(g_p, lp::lconstraint_kind::GT, aj));
auto new_ci = add(ci_a, ci_b);
init_occurs(new_ci);
if (m_constraints[new_ci].p.degree() < 3)
init_occurs(new_ci);
TRACE(arith, tout << "eliminate j" << x << ":\n";
display_constraint(tout << "ci: ", ci) << "\n";
display_constraint(tout << "other_ci: ", other_ci) << "\n";
display_constraint(tout << "ci_a: ", ci_a) << "\n";
display_constraint(tout << "ci_b: ", ci_b) << "\n";
display_constraint(tout << "new_ci: ", new_ci) << "\n");
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});
if (new_p.is_val()) {
lp::explanation ex;
lemma_builder new_lemma(c(), "stellensatz");
m_constraints_in_conflict.reset();
explain_constraint(new_lemma, new_ci, ex);
new_lemma &= ex;
IF_VERBOSE(2, verbose_stream() << "stellensatz unsat \n" << new_lemma << "\n");
TRACE(arith, tout << "unsat\n" << new_lemma << "\n");
c().lra_solver().settings().stats().m_nla_stellensatz++;
return l_false;
}
if (m_constraints[new_ci].p.degree() < 3) {
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});
}
}
}
}
return l_undef;
}
lp::lpvar stellensatz::select_variable_to_eliminate(lp::constraint_index ci) {
@ -319,29 +385,7 @@ namespace nla {
return {degree, lc, rest};
}
//
// convert a conflict from m_solver.lra()/lia() into
// a conflict for c().lra_solver()
// 1. constraints that are obtained by multiplication are explained from the original constraint
// 2. bounds justifications are added as justifications to the lemma.
//
void stellensatz::add_lemma() {
lp::explanation const &ex1 = m_solver.ex();
vector<ineq> const &ineqs = m_solver.ineqs();
TRACE(arith, display(tout); display_lemma(tout, ex1, ineqs));
auto& lra = c().lra_solver();
lp::explanation ex2;
lemma_builder new_lemma(c(), "stellensatz");
m_constraints_in_conflict.reset();
for (auto p : ex1)
explain_constraint(new_lemma, p.ci(), ex2);
new_lemma &= ex2;
for (auto const &ineq : ineqs)
new_lemma |= ineq;
IF_VERBOSE(2, verbose_stream() << "stellensatz unsat \n" << new_lemma << "\n");
TRACE(arith, tout << "unsat\n" << new_lemma << "\n");
c().lra_solver().settings().stats().m_nla_stellensatz++;
}
//
// a constraint can be explained by a set of bounds used as justifications
@ -354,7 +398,7 @@ namespace nla {
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);
c().lra_solver().push_explanation(dep.dep, ex);
}
else if (std::holds_alternative<multiplication_justification>(b)) {
auto& m = std::get<multiplication_justification>(b);
@ -366,8 +410,8 @@ namespace nla {
explain_constraint(new_lemma, m.left, ex);
explain_constraint(new_lemma, m.right, ex);
}
else if (std::holds_alternative<multiplication_const_justification>(b)) {
auto& m = std::get<multiplication_const_justification>(b);
else if (std::holds_alternative<multiplication_poly_justification>(b)) {
auto& m = std::get<multiplication_poly_justification>(b);
explain_constraint(new_lemma, m.ci, ex);
}
else if (std::holds_alternative<assumption_justification>(b)) {
@ -384,16 +428,16 @@ namespace nla {
rational stellensatz::cvalue(dd::pdd const& p) const {
dd::pdd_eval eval;
// eval.var2val() = [&](unsigned v) -> void { return c().val(v); };
eval.var2val() = [&](unsigned v) -> rational { return c().val(v); };
return eval(p);
}
lp::constraint_index stellensatz::gcd_normalize(lp::constraint_index ci) {
auto [p, k, j] = m_constraints[ci];
rational g(0);
bool is_int = all_of(p.free_vars(), [&](unsigned v) { return c().lra_solver().var_is_int(v); });
bool _is_int = is_int(p);
for (auto const& [c, is_constant] : p.coefficients())
if (!is_constant || !is_int)
if (!is_constant || !_is_int)
g = gcd(g, c);
if (g == 0 || g == 1)
return ci;
@ -405,7 +449,7 @@ namespace nla {
}
case lp::lconstraint_kind::GT: {
auto q = p;
if (is_int) {
if (_is_int) {
q += rational(1);
k = lp::lconstraint_kind::GE;
}
@ -439,11 +483,15 @@ namespace nla {
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);
// p >= 0 => a * p >= 0 if a > 0,
// p = 0 => p * q = 0 no matter what q
lp::constraint_index stellensatz::multiply(lp::constraint_index ci, dd::pdd q) {
auto const& [p, k, j] = m_constraints[ci];
return add_constraint(p * mul, mul > 0 ? k : swap_side(k), multiplication_const_justification{ci, mul});
auto k1 = k;
if (q.is_val() && q.val() < 0)
k1 = swap_side(k1);
SASSERT(q.is_val() || k1 == lp::lconstraint_kind::EQ);
return add_constraint(p * q, k1, multiplication_poly_justification{ci, q});
}
lp::constraint_index stellensatz::multiply(lp::constraint_index left, lp::constraint_index right) {
@ -453,13 +501,6 @@ namespace nla {
return add_constraint(p*q, k, multiplication_justification{left, right});
}
lp::constraint_index stellensatz::multiply_var(lp::constraint_index ci, lpvar x) {
auto const& [p, k, j] = m_constraints[ci];
SASSERT(k == lp::lconstriant_kind::EQ);
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());
@ -477,15 +518,11 @@ namespace nla {
}
bool stellensatz::is_int(svector<lp::lpvar> const& vars) const {
return all_of(vars, [&](lpvar v) { return m_solver.lra().var_is_int(v); });
return all_of(vars, [&](lpvar v) { return c().lra_solver().var_is_int(v); });
}
lpvar stellensatz::add_var(bool is_int) {
auto v = m_solver.lra().number_of_vars();
auto w = m_solver.lra().add_var(v, is_int);
SASSERT(v == w);
m_occurs.reserve(m_solver.lra().number_of_vars());
return w;
bool stellensatz::is_int(dd::pdd const &p) const {
return is_int(p.free_vars());
}
bool stellensatz::constraint_is_true(lp::constraint_index ci) const {
@ -503,12 +540,14 @@ namespace nla {
}
std::ostream& stellensatz::display(std::ostream& out) const {
m_solver.lra().display(out);
#if 0
// m_solver.lra().display(out);
for (auto const& [vars, v] : m_vars2mon) {
out << "j" << v << " := ";
display_product(out, vars);
out << "\n";
}
#endif
for (unsigned ci = 0; ci < m_constraints.size(); ++ci) {
out << "(" << ci << ") ";
display_constraint(out, ci);
@ -546,14 +585,10 @@ namespace nla {
}
std::ostream& stellensatz::display_var(std::ostream& out, lpvar j) const {
if (is_mon_var(j)) {
// display_product(out, c().emons()[j]);
}
else {
out << "j" << j;
}
return out;
if (c().is_monic_var(j))
return display_product(out, c().emon(j).vars());
else
return out << "j" << j;
}
std::ostream& stellensatz::display_constraint(std::ostream& out, lp::constraint_index ci) const {
@ -582,9 +617,9 @@ namespace nla {
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<multiplication_poly_justification>(j)) {
auto m = std::get<multiplication_poly_justification>(j);
out << m.p << " * (" << m.ci << ")";
}
else if (std::holds_alternative<assumption_justification>(j)) {
out << "assumption";
@ -598,8 +633,7 @@ namespace nla {
return out;
}
std::ostream& stellensatz::display_lemma(std::ostream& out, lp::explanation const& ex,
vector<ineq> const& ineqs) {
std::ostream &stellensatz::display_lemma(std::ostream &out, lp::explanation const &ex) {
m_constraints_in_conflict.reset();
svector<lp::constraint_index> todo;
for (auto c : ex)
@ -624,14 +658,14 @@ namespace nla {
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_poly_justification>(j)) {
auto m = std::get<multiplication_poly_justification>(j);
todo.push_back(m.ci);
}
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)) {
@ -642,10 +676,6 @@ namespace nla {
NOT_IMPLEMENTED_YET();
}
for (auto ineq : ineqs) {
term_offset t(ineq.term(), rational(0));
display(out, t) << " " << ineq.cmp() << " " << ineq.rs() << "\n";
}
return out;
}
@ -659,14 +689,63 @@ namespace nla {
// option: detect squares and add axioms for violated squares
// option: add NIA filters (non-linear divisbility axioms)
void stellensatz::solver::init() {
lra_solver = alloc(lp::lar_solver);
int_solver = alloc(lp::int_solver, *lra_solver);
m_ex.clear();
m_ineqs.reset();
m_internal2external_constraints.reset();
m_monomial_factory.reset();
auto &src = s.c().lra_solver();
auto &dst = *lra_solver;
for (unsigned v = 0; v < src.number_of_vars(); ++v)
dst.add_var(v, src.var_is_int(v));
for (lp::constraint_index ci = 0; ci < s.m_constraints.size(); ++ci) {
auto const &[p, k, j] = s.m_constraints[ci];
auto [t, offset] = to_term_offset(p);
auto coeffs = t.coeffs_as_vector();
if (coeffs.empty())
continue;
SASSERT(!coeffs.empty());
auto ti = dst.add_term(coeffs, dst.number_of_vars());
auto new_ci = dst.add_var_bound(ti, k, -offset);
m_internal2external_constraints.setx(new_ci, ci, ci);
}
}
stellensatz::term_offset stellensatz::solver::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, m_monomial_factory.mk_monomial(*lra_solver, vars));
}
return to;
}
//
// convert a conflict from m_solver.lra()/lia() into
// a conflict for c().lra_solver()
//
void stellensatz::solver::add_lemma() {
TRACE(arith, s.display(tout); s.display_lemma(tout, m_ex));
auto &src = s.c().lra_solver();
lp::explanation ex2;
lemma_builder new_lemma(s.c(), "stellensatz");
s.m_constraints_in_conflict.reset();
for (auto p : m_ex)
s.explain_constraint(new_lemma, m_internal2external_constraints[p.ci()], ex2);
new_lemma &= ex2;
IF_VERBOSE(2, verbose_stream() << "stellensatz unsat \n" << new_lemma << "\n");
TRACE(arith, tout << "unsat\n" << new_lemma << "\n");
s.c().lra_solver().settings().stats().m_nla_stellensatz++;
}
lbool stellensatz::solver::solve() {
init();
lbool r = solve_lra();
// IF_VERBOSE(0, verbose_stream() << "solve lra " << r << "\n");
if (r != l_true)
@ -675,6 +754,7 @@ namespace nla {
// IF_VERBOSE(0, verbose_stream() << "solve lia " << r << "\n");
if (r != l_true)
return r;
return l_undef;
if (update_values())
return l_true;
@ -687,6 +767,7 @@ namespace nla {
return l_true;
if (status == lp::lp_status::INFEASIBLE) {
lra_solver->get_infeasibility_explanation(m_ex);
add_lemma();
return l_false;
}
return l_undef;
@ -694,8 +775,11 @@ namespace nla {
lbool stellensatz::solver::solve_lia() {
switch (int_solver->check(&m_ex)) {
case lp::lia_move::sat: return l_true;
case lp::lia_move::conflict: return l_false;
case lp::lia_move::sat:
return l_true;
case lp::lia_move::conflict:
add_lemma();
return l_false;
default: // TODO: an option is to perform (bounded) search here to get an LIA verdict.
return l_undef;
}
@ -706,6 +790,8 @@ namespace nla {
// return true if the new assignment satisfies the products.
// return false if value constraints are not satisfied on monomials and there is a false constaint.
bool stellensatz::solver::update_values() {
return false;
#if 0
std::unordered_map<lpvar, rational> values;
lra_solver->get_model(values);
unsigned sz = lra_solver->number_of_vars();
@ -731,5 +817,6 @@ namespace nla {
s.c().lra_solver().set_column_value(v, lp::impq(values[v], rational(0)));
}
return satisfies_products;
#endif
}
}

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

@ -16,51 +16,17 @@ namespace nla {
class lar_solver;
class stellensatz : common {
class solver {
stellensatz &s;
scoped_ptr<lp::lar_solver> lra_solver;
scoped_ptr<lp::int_solver> int_solver;
lp::explanation m_ex;
vector<ineq> m_ineqs;
lbool solve_lra();
lbool solve_lia();
bool update_values();
vector<std::pair<lpvar, rational>> m_to_refine;
public:
solver(stellensatz &s) : s(s) {};
void init();
lbool solve();
lp::lar_solver &lra() { return *lra_solver; }
lp::lar_solver const &lra() const { return *lra_solver; }
lp::explanation &ex() { return m_ex; }
vector<ineq> &ineqs() { return m_ineqs; }
};
solver m_solver;
// factor t into x^degree*p + q, such that degree(x, q) < degree,
struct factorization {
unsigned degree;
dd::pdd p, q;
};
struct external_justification {
u_dependency *dep = nullptr;
external_justification(u_dependency *d) : dep(d) {}
};
struct assumption_justification {
};
struct assumption_justification {};
struct addition_justification {
lp::constraint_index left, right;
};
struct multiplication_const_justification {
struct multiplication_poly_justification {
lp::constraint_index ci;
rational mul;
};
struct multiplication_var_justification {
lp::constraint_index ci;
lpvar v;
dd::pdd p;
};
struct multiplication_justification {
lp::constraint_index left, right;
@ -72,14 +38,10 @@ namespace nla {
using justification = std::variant<
external_justification,
assumption_justification,
gcd_justification,
addition_justification,
multiplication_const_justification,
multiplication_var_justification,
multiplication_justification
>;
using term_offset = std::pair<lp::lar_term, rational>;
gcd_justification,
addition_justification,
multiplication_poly_justification,
multiplication_justification>;
struct constraint {
dd::pdd p;
@ -87,28 +49,74 @@ namespace nla {
justification j;
};
class monomial_factory {
struct eq {
bool operator()(unsigned_vector const &a, unsigned_vector const &b) const {
return a == b;
}
};
map<unsigned_vector, unsigned, svector_hash<unsigned_hash>, eq> m_vars2mon;
u_map<unsigned_vector> m_mon2vars;
bool is_mon_var(lpvar v) const {
return m_mon2vars.contains(v);
}
public:
void reset() {
m_vars2mon.reset();
m_mon2vars.reset();
}
lpvar mk_monomial(lp::lar_solver& lra, svector<lpvar> const &vars);
lpvar get_monomial(svector<lpvar> const &vars);
void init(lpvar v, svector<lpvar> const &_vars);
};
using term_offset = std::pair<lp::lar_term, rational>;
class solver {
stellensatz &s;
scoped_ptr<lp::lar_solver> lra_solver;
scoped_ptr<lp::int_solver> int_solver;
lp::explanation m_ex;
unsigned_vector m_internal2external_constraints;
monomial_factory m_monomial_factory;
lbool solve_lra();
lbool solve_lia();
bool update_values();
void init();
void add_lemma();
term_offset to_term_offset(dd::pdd const &p);
public:
solver(stellensatz &s) : s(s) {};
lbool solve();
};
solver m_solver;
// factor t into x^degree*p + q, such that degree(x, q) < degree,
struct factorization {
unsigned degree;
dd::pdd p, q;
};
coi m_coi;
dd::pdd_manager pddm;
vector<constraint> m_constraints;
monomial_factory m_monomial_factory;
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;
}
};
map<unsigned_vector, unsigned, svector_hash<unsigned_hash>, eq> m_vars2mon;
u_map<unsigned_vector> m_mon2vars;
bool is_mon_var(lpvar v) const { return m_mon2vars.contains(v); }
vector<svector<lp::constraint_index>> m_occurs; // map from variable to constraints they occur.
@ -118,7 +126,7 @@ namespace nla {
void init_occurs();
void init_occurs(lp::constraint_index ci);
void eliminate_variables();
lbool 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;
factorization factor(lpvar v, lp::constraint_index ci);
@ -128,16 +136,14 @@ namespace nla {
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 ci, dd::pdd p);
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;
bool is_int(dd::pdd const &p) const;
rational cvalue(dd::pdd const& p) const;
lpvar add_var(bool is_int);
// lemmas
void add_lemma();
indexed_uint_set m_constraints_in_conflict;
void explain_constraint(lemma_builder& new_lemma, lp::constraint_index ci, lp::explanation &ex);
@ -158,7 +164,7 @@ namespace nla {
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_lemma(std::ostream &out, lp::explanation const &ex);
std::ostream &display(std::ostream &out, term_offset const &t) const;
public: