mirror of
https://github.com/Z3Prover/z3
synced 2026-02-08 10:07:59 +00:00
1464 lines
64 KiB
C++
1464 lines
64 KiB
C++
#include "nlsat/levelwise.h"
|
|
#include "nlsat/nlsat_types.h"
|
|
#include "nlsat/nlsat_assignment.h"
|
|
#include "math/polynomial/algebraic_numbers.h"
|
|
#include "math/polynomial/polynomial.h"
|
|
#include "nlsat_common.h"
|
|
#include "util/z3_exception.h"
|
|
#include "util/vector.h"
|
|
#include "util/trace.h"
|
|
|
|
#include <algorithm>
|
|
#include <cstdint>
|
|
#include <set>
|
|
#include <unordered_map>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
static bool is_set(unsigned k) { return static_cast<int>(k) != -1; }
|
|
|
|
// Helper for sparse vector access with default value
|
|
template<typename T>
|
|
static T vec_get(std_vector<T> const& v, unsigned idx, T def) {
|
|
return idx < v.size() ? v[idx] : def;
|
|
}
|
|
|
|
// Helper for sparse vector write with auto-resize
|
|
template<typename T>
|
|
static void vec_setx(std_vector<T>& v, unsigned idx, T val, T def) {
|
|
if (idx >= v.size())
|
|
v.resize(idx + 1, def);
|
|
v[idx] = val;
|
|
}
|
|
|
|
namespace nlsat {
|
|
|
|
enum relation_mode {
|
|
biggest_cell,
|
|
spanning_tree
|
|
};
|
|
|
|
// 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 {
|
|
solver& m_solver;
|
|
polynomial_ref_vector const& m_P;
|
|
unsigned m_n; // maximal variable we project from
|
|
pmanager& m_pm;
|
|
anum_manager& m_am;
|
|
polynomial::cache& m_cache;
|
|
std_vector<root_function_interval> m_I; // intervals for variables 0..m_n-1
|
|
|
|
unsigned m_level = 0; // current level being processed
|
|
relation_mode m_rel_mode = biggest_cell;
|
|
unsigned m_spanning_tree_threshold = 3; // minimum both-side count for spanning tree
|
|
unsigned m_l_rf = UINT_MAX; // position of lower bound in m_rel.m_rfunc
|
|
unsigned m_u_rf = UINT_MAX; // position of upper bound in m_rel.m_rfunc (UINT_MAX in section case)
|
|
|
|
// Per-level state set by process_level/process_top_level
|
|
todo_set m_todo;
|
|
polynomial_ref_vector m_level_ps;
|
|
std_vector<polynomial_ref> m_witnesses;
|
|
std_vector<bool> m_poly_has_roots;
|
|
|
|
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.
|
|
std_vector<uint8_t> m_req;
|
|
|
|
// Vectors indexed by position in m_level_ps (more cache-friendly than maps)
|
|
mutable std_vector<uint8_t> m_side_mask; // bit0 = lower, bit1 = upper, 3 = both
|
|
mutable std_vector<bool> m_omit_lc;
|
|
mutable std_vector<bool> m_omit_disc;
|
|
mutable std_vector<unsigned> m_deg_in_order_graph; // degree of polynomial in resultant graph
|
|
mutable std_vector<unsigned> m_unique_neighbor; // UINT_MAX = not set, UINT_MAX-1 = multiple
|
|
|
|
assignment const& sample() const { return m_solver.sample(); }
|
|
|
|
// Utility: call fn for each distinct irreducible factor of poly
|
|
template <typename Func>
|
|
void for_each_distinct_factor(polynomial_ref const& poly_in, Func&& fn) {
|
|
if (!poly_in || is_zero(poly_in) || is_const(poly_in))
|
|
return;
|
|
polynomial_ref poly(poly_in);
|
|
polynomial_ref_vector factors(m_pm);
|
|
::nlsat::factor(poly, m_cache, factors);
|
|
for (unsigned i = 0; i < factors.size(); ++i) {
|
|
polynomial_ref f(factors.get(i), m_pm);
|
|
if (!f || is_zero(f) || is_const(f))
|
|
continue;
|
|
fn(f);
|
|
}
|
|
}
|
|
|
|
struct root_function {
|
|
scoped_anum val;
|
|
indexed_root_expr ire;
|
|
unsigned ps_idx; // index in m_level_ps
|
|
root_function(anum_manager& am, poly* p, unsigned idx, anum const& v, unsigned ps_idx)
|
|
: val(am), ire{ p, idx }, ps_idx(ps_idx) { am.set(val, v); }
|
|
root_function(root_function&& other) noexcept : val(other.val.m()), ire(other.ire), ps_idx(other.ps_idx) { val = other.val; }
|
|
root_function(root_function const&) = delete;
|
|
root_function& operator=(root_function const&) = delete;
|
|
root_function& operator=(root_function&& other) noexcept {
|
|
val = other.val;
|
|
ire = other.ire;
|
|
ps_idx = other.ps_idx;
|
|
return *this;
|
|
}
|
|
};
|
|
|
|
// Root functions (Theta) and the chosen relation (≼) on a given level.
|
|
struct relation_E {
|
|
std_vector<root_function> m_rfunc; // Θ: root functions on current level
|
|
std::set<std::pair<unsigned, unsigned>> m_pairs; // ≼ relation as unique m_level_ps index pairs
|
|
bool empty() const { return m_rfunc.empty() && m_pairs.empty(); }
|
|
void clear() {
|
|
m_pairs.clear();
|
|
m_rfunc.clear();
|
|
}
|
|
// Add pair by rfunc indices - canonicalizes and de-duplicates using ps_idx
|
|
void add_pair(unsigned j, unsigned k) {
|
|
unsigned a = m_rfunc[j].ps_idx;
|
|
unsigned b = m_rfunc[k].ps_idx;
|
|
if (a == b) return;
|
|
if (a > b) std::swap(a, b);
|
|
m_pairs.insert({a, b});
|
|
}
|
|
};
|
|
|
|
relation_E m_rel;
|
|
|
|
// Check that the relation forms a valid connected structure for order-invariance.
|
|
// Every root function on each side must be connected to the boundary through edges.
|
|
bool relation_invariant() const {
|
|
auto const& rfs = m_rel.m_rfunc;
|
|
unsigned n = rfs.size();
|
|
if (n == 0) return true;
|
|
|
|
// Build adjacency list from pairs (using ps_idx)
|
|
// Use vectors indexed by ps_idx for efficiency
|
|
unsigned max_ps_idx = 0;
|
|
for (auto const& rf : rfs)
|
|
if (rf.ps_idx > max_ps_idx) max_ps_idx = rf.ps_idx;
|
|
for (auto const& pr : m_rel.m_pairs) {
|
|
if (pr.first > max_ps_idx) max_ps_idx = pr.first;
|
|
if (pr.second > max_ps_idx) max_ps_idx = pr.second;
|
|
}
|
|
|
|
std_vector<std_vector<unsigned>> adj(max_ps_idx + 1);
|
|
for (auto const& pr : m_rel.m_pairs) {
|
|
adj[pr.first].push_back(pr.second);
|
|
adj[pr.second].push_back(pr.first);
|
|
}
|
|
|
|
// BFS to find all ps_idx reachable from a starting ps_idx
|
|
auto reachable_from = [&](unsigned start_ps_idx) -> std_vector<bool> {
|
|
std_vector<bool> visited(max_ps_idx + 1, false);
|
|
std_vector<unsigned> queue;
|
|
queue.push_back(start_ps_idx);
|
|
visited[start_ps_idx] = true;
|
|
while (!queue.empty()) {
|
|
unsigned curr = queue.back();
|
|
queue.pop_back();
|
|
for (unsigned neighbor : adj[curr]) {
|
|
if (!visited[neighbor]) {
|
|
visited[neighbor] = true;
|
|
queue.push_back(neighbor);
|
|
}
|
|
}
|
|
}
|
|
return visited;
|
|
};
|
|
|
|
// Check lower side: all rfuncs in [0, m_l_rf] must be connected to boundary
|
|
if (is_set(m_l_rf) && m_l_rf < n) {
|
|
unsigned boundary_ps_idx = rfs[m_l_rf].ps_idx;
|
|
std_vector<bool> reachable = reachable_from(boundary_ps_idx);
|
|
for (unsigned i = 0; i <= m_l_rf; ++i) {
|
|
unsigned ps_idx = rfs[i].ps_idx;
|
|
if (!reachable[ps_idx]) {
|
|
TRACE(lws, tout << "INVARIANT FAIL: lower side rfunc " << i
|
|
<< " (ps_idx=" << ps_idx << ") not connected to boundary "
|
|
<< m_l_rf << " (ps_idx=" << boundary_ps_idx << ")\n";);
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Check upper side: all rfuncs in [m_u_rf, n-1] must be connected to boundary
|
|
if (is_set(m_u_rf) && m_u_rf < n) {
|
|
unsigned boundary_ps_idx = rfs[m_u_rf].ps_idx;
|
|
std_vector<bool> reachable = reachable_from(boundary_ps_idx);
|
|
for (unsigned i = m_u_rf; i < n; ++i) {
|
|
unsigned ps_idx = rfs[i].ps_idx;
|
|
if (!reachable[ps_idx]) {
|
|
TRACE(lws, tout << "INVARIANT FAIL: upper side rfunc " << i
|
|
<< " (ps_idx=" << ps_idx << ") not connected to boundary "
|
|
<< m_u_rf << " (ps_idx=" << boundary_ps_idx << ")\n";);
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
// Weight function for estimating projection complexity.
|
|
// w(p, level) estimates the "cost" of polynomial p at the given level.
|
|
// w(disc(a)) ≈ 2*w(a), w(res(a,b)) ≈ w(a) + w(b)
|
|
unsigned w(poly* p, unsigned level) const {
|
|
if (!p) return 0;
|
|
return m_pm.degree(p, level);
|
|
}
|
|
|
|
// Estimate resultant weight for m_rel.m_pairs (pairs of m_level_ps indices)
|
|
unsigned estimate_resultant_weight() const {
|
|
unsigned total = 0;
|
|
for (auto const& pr : m_rel.m_pairs)
|
|
total += w(m_level_ps.get(pr.first), m_level) + w(m_level_ps.get(pr.second), m_level);
|
|
|
|
return total;
|
|
}
|
|
|
|
// Estimate leading coefficient weight for polynomials where omit_lc is false
|
|
unsigned estimate_lc_weight() const {
|
|
unsigned total = 0;
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
if (vec_get(m_omit_lc, i, false))
|
|
continue;
|
|
poly* p = m_level_ps.get(i);
|
|
unsigned deg = m_pm.degree(p, m_level);
|
|
if (deg > 0)
|
|
total += w(m_pm.coeff(p, m_level, deg), m_level - 1);
|
|
}
|
|
return total;
|
|
}
|
|
|
|
// Choose the best sector heuristic based on estimated weight.
|
|
// Also fills m_rel.m_pairs with the winning heuristic's pairs.
|
|
// Note: spanning_tree is handled at a higher level (process_level_sector)
|
|
impl(
|
|
solver& solver,
|
|
polynomial_ref_vector const& ps,
|
|
var max_x,
|
|
assignment const&,
|
|
pmanager& pm,
|
|
anum_manager& am,
|
|
polynomial::cache& cache)
|
|
: m_solver(solver),
|
|
m_P(ps),
|
|
m_n(max_x),
|
|
m_pm(pm),
|
|
m_am(am),
|
|
m_cache(cache),
|
|
m_todo(m_cache, true),
|
|
m_level_ps(m_pm),
|
|
m_psc_tmp(m_pm) {
|
|
m_I.reserve(m_n);
|
|
for (unsigned i = 0; i < m_n; ++i)
|
|
m_I.emplace_back(m_pm);
|
|
|
|
m_spanning_tree_threshold = m_solver.lws_spt_threshold();
|
|
}
|
|
|
|
void fail() { throw nullified_poly_exception(); }
|
|
|
|
static void reset_interval(root_function_interval& I) {
|
|
I.section = false;
|
|
I.l = nullptr;
|
|
I.u = nullptr;
|
|
I.l_index = 0;
|
|
I.u_index = 0;
|
|
}
|
|
|
|
// PSC-based discriminant: returns the first non-zero, non-constant element from the PSC chain.
|
|
polynomial_ref psc_discriminant(polynomial_ref const& p_in, unsigned x) {
|
|
if (!p_in || is_zero(p_in) || is_const(p_in))
|
|
return polynomial_ref(m_pm);
|
|
if (m_pm.degree(p_in, x) < 2)
|
|
return polynomial_ref(m_pm);
|
|
polynomial_ref p(p_in);
|
|
polynomial_ref d = derivative(p, x);
|
|
polynomial_ref_vector& chain = m_psc_tmp;
|
|
chain.reset();
|
|
m_cache.psc_chain(p, d, x, chain);
|
|
polynomial_ref disc(m_pm);
|
|
for (unsigned i = 0; i < chain.size(); ++i) {
|
|
disc = polynomial_ref(chain.get(i), m_pm);
|
|
if (!disc || is_zero(disc))
|
|
continue;
|
|
if (is_const(disc))
|
|
return polynomial_ref(m_pm); // constant means discriminant is non-zero constant
|
|
return disc;
|
|
}
|
|
return polynomial_ref(m_pm);
|
|
}
|
|
|
|
// PSC-based resultant: returns the first non-zero, non-constant element from the PSC chain.
|
|
polynomial_ref psc_resultant(poly* a, poly* b, unsigned x) {
|
|
if (!a || !b)
|
|
return polynomial_ref(m_pm);
|
|
polynomial_ref pa(a, m_pm);
|
|
polynomial_ref pb(b, m_pm);
|
|
polynomial_ref_vector& chain = m_psc_tmp;
|
|
chain.reset();
|
|
m_cache.psc_chain(pa, pb, x, chain);
|
|
polynomial_ref r(m_pm);
|
|
// Iterate forward: S[0] is the resultant (after reverse in psc_chain)
|
|
for (unsigned i = 0; i < chain.size(); ++i) {
|
|
r = polynomial_ref(chain.get(i), m_pm);
|
|
if (!r || is_zero(r))
|
|
continue;
|
|
if (is_const(r))
|
|
return polynomial_ref(m_pm); // constant means polynomials are coprime
|
|
return r;
|
|
}
|
|
return polynomial_ref(m_pm);
|
|
}
|
|
|
|
poly* request(poly* p, inv_req req) {
|
|
if (!p || req == inv_req::none)
|
|
return p;
|
|
p = m_todo.insert(p);
|
|
unsigned id = m_pm.id(p);
|
|
auto cur = static_cast<inv_req>(vec_get(m_req, id, static_cast<uint8_t>(inv_req::none)));
|
|
inv_req nxt = max_req(cur, req);
|
|
if (nxt != cur)
|
|
vec_setx(m_req, id, static_cast<uint8_t>(nxt), static_cast<uint8_t>(inv_req::none));
|
|
return p;
|
|
}
|
|
|
|
void request_factorized(polynomial_ref const& poly, inv_req req) {
|
|
for_each_distinct_factor(poly, [&](polynomial_ref const& f) {
|
|
TRACE(lws,
|
|
m_pm.display(tout << " request_factorized: factor=", f) << " at level " << m_pm.max_var(f) << "\n";
|
|
);
|
|
request(f.get(), req); // inherit tag across factorization (SMT-RAT appendOnCorrectLevel)
|
|
});
|
|
}
|
|
|
|
// Extract polynomials at max level from m_todo into m_level_ps.
|
|
// Sets m_level to the extracted level. Requirements remain in m_req.
|
|
void extract_max_tagged() {
|
|
m_level = m_todo.extract_max_polys(m_level_ps);
|
|
// Ensure all extracted polynomials have at least SIGN requirement
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
unsigned id = m_pm.id(m_level_ps.get(i));
|
|
if (vec_get(m_req, id, static_cast<uint8_t>(inv_req::none)) == static_cast<uint8_t>(inv_req::none))
|
|
vec_setx(m_req, id, static_cast<uint8_t>(inv_req::sign), static_cast<uint8_t>(inv_req::none));
|
|
}
|
|
}
|
|
|
|
// Get the requirement for polynomial at index i in m_level_ps
|
|
inv_req get_req(unsigned i) const {
|
|
unsigned id = m_pm.id(m_level_ps.get(i));
|
|
return static_cast<inv_req>(vec_get(m_req, id, static_cast<uint8_t>(inv_req::sign)));
|
|
}
|
|
|
|
// Set the requirement for polynomial at index i in m_level_ps
|
|
void set_req(unsigned i, inv_req req) {
|
|
unsigned id = m_pm.id(m_level_ps.get(i));
|
|
vec_setx(m_req, id, static_cast<uint8_t>(req), static_cast<uint8_t>(inv_req::none));
|
|
}
|
|
|
|
// Select a coefficient c of p (wrt x) such that c(s) != 0, or return null.
|
|
// 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) {
|
|
unsigned deg = m_pm.degree(p, x);
|
|
polynomial_ref coeff(m_pm);
|
|
|
|
// 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))
|
|
continue;
|
|
if (is_const(coeff))
|
|
return coeff;
|
|
}
|
|
|
|
// 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))
|
|
continue;
|
|
if (m_am.eval_sign_at(coeff, sample()) == 0)
|
|
continue;
|
|
|
|
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_projection_for_poly(polynomial_ref const& p, unsigned x, polynomial_ref const& nonzero_coeff, bool add_leading_coeff, bool add_discriminant) {
|
|
TRACE(lws,
|
|
m_pm.display(tout << " add_projection_for_poly: p=", p)
|
|
<< " x=" << x << " add_lc=" << add_leading_coeff << " add_disc=" << add_discriminant << "\n";
|
|
);
|
|
// Line 11 (non-null witness coefficient)
|
|
if (nonzero_coeff && !is_const(nonzero_coeff))
|
|
request_factorized(nonzero_coeff, inv_req::sign);
|
|
|
|
// Line 12 (disc + leading coefficient)
|
|
if (add_discriminant)
|
|
request_factorized(psc_discriminant(p, x), inv_req::ord);
|
|
if (add_leading_coeff) {
|
|
unsigned deg = m_pm.degree(p, x);
|
|
polynomial_ref lc(m_pm);
|
|
lc = m_pm.coeff(p, x, deg);
|
|
TRACE(lws, m_pm.display(tout << " adding lc: ", lc) << "\n";);
|
|
request_factorized(lc, inv_req::sign);
|
|
}
|
|
}
|
|
|
|
// ============================================================================
|
|
// noLdcf helpers - compute which leading coefficients can be omitted
|
|
// Inspired by SMT-RAT's "noLdcf" optimization in OneCellCAD.h.
|
|
// ============================================================================
|
|
|
|
// Compute side_mask: track which side(s) each polynomial appears on
|
|
// bit0 = lower (<= sample), bit1 = upper (> sample), 3 = both sides
|
|
void compute_side_mask() {
|
|
if (!is_set(m_l_rf))
|
|
return;
|
|
m_side_mask.resize(m_level_ps.size(), 0);
|
|
for (unsigned i = 0; i < m_rel.m_rfunc.size(); i ++) {
|
|
unsigned ps_idx = m_rel.m_rfunc[i].ps_idx;
|
|
if (i <= m_l_rf)
|
|
m_side_mask[ps_idx] |= 1;
|
|
else
|
|
m_side_mask[ps_idx] |= 2;
|
|
}
|
|
}
|
|
|
|
// Check if lower and upper bounds come from the same polynomial
|
|
bool same_boundary_poly() const {
|
|
if (!is_set(m_l_rf) || !is_set(m_u_rf))
|
|
return false;
|
|
return m_rel.m_rfunc[m_l_rf].ps_idx == m_rel.m_rfunc[m_u_rf].ps_idx;
|
|
}
|
|
|
|
// Get lower bound polynomial index (UINT_MAX if not set)
|
|
unsigned lower_bound_idx() const {
|
|
return is_set(m_l_rf) ? m_rel.m_rfunc[m_l_rf].ps_idx : UINT_MAX;
|
|
}
|
|
|
|
// Get upper bound polynomial index (UINT_MAX if not set)
|
|
unsigned upper_bound_idx() const {
|
|
return is_set(m_u_rf) ? m_rel.m_rfunc[m_u_rf].ps_idx : UINT_MAX;
|
|
}
|
|
|
|
// Compute deg: count distinct resultant-neighbors for each polynomial
|
|
// m_pairs contains index pairs into m_level_ps
|
|
void compute_resultant_graph_degree() {
|
|
m_deg_in_order_graph.clear();
|
|
m_deg_in_order_graph.resize(m_level_ps.size(), 0);
|
|
m_unique_neighbor.clear();
|
|
m_unique_neighbor.resize(m_level_ps.size(), UINT_MAX);
|
|
for (auto const& pr : m_rel.m_pairs) {
|
|
unsigned idx1 = pr.first;
|
|
unsigned idx2 = pr.second;
|
|
m_deg_in_order_graph[idx1]++;
|
|
m_deg_in_order_graph[idx2]++;
|
|
// Track unique neighbor
|
|
auto update_neighbor = [&](unsigned idx, unsigned other) {
|
|
if (m_unique_neighbor[idx] == UINT_MAX)
|
|
m_unique_neighbor[idx] = other;
|
|
else if (m_unique_neighbor[idx] != other)
|
|
m_unique_neighbor[idx] = UINT_MAX - 1; // multiple neighbors
|
|
};
|
|
update_neighbor(idx1, idx2);
|
|
update_neighbor(idx2, idx1);
|
|
}
|
|
}
|
|
|
|
// Compute noDisc for spanning_tree mode.
|
|
// A polynomial's discriminant can be omitted if it's a "leaf" in the resultant graph,
|
|
// connected only to a boundary polynomial.
|
|
//
|
|
// Relation diagram for noDisc-eligible polynomial p:
|
|
//
|
|
// lower side upper side
|
|
// ─────────────────────────────────────────────
|
|
// p ──Res── lb ═══════ ub
|
|
// ↑ ↑
|
|
// (boundary polynomials)
|
|
//
|
|
// Here p is a leaf with deg=1, connected only to lb (lower boundary).
|
|
// Its discriminant can be skipped because:
|
|
// - Root ordering of p relative to lb is fixed by Res(p, lb)
|
|
// - No other polynomial depends on p's internal root ordering
|
|
// - The boundary polynomial's discriminant handles its own root ordering
|
|
//
|
|
// Symmetrically, a leaf connected only to ub qualifies too.
|
|
void compute_omit_disc_for_spanning_tree() {
|
|
m_omit_disc.clear();
|
|
m_omit_disc.resize(m_level_ps.size(), false);
|
|
if (m_rel.m_rfunc.empty() || m_rel.m_pairs.empty())
|
|
return;
|
|
|
|
// Need graph info
|
|
compute_resultant_graph_degree();
|
|
|
|
unsigned l_bound_idx = lower_bound_idx();
|
|
unsigned u_bound_idx = upper_bound_idx();
|
|
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
// Skip boundary polynomials
|
|
if (i == l_bound_idx || i == u_bound_idx)
|
|
continue;
|
|
// Only skip if poly is not ORD-required by upstream projection
|
|
if (get_req(i) == inv_req::ord)
|
|
continue;
|
|
// Must be a leaf
|
|
if (m_deg_in_order_graph[i] != 1)
|
|
continue;
|
|
// Must be connected only to a boundary polynomial
|
|
unsigned neighbor = m_unique_neighbor[i];
|
|
if (neighbor != l_bound_idx && neighbor != u_bound_idx)
|
|
continue;
|
|
// If poly vanishes at sample, we need its discriminant for coinciding roots
|
|
poly* p = m_level_ps.get(i);
|
|
if (m_am.eval_sign_at(polynomial_ref(p, m_pm), sample()) == 0)
|
|
continue;
|
|
m_omit_disc[i] = true;
|
|
}
|
|
}
|
|
|
|
// Compute noDisc for same_boundary_poly case (sector with same lower/upper poly).
|
|
// When lower and upper bounds come from the same polynomial t, non-bound polynomials
|
|
// can omit their discriminant IF they don't vanish at sample AND their discriminant
|
|
// is non-zero at sample.
|
|
//
|
|
// Theory: If poly p doesn't vanish at sample, all its roots are on one side of the sample.
|
|
// Sign(p at sample) = sign(LC) * (-1)^(#roots below). If roots coincide (disc=0),
|
|
// two roots merge/disappear, parity changes by 2, sign unchanged. So disc not needed.
|
|
//
|
|
// However, if disc(p) = 0 at sample, p has a multiple root that can SPLIT when parameters
|
|
// change. When a double root splits, it creates two distinct roots, potentially one on
|
|
// each side of the sample, causing a sign change. So we must keep disc when disc=0.
|
|
//
|
|
// Example: p = (x-1)^2 at sample has disc=0. If it splits to (x-1)(x-9), the sample
|
|
// at x=5 has different sign than before.
|
|
void compute_omit_disc_for_same_boundary() {
|
|
m_omit_disc.clear();
|
|
m_omit_disc.resize(m_level_ps.size(), false);
|
|
if (!same_boundary_poly())
|
|
return;
|
|
|
|
unsigned bound_idx = lower_bound_idx();
|
|
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
// Skip boundary polynomial
|
|
if (i == bound_idx)
|
|
continue;
|
|
// Only skip if poly is not ORD-required
|
|
if (get_req(i) == inv_req::ord)
|
|
continue;
|
|
// Only skip if poly has roots (rootless needs disc to stay rootless)
|
|
if (!poly_has_roots(i))
|
|
continue;
|
|
poly* p = m_level_ps.get(i);
|
|
polynomial_ref p_ref(p, m_pm);
|
|
// Keep disc if poly vanishes at sample (root coincides with boundary)
|
|
if (m_am.eval_sign_at(p_ref, sample()) == 0)
|
|
continue;
|
|
// Keep disc if discriminant is zero at sample (multiple root can split)
|
|
polynomial_ref disc = psc_discriminant(p_ref, m_level);
|
|
if (disc && !is_const(disc) && m_am.eval_sign_at(disc, sample()) == 0)
|
|
continue;
|
|
m_omit_disc[i] = true;
|
|
}
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------
|
|
// noLdcf heuristic helpers
|
|
// ----------------------------------------------------------------------------
|
|
|
|
// Compute noLdcf: which leading coefficients can be omitted.
|
|
//
|
|
// Theory: A polynomial p with roots on BOTH sides of the sample has resultants
|
|
// with both bounds (Res(p, lb) and Res(p, ub)), ensuring its roots stay ordered.
|
|
// Even if p's LC becomes zero (degree drops), the existing roots remain correctly
|
|
// ordered via these resultant constraints.
|
|
//
|
|
// - biggest_cell mode (require_leaf=false): all non-bound polys with both-side roots
|
|
// - spanning_tree mode (require_leaf=true): only leaves (deg == 1) with both-side roots
|
|
//
|
|
void compute_omit_lc_both_sides(bool require_leaf) {
|
|
m_omit_lc.clear();
|
|
m_omit_lc.resize(m_level_ps.size(), false);
|
|
if (m_rel.m_rfunc.empty() || m_rel.m_pairs.empty() || m_side_mask.empty())
|
|
return;
|
|
|
|
if (require_leaf)
|
|
compute_resultant_graph_degree();
|
|
|
|
unsigned l_bound_idx = lower_bound_idx();
|
|
unsigned u_bound_idx = upper_bound_idx();
|
|
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
if (m_side_mask[i] != 3)
|
|
continue;
|
|
if (require_leaf && m_deg_in_order_graph[i] != 1)
|
|
continue;
|
|
if (i == l_bound_idx || i == u_bound_idx)
|
|
continue;
|
|
m_omit_lc[i] = true;
|
|
}
|
|
}
|
|
|
|
// Relation construction heuristics (same intent as previous implementation).
|
|
void fill_relation_sector_biggest_cell() {
|
|
TRACE(lws,
|
|
tout << " fill_biggest_cell: m_l_rf=" << m_l_rf << ", m_u_rf=" << m_u_rf << ", rfunc.size=" << m_rel.m_rfunc.size() << "\n";
|
|
);
|
|
if (is_set(m_l_rf))
|
|
for (unsigned j = 0; j < m_l_rf; ++j) {
|
|
TRACE(lws, tout << " add_pair(" << j << ", " << m_l_rf << ")\n";);
|
|
m_rel.add_pair(j, m_l_rf);
|
|
}
|
|
|
|
if (is_set(m_u_rf))
|
|
for (unsigned j = m_u_rf + 1; j < m_rel.m_rfunc.size(); ++j) {
|
|
TRACE(lws, tout << " add_pair(" << m_u_rf << ", " << j << ")\n";);
|
|
m_rel.add_pair(m_u_rf, j);
|
|
}
|
|
|
|
if (is_set(m_l_rf) && is_set(m_u_rf)) {
|
|
SASSERT(m_l_rf + 1 == m_u_rf);
|
|
TRACE(lws, tout << " add_pair(" << m_l_rf << ", " << m_u_rf << ")\n";);
|
|
m_rel.add_pair(m_l_rf, m_u_rf);
|
|
}
|
|
TRACE(lws,
|
|
tout << " fill_biggest_cell done: pairs.size=" << m_rel.m_pairs.size() << "\n";
|
|
);
|
|
}
|
|
|
|
void fill_relation_pairs_for_section_biggest_cell() {
|
|
auto const& rfs = m_rel.m_rfunc;
|
|
unsigned n = rfs.size();
|
|
if (n == 0)
|
|
return;
|
|
SASSERT(is_set(m_l_rf));
|
|
SASSERT(m_l_rf < n);
|
|
for (unsigned j = 0; j < m_l_rf; ++j)
|
|
m_rel.add_pair(j, m_l_rf);
|
|
for (unsigned j = m_l_rf + 1; j < n; ++j)
|
|
m_rel.add_pair(m_l_rf, j);
|
|
}
|
|
|
|
// ============================================================================
|
|
// Spanning tree heuristic based on the Reaching Orders Lemma.
|
|
// For polynomials appearing on both sides of the sample, we build a spanning
|
|
// tree that ensures order-invariance on both sides with n-1 edges.
|
|
// ============================================================================
|
|
|
|
// Build spanning tree on both-side polynomials using the lemma construction.
|
|
// Adds pairs directly to m_rel.m_pairs. Returns true if tree was built.
|
|
// K = lower rfunc positions, f = upper rfunc positions
|
|
void build_both_side_spanning_tree() {
|
|
auto const& rfs = m_rel.m_rfunc;
|
|
unsigned n = rfs.size();
|
|
SASSERT(n > 0 && is_set(m_l_rf) && is_set(m_u_rf));
|
|
SASSERT(!is_section());
|
|
SASSERT(!same_boundary_poly());
|
|
|
|
// Collect both-side polynomials with their rfunc indices on each side
|
|
struct both_info {
|
|
unsigned ps_idx;
|
|
unsigned lower_rf; // rfunc index on lower side
|
|
unsigned upper_rf; // rfunc index on upper side
|
|
};
|
|
std_vector<both_info> both;
|
|
|
|
// For sector: lower side is [0, m_l_rf], upper side is [m_u_rf, n)
|
|
|
|
// Build map from ps_idx to rfunc indices
|
|
std_vector<unsigned> lower_rf(m_level_ps.size(), UINT_MAX);
|
|
std_vector<unsigned> upper_rf(m_level_ps.size(), UINT_MAX);
|
|
|
|
for (unsigned i = 0; i <= m_l_rf; ++i)
|
|
lower_rf[rfs[i].ps_idx] = i;
|
|
for (unsigned i = m_u_rf; i < n; ++i)
|
|
upper_rf[rfs[i].ps_idx] = i;
|
|
// Collect both-side polynomials
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
if (m_side_mask[i] == 3) {
|
|
SASSERT(lower_rf[i] != UINT_MAX && upper_rf[i] != UINT_MAX);
|
|
both.push_back({i, lower_rf[i], upper_rf[i]});
|
|
}
|
|
}
|
|
|
|
SASSERT(both.size() >= m_spanning_tree_threshold);
|
|
|
|
// Sort by lower_rf (root position on lower side)
|
|
std::sort(both.begin(), both.end(),
|
|
[](both_info const& a, both_info const& b) { return a.lower_rf < b.lower_rf; });
|
|
|
|
// Build spanning tree using Reaching Orders Lemma (from the paper):
|
|
// Process elements in order of increasing lower_rf.
|
|
// For each new element a = min(remaining), connect to c where f(c) = min(f(K \ {a}))
|
|
// i.e., c has minimum upper_rf among all elements except a.
|
|
|
|
// Root of in-arborescence on lower side is max(K) = last element (max lower_rf)
|
|
// Root of out-arborescence on upper side is min(f(K)) = element with min upper_rf
|
|
unsigned upper_root_idx = 0;
|
|
for (unsigned i = 1; i < both.size(); ++i)
|
|
if (both[i].upper_rf < both[upper_root_idx].upper_rf)
|
|
upper_root_idx = i;
|
|
unsigned lower_root_idx = both.size() - 1;
|
|
|
|
// Process in order of lower_rf
|
|
// First element (index 0) has min lower_rf
|
|
for (unsigned a_idx = 0; a_idx < both.size() - 1; ++a_idx) {
|
|
// Find c = element with min upper_rf among all elements except a_idx
|
|
// Since we process in lower_rf order, elements [a_idx+1, n) are "K' = K \ {a}"
|
|
// and we need min upper_rf among them
|
|
unsigned c_idx = UINT_MAX;
|
|
for (unsigned i = a_idx + 1; i < both.size(); ++i)
|
|
if (c_idx == UINT_MAX || both[i].upper_rf < both[c_idx].upper_rf)
|
|
c_idx = i;
|
|
SASSERT(c_idx != UINT_MAX);
|
|
|
|
// Add edge {a, c}
|
|
unsigned ps_a = both[a_idx].ps_idx;
|
|
unsigned ps_c = both[c_idx].ps_idx;
|
|
m_rel.m_pairs.insert({std::min(ps_a, ps_c), std::max(ps_a, ps_c)});
|
|
}
|
|
|
|
|
|
// Check arborescence invariants
|
|
auto arb_invariant = [&]() {
|
|
// Reconstruct parent[] from the algorithm logic
|
|
std_vector<unsigned> parent(both.size(), UINT_MAX);
|
|
for (unsigned a_idx = 0; a_idx < both.size() - 1; ++a_idx) {
|
|
unsigned c_idx = UINT_MAX;
|
|
for (unsigned i = a_idx + 1; i < both.size(); ++i)
|
|
if (c_idx == UINT_MAX || both[i].upper_rf < both[c_idx].upper_rf)
|
|
c_idx = i;
|
|
parent[a_idx] = c_idx;
|
|
}
|
|
|
|
// Check it's a tree: exactly n-1 edges
|
|
unsigned edge_count = 0;
|
|
for (unsigned i = 0; i < both.size(); ++i)
|
|
if (parent[i] != UINT_MAX) edge_count++;
|
|
SASSERT(edge_count == both.size() - 1); // tree has n-1 edges
|
|
|
|
// Check roots are at extremes
|
|
// lower_root_idx should have max lower_rf
|
|
for (unsigned i = 0; i < both.size(); ++i)
|
|
SASSERT(both[i].lower_rf <= both[lower_root_idx].lower_rf);
|
|
// upper_root_idx should have min upper_rf
|
|
for (unsigned i = 0; i < both.size(); ++i)
|
|
SASSERT(both[i].upper_rf >= both[upper_root_idx].upper_rf);
|
|
|
|
// Root of in-arborescence (max lower_rf) has no parent
|
|
SASSERT(parent[lower_root_idx] == UINT_MAX);
|
|
for (unsigned i = 0; i < both.size(); ++i)
|
|
if (i != lower_root_idx)
|
|
SASSERT(parent[i] != UINT_MAX); // non-root has exactly 1 parent
|
|
|
|
// Check in-arborescence on left (lower) side:
|
|
// Each node can reach lower_root by following parent pointers
|
|
// and lower_rf increases (going right) at each step
|
|
for (unsigned i = 0; i < both.size(); ++i) {
|
|
unsigned curr = i;
|
|
unsigned steps = 0;
|
|
while (curr != lower_root_idx) {
|
|
unsigned p = parent[curr];
|
|
SASSERT(p != UINT_MAX);
|
|
// Going to parent increases lower_rf (toward root which has max lower_rf)
|
|
SASSERT(both[curr].lower_rf < both[p].lower_rf);
|
|
curr = p;
|
|
steps++;
|
|
SASSERT(steps < both.size()); // prevent infinite loop
|
|
}
|
|
}
|
|
|
|
// Check out-arborescence on right (upper) side:
|
|
// Re-orient edges based on upper_rf: edge goes from smaller to larger upper_rf
|
|
// Build adjacency from the edges in parent[]
|
|
std_vector<std_vector<unsigned>> out_children(both.size());
|
|
for (unsigned i = 0; i < both.size(); ++i) {
|
|
if (parent[i] != UINT_MAX) {
|
|
unsigned p = parent[i];
|
|
// Edge {i, p} exists. Orient based on upper_rf.
|
|
if (both[i].upper_rf < both[p].upper_rf) {
|
|
// i has smaller upper_rf, so edge is i → p in out-arborescence
|
|
out_children[i].push_back(p);
|
|
} else {
|
|
// p has smaller upper_rf, so edge is p → i in out-arborescence
|
|
out_children[p].push_back(i);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Each node can be reached from upper_root by following out_children
|
|
// and upper_rf increases (going away from root which has min upper_rf)
|
|
std_vector<bool> reached(both.size(), false);
|
|
std_vector<unsigned> stack;
|
|
stack.push_back(upper_root_idx);
|
|
reached[upper_root_idx] = true;
|
|
while (!stack.empty()) {
|
|
unsigned curr = stack.back();
|
|
stack.pop_back();
|
|
for (unsigned child : out_children[curr]) {
|
|
// Going to child increases upper_rf (away from root which has min upper_rf)
|
|
SASSERT(both[child].upper_rf > both[curr].upper_rf);
|
|
SASSERT(!reached[child]); // tree, no cycles
|
|
reached[child] = true;
|
|
stack.push_back(child);
|
|
}
|
|
}
|
|
// All nodes reachable from upper_root
|
|
for (unsigned i = 0; i < both.size(); ++i)
|
|
SASSERT(reached[i]);
|
|
|
|
return true;
|
|
};
|
|
SASSERT(arb_invariant());
|
|
}
|
|
|
|
// Sector spanning tree heuristic:
|
|
// 1. Build spanning tree on both-side polynomials
|
|
// 2. Extend with lowest_degree for single-side polynomials
|
|
// Helper: Connect non-tree (single-side) polynomials to their respective boundaries
|
|
void connect_non_tree_to_bounds() {
|
|
auto const& rfs = m_rel.m_rfunc;
|
|
unsigned n = rfs.size();
|
|
|
|
// Lower side: connect single-side polys to lower boundary
|
|
for (unsigned j = 0; j < m_l_rf; ++j)
|
|
if (m_side_mask[rfs[j].ps_idx] != 3)
|
|
m_rel.add_pair(j, m_l_rf);
|
|
|
|
// Upper side: connect single-side polys to upper boundary
|
|
for (unsigned j = m_u_rf + 1; j < n; ++j)
|
|
if (m_side_mask[rfs[j].ps_idx] != 3)
|
|
m_rel.add_pair(m_u_rf, j);
|
|
}
|
|
|
|
// Helper: Connect spanning tree extremes to boundaries (when boundaries are different polys)
|
|
void connect_tree_extremes_to_bounds() {
|
|
auto const& rfs = m_rel.m_rfunc;
|
|
unsigned n = rfs.size();
|
|
|
|
// Find max lower both-side poly (closest to lower boundary from below)
|
|
unsigned max_lower_both = UINT_MAX;
|
|
for (unsigned j = 0; j < m_l_rf; ++j)
|
|
if (m_side_mask[rfs[j].ps_idx] == 3)
|
|
max_lower_both = j;
|
|
|
|
// Find min upper both-side poly (closest to upper boundary from above)
|
|
unsigned min_upper_both = UINT_MAX;
|
|
for (unsigned j = m_u_rf + 1; j < n; ++j) {
|
|
if (m_side_mask[rfs[j].ps_idx] == 3) {
|
|
min_upper_both = j;
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Connect tree extremes to boundaries
|
|
if (max_lower_both != UINT_MAX)
|
|
m_rel.add_pair(max_lower_both, m_l_rf);
|
|
if (min_upper_both != UINT_MAX)
|
|
m_rel.add_pair(m_u_rf, min_upper_both);
|
|
}
|
|
|
|
void fill_relation_sector_spanning_tree() {
|
|
build_both_side_spanning_tree();
|
|
connect_non_tree_to_bounds();
|
|
connect_tree_extremes_to_bounds(); // otherwise the trees on both sides are connected to bounds already
|
|
|
|
// Connect lower and upper boundaries
|
|
SASSERT(m_l_rf + 1 == m_u_rf);
|
|
m_rel.add_pair(m_l_rf, m_u_rf);
|
|
}
|
|
|
|
// Extract roots of polynomial p around sample value v at the given level,
|
|
// partitioning them into lhalf (roots <= v) and uhalf (roots > v).
|
|
// ps_idx is the index of p in m_level_ps.
|
|
// Returns whether the polynomial has any roots.
|
|
bool isolate_roots_around_sample(poly* p, unsigned ps_idx,
|
|
anum const& v, std_vector<root_function>& lhalf,
|
|
std_vector<root_function>& uhalf) {
|
|
scoped_anum_vector roots(m_am);
|
|
|
|
|
|
// When the sample value is rational use a closest-root isolation
|
|
// returning at most closest to the sample roots
|
|
if (v.is_basic()) {
|
|
scoped_mpq s(m_am.qm());
|
|
m_am.to_rational(v, s);
|
|
svector<unsigned> idx_vec;
|
|
m_am.isolate_roots_closest(polynomial_ref(p, m_pm), undef_var_assignment(sample(), m_level), s, roots, idx_vec);
|
|
SASSERT(roots.size() == idx_vec.size());
|
|
SASSERT(roots.size() < 2 ||(roots.size() == 2 && m_am.compare(roots[0], v) < 0 && m_am.compare(roots[1], v) > 0));
|
|
for (unsigned k = 0; k < roots.size(); ++k) {
|
|
if (m_am.compare(roots[k], v) <= 0)
|
|
lhalf.emplace_back(m_am, p, idx_vec[k], roots[k], ps_idx);
|
|
else
|
|
uhalf.emplace_back(m_am, p, idx_vec[k], roots[k], ps_idx);
|
|
}
|
|
return roots.size();
|
|
}
|
|
|
|
m_am.isolate_roots(polynomial_ref(p, m_pm), undef_var_assignment(sample(), m_level), roots);
|
|
|
|
// SMT-RAT style reduction: keep only the closest to v roots: one below and one above v,
|
|
// or a single root at v
|
|
unsigned lower = UINT_MAX, upper = UINT_MAX;
|
|
bool section = false;
|
|
for (unsigned k = 0; k < roots.size(); ++k) {
|
|
int cmp = m_am.compare(roots[k], v);
|
|
if (cmp <= 0) {
|
|
lower = k;
|
|
if (cmp == 0) {
|
|
section = true;
|
|
break;
|
|
}
|
|
}
|
|
else {
|
|
upper = k;
|
|
break;
|
|
}
|
|
}
|
|
if (lower != UINT_MAX) {
|
|
lhalf.emplace_back(m_am, p, lower + 1, roots[lower], ps_idx);
|
|
if (section)
|
|
return roots.size(); // only keep the section root for this polynomial
|
|
}
|
|
if (upper != UINT_MAX)
|
|
uhalf.emplace_back(m_am, p, upper + 1, roots[upper], ps_idx);
|
|
return roots.size();
|
|
}
|
|
|
|
void init_poly_has_roots() {
|
|
m_poly_has_roots.clear();
|
|
m_poly_has_roots.resize(m_level_ps.size(), false);
|
|
}
|
|
|
|
bool collect_partitioned_root_functions_around_sample(anum const& v,
|
|
std_vector<root_function>& lhalf, std_vector<root_function>& uhalf) {
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i)
|
|
m_poly_has_roots[i] = isolate_roots_around_sample(m_level_ps.get(i), i, v, lhalf, uhalf);
|
|
return !lhalf.empty() || !uhalf.empty();
|
|
}
|
|
|
|
std_vector<root_function>::iterator set_relation_root_functions_from_partitions(
|
|
std_vector<root_function>& lhalf,
|
|
std_vector<root_function>& uhalf) {
|
|
auto& rfs = m_rel.m_rfunc;
|
|
size_t mid_pos = lhalf.size();
|
|
rfs.clear();
|
|
rfs.reserve(mid_pos + uhalf.size());
|
|
for (auto& rf : lhalf)
|
|
rfs.emplace_back(std::move(rf));
|
|
for (auto& rf : uhalf)
|
|
rfs.emplace_back(std::move(rf));
|
|
return rfs.begin() + mid_pos;
|
|
}
|
|
|
|
bool root_function_lt(root_function const& a, root_function const& b, bool degree_desc) {
|
|
if (a.ire.p == b.ire.p)
|
|
return a.ire.i < b.ire.i;
|
|
auto r = m_am.compare(a.val, b.val);
|
|
if (r)
|
|
return r == sign_neg;
|
|
unsigned dega = m_pm.degree(a.ire.p, m_level);
|
|
unsigned degb = m_pm.degree(b.ire.p, m_level);
|
|
if (dega != degb)
|
|
return degree_desc ? dega > degb : dega < degb;
|
|
return m_pm.id(a.ire.p) < m_pm.id(b.ire.p);
|
|
}
|
|
|
|
void sort_root_function_partitions(std_vector<root_function>::iterator mid) {
|
|
auto& rfs = m_rel.m_rfunc;
|
|
std::sort(rfs.begin(), mid,
|
|
[&](root_function const& a, root_function const& b) { return root_function_lt(a, b, true); });
|
|
std::sort(mid, rfs.end(),
|
|
[&](root_function const& a, root_function const& b) { return root_function_lt(a, b, false); });
|
|
}
|
|
|
|
// Populate Θ (root functions) around the sample, partitioned at `mid`, and sort each partition.
|
|
// Returns whether any roots were found.
|
|
bool build_sorted_root_functions_around_sample(anum const& v, std_vector<root_function>::iterator& mid) {
|
|
init_poly_has_roots();
|
|
|
|
std_vector<root_function> lhalf, uhalf;
|
|
if (!collect_partitioned_root_functions_around_sample(v, lhalf, uhalf))
|
|
return false;
|
|
|
|
mid = set_relation_root_functions_from_partitions(lhalf, uhalf);
|
|
|
|
// 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
|
|
sort_root_function_partitions(mid);
|
|
return true;
|
|
}
|
|
|
|
// Pick I_level around sample(m_level) from sorted Θ, split at `mid`.
|
|
// Sets m_l_rf/m_u_rf (positions in m_rfunc) and m_I[m_level] (interval with root indices).
|
|
void set_interval_from_root_partition(anum const& v, std_vector<root_function>::iterator mid) {
|
|
auto& rfs = m_rel.m_rfunc;
|
|
if (mid != rfs.begin()) {
|
|
m_l_rf = static_cast<unsigned>((mid - rfs.begin()) - 1);
|
|
auto& r = *(mid - 1);
|
|
if (m_am.eq(r.val, v)) {
|
|
// Section case: only m_l_rf is defined
|
|
m_I[m_level].section = true;
|
|
m_I[m_level].l = r.ire.p;
|
|
m_I[m_level].l_index = r.ire.i;
|
|
m_u_rf = m_l_rf;
|
|
}
|
|
else {
|
|
m_I[m_level].l = r.ire.p;
|
|
m_I[m_level].l_index = r.ire.i;
|
|
if (mid != rfs.end()) {
|
|
m_u_rf = m_l_rf + 1;
|
|
m_I[m_level].u = mid->ire.p;
|
|
m_I[m_level].u_index = mid->ire.i;
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
// sample is below all roots: I = (-oo, theta_1)
|
|
m_l_rf = UINT_MAX;
|
|
m_u_rf = 0;
|
|
auto& r = *mid;
|
|
m_I[m_level].u = r.ire.p;
|
|
m_I[m_level].u_index = r.ire.i;
|
|
}
|
|
}
|
|
|
|
// Build Θ (root functions) and pick I_level around sample(level).
|
|
// Sets m_l_rf/m_u_rf and m_I[level].
|
|
// Returns whether any roots were found (i.e., whether a relation can be built).
|
|
bool build_interval() {
|
|
m_rel.clear();
|
|
reset_interval(m_I[m_level]);
|
|
m_l_rf = UINT_MAX;
|
|
m_u_rf = UINT_MAX;
|
|
|
|
anum const& v = sample().value(m_level);
|
|
std_vector<root_function>::iterator mid;
|
|
if (!build_sorted_root_functions_around_sample(v, mid))
|
|
return false;
|
|
|
|
set_interval_from_root_partition(v, mid);
|
|
compute_side_mask();
|
|
return true;
|
|
}
|
|
|
|
|
|
void add_relation_resultants() {
|
|
TRACE(lws, tout << " add_relation_resultants: " << m_rel.m_pairs.size() << " pairs\n";);
|
|
for (auto const& pr : m_rel.m_pairs) {
|
|
poly* p1 = m_level_ps.get(pr.first);
|
|
poly* p2 = m_level_ps.get(pr.second);
|
|
TRACE(lws,
|
|
m_pm.display(m_pm.display(tout << " resultant(" << pr.first << ", " << pr.second << "): ", p1) << " and ", p2) << "\n";
|
|
);
|
|
polynomial_ref res = psc_resultant(p1, p2, m_level);
|
|
TRACE(lws,
|
|
tout << " resultant poly: ";
|
|
if (res)
|
|
m_pm.display(tout, res) << "\n resultant sign at sample: " << m_am.eval_sign_at(res, sample());
|
|
else
|
|
tout << "(null)";
|
|
tout << "\n";
|
|
);
|
|
request_factorized(res, inv_req::ord);
|
|
}
|
|
}
|
|
|
|
// Top level (m_n): add resultants between adjacent roots of different polynomials.
|
|
// Fills m_poly_has_roots as a side effect.
|
|
void add_adjacent_root_resultants() {
|
|
m_poly_has_roots.clear();
|
|
m_poly_has_roots.resize(m_level_ps.size(), false);
|
|
|
|
std_vector<std::pair<scoped_anum, poly*>> root_vals;
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
poly* p = m_level_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);
|
|
m_poly_has_roots[i] = !roots.empty();
|
|
TRACE(lws,
|
|
tout << " poly[" << i << "] has " << roots.size() << " roots: ";
|
|
for (unsigned k = 0; k < roots.size(); ++k) {
|
|
if (k > 0) tout << ", ";
|
|
m_am.display_decimal(tout, roots[k], 5);
|
|
}
|
|
tout << "\n";
|
|
);
|
|
for (unsigned k = 0; k < roots.size(); ++k) {
|
|
scoped_anum root_v(m_am);
|
|
m_am.set(root_v, roots[k]);
|
|
root_vals.emplace_back(std::move(root_v), p);
|
|
}
|
|
}
|
|
if (root_vals.size() < 2)
|
|
return;
|
|
|
|
std::sort(root_vals.begin(), root_vals.end(), [&](auto const& a, auto const& b) {
|
|
return m_am.lt(a.first, b.first);
|
|
});
|
|
|
|
TRACE(lws,
|
|
tout << " Sorted roots:\n";
|
|
for (unsigned j = 0; j < root_vals.size(); ++j)
|
|
m_pm.display(m_am.display_decimal(tout << " [" << j << "] val=", root_vals[j].first, 5) << " poly=", root_vals[j].second) << "\n";
|
|
);
|
|
|
|
std::set<std::pair<poly*, poly*>> added_pairs;
|
|
for (unsigned j = 0; j + 1 < root_vals.size(); ++j) {
|
|
poly* p1 = root_vals[j].second;
|
|
poly* p2 = root_vals[j + 1].second;
|
|
if (!p1 || !p2 || p1 == p2)
|
|
continue;
|
|
if (p1 > p2) std::swap(p1, p2);
|
|
if (!added_pairs.insert({p1, p2}).second)
|
|
continue;
|
|
TRACE(lws,
|
|
m_pm.display(m_pm.display(tout << " Adjacent resultant pair: ", p1) << " and ", p2) << "\n";
|
|
);
|
|
request_factorized(psc_resultant(p1, p2, m_n), inv_req::ord);
|
|
}
|
|
}
|
|
|
|
void upgrade_bounds_to_ord() {
|
|
poly* lb = m_I[m_level].l;
|
|
poly* ub = m_I[m_level].u;
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
poly* p = m_level_ps.get(i);
|
|
if (p == lb || p == ub)
|
|
set_req(i, inv_req::ord);
|
|
}
|
|
}
|
|
|
|
bool is_section() { return m_I[m_level].is_section();}
|
|
|
|
// Common projection loop used by both section and sector cases.
|
|
// m_omit_lc and m_omit_disc are precomputed (empty for section case).
|
|
void add_level_projections() {
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
polynomial_ref p(m_level_ps.get(i), m_pm);
|
|
polynomial_ref lc(m_pm);
|
|
unsigned deg = m_pm.degree(p, m_level);
|
|
lc = m_pm.coeff(p, m_level, deg);
|
|
|
|
bool add_lc = !vec_get(m_omit_lc, i, false);
|
|
bool add_disc = !vec_get(m_omit_disc, i, false);
|
|
|
|
TRACE(lws,
|
|
tout << " poly[" << i << "]";
|
|
tout << " omit_lc=" << vec_get(m_omit_lc, i, false);
|
|
tout << " omit_disc=" << vec_get(m_omit_disc, i, false);
|
|
tout << "\n";
|
|
);
|
|
|
|
// SMT-RAT-style coeffNonNull: if the leading coefficient is already non-zero at the sample
|
|
// AND we're adding lc, we do not need to project an additional non-null coefficient witness.
|
|
polynomial_ref witness = m_witnesses[i];
|
|
if (add_lc && witness && !is_const(witness))
|
|
if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0)
|
|
witness = polynomial_ref(m_pm);
|
|
|
|
add_projection_for_poly(p, m_level, witness, add_lc, add_disc);
|
|
}
|
|
}
|
|
|
|
// Check if spanning tree should be used based on both_count threshold.
|
|
// Note: This is only called from process_level_sector() which handles non-section cases.
|
|
// Spanning tree requires: distinct lower/upper bounds (!same_boundary_poly) and enough both-side polys.
|
|
bool should_use_spanning_tree() {
|
|
SASSERT(!is_section()); // spanning_tree is for sector case only
|
|
// Need at least 2 polynomials for a spanning tree to have edges
|
|
if (m_spanning_tree_threshold < 2)
|
|
return false;
|
|
if (m_rel.m_rfunc.size() < 4) // 2 different bounds + at least 2 same out of bounds
|
|
return false;
|
|
if (m_side_mask.size() == 0)
|
|
return false;
|
|
|
|
if (!is_set(m_l_rf) || !is_set(m_u_rf))
|
|
return false;
|
|
|
|
if (same_boundary_poly())
|
|
return false; // spanning tree requires distinct lower/upper bounds
|
|
|
|
unsigned both_count = 0;
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i)
|
|
if (m_side_mask[i] == 3)
|
|
if (++both_count >= m_spanning_tree_threshold)
|
|
return true;
|
|
return false;
|
|
}
|
|
|
|
// Clear all per-level state that could be stale from previous levels.
|
|
// This ensures no leftover heuristic decisions affect the current level.
|
|
void clear_level_state() {
|
|
m_omit_lc.clear();
|
|
m_omit_disc.clear();
|
|
m_rel.m_pairs.clear();
|
|
m_side_mask.clear();
|
|
m_deg_in_order_graph.clear();
|
|
m_unique_neighbor.clear();
|
|
}
|
|
|
|
void process_level_section(bool have_interval) {
|
|
SASSERT(m_I[m_level].section);
|
|
m_rel_mode = biggest_cell; // default
|
|
if (have_interval) {
|
|
fill_relation_pairs_for_section_biggest_cell();
|
|
SASSERT(relation_invariant());
|
|
}
|
|
upgrade_bounds_to_ord();
|
|
// For section case, m_omit_lc and m_omit_disc are empty (cleared by clear_level_state)
|
|
// so all LCs and discriminants are added
|
|
add_level_projections();
|
|
add_relation_resultants();
|
|
}
|
|
|
|
void process_level_sector(bool have_interval) {
|
|
SASSERT(!m_I[m_level].section);
|
|
m_rel_mode = biggest_cell; // default
|
|
if (have_interval) {
|
|
// Check spanning tree threshold first
|
|
if (should_use_spanning_tree()) {
|
|
fill_relation_sector_spanning_tree();
|
|
compute_omit_disc_for_spanning_tree();
|
|
m_rel_mode = spanning_tree;
|
|
} else {
|
|
fill_relation_sector_biggest_cell();
|
|
// For same_boundary_poly, compute disc omission
|
|
if (same_boundary_poly())
|
|
compute_omit_disc_for_same_boundary();
|
|
// TODO: Could also drop discriminants for single-side leaves in biggest_cell mode.
|
|
// A leaf p with all roots on one side (e.g., below lb) has sign at sample determined
|
|
// solely by being above/below all roots - internal root ordering doesn't matter.
|
|
// SMT-RAT doesn't implement this, but it should be theoretically sound.
|
|
}
|
|
compute_omit_lc_both_sides(m_rel_mode == spanning_tree);
|
|
|
|
SASSERT(relation_invariant());
|
|
}
|
|
upgrade_bounds_to_ord();
|
|
add_level_projections();
|
|
add_relation_resultants();
|
|
}
|
|
|
|
void collect_non_null_witnesses() {
|
|
// Line 10/11: detect nullification + pick a non-zero coefficient witness per p.
|
|
m_witnesses.clear();
|
|
m_witnesses.reserve(m_level_ps.size());
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
polynomial_ref p(m_level_ps.get(i), m_pm);
|
|
polynomial_ref w = choose_nonzero_coeff(p, m_level);
|
|
if (!w)
|
|
fail();
|
|
m_witnesses.push_back(w);
|
|
}
|
|
}
|
|
|
|
void display_polys_at_level(std::ostream& out) {
|
|
TRACE(lws, out << "Polynomials at level " << m_level << "\n";
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i)
|
|
m_pm.display(out, m_level_ps.get(i)) << "\n";);
|
|
}
|
|
|
|
void process_level() {
|
|
TRACE(lws, display_polys_at_level(tout););
|
|
clear_level_state(); // Clear stale state from previous level
|
|
collect_non_null_witnesses();
|
|
// Lines 3-8: Θ + I_level + relation ≼
|
|
bool have_interval = build_interval();
|
|
|
|
TRACE(lws,
|
|
display(tout << "Interval: ", m_solver, m_I[m_level])
|
|
<< "\nSection? " << (m_I[m_level].section ? "yes" : "no")
|
|
<< "\nhave_interval=" << have_interval << ", rfunc.size=" << m_rel.m_rfunc.size() << "\n";
|
|
for (unsigned i = 0; i < m_rel.m_rfunc.size(); ++i)
|
|
m_am.display(tout << " rfunc[" << i << "]: ps_idx=" << m_rel.m_rfunc[i].ps_idx << ", val=", m_rel.m_rfunc[i].val) << "\n";
|
|
);
|
|
|
|
if (m_I[m_level].section)
|
|
process_level_section(have_interval);
|
|
else
|
|
process_level_sector(have_interval);
|
|
}
|
|
|
|
bool poly_has_roots(unsigned i) { return i < vec_get(m_poly_has_roots, i, false); }
|
|
|
|
void process_top_level() {
|
|
TRACE(lws, display_polys_at_level(tout););
|
|
collect_non_null_witnesses();
|
|
add_adjacent_root_resultants();
|
|
|
|
// Projections (coeff witness, disc, leading coeff).
|
|
for (unsigned i = 0; i < m_level_ps.size(); ++i) {
|
|
polynomial_ref p(m_level_ps.get(i), m_pm);
|
|
polynomial_ref lc(m_pm);
|
|
unsigned deg = m_pm.degree(p, m_n);
|
|
lc = m_pm.coeff(p, m_n, deg);
|
|
|
|
bool add_lc = true;
|
|
if (!poly_has_roots(i))
|
|
if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0)
|
|
add_lc = false;
|
|
|
|
// if the leading coefficient is already non-zero at the sample
|
|
// AND we're adding lc, we do not need to project an additional non-null coefficient witness.
|
|
polynomial_ref witness = m_witnesses[i];
|
|
if (add_lc && witness && !is_const(witness))
|
|
if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0)
|
|
witness = polynomial_ref(m_pm); // zero the witnsee as lc will be the witness
|
|
add_projection_for_poly(p, m_n, witness, add_lc, true); //true to add the discriminant
|
|
}
|
|
}
|
|
|
|
std_vector<root_function_interval> single_cell_work() {
|
|
TRACE(lws,
|
|
tout << "Input polynomials (" << m_P.size() << "):\n";
|
|
for (unsigned i = 0; i < m_P.size(); ++i)
|
|
m_pm.display(tout << " p[" << i << "]: ", m_P.get(i)) << "\n";
|
|
tout << "Sample values:\n";
|
|
for (unsigned j = 0; j < m_n; ++j)
|
|
m_am.display(tout << " x" << j << " = ", sample().value(j)) << "\n";
|
|
);
|
|
|
|
if (m_n == 0) return m_I;
|
|
|
|
m_todo.reset();
|
|
|
|
// Initialize m_todo with distinct irreducible factors of the input polynomials.
|
|
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) {
|
|
request(f.get(), inv_req::sign);
|
|
});
|
|
}
|
|
|
|
if (m_todo.empty()) return m_I;
|
|
|
|
// Process top level m_n (we project from x_{m_n} and do not return an interval for it).
|
|
if (m_todo.max_var() == m_n) {
|
|
extract_max_tagged();
|
|
process_top_level();
|
|
}
|
|
|
|
// Now iteratively process remaining levels (descending).
|
|
while (!m_todo.empty()) {
|
|
extract_max_tagged();
|
|
SASSERT(m_level < m_n);
|
|
process_level();
|
|
TRACE(lws,
|
|
tout << "After level " << m_level << ": m_todo.empty()=" << m_todo.empty();
|
|
if (!m_todo.empty()) tout << ", m_todo.max_var()=" << m_todo.max_var();
|
|
tout << "\n";
|
|
);
|
|
}
|
|
|
|
TRACE(lws,
|
|
for (unsigned i = 0; i < m_I.size(); ++i)
|
|
display(tout << "I[" << i << "]: ", m_solver, m_I[i]) << "\n";
|
|
);
|
|
|
|
return m_I;
|
|
}
|
|
|
|
std_vector<root_function_interval> single_cell() {
|
|
try {
|
|
return single_cell_work();
|
|
}
|
|
catch (nullified_poly_exception&) {
|
|
m_fail = true;
|
|
return m_I;
|
|
}
|
|
}
|
|
};
|
|
|
|
levelwise::levelwise(
|
|
nlsat::solver& solver,
|
|
polynomial_ref_vector const& ps,
|
|
var n,
|
|
assignment const& s,
|
|
pmanager& pm,
|
|
anum_manager& am,
|
|
polynomial::cache& cache)
|
|
: m_impl(new impl(solver, ps, n, s, pm, am, cache)) {}
|
|
|
|
levelwise::~levelwise() { delete m_impl; }
|
|
|
|
std_vector<levelwise::root_function_interval> levelwise::single_cell() {
|
|
return m_impl->single_cell();
|
|
}
|
|
|
|
bool levelwise::failed() const { return m_impl->m_fail; }
|
|
|
|
} // namespace nlsat
|
|
|
|
// Free pretty-printer for symbolic_interval
|
|
std::ostream& nlsat::display(std::ostream& out, solver& s, levelwise::root_function_interval const& I) {
|
|
if (I.section) {
|
|
out << "Section: ";
|
|
if (I.l == nullptr)
|
|
out << "(undef)";
|
|
else {
|
|
::nlsat::display(out, s, I.l);
|
|
out << "[root_index:" << I.l_index << "]";
|
|
}
|
|
}
|
|
else {
|
|
out << "Sector: (";
|
|
if (I.l_inf())
|
|
out << "-oo";
|
|
else {
|
|
::nlsat::display(out, s, I.l);
|
|
out << "[root_index:" << I.l_index << "]";
|
|
}
|
|
out << ", ";
|
|
if (I.u_inf())
|
|
out << "+oo";
|
|
else {
|
|
::nlsat::display(out, s, I.u);
|
|
out << "[root_index:" << I.u_index << "]";
|
|
}
|
|
out << ")";
|
|
}
|
|
return out;
|
|
}
|