diff --git a/.gitattributes b/.gitattributes index 647728dcc..3fcf2b5d2 100644 --- a/.gitattributes +++ b/.gitattributes @@ -3,4 +3,7 @@ src/api/dotnet/Properties/AssemblyInfo.cs text eol=crlf -.github/workflows/*.lock.yml linguist-generated=true merge=ours \ No newline at end of file +.github/workflows/*.lock.yml linguist-generated=true merge=ours + +# Use bd merge for beads JSONL files +.beads/issues.jsonl merge=beads diff --git a/.gitignore b/.gitignore index 8104503a8..ee5df58ff 100644 --- a/.gitignore +++ b/.gitignore @@ -115,3 +115,5 @@ genaisrc/genblogpost.genai.mts *.mts # Bazel generated files bazel-* +# Local issue tracking +.beads diff --git a/levelwise.pdf b/levelwise.pdf new file mode 100644 index 000000000..ecfb982f5 Binary files /dev/null and b/levelwise.pdf differ diff --git a/src/math/polynomial/algebraic_numbers.cpp b/src/math/polynomial/algebraic_numbers.cpp index 24cb3c782..7c72ffd63 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/math/polynomial/polynomial_cache.cpp b/src/math/polynomial/polynomial_cache.cpp index ef407023d..7d3a15ded 100644 --- a/src/math/polynomial/polynomial_cache.cpp +++ b/src/math/polynomial/polynomial_cache.cpp @@ -128,7 +128,7 @@ namespace polynomial { m_factor_cache.reset(); } - unsigned pid(polynomial * p) const { return m.id(p); } + unsigned pid(const polynomial * p) const { return m.id(p); } polynomial * mk_unique(polynomial * p) { if (m_in_cache.get(pid(p), false)) @@ -141,6 +141,28 @@ namespace polynomial { return p_prime; } + bool contains(const polynomial * p) const { + return m_in_cache.get(pid(p), false); + } + + bool contains_chain(polynomial * p, polynomial * q, var x) const { + if (!m_in_cache.get(pid(p), false)) { + polynomial * const * p2 = m_poly_table.find_core(p); + if (!p2) + return false; + p = *p2; + } + if (!m_in_cache.get(pid(q), false)) { + polynomial * const * q2 = m_poly_table.find_core(q); + if (!q2) + return false; + q = *q2; + } + unsigned h = hash_u_u(pid(p), pid(q)); + psc_chain_entry key(p, q, x, h); + return m_psc_chain_cache.contains(&key); + } + void psc_chain(polynomial * p, polynomial * q, var x, polynomial_ref_vector & S) { p = mk_unique(p); q = mk_unique(q); @@ -213,6 +235,14 @@ namespace polynomial { return m_imp->mk_unique(p); } + bool cache::contains(const polynomial * p) const { + return m_imp->contains(p); + } + + bool cache::contains_chain(polynomial const * p, polynomial const * q, var x) const { + return m_imp->contains_chain(const_cast(p), const_cast(q), x); + } + void cache::psc_chain(polynomial const * p, polynomial const * q, var x, polynomial_ref_vector & S) { m_imp->psc_chain(const_cast(p), const_cast(q), x, S); } diff --git a/src/math/polynomial/polynomial_cache.h b/src/math/polynomial/polynomial_cache.h index fbd70e500..a27abf42d 100644 --- a/src/math/polynomial/polynomial_cache.h +++ b/src/math/polynomial/polynomial_cache.h @@ -34,9 +34,10 @@ namespace polynomial { manager & m() const; manager & pm() const { return m(); } polynomial * mk_unique(polynomial * p); + bool contains(const polynomial * p) const; + bool contains_chain(polynomial const * p, polynomial const * q, var x) const; void psc_chain(polynomial const * p, polynomial const * q, var x, polynomial_ref_vector & S); void factor(polynomial const * p, polynomial_ref_vector & distinct_factors); void reset(); }; }; - diff --git a/src/nlsat/CMakeLists.txt b/src/nlsat/CMakeLists.txt index 077c420f9..7faae8b19 100644 --- a/src/nlsat/CMakeLists.txt +++ b/src/nlsat/CMakeLists.txt @@ -1,9 +1,11 @@ z3_add_component(nlsat SOURCES nlsat_clause.cpp + nlsat_common.cpp nlsat_evaluator.cpp nlsat_explain.cpp nlsat_interval_set.cpp + levelwise.cpp nlsat_simplify.cpp nlsat_solver.cpp nlsat_types.cpp diff --git a/src/nlsat/levelwise.cpp b/src/nlsat/levelwise.cpp new file mode 100644 index 000000000..1a6583c9f --- /dev/null +++ b/src/nlsat/levelwise.cpp @@ -0,0 +1,1428 @@ +#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 +#include +#include +#include +#include +#include + +static bool is_set(unsigned k) { return static_cast(k) != -1; } + +// Helper for sparse vector access with default value +template +static T vec_get(std_vector const& v, unsigned idx, T def) { + return idx < v.size() ? v[idx] : def; +} + +// Helper for sparse vector write with auto-resize +template +static void vec_setx(std_vector& v, unsigned idx, T val, T def) { + if (idx >= v.size()) + v.resize(idx + 1, def); + v[idx] = val; +} + +namespace nlsat { + + // The three projection modes for a level: + // 1. section_biggest_cell: Sample is on a root. All disc/lc added. + // 2. sector_biggest_cell: Sample between roots. noLdcf optimization only. + // 3. sector_spanning_tree: Sample between roots with many both-side polys. + // Both noLdcf and noDisc optimizations. + enum class projection_mode { + section_biggest_cell, + sector_biggest_cell, + sector_spanning_tree + }; + +// Tag indicating what kind of invariance we need to preserve for a polynomial on the cell: +// - SIGN: sign-invariance is sufficient (polynomial doesn't change sign within the cell) +// - ORD: order-invariance is required (root multiplicities are constant within the cell) +// +// Order-invariance is stronger than sign-invariance and is needed for: +// - Discriminants: to ensure root functions remain continuous and ordered (Theorem 2.1 in [1]) +// - Resultants: to ensure relative ordering of roots from different polynomials (Theorem 2.2 in [1]) +// Sign-invariance suffices for leading coefficients (ensures polynomial degree doesn't drop). +// +// [1] Nalbach et al., "Projective Delineability for Single Cell Construction", SC² 2025 + 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 { + 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 m_I; // intervals for variables 0..m_n-1 + + unsigned m_level = 0; // current level being processed + 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 m_witnesses; + std_vector 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 m_req; + + // Vectors indexed by position in m_level_ps (more cache-friendly than maps) + mutable std_vector m_side_mask; // bit0 = lower, bit1 = upper, 3 = both + mutable std_vector m_omit_lc; + mutable std_vector m_omit_disc; + mutable std_vector m_deg_in_order_graph; // degree of polynomial in resultant graph + mutable std_vector 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 + 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 m_rfunc; // Θ: root functions on current level + std::set> 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> 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 { + std_vector visited(max_ps_idx + 1, false); + std_vector 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 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 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(vec_get(m_req, id, static_cast(inv_req::none))); + inv_req nxt = max_req(cur, req); + if (nxt != cur) + vec_setx(m_req, id, static_cast(nxt), static_cast(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"; + ); + // Each irreducible factor inherits the invariance requirement. + // If already requested with a weaker tag, upgrade to the stronger one. + request(f.get(), req); + }); + } + + // 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(inv_req::none)) == static_cast(inv_req::none)) + vec_setx(m_req, id, static_cast(inv_req::sign), static_cast(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(vec_get(m_req, id, static_cast(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(req), static_cast(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 optimization - compute which leading coefficients can be omitted + // + // This optimization is justified by PROJECTIVE DELINEABILITY theory [1,2]. + // + // Regular delineability (Theorem 2.1 in [2]) requires the leading coefficient + // to be sign-invariant to ensure the polynomial's degree doesn't drop. However, + // projective delineability (Theorem 3.1 in [2]) shows that order-invariance of + // the discriminant alone (without LC sign-invariance) is sufficient to ensure + // the polynomial's roots behave continuously - they may "go to infinity" when + // the LC vanishes, but won't suddenly appear or disappear within the cell. + // + // For a polynomial p with roots on BOTH sides of the sample: + // - Resultants with both boundary polynomials are computed anyway + // - These resultants ensure p's roots don't cross the cell boundaries (Thm 3.2) + // - Even if p's degree drops (LC becomes zero), existing roots remain ordered + // - New roots can only appear "at infinity", outside the bounded cell + // + // [1] Michel et al., "On Projective Delineability", arXiv:2411.13300, 2024 + // [2] Nalbach et al., "Projective Delineability for Single Cell Construction", SC² 2025 + // ============================================================================ + + // 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, or 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, or 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); + } + } + + // TODO: Investigate noDisc optimization for spanning_tree mode when lb and ub are the same + // polynomial. SMT-RAT does noDisc for leaves connected to the section-defining polynomial + // in section case, which might extend to sectors where lb == ub. The previous implementation + // here was unsound because it applied to general leaves connected to any boundary, but + // discriminants are needed to ensure a polynomial's own roots don't merge/split. + + // ---------------------------------------------------------------------------- + // noLdcf heuristic helpers + // ---------------------------------------------------------------------------- + + // Compute noLdcf: which leading coefficients can be omitted using projective delineability. + // + // A polynomial p with roots on BOTH sides of the sample has resultants with both bounds. + // By projective delineability (Theorem 3.1 in [2]), we only need the discriminant to be + // order-invariant - the LC can be omitted because roots "going to infinity" don't affect + // sign-invariance within the bounded cell. + // + // - biggest_cell mode, require_leaf=false: all non-bound polys with both-side roots + // - spanning_tree mode, require_leaf=true: only leaves with deg == 1 and both-side roots + // + // [2] Nalbach et al., "Projective Delineability for Single Cell Construction", SC² 2025 + 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; + + // 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 lower_rf(m_level_ps.size(), UINT_MAX); + std_vector 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; + (void)upper_root_idx; // used in DEBUG_CODE below + + // 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 (used in debug via SASSERT) + DEBUG_CODE( + unsigned lower_root_idx = both.size() - 1; + auto arb_invariant = [&]() { + // Reconstruct parent[] from the algorithm logic + std_vector 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> 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) + out_children[i].push_back(p); // i has smaller upper_rf + else + out_children[p].push_back(i); // p has smaller upper_rf + } + } + + // 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 reached(both.size(), false); + std_vector 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& lhalf, + std_vector& 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 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); + + // For single cell construction, we only need the closest roots to the sample point: + // at most one root below v and one root above v, or a single root at v for sections. + // Other roots are irrelevant for bounding the cell containing the sample. + 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& lhalf, std_vector& 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::iterator set_relation_root_functions_from_partitions( + std_vector& lhalf, + std_vector& 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::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::iterator& mid) { + init_poly_has_roots(); + + std_vector 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::iterator mid) { + auto& rfs = m_rel.m_rfunc; + if (mid != rfs.begin()) { + m_l_rf = static_cast((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::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> 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> 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();} + + // Section case: the section-defining polynomial needs disc and lc projections. + // For polynomials WITH roots: resultants with section polynomial ensure root ordering. + // For polynomials WITHOUT roots: they need LC + disc to ensure they don't gain roots. + // + // Theory: For a section (sample on a root), sign-invariance of polynomials with roots + // is ensured by resultants with the section-defining polynomial (Thm 2.2 in [1]). + // But polynomials without roots need delineability (LC + disc) to stay root-free. + // + // [1] Nalbach et al., "Projective Delineability for Single Cell Construction", SC² 2025 + void add_section_projections() { + poly* section_poly = m_I[m_level].l; + SASSERT(section_poly); + + // Build set of polynomial indices that have roots at this level + std::set has_roots; + for (auto const& rf : m_rel.m_rfunc) + has_roots.insert(rf.ps_idx); + + for (unsigned i = 0; i < m_level_ps.size(); ++i) { + polynomial_ref p(m_level_ps.get(i), m_pm); + polynomial_ref witness = m_witnesses[i]; + + if (m_level_ps.get(i) == section_poly) + add_projection_for_poly(p, m_level, witness, true, true); // section poly: full projection + else if (has_roots.find(i) == has_roots.end()) + add_projection_for_poly(p, m_level, witness, true, true); // no roots: need LC+disc for delineability + else if (witness && !is_const(witness)) + request_factorized(witness, inv_req::sign); // has roots: witness only + } + } + + // Sector case: projection loop using m_omit_lc and m_omit_disc. + void add_sector_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"; + ); + + // coeffNonNull optimization: if the leading coefficient is already non-zero at the + // sample point AND we're projecting it, the LC itself serves as the non-null witness. + // No need for an additional coefficient witness in this case. + 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 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_with_mode(projection_mode mode, bool have_interval) { + if (have_interval) { + switch (mode) { + case projection_mode::section_biggest_cell: + fill_relation_pairs_for_section_biggest_cell(); + break; + + case projection_mode::sector_biggest_cell: + fill_relation_sector_biggest_cell(); + compute_side_mask(); + compute_omit_lc_both_sides(/*require_leaf=*/false); + // m_omit_disc stays empty - all discriminants added + break; + + case projection_mode::sector_spanning_tree: + fill_relation_sector_spanning_tree(); + compute_side_mask(); + compute_omit_lc_both_sides(/*require_leaf=*/true); + // m_omit_disc stays empty - all discriminants added + break; + } + SASSERT(relation_invariant()); + } + + upgrade_bounds_to_ord(); + + if (mode == projection_mode::section_biggest_cell) + add_section_projections(); + else + add_sector_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"; + ); + + projection_mode mode; + if (m_I[m_level].section) + mode = projection_mode::section_biggest_cell; + else if (should_use_spanning_tree()) + mode = projection_mode::sector_spanning_tree; + else + mode = projection_mode::sector_biggest_cell; + + process_level_with_mode(mode, have_interval); + } + + bool poly_has_roots(unsigned i) { return 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 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 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::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; +} diff --git a/src/nlsat/levelwise.h b/src/nlsat/levelwise.h new file mode 100644 index 000000000..27e379374 --- /dev/null +++ b/src/nlsat/levelwise.h @@ -0,0 +1,56 @@ +#pragma once +#include "nlsat_types.h" +#include "math/polynomial/polynomial_cache.h" + +namespace nlsat { + + class assignment; // forward declared in nlsat_types.h + + class levelwise { + public: + struct indexed_root_expr { + poly* p; + unsigned i; + }; + struct root_function_interval { + bool section = false; + polynomial_ref l; + unsigned l_index; // the low bound root index + polynomial_ref u; + unsigned u_index; // the upper bound root index + bool l_inf() const { return l == nullptr; } + bool u_inf() const { return u == nullptr; } + bool is_section() const { return section; } + bool is_sector() const { return !section; } + polynomial_ref& section_poly() { + SASSERT(is_section()); + return l; + } + root_function_interval(polynomial::manager & pm):l(pm), u(pm) {} + }; + // Free pretty-printer declared here so external modules (e.g., nlsat_explain) can + // display intervals without depending on levelwise internals. + // Implemented in levelwise.cpp + friend std::ostream& display(std::ostream& out, solver& s, root_function_interval const& I); + + private: + + + struct impl; + impl* m_impl; + public: + // Construct with polynomials ps, maximal variable max_x, current sample s, polynomial manager pm, and algebraic-number manager am + levelwise(nlsat::solver& solver, polynomial_ref_vector const& ps, var max_x, assignment const& s, pmanager& pm, anum_manager& am, polynomial::cache & cache); + ~levelwise(); + + levelwise(levelwise const&) = delete; + levelwise& operator=(levelwise const&) = delete; + std_vector single_cell(); + bool failed() const; + }; + + // + // Free pretty-printer (non-member) for levelwise::symbolic_interval + std::ostream& display(std::ostream& out, solver& s, levelwise::root_function_interval const& I); + +} // namespace nlsat diff --git a/src/nlsat/nlsat_common.cpp b/src/nlsat/nlsat_common.cpp new file mode 100644 index 000000000..547d2660d --- /dev/null +++ b/src/nlsat/nlsat_common.cpp @@ -0,0 +1,123 @@ +/*++ +Copyright (c) 2012 Microsoft Corporation + +Module Name: + + nlsat_common.cpp + +Abstract: + + some common routines from nlsat + +Author: + + Lev Nachmanson(levnach@hotmail.com) 2025-October. + +Revision History: + +--*/ +#include "nlsat/nlsat_common.h" +namespace nlsat { + + todo_set::todo_set(polynomial::cache& u, bool canonicalize): m_cache(u), m_set(u.pm()), m_canonicalize(canonicalize) {} + + void todo_set::reset() { + pmanager& pm = m_set.m(); + unsigned sz = m_set.size(); + for (unsigned i = 0; i < sz; i++) { + m_in_set[pm.id(m_set.get(i))] = false; + } + m_set.reset(); + } + + poly* todo_set::insert(poly* p) { + pmanager& pm = m_set.m(); + if (m_in_set.get(pm.id(p), false)) + return p; + if (m_cache.contains(p)) { + // still have to insert in the set + m_in_set.setx(pm.id(p), true, false); + m_set.push_back(p); + return p; + } + polynomial_ref pinned(pm); // keep canonicalized polynomial alive until it is stored + if (m_canonicalize) { + // Canonicalize content+sign so scalar multiples share the same representative. + if (!pm.is_zero(p) && !pm.is_const(p)) { + var x = pm.max_var(p); + pm.primitive(p, x, pinned); + p = pinned.get(); + } + else + pinned = p; + p = pm.flip_sign_if_lm_neg(p); + pinned = p; + } + p = m_cache.mk_unique(p); + unsigned pid = pm.id(p); + if (!m_in_set.get(pid, false)) { + m_in_set.setx(pid, true, false); + m_set.push_back(p); + } + return p; + } + + bool todo_set::contains(poly* p) const { + if (!p) + return false; + pmanager& pm = m_set.m(); + return m_in_set.get(pm.id(p), false); + } + + bool todo_set::empty() const { return m_set.empty(); } + + // Return max variable in todo_set + var todo_set::max_var() const { + pmanager& pm = m_set.m(); + var max = null_var; + unsigned sz = m_set.size(); + for (unsigned i = 0; i < sz; i++) { + var x = pm.max_var(m_set.get(i)); + SASSERT(x != null_var); + if (max == null_var || x > max) + max = x; + } + return max; + } + + /** + \brief Remove the maximal polynomials from the set and store + them in max_polys. Return the maximal variable + */ + var todo_set::extract_max_polys(polynomial_ref_vector& max_polys) { + max_polys.reset(); + var x = max_var(); + pmanager& pm = m_set.m(); + unsigned sz = m_set.size(); + unsigned j = 0; + for (unsigned i = 0; i < sz; i++) { + poly* p = m_set.get(i); + var y = pm.max_var(p); + SASSERT(y <= x); + if (y == x) { + max_polys.push_back(p); + m_in_set[pm.id(p)] = false; + } + else { + m_set.set(j, p); + j++; + } + } + m_set.shrink(j); + return x; + } + + /** + \brief Wrapper for factorization + */ + void factor(polynomial_ref & p, polynomial::cache& cache, polynomial_ref_vector & fs) { + TRACE(nlsat_factor, tout << "factor\n" << p << "\n";); + fs.reset(); + cache.factor(p.get(), fs); + } +} diff --git a/src/nlsat/nlsat_common.h b/src/nlsat/nlsat_common.h new file mode 100644 index 000000000..c7864e991 --- /dev/null +++ b/src/nlsat/nlsat_common.h @@ -0,0 +1,156 @@ +/*++ + Copyright (c) 2025 + + Module Name: + + nlsat_common.h + + Abstract: + + Pretty-print helpers for NLSAT components as free functions. + These forward to existing solver/pmanager display facilities. + + --*/ +#pragma once + +#include "nlsat/nlsat_assignment.h" +#include "nlsat/nlsat_solver.h" +#include "nlsat/nlsat_scoped_literal_vector.h" +#include "math/polynomial/polynomial_cache.h" + +namespace nlsat { + + // Lightweight set of polynomials that keeps elements unique (by cache id) and + // supports extracting all polynomials whose maximal variable is maximal. + // Optional canonicalization (primitive + sign) can be enabled per instance. + struct todo_set { + polynomial::cache& m_cache; + polynomial_ref_vector m_set; + svector m_in_set; + bool m_canonicalize; + + todo_set(polynomial::cache& u, bool canonicalize); + + void reset(); + // Insert polynomial (canonicalizing if requested) and return the cached representative. + poly* insert(poly* p); + bool contains(poly* p) const; + bool contains(polynomial_ref const& p) const { return contains(p.get()); } + bool empty() const; + // Return max variable in todo_set + var max_var() const; + /** + \brief Remove the maximal polynomials from the set and store + them in max_polys. Return the maximal variable + */ + var extract_max_polys(polynomial_ref_vector& max_polys); + }; + + inline std::ostream& display(std::ostream& out, pmanager& pm, polynomial_ref const& p, display_var_proc const& proc) { + pm.display(out, p, proc); + return out; + } + + inline std::ostream& display(std::ostream& out, pmanager& pm, polynomial_ref_vector const& ps, display_var_proc const& proc, char const* delim = "\n") { + for (unsigned i = 0; i < ps.size(); ++i) { + if (i > 0) out << delim; + pm.display(out, ps.get(i), proc); + } + return out; + } + + inline std::ostream& display(std::ostream& out, solver& s, polynomial_ref const& p) { + return display(out, s.pm(), p, s.display_proc()); + } + + inline std::ostream& display(std::ostream& out, solver& s, poly* p) { + return display(out, s.pm(), polynomial_ref(p, s.pm()), s.display_proc()); + } + + inline std::ostream& display(std::ostream& out, solver& s, polynomial_ref_vector const& ps, char const* delim = "\n") { + return display(out, s.pm(), ps, s.display_proc(), delim); + } + + inline std::ostream& display(std::ostream& out, solver& s, literal l) { + return s.display(out, l); + } + + inline std::ostream& display(std::ostream& out, solver& s, unsigned n, literal const* ls) { + return s.display(out, n, ls); + } + + inline std::ostream& display(std::ostream& out, solver& s, literal_vector const& ls) { + return s.display(out, ls); + } + + inline std::ostream& display(std::ostream& out, solver& s, scoped_literal_vector const& ls) { + return s.display(out, ls.size(), ls.data()); + } + + inline std::ostream& display_var(std::ostream& out, solver& s, var x) { + return s.display(out, x); + } + /** + \brief evaluate the given polynomial in the current interpretation. + max_var(p) must be assigned in the current interpretation. + */ + inline ::sign sign(polynomial_ref const & p, assignment & x2v, anum_manager& am) { + SASSERT(max_var(p) == null_var || x2v.is_assigned(max_var(p))); + auto s = am.eval_sign_at(p, x2v); + return s; + } + + /** + * Check whether all coefficients of the polynomial `s` (viewed as a polynomial + * in its main variable) evaluate to zero under the given assignment `x2v`. + * This is exactly the logic used in several places in the nlsat codebase + * (e.g. coeffs_are_zeroes_in_factor in nlsat_explain.cpp). + */ + inline bool coeffs_are_zeroes_on_sample(polynomial_ref const & s, pmanager & pm, assignment & x2v, anum_manager & am) { + polynomial_ref c(pm); + var x = pm.max_var(s); + unsigned n = pm.degree(s, x); + for (unsigned j = 0; j <= n; ++j) { + c = pm.coeff(s, x, j); + if (sign(c, x2v, am) != 0) + return false; + } + return true; + } + + /** + * \brief Display a vector of algebraic numbers in several commonly useful formats. + * + * This mirrors the ad-hoc helper that existed in `src/test/algebraic.cpp` so that + * solver / explanation code can conveniently dump root sets while debugging. + * + * For each algebraic number it prints (in order): + * - a decimal approximation (10 digits) + * - the root object representation (defining polynomial & isolating interval) + * - the isolating interval alone + * + */ + inline void display(std::ostream & out, scoped_anum_vector const & rs) { + algebraic_numbers::manager & m = rs.m(); + out << "numbers in decimal:\n"; + for (const auto & r : rs) { + m.display_decimal(out, r, 10); + out << '\n'; + } + out << "numbers as root objects\n"; + for (const auto & r : rs) { + m.display_root(out, r); + out << '\n'; + } + out << "numbers as intervals\n"; + for (const auto & r : rs) { + m.display_interval(out, r); + out << '\n'; + } + } + /** + \brief Wrapper for factorization + */ + void factor(polynomial_ref & p, polynomial::cache& cache, polynomial_ref_vector & fs); + +} // namespace nlsat diff --git a/src/nlsat/nlsat_explain.cpp b/src/nlsat/nlsat_explain.cpp index bceb46dc0..714e49c3f 100644 --- a/src/nlsat/nlsat_explain.cpp +++ b/src/nlsat/nlsat_explain.cpp @@ -17,9 +17,11 @@ Revision History: --*/ #include "nlsat/nlsat_explain.h" +#include "nlsat/levelwise.h" #include "nlsat/nlsat_assignment.h" #include "nlsat/nlsat_evaluator.h" #include "math/polynomial/algebraic_numbers.h" +#include "nlsat/nlsat_common.h" #include "util/ref_buffer.h" #include "util/mpq.h" @@ -32,7 +34,6 @@ namespace nlsat { struct explain::imp { solver & m_solver; - assignment const & m_assignment; atom_vector const & m_atoms; atom_vector const & m_x2eq; anum_manager & m_am; @@ -49,85 +50,25 @@ namespace nlsat { bool m_factor; bool m_add_all_coeffs; bool m_add_zero_disc; - bool m_cell_sample; + assignment const & sample() const { return m_solver.sample(); } + assignment & sample() { return m_solver.sample(); } - struct todo_set { - polynomial::cache & m_cache; - polynomial_ref_vector m_set; - svector m_in_set; - - todo_set(polynomial::cache & u):m_cache(u), m_set(u.pm()) {} - - void reset() { - pmanager & pm = m_set.m(); - unsigned sz = m_set.size(); - for (unsigned i = 0; i < sz; ++i) { - m_in_set[pm.id(m_set.get(i))] = false; - } - m_set.reset(); - } - - void insert(poly * p) { - pmanager & pm = m_set.m(); - p = m_cache.mk_unique(p); - unsigned pid = pm.id(p); - if (m_in_set.get(pid, false)) - return; - m_in_set.setx(pid, true, false); - m_set.push_back(p); - } - - bool empty() const { return m_set.empty(); } - - // Return max variable in todo_set - var max_var() const { - pmanager & pm = m_set.m(); - var max = null_var; - unsigned sz = m_set.size(); - for (unsigned i = 0; i < sz; ++i) { - var x = pm.max_var(m_set.get(i)); - SASSERT(x != null_var); - if (max == null_var || x > max) - max = x; - } - return max; - } - - /** - \brief Remove the maximal polynomials from the set and store - them in max_polys. Return the maximal variable - */ - var extract_max_polys(polynomial_ref_vector & max_polys) { - max_polys.reset(); - var x = max_var(); - pmanager & pm = m_set.m(); - unsigned sz = m_set.size(); - unsigned j = 0; - for (unsigned i = 0; i < sz; ++i) { - poly * p = m_set.get(i); - var y = pm.max_var(p); - SASSERT(y <= x); - if (y == x) { - max_polys.push_back(p); - m_in_set[pm.id(p)] = false; - } - else { - m_set.set(j, p); - j++; - } - } - m_set.shrink(j); - return x; - } - }; - // temporary field for store todo set of polynomials todo_set m_todo; + + // Track polynomials already processed in current projection to avoid cycles + todo_set m_processed; // temporary fields for preprocessing core scoped_literal_vector m_core1; scoped_literal_vector m_core2; + + // Lower-stage polynomials encountered during normalization that need to be projected + polynomial_ref_vector m_lower_stage_polys; + + // Store last levelwise input for debugging unsound lemmas + polynomial_ref_vector m_last_lws_input_polys; // temporary fields for storing the result scoped_literal_vector * m_result = nullptr; @@ -136,9 +77,8 @@ namespace nlsat { evaluator & m_evaluator; imp(solver & s, assignment const & x2v, polynomial::cache & u, atom_vector const & atoms, atom_vector const & x2eq, - evaluator & ev, bool is_sample): + evaluator & ev, bool canonicalize): m_solver(s), - m_assignment(x2v), m_atoms(atoms), m_x2eq(x2eq), m_am(x2v.am()), @@ -150,10 +90,12 @@ namespace nlsat { m_factors(m_pm), m_factors_save(m_pm), m_roots_tmp(m_am), - m_cell_sample(is_sample), - m_todo(u), + m_todo(u, canonicalize), + m_processed(u, canonicalize), m_core1(s), m_core2(s), + m_lower_stage_polys(m_pm), + m_last_lws_input_polys(m_pm), m_evaluator(ev) { m_simplify_cores = false; m_full_dimensional = false; @@ -162,25 +104,8 @@ namespace nlsat { m_add_zero_disc = true; } - std::ostream& display(std::ostream & out, polynomial_ref const & p) const { - m_pm.display(out, p, m_solver.display_proc()); - return out; - } - - std::ostream& display(std::ostream & out, polynomial_ref_vector const & ps, char const * delim = "\n") const { - for (unsigned i = 0; i < ps.size(); ++i) { - if (i > 0) - out << delim; - m_pm.display(out, ps.get(i), m_solver.display_proc()); - } - return out; - } - - std::ostream& display(std::ostream & out, literal l) const { return m_solver.display(out, l); } - std::ostream& display_var(std::ostream & out, var x) const { return m_solver.display(out, x); } - std::ostream& display(std::ostream & out, unsigned sz, literal const * ls) const { return m_solver.display(out, sz, ls); } - std::ostream& display(std::ostream & out, literal_vector const & ls) const { return display(out, ls.size(), ls.data()); } - std::ostream& display(std::ostream & out, scoped_literal_vector const & ls) const { return display(out, ls.size(), ls.data()); } + // display helpers moved to free functions in nlsat_common.h + /** \brief Add literal to the result vector. @@ -208,28 +133,7 @@ namespace nlsat { SASSERT(check_already_added()); } - - /** - \brief evaluate the given polynomial in the current interpretation. - max_var(p) must be assigned in the current interpretation. - */ - ::sign sign(polynomial_ref const & p) { - SASSERT(max_var(p) == null_var || m_assignment.is_assigned(max_var(p))); - auto s = m_am.eval_sign_at(p, m_assignment); - TRACE(nlsat_explain, tout << "p: " << p << " var: " << max_var(p) << " sign: " << s << "\n";); - return s; - } - /** - \brief Wrapper for factorization - */ - void factor(polynomial_ref & p, polynomial_ref_vector & fs) { - // TODO: add params, caching - TRACE(nlsat_factor, tout << "factor\n" << p << "\n";); - fs.reset(); - m_cache.factor(p.get(), fs); - } - /** \brief Wrapper for psc chain computation */ @@ -306,14 +210,14 @@ namespace nlsat { void find_cell_roots(polynomial_ref_vector & ps, var y, cell_root_info & info) { info.reset(); - SASSERT(m_assignment.is_assigned(y)); + SASSERT(m_solver.sample().is_assigned(y)); bool lower_inf = true; bool upper_inf = true; scoped_anum_vector & roots = m_roots_tmp; scoped_anum lower(m_am); scoped_anum upper(m_am); - anum const & y_val = m_assignment.value(y); - TRACE(nlsat_explain, tout << "adding literals for "; display_var(tout, y); tout << " -> "; + anum const & y_val = m_solver.sample().value(y); + TRACE(nlsat_explain, tout << "adding literals for "; m_solver.display_var(tout, y); tout << " -> "; m_am.display_decimal(tout, y_val); tout << "\n";); polynomial_ref p(m_pm); unsigned sz = ps.size(); @@ -322,12 +226,12 @@ namespace nlsat { if (max_var(p) != y) continue; roots.reset(); - // Variable y is assigned in m_assignment. We must temporarily unassign it. + // Variable y is assigned in m_solver.sample(). We must temporarily unassign it. // Otherwise, the isolate_roots procedure will assume p is a constant polynomial. - m_am.isolate_roots(p, undef_var_assignment(m_assignment, y), roots); + m_am.isolate_roots(p, undef_var_assignment(m_solver.sample(), y), roots); unsigned num_roots = roots.size(); TRACE(nlsat_explain, - tout << "isolated roots for "; display_var(tout, y); + tout << "isolated roots for "; m_solver.display_var(tout, y); tout << " with polynomial: " << p << "\n"; for (unsigned ri = 0; ri < num_roots; ++ri) { m_am.display_decimal(tout << " root[" << (ri+1) << "] = ", roots[ri]); @@ -386,6 +290,11 @@ namespace nlsat { } } + ::sign sign(polynomial_ref const & p) { + return ::nlsat::sign(p, m_solver.sample(), m_am); + } + + void add_zero_assumption(polynomial_ref& p) { // Build a square-free representative of p so that we can speak about // a specific root that coincides with the current assignment. @@ -400,7 +309,7 @@ namespace nlsat { SASSERT(maxx != null_var); if (maxx == null_var) return; - SASSERT(m_assignment.is_assigned(maxx)); + SASSERT(m_solver.sample().is_assigned(maxx)); polynomial_ref_vector singleton(m_pm); singleton.push_back(q); @@ -411,7 +320,7 @@ namespace nlsat { TRACE(nlsat_explain, tout << "adding zero-assumption root literal for "; - display_var(tout, maxx); tout << " using root[" << info.m_eq_idx << "] of " << q << "\n";); + m_solver.display_var(tout, maxx); tout << " using root[" << info.m_eq_idx << "] of " << q << "\n";); add_root_literal(atom::ROOT_EQ, maxx, info.m_eq_idx, info.m_eq); } @@ -456,23 +365,23 @@ namespace nlsat { TRACE(nlsat_explain, tout << "lc: " << lc << " reduct: " << reduct << "\n";); insert_fresh_factors_in_todo(lc); if (!is_zero(lc) && sign(lc)) { - insert_fresh_factors_in_todo(lc); TRACE(nlsat_explain, tout << "lc does no vaninsh\n";); return; } - add_zero_assumption(lc); if (k == 0) { // all coefficients of p vanished in the current interpretation, // and were added as assumptions. p = m_pm.mk_zero(); TRACE(nlsat_explain, tout << "all coefficients of p vanished\n";); if (m_add_all_coeffs) { + add_zero_assumption(lc); return; } TRACE(nlsat_explain, tout << "falling back to add-all-coeffs projection\n";); m_add_all_coeffs = true; throw add_all_coeffs_restart(); } + add_zero_assumption(lc); k--; p = reduct; } @@ -524,6 +433,7 @@ namespace nlsat { Remark: root atoms are not normalized */ literal normalize(literal l, var max) { + TRACE(nlsat_explain, display(tout << "l:", m_solver, l) << '\n';); bool_var b = l.var(); if (b == true_bool_var) return l; @@ -546,6 +456,8 @@ namespace nlsat { SASSERT(max_var(p) != null_var); SASSERT(max_var(p) < max); // factor p is a lower stage polynomial, so we should add assumption to justify p being eliminated + // Also collect it for projection to ensure proper cell construction + m_lower_stage_polys.push_back(p); if (s == 0) add_assumption(atom::EQ, p); // add assumption p = 0 else if (a->is_even(i)) @@ -613,6 +525,7 @@ namespace nlsat { stages) from (arithmetic) literals, */ void normalize(scoped_literal_vector & C, var max) { + TRACE(nlsat_explain, display(tout << "C:\n", m_solver, C) << '\n';); unsigned sz = C.size(); unsigned j = 0; for (unsigned i = 0; i < sz; ++i) { @@ -699,6 +612,10 @@ namespace nlsat { \brief Add factors of p to todo */ void insert_fresh_factors_in_todo(polynomial_ref & p) { + // Skip if already processed in this projection (prevents cycles) + if (m_processed.contains(p)) + return; + if (is_const(p)) return; elim_vanishing(p); @@ -706,14 +623,14 @@ namespace nlsat { return; if (m_factor) { restore_factors _restore(m_factors, m_factors_save); - factor(p, m_factors); - TRACE(nlsat_explain, display(tout << "adding factors of\n", p); tout << "\n" << m_factors << "\n";); + factor(p, m_cache, m_factors); + TRACE(nlsat_explain, display(tout << "adding factors of\n", m_solver, p); tout << "\n" << m_factors << "\n";); polynomial_ref f(m_pm); for (unsigned i = 0; i < m_factors.size(); ++i) { f = m_factors.get(i); elim_vanishing(f); - if (!is_const(f)) { - TRACE(nlsat_explain, tout << "adding factor:\n"; display(tout, f); tout << "\n";); + if (!is_const(f) && !m_processed.contains(f)) { + TRACE(nlsat_explain, tout << "adding factor:\n"; display(tout, m_solver, f); tout << "\n";); m_todo.insert(f); } } @@ -736,10 +653,10 @@ namespace nlsat { unsigned k_deg = m_pm.degree(p, x); if (k_deg == 0) continue; // p depends on x - TRACE(nlsat_explain, tout << "processing poly of degree " << k_deg << " w.r.t x" << x << ": "; display(tout, p) << "\n";); + TRACE(nlsat_explain, tout << "processing poly of degree " << k_deg << " w.r.t x" << x << ": "; m_solver.display(tout, p) << "\n";); for (int j_coeff_deg = k_deg; j_coeff_deg >= 0; j_coeff_deg--) { coeff = m_pm.coeff(p, x, j_coeff_deg); - TRACE(nlsat_explain, tout << " coeff deg " << j_coeff_deg << ": "; display(tout, coeff) << "\n";); + TRACE(nlsat_explain, tout << " coeff deg " << j_coeff_deg << ": "; display(tout, m_solver, coeff) << "\n";); insert_fresh_factors_in_todo(coeff); if (!m_add_all_coeffs) break; @@ -766,7 +683,7 @@ namespace nlsat { // this function also explains the value 0, if met bool coeffs_are_zeroes(polynomial_ref &s) { restore_factors _restore(m_factors, m_factors_save); - factor(s, m_factors); + factor(s, m_cache, m_factors); unsigned num_factors = m_factors.size(); m_zero_fs.reset(); m_is_even.reset(); @@ -774,7 +691,7 @@ namespace nlsat { bool have_zero = false; for (unsigned i = 0; i < num_factors; ++i) { f = m_factors.get(i); - if (coeffs_are_zeroes_in_factor(f)) { + if (coeffs_are_zeroes_on_sample(f, m_pm, sample(), m_am)) { have_zero = true; break; } @@ -799,7 +716,7 @@ namespace nlsat { auto c = polynomial_ref(this->m_pm); for (unsigned j = 0; j <= n; ++j) { c = m_pm.coeff(s, x, j); - if (sign(c) != 0) + if (nlsat::sign(c, sample(), m_am) != 0) return false; } return true; @@ -814,7 +731,7 @@ namespace nlsat { psc_chain(p, q, x, S); unsigned sz = S.size(); - TRACE(nlsat_explain, tout << "computing psc of\n"; display(tout, p); tout << "\n"; display(tout, q); tout << "\n"; + TRACE(nlsat_explain, tout << "computing psc of\n"; display(tout, m_solver, p); tout << "\n"; display(tout, m_solver, q); tout << "\n"; for (unsigned i = 0; i < sz; ++i) { s = S.get(i); tout << "psc: " << s << "\n"; @@ -822,10 +739,16 @@ namespace nlsat { for (unsigned i = 0; i < sz; ++i) { s = S.get(i); - TRACE(nlsat_explain, display(tout << "processing psc(" << i << ")\n", s) << "\n";); + TRACE(nlsat_explain, display(tout << "processing psc(" << i << ")\n", m_solver, s) << "\n";); if (is_zero(s)) { - TRACE(nlsat_explain, tout << "skipping psc is the zero polynomial\n";); - continue; + // PSC is identically zero - polynomials share a common factor. + // This can cause unsound lemmas. Fall back to add-all-coeffs projection. + TRACE(nlsat_explain, tout << "psc is zero polynomial - polynomials share common factor\n";); + if (m_add_all_coeffs) + continue; + TRACE(nlsat_explain, tout << "falling back to add-all-coeffs projection\n";); + m_add_all_coeffs = true; + throw add_all_coeffs_restart(); } if (is_const(s)) { TRACE(nlsat_explain, tout << "done, psc is a constant\n";); @@ -836,11 +759,11 @@ namespace nlsat { } TRACE(nlsat_explain, tout << "adding v-psc of\n"; - display(tout, p); + display(tout, m_solver, p); tout << "\n"; - display(tout, q); + display(tout, m_solver, q); tout << "\n---->\n"; - display(tout, s); + display(tout, m_solver, s); tout << "\n";); // s did not vanish completely, but its leading coefficient may have vanished insert_fresh_factors_in_todo(s); @@ -900,14 +823,14 @@ namespace nlsat { void add_root_literal(atom::kind k, var y, unsigned i, poly * p) { polynomial_ref pr(p, m_pm); - TRACE(nlsat_explain, - display(tout << "x" << y << " " << k << "[" << i << "](", pr); tout << ")\n";); + TRACE(nlsat_explain, + display(tout << "x" << y << " " << k << "[" << i << "](", m_solver, pr); tout << ")\n";); if (!mk_linear_root(k, y, i, p) && !mk_quadratic_root(k, y, i, p)) { bool_var b = m_solver.mk_root_atom(k, y, i, p); literal l(b, true); - TRACE(nlsat_explain, tout << "adding literal\n"; display(tout, l); tout << "\n";); + TRACE(nlsat_explain, tout << "adding literal\n"; display(tout, m_solver, l); tout << "\n";); add_literal(l); } } @@ -968,7 +891,7 @@ namespace nlsat { return false; } - SASSERT(m_assignment.is_assigned(y)); + SASSERT(sample().is_assigned(y)); polynomial_ref A(m_pm), B(m_pm), C(m_pm), q(m_pm), p_diff(m_pm), yy(m_pm); A = m_pm.coeff(p, y, 2); B = m_pm.coeff(p, y, 1); @@ -1000,6 +923,8 @@ namespace nlsat { if (!is_const(p)) { TRACE(nlsat_explain, tout << p << "\n";); add_assumption(s == 0 ? atom::EQ : (s < 0 ? atom::LT : atom::GT), p); + // Also add p to the projection to ensure proper cell construction + insert_fresh_factors_in_todo(p); } return s; } @@ -1013,7 +938,7 @@ namespace nlsat { */ void mk_linear_root(atom::kind k, var y, unsigned i, poly * p, bool mk_neg) { - TRACE(nlsat_explain, display_var(tout, y); m_pm.display(tout << ": ", p, m_solver.display_proc()); tout << "\n"); + TRACE(nlsat_explain, m_solver.display_var(tout, y); m_pm.display(tout << ": ", p, m_solver.display_proc()); tout << "\n"); polynomial_ref p_prime(m_pm); p_prime = p; bool lsign = false; @@ -1104,33 +1029,37 @@ namespace nlsat { \brief Apply model-based projection operation defined in our paper. */ - void project_original(polynomial_ref_vector & ps, var max_x) { - if (ps.empty()) - return; - m_todo.reset(); - for (poly* p : ps) { - m_todo.insert(p); - } - var x = m_todo.extract_max_polys(ps); - // Remark: after vanishing coefficients are eliminated, ps may not contain max_x anymore - if (x < max_x) - add_cell_lits(ps, x); - while (true) { - if (all_univ(ps, x) && m_todo.empty()) { - m_todo.reset(); - break; + bool levelwise_single_cell(polynomial_ref_vector & ps, var max_x, polynomial::cache & cache) { + // Store polynomials for debugging unsound lemmas + m_last_lws_input_polys.reset(); + for (unsigned i = 0; i < ps.size(); i++) + m_last_lws_input_polys.push_back(ps.get(i)); + + levelwise lws(m_solver, ps, max_x, sample(), m_pm, m_am, cache); + auto cell = lws.single_cell(); + if (lws.failed()) + return false; + + TRACE(lws, for (unsigned i = 0; i < cell.size(); i++) + display(tout << "I[" << i << "]:", m_solver, cell[i]) << "\n";); + // Enumerate all intervals in the computed cell and add literals for each non-trivial interval. + // Non-trivial = section, or sector with at least one finite bound (ignore (-oo,+oo)). + for (auto const & I : cell) { + if (I.is_section()) { + SASSERT(I.l && I.l_index); + add_root_literal(atom::ROOT_EQ, max_var(I.l.get()), I.l_index, I.l.get()); + continue; } - TRACE(nlsat_explain, tout << "project loop, processing var "; display_var(tout, x); - tout << "\npolynomials\n"; - display(tout, ps); tout << "\n";); - add_lcs(ps, x); - psc_discriminant(ps, x); - psc_resultant(ps, x); - if (m_todo.empty()) - break; - x = m_todo.extract_max_polys(ps); - add_cell_lits(ps, x); + if (I.l_inf() && I.u_inf()) + continue; // skip whole-line sector + if (!I.l_inf()) + add_root_literal(m_full_dimensional ? atom::ROOT_GE : + atom::ROOT_GT, max_var(I.l.get()), I.l_index, I.l.get()); + if (!I.u_inf()) + add_root_literal(m_full_dimensional ? atom::ROOT_LE : + atom::ROOT_LT, max_var(I.u.get()), I.u_index, I.u.get()); } + return true; } /** @@ -1140,12 +1069,17 @@ namespace nlsat { * "Solving Satisfiability of Polynomial Formulas By Sample - Cell Projection" * https://arxiv.org/abs/2003.00409 */ + void collect_to_processed(polynomial_ref_vector & ps) { + for (unsigned i = 0; i < ps.size(); ++i) + m_processed.insert(ps.get(i)); + } + void project_cdcac(polynomial_ref_vector & ps, var max_x) { bool first = true; if (ps.empty()) return; - m_todo.reset(); + m_processed.reset(); for (unsigned i = 0; i < ps.size(); ++i) { polynomial_ref p(m_pm); p = ps.get(i); @@ -1155,10 +1089,22 @@ namespace nlsat { ps.reset(); for (auto p: m_todo.m_set) ps.push_back(p); + + bool use_all_coeffs = false; + + if (m_solver.apply_levelwise()) { + bool levelwise_ok = levelwise_single_cell(ps, max_x, m_cache); + m_solver.record_levelwise_result(levelwise_ok); + if (levelwise_ok) + return; + // Levelwise failed, use add_all_coeffs mode for fallback projection + use_all_coeffs = true; + } + // Set m_add_all_coeffs for the rest of the projection, restore on function exit + flet _set_all_coeffs(m_add_all_coeffs, use_all_coeffs || m_add_all_coeffs); var x = m_todo.extract_max_polys(ps); - // Remark: after vanishing coefficients are eliminated, ps may not contain max_x anymore - + collect_to_processed(ps); polynomial_ref_vector samples(m_pm); if (x < max_x) cac_add_cell_lits(ps, x, samples); @@ -1168,8 +1114,8 @@ namespace nlsat { m_todo.reset(); break; } - TRACE(nlsat_explain, tout << "project loop, processing var "; display_var(tout, x); tout << "\npolynomials\n"; - display(tout, ps); tout << "\n";); + TRACE(nlsat_explain, tout << "project loop, processing var "; m_solver.display_var(tout, x); tout << "\npolynomials\n"; + display(tout, m_solver, ps); tout << "\n";); if (first) { // The first run is special because x is not constrained by the sample, we cannot surround it by the root functions. // we make the polynomials in ps delinable add_lcs(ps, x); @@ -1186,17 +1132,14 @@ namespace nlsat { if (m_todo.empty()) break; x = m_todo.extract_max_polys(ps); + collect_to_processed(ps); cac_add_cell_lits(ps, x, samples); } + } void project(polynomial_ref_vector & ps, var max_x) { - if (m_cell_sample) { - project_cdcac(ps, max_x); - } - else { - project_original(ps, max_x); - } + project_cdcac(ps, max_x); } bool check_already_added() const { @@ -1291,7 +1234,7 @@ namespace nlsat { new_lit = l; return; } - TRACE(nlsat_simplify_core, display(tout << "trying to simplify literal\n", l) << "\nusing equation\n"; + TRACE(nlsat_simplify_core, display(tout << "trying to simplify literal\n", m_solver, l) << "\nusing equation\n"; m_pm.display(tout, info.m_eq, m_solver.display_proc()); tout << "\n";); int atom_sign = 1; bool modified_lit = false; @@ -1374,14 +1317,14 @@ namespace nlsat { new_lit = m_solver.mk_ineq_literal(new_k, new_factors.size(), new_factors.data(), new_factors_even.data()); if (l.sign()) new_lit.neg(); - TRACE(nlsat_simplify_core, tout << "simplified literal:\n"; display(tout, new_lit) << " " << m_solver.value(new_lit) << "\n";); + TRACE(nlsat_simplify_core, tout << "simplified literal:\n"; display(tout, m_solver, new_lit) << " " << m_solver.value(new_lit) << "\n";); if (max_var(new_lit) < max) { if (m_solver.value(new_lit) == l_true) { TRACE(nlsat_simplify_bug, tout << "literal normalized away because it is already true after rewriting:\n"; - display(tout << " original: ", l) << "\n"; - display(tout << " rewritten: ", new_lit) << "\n"; + m_solver.display(tout << " original: ", l) << "\n"; + m_solver.display(tout << " rewritten: ", new_lit) << "\n"; if (info.m_eq) { polynomial_ref eq_ref(const_cast(info.m_eq), m_pm); m_pm.display(tout << " equation used: ", eq_ref, m_solver.display_proc()); @@ -1393,8 +1336,8 @@ namespace nlsat { literal assumption = new_lit; TRACE(nlsat_simplify_bug, tout << "literal replaced by lower-stage assumption due to rewriting:\n"; - display(tout << " original: ", l) << "\n"; - display(tout << " assumption: ", assumption) << "\n"; + m_solver.display(tout << " original: ", l) << "\n"; + m_solver.display(tout << " assumption: ", assumption) << "\n"; if (info.m_eq) { polynomial_ref eq_ref(const_cast(info.m_eq), m_pm); m_pm.display(tout << " equation used: ", eq_ref, m_solver.display_proc()); @@ -1408,12 +1351,12 @@ namespace nlsat { literal before = new_lit; (void)before; new_lit = normalize(new_lit, max); - TRACE(nlsat_simplify_core, tout << "simplified literal after normalization:\n"; display(tout, new_lit); tout << " " << m_solver.value(new_lit) << "\n";); + TRACE(nlsat_simplify_core, tout << "simplified literal after normalization:\n"; m_solver.display(tout, new_lit); tout << " " << m_solver.value(new_lit) << "\n";); if (new_lit == true_literal || new_lit == false_literal) { TRACE(nlsat_simplify_bug, tout << "normalize() turned rewritten literal into constant value:\n"; - display(tout << " original: ", l) << "\n"; - display(tout << " rewritten: ", before) << "\n"; + m_solver.display(tout << " original: ", l) << "\n"; + m_solver.display(tout << " rewritten: ", before) << "\n"; tout << " result: " << (new_lit == true_literal ? "true" : "false") << "\n"; if (info.m_eq) { polynomial_ref eq_ref(const_cast(info.m_eq), m_pm); @@ -1577,7 +1520,7 @@ namespace nlsat { poly * eq_p = eq->p(0); VERIFY(simplify(C, eq_p, max)); // add equation as an assumption - TRACE(nlsat_simpilfy_core, display(tout << "adding equality as assumption ", literal(eq->bvar(), true)); tout << "\n";); + TRACE(nlsat_simpilfy_core, display(tout << "adding equality as assumption ", m_solver, literal(eq->bvar(), true)); tout << "\n";); add_literal(literal(eq->bvar(), true)); } } @@ -1586,33 +1529,52 @@ namespace nlsat { \brief Main procedure. The explain the given unsat core, and store the result in m_result */ void main(unsigned num, literal const * ls) { - if (num == 0) + if (num == 0) { + TRACE(nlsat_explain, tout << "num:" << num << "\n";); return; + } collect_polys(num, ls, m_ps); + + // Add lower-stage polynomials collected during normalization + // These polynomials had assumptions added but need to be projected as well + for (unsigned i = 0; i < m_lower_stage_polys.size(); i++) { + m_ps.push_back(m_lower_stage_polys.get(i)); + } + var max_x = max_var(m_ps); - TRACE(nlsat_explain, tout << "polynomials in the conflict:\n"; display(tout, m_ps); tout << "\n";); + TRACE(nlsat_explain, tout << "polynomials in the conflict:\n"; display(tout, m_solver, m_ps); tout << "\n";); + + // Note: levelwise is now called in process2() before normalize() + // to ensure it receives the original polynomials + elim_vanishing(m_ps); - TRACE(nlsat_explain, tout << "after elim_vanishing m_ps:\n"; display(tout, m_ps); tout << "\n";); + TRACE(nlsat_explain, tout << "after elim_vanishing m_ps:\n"; display(tout, m_solver, m_ps); tout << "\n";); project(m_ps, max_x); - TRACE(nlsat_explain, tout << "after projection\n"; display(tout, m_ps); tout << "\n";); + TRACE(nlsat_explain, tout << "after projection\n"; display(tout, m_solver, m_ps); tout << "\n";); } void process2(unsigned num, literal const * ls) { + // Reset lower-stage polynomials collection + m_lower_stage_polys.reset(); + + // Try levelwise with ORIGINAL polynomials BEFORE any normalization + if (m_simplify_cores) { + TRACE(nlsat_explain, tout << "m_simplify_cores is true\n";); m_core2.reset(); m_core2.append(num, ls); var max = max_var(num, ls); SASSERT(max != null_var); - TRACE(nlsat_explain, display(tout << "core before normalization\n", m_core2) << "\n";); + TRACE(nlsat_explain, display(tout << "core before normalization\n", m_solver, m_core2) << "\n";); normalize(m_core2, max); - TRACE(nlsat_explain, display(tout << "core after normalization\n", m_core2) << "\n";); + TRACE(nlsat_explain, display(tout << "core after normalization\n", m_solver, m_core2) << "\n";); simplify(m_core2, max); - TRACE(nlsat_explain, display(tout << "core after simplify\n", m_core2) << "\n";); + TRACE(nlsat_explain, display(tout << "core after simplify\n", m_solver, m_core2) << "\n";); main(m_core2.size(), m_core2.data()); m_core2.reset(); } else { - TRACE(nlsat_explain, display(tout << "core befor normalization\n", m_core2) << "\n";); + TRACE(nlsat_explain, display(tout << "core befor normalization\n", m_solver, m_core2) << "\n";); main(num, ls); } } @@ -1674,39 +1636,42 @@ namespace nlsat { todo.reset(); core.reset(); todo.append(num, ls); while (true) { - TRACE(nlsat_minimize, tout << "core minimization:\n"; display(tout, todo); tout << "\nCORE:\n"; display(tout, core) << "\n";); + TRACE(nlsat_minimize, tout << "core minimization:\n"; display(tout, m_solver, todo); tout << "\nCORE:\n"; display(tout, m_solver, core) << "\n";); if (!minimize_core(todo, core)) break; std::reverse(todo.begin(), todo.end()); - TRACE(nlsat_minimize, tout << "core minimization:\n"; display(tout, todo); tout << "\nCORE:\n"; display(tout, core) << "\n";); + TRACE(nlsat_minimize, tout << "core minimization:\n"; display(tout, m_solver, todo); tout << "\nCORE:\n"; display(tout, m_solver, core) << "\n";); if (!minimize_core(todo, core)) break; } - TRACE(nlsat_minimize, tout << "core:\n"; display(tout, core);); + TRACE(nlsat_minimize, tout << "core:\n"; display(tout, m_solver, core);); r.append(core.size(), core.data()); } void process(unsigned num, literal const * ls) { if (m_minimize_cores && num > 1) { + TRACE(nlsat_explain, tout << "m_minimize_cores:" << m_minimize_cores << ", num:" << num;); m_core1.reset(); minimize(num, ls, m_core1); process2(m_core1.size(), m_core1.data()); m_core1.reset(); } else { + TRACE(nlsat_explain, tout << "directly to process2\n";); process2(num, ls); } } - void operator()(unsigned num, literal const * ls, scoped_literal_vector & result) { + void compute_conflict_explanation(unsigned num, literal const * ls, scoped_literal_vector & result) { SASSERT(check_already_added()); SASSERT(num > 0); TRACE(nlsat_explain, - tout << "[explain] set of literals is infeasible in the current interpretation\n"; - display(tout, num, ls) << "\n"; - m_solver.display_assignment(tout); + tout << "the infeasible clause:\n"; + display(tout, m_solver, num, ls) << "\n"; + + m_solver.display_assignment(tout << "assignment:\n"); ); - if (max_var(num, ls) == 0 && !m_assignment.is_assigned(0)) { + if (max_var(num, ls) == 0 && !m_solver.sample().is_assigned(0)) { TRACE(nlsat_explain, tout << "all literals use unassigned max var; returning justification\n";); result.reset(); return; @@ -1718,7 +1683,7 @@ namespace nlsat { process(num, ls); reset_already_added(); m_result = nullptr; - TRACE(nlsat_explain, display(tout << "[explain] result\n", result) << "\n";); + TRACE(nlsat_explain, display(tout << "[explain] result\n", m_solver, result) << "\n";); CASSERT("nlsat", check_already_added()); break; } @@ -1848,12 +1813,12 @@ namespace nlsat { collect_polys(lits.size(), lits.data(), m_ps); unbounded = true; scoped_anum x_val(m_am); - x_val = m_assignment.value(x); + x_val = sample().value(x); for (unsigned i = 0; i < m_ps.size(); ++i) { p = m_ps.get(i); scoped_anum_vector & roots = m_roots_tmp; roots.reset(); - m_am.isolate_roots(p, undef_var_assignment(m_assignment, x), roots); + m_am.isolate_roots(p, undef_var_assignment(sample(), x), roots); for (unsigned j = 0; j < roots.size(); ++j) { int s = m_am.compare(x_val, roots[j]); if (s <= 0 && (unbounded || m_am.compare(roots[j], val) <= 0)) { @@ -1867,8 +1832,8 @@ namespace nlsat { }; explain::explain(solver & s, assignment const & x2v, polynomial::cache & u, - atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev, bool use_cell_sample) { - m_imp = alloc(imp, s, x2v, u, atoms, x2eq, ev, use_cell_sample); + atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev, bool canonicalize) { + m_imp = alloc(imp, s, x2v, u, atoms, x2eq, ev, canonicalize); } explain::~explain() { @@ -1904,8 +1869,8 @@ namespace nlsat { m_imp->m_add_zero_disc = f; } - void explain::main_operator(unsigned n, literal const * ls, scoped_literal_vector & result) { - (*m_imp)(n, ls, result); + void explain::compute_conflict_explanation(unsigned n, literal const * ls, scoped_literal_vector & result) { + m_imp->compute_conflict_explanation(n, ls, result); } void explain::project(var x, unsigned n, literal const * ls, scoped_literal_vector & result) { @@ -1916,6 +1881,16 @@ namespace nlsat { m_imp->maximize(x, n, ls, val, unbounded); } + void explain::display_last_lws_input(std::ostream& out) { + out << "=== POLYNOMIALS PASSED TO LEVELWISE ===\n"; + for (unsigned i = 0; i < m_imp->m_last_lws_input_polys.size(); i++) { + out << " p[" << i << "]: "; + m_imp->m_pm.display(out, m_imp->m_last_lws_input_polys.get(i)); + out << "\n"; + } + out << "=== END LEVELWISE INPUT (" << m_imp->m_last_lws_input_polys.size() << " polynomials) ===\n"; + } + void explain::test_root_literal(atom::kind k, var y, unsigned i, poly* p, scoped_literal_vector & result) { m_imp->test_root_literal(k, y, i, p, result); } @@ -1924,29 +1899,29 @@ namespace nlsat { #ifdef Z3DEBUG #include void pp(nlsat::explain::imp & ex, unsigned num, nlsat::literal const * ls) { - ex.display(std::cout, num, ls); + display(std::cout, ex.m_solver, num, ls); } void pp(nlsat::explain::imp & ex, nlsat::scoped_literal_vector & ls) { - ex.display(std::cout, ls); + display(std::cout, ex.m_solver, ls); } void pp(nlsat::explain::imp & ex, polynomial_ref const & p) { - ex.display(std::cout, p); + display(std::cout, ex.m_solver, p); std::cout << std::endl; } void pp(nlsat::explain::imp & ex, polynomial::polynomial * p) { polynomial_ref _p(p, ex.m_pm); - ex.display(std::cout, _p); + display(std::cout, ex.m_solver, _p); std::cout << std::endl; } void pp(nlsat::explain::imp & ex, polynomial_ref_vector const & ps) { - ex.display(std::cout, ps); + display(std::cout, ex.m_solver, ps); } void pp_var(nlsat::explain::imp & ex, nlsat::var x) { - ex.display(std::cout, x); + display_var(std::cout, ex.m_solver, x); std::cout << std::endl; } void pp_lit(nlsat::explain::imp & ex, nlsat::literal l) { - ex.display(std::cout, l); + display(std::cout, ex.m_solver, l); std::cout << std::endl; } #endif diff --git a/src/nlsat/nlsat_explain.h b/src/nlsat/nlsat_explain.h index 6ca08e699..3ff4fb982 100644 --- a/src/nlsat/nlsat_explain.h +++ b/src/nlsat/nlsat_explain.h @@ -35,7 +35,7 @@ namespace nlsat { imp * m_imp; public: explain(solver & s, assignment const & x2v, polynomial::cache & u, - atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev, bool use_cell_sample_proj); + atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev, bool canonicalize); ~explain(); @@ -64,7 +64,7 @@ namespace nlsat { - s_1, ..., s_m do not contain variable x. - s_1, ..., s_m are false in the current interpretation */ - void main_operator(unsigned n, literal const * ls, scoped_literal_vector & result); + void compute_conflict_explanation(unsigned n, literal const * ls, scoped_literal_vector & result); /** @@ -103,6 +103,11 @@ namespace nlsat { */ void maximize(var x, unsigned n, literal const * ls, scoped_anum& val, bool& unbounded); + /** + Print the polynomials that were passed to levelwise in the last call (for debugging). + */ + void display_last_lws_input(std::ostream& out); + /** Unit test routine. */ diff --git a/src/nlsat/nlsat_interval_set.h b/src/nlsat/nlsat_interval_set.h index 297318912..263ab8a10 100644 --- a/src/nlsat/nlsat_interval_set.h +++ b/src/nlsat/nlsat_interval_set.h @@ -35,6 +35,7 @@ namespace nlsat { interval_set_manager(anum_manager & m, small_object_allocator & a); void set_seed(unsigned s) { m_rand.set_seed(s); } + unsigned get_seed() const { return m_rand.get_seed(); } /** \brief Return the empty set. diff --git a/src/nlsat/nlsat_params.pyg b/src/nlsat/nlsat_params.pyg index 2403f94b2..1ab305984 100644 --- a/src/nlsat/nlsat_params.pyg +++ b/src/nlsat/nlsat_params.pyg @@ -5,7 +5,6 @@ def_module_params('nlsat', params=(max_memory_param(), ('simple_check', BOOL, False, "precheck polynomials using variables sign"), ('variable_ordering_strategy', UINT, 0, "Variable Ordering Strategy, 0 for none, 1 for BROWN, 2 for TRIANGULAR, 3 for ONLYPOLY"), - ('cell_sample', BOOL, True, "cell sample projection"), ('lazy', UINT, 0, "how lazy the solver is."), ('reorder', BOOL, True, "reorder variables."), ('log_lemmas', BOOL, False, "display lemmas as self-contained SMT formulas"), @@ -22,6 +21,8 @@ def_module_params('nlsat', ('factor', BOOL, True, "factor polynomials produced during conflict resolution."), ('add_all_coeffs', BOOL, False, "add all polynomial coefficients during projection."), ('zero_disc', BOOL, False, "add_zero_assumption to the vanishing discriminant."), - ('known_sat_assignment_file_name', STRING, "", "the file name of a known solution: used for debugging only") - - )) + ('known_sat_assignment_file_name', STRING, "", "the file name of a known solution: used for debugging only"), + ('lws', BOOL, True, "apply levelwise."), + ('lws_spt_threshold', UINT, 2, "minimum both-side polynomial count to apply spanning tree optimization; < 2 disables spanning tree"), + ('canonicalize', BOOL, True, "canonicalize polynomials.") + )) diff --git a/src/nlsat/nlsat_solver.cpp b/src/nlsat/nlsat_solver.cpp index a8497aae5..2e1f7df48 100644 --- a/src/nlsat/nlsat_solver.cpp +++ b/src/nlsat/nlsat_solver.cpp @@ -1,3 +1,4 @@ +// int tttt = 0; /*++ Copyright (c) 2012 Microsoft Corporation @@ -41,6 +42,7 @@ Revision History: #include "nlsat/nlsat_simplify.h" #include "nlsat/nlsat_simple_checker.h" #include "nlsat/nlsat_variable_ordering_strategy.h" +#include "nlsat_solver.h" #define NLSAT_EXTRA_VERBOSE @@ -52,7 +54,6 @@ Revision History: namespace nlsat { - typedef chashtable ineq_atom_table; typedef chashtable root_atom_table; @@ -227,9 +228,8 @@ namespace nlsat { unsigned m_max_conflicts; unsigned m_lemma_rlimit; unsigned m_lemma_count; - unsigned m_variable_ordering_strategy; + unsigned m_variable_ordering_strategy; bool m_set_0_more; - bool m_cell_sample; struct stats { unsigned m_simplifications; @@ -239,13 +239,18 @@ namespace nlsat { unsigned m_decisions; unsigned m_stages; unsigned m_irrational_assignments; // number of irrational witnesses + unsigned m_levelwise_calls; + unsigned m_levelwise_failures; + unsigned m_lws_initial_fail; void reset() { memset(this, 0, sizeof(*this)); } stats() { reset(); } }; // statistics stats m_stats; std::string m_debug_known_solution_file_name; - + bool m_apply_lws; + bool m_last_conflict_used_lws = false; // Track if last conflict explanation used levelwise + unsigned m_lws_spt_threshold = 3; imp(solver& s, ctx& c): m_ctx(c), m_solver(s), @@ -264,7 +269,7 @@ namespace nlsat { m_simplify(s, m_atoms, m_clauses, m_learned, m_pm), m_display_var(m_perm), m_display_assumption(nullptr), - m_explain(s, m_assignment, m_cache, m_atoms, m_var2eq, m_evaluator, nlsat_params(c.m_params).cell_sample()), + m_explain(s, m_assignment, m_cache, m_atoms, m_var2eq, m_evaluator, nlsat_params(c.m_params).canonicalize()), m_scope_lvl(0), m_lemma(s), m_lazy_clause(s), @@ -305,8 +310,9 @@ namespace nlsat { m_check_lemmas = p.check_lemmas(); m_variable_ordering_strategy = p.variable_ordering_strategy(); m_debug_known_solution_file_name = p.known_sat_assignment_file_name(); + m_apply_lws = p.lws(); + m_lws_spt_threshold = p.lws_spt_threshold(); // 0 disables spanning tree m_check_lemmas |= !(m_debug_known_solution_file_name.empty()); - m_cell_sample = p.cell_sample(); m_ism.set_seed(m_random_seed); m_explain.set_simplify_cores(m_simplify_cores); @@ -1016,7 +1022,8 @@ namespace nlsat { } // Helper: Setup checker solver and translate atoms/clauses - void setup_checker(imp& checker, scoped_bool_vars& tr, unsigned n, literal const* cls, assumption_set a) { + // Returns false if the lemma cannot be properly translated for checking + bool setup_checker(imp& checker, scoped_bool_vars& tr, unsigned n, literal const* cls, assumption_set a) { auto pconvert = [&](poly* p) { return convert(m_pm, p, checker.m_pm); }; @@ -1052,7 +1059,9 @@ namespace nlsat { } else { // root atom cannot be translated due to variable ordering - bv = checker.mk_bool_var(); + // Skip lemma check in this case + TRACE(nlsat, tout << "check_lemma: skipping due to untranslatable root atom\n";); + return false; } } else { @@ -1079,20 +1088,32 @@ namespace nlsat { literal nlit(tr[lit.var()], !lit.sign()); checker.mk_external_clause(1, &nlit, nullptr); } + return true; // Successfully set up the checker } // Helper: Display unsound lemma failure information void display_unsound_lemma(imp& checker, scoped_bool_vars& tr, unsigned n, literal const* cls) { verbose_stream() << "\n"; verbose_stream() << "========== UNSOUND LEMMA DETECTED ==========\n"; + verbose_stream() << "Levelwise used for this conflict: " << (m_last_conflict_used_lws ? "YES" : "NO") << "\n"; + + // Print polynomials passed to levelwise + if (m_last_conflict_used_lws) { + m_explain.display_last_lws_input(verbose_stream()); + } + verbose_stream() << "The following lemma is NOT implied by the original clauses:\n"; display(verbose_stream() << " Lemma: ", n, cls) << "\n\n"; verbose_stream() << "Reason: Found a satisfying assignment where:\n"; verbose_stream() << " - The original clauses are satisfied\n"; verbose_stream() << " - But ALL literals in the lemma are FALSE\n\n"; + // Display sample point (original solver's assignment) + verbose_stream() << "Variable values at SAMPLE point:\n"; + display_num_assignment(verbose_stream()); + // Display variable values used in the lemma - verbose_stream() << "Variable values in counterexample:\n"; + verbose_stream() << "\nVariable values in counterexample:\n"; auto lemma_vars = collect_vars_on_clause(n, cls); for (var x : lemma_vars) { if (checker.m_assignment.is_assigned(x)) { @@ -1120,11 +1141,15 @@ namespace nlsat { TRACE(nlsat, display(tout << "check lemma: ", n, cls) << "\n"; display(tout);); + // Save RNG state before check_lemma to ensure determinism + unsigned saved_random_seed = m_random_seed; + unsigned saved_ism_seed = m_ism.get_seed(); + try { // Create a separate reslimit for the checker with 10 second timeout reslimit checker_rlimit; cancel_eh eh(checker_rlimit); - scoped_timer timer(10000, &eh); + scoped_timer timer(1000, &eh); // one second ctx c(checker_rlimit, m_ctx.m_params, m_ctx.m_incremental); solver solver2(c); @@ -1133,14 +1158,28 @@ namespace nlsat { checker.m_log_lemmas = false; checker.m_dump_mathematica = false; checker.m_inline_vars = false; + checker.m_apply_lws = false; // Disable levelwise for checker to avoid recursive issues scoped_bool_vars tr(checker); - setup_checker(checker, tr, n, cls, a); + if (!setup_checker(checker, tr, n, cls, a)) { + // Restore RNG state + m_random_seed = saved_random_seed; + m_ism.set_seed(saved_ism_seed); + return; // Lemma contains untranslatable atoms, skip check + } lbool r = checker.check(); - if (r == l_undef) // Checker timed out - skip this lemma check - return; + if (r == l_undef) { + // Restore RNG state + m_random_seed = saved_random_seed; + m_ism.set_seed(saved_ism_seed); + return; // Checker timed out - skip this lemma check + } if (r == l_true) { + // Before reporting unsound, dump the lemma to see what we're checking + verbose_stream() << "Dumping lemma that internal checker thinks is not a tautology:\n"; + verbose_stream() << "Checker levelwise calls: " << checker.m_stats.m_levelwise_calls << "\n"; + log_lemma(verbose_stream(), n, cls, true, "internal-check-fail"); display_unsound_lemma(checker, tr, n, cls); exit(1); } @@ -1148,6 +1187,10 @@ namespace nlsat { catch (...) { // Ignore exceptions from the checker - just skip this lemma check } + + // Restore RNG state after check_lemma + m_random_seed = saved_random_seed; + m_ism.set_seed(saved_ism_seed); } void log_lemma(std::ostream& out, clause const& cls, std::string annotation) { @@ -1165,7 +1208,7 @@ namespace nlsat { used_bools[b] = true; vars.reset(); this->vars(lit, vars); - for (var v : vars) + for (var v : vars) used_vars[v] = true; } display(out << "(echo \"#" << m_lemma_count++ << ":" << annotation << ":", n, cls) << "\")\n"; @@ -1182,7 +1225,6 @@ namespace nlsat { for (unsigned i = 0; i < n; ++i) display_smt2(out << "(assert ", ~cls[i]) << ")\n"; out << "(check-sat)\n(reset)\n"; - } clause * mk_clause_core(unsigned num_lits, literal const * lits, bool learned, _assumption_set a) { @@ -2120,9 +2162,8 @@ namespace nlsat { collect(assumptions, m_learned); del_clauses(m_valids); - if (m_check_lemmas) - for (clause* c : m_learned) - check_lemma(c->size(), c->data(), nullptr); + // Note: Don't check learned clauses here - they are the result of resolution + // and may not be tautologies. Conflict lemmas are checked in resolve_lazy_justification. assumptions.reset(); assumptions.append(result); @@ -2338,8 +2379,8 @@ namespace nlsat { m_lazy_clause.reset(); - m_explain.main_operator(jst.num_lits(), jst.lits(), m_lazy_clause); - for (unsigned i = 0; i < sz; ++i) + m_explain.compute_conflict_explanation(jst.num_lits(), jst.lits(), m_lazy_clause); + for (unsigned i = 0; i < sz; i++) m_lazy_clause.push_back(~jst.lit(i)); // lazy clause is a valid clause @@ -2359,6 +2400,8 @@ namespace nlsat { } if (m_check_lemmas) { + TRACE(nlsat, tout << "Checking lazy clause with " << m_lazy_clause.size() << " literals:\n"; + display(tout, m_lazy_clause.size(), m_lazy_clause.data()) << "\n";); check_lemma(m_lazy_clause.size(), m_lazy_clause.data(), nullptr); m_valids.push_back(mk_clause_core(m_lazy_clause.size(), m_lazy_clause.data(), false, nullptr)); } @@ -2547,7 +2590,8 @@ namespace nlsat { resolve_clause(b, *(jst.get_clause())); break; case justification::LAZY: - resolve_lazy_justification(b, *(jst.get_lazy())); + + resolve_lazy_justification(b, *(jst.get_lazy())); break; case justification::DECISION: SASSERT(m_num_marks == 0); @@ -2603,9 +2647,8 @@ namespace nlsat { TRACE(nlsat, tout << "new lemma:\n"; display(tout, m_lemma.size(), m_lemma.data()); tout << "\n"; tout << "found_decision: " << found_decision << "\n";); - if (m_check_lemmas) { - check_lemma(m_lemma.size(), m_lemma.data(), m_lemma_assumptions.get()); - } + // Note: Don't check m_lemma here - it's the result of resolution + // and may not be a tautology. Conflict lemmas are checked in resolve_lazy_justification. // if (m_log_lemmas) // log_lemma(std::cout, m_lemma.size(), m_lemma.data(), false); @@ -2772,6 +2815,8 @@ namespace nlsat { st.update("nlsat stages", m_stats.m_stages); st.update("nlsat simplifications", m_stats.m_simplifications); st.update("nlsat irrational assignments", m_stats.m_irrational_assignments); + st.update("levelwise calls", m_stats.m_levelwise_calls); + st.update("levelwise failures", m_stats.m_levelwise_failures); } void reset_statistics() { @@ -4373,10 +4418,19 @@ namespace nlsat { nlsat_params::collect_param_descrs(d); } - unsynch_mpq_manager & solver::qm() { + const assignment &solver::sample() const { + return m_imp->m_assignment; + } + + assignment &solver::sample() { + return m_imp->m_assignment; + } + + unsynch_mpq_manager &solver::qm() + { return m_imp->m_qm; } - + anum_manager & solver::am() { return m_imp->m_am; } @@ -4626,6 +4680,15 @@ namespace nlsat { m_imp->m_stats.m_simplifications++; } + void solver::record_levelwise_result(bool success) { + m_imp->m_stats.m_levelwise_calls++; + m_imp->m_last_conflict_used_lws = success; // Track for unsound lemma reporting + if (!success) { + m_imp->m_stats.m_levelwise_failures++; + // m_imp->m_apply_lws = false; // is it useful to throttle + } + } + bool solver::has_root_atom(clause const& c) const { return m_imp->has_root_atom(c); } @@ -4637,6 +4700,9 @@ namespace nlsat { assumption solver::join(assumption a, assumption b) { return (m_imp->m_asm.mk_join(static_cast(a), static_cast(b))); } - - + + bool solver::apply_levelwise() const { return m_imp->m_apply_lws; } + + unsigned solver::lws_spt_threshold() const { return m_imp->m_lws_spt_threshold; } + }; diff --git a/src/nlsat/nlsat_solver.h b/src/nlsat/nlsat_solver.h index ae403fc83..566a7c09e 100644 --- a/src/nlsat/nlsat_solver.h +++ b/src/nlsat/nlsat_solver.h @@ -194,6 +194,7 @@ namespace nlsat { assumption join(assumption a, assumption b); void inc_simplify(); + void record_levelwise_result(bool success); void add_bound(bound_constraint const& c); /** @@ -244,7 +245,10 @@ namespace nlsat { // ----------------------- void updt_params(params_ref const & p); static void collect_param_descrs(param_descrs & d); - + const assignment& sample() const; + assignment& sample(); + bool apply_levelwise() const; + unsigned lws_spt_threshold() const; void reset(); void collect_statistics(statistics & st); void reset_statistics(); @@ -294,8 +298,7 @@ namespace nlsat { std::ostream& display_assignment(std::ostream& out) const; std::ostream& display_var(std::ostream& out, unsigned j) const; - + }; }; - diff --git a/src/nlsat/nlsat_types.h b/src/nlsat/nlsat_types.h index 720803861..be8f625c9 100644 --- a/src/nlsat/nlsat_types.h +++ b/src/nlsat/nlsat_types.h @@ -22,7 +22,7 @@ Revision History: #include "util/buffer.h" #include "sat/sat_types.h" #include "util/z3_exception.h" - +#include namespace algebraic_numbers { class anum; class manager; @@ -143,6 +143,50 @@ namespace nlsat { struct eq_proc { bool operator()(root_atom const * a1, root_atom const * a2) const; }; }; + /** + \brief An indexed root expression designates the i-th real root of a (square-free) polynomial p when regarded as + a univariate over its maximal variable x after substituting the current values for variables < x. + + It is a lightweight value object used by projection / sample-cell algorithms. It does not own the polynomial. + */ + struct indexed_root { + poly * m_p; // polynomial whose root is referenced (assumed non-null) + unsigned m_index; // root index (0-based internally; corresponds to paper's i) + var m_var; // the main variable of m_p when treated as univariate + indexed_root(): m_p(nullptr), m_index(0), m_var(null_var) {} + indexed_root(poly* p, unsigned i, var x): m_p(p), m_index(i), m_var(x) {} + poly * p() const { return m_p; } + unsigned index() const { return m_index; } + var x() const { return m_var; } + bool is_null() const { return m_p == nullptr; } + // ordering helper (by variable then polynomial id then index) - total but arbitrary + struct lt { + pmanager & m_pm; + lt(pmanager & pm): m_pm(pm) {} + bool operator()(indexed_root const & a, indexed_root const & b) const { + if (a.m_var != b.m_var) return a.m_var < b.m_var; + if (a.m_p != b.m_p) return m_pm.id(a.m_p) < m_pm.id(b.m_p); + return a.m_index < b.m_index; + } + }; + struct hash_proc { + pmanager & m_pm; + hash_proc(pmanager & pm): m_pm(pm) {} + unsigned operator()(indexed_root const & r) const { + return combine(combine(m_pm.id(r.m_p), r.m_index), r.m_var); + } + static unsigned combine(unsigned a, unsigned b) { return a * 1315423911u + b + (a<<5) + (b>>2); } + }; + struct eq_proc { + pmanager & m_pm; + eq_proc(pmanager & pm): m_pm(pm) {} + bool operator()(indexed_root const & a, indexed_root const & b) const { + return a.m_var == b.m_var && a.m_index == b.m_index && a.m_p == b.m_p; } + }; + }; + + typedef std::vector indexed_root_vector; + inline ineq_atom * to_ineq_atom(atom * a) { SASSERT(a->is_ineq_atom()); return static_cast(a); } inline root_atom * to_root_atom(atom * a) { SASSERT(a->is_root_atom()); return static_cast(a); } inline ineq_atom const * to_ineq_atom(atom const * a) { SASSERT(a->is_ineq_atom()); return static_cast(a); } diff --git a/src/test/nlsat.cpp b/src/test/nlsat.cpp index 35e9278ce..f26995462 100644 --- a/src/test/nlsat.cpp +++ b/src/test/nlsat.cpp @@ -20,6 +20,7 @@ Notes: #include "nlsat/nlsat_interval_set.h" #include "nlsat/nlsat_evaluator.h" #include "nlsat/nlsat_solver.h" +#include "nlsat/levelwise.h" #include "util/util.h" #include "nlsat/nlsat_explain.h" #include "math/polynomial/polynomial_cache.h" @@ -27,6 +28,112 @@ Notes: #include #include +// Helper function to check if a point (given by counter_as) is inside a levelwise cell. +// Returns true if the point is inside ALL cell intervals, false otherwise. +// Prints diagnostic info about which constraint was violated (if any). +static bool is_point_inside_cell( + anum_manager& am, + polynomial::manager& pm, + std_vector const& cell, + nlsat::assignment const& counter_as, + bool verbose = true) +{ + for (unsigned level = 0; level < cell.size(); ++level) { + auto const& interval = cell[level]; + + // Skip whole-line sectors - no constraint + if (!interval.is_section() && interval.l_inf() && interval.u_inf()) + continue; + + // Get the value at this level from counter_as + if (!counter_as.is_assigned(level)) + continue; // Variable not assigned in counterexample + + scoped_anum counter_val(am); + am.set(counter_val, counter_as.value(level)); + + // Build partial assignment up to this level (exclusive) for root isolation + nlsat::assignment partial_as(am); + for (unsigned j = 0; j < level; ++j) { + if (counter_as.is_assigned(j)) { + scoped_anum v(am); + am.set(v, counter_as.value(j)); + partial_as.set(j, v); + } + } + + if (interval.is_section()) { + // Section: counter must equal the root + polynomial_ref section_p(interval.l, pm); + scoped_anum_vector roots(am); + am.isolate_roots(section_p, nlsat::undef_var_assignment(partial_as, level), roots); + + if (roots.size() >= interval.l_index) { + bool at_root = am.eq(counter_val, roots[interval.l_index - 1]); + if (!at_root) { + if (verbose) { + std::cout << " Level " << level << ": OUTSIDE (not at section root)\n"; + std::cout << " Section root = "; + am.display_decimal(std::cout, roots[interval.l_index - 1], 6); + std::cout << ", counter = "; + am.display_decimal(std::cout, counter_val, 6); + std::cout << "\n"; + } + return false; + } + } + } else { + // Sector: check lower and upper bounds + if (!interval.l_inf()) { + polynomial_ref lower_p(interval.l, pm); + scoped_anum_vector lower_roots(am); + am.isolate_roots(lower_p, nlsat::undef_var_assignment(partial_as, level), lower_roots); + + if (lower_roots.size() >= interval.l_index) { + int cmp = am.compare(counter_val, lower_roots[interval.l_index - 1]); + if (cmp <= 0) { // counter <= lower bound, violates (lower, ...) + if (verbose) { + std::cout << " Level " << level << ": OUTSIDE (at or below lower bound)\n"; + std::cout << " Lower bound = "; + am.display_decimal(std::cout, lower_roots[interval.l_index - 1], 6); + std::cout << ", counter = "; + am.display_decimal(std::cout, counter_val, 6); + std::cout << "\n"; + } + return false; + } + } + } + + if (!interval.u_inf()) { + polynomial_ref upper_p(interval.u, pm); + scoped_anum_vector upper_roots(am); + am.isolate_roots(upper_p, nlsat::undef_var_assignment(partial_as, level), upper_roots); + + if (upper_roots.size() >= interval.u_index) { + int cmp = am.compare(counter_val, upper_roots[interval.u_index - 1]); + if (cmp >= 0) { // counter >= upper bound, violates (..., upper) + if (verbose) { + std::cout << " Level " << level << ": OUTSIDE (at or above upper bound)\n"; + std::cout << " Upper bound = "; + am.display_decimal(std::cout, upper_roots[interval.u_index - 1], 6); + std::cout << ", counter = "; + am.display_decimal(std::cout, counter_val, 6); + std::cout << "\n"; + } + return false; + } + } + } + } + } + + if (verbose) { + std::cout << " Point is INSIDE the cell (all constraints satisfied)\n"; + } + return true; +} + nlsat::interval_set_ref tst_interval(nlsat::interval_set_ref const & s1, nlsat::interval_set_ref const & s2, unsigned expected_num_intervals, @@ -320,7 +427,7 @@ static void project(nlsat::solver& s, nlsat::explain& ex, nlsat::var x, unsigned static void project_fa(nlsat::solver& s, nlsat::explain& ex, nlsat::var x, unsigned num, nlsat::literal const* lits) { std::cout << "Project "; nlsat::scoped_literal_vector result(s); - ex.main_operator(num, lits, result); + ex.compute_conflict_explanation(num, lits, result); std::cout << "(or"; for (auto l : result) { s.display(std::cout << " ", l); @@ -360,6 +467,11 @@ static nlsat::literal mk_eq(nlsat::solver& s, nlsat::poly* p) { return s.mk_ineq_literal(nlsat::atom::EQ, 1, _p, is_even); } +static nlsat::literal mk_root_eq(nlsat::solver& s, nlsat::poly* p, nlsat::var x, unsigned i) { + nlsat::bool_var b = s.mk_root_atom(nlsat::atom::ROOT_EQ, x, i, p); + return nlsat::literal(b, false); +} + static void set_assignment_value(nlsat::assignment& as, anum_manager& am, nlsat::var v, rational const& val) { scoped_anum tmp(am); am.set(tmp, val.to_mpq()); @@ -849,7 +961,7 @@ void tst_nlsat_mv() { s.set_rvalues(assignment); nlsat::scoped_literal_vector result(s); - ex.main_operator(lits.size(), lits.data(), result); + ex.compute_conflict_explanation(lits.size(), lits.data(), result); std::cout << "main_operator root regression core:\n"; s.display(std::cout, lits.size(), lits.data()); @@ -959,11 +1071,1484 @@ x7 := 1 } +static void tst_polynomial_cache_mk_unique() { + params_ref ps; + reslimit rlim; + nlsat::solver s(rlim, ps, false); + nlsat::pmanager& pm = s.pm(); + polynomial::cache cache(pm); + nlsat::var x = s.mk_var(false); + polynomial_ref x_poly(pm); + x_poly = pm.mk_polynomial(x); + + polynomial_ref p(pm); + p = 2 * x_poly; + pm.display(std::cout, p); + std::cout << "\n"; + + ENSURE(!cache.contains(p.get())); + polynomial::polynomial* unique_p = cache.mk_unique(p.get()); + ENSURE(unique_p != nullptr); + ENSURE(cache.contains(unique_p)); + pm.display(std::cout, unique_p); + std::cout << "\n"; + + polynomial_ref p_again(pm); + p_again = 2 * x_poly; + ENSURE(cache.mk_unique(p_again.get()) == unique_p); + + polynomial_ref p_neg(pm); + p_neg = -2 * x_poly; + ENSURE(!cache.contains(p_neg.get())); + polynomial::polynomial* unique_p_neg = cache.mk_unique(p_neg.get()); + ENSURE(unique_p_neg != nullptr); + ENSURE(unique_p_neg != unique_p); + ENSURE(cache.contains(unique_p_neg)); + + polynomial_ref p_neg_again(pm); + p_neg_again = -2 * x_poly; + ENSURE(cache.mk_unique(p_neg_again.get()) == unique_p_neg); + + polynomial_ref ca(p, pm), cb(p_neg, pm); + pm.primitive(ca, x, ca); + pm.primitive(cb, x, cb); + pm.display(std::cout << "ca:", ca) << "\n"; + pm.display(std::cout << "cb:", cb) << "\n"; + +} + +static void tst_nullified_polynomial() { + params_ref ps; + reslimit rlim; + nlsat::solver s(rlim, ps, false); + nlsat::pmanager & pm = s.pm(); + anum_manager & am = s.am(); + polynomial::cache cache(pm); + + nlsat::var x0 = s.mk_var(false); + nlsat::var x1 = s.mk_var(false); + nlsat::var x2 = s.mk_var(false); + nlsat::var x3 = s.mk_var(false); + ENSURE(x0 < x1 && x1 < x2); + + polynomial_ref _x0(pm), _x1(pm), _x2(pm), _x3(pm); + _x0 = pm.mk_polynomial(x0); + _x1 = pm.mk_polynomial(x1); + _x2 = pm.mk_polynomial(x2); + _x3 = pm.mk_polynomial(x3); + + polynomial_ref p1(pm), p2(pm); + p1 = _x2 * _x1; + p2 = _x0; + p1 = p1 + p2; + p2 = _x3; + polynomial_ref_vector polys(pm); + polys.push_back(p1); + polys.push_back(p2); + nlsat::assignment as(am); + scoped_anum zero(am); + am.set(zero, 0); + as.set(x0, zero); + as.set(x1, zero); + as.set(x2, zero); + as.set(x3, zero); + s.set_rvalues(as); + + unsigned max_x = 0; + for (unsigned i = 0; i < polys.size(); ++i) { + unsigned lvl = pm.max_var(polys.get(i)); + if (lvl > max_x) + max_x = lvl; + } + ENSURE(max_x == x3); + + nlsat::levelwise lws(s, polys, max_x, s.sample(), pm, am, cache); + auto cell = lws.single_cell(); + ENSURE(lws.failed()); +} + +// Test case for unsound lemma lws2380 - comparing standard projection vs levelwise +// The issue: x7 is unconstrained in levelwise output but affects the section polynomial +static void tst_unsound_lws2380() { + enable_trace("nlsat_explain"); + + auto run_test = [](bool use_lws) { + std::cout << "=== tst_unsound_lws2380: " << (use_lws ? "Levelwise" : "Standard") << " projection (lws=" << use_lws << ") ===\n"; + params_ref ps; + ps.set_bool("lws", use_lws); + reslimit rlim; + nlsat::solver s(rlim, ps, false); + anum_manager & am = s.am(); + nlsat::pmanager & pm = s.pm(); + nlsat::assignment as(am); + nlsat::explain& ex = s.get_explain(); + + // Create 14 variables x0-x13 + nlsat::var x0 = s.mk_var(false); + nlsat::var x1 = s.mk_var(false); + nlsat::var x2 = s.mk_var(false); + nlsat::var x3 = s.mk_var(false); + nlsat::var x4 = s.mk_var(false); + nlsat::var x5 = s.mk_var(false); + nlsat::var x6 = s.mk_var(false); + nlsat::var x7 = s.mk_var(false); + nlsat::var x8 = s.mk_var(false); + nlsat::var x9 = s.mk_var(false); + nlsat::var x10 = s.mk_var(false); + nlsat::var x11 = s.mk_var(false); + nlsat::var x12 = s.mk_var(false); + nlsat::var x13 = s.mk_var(false); + + polynomial_ref _x0(pm), _x1(pm), _x2(pm), _x3(pm), _x4(pm), _x5(pm), _x6(pm); + polynomial_ref _x7(pm), _x8(pm), _x9(pm), _x10(pm), _x11(pm), _x12(pm), _x13(pm); + _x0 = pm.mk_polynomial(x0); + _x1 = pm.mk_polynomial(x1); + _x2 = pm.mk_polynomial(x2); + _x3 = pm.mk_polynomial(x3); + _x4 = pm.mk_polynomial(x4); + _x5 = pm.mk_polynomial(x5); + _x6 = pm.mk_polynomial(x6); + _x7 = pm.mk_polynomial(x7); + _x8 = pm.mk_polynomial(x8); + _x9 = pm.mk_polynomial(x9); + _x10 = pm.mk_polynomial(x10); + _x11 = pm.mk_polynomial(x11); + _x12 = pm.mk_polynomial(x12); + _x13 = pm.mk_polynomial(x13); + + // p[0]: x13 + polynomial_ref p0(pm); + p0 = _x13; + + // p[1]: 170*x8*x13 + x10*x11*x12 - x11*x12 - x7*x8*x12 + x5*x10*x11 + 184*x1*x10*x11 + // - x0*x10*x11 - x5*x11 - 184*x1*x11 + x0*x11 - x3*x8*x10 + x2*x8*x10 + // - 2*x10 - 184*x1*x7*x8 - x2*x8 + 2 + polynomial_ref p1(pm); + p1 = (170 * _x8 * _x13) + + (_x10 * _x11 * _x12) - + (_x11 * _x12) - + (_x7 * _x8 * _x12) + + (_x5 * _x10 * _x11) + + (184 * _x1 * _x10 * _x11) - + (_x0 * _x10 * _x11) - + (_x5 * _x11) - + (184 * _x1 * _x11) + + (_x0 * _x11) - + (_x3 * _x8 * _x10) + + (_x2 * _x8 * _x10) - + (2 * _x10) - + (184 * _x1 * _x7 * _x8) - + (_x2 * _x8) + + 2; + + std::cout << "p0: " << p0 << "\n"; + std::cout << "p1: " << p1 << "\n"; + + // Set sample: x0=-1, x1=-1, x2=0, x3=-1, x5=0, x7=0, x8=2, x10=0, x11=0, x12=2 + scoped_anum val(am); + am.set(val, -1); as.set(x0, val); + am.set(val, -1); as.set(x1, val); + am.set(val, 0); as.set(x2, val); + am.set(val, -1); as.set(x3, val); + am.set(val, 0); as.set(x4, val); + am.set(val, 0); as.set(x5, val); + am.set(val, 0); as.set(x6, val); + am.set(val, 0); as.set(x7, val); + am.set(val, 2); as.set(x8, val); + am.set(val, 0); as.set(x9, val); + am.set(val, 0); as.set(x10, val); + am.set(val, 0); as.set(x11, val); + am.set(val, 2); as.set(x12, val); + am.set(val, 1); as.set(x13, val); + s.set_rvalues(as); + + // Create literals for the two polynomials + nlsat::scoped_literal_vector lits(s); + lits.push_back(mk_gt(s, p0.get())); // x13 > 0 + lits.push_back(mk_gt(s, p1.get())); // p1 > 0 + + project_fa(s, ex, x13, lits.size(), lits.data()); + std::cout << "\n"; + }; + + run_test(false); // Standard projection + run_test(true); // Levelwise projection +} + +// Test case for unsound lemma - levelwise produces cell that's too large +// Input: 5 polynomials with max_var=x3, sample x0=-7, x1=-1, x2=1 +// Counterexample: x0=-4, x1=-8, x2=5, x3=6 +static void tst_unsound_lws_x3() { + std::cout << "=== tst_unsound_lws_x3 ===\n"; + params_ref ps; + ps.set_bool("lws", true); + reslimit rlim; + nlsat::solver s(rlim, ps, false); + anum_manager & am = s.am(); + nlsat::pmanager & pm = s.pm(); + nlsat::assignment sample_as(am); + nlsat::assignment counter_as(am); + polynomial::cache cache(pm); + + // Create 4 variables x0, x1, x2, x3 + nlsat::var x0 = s.mk_var(false); + nlsat::var x1 = s.mk_var(false); + nlsat::var x2 = s.mk_var(false); + nlsat::var x3 = s.mk_var(false); + + polynomial_ref _x0(pm), _x1(pm), _x2(pm), _x3(pm); + _x0 = pm.mk_polynomial(x0); + _x1 = pm.mk_polynomial(x1); + _x2 = pm.mk_polynomial(x2); + _x3 = pm.mk_polynomial(x3); + + // p[0]: x3 + x0 + polynomial_ref p0(pm); + p0 = _x3 + _x0; + + // p[1]: x1 + polynomial_ref p1(pm); + p1 = _x1; + + // p[2]: 9 x1^2 x3^2 - 6 x0 x1 x2 x3 + 3 x0 x1^2 x3 - 12 x1^2 x3 + x0^2 x2^2 - x0^2 x1 x2 - 4 x0 x1 x2 + 2 x0 x1^2 + 4 x1^2 + polynomial_ref p2(pm); + p2 = 9 * (_x1^2) * (_x3^2) + - 6 * _x0 * _x1 * _x2 * _x3 + + 3 * _x0 * (_x1^2) * _x3 + - 12 * (_x1^2) * _x3 + + (_x0^2) * (_x2^2) + - (_x0^2) * _x1 * _x2 + - 4 * _x0 * _x1 * _x2 + + 2 * _x0 * (_x1^2) + + 4 * (_x1^2); + + // p[3]: 9 x1^2 x3^2 - 6 x0 x1 x2 x3 + 6 x0 x1^2 x3 - 12 x1^2 x3 + x0^2 x2^2 - 2 x0^2 x1 x2 - 4 x0 x1 x2 + x0^2 x1^2 + 4 x0 x1^2 + 4 x1^2 + polynomial_ref p3(pm); + p3 = 9 * (_x1^2) * (_x3^2) + - 6 * _x0 * _x1 * _x2 * _x3 + + 6 * _x0 * (_x1^2) * _x3 + - 12 * (_x1^2) * _x3 + + (_x0^2) * (_x2^2) + - 2 * (_x0^2) * _x1 * _x2 + - 4 * _x0 * _x1 * _x2 + + (_x0^2) * (_x1^2) + + 4 * _x0 * (_x1^2) + + 4 * (_x1^2); + + // p[4]: x3 + x1 + x0 + polynomial_ref p4(pm); + p4 = _x3 + _x1 + _x0; + + std::cout << "p0: " << p0 << "\n"; + std::cout << "p1: " << p1 << "\n"; + std::cout << "p2: " << p2 << "\n"; + std::cout << "p3: " << p3 << "\n"; + std::cout << "p4: " << p4 << "\n\n"; + + // Set sample: x0=-7, x1=-1, x2=1, x3=? (need to pick a value in the cell) + // For the sample, we need an x3 value. Let's use x3=8 (which is > -x0 = 7, so p0 > 0) + scoped_anum val(am); + am.set(val, -7); sample_as.set(x0, val); + am.set(val, -1); sample_as.set(x1, val); + am.set(val, 1); sample_as.set(x2, val); + am.set(val, 8); sample_as.set(x3, val); + + // Counterexample: x0=-4, x1=-8, x2=5, x3=6 + am.set(val, -4); counter_as.set(x0, val); + am.set(val, -8); counter_as.set(x1, val); + am.set(val, 5); counter_as.set(x2, val); + am.set(val, 6); counter_as.set(x3, val); + + std::cout << "Sample point: x0=-7, x1=-1, x2=1, x3=8\n"; + std::cout << "Counterexample: x0=-4, x1=-8, x2=5, x3=6\n\n"; + + // Evaluate polynomials at sample + std::cout << "Polynomial signs at SAMPLE:\n"; + std::cout << " p0 sign: " << am.eval_sign_at(p0, sample_as) << "\n"; + std::cout << " p1 sign: " << am.eval_sign_at(p1, sample_as) << "\n"; + std::cout << " p2 sign: " << am.eval_sign_at(p2, sample_as) << "\n"; + std::cout << " p3 sign: " << am.eval_sign_at(p3, sample_as) << "\n"; + std::cout << " p4 sign: " << am.eval_sign_at(p4, sample_as) << "\n\n"; + + // Evaluate polynomials at counterexample + std::cout << "Polynomial signs at COUNTEREXAMPLE:\n"; + std::cout << " p0 sign: " << am.eval_sign_at(p0, counter_as) << "\n"; + std::cout << " p1 sign: " << am.eval_sign_at(p1, counter_as) << "\n"; + std::cout << " p2 sign: " << am.eval_sign_at(p2, counter_as) << "\n"; + std::cout << " p3 sign: " << am.eval_sign_at(p3, counter_as) << "\n"; + std::cout << " p4 sign: " << am.eval_sign_at(p4, counter_as) << "\n\n"; + + // Set solver assignment for levelwise (without x3) + am.set(val, -7); sample_as.set(x0, val); + am.set(val, -1); sample_as.set(x1, val); + am.set(val, 1); sample_as.set(x2, val); + s.set_rvalues(sample_as); + + // Build polynomial vector + polynomial_ref_vector polys(pm); + polys.push_back(p0); + polys.push_back(p1); + polys.push_back(p2); + polys.push_back(p3); + polys.push_back(p4); + + unsigned max_x = x3; + + // Print roots of each polynomial at sample + std::cout << "Roots of polynomials at sample (in x3):\n"; + for (unsigned i = 0; i < polys.size(); ++i) { + polynomial_ref p(polys.get(i), pm); + if (pm.max_var(p) != x3) { + std::cout << " p" << i << ": max_var is not x3, skipping\n"; + continue; + } + scoped_anum_vector roots(am); + am.isolate_roots(p, nlsat::undef_var_assignment(sample_as, x3), roots); + std::cout << " p" << i << " roots: "; + if (roots.empty()) { + std::cout << "(none)"; + } else { + for (unsigned j = 0; j < roots.size(); ++j) { + if (j > 0) std::cout << ", "; + am.display_decimal(std::cout, roots[j], 5); + } + } + std::cout << "\n"; + } + std::cout << "\n"; + + // Compute and evaluate resultant of p3 and p4 + std::cout << "Resultant of p3 and p4 (in x3):\n"; + polynomial_ref res_p3_p4(pm); + { + pm.resultant(p3, p4, x3, res_p3_p4); + std::cout << " Res(p3, p4) = "; + pm.display(std::cout, res_p3_p4); + std::cout << "\n"; + std::cout << " Sign at sample (x0=-7, x1=-1, x2=1): " << am.eval_sign_at(res_p3_p4, sample_as) << "\n"; + std::cout << " Sign at counter (x0=-4, x1=-8, x2=5): " << am.eval_sign_at(res_p3_p4, counter_as) << "\n"; + + // Check roots of the resultant at x2 level (parametric in x0, x1) + std::cout << " Roots at sample x0,x1 (in x2): "; + scoped_anum_vector res_roots(am); + nlsat::assignment partial_sample(am); + scoped_anum val(am); + am.set(val, -7); partial_sample.set(x0, val); + am.set(val, -1); partial_sample.set(x1, val); + am.isolate_roots(res_p3_p4, nlsat::undef_var_assignment(partial_sample, x2), res_roots); + for (unsigned j = 0; j < res_roots.size(); ++j) { + if (j > 0) std::cout << ", "; + am.display_decimal(std::cout, res_roots[j], 5); + } + std::cout << "\n"; + + // Check roots at counterexample x0,x1 + std::cout << " Roots at counter x0,x1 (in x2): "; + nlsat::assignment partial_counter(am); + am.set(val, -4); partial_counter.set(x0, val); + am.set(val, -8); partial_counter.set(x1, val); + scoped_anum_vector res_roots_counter(am); + am.isolate_roots(res_p3_p4, nlsat::undef_var_assignment(partial_counter, x2), res_roots_counter); + for (unsigned j = 0; j < res_roots_counter.size(); ++j) { + if (j > 0) std::cout << ", "; + am.display_decimal(std::cout, res_roots_counter[j], 5); + } + std::cout << "\n"; + + // Compute and check discriminant of Res(p3,p4) in x2 + std::cout << "\n Discriminant of Res(p3,p4) in x2:\n"; + polynomial_ref disc_res(pm); + pm.discriminant(res_p3_p4, x2, disc_res); + std::cout << " Disc = "; + pm.display(std::cout, disc_res); + std::cout << "\n"; + std::cout << " Sign at sample (x0=-7, x1=-1): " << am.eval_sign_at(disc_res, sample_as) << "\n"; + std::cout << " Sign at counter (x0=-4, x1=-8): " << am.eval_sign_at(disc_res, counter_as) << "\n"; + } + std::cout << "\n"; + + std::cout << "Running levelwise with max_x = x3\n"; + + // Run levelwise + nlsat::levelwise lws(s, polys, max_x, s.sample(), pm, am, cache); + auto cell = lws.single_cell(); + + std::cout << "Levelwise " << (lws.failed() ? "FAILED" : "succeeded") << "\n"; + std::cout << "Cell intervals (count=" << cell.size() << "):\n"; + for (auto const& interval : cell) { + nlsat::display(std::cout << " ", s, interval) << "\n"; + } + + // Evaluate cell bounds at counterexample to check if counterexample is in cell + std::cout << "\n--- Checking if counterexample is in cell ---\n"; + std::cout << "For a SECTOR (lower_root, upper_root), variable x satisfies:\n"; + std::cout << " x > lower_root AND x < upper_root\n\n"; + + // For univariate evaluation, we need to substitute lower vars + // Level 0: x0 interval, evaluate at x0=-4 + // Level 1: x1 interval (parametric in x0), evaluate at (x0=-4, x1=-8) + // Level 2: x2 interval (parametric in x0,x1), evaluate at (x0=-4,x1=-8,x2=5) + + bool counterexample_outside_cell = false; + + for (unsigned i = 0; i < cell.size(); ++i) { + auto const& interval = cell[i]; + nlsat::var level = i; + std::cout << "Level " << level << ":\n"; + + // Build assignment up to this level (exclusive) for root isolation + nlsat::assignment partial_as(am); + scoped_anum val(am); + if (level > 0) { am.set(val, -4); partial_as.set(x0, val); } + if (level > 1) { am.set(val, -8); partial_as.set(x1, val); } + if (level > 2) { am.set(val, 5); partial_as.set(x2, val); } + + scoped_anum counter_val(am); + if (level == 0) am.set(counter_val, -4); + else if (level == 1) am.set(counter_val, -8); + else if (level == 2) am.set(counter_val, 5); + + if (interval.is_section()) { + std::cout << " Section case\n"; + } else { + // Isolate roots and check bounds + if (!interval.l_inf()) { + polynomial_ref lower_p(interval.l, pm); + scoped_anum_vector lower_roots(am); + am.isolate_roots(lower_p, nlsat::undef_var_assignment(partial_as, level), lower_roots); + if (lower_roots.size() >= interval.l_index) { + std::cout << " Lower root (root[" << interval.l_index << "]): "; + am.display_decimal(std::cout, lower_roots[interval.l_index - 1], 10); + std::cout << "\n"; + std::cout << " Counter x" << level << " = "; + am.display_decimal(std::cout, counter_val, 10); + int cmp = am.compare(counter_val, lower_roots[interval.l_index - 1]); + std::cout << " -> " << (cmp > 0 ? "ABOVE" : (cmp < 0 ? "BELOW" : "EQUAL")) << " lower bound\n"; + if (cmp <= 0) counterexample_outside_cell = true; + } + } + if (!interval.u_inf()) { + polynomial_ref upper_p(interval.u, pm); + scoped_anum_vector upper_roots(am); + am.isolate_roots(upper_p, nlsat::undef_var_assignment(partial_as, level), upper_roots); + if (upper_roots.size() >= interval.u_index) { + std::cout << " Upper root (root[" << interval.u_index << "]): "; + am.display_decimal(std::cout, upper_roots[interval.u_index - 1], 10); + std::cout << "\n"; + std::cout << " Counter x" << level << " = "; + am.display_decimal(std::cout, counter_val, 10); + int cmp = am.compare(counter_val, upper_roots[interval.u_index - 1]); + std::cout << " -> " << (cmp > 0 ? "ABOVE" : (cmp < 0 ? "BELOW" : "EQUAL")) << " upper bound\n"; + if (cmp >= 0) counterexample_outside_cell = true; + } + } + } + std::cout << "\n"; + } + + // The counterexample has different polynomial signs than the sample. + // For a sound cell, the counterexample must be OUTSIDE the cell. + ENSURE(counterexample_outside_cell); + std::cout << "SUCCESS: Counterexample is OUTSIDE the cell (cell is sound)\n"; + + std::cout << "=== END tst_unsound_lws_x3 ===\n\n"; +} + +// Test case for unsound lemma from From_T2__n-46.t2__p4432_terminationG_0.smt2 +// This test calls levelwise with the input polynomials and analyzes what was missed. +// +// Input polynomials passed to levelwise: +// p[0]: x2 +// p[1]: x2 x6^2 + x0 x5 x6 + x2 x4 x6 + x2 x3 x6 - x0 x6 - x0 x4 +// p[2]: x2 x6^2 x7 + x0 x5 x6 x7 + x2 x4 x6 x7 + x2 x3 x6 x7 - x0 x6 x7 - x0 x4 x7 + 2 x6 + 2 x4 +// p[3]: x6 +// p[4]: x2 x6 x7 - 2 +// +// Sample: x0=1, x1=2, x2=1, x3=-1, x4=-1, x5=1, x6=7/8 +// Counterexample: x0=1, x2=1, x3=0, x4=-9, x5=0, x6=5, x7=0 +// +// Unsound lemma produced: +// !(x0 > 0) or !(x2 > 0) or !(x4 < root[1](2 x2 x4 + x0)) or +// !(x5 = root[1](x0 x5 + x2 x3)) or +// !(x0^2 x5^2 + ... > 0) or !(2 x2 > 0) or +// !(4 x2 x6 + x0 x5 + 2 x2 x4 + x2 x3 - x0 > 0) or +// !(2 x2 x6^2 + x0 x5 x6 + 2 x2 x4 x6 + x2 x3 x6 - x0 x6 - x0 x4 < 0) or +// x7 < root[1](x2 x6^2 x7 + ... + 2 x6 + 2 x4) or +// !(x7 < root[1](x2 x6 x7 - 2)) +static void tst_unsound_lws_n46() { + std::cout << "=== tst_unsound_lws_n46 ===\n"; + + params_ref ps; + ps.set_bool("lws", true); + reslimit rlim; + nlsat::solver s(rlim, ps, false); + anum_manager & am = s.am(); + nlsat::pmanager & pm = s.pm(); + nlsat::assignment sample_as(am); + nlsat::assignment counter_as(am); + polynomial::cache cache(pm); + + // Create 8 variables x0-x7 + nlsat::var x0 = s.mk_var(false); + nlsat::var x1 = s.mk_var(false); + nlsat::var x2 = s.mk_var(false); + nlsat::var x3 = s.mk_var(false); + nlsat::var x4 = s.mk_var(false); + nlsat::var x5 = s.mk_var(false); + nlsat::var x6 = s.mk_var(false); + nlsat::var x7 = s.mk_var(false); + + polynomial_ref _x0(pm), _x1(pm), _x2(pm), _x3(pm), _x4(pm), _x5(pm), _x6(pm), _x7(pm); + _x0 = pm.mk_polynomial(x0); + _x1 = pm.mk_polynomial(x1); + _x2 = pm.mk_polynomial(x2); + _x3 = pm.mk_polynomial(x3); + _x4 = pm.mk_polynomial(x4); + _x5 = pm.mk_polynomial(x5); + _x6 = pm.mk_polynomial(x6); + _x7 = pm.mk_polynomial(x7); + (void)_x1; // unused + + // Input polynomials passed to levelwise (from unsound_log) + // p[0]: x2 + polynomial_ref p0(pm); + p0 = _x2; + + // p[1]: x2 x6^2 + x0 x5 x6 + x2 x4 x6 + x2 x3 x6 - x0 x6 - x0 x4 + polynomial_ref p1(pm); + p1 = _x2 * (_x6^2) + _x0 * _x5 * _x6 + _x2 * _x4 * _x6 + _x2 * _x3 * _x6 - _x0 * _x6 - _x0 * _x4; + + // p[2]: x2 x6^2 x7 + x0 x5 x6 x7 + x2 x4 x6 x7 + x2 x3 x6 x7 - x0 x6 x7 - x0 x4 x7 + 2 x6 + 2 x4 + polynomial_ref p2(pm); + p2 = _x2 * (_x6^2) * _x7 + _x0 * _x5 * _x6 * _x7 + _x2 * _x4 * _x6 * _x7 + + _x2 * _x3 * _x6 * _x7 - _x0 * _x6 * _x7 - _x0 * _x4 * _x7 + + 2 * _x6 + 2 * _x4; + + // p[3]: x6 + polynomial_ref p3(pm); + p3 = _x6; + + // p[4]: x2 x6 x7 - 2 + polynomial_ref p4(pm); + p4 = _x2 * _x6 * _x7 - 2; + + std::cout << "Input polynomials:\n"; + std::cout << " p0: " << p0 << "\n"; + std::cout << " p1: " << p1 << "\n"; + std::cout << " p2: " << p2 << "\n"; + std::cout << " p3: " << p3 << "\n"; + std::cout << " p4: " << p4 << "\n\n"; + + // Set sample point: x0=1, x1=2, x2=1, x3=-1, x4=-1, x5=1, x6=7/8 + scoped_anum val(am); + rational q(7, 8); // 0.875 = 7/8 + am.set(val, 1); sample_as.set(x0, val); + am.set(val, 2); sample_as.set(x1, val); + am.set(val, 1); sample_as.set(x2, val); + am.set(val, -1); sample_as.set(x3, val); + am.set(val, -1); sample_as.set(x4, val); + am.set(val, 1); sample_as.set(x5, val); + am.set(val, q.to_mpq()); sample_as.set(x6, val); + + std::cout << "Sample point: x0=1, x1=2, x2=1, x3=-1, x4=-1, x5=1, x6=7/8\n"; + + // Set counterexample: x0=1, x2=1, x3=0, x4=-9, x5=0, x6=5, x7=0 + am.set(val, 1); counter_as.set(x0, val); + am.set(val, 0); counter_as.set(x1, val); + am.set(val, 1); counter_as.set(x2, val); + am.set(val, 0); counter_as.set(x3, val); + am.set(val, -9); counter_as.set(x4, val); + am.set(val, 0); counter_as.set(x5, val); + am.set(val, 5); counter_as.set(x6, val); + am.set(val, 0); counter_as.set(x7, val); + + std::cout << "Counterexample: x0=1, x2=1, x3=0, x4=-9, x5=0, x6=5, x7=0\n\n"; + + // Set solver assignment for levelwise + s.set_rvalues(sample_as); + + // Build polynomial vector + polynomial_ref_vector polys(pm); + polys.push_back(p0); + polys.push_back(p1); + polys.push_back(p2); + polys.push_back(p3); + polys.push_back(p4); + + nlsat::var max_x = x7; + + // Run levelwise + std::cout << "Running levelwise with max_x = x7...\n"; + nlsat::levelwise lws(s, polys, max_x, s.sample(), pm, am, cache); + auto cell = lws.single_cell(); + + std::cout << "Levelwise " << (lws.failed() ? "FAILED" : "succeeded") << "\n"; + std::cout << "Cell intervals:\n"; + for (unsigned i = 0; i < cell.size(); ++i) { + std::cout << " Level " << i << ": "; + nlsat::display(std::cout, s, cell[i]) << "\n"; + } + + // Print the lemma produced by levelwise + std::cout << "\n--- LEMMA from levelwise ---\n"; + for (unsigned i = 0; i < cell.size(); ++i) { + auto const& interval = cell[i]; + if (interval.section) { + std::cout << "!(x" << i << " = root[" << interval.l_index << "]("; + pm.display(std::cout, interval.l) << "))\n"; + } else { + if (!interval.l_inf()) { + std::cout << "!(x" << i << " > root[" << interval.l_index << "]("; + pm.display(std::cout, interval.l) << "))\n"; + } + if (!interval.u_inf()) { + std::cout << "!(x" << i << " < root[" << interval.u_index << "]("; + pm.display(std::cout, interval.u) << "))\n"; + } + } + } + std::cout << "--- END LEMMA ---\n\n"; + + // Test for the discriminant projection fix: + // + // BUG: When p1 = (x6-1)^2 at the sample (a double root), the discriminant of p1 + // is zero. The function compute_omit_disc_for_same_boundary() incorrectly omitted + // this discriminant because it only checked if p1(sample) != 0, not if disc(p1) = 0. + // + // FIX: Now we also check if disc(p) = 0 at sample, and if so, we keep the discriminant. + // This causes the discriminant's root polynomial (x2*x4 + x0) to be projected, + // creating a section at level 4 that excludes the counterexample. + + std::cout << "=== Verifying discriminant projection fix ===\n"; + + // Check 1: Level 4 should now be a SECTION (not a sector as before the fix) + // The discriminant of p1 w.r.t. x6 has a factor (x2*x4 + x0) that becomes the section + ENSURE(cell.size() > 4); + ENSURE(cell[4].section); // Level 4 must be a section after the fix + std::cout << "Level 4 is a section: " << (cell[4].section ? "YES (FIX WORKING)" : "NO (BUG!)") << "\n"; + + // Check 2: The section polynomial at level 4 should be x2*x4 + x0 (or equivalent) + // At sample: x2=1, x0=1, so root is x4 = -x0/x2 = -1 (matches sample x4=-1) + polynomial_ref section_poly(cell[4].l, pm); + std::cout << "Level 4 section polynomial: " << section_poly << "\n"; + + // Check 3: Verify the counterexample is OUTSIDE the cell + // At counterexample: x2=1, x0=1, so section root is x4 = -1 + // But counterexample has x4 = -9, which is NOT equal to -1 + // Therefore the literal !(x4 = root[1](...)) is TRUE, making the lemma sound + + polynomial_ref x4_section(pm); + x4_section = _x2 * _x4 + _x0; // Expected section polynomial + scoped_anum_vector roots_x4(am); + am.isolate_roots(x4_section, nlsat::undef_var_assignment(counter_as, x4), roots_x4); + + std::cout << "At counterexample:\n"; + std::cout << " Section polynomial: x2*x4 + x0 = x4 + 1\n"; + std::cout << " Section root: x4 = "; + if (!roots_x4.empty()) am.display_decimal(std::cout, roots_x4[0], 6); + std::cout << "\n"; + std::cout << " Counterexample x4 = -9\n"; + + bool x4_at_section = !roots_x4.empty() && am.eq(counter_as.value(x4), roots_x4[0]); + std::cout << " Is x4=-9 equal to section root? " << (x4_at_section ? "YES" : "NO") << "\n"; + + // The fix ensures x4_at_section is FALSE, meaning the counterexample is OUTSIDE the cell + ENSURE(!x4_at_section); // Counterexample must NOT satisfy the section constraint + + std::cout << "\n=== FIX VERIFIED: Counterexample is outside the cell ===\n"; + std::cout << "The lemma literal !(x4 = root[1](x2*x4 + x0)) is TRUE at counterexample.\n"; + std::cout << "Therefore the lemma is SOUND (disjunction has a true literal).\n"; + + std::cout << "=== END tst_unsound_lws_n46 ===\n\n"; +} + +// Test case for unsound lemma from From_AProVE_2014__Et4-rec.jar-obl-8__p28996_terminationG_0.smt2 +// Sample: x0=4, x1=5, x2=3.5, x3=-8, x4=5 +// Counterexample: x0=5, x3=3, x4=0, x5=0 +// Polynomials: +// p[0]: x0 +// p[1]: 2 x0 x4^2 + 2 x3 x4 - x0 x4 - 2 x3 +// p[2]: 2 x0 x4^2 x5 + 2 x3 x4 x5 - x0 x4 x5 - 2 x3 x5 + 4 x3 x4^2 + 9 x0 x3 x4 - 26 x3 x4 - 3 x0 x4 +// p[3]: x5 - 9 +static void tst_unsound_lws_et4() { + std::cout << "=== tst_unsound_lws_et4 ===\n"; + params_ref ps; + ps.set_bool("lws", true); + reslimit rlim; + nlsat::solver s(rlim, ps, false); + anum_manager & am = s.am(); + nlsat::pmanager & pm = s.pm(); + nlsat::assignment sample_as(am); + nlsat::assignment counter_as(am); + polynomial::cache cache(pm); + + // Create 6 variables x0, x1, x2, x3, x4, x5 + nlsat::var x0 = s.mk_var(false); + nlsat::var x1 = s.mk_var(false); + nlsat::var x2 = s.mk_var(false); + nlsat::var x3 = s.mk_var(false); + nlsat::var x4 = s.mk_var(false); + nlsat::var x5 = s.mk_var(false); + + polynomial_ref _x0(pm), _x1(pm), _x2(pm), _x3(pm), _x4(pm), _x5(pm); + _x0 = pm.mk_polynomial(x0); + _x1 = pm.mk_polynomial(x1); + _x2 = pm.mk_polynomial(x2); + _x3 = pm.mk_polynomial(x3); + _x4 = pm.mk_polynomial(x4); + _x5 = pm.mk_polynomial(x5); + + // p[0]: x0 + polynomial_ref p0(pm); + p0 = _x0; + + // p[1]: 2 x0 x4^2 + 2 x3 x4 - x0 x4 - 2 x3 + polynomial_ref p1(pm); + p1 = 2 * _x0 * (_x4^2) + 2 * _x3 * _x4 - _x0 * _x4 - 2 * _x3; + + // p[2]: 2 x0 x4^2 x5 + 2 x3 x4 x5 - x0 x4 x5 - 2 x3 x5 + 4 x3 x4^2 + 9 x0 x3 x4 - 26 x3 x4 - 3 x0 x4 + polynomial_ref p2(pm); + p2 = 2 * _x0 * (_x4^2) * _x5 + + 2 * _x3 * _x4 * _x5 + - _x0 * _x4 * _x5 + - 2 * _x3 * _x5 + + 4 * _x3 * (_x4^2) + + 9 * _x0 * _x3 * _x4 + - 26 * _x3 * _x4 + - 3 * _x0 * _x4; + + // p[3]: x5 - 9 + polynomial_ref p3(pm); + p3 = _x5 - 9; + + std::cout << "p0: " << p0 << "\n"; + std::cout << "p1: " << p1 << "\n"; + std::cout << "p2: " << p2 << "\n"; + std::cout << "p3: " << p3 << "\n\n"; + + // Sample: x0=4, x1=5, x2=3.5, x3=-8, x4=5 + scoped_anum val(am); + am.set(val, 4); sample_as.set(x0, val); + am.set(val, 5); sample_as.set(x1, val); + rational q(7, 2); // 3.5 + am.set(val, q.to_mpq()); sample_as.set(x2, val); + am.set(val, -8); sample_as.set(x3, val); + am.set(val, 5); sample_as.set(x4, val); + // x5 not assigned (top level) + + // Counterexample: x0=5, x3=3, x4=0, x5=0 + am.set(val, 5); counter_as.set(x0, val); + am.set(val, 5); counter_as.set(x1, val); // use same as sample + am.set(val, q.to_mpq()); counter_as.set(x2, val); // use same as sample + am.set(val, 3); counter_as.set(x3, val); + am.set(val, 0); counter_as.set(x4, val); + am.set(val, 0); counter_as.set(x5, val); + + std::cout << "Sample point: x0=4, x1=5, x2=3.5, x3=-8, x4=5\n"; + std::cout << "Counterexample: x0=5, x3=3, x4=0, x5=0\n\n"; + + // Evaluate polynomials at sample (need to set x5 for evaluation) + scoped_anum sample_x5(am); + am.set(sample_x5, 0); // pick some value in the cell + nlsat::assignment sample_full(am); + am.set(val, 4); sample_full.set(x0, val); + am.set(val, 5); sample_full.set(x1, val); + am.set(val, q.to_mpq()); sample_full.set(x2, val); + am.set(val, -8); sample_full.set(x3, val); + am.set(val, 5); sample_full.set(x4, val); + am.set(val, 0); sample_full.set(x5, val); + + std::cout << "Polynomial signs at SAMPLE (with x5=0):\n"; + std::cout << " p0 sign: " << am.eval_sign_at(p0, sample_full) << "\n"; + std::cout << " p1 sign: " << am.eval_sign_at(p1, sample_full) << "\n"; + std::cout << " p2 sign: " << am.eval_sign_at(p2, sample_full) << "\n"; + std::cout << " p3 sign: " << am.eval_sign_at(p3, sample_full) << "\n\n"; + + // Evaluate polynomials at counterexample + std::cout << "Polynomial signs at COUNTEREXAMPLE:\n"; + std::cout << " p0 sign: " << am.eval_sign_at(p0, counter_as) << "\n"; + std::cout << " p1 sign: " << am.eval_sign_at(p1, counter_as) << "\n"; + std::cout << " p2 sign: " << am.eval_sign_at(p2, counter_as) << "\n"; + std::cout << " p3 sign: " << am.eval_sign_at(p3, counter_as) << "\n\n"; + + // Set solver assignment for levelwise (without x5) + s.set_rvalues(sample_as); + + // Build polynomial vector + polynomial_ref_vector polys(pm); + polys.push_back(p0); + polys.push_back(p1); + polys.push_back(p2); + polys.push_back(p3); + + unsigned max_x = x5; + + // Print roots of each polynomial at sample + std::cout << "Roots of polynomials at sample (in x5):\n"; + for (unsigned i = 0; i < polys.size(); ++i) { + polynomial_ref p(polys.get(i), pm); + if (pm.max_var(p) != x5) { + std::cout << " p" << i << ": max_var is not x5, skipping\n"; + continue; + } + scoped_anum_vector roots(am); + am.isolate_roots(p, nlsat::undef_var_assignment(sample_as, x5), roots); + std::cout << " p" << i << " roots: "; + if (roots.empty()) { + std::cout << "(none)"; + } else { + for (unsigned j = 0; j < roots.size(); ++j) { + if (j > 0) std::cout << ", "; + am.display_decimal(std::cout, roots[j], 5); + } + } + std::cout << "\n"; + } + std::cout << "\n"; + + std::cout << "Running levelwise with max_x = x5\n"; + + // Run levelwise + nlsat::levelwise lws(s, polys, max_x, s.sample(), pm, am, cache); + auto cell = lws.single_cell(); + + std::cout << "Levelwise " << (lws.failed() ? "FAILED" : "succeeded") << "\n"; + std::cout << "Cell intervals (count=" << cell.size() << "):\n"; + for (auto const& interval : cell) { + nlsat::display(std::cout << " ", s, interval) << "\n"; + } + + // Evaluate cell bounds at counterexample to check if counterexample is in cell + std::cout << "\n--- Checking if counterexample is in cell ---\n"; + + bool counterexample_outside_cell = false; + + for (unsigned i = 0; i < cell.size(); ++i) { + auto const& interval = cell[i]; + nlsat::var level = i; + std::cout << "Level " << level << " (x" << level << "):\n"; + + // Build assignment up to this level (exclusive) for root isolation + nlsat::assignment partial_as(am); + scoped_anum val(am); + if (level > 0) { am.set(val, 5); partial_as.set(x0, val); } // counter x0 + if (level > 1) { am.set(val, 5); partial_as.set(x1, val); } + if (level > 2) { am.set(val, q.to_mpq()); partial_as.set(x2, val); } + if (level > 3) { am.set(val, 3); partial_as.set(x3, val); } // counter x3 + if (level > 4) { am.set(val, 0); partial_as.set(x4, val); } // counter x4 + + scoped_anum counter_val(am); + if (level == 0) am.set(counter_val, 5); // x0 + else if (level == 1) am.set(counter_val, 5); + else if (level == 2) am.set(counter_val, q.to_mpq()); + else if (level == 3) am.set(counter_val, 3); // x3 + else if (level == 4) am.set(counter_val, 0); // x4 + else if (level == 5) am.set(counter_val, 0); // x5 + + if (interval.is_section()) { + std::cout << " Section case\n"; + } else { + // Isolate roots and check bounds + if (!interval.l_inf()) { + polynomial_ref lower_p(interval.l, pm); + scoped_anum_vector lower_roots(am); + am.isolate_roots(lower_p, nlsat::undef_var_assignment(partial_as, level), lower_roots); + if (lower_roots.size() >= interval.l_index) { + std::cout << " Lower root (root[" << interval.l_index << "]): "; + am.display_decimal(std::cout, lower_roots[interval.l_index - 1], 10); + std::cout << "\n"; + std::cout << " Counter x" << level << " = "; + am.display_decimal(std::cout, counter_val, 10); + int cmp = am.compare(counter_val, lower_roots[interval.l_index - 1]); + std::cout << " -> " << (cmp > 0 ? "ABOVE" : (cmp < 0 ? "BELOW" : "EQUAL")) << " lower bound\n"; + if (cmp <= 0) counterexample_outside_cell = true; + } + } + if (!interval.u_inf()) { + polynomial_ref upper_p(interval.u, pm); + scoped_anum_vector upper_roots(am); + am.isolate_roots(upper_p, nlsat::undef_var_assignment(partial_as, level), upper_roots); + if (upper_roots.size() >= interval.u_index) { + std::cout << " Upper root (root[" << interval.u_index << "]): "; + am.display_decimal(std::cout, upper_roots[interval.u_index - 1], 10); + std::cout << "\n"; + std::cout << " Counter x" << level << " = "; + am.display_decimal(std::cout, counter_val, 10); + int cmp = am.compare(counter_val, upper_roots[interval.u_index - 1]); + std::cout << " -> " << (cmp > 0 ? "ABOVE" : (cmp < 0 ? "BELOW" : "EQUAL")) << " upper bound\n"; + if (cmp >= 0) counterexample_outside_cell = true; + } + } + } + std::cout << "\n"; + } + + // The counterexample has different polynomial signs than the sample. + // For a sound cell, the counterexample must be OUTSIDE the cell. + ENSURE(counterexample_outside_cell); + std::cout << "SUCCESS: Counterexample is OUTSIDE the cell (cell is sound)\n"; + + std::cout << "=== END tst_unsound_lws_et4 ===\n\n"; +} + +// Test case for unsound lemma with disc=0 at sample for same_boundary_poly sector case +// Polynomials: +// p[0]: - x2^3 + x1^3 x2 - x1^4 +// p[1]: - x2^3 x3 + x1^3 x2 x3 - x1^4 x3 - x2^3 + x1^3 x2^2 + 4 x1^4 x2 - x1^3 x2 - x1^5 - x1^4 +// p[2]: x3 +// Sample: x0=1, x1=1, x2=1 +// Counterexample: x1=12, x2=16, x3=0 +static void tst_unsound_lws_disc_zero() { + std::cout << "=== tst_unsound_lws_disc_zero ===\n"; + + params_ref ps; + ps.set_bool("lws", true); + reslimit rlim; + nlsat::solver s(rlim, ps, false); + anum_manager & am = s.am(); + nlsat::pmanager & pm = s.pm(); + nlsat::assignment sample_as(am); + nlsat::assignment counter_as(am); + polynomial::cache cache(pm); + + // Create 4 variables x0-x3 + nlsat::var x0 = s.mk_var(false); + nlsat::var x1 = s.mk_var(false); + nlsat::var x2 = s.mk_var(false); + nlsat::var x3 = s.mk_var(false); + + polynomial_ref _x0(pm), _x1(pm), _x2(pm), _x3(pm); + _x0 = pm.mk_polynomial(x0); + _x1 = pm.mk_polynomial(x1); + _x2 = pm.mk_polynomial(x2); + _x3 = pm.mk_polynomial(x3); + (void)_x0; // unused + + // p[0]: - x2^3 + x1^3 x2 - x1^4 + polynomial_ref p0(pm); + p0 = -(_x2^3) + (_x1^3) * _x2 - (_x1^4); + + // p[1]: - x2^3 x3 + x1^3 x2 x3 - x1^4 x3 - x2^3 + x1^3 x2^2 + 4 x1^4 x2 - x1^3 x2 - x1^5 - x1^4 + polynomial_ref p1(pm); + p1 = -(_x2^3) * _x3 + (_x1^3) * _x2 * _x3 - (_x1^4) * _x3 + - (_x2^3) + (_x1^3) * (_x2^2) + 4 * (_x1^4) * _x2 - (_x1^3) * _x2 - (_x1^5) - (_x1^4); + + // p[2]: x3 + polynomial_ref p2(pm); + p2 = _x3; + + std::cout << "Input polynomials:\n"; + std::cout << " p0: " << p0 << "\n"; + std::cout << " p1: " << p1 << "\n"; + std::cout << " p2: " << p2 << "\n\n"; + + // Set sample point: x0=1, x1=1, x2=1 + scoped_anum val(am); + am.set(val, 1); sample_as.set(x0, val); + am.set(val, 1); sample_as.set(x1, val); + am.set(val, 1); sample_as.set(x2, val); + + std::cout << "Sample point: x0=1, x1=1, x2=1\n"; + + // Set counterexample: x1=12, x2=16, x3=0 + am.set(val, 1); counter_as.set(x0, val); + am.set(val, 12); counter_as.set(x1, val); + am.set(val, 16); counter_as.set(x2, val); + am.set(val, 0); counter_as.set(x3, val); + + std::cout << "Counterexample: x1=12, x2=16, x3=0\n\n"; + + // Set solver assignment for levelwise + s.set_rvalues(sample_as); + + // Build polynomial vector + polynomial_ref_vector polys(pm); + polys.push_back(p0); + polys.push_back(p1); + polys.push_back(p2); + + nlsat::var max_x = x3; + + // Run levelwise + std::cout << "Running levelwise with max_x = x3...\n"; + nlsat::levelwise lws(s, polys, max_x, s.sample(), pm, am, cache); + auto cell = lws.single_cell(); + + std::cout << "Levelwise " << (lws.failed() ? "FAILED" : "succeeded") << "\n"; + std::cout << "Cell intervals:\n"; + for (unsigned i = 0; i < cell.size(); ++i) { + std::cout << " Level " << i << ": "; + nlsat::display(std::cout, s, cell[i]) << "\n"; + } + + // Print the lemma produced by levelwise + std::cout << "\n--- LEMMA from levelwise ---\n"; + for (unsigned i = 0; i < cell.size(); ++i) { + auto const& interval = cell[i]; + if (interval.section) { + std::cout << "!(x" << i << " = root[" << interval.l_index << "]("; + pm.display(std::cout, interval.l) << "))\n"; + } else { + if (!interval.l_inf()) { + std::cout << "!(x" << i << " > root[" << interval.l_index << "]("; + pm.display(std::cout, interval.l) << "))\n"; + } + if (!interval.u_inf()) { + std::cout << "!(x" << i << " < root[" << interval.u_index << "]("; + pm.display(std::cout, interval.u) << "))\n"; + } + } + } + std::cout << "--- END LEMMA ---\n\n"; + + // Check sign of p0 at sample vs counterexample + std::cout << "=== Checking sign of p0 (projection poly) ===\n"; + sign p0_sample = am.eval_sign_at(p0, sample_as); + sign p0_counter = am.eval_sign_at(p0, counter_as); + std::cout << " p0 at sample (x1=1, x2=1): " << p0_sample << "\n"; + std::cout << " p0 at counter (x1=12, x2=16): " << p0_counter << "\n"; + std::cout << " Signs " << (p0_sample == p0_counter ? "match" : "DIFFER") << "\n"; + + // Check if counterexample is inside the cell at each level + // For soundness: if counter has different poly sign, it MUST be outside the cell + std::cout << "\n=== Checking if counterexample is in the cell ===\n"; + bool counter_outside = false; + + // Level 0: (-oo, +oo) - always inside + std::cout << "Level 0: (-oo, +oo) -> counter inside\n"; + + // Level 1: check if x1=12 is in the cell + { + auto const& interval = cell[1]; + if (!interval.l_inf()) { + polynomial_ref lower_p(interval.l, pm); + scoped_anum_vector lower_roots(am); + // Partial assignment for level 1 (just x0) + nlsat::assignment partial_as(am); + scoped_anum val0(am); + am.set(val0, 1); // counter x0 = 1 + partial_as.set(x0, val0); + am.isolate_roots(lower_p, nlsat::undef_var_assignment(partial_as, 1), lower_roots); + if (lower_roots.size() >= interval.l_index) { + scoped_anum counter_x1(am); + am.set(counter_x1, 12); + int cmp = am.compare(counter_x1, lower_roots[interval.l_index - 1]); + std::cout << "Level 1 lower bound: root[" << interval.l_index << "] of poly, counter x1=12 is " + << (cmp > 0 ? "ABOVE" : (cmp < 0 ? "BELOW" : "AT")) << " lower bound\n"; + if (cmp <= 0) counter_outside = true; + } + } + if (!interval.u_inf()) { + polynomial_ref upper_p(interval.u, pm); + scoped_anum_vector upper_roots(am); + // Partial assignment for level 1 (just x0) + nlsat::assignment partial_as(am); + scoped_anum val0(am); + am.set(val0, 1); // counter x0 = 1 + partial_as.set(x0, val0); + am.isolate_roots(upper_p, nlsat::undef_var_assignment(partial_as, 1), upper_roots); + if (upper_roots.size() >= interval.u_index) { + scoped_anum counter_x1(am); + am.set(counter_x1, 12); + int cmp = am.compare(counter_x1, upper_roots[interval.u_index - 1]); + std::cout << "Level 1 upper bound: root[" << interval.u_index << "] of poly, counter x1=12 is " + << (cmp > 0 ? "ABOVE" : (cmp < 0 ? "BELOW" : "AT")) << " upper bound\n"; + if (cmp >= 0) counter_outside = true; + } + } + } + + // For a sound cell, if polynomial signs differ, counter MUST be outside the cell + // The fix (projecting p0's discriminant) should create bounds that exclude the counterexample + if (p0_sample != p0_counter) { + std::cout << "\nPoly signs differ between sample and counter.\n"; + std::cout << "For cell to be sound, counter must be OUTSIDE the cell.\n"; + std::cout << "Counter is " << (counter_outside ? "OUTSIDE" : "INSIDE") << " the cell.\n"; + ENSURE(counter_outside); // Cell is sound if counter is outside + } else { + std::cout << "\nPoly signs match - cell is trivially sound.\n"; + } + + std::cout << "\n=== END tst_unsound_lws_disc_zero ===\n\n"; +} + +// Test case for unsound lemma from ppblockterm.t2__p7867_terminationG_0.smt2 +// Issue z3-76w: levelwise produces unsound cell +static void tst_unsound_lws_ppblockterm() { + std::cout << "=== tst_unsound_lws_ppblockterm ===\n"; + params_ref ps; + ps.set_bool("lws", true); + reslimit rlim; + nlsat::solver s(rlim, ps, false); + anum_manager & am = s.am(); + nlsat::pmanager & pm = s.pm(); + polynomial::cache cache(pm); + + // Create 25 variables x0 to x24 + nlsat::var vars[25]; + polynomial_ref xvars[25] = { + polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), + polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), + polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), + polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), + polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm), polynomial_ref(pm) + }; + + for (unsigned i = 0; i < 25; i++) { + vars[i] = s.mk_var(false); + xvars[i] = pm.mk_polynomial(vars[i]); + } + + // Aliases for readability + polynomial_ref &x2 = xvars[2], &x4 = xvars[4], &x5 = xvars[5]; + polynomial_ref &x9 = xvars[9], &x10 = xvars[10], &x11 = xvars[11]; + polynomial_ref &x12 = xvars[12], &x13 = xvars[13], &x15 = xvars[15]; + polynomial_ref &x16 = xvars[16], &x19 = xvars[19], &x24 = xvars[24]; + + // p[0]: x16 + x9 + polynomial_ref p0(pm); + p0 = x16 + x9; + + // p[1]: x16 x24 + x9 x24 + x13 x19 + x5 x19 + x12 x15 + x10 x15 - 2 + polynomial_ref p1(pm); + p1 = x16 * x24 + x9 * x24 + x13 * x19 + x5 * x19 + x12 * x15 + x10 * x15 - 2; + + // p[2]: x4 + polynomial_ref p2(pm); + p2 = x4; + + // p[3]: x4 x24 + x2 x19 + x11 x15 + polynomial_ref p3(pm); + p3 = x4 * x24 + x2 * x19 + x11 * x15; + + std::cout << "Input polynomials:\n"; + std::cout << " p0: " << p0 << "\n"; + std::cout << " p1: " << p1 << "\n"; + std::cout << " p2: " << p2 << "\n"; + std::cout << " p3: " << p3 << "\n\n"; + + // Set sample point values + // x0->1, x1->1, x2->-1, x3->-1, x4->1, x5->-1, x6->1, x7->2, x8->1 + // x9->1, x10->-1, x11->1, x12->-2, x13->-2, x14->-2, x15->0, x16->2 + // x17->0, x18->0, x19->0, x20->0, x21->0, x22->0, x23->0 + int sample_vals[24] = {1, 1, -1, -1, 1, -1, 1, 2, 1, 1, -1, 1, -2, -2, -2, 0, 2, 0, 0, 0, 0, 0, 0, 0}; + + scoped_anum val(am); + nlsat::assignment sample_as(am); + for (unsigned i = 0; i < 24; i++) { + am.set(val, sample_vals[i]); + sample_as.set(vars[i], val); + } + s.set_rvalues(sample_as); + + std::cout << "Sample point (x0..x23):\n"; + for (unsigned i = 0; i < 24; i++) { + std::cout << " x" << i << " -> " << sample_vals[i] << "\n"; + } + std::cout << "\n"; + + // Evaluate polynomials at sample (with x24 = 0 for checking) + am.set(val, 0); + sample_as.set(vars[24], val); + std::cout << "Polynomial signs at sample (x24=0):\n"; + std::cout << " p0 sign: " << am.eval_sign_at(p0, sample_as) << "\n"; + std::cout << " p1 sign: " << am.eval_sign_at(p1, sample_as) << "\n"; + std::cout << " p2 sign: " << am.eval_sign_at(p2, sample_as) << "\n"; + std::cout << " p3 sign: " << am.eval_sign_at(p3, sample_as) << "\n\n"; + + // Counterexample from the unsound lemma: + // x2=-1, x4=1, x5=0, x9=0, x10=0, x11=-1, x12=-1, x13=-2, x15=3, x16=2, x19=0, x24=3 + nlsat::assignment counter_as(am); + am.set(val, -1); counter_as.set(vars[2], val); + am.set(val, 1); counter_as.set(vars[4], val); + am.set(val, 0); counter_as.set(vars[5], val); + am.set(val, 0); counter_as.set(vars[9], val); + am.set(val, 0); counter_as.set(vars[10], val); + am.set(val, -1); counter_as.set(vars[11], val); + am.set(val, -1); counter_as.set(vars[12], val); + am.set(val, -2); counter_as.set(vars[13], val); + am.set(val, 3); counter_as.set(vars[15], val); + am.set(val, 2); counter_as.set(vars[16], val); + am.set(val, 0); counter_as.set(vars[19], val); + am.set(val, 3); counter_as.set(vars[24], val); + + std::cout << "Counterexample point:\n"; + std::cout << " x2=-1, x4=1, x5=0, x9=0, x10=0, x11=-1, x12=-1, x13=-2, x15=3, x16=2, x19=0, x24=3\n\n"; + + std::cout << "Polynomial signs at counterexample:\n"; + std::cout << " p0 sign: " << am.eval_sign_at(p0, counter_as) << "\n"; + std::cout << " p1 sign: " << am.eval_sign_at(p1, counter_as) << "\n"; + std::cout << " p2 sign: " << am.eval_sign_at(p2, counter_as) << "\n"; + std::cout << " p3 sign: " << am.eval_sign_at(p3, counter_as) << "\n\n"; + + // Build polynomial vector + polynomial_ref_vector polys(pm); + polys.push_back(p0); + polys.push_back(p1); + polys.push_back(p2); + polys.push_back(p3); + + unsigned max_x = vars[24]; // x24 is the highest variable + + std::cout << "Running levelwise with max_x = x24...\n"; + nlsat::levelwise lws(s, polys, max_x, s.sample(), pm, am, cache); + auto cell = lws.single_cell(); + + if (lws.failed()) { + std::cout << "Levelwise FAILED\n"; + } else { + std::cout << "Levelwise succeeded\n"; + std::cout << "--- LEMMA from levelwise ---\n"; + for (unsigned i = 0; i < cell.size(); i++) { + auto const& interval = cell[i]; + std::cout << "Level x" << i << ": "; + if (interval.is_section()) { + std::cout << "section at root[" << interval.l_index << "] of " << interval.l << "\n"; + } else { + if (interval.l_inf()) + std::cout << "(-oo, "; + else + std::cout << "(root[" << interval.l_index << "] of " << interval.l << ", "; + if (interval.u_inf()) + std::cout << "+oo)"; + else + std::cout << "root[" << interval.u_index << "] of " << interval.u << ")"; + std::cout << "\n"; + } + } + std::cout << "--- END LEMMA ---\n\n"; + + // Check polynomial signs at sample and counterexample + int p0_sample = am.eval_sign_at(p0, sample_as); + int p1_sample = am.eval_sign_at(p1, sample_as); + int p2_sample = am.eval_sign_at(p2, sample_as); + int p3_sample = am.eval_sign_at(p3, sample_as); + + int p0_counter = am.eval_sign_at(p0, counter_as); + int p1_counter = am.eval_sign_at(p1, counter_as); + int p2_counter = am.eval_sign_at(p2, counter_as); + int p3_counter = am.eval_sign_at(p3, counter_as); + + bool signs_differ = (p0_sample != p0_counter) || (p1_sample != p1_counter) || + (p2_sample != p2_counter) || (p3_sample != p3_counter); + + if (signs_differ) { + std::cout << "Polynomial signs DIFFER between sample and counterexample.\n"; + } else { + std::cout << "Polynomial signs match between sample and counterexample.\n"; + } + + // Verify that the counterexample is OUTSIDE the cell (cell is sound) + std::cout << "\nChecking if counterexample is inside cell:\n"; + bool inside = is_point_inside_cell(am, pm, cell, counter_as); + + // For a sound cell with differing signs, counterexample must be OUTSIDE + ENSURE(!inside); + std::cout << "SUCCESS: Counterexample is OUTSIDE the cell (cell is sound)\n"; + } + + std::cout << "=== END tst_unsound_lws_ppblockterm ===\n\n"; +} + +// Test case for gh-8397: unsound lemma with lws=false +// Unsound lemma detected by nlsat.check_lemmas=true: +// !(1024 x1 = 0) or !(x1 + 1 > 0) or !(2048 x1^2 + 4096 x1 = 0) or +// !(x1 = root[3](1024 x1^3 + 6144 x1^2 + 6144 x1)) or +// !(x1 = root[3](2048 x1^3 + 6144 x1^2 + 4096 x1)) or !(x1 = 0) or +// 4 x6^3 + 4 x5 x6^2 - 4 x1 x6^2 + 4 x6^2 - 4 x1 x5 x6 - 4 x1 x6 < 0 or +// x6 - x1 = 0 or x6 > 0 +// Counterexample: x1 = 0, x5 = 0, x6 = -0.5 +static void tst_unsound_gh8397() { + // Reproduce exact unsound lemma from gh-8397 + // Unsound lemma: !(1024 x1 = 0) or !(x1 + 1 > 0) or !(2048 x1^2 + 4096 x1 = 0) or + // !(x1 = root[3](1024 x1^3 + 6144 x1^2 + 6144 x1)) or + // !(x1 = root[3](2048 x1^3 + 6144 x1^2 + 4096 x1)) or !(x1 = 0) or + // 4 x6^3 + 4 x5 x6^2 - 4 x1 x6^2 + 4 x6^2 - 4 x1 x5 x6 - 4 x1 x6 < 0 or x6 - x1 = 0 or x6 > 0 + // Counterexample: x1=0, x5=0, x6=-0.5 makes ALL literals FALSE + // Sample point: x0=-1, x1=0, x2=0, x3=-1, x4=0, x5=-1 (x6 is conflict var) + std::cout << "=== tst_unsound_gh8397 ===\n"; + + auto run_test = [](bool use_lws) { + std::cout << "\n--- Running with lws=" << (use_lws ? "true" : "false") << " ---\n"; + + params_ref ps; + ps.set_bool("lws", use_lws); + reslimit rlim; + nlsat::solver s(rlim, ps, false); + anum_manager& am = s.am(); + nlsat::pmanager& pm = s.pm(); + nlsat::explain& ex = s.get_explain(); + + // Create variables in order: x0, x1, x2, x3, x4, x5, x6 + nlsat::var x0 = s.mk_var(false); + nlsat::var x1 = s.mk_var(false); + nlsat::var x2 = s.mk_var(false); + nlsat::var x3 = s.mk_var(false); + nlsat::var x4 = s.mk_var(false); + nlsat::var x5 = s.mk_var(false); + nlsat::var x6 = s.mk_var(false); + (void)x0; (void)x2; (void)x3; (void)x4; // suppress unused warnings + + polynomial_ref _x1(pm), _x5(pm), _x6(pm); + _x1 = pm.mk_polynomial(x1); + _x5 = pm.mk_polynomial(x5); + _x6 = pm.mk_polynomial(x6); + + // Polynomials from the unsound lemma + // p_main: 4 x6^3 + 4 x5 x6^2 - 4 x1 x6^2 + 4 x6^2 - 4 x1 x5 x6 - 4 x1 x6 + polynomial_ref p_main(pm); + p_main = 4*_x6*_x6*_x6 + 4*_x5*_x6*_x6 - 4*_x1*_x6*_x6 + 4*_x6*_x6 - 4*_x1*_x5*_x6 - 4*_x1*_x6; + + // p_diff: x6 - x1 + polynomial_ref p_diff(pm); + p_diff = _x6 - _x1; + + // Additional polynomials for x1 constraints (from root literals) + // 1024 x1^3 + 6144 x1^2 + 6144 x1 = 1024*x1*(x1^2 + 6*x1 + 6) + polynomial_ref p_root1(pm); + p_root1 = 1024*_x1*_x1*_x1 + 6144*_x1*_x1 + 6144*_x1; + + // 2048 x1^3 + 6144 x1^2 + 4096 x1 = 1024*x1*(2*x1^2 + 6*x1 + 4) + polynomial_ref p_root2(pm); + p_root2 = 2048*_x1*_x1*_x1 + 6144*_x1*_x1 + 4096*_x1; + + // 1024 x1 + polynomial_ref p_x1_eq(pm); + p_x1_eq = 1024*_x1; + + // x1 + 1 + polynomial_ref p_x1_gt(pm); + p_x1_gt = _x1 + 1; + + // 2048 x1^2 + 4096 x1 + polynomial_ref p_x1_eq2(pm); + p_x1_eq2 = 2048*_x1*_x1 + 4096*_x1; + + std::cout << "p_main: " << p_main << "\n"; + std::cout << "p_diff: " << p_diff << "\n"; + std::cout << "p_root1: " << p_root1 << "\n"; + std::cout << "p_root2: " << p_root2 << "\n"; + + // Set sample point: x0=-1, x1=0, x2=0, x3=-1, x4=0, x5=-1 + nlsat::assignment sample_as(am); + scoped_anum val(am); + am.set(val, -1); sample_as.set(x0, val); + am.set(val, 0); sample_as.set(x1, val); + am.set(val, 0); sample_as.set(x2, val); + am.set(val, -1); sample_as.set(x3, val); + am.set(val, 0); sample_as.set(x4, val); + am.set(val, -1); sample_as.set(x5, val); + // x6 is the conflict variable, set to 0 for sample + am.set(val, 0); sample_as.set(x6, val); + s.set_rvalues(sample_as); + + // Input literals for the conflict (constraints on x6): + // The conflict is that x6 > x1 (i.e. x6 > 0) AND p_main < 0 + // But at sample point x1=0, x5=-1, these are infeasible at x6=0 + nlsat::scoped_literal_vector lits(s); + lits.push_back(mk_gt(s, p_diff.get())); // x6 - x1 > 0 (i.e. x6 > 0 since x1=0) + lits.push_back(mk_lt(s, p_main.get())); // p_main < 0 + + // Also add x1 constraints that were part of the conflict core + lits.push_back(mk_eq(s, p_x1_eq.get())); // 1024 x1 = 0 + lits.push_back(mk_gt(s, p_x1_gt.get())); // x1 + 1 > 0 + lits.push_back(mk_eq(s, p_x1_eq2.get())); // 2048 x1^2 + 4096 x1 = 0 + // Root literals: x1 = root[3](p_root1), x1 = root[3](p_root2) + lits.push_back(mk_root_eq(s, p_root1.get(), x1, 3)); // x1 = root[3](p_root1) + lits.push_back(mk_root_eq(s, p_root2.get(), x1, 3)); // x1 = root[3](p_root2) + lits.push_back(mk_eq(s, _x1.get())); // x1 = 0 + + std::cout << "Input literals:\n"; + s.display(std::cout, lits.size(), lits.data()); + std::cout << "\n"; + + // Call compute_conflict_explanation + nlsat::scoped_literal_vector result(s); + ex.compute_conflict_explanation(lits.size(), lits.data(), result); + + std::cout << "Result lemma (lws=" << (use_lws ? "true" : "false") << "):\n"; + std::cout << "(or"; + for (auto l : result) { + s.display(std::cout << "\n ", l); + } + for (unsigned i = 0; i < lits.size(); ++i) { + s.display(std::cout << "\n ", ~lits[i]); + } + std::cout << "\n)\n"; + + // Counterexample from checker: x1=0, x5=0, x6=-0.5 + nlsat::assignment counter_as(am); + am.set(val, -1); counter_as.set(x0, val); + am.set(val, 0); counter_as.set(x1, val); + am.set(val, 0); counter_as.set(x2, val); + am.set(val, -1); counter_as.set(x3, val); + am.set(val, 0); counter_as.set(x4, val); + am.set(val, 0); counter_as.set(x5, val); // changed from -1 to 0 + mpq half; am.qm().set(half, -1, 2); // -1/2 + am.set(val, half); counter_as.set(x6, val); + + std::cout << "\nEvaluating at counterexample (x1=0, x5=0, x6=-0.5):\n"; + + // Evaluate each result literal at counterexample + small_object_allocator allocator; + nlsat::evaluator eval(s, counter_as, pm, allocator); + + bool all_false = true; + for (auto l : result) { + nlsat::atom* a = s.bool_var2atom(l.var()); + if (a) { + bool value = eval.eval(a, l.sign()); + s.display(std::cout << " ", l); + std::cout << " -> " << (value ? "TRUE" : "FALSE") << "\n"; + if (value) all_false = false; + } + } + + // Also check negated input literals + for (unsigned i = 0; i < lits.size(); ++i) { + nlsat::literal nl = ~lits[i]; + nlsat::atom* a = s.bool_var2atom(nl.var()); + if (a) { + bool value = eval.eval(a, nl.sign()); + s.display(std::cout << " ", nl); + std::cout << " -> " << (value ? "TRUE" : "FALSE") << "\n"; + if (value) all_false = false; + } + } + + if (all_false) { + std::cout << "*** ALL literals FALSE at counterexample - LEMMA IS UNSOUND! ***\n"; + } else { + std::cout << "At least one literal is TRUE - lemma is sound at this point\n"; + } + }; + + run_test(false); // lws=false (buggy) + run_test(true); // lws=true (should be correct) + + std::cout << "\n=== END tst_unsound_gh8397 ===\n\n"; +} + void tst_nlsat() { + tst_unsound_gh8397(); + return; + std::cout << "------------------\n"; + tst_unsound_lws_ppblockterm(); + std::cout << "------------------\n"; + tst_unsound_lws_n46(); + std::cout << "------------------\n"; + tst_unsound_lws_et4(); + std::cout << "------------------\n"; + tst_unsound_lws_x3(); + std::cout << "------------------\n"; + tst_unsound_lws2380(); + std::cout << "------------------\n"; + tst_polynomial_cache_mk_unique(); + std::cout << "------------------\n"; + tst_nullified_polynomial(); std::cout << "------------------\n"; tst11(); std::cout << "------------------\n"; - return; tst10(); std::cout << "------------------\n"; tst9(); diff --git a/src/util/trace_tags.def b/src/util/trace_tags.def index 7d8c0928a..8eefa1d04 100644 --- a/src/util/trace_tags.def +++ b/src/util/trace_tags.def @@ -551,6 +551,8 @@ X(Global, linear_equation_mk, "linear equation mk") X(Global, list, "list") X(Global, literal_occ, "literal occ") X(Global, lp_core, "lp core") +X(Global, lws, "levelwise") +X(Global, lws_norm, "levelwise normalize") X(Global, macro_bug, "macro bug") X(Global, macro_finder, "macro finder") X(Global, macro_insert, "macro insert") diff --git a/src/util/util.h b/src/util/util.h index 0056c132f..84584f187 100644 --- a/src/util/util.h +++ b/src/util/util.h @@ -345,6 +345,7 @@ public: } void set_seed(unsigned s) { m_data = s; } + unsigned get_seed() const { return m_data; } int operator()() { return ((m_data = m_data * 214013L + 2531011L) >> 16) & 0x7fff;