From 66ef8781624f15a80e47fc4a207ce3bb7e3722b9 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Sat, 10 Jan 2026 19:06:52 -1000 Subject: [PATCH] introduce isolate_root_closest Signed-off-by: Lev Nachmanson --- src/math/polynomial/algebraic_numbers.cpp | 238 +++++++++++++++++ src/math/polynomial/algebraic_numbers.h | 13 + src/nlsat/levelwise.cpp | 309 ++++++++++++++++++---- 3 files changed, 512 insertions(+), 48 deletions(-) diff --git a/src/math/polynomial/algebraic_numbers.cpp b/src/math/polynomial/algebraic_numbers.cpp index 4cdf9c4b8..ae78326de 100644 --- a/src/math/polynomial/algebraic_numbers.cpp +++ b/src/math/polynomial/algebraic_numbers.cpp @@ -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(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 & 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 & 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 & 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 & signs) { m_imp->isolate_roots(p, x2v, roots, signs); } diff --git a/src/math/polynomial/algebraic_numbers.h b/src/math/polynomial/algebraic_numbers.h index 88792bbc2..46cb3c6da 100644 --- a/src/math/polynomial/algebraic_numbers.h +++ b/src/math/polynomial/algebraic_numbers.h @@ -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 & indices); + /** \brief Isolate the roots of the given polynomial, and compute its sign between them. */ diff --git a/src/nlsat/levelwise.cpp b/src/nlsat/levelwise.cpp index cd494db01..bad826a5f 100644 --- a/src/nlsat/levelwise.cpp +++ b/src/nlsat/levelwise.cpp @@ -5,8 +5,10 @@ #include "math/polynomial/polynomial.h" #include "nlsat_common.h" #include "util/z3_exception.h" +#include "util/vector.h" #include +#include #include #include #include @@ -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(a) < static_cast(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 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(m_req.get(id, static_cast(inv_req::none))); + inv_req nxt = max_req(cur, req); + if (nxt != cur) + m_req.setx(id, static_cast(nxt), static_cast(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; // + + var extract_max_tagged(todo_set& P, polynomial_ref_vector& max_ps, std::vector& 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(m_req.get(id, static_cast(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(inv_req::none), static_cast(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(deg); j >= 0; --j) { - coeff = m_pm.coeff(p, x, static_cast(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& 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 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(-1); unsigned u_index = static_cast(-1); @@ -448,17 +603,20 @@ struct levelwise::impl { std::pair 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& poly_has_roots) { + poly_has_roots.reset(); + poly_has_roots.resize(top_ps.size(), false); std::vector> 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 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& 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& level_tags) { + SASSERT(level_ps.size() == level_tags.size()); // Line 10/11: detect nullification + pick a non-zero coefficient witness per p. std::vector 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 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& top_tags) { + SASSERT(top_ps.size() == top_tags.size()); std::vector 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 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 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 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;