mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2026-01-13 14:16:14 +00:00
Code Activity

introduce isolate_root_closest

Browse source 
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2026-01-10 19:06:52 -10:00
parent 163977506c
commit 66ef878162
3 changed files with 512 additions and 48 deletions
Download patch file Download diff file Expand all files Collapse all files

238
src/math/polynomial/algebraic_numbers.cpp
View file

@ -678,6 +678,169 @@ namespace algebraic_numbers {
isolate_roots(up, roots);
}
unsigned sign_variations_at_mpq(upolynomial::upolynomial_sequence const & seq, mpq const & b) {
unsigned sz = seq.size();
if (sz <= 1)
return 0;
unsigned r = 0;
int sign = 0, prev_sign = 0;
for (unsigned i = 0; i < sz; ++i) {
unsigned psz = seq.size(i);
mpz const * p = seq.coeffs(i);
sign = static_cast<int>(upm().eval_sign_at(psz, p, b));
if (sign == 0)
continue;
if (prev_sign != 0 && sign != prev_sign)
r++;
prev_sign = sign;
}
return r;
}
// Isolate the i-th real root of sqf_p (1-based), where sqf_p is square-free and univariate.
// Return the root as an algebraic number in r. The polynomial stored in the result is sqf_p.
void isolate_kth_root(scoped_upoly const & sqf_p, upolynomial::upolynomial_sequence const & seq, unsigned i, numeral & r) {
SASSERT(i > 0);
unsigned sz = sqf_p.size();
mpz const * p = sqf_p.data();
SASSERT(sz > 0);
if (sz == 2) {
// Linear polynomial ax + b: root is -b/a (always rational).
scoped_mpq q(qm());
qm().set(q, p[0], p[1]);
qm().neg(q);
set(r, q);
return;
}
unsigned pos_k = upm().knuth_positive_root_upper_bound(sz, p);
unsigned neg_k = upm().knuth_negative_root_upper_bound(sz, p);
scoped_mpbq lo(bqm()), hi(bqm()), mid(bqm());
unsigned vminus = upm().sign_variations_at_minus_inf(seq);
unsigned v0 = upm().sign_variations_at_zero(seq);
unsigned vplus = upm().sign_variations_at_plus_inf(seq);
unsigned le0_cnt = vminus - v0; // roots in (-oo, 0]
unsigned vlo, vhi;
unsigned target;
if (i <= le0_cnt) {
// Isolate within (-2^neg_k, 0] to keep the interval on the non-positive side.
bqm().power(mpbq(2), neg_k, lo);
bqm().neg(lo);
bqm().set(hi, 0);
vlo = vminus;
vhi = v0;
target = i;
}
else {
// Isolate within (0, 2^pos_k] to keep the interval on the non-negative side.
bqm().set(lo, 0);
bqm().power(mpbq(2), pos_k, hi);
vlo = v0;
vhi = vplus;
target = i - le0_cnt;
}
// Sanity: sqf_p has at least i roots.
SASSERT(vlo >= vhi);
SASSERT(i <= vminus - vplus);
SASSERT(target > 0);
SASSERT(target <= vlo - vhi);
while (vlo > vhi + 1) {
checkpoint();
bqm().add(lo, hi, mid);
bqm().div2(mid);
unsigned vmid = upm().sign_variations_at(seq, mid);
unsigned left_cnt = vlo - vmid; // roots in (lo, mid]
if (target <= left_cnt) {
bqm().set(hi, mid);
vhi = vmid;
}
else {
bqm().set(lo, mid);
vlo = vmid;
target -= left_cnt;
}
}
SASSERT(vlo == vhi + 1);
// If the upper endpoint is exactly a dyadic root, return it as a basic number.
if (upm().eval_sign_at(sz, p, hi) == 0) {
scoped_mpq q(qm());
to_mpq(qm(), hi, q);
set(r, q);
return;
}
// Convert the isolating interval into a refinable one (or discover a dyadic root on the way).
scoped_mpbq a(bqm()), b(bqm());
bqm().set(a, lo);
bqm().set(b, hi);
if (!upm().isolating2refinable(sz, p, bqm(), a, b)) {
scoped_mpq q(qm());
to_mpq(qm(), a, q);
set(r, q);
return;
}
del(r);
r = mk_algebraic_cell(sz, p, a, b, false /* minimal */);
SASSERT(acell_inv(*r.to_algebraic()));
}
// Closest-root isolation for an (integer) univariate polynomial.
void isolate_roots_closest_univariate(polynomial_ref const & p, mpq const & s, numeral_vector & roots, svector<unsigned> & indices) {
SASSERT(is_univariate(p));
SASSERT(roots.empty());
indices.reset();
if (::is_zero(p) || ::is_const(p))
return;
// Convert to dense univariate form and take the square-free part.
scoped_upoly & up = m_isolate_tmp1;
scoped_upoly & sqf_p = m_isolate_tmp3;
up.reset();
sqf_p.reset();
upm().to_numeral_vector(p, up);
if (up.empty())
return;
upm().square_free(up.size(), up.data(), sqf_p);
if (sqf_p.empty() || upm().degree(sqf_p) == 0)
return;
upolynomial::scoped_upolynomial_sequence seq(upm());
upm().sturm_seq(sqf_p.size(), sqf_p.data(), seq);
unsigned vminus = upm().sign_variations_at_minus_inf(seq);
unsigned vplus = upm().sign_variations_at_plus_inf(seq);
if (vminus <= vplus)
return;
unsigned vs = sign_variations_at_mpq(seq, s);
unsigned total = vminus - vplus;
unsigned k = vminus - vs; // #roots in (-oo, s]
if (upm().eval_sign_at(sqf_p.size(), sqf_p.data(), s) == 0) {
roots.push_back(numeral());
set(roots.back(), s);
indices.push_back(k);
return;
}
// predecessor (<= s)
if (k > 0) {
roots.push_back(numeral());
isolate_kth_root(sqf_p, seq, k, roots.back());
indices.push_back(k);
}
// successor (> s)
if (k < total) {
roots.push_back(numeral());
isolate_kth_root(sqf_p, seq, k + 1, roots.back());
indices.push_back(k + 1);
}
}
void mk_root(scoped_upoly const & up, unsigned i, numeral & r) {
// TODO: implement version that finds i-th root without isolating all roots.
if (i == 0)
@ -2547,6 +2710,77 @@ namespace algebraic_numbers {
}
}
void isolate_roots_closest(polynomial_ref const & p, polynomial::var2anum const & x2v, mpq const & s, numeral_vector & roots, svector<unsigned> & indices) {
TRACE(isolate_roots, tout << "isolating closest roots of: " << p << " around " << m_qmanager.to_string(s) << "\n";);
SASSERT(roots.empty());
indices.reset();
polynomial::manager & ext_pm = p.m();
if (ext_pm.is_zero(p) || ext_pm.is_const(p))
return;
if (ext_pm.is_univariate(p)) {
isolate_roots_closest_univariate(p, s, roots, indices);
return;
}
// eliminate rationals
polynomial_ref p_prime(ext_pm);
var2basic x2v_basic(*this, x2v);
p_prime = ext_pm.substitute(p, x2v_basic);
if (ext_pm.is_zero(p_prime) || ext_pm.is_const(p_prime))
return;
if (ext_pm.is_univariate(p_prime)) {
polynomial::var x = ext_pm.max_var(p_prime);
if (x2v.contains(x)) {
// The remaining variable is assigned, the actual unassigned variable vanished when we replaced rational values.
// So, the polynomial does not have any roots.
return;
}
isolate_roots_closest_univariate(p_prime, s, roots, indices);
return;
}
// Fallback: isolate all roots then select closest ones.
scoped_numeral_vector all(m_wrapper);
isolate_roots(p, x2v, all);
if (all.empty())
return;
scoped_numeral sv(m_wrapper);
set(sv, s);
unsigned lower = UINT_MAX, upper = UINT_MAX;
for (unsigned k = 0; k < all.size(); ++k) {
auto cmp = compare(all[k], sv);
if (cmp <= 0)
lower = k;
else {
upper = k;
break;
}
}
if (lower != UINT_MAX && eq(all[lower], s)) {
roots.push_back(numeral());
set(roots.back(), all[lower]);
indices.push_back(lower + 1);
return;
}
if (lower != UINT_MAX) {
roots.push_back(numeral());
set(roots.back(), all[lower]);
indices.push_back(lower + 1);
}
if (upper != UINT_MAX) {
roots.push_back(numeral());
set(roots.back(), all[upper]);
indices.push_back(upper + 1);
}
}
sign eval_at_mpbq(polynomial_ref const & p, polynomial::var2anum const & x2v, polynomial::var x, mpbq const & v) {
scoped_mpq qv(qm());
to_mpq(qm(), v, qv);
@ -3011,6 +3245,10 @@ namespace algebraic_numbers {
m_imp->isolate_roots(p, x2v, roots);
}
void manager::isolate_roots_closest(polynomial_ref const & p, polynomial::var2anum const & x2v, mpq const & s, numeral_vector & roots, svector<unsigned> & indices) {
m_imp->isolate_roots_closest(p, x2v, s, roots, indices);
}
void manager::isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots, svector<sign> & signs) {
m_imp->isolate_roots(p, x2v, roots, signs);
}

13
src/math/polynomial/algebraic_numbers.h
View file

@ -169,6 +169,19 @@ namespace algebraic_numbers {
*/
void isolate_roots(polynomial_ref const & p, polynomial::var2anum const & x2v, numeral_vector & roots);
/**
\brief Isolate the closest real roots of a multivariate polynomial p around the rational point s.
The method behaves like isolate_roots(p, x2v, roots) but only returns:
- the last root r such that r <= s (if it exists), and
- the first root r such that r > s (if it exists),
or a single root if s itself is a root.
The returned roots are sorted increasingly. The associated 1-based root indices
(with respect to the full increasing root list) are stored in \c indices.
*/
void isolate_roots_closest(polynomial_ref const & p, polynomial::var2anum const & x2v, mpq const & s, numeral_vector & roots, svector<unsigned> & indices);
/**
\brief Isolate the roots of the given polynomial, and compute its sign between them.
*/

309
src/nlsat/levelwise.cpp
View file

@ -5,8 +5,10 @@
#include "math/polynomial/polynomial.h"
#include "nlsat_common.h"
#include "util/z3_exception.h"
#include "util/vector.h"
#include <algorithm>
#include <cstdint>
#include <set>
#include <utility>
#include <vector>
@ -21,6 +23,18 @@ enum relation_mode {
lowest_degree
};
// Tag indicating what kind of invariance we need to preserve for a polynomial on the cell:
// - SIGN: sign-invariance is sufficient
// - ORD: order-invariance is required (a stronger requirement)
//
// This is inspired by SMT-RAT's InvarianceType and is used together with a
// de-dup/upgrade worklist discipline (appendOnCorrectLevel()).
enum class inv_req : uint8_t { none = 0, sign = 1, ord = 2 };
static inline inv_req max_req(inv_req a, inv_req b) {
return static_cast<uint8_t>(a) < static_cast<uint8_t>(b) ? b : a;
}
struct nullified_poly_exception {};
struct levelwise::impl {
@ -37,6 +51,9 @@ struct levelwise::impl {
polynomial_ref_vector m_psc_tmp; // scratch for PSC chains
bool m_fail = false;
// Current requirement tag for polynomials stored in the todo_set, keyed by pm.id(poly*).
// Entries are upgraded SIGN -> ORD as needed and cleared when a polynomial is extracted.
svector<uint8_t> m_req;
assignment const& sample() const { return m_solver.sample(); }
@ -175,39 +192,104 @@ struct levelwise::impl {
return polynomial_ref(m_pm);
}
void insert_factorized(todo_set& P, polynomial_ref const& poly) {
poly* request(todo_set& P, poly* p, inv_req req) {
if (!p || req == inv_req::none)
return p;
p = P.insert(p);
unsigned id = m_pm.id(p);
auto cur = static_cast<inv_req>(m_req.get(id, static_cast<uint8_t>(inv_req::none)));
inv_req nxt = max_req(cur, req);
if (nxt != cur)
m_req.setx(id, static_cast<uint8_t>(nxt), static_cast<uint8_t>(inv_req::none));
return p;
}
void request_factorized(todo_set& P, polynomial_ref const& poly, inv_req req) {
for_each_distinct_factor(poly, [&](polynomial_ref const& f) {
P.insert(f.get());
request(P, f.get(), req); // inherit tag across factorization (SMT-RAT appendOnCorrectLevel)
});
}
using tagged_id = std::pair<unsigned, inv_req>; // <pm.id(poly), tag>
var extract_max_tagged(todo_set& P, polynomial_ref_vector& max_ps, std::vector<tagged_id>& tagged) {
var level = P.extract_max_polys(max_ps);
tagged.clear();
tagged.reserve(max_ps.size());
for (unsigned i = 0; i < max_ps.size(); ++i) {
poly* p = max_ps.get(i);
unsigned id = m_pm.id(p);
inv_req req = static_cast<inv_req>(m_req.get(id, static_cast<uint8_t>(inv_req::sign)));
if (req == inv_req::none)
req = inv_req::sign;
tagged.emplace_back(id, req);
// Clear: extracted polynomials are no longer part of the worklist.
m_req.setx(id, static_cast<uint8_t>(inv_req::none), static_cast<uint8_t>(inv_req::none));
}
return level;
}
// Select a coefficient c of p (wrt x) such that c(s) != 0, or return null.
polynomial_ref choose_nonzero_coeff(polynomial_ref const& p, unsigned x) {
// The coefficients are polynomials in lower variables; we prefer "simpler" ones
// to reduce projection blow-up.
polynomial_ref choose_nonzero_coeff(polynomial_ref const& p, unsigned x, poly* exclude = nullptr) {
unsigned deg = m_pm.degree(p, x);
polynomial_ref coeff(m_pm);
for (int j = static_cast<int>(deg); j >= 0; --j) {
coeff = m_pm.coeff(p, x, static_cast<unsigned>(j));
if (!coeff || is_zero(coeff))
continue;
// First pass: any non-zero constant coefficient is ideal (no need to project it).
for (unsigned j = 0; j <= deg; ++j) {
coeff = m_pm.coeff(p, x, j);
if (!coeff || is_zero(coeff) || (exclude && coeff.get() == exclude))
continue;
if (!is_const(coeff))
continue;
if (m_am.eval_sign_at(coeff, sample()) != 0)
return coeff;
}
return polynomial_ref(m_pm);
// Second pass: pick the non-constant coefficient with minimal (total_degree, size, max_var).
polynomial_ref best(m_pm);
unsigned best_td = 0;
unsigned best_sz = 0;
var best_mv = null_var;
for (unsigned j = 0; j <= deg; ++j) {
coeff = m_pm.coeff(p, x, j);
if (!coeff || is_zero(coeff) || (exclude && coeff.get() == exclude))
continue;
if (m_am.eval_sign_at(coeff, sample()) == 0)
continue;
if (is_const(coeff))
continue; // handled above
unsigned td = total_degree(coeff);
unsigned sz = m_pm.size(coeff);
var mv = m_pm.max_var(coeff);
if (!best ||
td < best_td ||
(td == best_td && (sz < best_sz ||
(sz == best_sz && (best_mv == null_var || mv < best_mv))))) {
best = coeff;
best_td = td;
best_sz = sz;
best_mv = mv;
}
}
return best;
}
void add_projections_for(todo_set& P, polynomial_ref const& p, unsigned x, polynomial_ref const& nonzero_coeff) {
void add_projections_for(todo_set& P, polynomial_ref const& p, unsigned x, polynomial_ref const& nonzero_coeff, bool add_leading_coeff) {
// Line 11 (non-null witness coefficient)
if (nonzero_coeff && !is_const(nonzero_coeff))
insert_factorized(P, nonzero_coeff);
request_factorized(P, nonzero_coeff, inv_req::sign);
// Line 12 (disc + leading coefficient)
insert_factorized(P, psc_discriminant(p, x));
request_factorized(P, psc_discriminant(p, x), inv_req::ord);
unsigned deg = m_pm.degree(p, x);
polynomial_ref lc(m_pm);
lc = m_pm.coeff(p, x, deg);
insert_factorized(P, lc);
if (add_leading_coeff && lc.get() != nonzero_coeff.get())
request_factorized(P, lc, inv_req::sign);
}
// Relation construction heuristics (same intent as previous implementation).
@ -362,45 +444,118 @@ struct levelwise::impl {
}
// Build Θ (root functions), pick I_level around sample(level), and build relation ≼.
void build_interval_and_relation(unsigned level, polynomial_ref_vector const& level_ps) {
void build_interval_and_relation(unsigned level, polynomial_ref_vector const& level_ps, svector<char>& poly_has_roots) {
m_level = level;
m_rel.clear();
reset_interval(m_I[level]);
poly_has_roots.reset();
poly_has_roots.resize(level_ps.size(), false);
anum const& v = sample().value(level);
for (unsigned i = 0; i < level_ps.size(); ++i) {
poly* p = level_ps.get(i);
scoped_anum_vector roots(m_am);
// SMT-RAT caches isolated bound roots inside the constructed cell components
// (e.g., Section::isolatedRoot in OneCellCAD.h). Z3's levelwise currently
// recomputes root isolations on demand.
//
// Optimization: when the sample value is rational, use a closest-root isolation
// routine to avoid isolating all roots.
if (v.is_basic()) {
scoped_mpq s(m_am.qm());
m_am.to_rational(v, s);
svector<unsigned> idx;
m_am.isolate_roots_closest(polynomial_ref(p, m_pm), undef_var_assignment(sample(), level), s, roots, idx);
poly_has_roots[i] = !roots.empty();
SASSERT(roots.size() == idx.size());
for (unsigned k = 0; k < roots.size(); ++k)
m_rel.m_rfunc.emplace_back(m_am, p, idx[k], roots[k]);
continue;
}
m_am.isolate_roots(polynomial_ref(p, m_pm), undef_var_assignment(sample(), level), roots);
for (unsigned k = 0; k < roots.size(); ++k)
m_rel.m_rfunc.emplace_back(m_am, p, k + 1, roots[k]);
poly_has_roots[i] = !roots.empty();
// SMT-RAT style reduction: keep only the closest root below (<= v) and above (> v),
// or a single root if v is a root of p.
unsigned lower = UINT_MAX, upper = UINT_MAX;
for (unsigned k = 0; k < roots.size(); ++k) {
auto cmp = m_am.compare(roots[k], v);
if (cmp <= 0)
lower = k;
else {
upper = k;
break;
}
}
if (lower != UINT_MAX) {
m_rel.m_rfunc.emplace_back(m_am, p, lower + 1, roots[lower]);
if (m_am.eq(roots[lower], v))
continue; // only keep the section root for this polynomial
}
if (upper != UINT_MAX)
m_rel.m_rfunc.emplace_back(m_am, p, upper + 1, roots[upper]);
}
if (m_rel.m_rfunc.empty())
return;
anum const& v = sample().value(level);
// Comparator for sorting roots by value (with deterministic tie-breaking)
auto cmp = [&](root_function const& a, root_function const& b) {
if (a.ire.p == b.ire.p) // compare pointers
return a.ire.i < b.ire.i; // compare root indices in the same polynomial
auto r = m_am.compare(a.val, b.val);
if (r)
return r == sign_neg;
// Tie-break: same value - use polynomial id
unsigned ida = m_pm.id(a.ire.p);
unsigned idb = m_pm.id(b.ire.p);
return ida < idb;
};
// Partition: roots <= v come first, then roots > v
auto& rfs = m_rel.m_rfunc;
auto mid = std::partition(rfs.begin(), rfs.end(), [&](root_function const& f) {
return m_am.compare(f.val, v) <= 0;
});
// Sort each partition separately - O(n log n) comparisons between roots only
std::sort(rfs.begin(), mid, cmp);
std::sort(mid, rfs.end(), cmp);
// Sort each partition separately. For ties (equal root values), prefer
// low-degree defining polynomials as interval boundaries:
// - lower bound comes from the last element in the <= partition, so sort ties by degree descending
// - upper bound comes from the first element in the > partition, so sort ties by degree ascending
auto cmp_value = [&](root_function const& a, root_function const& b) {
auto r = m_am.compare(a.val, b.val);
if (r)
return r == sign_neg;
if (a.ire.p == b.ire.p)
return a.ire.i < b.ire.i;
return false;
};
auto tie_id = [&](root_function const& a, root_function const& b) {
unsigned ida = m_pm.id(a.ire.p);
unsigned idb = m_pm.id(b.ire.p);
if (ida != idb)
return ida < idb;
return a.ire.i < b.ire.i;
};
auto cmp_lo = [&](root_function const& a, root_function const& b) {
if (cmp_value(a, b))
return true;
if (cmp_value(b, a))
return false;
if (a.ire.p == b.ire.p)
return a.ire.i < b.ire.i;
unsigned dega = m_pm.degree(a.ire.p, level);
unsigned degb = m_pm.degree(b.ire.p, level);
if (dega != degb)
return dega > degb; // descending
return tie_id(a, b);
};
auto cmp_hi = [&](root_function const& a, root_function const& b) {
if (cmp_value(a, b))
return true;
if (cmp_value(b, a))
return false;
if (a.ire.p == b.ire.p)
return a.ire.i < b.ire.i;
unsigned dega = m_pm.degree(a.ire.p, level);
unsigned degb = m_pm.degree(b.ire.p, level);
if (dega != degb)
return dega < degb; // ascending
return tie_id(a, b);
};
std::sort(rfs.begin(), mid, cmp_lo);
std::sort(mid, rfs.end(), cmp_hi);
unsigned l_index = static_cast<unsigned>(-1);
unsigned u_index = static_cast<unsigned>(-1);
@ -448,17 +603,20 @@ struct levelwise::impl {
std::pair<unsigned, unsigned> key = id1 < id2 ? std::make_pair(id1, id2) : std::make_pair(id2, id1);
if (!added_pairs.insert(key).second)
continue;
insert_factorized(P, psc_resultant(a, b, level));
request_factorized(P, psc_resultant(a, b, level), inv_req::ord);
}
}
// Top level (m_n): add resultants between adjacent roots of different polynomials.
void add_adjacent_root_resultants(todo_set& P, polynomial_ref_vector const& top_ps) {
void add_adjacent_root_resultants(todo_set& P, polynomial_ref_vector const& top_ps, svector<char>& poly_has_roots) {
poly_has_roots.reset();
poly_has_roots.resize(top_ps.size(), false);
std::vector<std::pair<scoped_anum, poly*>> root_vals;
for (unsigned i = 0; i < top_ps.size(); ++i) {
poly* p = top_ps.get(i);
scoped_anum_vector roots(m_am);
m_am.isolate_roots(polynomial_ref(p, m_pm), undef_var_assignment(sample(), m_n), roots);
poly_has_roots[i] = !roots.empty();
for (unsigned k = 0; k < roots.size(); ++k) {
scoped_anum v(m_am);
m_am.set(v, roots[k]);
@ -483,16 +641,28 @@ struct levelwise::impl {
std::pair<unsigned, unsigned> key = id1 < id2 ? std::make_pair(id1, id2) : std::make_pair(id2, id1);
if (!added_pairs.insert(key).second)
continue;
insert_factorized(P, psc_resultant(p1, p2, m_n));
request_factorized(P, psc_resultant(p1, p2, m_n), inv_req::ord);
}
}
void process_level(todo_set& P, unsigned level, polynomial_ref_vector const& level_ps) {
void upgrade_bounds_to_ord(unsigned level, polynomial_ref_vector const& level_ps, std::vector<tagged_id>& tags) {
poly* lb = m_I[level].l;
poly* ub = m_I[level].u;
for (unsigned i = 0; i < level_ps.size(); ++i) {
poly* p = level_ps.get(i);
if (p == lb || p == ub)
tags[i].second = inv_req::ord;
}
}
void process_level(todo_set& P, unsigned level, polynomial_ref_vector const& level_ps, std::vector<tagged_id>& level_tags) {
SASSERT(level_ps.size() == level_tags.size());
// Line 10/11: detect nullification + pick a non-zero coefficient witness per p.
std::vector<polynomial_ref> witnesses;
witnesses.reserve(level_ps.size());
for (unsigned i = 0; i < level_ps.size(); ++i) {
polynomial_ref p(level_ps.get(i), m_pm);
SASSERT(level_tags[i].first == m_pm.id(p.get()));
polynomial_ref w = choose_nonzero_coeff(p, level);
if (!w)
fail();
@ -500,25 +670,48 @@ struct levelwise::impl {
}
// Lines 3-8: Θ + I_level + relation ≼
build_interval_and_relation(level, level_ps);
svector<char> poly_has_roots;
build_interval_and_relation(level, level_ps, poly_has_roots);
// SMT-RAT (levelwise/OneCellCAD.h) upgrades the polynomials defining the current
// bounds/sections to ORD_INV. Z3's levelwise does not currently implement SMT-RAT's
// dedicated section/sector heuristics (sectionHeuristic/sectorHeuristic) for choosing
// additional resultants/leading-coefficients; it instead uses relation_mode and the
// closest-root reduction when building Θ.
upgrade_bounds_to_ord(level, level_ps, level_tags);
// Lines 11-12 (Algorithm 1): add projections for each p
// Note: Algorithm 1 adds disc + ldcf for ALL polynomials (classical delineability)
// Algorithm 2 (projective delineability) would only add ldcf for bound-defining polys
// We additionally omit leading coefficients for rootless polynomials when possible
// (cf. projective delineability, Lemma 3.2).
for (unsigned i = 0; i < level_ps.size(); ++i) {
polynomial_ref p(level_ps.get(i), m_pm);
add_projections_for(P, p, level, witnesses[i]);
bool add_lc = true;
polynomial_ref lc(m_pm);
if (i < poly_has_roots.size() && !poly_has_roots[i]) {
unsigned deg = m_pm.degree(p, level);
lc = m_pm.coeff(p, level, deg);
if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0)
add_lc = false;
}
if (!add_lc && lc && witnesses[i].get() == lc.get() && !is_const(lc)) {
polynomial_ref alt = choose_nonzero_coeff(p, level, lc.get());
if (alt)
witnesses[i] = alt;
}
add_projections_for(P, p, level, witnesses[i], add_lc);
}
// Line 14 (Algorithm 1): add resultants for chosen relation pairs
add_relation_resultants(P, level);
}
void process_top_level(todo_set& P, polynomial_ref_vector const& top_ps) {
void process_top_level(todo_set& P, polynomial_ref_vector const& top_ps, std::vector<tagged_id>& top_tags) {
SASSERT(top_ps.size() == top_tags.size());
std::vector<polynomial_ref> witnesses;
witnesses.reserve(top_ps.size());
for (unsigned i = 0; i < top_ps.size(); ++i) {
polynomial_ref p(top_ps.get(i), m_pm);
SASSERT(top_tags[i].first == m_pm.id(p.get()));
polynomial_ref w = choose_nonzero_coeff(p, m_n);
if (!w)
fail();
@ -526,12 +719,30 @@ struct levelwise::impl {
}
// Resultants between adjacent root functions (a lightweight ordering for the top level).
add_adjacent_root_resultants(P, top_ps);
svector<char> poly_has_roots;
add_adjacent_root_resultants(P, top_ps, poly_has_roots);
// Projections (coeff witness, disc, leading coeff).
// Note: SMT-RAT's levelwise implementation additionally has dedicated heuristics for
// selecting resultants and selectively adding leading coefficients (see OneCellCAD.h,
// sectionHeuristic/sectorHeuristic). Z3's levelwise currently uses the relation_mode
// ordering heuristics instead of these specialized cases.
for (unsigned i = 0; i < top_ps.size(); ++i) {
polynomial_ref p(top_ps.get(i), m_pm);
add_projections_for(P, p, m_n, witnesses[i]);
bool add_lc = true;
polynomial_ref lc(m_pm);
if (i < poly_has_roots.size() && !poly_has_roots[i]) {
unsigned deg = m_pm.degree(p, m_n);
lc = m_pm.coeff(p, m_n, deg);
if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0)
add_lc = false;
}
if (!add_lc && lc && witnesses[i].get() == lc.get() && !is_const(lc)) {
polynomial_ref alt = choose_nonzero_coeff(p, m_n, lc.get());
if (alt)
witnesses[i] = alt;
}
add_projections_for(P, p, m_n, witnesses[i], add_lc);
}
}
@ -545,7 +756,7 @@ struct levelwise::impl {
for (unsigned i = 0; i < m_P.size(); ++i) {
polynomial_ref pi(m_P.get(i), m_pm);
for_each_distinct_factor(pi, [&](polynomial_ref const& f) {
P.insert(f.get());
request(P, f.get(), inv_req::sign);
});
}
@ -555,16 +766,18 @@ struct levelwise::impl {
// Process top level m_n (we project from x_{m_n} and do not return an interval for it).
if (P.max_var() == m_n) {
polynomial_ref_vector top_ps(m_pm);
P.extract_max_polys(top_ps);
process_top_level(P, top_ps);
std::vector<tagged_id> top_tags;
extract_max_tagged(P, top_ps, top_tags);
process_top_level(P, top_ps, top_tags);
}
// Now iteratively process remaining levels (descending).
while (!P.empty()) {
polynomial_ref_vector level_ps(m_pm);
var level = P.extract_max_polys(level_ps);
std::vector<tagged_id> level_tags;
var level = extract_max_tagged(P, level_ps, level_tags);
SASSERT(level < m_n);
process_level(P, level, level_ps);
process_level(P, level, level_ps, level_tags);
}
return m_I;