diff --git a/.clang-format b/.clang-format index 6bfa48374..f8fd218e4 100644 --- a/.clang-format +++ b/.clang-format @@ -28,7 +28,7 @@ AccessModifierOffset: -4 # Function definitions AlwaysBreakAfterReturnType: None -AllowShortFunctionsOnASingleLine: Empty +AllowShortFunctionsOnASingleLine: true AllowShortIfStatementsOnASingleLine: false AllowShortLoopsOnASingleLine: false # Ensure function-opening brace is attached to the signature diff --git a/src/math/lp/nla_stellensatz.cpp b/src/math/lp/nla_stellensatz.cpp index c2f399485..a7620344b 100644 --- a/src/math/lp/nla_stellensatz.cpp +++ b/src/math/lp/nla_stellensatz.cpp @@ -117,36 +117,35 @@ namespace nla { common(core), m_solver(*this), m_coi(*core), - pddm(core->lra_solver().number_of_vars()) { + pddm(core->lra_solver().number_of_vars()), + m_di(m_dm, core->reslim()) { } lbool stellensatz::saturate() { + // TODO: set up initial bounds unsigned num_conflicts = 0, num_constraints = 0; init_solver(); TRACE(arith, display(tout << "stellensatz before saturation\n")); start_saturate: - m_last_constraint = lp::null_ci; num_constraints = m_constraints.size(); - lbool r = model_repair(); - - if (num_constraints == m_constraints.size()) - return l_undef; - + lbool r = search(); + if (r == l_undef) r = m_solver.solve(m_core); TRACE(arith, display(tout << "stellensatz after saturation " << r << "\n")); if (r == l_false) { - ++num_conflicts; - IF_VERBOSE(2, verbose_stream() << "(nla.stellensatz :conflicts " << num_conflicts << ":constraints " << m_constraints.size() << ")\n"); + ++num_conflicts; + IF_VERBOSE(2, verbose_stream() << "(nla.stellensatz :conflicts " << num_conflicts << ":constraints " + << m_constraints.size() << ")\n"); if (num_conflicts >= m_config.max_conflicts) return l_undef; while (backtrack(m_core)) { - ++c().lp_settings().stats().m_stellensatz.m_num_backtracks; + ++c().lp_settings().stats().m_stellensatz.m_num_backtracks; lbool rb = m_solver.solve(m_core); if (rb == l_false) continue; - if (rb == l_undef) - return rb; + if (rb == l_undef) + return rb; m_solver.update_values(m_values); goto start_saturate; } @@ -158,6 +157,7 @@ namespace nla { return r; } + void stellensatz::pop_constraint() { auto const &[p, k, j] = m_constraints.back(); m_constraint_index.erase({p.index(), k}); @@ -165,12 +165,14 @@ namespace nla { } void stellensatz::init_solver() { - m_trail.reset(); + m_ctrail.reset(); + m_vtrail.reset(); m_monomial_factory.reset(); m_coi.init(); m_values.reset(); init_vars(); init_occurs(); + init_bounds(); } void stellensatz::init_vars() { @@ -203,6 +205,42 @@ namespace nla { m_occurs.reserve(sz); } + void stellensatz::init_bounds() { + for (unsigned ci = 0; ci < m_constraints.size(); ++ci) { + auto [p, k, j] = m_constraints[ci]; + if (!p.is_unilinear()) + continue; + // ax + b >= 0 + auto b = p.lo().val(); + auto a = p.hi().val(); + auto x = p.var(); + if (a > 0) { + auto lb = -b / a; + if (var_is_int(x)) + lb = floor(lb); + if (!has_lo(x) || + lo_val(x) < lb || (!lo_is_strict(x) && k == lp::lconstraint_kind::GT && lo_val(x) == lb)) { + bound_info bi(x, k, lb, m_level, m_lower[x]); + m_lower[x] = m_bounds.size(); + m_bounds.push_back(bi); + } + } + else if (a < 0) { + auto ub = -b / a; + if (var_is_int(x)) + ub = ceil(ub); + if (!has_hi(x) || hi_val(x) > ub || + (!hi_is_strict(x) && k == lp::lconstraint_kind::GT && hi_val(x) == ub)) { + bound_info bi(x, + k == lp::lconstraint_kind::GT ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE, + ub, m_level, m_upper[x]); + m_upper[x] = m_bounds.size(); + m_bounds.push_back(bi); + } + } + } + } + // set the model based on m_values computed by the solver bool stellensatz::set_model() { if (any_of(m_constraints, [&](auto const &c) { return !constraint_is_true(c); })) @@ -318,54 +356,45 @@ namespace nla { s.pop_constraint(); } }; - m_trail.push_scope(); - m_trail.push(undo_constraint(*this)); + m_ctrail.push_scope(); + m_ctrail.push(undo_constraint(*this)); return ci; } - void stellensatz::add_active(lp::constraint_index ci) { - if (m_active.contains(ci)) - return; - if (m_constraints[ci].p.degree() > m_config.max_degree) - return; - m_active.insert(ci); - } - - lp::constraint_index stellensatz::resolve_variable(lpvar x, lp::constraint_index ci, lp::constraint_index other_ci, - rational const &p_value, factorization const &f, unsigned_vector const &_m1, - dd::pdd _f_p) { + lp::constraint_index stellensatz::resolve_variable(lpvar x, lp::constraint_index ci, + lp::constraint_index other_ci) { TRACE(arith, tout << "resolve j" << x << " between ci " << ci << " and other_ci " << other_ci << "\n"); - if (ci == other_ci) + if (ci == lp::null_ci || other_ci == lp::null_ci) return lp::null_ci; + if (constraint_is_true(ci) && constraint_is_true(other_ci)) + return lp::null_ci; + auto f = factor(x, ci); if (f.degree != 1) - return lp::null_ci; // future - could handle this - auto const &[other_p, other_k, other_j] = m_constraints[other_ci]; - auto const &[p, k, j] = m_constraints[ci]; + return lp::null_ci; // future - could handle this auto g = factor(x, other_ci); if (g.degree != 1) - return lp::null_ci; // not handling degree > 1 + return lp::null_ci; // not handling degree > 1 + auto p_value1 = value(f.p); auto p_value2 = value(g.p); // - // check that p_value and p_value2 have different signs + // check that p_value1 and p_value2 have different signs // check that other_p is disjoint from tabu // compute overlaps extending x // - SASSERT(g.degree > 0); - SASSERT(p_value != 0); - if (g.degree > f.degree) // future: could handle this too by considering tabu to be a map into degrees. + SASSERT(p_value1 != 0); + if (p_value1 > 0 && p_value2 > 0) return lp::null_ci; - if (p_value > 0 && p_value2 > 0) - return lp::null_ci; - if (p_value < 0 && p_value2 < 0) + if (p_value1 < 0 && p_value2 < 0) return lp::null_ci; + 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; - auto m1(_m1); - auto f_p(_f_p); dd::pdd::merge(m1, m2, m1m2); SASSERT(m1m2.contains(x)); f_p = f_p.mul(m1); @@ -380,6 +409,8 @@ namespace nla { // m1m2 * g_p + g.q >= 0 // lp::constraint_index ci_a, ci_b; + auto const &[other_p, other_k, other_j] = m_constraints[other_ci]; + auto const &[p, k, j] = m_constraints[ci]; bool g_strict = other_k == lp::lconstraint_kind::GT; bool f_strict = k == lp::lconstraint_kind::GT; @@ -399,26 +430,20 @@ namespace nla { auto new_ci = add(ci_a, ci_b); + move_up(x); ++c().lp_settings().stats().m_stellensatz.m_num_resolutions; - - TRACE(arith, tout << "eliminate j" << x << ":\n"; - display_constraint(tout, ci) << "\n"; - display_constraint(tout, other_ci) << "\n"; - display_constraint(tout, ci_a) << "\n"; - display_constraint(tout, ci_b) << "\n"; - display_constraint(tout, new_ci) << "\n"); - - if (constraint_is_conflict(new_ci)) - return new_ci; + TRACE(arith, tout << "eliminate j" << x << ":\n"; display_constraint(tout, ci) << "\n"; + display_constraint(tout, other_ci) << "\n"; display_constraint(tout, ci_a) << "\n"; + display_constraint(tout, ci_b) << "\n"; display_constraint(tout, new_ci) << "\n"); new_ci = factor(new_ci); - init_occurs(new_ci); + init_occurs(new_ci); - CTRACE(arith, !is_new_constraint(new_ci), display_constraint(tout << "not new ", new_ci) << "\n"); - CTRACE(arith, is_new_constraint(new_ci), display_constraint(tout << "new ", new_ci) << "\n"); - return new_ci; + CTRACE(arith, new_ci == m_constraints.size() - 1, display_constraint(tout << "not new ", new_ci) << "\n"); + CTRACE(arith, new_ci != m_constraints.size() - 1, display_constraint(tout << "new ", new_ci) << "\n"); + return new_ci; } bool stellensatz::conflict(lp::constraint_index ci) { @@ -441,180 +466,6 @@ namespace nla { c().lra_solver().settings().stats().m_nla_stellensatz++; } - lbool stellensatz::model_repair() { - m_active.reset(); - m_processed.reset(); - for (lp::constraint_index ci = 0; ci < m_constraints.size(); ++ci) - m_active.insert(ci); - m_level2var.reset(); - m_var2level.reset(); - for (unsigned v = 0; v < m_values.size(); ++v) - m_level2var.push_back(v); - random_gen rand(c().random()); - shuffle(m_level2var.size(), m_level2var.begin(), rand); - m_var2level.resize(m_values.size(), m_level2var.size()); - for (unsigned i = 0; i < m_level2var.size(); ++i) - m_var2level[m_level2var[i]] = i; - return model_repair_loop(); - } - - lbool stellensatz::model_repair_loop() { - unsigned level = 0; - while (c().reslim().inc()) { - if (level >= m_level2var.size()) - break; - m_last_constraint = lp::null_ci; - auto ci = model_repair(m_level2var[level]); - TRACE(arith, - display_constraint(tout << "after repairing j" << m_level2var[level] << " @ " << level << " ", ci) - << "\nconflict " - << conflict(ci) << "\n"); - if (ci == lp::null_ci) - ++level; - else if (conflict(ci)) - return l_false; - else if (m_processed.contains(ci) || constraint_is_trivial(ci)) - ++level; - else - level = max_level(m_constraints[ci]); - } - return all_of(m_constraints, [&](constraint const &c) { return constraint_is_true(c); }) ? l_true : l_undef; - } - - unsigned stellensatz::max_level(constraint const &c) const { - unsigned level = 0; - auto const &vars = c.p.free_vars(); - for (auto v : vars) - if (true || c.p.degree(v) == 1) - level = std::max(level, m_var2level[v]); - // display_constraint(verbose_stream(), c) << " " << level << "\n"; - return level; - } - - lp::constraint_index stellensatz::model_repair(lp::lpvar v) { - auto [lo, hi, inf, sup, vanishing] = find_bounds(v); - TRACE(arith, tout << "bounds for v" << v << " @ " << m_var2level[v] << " : " << lo << " to " << hi << "\n"; - display_constraint(tout << " inf ", inf) << "\n"; - display_constraint(tout << " sup ", sup) << "\n"; - display_constraint(tout << " vanishing ", vanishing) << "\n";); - if (vanishing != lp::null_ci) - return vanishing; - - if (inf == lp::null_ci && sup == lp::null_ci) - return lp::null_ci; - - if (!constraint_is_false(inf) && !constraint_is_false(sup)) - return lp::null_ci; - - TRACE(arith, tout << "v" << v << " @ " << m_var2level[v] << "\n"); - - if (inf == lp::null_ci) { - SASSERT(sup != lp::null_ci); - // repair v by setting it below sup - auto f = factor(v, sup); - m_values[v] = floor(-value(f.q) / value(f.p)); - if (lp::lconstraint_kind::GT == m_constraints[sup].k) - m_values[v]--; - CTRACE(arith, constraint_is_false(sup), display_constraint(tout << m_values[v] << " ", sup) << "\n"); - return lp::null_ci; - } - if (sup == lp::null_ci) { - SASSERT(inf != lp::null_ci); - // repair v by setting it above inf - auto f = factor(v, inf); - m_values[v] = ceil(-value(f.q) / value(f.p)); - if (lp::lconstraint_kind::GT == m_constraints[inf].k) - m_values[v]++; - CTRACE(arith, constraint_is_false(inf), display_constraint(tout << m_values[v] << " ", inf) << "\n"); - return lp::null_ci; - } - - bool strict_lo = m_constraints[inf].k == lp::lconstraint_kind::GT; - bool strict_hi = m_constraints[sup].k == lp::lconstraint_kind::GT; - bool strict = strict_lo || strict_hi; - if (((!strict && lo <= hi) || lo < hi) && (!var_is_int(v) || ceil(lo) <= floor(hi))) { - // repair v by setting it between lo and hi assuming it is integral when v is integer - m_values[v] = var_is_int(v) ? ceil(lo) : ((lo + hi) / 2); - CTRACE(arith, constraint_is_false(sup) || constraint_is_false(inf), - tout << m_values[v] << " " << " between " << lo << " and " << hi << "\n"; - display_constraint(tout, sup) << "\n"; - display_constraint(tout, inf) << "\n";); - return lp::null_ci; - } - - TRACE(arith, tout << "cannot repair v" << v << " between " << lo << " and " << hi << " " << (lo > hi) << " is int " << var_is_int(v) << "\n"; - display_constraint(tout, inf) << "\n"; display_constraint(tout, sup) << "\n";); - auto f = factor(v, inf); - SASSERT(f.degree == 1); - auto p_value = value(f.p); - f.p *= pddm.mk_var(v); - auto [m, f_p] = f.p.var_factors(); - return resolve_variable(v, inf, sup, p_value, f, m, f_p); - } - - stellensatz::bound_info stellensatz::find_bounds(lpvar v) { - bound_info result; - auto &[lo, hi, inf, sup, van] = result; - for (auto ci : m_occurs[v]) { - CTRACE(arith_verbose, !m_active.contains(ci), display_constraint(tout << "skip inactive ", ci) << "\n"); - if (!m_active.contains(ci)) - continue; - - // consider only constraints where v is maximal - auto const &p = m_constraints[ci].p; - auto const &vars = p.free_vars(); - if (any_of(vars, [&](unsigned j) { return p.degree(j) == 1 && m_var2level[j] > m_var2level[v]; })) { - m_processed.insert(ci); - TRACE(arith_verbose, display_constraint(tout << "skip non-maximal ", ci) << "\n"); - continue; - } - auto f = factor(v, ci); - auto q_ge_0 = vanishing(v, f, ci); - if (q_ge_0 != lp::null_ci) { - if (m_processed.contains(ci)) - continue; - m_processed.insert(ci); - if (!constraint_is_true(q_ge_0)) { - van = q_ge_0; - return result; - } - TRACE(arith_verbose, display_constraint(tout << "vanished j" << v << " in ", ci) << "\n"; - display_constraint(tout << " to ", q_ge_0) << "\n";); - continue; - } - m_processed.insert(ci); - if (f.degree > 1) - continue; - TRACE(arith_verbose, display_constraint(tout << "process maximal ", ci) << "\n"); - auto p_value = value(f.p); - auto q_value = value(f.q); - SASSERT(f.degree == 1); - SASSERT(p_value != 0); - auto quot = -q_value / p_value; - if (p_value < 0) { - // p*x + q >= 0 => x <= -q / p if p < 0 - if (sup == lp::null_ci || hi > quot) { - hi = quot; - sup = ci; - } - else if (hi == quot && m_constraints[ci].k == lp::lconstraint_kind::GT) { - sup = ci; - } - } - else { - // p*x + q >= 0 => x >= -q / p if p > 0 - if (inf == lp::null_ci || lo < quot) { - lo = quot; - inf = ci; - } - else if (lo == quot && m_constraints[ci].k == lp::lconstraint_kind::GT) { - inf = ci; - } - } - } - return result; - } - lp::constraint_index stellensatz::vanishing(lpvar x, factorization const &f, lp::constraint_index ci) { if (f.p.is_zero()) return lp::null_ci; @@ -679,19 +530,6 @@ namespace nla { return new_ci; } - bool stellensatz::is_new_constraint(lp::constraint_index ci) const { - return ci == m_last_constraint; - } - - lp::lpvar stellensatz::select_variable_to_eliminate(lp::constraint_index ci) { - auto const& [p, k, j] = m_constraints[ci]; - lpvar best_var = lp::null_lpvar; - for (auto v : p.free_vars()) - if (best_var > v) - best_var = v; - return best_var; - } - unsigned stellensatz::degree_of_var_in_constraint(lpvar var, lp::constraint_index ci) const { return m_constraints[ci].p.degree(var); } @@ -736,7 +574,7 @@ namespace nla { external.shrink(i); auto [p, k, j] = m_constraints[max_ci]; while (m_constraints.size() > max_ci) - m_trail.pop_scope(1); + m_ctrail.pop_scope(1); for (auto const &[_p, _k, _j] : replay) { auto ci = add_constraint(_p, _k, _j); @@ -975,10 +813,8 @@ namespace nla { s.remove_occurs(ci); } }; - add_active(ci); - m_trail.push(undo_occurs(*this, ci)); + m_ctrail.push(undo_occurs(*this, ci)); m_has_occurs[ci] = true; - m_last_constraint = ci; auto const &con = m_constraints[ci]; for (auto v : con.p.free_vars()) m_occurs[v].push_back(ci); @@ -1010,7 +846,7 @@ namespace nla { } bool stellensatz::constraint_is_false(lp::constraint_index ci) const { - return ci != lp::null_ci && !constraint_is_true(m_constraints[ci]); + return ci != lp::null_ci && constraint_is_false(m_constraints[ci]); } bool stellensatz::constraint_is_true(constraint const &c) const { @@ -1054,20 +890,411 @@ namespace nla { || (p.val() <= 0 && k == lp::lconstraint_kind::GT)); } + bool stellensatz::constraint_is_bound_conflict(constraint const& c, u_dependency*& d) { + auto const& [p, k, j] = c; + scoped_dep_interval iv(m_di); + interval(p, iv); + if (iv->m_upper_inf) + return false; + if (iv->m_upper > 0) + return false; + bool is_negative = iv->m_upper < 0 || iv->m_upper_open; + bool is_conflict = k == lp::lconstraint_kind::GT || (k == lp::lconstraint_kind::GE && is_negative); + if (!is_conflict) + return false; + d = iv->m_upper_dep; + return true; + } + + bool stellensatz::var_is_bound_conflict(lpvar v) const { + if (!has_lo(v)) + return false; + if (!has_hi(v)) + return false; + if (lo_val(v) > hi_val(v)) + return true; + if (lo_val(v) < hi_val(v)) + return false; + return lo_is_strict(v) || hi_is_strict(v); + } + + void stellensatz::interval(dd::pdd p, scoped_dep_interval& iv) { + auto &m = iv.m(); + if (p.is_val()) { + m.set_lower(iv, p.val()); + m.set_lower_dep(iv, nullptr); + m.set_lower_is_open(iv, false); + m.set_lower_is_inf(iv, false); + m.set_upper(iv, p.val()); + m.set_upper_dep(iv, nullptr); + m.set_upper_is_open(iv, false); + m.set_upper_is_inf(iv, false); + return; + } + scoped_dep_interval x(m), lo(m), hi(m); + auto v = p.var(); + auto lb = m_lower[v]; + if (lb == UINT_MAX) { + m.set_lower_is_inf(x, true); + } + else { + auto lo = m_bounds[lb]; + m.set_lower_is_open(x, lo.m_kind == lp::lconstraint_kind::GT); + m.set_lower(x, lo.m_value); + m.set_lower_dep(x, lo.m_bound_justifications); + } + auto ub = m_upper[v]; + if (ub == UINT_MAX) { + m.set_upper_is_inf(x, true); + } + else { + auto hi = m_bounds[ub]; + m.set_upper_is_open(x, hi.m_kind == lp::lconstraint_kind::LT); + m.set_upper(x, hi.m_value); + m.set_upper_dep(x, hi.m_bound_justifications); + } + interval(p.hi(), hi); + interval(p.lo(), lo); + m.mul(x, hi, hi); + m.add(lo, hi, iv); + } + + lbool stellensatz::search() { + init_search(); + lbool is_sat = l_undef; + while (is_sat == l_undef && c().reslim().inc()) { + if (is_conflict()) + is_sat = resolve_conflict(); + else if (should_propagate()) + propagate(); + else if (!decide()) + is_sat = l_true; + } + return is_sat; + } + + void stellensatz::init_search() { + m_processed.reset(); + m_level2var.reset(); + m_var2level.reset(); + for (unsigned v = 0; v < m_values.size(); ++v) + m_level2var.push_back(v); + random_gen rand(c().random()); + shuffle(m_level2var.size(), m_level2var.begin(), rand); + // move up (or down?) variables that live in unsat constraints under monomials + for (auto const &c : m_constraints) + if (constraint_is_false(c)) + for (auto v : c.p.free_vars()) + move_up(v); + m_var2level.resize(m_values.size(), m_level2var.size()); + for (unsigned i = 0; i < m_level2var.size(); ++i) + m_var2level[m_level2var[i]] = i; + } + + + // walk m_bounds based on the propagated bounds + lbool stellensatz::resolve_conflict() { + SASSERT(is_conflict()); + m_conflict_marked.reset(); + mark_dependencies(m_conflict_dep); + unsigned conflict_level = 0; + m_core.reset(); + for (auto d : m_conflict_marked) { + auto &b = m_bounds[d]; + conflict_level = std::max(conflict_level, b.m_level); + } + bool found_decision = false; + while (!found_decision && !m_bounds.empty() && !m_conflict_marked.empty()) { + while (!m_conflict_marked.contains(m_bounds.size() - 1)) + pop_bound(); + + auto const &[v, k, value, level, last_bound, is_decision, dep, ci] = m_bounds.back(); + + mark_dependencies(dep); + m_conflict_marked.remove(m_bounds.size() - 1); + m_conflict = resolve_variable(v, m_conflict, ci); // adds assumptions.. + if (conflict_level == 0 && ci != lp::null_ci) + m_core.push_back(ci); + found_decision = is_decision; + } + SASSERT(found_decision || conflict_level == 0); + SASSERT(!found_decision || conflict_level != 0); + if (conflict_level == 0) + return l_false; + + auto [v, k, value, level, last_bound, is_decision, dep, ci] = m_bounds.back(); + SASSERT(is_decision); + m_conflict_marked.remove(m_bounds.size() - 1); + while (!m_bounds.empty() && !m_conflict_marked.contains(m_bounds.size() - 1)) + pop_bound(); + switch (k) { + case lp::lconstraint_kind::GE: + if (var_is_int(v)) { + k = lp::lconstraint_kind::LE; + value = value - 1; + } + else + k = lp::lconstraint_kind::LT; + break; + case lp::lconstraint_kind::GT: k = lp::lconstraint_kind::LE; break; + case lp::lconstraint_kind::LT: k = lp::lconstraint_kind::GE; break; + case lp::lconstraint_kind::LE: + if (var_is_int(v)) { + k = lp::lconstraint_kind::GE; + value = value + 1; + } + else + k = lp::lconstraint_kind::GT; + break; + } + dep = nullptr; + for (auto d : m_conflict_marked) + dep = m_dm.mk_join(m_dm.mk_leaf(d), dep); + bool is_upper = (k == lp::lconstraint_kind::LE) || k == lp::lconstraint_kind::LT; + last_bound = is_upper ? m_upper[v] : m_lower[v]; + (is_upper ? m_upper[v] : m_lower[v]) = m_bounds.size(); + m_bounds.push_back(bound_info(v, k, value, m_level, last_bound, dep, m_conflict)); + // set the value within bounds if the bounds are consistent. + set_in_bounds(v); + return l_undef; + } + + void stellensatz::mark_dependencies(u_dependency* d) { + unsigned_vector deps; + m_dm.linearize(d, deps); + for (auto d : deps) + m_conflict_marked.insert(d); + } + + void stellensatz::pop_bound() { + auto const &[v, k, value, level, last_bound, is_decision, dep, ci] = m_bounds.back(); + bool is_upper = k == lp::lconstraint_kind::LE || k == lp::lconstraint_kind::LT; + (is_upper ? m_upper[v] : m_lower[v]) = last_bound; + if (is_decision) { + m_dm.pop_scope(1); + m_vtrail.pop_scope(1); + } + m_bounds.pop_back(); + } + + void stellensatz::propagate() { + if (m_prop_qhead == m_bounds.size()) + return; + m_vtrail.push(value_trail(m_prop_qhead)); + for (; m_prop_qhead < m_bounds.size(); ++m_prop_qhead) { + if (!c().reslim().inc()) + return; + auto v = m_bounds[m_prop_qhead].m_var; + if (var_is_bound_conflict(v)) { + set_conflict(v); + return; + } + for (auto ci : m_occurs[v]) { + if (constraint_is_true(ci)) + continue; + + // detect conflict or propagate + u_dependency *d = nullptr; + if (constraint_is_bound_conflict(ci, d)) { + set_conflict(ci, d); + return; + } + lpvar w; + rational value; + lp::lconstraint_kind k; + unsigned_vector deps; + if (constraint_is_propagating(ci, d, w, value, k)) { + TRACE(arith, display_constraint(tout, ci) << "\n"; + tout << "v" << w << " " << k << " " << value << "\n"; m_dm.linearize(d, deps); + for (auto d : deps) display_bound(tout, d);); + + bool is_upper = k == lp::lconstraint_kind::LE || k == lp::lconstraint_kind::LT; + bool is_strict = k == lp::lconstraint_kind::LT || k == lp::lconstraint_kind::GT; + auto last_bound = is_upper ? m_upper[w] : m_lower[w]; + (is_upper ? m_upper[w] : m_lower[w]) = m_bounds.size(); + m_bounds.push_back(bound_info(w, k, value, m_level, last_bound, d, ci)); + if (!is_strict) + m_values[w] = value; + else { + SASSERT(has_lo(w) || has_hi(w)); + if (!has_lo(w)) + m_values[w] = hi_val(w) - 1; + else if (!has_hi(w)) + m_values[w] = lo_val(w) + 1; + else + m_values[w] = (lo_val(w) + hi_val(w)) / 2; + } + SASSERT(well_formed(w)); + continue; + } + // constraint is false, but not propagating + } + } + } + + bool stellensatz::decide() { + rational value; + lp::lconstraint_kind k; + // first attempt to repair variables in some order + // then if hitting a variable that cannot be repaired, create a decision based on the value conflict. + // resume from some level? + for (unsigned level = 0; level < m_level2var.size(); ++level) { + auto w = m_level2var[level]; + lp::constraint_index ci = lp::null_ci; + bool is_sat = repair_variable(w, value, k, ci); + if (is_sat) + continue; + if (ci != lp::null_ci) { + // backtrack decision to max variable in ci + level = max_level(m_constraints[ci]) - 1; + continue; + } + SASSERT(k == lp::lconstraint_kind::LE || k == lp::lconstraint_kind::GE); + bool is_upper = k == lp::lconstraint_kind::LE; + m_vtrail.push_scope(); + m_vtrail.push(value_trail(m_level)); + m_dm.push_scope(); + m_level++; + auto last_bound = is_upper ? m_upper[w] : m_lower[w]; + (is_upper ? m_upper[w] : m_lower[w]) = m_bounds.size(); + m_bounds.push_back(bound_info(w, k, value, m_level, last_bound)); + m_values[w] = value; + SASSERT(well_formed(w)); + return true; + } + return false; + } + + unsigned stellensatz::max_level(constraint const &c) const { + unsigned level = 0; + auto const &vars = c.p.free_vars(); + for (auto v : vars) + level = std::max(level, m_var2level[v]); + return level; + } + + bool stellensatz::constraint_is_propagating(lp::constraint_index ci, u_dependency*& d, lpvar& v, rational& value, lp::lconstraint_kind& k) { + auto [p, ck, j] = m_constraints[ci]; + for (auto x : p.free_vars()) { + auto f = factor(x, ci); + if (f.degree > 1) + continue; + scoped_dep_interval ivp(m_di), ivq(m_di); + interval(f.p, ivp); + interval(-f.q, ivq); + // p > 0: + if (!ivq->m_lower_inf && !ivp->m_lower_inf && ivp->m_lower > 0) { + // v >= -q / p + auto quot = ivq->m_lower / ivp->m_lower; + if (var_is_int(x)) + quot = ceil(quot); + bool is_strict = !var_is_int(x) && (ck == lp::lconstraint_kind::GT || ivq->m_lower_open || ivp->m_lower_open); + if (!has_lo(x) || quot > lo_val(x) || (!lo_is_strict(x) && quot == lo_val(x) && is_strict)) { + d = m_dm.mk_join(ivq->m_lower_dep, ivp->m_lower_dep); + k = is_strict ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE; + v = x; + return true; + } + } + // p < 0 + if (ivp->m_upper_inf && ivp->m_upper < 0 && !ivq->m_upper_inf) { + // v <= -q / p + auto quot = ivq->m_upper / ivp->m_upper; + if (var_is_int(x)) + quot = floor(quot); + bool is_strict = + !var_is_int(x) && (ck == lp::lconstraint_kind::GT || ivq->m_upper_open || ivp->m_upper_open); + if (!has_hi(x) || quot < hi_val(x) || (!hi_is_strict(x) && quot == hi_val(x) && is_strict)) { + d = m_dm.mk_join(ivq->m_upper_dep, ivp->m_upper_dep); + k = is_strict ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE; + v = x; + return true; + } + } + } + return false; + } + + bool stellensatz::well_formed() const { + for (unsigned v = 0; v < num_vars(); ++v) + if (!well_formed(v)) + return false; + return true; + } + + bool stellensatz::well_formed(lpvar v) const { + if (var_is_bound_conflict(v)) + return true; + auto const &value = m_values[v]; + if (var_is_int(v) && value.is_int()) + return false; + // values of variables are within bounds + if (has_lo(v) && value < lo_val(v)) + return false; + if (has_lo(v) && lo_is_strict(v) && value == lo_val(v)) + return false; + if (has_lo(v) && lo_kind(v) != lp::lconstraint_kind::GE && lo_kind(v) != lp::lconstraint_kind::GT) + return false; + if (has_hi(v) && value > hi_val(v)) + return false; + if (has_hi(v) && hi_is_strict(v) && value == hi_val(v)) + return false; + if (has_hi(v) && hi_kind(v) != lp::lconstraint_kind::LE && hi_kind(v) != lp::lconstraint_kind::LT) + return false; + return true; + } + + std::ostream& stellensatz::display_bound(std::ostream& out, unsigned i) const { + auto const &[v, k, val, last_bound, level, is_decision, d, ci] = m_bounds[i]; + out << level << ":" << i << ": v" << v << " " << k << " " + << val << " " + << ((last_bound == UINT_MAX) ? -1 : last_bound) + << (is_decision ? " d " : ""); + unsigned_vector deps; + m_dm.linearize(d, deps); + out << deps; + if (ci != lp::null_ci) + out << " (" << ci << ")"; + out << "\n"; + return out; + } + std::ostream& stellensatz::display(std::ostream& out) const { - #if 0 -// m_solver.lra().display(out); - for (auto const& [vars, v] : m_vars2mon) { - out << "j" << v << " := "; - display_product(out, vars); + + for (unsigned i = 0; i < m_bounds.size(); ++i) + display_bound(out, i); + + for (unsigned v = 0; v < num_vars(); ++v) { + out << "v" << v << " " << m_values[v] << " "; + if (has_lo(v)) { + if (lo_is_strict(v)) + out << "("; + else + out << "["; + out << lo_val(v); + } + else + out << "(-oo"; + out << ", "; + if (has_hi(v)) { + out << hi_val(v); + if (hi_is_strict(v)) + out << ")"; + else + out << "]"; + } + else + out << "+oo)"; + if (has_lo(v)) + out << " " << m_lower[v]; + if (has_hi(v)) + out << " " << m_upper[v]; + out << "\n"; } - #endif for (unsigned ci = 0; ci < m_constraints.size(); ++ci) { display_constraint(out, ci) << "\n"; - out << "\t<- "; - display(out, m_constraints[ci].j); - out << "\n"; + display(out << "\t<- ", m_constraints[ci].j) << "\n"; } return out; } @@ -1107,7 +1334,9 @@ namespace nla { if (ci == lp::null_ci) return out << "(null) "; bool is_true = constraint_is_true(ci); - return display_constraint(out << "(" << ci << ") ", m_constraints[ci]) << (is_true ? " [true]" : " [false]"); + auto const& [p, k, j] = m_constraints[ci]; + return display_constraint(out << "(" << ci << ") ", m_constraints[ci]) + << (is_true ? " [true] " : " [false] ") << "(" << value(p) << " " << k << " 0)"; } std::ostream& stellensatz::display_constraint(std::ostream& out, constraint const &c) const { @@ -1331,4 +1560,183 @@ namespace nla { values[i] = _values[i]; TRACE(arith, for (unsigned i = 0; i < sz; ++i) tout << "j" << i << " := " << values[i] << "\n";); } + + // + // has_lo(v), has_hi(v), bounds by constraints + // new bound should be consistent with lo_val/hi_val + // such conflicts should be caught by bounds propagation already + // + bool stellensatz::repair_variable(lp::lpvar &v, rational &r, lp::lconstraint_kind &k, lp::constraint_index &ci) { + ci = lp::null_ci; + auto [lo, hi, inf, sup, vanishing] = find_bounds(v); + TRACE(arith, tout << "bounds for v" << v << " @ " << m_var2level[v] << " : " << lo << " to " << hi << "\n"; + display_constraint(tout << " inf ", inf) << "\n"; display_constraint(tout << " sup ", sup) << "\n"; + display_constraint(tout << " vanishing ", vanishing) << "\n";); + if (vanishing != lp::null_ci) { + ci = vanishing; + return false; + } + + if (inf == lp::null_ci && sup == lp::null_ci) + return true; + + if (!constraint_is_false(inf) && !constraint_is_false(sup)) + return true; + + TRACE(arith, tout << "v" << v << " @ " << m_var2level[v] << "\n"); + + if (constraint_is_false(sup)) { + SASSERT(sup != lp::null_ci); + // repair v by setting it below sup + auto f = factor(v, sup); + auto hi = -value(f.q) / value(f.p); + r = floor(hi); + if (lp::lconstraint_kind::GT == m_constraints[sup].k) { + // todo - be more refined for adusting bound for reals + r--; + } + k = lp::lconstraint_kind::LE; + if (in_bounds(v, r)) { + m_values[v] = r; + if (!constraint_is_false(inf)) + return true; + } + find_split(v, r, k, sup); + } + else { + SASSERT(inf != lp::null_ci); + // repair v by setting it below sup + auto f = factor(v, inf); + auto lo = -value(f.q) / value(f.p); + r = ceil(hi); + if (lp::lconstraint_kind::GT == m_constraints[inf].k) { + // todo - be more refined for adjusting bound split for reals + r++; + } + k = lp::lconstraint_kind::LE; + if (in_bounds(v, r)) { + m_values[v] = r; + if (!constraint_is_false(sup)) + return true; + } + find_split(v, r, k, inf); + } + return false; + } + + + void stellensatz::find_split(lpvar & v, rational & r, lp::lconstraint_kind & k, lp::constraint_index ci) { + unsigned n = 0; + for (auto w : m_constraints[ci].p.free_vars()) { + if (is_fixed(w)) + continue; + if (c().random(++n) == 0) { + r = m_values[w]; + if (has_hi(w) && hi_val(w) == r) + k = lp::lconstraint_kind::GE; + else + k = lp::lconstraint_kind::LE; + v = w; + } + } + } + + void stellensatz::set_in_bounds(lpvar v) { + if (in_bounds(v)) + return; + if (!has_lo(v)) + m_values[v] = hi_is_strict(v) ? hi_val(v) - 1 : hi_val(v); + else if (!has_hi(v)) + m_values[v] = lo_is_strict(v) ? lo_val(v) + 1 : lo_val(v); + else if (lo_is_strict(v) || hi_is_strict(v)) + m_values[v] = (lo_val(v) + hi_val(v)) / 2; + else + m_values[v] = lo_val(v); + } + + bool stellensatz::in_bounds(lpvar v, rational const& value) const { + if (has_lo(v)) { + if (value < lo_val(v)) + return false; + if (lo_is_strict(v) && value <= lo_val(v)) + return false; + } + if (has_hi(v)) { + if (value > hi_val(v)) + return false; + if (hi_is_strict(v) && value >= hi_val(v)) + return false; + } + return true; + } + + stellensatz::repair_var_info stellensatz::find_bounds(lpvar v) { + repair_var_info result; + auto &[lo, hi, inf, sup, van] = result; + for (auto ci : m_occurs[v]) { + + // consider only constraints where v is maximal + auto const &p = m_constraints[ci].p; + auto const &vars = p.free_vars(); + if (any_of(vars, [&](unsigned j) { return p.degree(j) == 1 && m_var2level[j] > m_var2level[v]; })) { + TRACE(arith_verbose, display_constraint(tout << "skip non-maximal ", ci) << "\n"); + continue; + } + auto f = factor(v, ci); + auto q_ge_0 = vanishing(v, f, ci); + if (q_ge_0 != lp::null_ci) { + if (m_processed.contains(ci)) + continue; + if (!constraint_is_true(q_ge_0)) { + m_processed.insert(ci); + van = q_ge_0; + return result; + } + TRACE(arith_verbose, display_constraint(tout << "vanished j" << v << " in ", ci) << "\n"; + display_constraint(tout << " to ", q_ge_0) << "\n";); + continue; + } + if (f.degree > 1) + continue; + TRACE(arith_verbose, display_constraint(tout << "process maximal ", ci) << "\n"); + auto p_value = value(f.p); + auto q_value = value(f.q); + SASSERT(f.degree == 1); + SASSERT(p_value != 0); + auto quot = -q_value / p_value; + if (p_value < 0) { + // p*x + q >= 0 => x <= -q / p if p < 0 + if (sup == lp::null_ci || hi > quot) { + hi = quot; + sup = ci; + } + else if (hi == quot && m_constraints[ci].k == lp::lconstraint_kind::GT) { + sup = ci; + } + } + else { + // p*x + q >= 0 => x >= -q / p if p > 0 + if (inf == lp::null_ci || lo < quot) { + lo = quot; + inf = ci; + } + else if (lo == quot && m_constraints[ci].k == lp::lconstraint_kind::GT) { + inf = ci; + } + } + } + return result; + } + + void stellensatz::move_up(lpvar v) { + auto l = m_var2level[v]; + if (l == 0) + return; + auto l0 = c().random(l); + auto w = m_level2var[l0]; + m_var2level[v] = l0; + m_var2level[w] = l; + m_level2var[l0] = v; + m_level2var[l] = w; + } } diff --git a/src/math/lp/nla_stellensatz.cpp.2 b/src/math/lp/nla_stellensatz.cpp.2 new file mode 100644 index 000000000..35793d800 --- /dev/null +++ b/src/math/lp/nla_stellensatz.cpp.2 @@ -0,0 +1,1336 @@ +/*++ + Copyright (c) 2025 Microsoft Corporation + + + vanishing(v, p) = (q, r) with q = 0, r >= 0 if p := q*v + r, M(q) = 0 + + + --*/ + +#include "params/smt_params_helper.hpp" +#include "math/dd/pdd_eval.h" +#include "math/lp/nla_core.h" +#include "math/lp/nla_coi.h" +#include "math/lp/nla_stellensatz.h" +#include +#include +#include +#include +#include +#include +#include "explanation.h" +#include "int_solver.h" +#include "lar_solver.h" +#include "lia_move.h" +#include "lp_settings.h" +#include "lp_types.h" +#include "nla_common.h" +#include "nla_defs.h" +#include "nla_types.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace nla { + + lp::lconstraint_kind join(lp::lconstraint_kind k1, lp::lconstraint_kind k2) { + if (k1 == k2) + return k1; + if (k1 == lp::lconstraint_kind::EQ) + return k2; + if (k2 == lp::lconstraint_kind::EQ) + return k1; + 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::LT; + 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::GT; + UNREACHABLE(); + 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; + } + + lpvar stellensatz::monomial_factory::mk_monomial(lp::lar_solver &lra, svector const &vars) { + lpvar v = lp::null_lpvar; + if (vars.empty()) + return v; + if (vars.size() == 1) + return vars[0]; + svector _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 const &vars) { + lpvar v = lp::null_lpvar; + if (vars.empty()) + return v; + if (vars.size() == 1) + return vars[0]; + svector _vars(vars); + std::sort(_vars.begin(), _vars.end()); + if (m_vars2mon.find(_vars, v)) + return v; + return lp::null_lpvar; + } + + void stellensatz::monomial_factory::init(lpvar v, svector const &_vars) { + svector 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(core->lra_solver().number_of_vars()) { + } + + lbool stellensatz::saturate() { + unsigned num_conflicts = 0, num_constraints = 0; + init_solver(); + TRACE(arith, display(tout << "stellensatz before saturation\n")); + start_saturate: + m_last_constraint = lp::null_ci; + num_constraints = m_constraints.size(); + lbool r = model_repair(); + + if (num_constraints == m_constraints.size()) + return l_undef; + + if (r == l_undef) + r = m_solver.solve(m_core); + TRACE(arith, display(tout << "stellensatz after saturation " << r << "\n")); + if (r == l_false) { + ++num_conflicts; + IF_VERBOSE(2, verbose_stream() << "(nla.stellensatz :conflicts " << num_conflicts << ":constraints " << m_constraints.size() << ")\n"); + if (num_conflicts >= m_config.max_conflicts) + return l_undef; + while (backtrack(m_core)) { + ++c().lp_settings().stats().m_stellensatz.m_num_backtracks; + lbool rb = m_solver.solve(m_core); + if (rb == l_false) + continue; + if (rb == l_undef) + return rb; + m_solver.update_values(m_values); + goto start_saturate; + } + ++c().lp_settings().stats().m_stellensatz.m_num_conflicts; + conflict(m_core); + } + if (r == l_true && !set_model()) + r = l_undef; + return r; + } + + void stellensatz::pop_constraint() { + auto const &[p, k, j] = m_constraints.back(); + m_constraint_index.erase({p.index(), k}); + m_constraints.pop_back(); + } + + void stellensatz::init_solver() { + m_trail.reset(); + m_monomial_factory.reset(); + m_coi.init(); + m_values.reset(); + init_vars(); + init_occurs(); + } + + void stellensatz::init_vars() { + auto const& src = c().lra_solver(); + auto sz = src.number_of_vars(); + for (unsigned v = 0; v < sz; ++v) { + m_values.push_back(c().val(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)) { + 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); + add_var_bound(v, k, rhs, external_justification(dep)); + } + if (src.column_has_upper_bound(v)) { + auto hi_bound = src.get_upper_bound(v); + SASSERT(hi_bound.y <= 0); + 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); + add_var_bound(v, k, rhs, external_justification(dep)); + } + } + m_occurs.reserve(sz); + } + + // set the model based on m_values computed by the solver + bool stellensatz::set_model() { + if (any_of(m_constraints, [&](auto const &c) { return !constraint_is_true(c); })) + return false; + auto &src = c().lra_solver(); + c().m_use_nra_model = true; + m_values.reserve(src.number_of_vars()); + for (unsigned v = 0; v < src.number_of_vars(); ++v) { + if (c().is_monic_var(v)) { + auto& mon = c().emons()[v]; + rational val(1); + for (auto w : mon.vars()) + val *= m_values[w]; + m_values[v] = val; + } + else if (src.column_has_term(v)) { + auto const& t = src.get_term(v); + rational val(0); + for (auto cv : t) + val += m_values[cv.j()] * cv.coeff(); + m_values[v] = val; + } + else if (m_coi.vars().contains(v)) { + // variable is in coi, m_values[v] is set. + } + else { + m_values[v] = c().val(v); + } + c().m_nra.set_value(v, m_values[v]); + } + return true; + } + + 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, m_monomial_factory.get_monomial(vars)); + } + return to; + } + + bool stellensatz::has_term_offset(dd::pdd const &p) { + for (auto const &[coeff, vars] : p) + if (!vars.empty() && lp::null_lpvar == m_monomial_factory.get_monomial(vars)) + return false; + return true; + } + + void stellensatz::init_monomial(unsigned mon_var) { + auto &mon = c().emons()[mon_var]; + 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) { + 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; + } + if (k == lp::lconstraint_kind::GT && is_int(p)) { + k = lp::lconstraint_kind::GE; + p -= rational(1); + } + lp::constraint_index ci = lp::null_ci; + if (m_constraint_index.find({p.index(), k}, ci)) + return ci; + ci = m_constraints.size(); + m_constraints.push_back({p, k, j }); + m_constraint_index.insert({p.index(), k}, ci); + ++c().lp_settings().stats().m_stellensatz.m_num_constraints; + + m_has_occurs.reserve(ci + 1); + struct undo_constraint : public trail { + stellensatz &s; + undo_constraint(stellensatz& s): s(s) {} + void undo() override { + s.pop_constraint(); + } + }; + m_trail.push_scope(); + m_trail.push(undo_constraint(*this)); + return ci; + } + + void stellensatz::add_active(lp::constraint_index ci) { + if (m_active.contains(ci)) + return; + if (m_constraints[ci].p.degree() > m_config.max_degree) + return; + m_active.insert(ci); + } + + lp::constraint_index stellensatz::resolve_variable(lpvar x, lp::constraint_index ci, lp::constraint_index other_ci, + rational const &p_value, factorization const &f, unsigned_vector const &_m1, + dd::pdd _f_p) { + TRACE(arith, tout << "resolve j" << x << " between ci " << ci << " and other_ci " << other_ci << "\n"); + if (ci == other_ci) + return lp::null_ci; + if (f.degree != 1) + return lp::null_ci; // future - could handle this + auto const &[other_p, other_k, other_j] = m_constraints[other_ci]; + auto const &[p, k, j] = m_constraints[ci]; + auto g = factor(x, other_ci); + if (g.degree != 1) + return lp::null_ci; // not handling degree > 1 + auto p_value2 = value(g.p); + // + // check that p_value and p_value2 have different signs + // check that other_p is disjoint from tabu + // compute overlaps extending x + // + SASSERT(g.degree > 0); + SASSERT(p_value != 0); + if (g.degree > f.degree) // future: could handle this too by considering tabu to be a map into degrees. + return lp::null_ci; + if (p_value > 0 && p_value2 > 0) + return lp::null_ci; + if (p_value < 0 && p_value2 < 0) + return lp::null_ci; + + for (unsigned i = 0; i < g.degree; ++i) + g.p *= to_pdd(x); + + 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 = value(f_p) < 0; + SASSERT(value(f_p) != 0); + + // m1m2 * f_p + f.q >= 0 + // m1m2 * g_p + g.q >= 0 + // + lp::constraint_index ci_a, ci_b; + + bool g_strict = other_k == lp::lconstraint_kind::GT; + bool f_strict = k == lp::lconstraint_kind::GT; + if (g_p.is_val()) + ci_a = multiply(ci, pddm.mk_val(g_p.val())); + else if (!sign_f) + ci_a = multiply(ci, assume(g_p, f_strict ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE)); + else + ci_a = multiply(ci, assume(g_p, f_strict ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE)); + + if (f_p.is_val()) + ci_b = multiply(other_ci, pddm.mk_val(f_p.val())); + else if (sign_f) + ci_b = multiply(other_ci, assume(f_p, g_strict ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE)); + else + ci_b = multiply(other_ci, assume(f_p, g_strict ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE)); + + auto new_ci = add(ci_a, ci_b); + + ++c().lp_settings().stats().m_stellensatz.m_num_resolutions; + + + TRACE(arith, tout << "eliminate j" << x << ":\n"; + display_constraint(tout, ci) << "\n"; + display_constraint(tout, other_ci) << "\n"; + display_constraint(tout, ci_a) << "\n"; + display_constraint(tout, ci_b) << "\n"; + display_constraint(tout, new_ci) << "\n"); + + if (constraint_is_conflict(new_ci)) + return new_ci; + + new_ci = factor(new_ci); + + init_occurs(new_ci); + + CTRACE(arith, !is_new_constraint(new_ci), display_constraint(tout << "not new ", new_ci) << "\n"); + CTRACE(arith, is_new_constraint(new_ci), display_constraint(tout << "new ", new_ci) << "\n"); + return new_ci; + } + + bool stellensatz::conflict(lp::constraint_index ci) { + if (!constraint_is_conflict(ci)) + return false; + m_core.reset(); + m_core.push_back(ci); + return true; + } + + void stellensatz::conflict(svector const& core) { + lp::explanation ex; + lemma_builder new_lemma(c(), "stellensatz"); + m_constraints_in_conflict.reset(); + for (auto ci : core) + explain_constraint(new_lemma, 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++; + } + + lbool stellensatz::model_repair() { + m_active.reset(); + m_processed.reset(); + for (lp::constraint_index ci = 0; ci < m_constraints.size(); ++ci) + m_active.insert(ci); + m_level2var.reset(); + m_var2level.reset(); + for (unsigned v = 0; v < m_values.size(); ++v) + m_level2var.push_back(v); + random_gen rand(c().random()); + shuffle(m_level2var.size(), m_level2var.begin(), rand); + m_var2level.resize(m_values.size(), m_level2var.size()); + for (unsigned i = 0; i < m_level2var.size(); ++i) + m_var2level[m_level2var[i]] = i; + return model_repair_loop(); + } + + lbool stellensatz::model_repair_loop() { + unsigned level = 0; + while (c().reslim().inc()) { + if (level >= m_level2var.size()) + break; + m_last_constraint = lp::null_ci; + auto ci = model_repair(m_level2var[level]); + TRACE(arith, + display_constraint(tout << "after repairing j" << m_level2var[level] << " @ " << level << " ", ci) + << "\nconflict " + << conflict(ci) << "\n"); + if (ci == lp::null_ci) + ++level; + else if (conflict(ci)) + return l_false; + else if (m_processed.contains(ci) || constraint_is_trivial(ci)) + ++level; + else + level = max_level(m_constraints[ci]); + } + return all_of(m_constraints, [&](constraint const &c) { return constraint_is_true(c); }) ? l_true : l_undef; + } + + unsigned stellensatz::max_level(constraint const &c) const { + unsigned level = 0; + auto const &vars = c.p.free_vars(); + for (auto v : vars) + if (true || c.p.degree(v) == 1) + level = std::max(level, m_var2level[v]); + // display_constraint(verbose_stream(), c) << " " << level << "\n"; + return level; + } + + lp::constraint_index stellensatz::model_repair(lp::lpvar v) { + auto [lo, hi, inf, sup, vanishing] = find_bounds(v); + TRACE(arith, tout << "bounds for v" << v << " @ " << m_var2level[v] << " : " << lo << " to " << hi << "\n"; + display_constraint(tout << " inf ", inf) << "\n"; + display_constraint(tout << " sup ", sup) << "\n"; + display_constraint(tout << " vanishing ", vanishing) << "\n";); + if (vanishing != lp::null_ci) + return vanishing; + + if (inf == lp::null_ci && sup == lp::null_ci) + return lp::null_ci; + + if (!constraint_is_false(inf) && !constraint_is_false(sup)) + return lp::null_ci; + + TRACE(arith, tout << "v" << v << " @ " << m_var2level[v] << "\n"); + + if (inf == lp::null_ci) { + SASSERT(sup != lp::null_ci); + // repair v by setting it below sup + auto f = factor(v, sup); + m_values[v] = floor(-value(f.q) / value(f.p)); + if (lp::lconstraint_kind::GT == m_constraints[sup].k) + m_values[v]--; + CTRACE(arith, constraint_is_false(sup), display_constraint(tout << m_values[v] << " ", sup) << "\n"); + return lp::null_ci; + } + if (sup == lp::null_ci) { + SASSERT(inf != lp::null_ci); + // repair v by setting it above inf + auto f = factor(v, inf); + m_values[v] = ceil(-value(f.q) / value(f.p)); + if (lp::lconstraint_kind::GT == m_constraints[inf].k) + m_values[v]++; + CTRACE(arith, constraint_is_false(inf), display_constraint(tout << m_values[v] << " ", inf) << "\n"); + return lp::null_ci; + } + + bool strict_lo = m_constraints[inf].k == lp::lconstraint_kind::GT; + bool strict_hi = m_constraints[sup].k == lp::lconstraint_kind::GT; + bool strict = strict_lo || strict_hi; + if (((!strict && lo <= hi) || lo < hi) && (!var_is_int(v) || ceil(lo) <= floor(hi))) { + // repair v by setting it between lo and hi assuming it is integral when v is integer + m_values[v] = var_is_int(v) ? ceil(lo) : ((lo + hi) / 2); + CTRACE(arith, constraint_is_false(sup) || constraint_is_false(inf), + tout << m_values[v] << " " << " between " << lo << " and " << hi << "\n"; + display_constraint(tout, sup) << "\n"; + display_constraint(tout, inf) << "\n";); + return lp::null_ci; + } + + TRACE(arith, tout << "cannot repair v" << v << " between " << lo << " and " << hi << " " << (lo > hi) << " is int " << var_is_int(v) << "\n"; + display_constraint(tout, inf) << "\n"; display_constraint(tout, sup) << "\n";); + auto f = factor(v, inf); + SASSERT(f.degree == 1); + auto p_value = value(f.p); + f.p *= pddm.mk_var(v); + auto [m, f_p] = f.p.var_factors(); + return resolve_variable(v, inf, sup, p_value, f, m, f_p); + } + + stellensatz::bound_info stellensatz::find_bounds(lpvar v) { + bound_info result; + auto &[lo, hi, inf, sup, van] = result; + for (auto ci : m_occurs[v]) { + CTRACE(arith_verbose, !m_active.contains(ci), display_constraint(tout << "skip inactive ", ci) << "\n"); + if (!m_active.contains(ci)) + continue; + + // consider only constraints where v is maximal + auto const &p = m_constraints[ci].p; + auto const &vars = p.free_vars(); + if (any_of(vars, [&](unsigned j) { return p.degree(j) == 1 && m_var2level[j] > m_var2level[v]; })) { + m_processed.insert(ci); + TRACE(arith_verbose, display_constraint(tout << "skip non-maximal ", ci) << "\n"); + continue; + } + auto f = factor(v, ci); + auto q_ge_0 = vanishing(v, f, ci); + if (q_ge_0 != lp::null_ci) { + if (m_processed.contains(ci)) + continue; + m_processed.insert(ci); + if (!constraint_is_true(q_ge_0)) { + van = q_ge_0; + return result; + } + TRACE(arith_verbose, display_constraint(tout << "vanished j" << v << " in ", ci) << "\n"; + display_constraint(tout << " to ", q_ge_0) << "\n";); + continue; + } + m_processed.insert(ci); + if (f.degree > 1) + continue; + TRACE(arith_verbose, display_constraint(tout << "process maximal ", ci) << "\n"); + auto p_value = value(f.p); + auto q_value = value(f.q); + SASSERT(f.degree == 1); + SASSERT(p_value != 0); + auto quot = -q_value / p_value; + if (p_value < 0) { + // p*x + q >= 0 => x <= -q / p if p < 0 + if (sup == lp::null_ci || hi > quot) { + hi = quot; + sup = ci; + } + else if (hi == quot && m_constraints[ci].k == lp::lconstraint_kind::GT) { + sup = ci; + } + } + else { + // p*x + q >= 0 => x >= -q / p if p > 0 + if (inf == lp::null_ci || lo < quot) { + lo = quot; + inf = ci; + } + else if (lo == quot && m_constraints[ci].k == lp::lconstraint_kind::GT) { + inf = ci; + } + } + } + return result; + } + + lp::constraint_index stellensatz::vanishing(lpvar x, factorization const &f, lp::constraint_index ci) { + if (f.p.is_zero()) + return lp::null_ci; + auto p_value = value(f.p); + if (p_value != 0) + return lp::null_ci; + + ++c().lp_settings().stats().m_stellensatz.m_num_vanishings; + // add p = 0 as assumption and reduce to q. + auto p_is_0 = assume(f.p, lp::lconstraint_kind::EQ); + // ci & -p_is_0*x^f.degree => new_ci + dd::pdd r = pddm.mk_val(rational(-1)); + for (unsigned i = 0; i < f.degree; ++i) + 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 << "j" << x << " ", ci) << "\n"; + display_constraint(tout << "reduced ", new_ci) << "\n"; + display_constraint(tout, p_is_0) << "\n"); + init_occurs(new_ci); + return new_ci; + } + + lp::constraint_index stellensatz::factor(lp::constraint_index ci) { + auto const &[p, k, j] = m_constraints[ci]; + auto [vars, q] = p.var_factors(); // p = vars * q + + bool is_gt = k == lp::lconstraint_kind::GT; + unsigned i = 0; + for (auto v : vars) + if (is_gt || value(v) != 0) + vars[i++] = v; + else + q *= pddm.mk_var(v); + vars.shrink(i); + if (vars.empty()) + return ci; + if (vars.size() == 1 && q.is_val()) + return ci; + if (q.is_val()) { + q *= pddm.mk_var(vars.back()); + vars.pop_back(); + } + + vector muls; + muls.push_back(q); + for (auto v : vars) + muls.push_back(muls.back() * pddm.mk_var(v)); + auto new_ci = ci; + SASSERT(muls.back() == p); + int sign = 1; + for (unsigned i = vars.size(); i-- > 0;) { + auto pv = pddm.mk_var(vars[i]); + auto val = value(vars[i]); + auto k1 = val == 0 ? lp::lconstraint_kind::GE : (val > 0 ? lp::lconstraint_kind::GT : lp::lconstraint_kind::LT); + if (val < 0) + sign = -sign; + new_ci = divide(new_ci, assume(pv, k1), sign * muls[i]); + TRACE(arith_verbose, display_constraint(tout, new_ci) << "\n"; display_constraint(tout, assume(pv, k1)) << "\n";); + } + TRACE(arith, display_constraint(tout << "factor ", new_ci) << "\n"); + return new_ci; + } + + bool stellensatz::is_new_constraint(lp::constraint_index ci) const { + return ci == m_last_constraint; + } + + lp::lpvar stellensatz::select_variable_to_eliminate(lp::constraint_index ci) { + auto const& [p, k, j] = m_constraints[ci]; + lpvar best_var = lp::null_lpvar; + for (auto v : p.free_vars()) + if (best_var > v) + best_var = v; + return best_var; + } + + unsigned stellensatz::degree_of_var_in_constraint(lpvar var, lp::constraint_index ci) const { + return m_constraints[ci].p.degree(var); + } + + stellensatz::factorization stellensatz::factor(lpvar v, lp::constraint_index ci) { + auto const& [p, k, j] = m_constraints[ci]; + auto degree = degree_of_var_in_constraint(v, ci); + dd::pdd lc(pddm), rest(pddm); + p.factor(v, degree, lc, rest); + return {degree, lc, rest}; + } + + // + // check if core depends on an assumption + // identify the maximal assumption + // undo m_constraints down to max_ci. + // negate max_ci + // propagate it using remaining external and assumptions. + // find a new satisfying assignment (if any) before continuing. + // + bool stellensatz::backtrack(svector const &core) { + m_constraints_in_conflict.reset(); + svector external, assumptions; + for (auto ci : core) + explain_constraint(ci, external, assumptions); + if (assumptions.empty()) + return false; + lp::constraint_index max_ci = 0; + for (auto ci : assumptions) + max_ci = std::max(ci, max_ci); + TRACE(arith, tout << "max " << max_ci << " external " << external << " assumptions " << assumptions << "\n"; + display_constraint(tout, max_ci) << "\n";); + // TODO, we can in reality replay all constraints that don't depend on max_ci + vector replay; + unsigned i = 0; + for (auto ci : external) { + if (ci > max_ci) + replay.push_back(m_constraints[ci]); + else + external[i++] = ci; + } + external.shrink(i); + auto [p, k, j] = m_constraints[max_ci]; + while (m_constraints.size() > max_ci) + m_trail.pop_scope(1); + + for (auto const &[_p, _k, _j] : replay) { + auto ci = add_constraint(_p, _k, _j); + init_occurs(ci); + external.push_back(ci); + } + assumptions.erase(max_ci); + external.append(assumptions); + propagation_justification prop{external}; + auto new_ci = add_constraint(p, negate(k), prop); + TRACE(arith, display_constraint(tout << "backtrack ", new_ci) << "\n"); + init_occurs(new_ci); + return true; + } + + bool stellensatz::core_is_linear(svector const &core) { + m_constraints_in_conflict.reset(); + svector external, assumptions; + for (auto ci : core) + explain_constraint(ci, external, assumptions); + return all_of(assumptions, [&](auto ci) { return has_term_offset(m_constraints[ci].p); }); + } + + void stellensatz::explain_constraint(lp::constraint_index ci, svector &external, + svector &assumptions) { + if (m_constraints_in_conflict.contains(ci)) + return; + m_constraints_in_conflict.insert(ci); + auto &[p, k, b] = m_constraints[ci]; + if (std::holds_alternative(b)) { + external.push_back(ci); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.left, external, assumptions); + explain_constraint(m.right, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.left, external, assumptions); + explain_constraint(m.right, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + explain_constraint(m.divisor, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + explain_constraint(m.ci_eq, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + for (auto c : m.cs) + explain_constraint(c, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.left, external, assumptions); + explain_constraint(m.right, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + } + else if (std::holds_alternative(b)) { + assumptions.push_back(ci); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + } + else + UNREACHABLE(); + } + + // + // a constraint can be explained by a set of bounds used as justifications + // and by an original constraint. + // + void stellensatz::explain_constraint(lemma_builder &new_lemma, lp::constraint_index ci, lp::explanation &ex) { + svector external, assumptions; + explain_constraint(ci, external, assumptions); + for (auto ci : external) { + auto &[p, k, b] = m_constraints[ci]; + auto dep = std::get(b); + c().lra_solver().push_explanation(dep.dep, ex); + } + for (auto ci : assumptions) { + auto &[p, k, b] = m_constraints[ci]; + auto [t, coeff] = to_term_offset(p); + new_lemma |= ineq(t, negate(k), -coeff); + } + } + + rational stellensatz::value(dd::pdd const& p) const { + dd::pdd_eval eval; + eval.var2val() = [&](unsigned v) -> rational { return value(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 = is_int(p); + for (auto const& [c, is_constant] : p.coefficients()) + if (!is_constant || !_is_int) + g = gcd(g, c); + if (g == 0 || g == 1) + return ci; + switch (k) { + 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: + case lp::lconstraint_kind::LE: + UNREACHABLE(); + case lp::lconstraint_kind::EQ: + case lp::lconstraint_kind::NE: + 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::assume(dd::pdd const& p, lp::lconstraint_kind k) { + if (k == lp::lconstraint_kind::EQ) { + auto left = assume(p, lp::lconstraint_kind::GE); + auto right = assume(-p, lp::lconstraint_kind::GE); + return add_constraint(p, k, eq_justification{left, right}); + } + u_dependency *d = nullptr; + auto has_bound = [&](rational const& a, lpvar x, rational const& b) { + rational bound; + bool is_strict = false; + if (a == 1 && k == lp::lconstraint_kind::GE && c().lra_solver().has_lower_bound(x, d, bound, is_strict) && + bound >= -b) { + return true; + } + if (a == 1 && k == lp::lconstraint_kind::GT && c().lra_solver().has_lower_bound(x, d, bound, is_strict) && + (bound > -b || (is_strict && bound >= -b))) { + return true; + } + if (a == -1 && k == lp::lconstraint_kind::GE && c().lra_solver().has_upper_bound(x, d, bound, is_strict) && + bound <= -b) { + return true; + } + if (a == -1 && k == lp::lconstraint_kind::GT && c().lra_solver().has_upper_bound(x, d, bound, is_strict) && + (bound < -b || (is_strict && bound <= -b))) { + return true; + } + return false; + }; + + if (p.is_unilinear() && has_bound(p.hi().val(), p.var(), p.lo().val())) + // ax + b ~k~ 0 + return add_constraint(p, k, external_justification(d)); + return add_constraint(p, k, assumption_justification()); + } + + // p1 >= lo1, p2 >= lo2 => p1 + p2 >= lo1 + lo2 + lp::constraint_index stellensatz::add(lp::constraint_index left, lp::constraint_index right) { + 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 >= 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]; + 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) { + 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}); + } + + // p k 0, d > 0 -> p/d k 0, where q := d / d + // q is the quotient obtained by dividing the divisor constraint, which is of the form d - 1 >= 0 or d > 0 + lp::constraint_index stellensatz::divide(lp::constraint_index ci, lp::constraint_index divisor, dd::pdd q) { + auto const &[p, k, j] = m_constraints[ci]; + return add_constraint(q, k, division_justification{ci, divisor}); + } + + lp::constraint_index stellensatz::substitute(lp::constraint_index ci, lp::constraint_index ci_eq, lpvar v, + dd::pdd p) { + auto const &[p1, k1, j1] = m_constraints[ci]; + auto const &[p2, k2, j2] = m_constraints[ci_eq]; + SASSERT(k2 == lp::lconstraint_kind::EQ); + auto q = p1.subst_pdd(v, p); + return add_constraint(q, k1, substitute_justification{ci, ci_eq, v, p}); + } + + void stellensatz::init_occurs() { + m_occurs.reset(); + m_occurs.reserve(c().lra_solver().number_of_vars()); + 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; + if (m_has_occurs[ci]) + return; + struct undo_occurs : public trail { + stellensatz & s; + lp::constraint_index ci; + undo_occurs(stellensatz & s, lp::constraint_index ci) : s(s), ci(ci) {} + void undo() override { + s.remove_occurs(ci); + } + }; + add_active(ci); + m_trail.push(undo_occurs(*this, ci)); + m_has_occurs[ci] = true; + m_last_constraint = ci; + auto const &con = m_constraints[ci]; + for (auto v : con.p.free_vars()) + m_occurs[v].push_back(ci); + + } + + void stellensatz::remove_occurs(lp::constraint_index ci) { + SASSERT(m_has_occurs[ci]); + m_has_occurs[ci] = false; + auto const &con = m_constraints[ci]; + for (auto v : con.p.free_vars()) + m_occurs[v].pop_back(); + } + + bool stellensatz::is_int(svector const& vars) const { + return all_of(vars, [&](lpvar v) { return var_is_int(v); }); + } + + bool stellensatz::is_int(dd::pdd const &p) const { + return is_int(p.free_vars()); + } + + bool stellensatz::var_is_int(lp::lpvar v) const { + return c().lra_solver().var_is_int(v); + } + + bool stellensatz::constraint_is_true(lp::constraint_index ci) const { + return ci != lp::null_ci && constraint_is_true(m_constraints[ci]); + } + + bool stellensatz::constraint_is_false(lp::constraint_index ci) const { + return ci != lp::null_ci && !constraint_is_true(m_constraints[ci]); + } + + bool stellensatz::constraint_is_true(constraint const &c) const { + auto const& [p, k, j] = c; + auto lhs = value(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; + case lp::lconstraint_kind::LT: return lhs < 0; + case lp::lconstraint_kind::EQ: return lhs == 0; + case lp::lconstraint_kind::NE: return lhs != 0; + default: UNREACHABLE(); return false; + } + } + + bool stellensatz::constraint_is_trivial(lp::constraint_index ci) const { + return ci != lp::null_ci && constraint_is_trivial(m_constraints[ci]); + } + + bool stellensatz::constraint_is_trivial(constraint const &c) const { + auto const &[p, k, j] = c; + if (!p.is_val()) + return false; + if (p.val() > 0) + return k == lp::lconstraint_kind::GE || k == lp::lconstraint_kind::GT || k == lp::lconstraint_kind::NE; + else if (p.val() == 0) + return k == lp::lconstraint_kind::EQ || k == lp::lconstraint_kind::GE || k == lp::lconstraint_kind::LE; + else // p.val() < 0 + return k == lp::lconstraint_kind::LE || k == lp::lconstraint_kind::LT || k == lp::lconstraint_kind::NE; + } + + bool stellensatz::constraint_is_conflict(lp::constraint_index ci) const { + return ci != lp::null_ci && constraint_is_conflict(m_constraints[ci]); + } + + bool stellensatz::constraint_is_conflict(constraint const& c) const { + auto const &[p, k, j] = c; + return p.is_val() && + ((p.val() < 0 && k == lp::lconstraint_kind::GE) + || (p.val() <= 0 && k == lp::lconstraint_kind::GT)); + } + + std::ostream& stellensatz::display(std::ostream& out) const { + #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) { + display_constraint(out, ci) << "\n"; + out << "\t<- "; + display(out, m_constraints[ci].j); + out << "\n"; + } + return out; + } + + std::ostream& stellensatz::display_product(std::ostream& out, svector const& vars) const { + bool first = true; + for (auto v : vars) { + if (first) + first = false; + else + out << " * "; + out << "j" << v; + } + return out; + } + + std::ostream& stellensatz::display(std::ostream& out, term_offset const& t) const { + bool first = true; + for (auto [v, coeff] : t.first) { + c().display_coeff(out, first, coeff); + first = false; + out << pp_j(*this, v); + } + if (t.second != 0) + out << " + " << t.second; + return out; + } + + std::ostream& stellensatz::display_var(std::ostream& out, lpvar j) const { + 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 { + if (ci == lp::null_ci) + return out << "(null) "; + bool is_true = constraint_is_true(ci); + auto const& [p, k, j] = m_constraints[ci]; + return display_constraint(out << "(" << ci << ") ", m_constraints[ci]) + << (is_true ? " [true] " : " [false] ") << "(" << value(p) << " " << k << " 0)"; + } + + 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"; + } + + std::ostream &stellensatz::display(std::ostream &out, justification const &j) const { + if (std::holds_alternative(j)) { + auto dep = std::get(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(j)) { + auto m = std::get(j); + out << "(" << m.left << ") + (" << m.right << ")"; + } + else if (std::holds_alternative(j)) { + auto &m = std::get(j); + out << "(" << m.left << ") & (" << m.right << ")"; + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + out << "(" << m.ci << ") (" << m.ci_eq << ") by j" << m.v << " := " << m.p; + } + else if (std::holds_alternative(j)) { + auto &m = std::get(j); + out << "propagate "; + for (auto c : m.cs) + out << "(" << c << ") "; + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + out << "(" << m.left << ") * (" << m.right << ")"; + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + out << m.p << " * (" << m.ci << ")"; + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + out << "(" << m.ci << ") / (" << m.divisor << ")"; + } + else if (std::holds_alternative(j)) { + out << "assumption"; + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + out << " gcd (" << m.ci << ")"; + } + else + UNREACHABLE(); + return out; + } + + std::ostream &stellensatz::display_lemma(std::ostream &out, lp::explanation const &ex) { + m_constraints_in_conflict.reset(); + svector todo; + for (auto c : ex) + todo.push_back(c.ci()); + for (unsigned i = 0; i < todo.size(); ++i) { + auto ci = todo[i]; + if (m_constraints_in_conflict.contains(ci)) + continue; + m_constraints_in_conflict.insert(ci); + out << "(" << ci << ") "; + display_constraint(out, ci) << " "; + auto const& j = m_constraints[ci].j; + + display(out, j) << "\n"; + if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.left); + todo.push_back(m.right); + } + if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.left); + todo.push_back(m.right); + } + if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.append(m.cs); + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.ci); + todo.push_back(m.ci_eq); + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.ci); + todo.push_back(m.divisor); + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.left); + todo.push_back(m.right); + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.ci); + } + else if (std::holds_alternative(j)) { + // do nothing + } + else if (std::holds_alternative(j)) { + // do nothing + } + else if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.ci); + } + else + NOT_IMPLEMENTED_YET(); + + } + return out; + } + + void stellensatz::updt_params(params_ref const& p) { + smt_params_helper sp(p); + m_config.max_degree = sp.arith_nl_stellensatz_max_degree(); + m_config.max_conflicts = sp.arith_nl_stellensatz_max_conflicts(); + m_config.max_constraints = sp.arith_nl_stellensatz_max_constraints(); + m_config.strategy = sp.arith_nl_stellensatz_strategy(); + } + + + // Solver component + // check LRA feasibilty + // (partial) check LIA feasibility + // try patch LIA model + // option: iterate by continuing saturation based on partial results + // option: add incremental linearization axioms + // 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_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; + } + + lbool stellensatz::solver::solve(svector& core) { + init(); + lbool r = solve_lra(); + if (r == l_true) + r = solve_lia(); + + if (r == l_false) { + core.reset(); + for (auto p : m_ex) + core.push_back(m_internal2external_constraints[p.ci()]); + } + return r; + } + + lbool stellensatz::solver::solve_lra() { + auto status = lra_solver->find_feasible_solution();; + if (lra_solver->is_feasible()) + return l_true; + if (status == lp::lp_status::INFEASIBLE) { + lra_solver->get_infeasibility_explanation(m_ex); + return l_false; + } + return l_undef; + } + + 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; + default: // TODO: an option is to perform (bounded) search here to get an LIA verdict. + return l_undef; + } + return l_undef; + } + + // update values using the model + void stellensatz::solver::update_values(vector& values) { + std::unordered_map _values; + lra_solver->get_model(_values); + unsigned sz = values.size(); + for (unsigned i = 0; i < sz; ++i) + values[i] = _values[i]; + TRACE(arith, for (unsigned i = 0; i < sz; ++i) tout << "j" << i << " := " << values[i] << "\n";); + } +} diff --git a/src/math/lp/nla_stellensatz.h b/src/math/lp/nla_stellensatz.h index 0583c95b4..b6dd27e24 100644 --- a/src/math/lp/nla_stellensatz.h +++ b/src/math/lp/nla_stellensatz.h @@ -5,7 +5,9 @@ #pragma once #include "util/scoped_ptr_vector.h" +#include "util/dependency.h" #include "math/dd/dd_pdd.h" +#include "math/interval/dep_intervals.h" #include "math/lp/nla_coi.h" #include "math/lp/int_solver.h" #include @@ -131,23 +133,6 @@ namespace nla { unsigned strategy = 0; }; - - - trail_stack m_trail; - coi m_coi; - dd::pdd_manager pddm; - config m_config; - vector m_constraints; - monomial_factory m_monomial_factory; - indexed_uint_set m_active, m_processed; - vector m_values; - svector m_core; - vector> m_occurs; // map from variable to constraints they occur. - bool_vector m_has_occurs; - lp::constraint_index m_last_constraint = lp::null_ci; - - unsigned_vector m_var2level, m_level2var; - struct constraint_key { unsigned pdd; lp::lconstraint_kind k; @@ -159,28 +144,171 @@ namespace nla { struct hash { unsigned operator()(constraint_key const &c) const { return hash_u_u(c.pdd, static_cast(c.k)); - } + } }; }; + + struct bound_info { + unsigned m_var; + lp::lconstraint_kind m_kind; + rational m_value; + unsigned m_level = 0; + unsigned m_last_bound = UINT_MAX; + bool m_is_decision = true; + u_dependency* m_bound_justifications = nullptr; // index into bounds + lp::constraint_index m_constraint_justification = lp::null_ci; + bound_info(lpvar v, lp::lconstraint_kind k, rational const &value, unsigned level, unsigned last_bound, u_dependency *d, + lp::constraint_index ci) + : m_var(v), m_kind(k), m_value(value), m_level(level), m_last_bound(last_bound), m_is_decision(false), + m_bound_justifications(d), m_constraint_justification(ci) {} + bound_info(lpvar v, lp::lconstraint_kind k, rational const &value, unsigned level, unsigned last_bound) + : m_var(v), m_kind(k), m_value(value), m_level(level), m_last_bound(last_bound), m_is_decision(true) {} + }; + + struct assignment { + lpvar m_var; + bool m_is_upper; + }; + + trail_stack m_ctrail, m_vtrail; // constraint and variable trail + coi m_coi; + dd::pdd_manager pddm; + config m_config; + vector m_constraints; + monomial_factory m_monomial_factory; + vector m_values; + svector m_core; + vector> m_occurs; // map from variable to constraints they occur. + bool_vector m_has_occurs; // is the constraint indexed already map m_constraint_index; + // + // there is a default interpretation of variables. + // The default interpretation has the invariant that the values are within the asserted bounds of variables. + // Every time the default interpretation changes, the set of false constraints is updated. + // A false constraint is conflicting if the lower and upper bounds render the constraint unfixable. + // The solver enters the conflict stage and performs backjumping. + // If none of the constraints are unfixable, then + // propagation loops over false constraints and performs bounds propagation if it merits fixes. + // if there is no propagation, perform a decision to fix the default interpretation of a constraint. + // the side effect of decisions and propagations is that the set of false constraints are updated. + // + // Extensions: + // - we can assign priorities to variables to choose one variable to fix over another when repairing a + // constraint that is false. + // - we can incorporate backtracking instead of backjumping by replaying propagations + // + unsigned_vector m_lower, m_upper; // var -> index into m_bounds + vector m_bounds; // bound index -> bound meta-data + unsigned m_level = 0; + + unsigned m_prop_qhead = 0; // head into propagation queue + lp::constraint_index m_conflict = lp::null_ci; + u_dependency *m_conflict_dep = nullptr; + + u_dependency_manager m_dm; + dep_intervals m_di; + indexed_uint_set m_conflict_marked; + + void propagate(); + bool decide(); + lbool search(); + lbool resolve_conflict(); + void init_search(); + void pop_bound(); + void mark_dependencies(u_dependency *d); + bool should_propagate() const { return m_prop_qhead < m_bounds.size(); } + // assuming variables have bounds determine if polynomial has lower/upper bounds + void interval(dd::pdd p, scoped_dep_interval &iv); + + void set_conflict(lp::constraint_index ci, u_dependency *d) { + m_conflict = ci; + m_conflict_dep = d; + } + void set_conflict(lpvar v) { + m_conflict_dep = m_dm.mk_join(lo_dep(v), hi_dep(v)); + m_conflict = resolve_variable(v, lo_constraint(v), hi_constraint(v)); + } + void reset_conflict() { m_conflict = lp::null_ci; } + bool is_conflict() const { return m_conflict != lp::null_ci; } + + indexed_uint_set m_processed; + unsigned_vector m_var2level, m_level2var; + bool has_lo(lpvar v) const { + return m_lower[v] != UINT_MAX; + } + bool has_hi(lpvar v) const { + return m_upper[v] != UINT_MAX; + } + rational const& lo_val(lpvar v) const { + SASSERT(has_lo(v)); + return m_bounds[m_lower[v]].m_value; + } + rational const& hi_val(lpvar v) const { + SASSERT(has_hi(v)); + return m_bounds[m_upper[v]].m_value; + } + lp::lconstraint_kind lo_kind(lpvar v) const { + SASSERT(has_lo(v)); + return m_bounds[m_lower[v]].m_kind; + } + lp::lconstraint_kind hi_kind(lpvar v) const { + SASSERT(has_hi(v)); + return m_bounds[m_upper[v]].m_kind; + } + bool lo_is_strict(lpvar v) const { + SASSERT(has_lo(v)); + return lo_kind(v) == lp::lconstraint_kind::GT; + } + bool hi_is_strict(lpvar v) const { + SASSERT(has_hi(v)); + return hi_kind(v) == lp::lconstraint_kind::LT; + } + u_dependency *lo_dep(lpvar v) const { + SASSERT(has_lo(v)); + return m_bounds[m_lower[v]].m_bound_justifications; + } + u_dependency *hi_dep(lpvar v) const { + SASSERT(has_hi(v)); + return m_bounds[m_upper[v]].m_bound_justifications; + } + lp::constraint_index lo_constraint(lpvar v) const { return m_bounds[m_lower[v]].m_constraint_justification; } + lp::constraint_index hi_constraint(lpvar v) const { return m_bounds[m_upper[v]].m_constraint_justification; } + bool is_fixed(lpvar v) const { return has_lo(v) && has_hi(v) && lo_val(v) == hi_val(v); } + void move_up(lpvar v); + + + + struct repair_var_info { + rational lo_value, hi_value; + lp::constraint_index lo = lp::null_ci, hi = lp::null_ci, vanishing = lp::null_ci; + }; + repair_var_info find_bounds(lpvar v); + unsigned max_level(constraint const &c) const; + bool repair_variable(lpvar& v, rational &r, lp::lconstraint_kind& k, lp::constraint_index& ci); + void find_split(lpvar& v, rational& r, lp::lconstraint_kind& k, lp::constraint_index ci); + void set_in_bounds(lpvar v); + bool in_bounds(lpvar v) { return in_bounds(v, m_values[v]); } + + + + bool in_bounds(lpvar v, rational const &value) const; dd::pdd to_pdd(lpvar v); void init_monomial(unsigned mon_var); term_offset to_term_offset(dd::pdd const &p); bool has_term_offset(dd::pdd const &p); - lp::constraint_index add_constraint(dd::pdd p, lp::lconstraint_kind k, justification j); - + 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); - // initialization void init_solver(); void init_vars(); void init_occurs(); void init_occurs(lp::constraint_index ci); + void init_bounds(); void pop_constraint(); void remove_occurs(lp::constraint_index ci); @@ -188,30 +316,23 @@ namespace nla { bool conflict(lp::constraint_index ci); void conflict(svector const& core); lp::constraint_index vanishing(lpvar x, factorization const& f, lp::constraint_index ci); - 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); - lp::constraint_index resolve_variable(lpvar x, lp::constraint_index ci, lp::constraint_index other_ci, rational const& p_value, - factorization const& f, unsigned_vector const& m1, dd::pdd _f_p); - - lbool model_repair(); - lbool model_repair_loop(); - lp::constraint_index model_repair(lp::lpvar v); - struct bound_info { - rational lo_value, hi_value; - lp::constraint_index lo = lp::null_ci, hi = lp::null_ci, vanishing = lp::null_ci; - }; - bound_info find_bounds(lpvar v); - unsigned max_level(constraint const &c) const; + lp::constraint_index resolve_variable(lpvar x, lp::constraint_index ci, lp::constraint_index other_ci); bool constraint_is_true(lp::constraint_index ci) const; - bool constraint_is_false(lp::constraint_index ci) const; bool constraint_is_true(constraint const &c) const; + bool constraint_is_false(lp::constraint_index ci) const; + bool constraint_is_false(constraint const &c) const { return !constraint_is_true(c); } bool constraint_is_conflict(lp::constraint_index ci) const; bool constraint_is_conflict(constraint const &c) const; bool constraint_is_trivial(lp::constraint_index ci) const; bool constraint_is_trivial(constraint const& c) const; - bool is_new_constraint(lp::constraint_index ci) const; + bool constraint_is_bound_conflict(constraint const &c, u_dependency*& d); + bool constraint_is_bound_conflict(lp::constraint_index ci, u_dependency*& d) { return constraint_is_bound_conflict(m_constraints[ci], d); } + bool var_is_bound_conflict(lpvar v) const; + + bool constraint_is_propagating(lp::constraint_index ci, u_dependency *&d, lpvar &v, rational &value, lp::lconstraint_kind& k); lp::constraint_index gcd_normalize(lp::constraint_index ci); lp::constraint_index assume(dd::pdd const& p, lp::lconstraint_kind k); @@ -226,10 +347,9 @@ namespace nla { bool var_is_int(lp::lpvar v) const; rational value(dd::pdd const& p) const; rational value(lp::lpvar v) const { return m_values[v]; } + unsigned num_vars() const { return m_values.size(); } bool set_model(); - void add_active(lp::constraint_index ci); - // lemmas indexed_uint_set m_constraints_in_conflict; void explain_constraint(lemma_builder& new_lemma, lp::constraint_index ci, lp::explanation &ex); @@ -238,6 +358,9 @@ namespace nla { bool backtrack(svector const& core); bool core_is_linear(svector const &core); + bool well_formed() const; + bool well_formed(lpvar v) const; + struct pp_j { stellensatz const &s; lpvar j; @@ -253,6 +376,7 @@ namespace nla { std::ostream& display_product(std::ostream& out, svector const& vars) const; std::ostream& display_constraint(std::ostream& out, lp::constraint_index ci) const; std::ostream& display_constraint(std::ostream& out, constraint const& c) const; + std::ostream &display_bound(std::ostream &out, unsigned bound_index) 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); diff --git a/src/math/lp/nla_stellensatz.h.2 b/src/math/lp/nla_stellensatz.h.2 new file mode 100644 index 000000000..0583c95b4 --- /dev/null +++ b/src/math/lp/nla_stellensatz.h.2 @@ -0,0 +1,268 @@ +/*++ + Copyright (c) 2025 Microsoft Corporation + + --*/ +#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 + +namespace nla { + + class core; + class lar_solver; + class stellensatz : common { + + struct external_justification { + u_dependency *dep = nullptr; + external_justification(u_dependency *d) : dep(d) {} + }; + struct assumption_justification {}; + struct addition_justification { + lp::constraint_index left, right; + }; + struct multiplication_poly_justification { + lp::constraint_index ci; + dd::pdd p; + }; + struct multiplication_justification { + lp::constraint_index left, right; + }; + struct division_justification { + lp::constraint_index ci; + lp::constraint_index divisor; + }; + struct substitute_justification { + lp::constraint_index ci; + lp::constraint_index ci_eq; + lpvar v; + dd::pdd p; + }; + struct eq_justification { + lp::constraint_index left; + lp::constraint_index right; + }; + struct gcd_justification { + lp::constraint_index ci; + }; + struct propagation_justification { + svector cs; + }; + + using justification = std::variant< + external_justification, + assumption_justification, + gcd_justification, + substitute_justification, + addition_justification, + division_justification, + eq_justification, + propagation_justification, + multiplication_poly_justification, + multiplication_justification>; + + struct constraint { + dd::pdd p; + lp::lconstraint_kind k; + justification j; + }; + + class monomial_factory { + struct eq { + bool operator()(unsigned_vector const &a, unsigned_vector const &b) const { + return a == b; + } + }; + map, eq> m_vars2mon; + u_map 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 const &vars); + + lpvar get_monomial(svector const &vars); + + void init(lpvar v, svector const &_vars); + }; + + using term_offset = std::pair; + + class solver { + stellensatz &s; + scoped_ptr lra_solver; + scoped_ptr int_solver; + lp::explanation m_ex; + unsigned_vector m_internal2external_constraints; + monomial_factory m_monomial_factory; + lbool solve_lra(); + lbool solve_lia(); + + void init(); + term_offset to_term_offset(dd::pdd const &p); + public: + solver(stellensatz &s) : s(s) {}; + lbool solve(svector& core); + void update_values(vector& values); + }; + + 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 config { + unsigned max_degree = 3; + unsigned max_conflicts = UINT_MAX; + unsigned max_constraints = UINT_MAX; + unsigned strategy = 0; + }; + + + + trail_stack m_trail; + coi m_coi; + dd::pdd_manager pddm; + config m_config; + vector m_constraints; + monomial_factory m_monomial_factory; + indexed_uint_set m_active, m_processed; + vector m_values; + svector m_core; + vector> m_occurs; // map from variable to constraints they occur. + bool_vector m_has_occurs; + lp::constraint_index m_last_constraint = lp::null_ci; + + unsigned_vector m_var2level, m_level2var; + + struct constraint_key { + unsigned pdd; + lp::lconstraint_kind k; + struct eq { + bool operator()(constraint_key const &a, constraint_key const &b) const { + return a.pdd == b.pdd && a.k == b.k; + } + }; + struct hash { + unsigned operator()(constraint_key const &c) const { + return hash_u_u(c.pdd, static_cast(c.k)); + } + }; + }; + map m_constraint_index; + + + dd::pdd to_pdd(lpvar v); + void init_monomial(unsigned mon_var); + term_offset to_term_offset(dd::pdd const &p); + bool has_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); + + + + // initialization + void init_solver(); + void init_vars(); + void init_occurs(); + void init_occurs(lp::constraint_index ci); + void pop_constraint(); + void remove_occurs(lp::constraint_index ci); + + lp::constraint_index factor(lp::constraint_index ci); + bool conflict(lp::constraint_index ci); + void conflict(svector const& core); + lp::constraint_index vanishing(lpvar x, factorization const& f, lp::constraint_index ci); + 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); + lp::constraint_index resolve_variable(lpvar x, lp::constraint_index ci, lp::constraint_index other_ci, rational const& p_value, + factorization const& f, unsigned_vector const& m1, dd::pdd _f_p); + + lbool model_repair(); + lbool model_repair_loop(); + lp::constraint_index model_repair(lp::lpvar v); + struct bound_info { + rational lo_value, hi_value; + lp::constraint_index lo = lp::null_ci, hi = lp::null_ci, vanishing = lp::null_ci; + }; + bound_info find_bounds(lpvar v); + unsigned max_level(constraint const &c) const; + + bool constraint_is_true(lp::constraint_index ci) const; + bool constraint_is_false(lp::constraint_index ci) const; + bool constraint_is_true(constraint const &c) const; + bool constraint_is_conflict(lp::constraint_index ci) const; + bool constraint_is_conflict(constraint const &c) const; + bool constraint_is_trivial(lp::constraint_index ci) const; + bool constraint_is_trivial(constraint const& c) const; + bool is_new_constraint(lp::constraint_index ci) const; + + 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, dd::pdd p); + lp::constraint_index multiply(lp::constraint_index left, lp::constraint_index right); + lp::constraint_index divide(lp::constraint_index ci, lp::constraint_index divisor, dd::pdd d); + lp::constraint_index substitute(lp::constraint_index ci, lp::constraint_index ci_eq, lpvar v, dd::pdd p); + + bool is_int(svector const& vars) const; + bool is_int(dd::pdd const &p) const; + bool var_is_int(lp::lpvar v) const; + rational value(dd::pdd const& p) const; + rational value(lp::lpvar v) const { return m_values[v]; } + bool set_model(); + + void add_active(lp::constraint_index ci); + + // lemmas + indexed_uint_set m_constraints_in_conflict; + void explain_constraint(lemma_builder& new_lemma, lp::constraint_index ci, lp::explanation &ex); + void explain_constraint(lp::constraint_index ci, svector &external, + svector &assumptions); + bool backtrack(svector const& core); + bool core_is_linear(svector const &core); + + struct pp_j { + stellensatz const &s; + lpvar j; + pp_j(stellensatz const&s, lpvar j) : s(s), j(j) {} + std::ostream &display(std::ostream &out) const { + return s.display_var(out, j); + } + }; + friend std::ostream &operator<<(std::ostream &out, pp_j const &p) { + return p.display(out); + } + std::ostream& display(std::ostream& out) const; + std::ostream& display_product(std::ostream& out, svector const& vars) const; + std::ostream& display_constraint(std::ostream& out, lp::constraint_index ci) 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); + std::ostream &display(std::ostream &out, term_offset const &t) const; + + public: + stellensatz(core* core); + lbool saturate(); + + void updt_params(params_ref const &p); + }; + +}