diff --git a/src/nlsat/levelwise.cpp b/src/nlsat/levelwise.cpp index c890b514c..4f006123c 100644 --- a/src/nlsat/levelwise.cpp +++ b/src/nlsat/levelwise.cpp @@ -105,9 +105,9 @@ namespace nlsat { } struct root_function { - scoped_anum val; + scoped_anum val; indexed_root_expr ire; - unsigned ps_idx; // index in m_level_ps + 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; } @@ -345,10 +345,10 @@ namespace nlsat { void request_factorized(polynomial_ref const& poly, inv_req req) { for_each_distinct_factor(poly, [&](polynomial_ref const& f) { TRACE(lws, - tout << " request_factorized: factor="; - m_pm.display(tout, f); - tout << " at level " << m_pm.max_var(f) << "\n"; - ); + tout << " request_factorized: factor="; + m_pm.display(tout, f); + tout << " at level " << m_pm.max_var(f) << "\n"; + ); request(f.get(), req); // inherit tag across factorization (SMT-RAT appendOnCorrectLevel) }); } @@ -421,12 +421,12 @@ namespace nlsat { return best; } - void add_projections_for(polynomial_ref const& p, unsigned x, polynomial_ref const& nonzero_coeff, bool add_leading_coeff, bool add_discriminant) { + 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, - tout << " add_projections_for: p="; - m_pm.display(tout, p); - tout << " x=" << x << " add_lc=" << add_leading_coeff << " add_disc=" << add_discriminant << "\n"; - ); + tout << " add_projection_for_poly: p="; + m_pm.display(tout, p); + tout << " 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); @@ -439,10 +439,10 @@ namespace nlsat { polynomial_ref lc(m_pm); lc = m_pm.coeff(p, x, deg); TRACE(lws, - tout << " adding lc: "; - m_pm.display(tout, lc); - tout << "\n"; - ); + tout << " adding lc: "; + m_pm.display(tout, lc); + tout << "\n"; + ); request_factorized(lc, inv_req::sign); } } @@ -455,14 +455,15 @@ namespace nlsat { // Compute side_mask: track which side(s) each polynomial appears on // bit0 = lower (<= sample), bit1 = upper (> sample), 3 = both sides void compute_side_mask() { - m_side_mask.clear(); + if (!is_set(m_l_rf)) + return; m_side_mask.resize(m_level_ps.size(), 0); - anum const& v = sample().value(m_level); - for (auto const& rf : m_rel.m_rfunc) { - if (m_am.compare(rf.val, v) <= 0) - m_side_mask[rf.ps_idx] |= 1; + 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[rf.ps_idx] |= 2; + m_side_mask[ps_idx] |= 2; } } @@ -545,7 +546,6 @@ namespace nlsat { if (m_rel.m_rfunc.empty() || m_rel.m_pairs.empty()) return; - compute_side_mask(); if (require_leaf) compute_resultant_graph_degree(); @@ -564,28 +564,22 @@ namespace nlsat { } } - // Compute noLdcf for sector case void compute_omit_lc_sector() { if (!is_set(m_l_rf) || !is_set(m_u_rf)) return; - if (m_rel_mode == spanning_tree) - compute_omit_lc_both_sides(true); - else - compute_omit_lc_both_sides(false); + // m_rel == spanning_tree, only leaves can drop ldcfs, but in the biggest cell all non-bounds are leaves, + compute_omit_lc_both_sides(m_rel_mode == spanning_tree); } - // Compute noLdcf for section case void compute_omit_lc_section() { - if (m_rel_mode == biggest_cell) - m_omit_lc.clear(); // no omit_lc for biggest_cell - else // spanning_tree + if (m_rel_mode == spanning_tree) compute_omit_lc_both_sides(true); } // Relation construction heuristics (same intent as previous implementation). - void fill_relation_with_biggest_cell_heuristic() { + 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"; - ); + 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";); @@ -604,8 +598,8 @@ namespace nlsat { 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"; - ); + tout << " fill_biggest_cell done: pairs.size=" << m_rel.m_pairs.size() << "\n"; + ); } void fill_relation_pairs_for_section_biggest_cell() { @@ -630,229 +624,146 @@ namespace nlsat { // 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 - bool build_both_side_spanning_tree() { + void build_both_side_spanning_tree() { auto const& rfs = m_rel.m_rfunc; unsigned n = rfs.size(); - if (n == 0) return false; + SASSERT(n > 0 && is_set(m_l_rf) && is_set(m_u_rf)); + SASSERT(!is_section()); + SASSERT(!same_boundary_poly()); - // Map ps_idx -> rfunc index on lower/upper side - std_vector lower_pos(m_level_ps.size(), UINT_MAX); - std_vector upper_pos(m_level_ps.size(), UINT_MAX); + // 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; - if (is_set(m_l_rf)) { - for (unsigned i = 0; i <= m_l_rf; ++i) - lower_pos[rfs[i].ps_idx] = i; - } - // For sector case: upper side is [m_u_rf, n) - // For section case: upper side is (m_l_rf, n) since m_u_rf is not set - unsigned upper_start = is_set(m_u_rf) ? m_u_rf : (is_set(m_l_rf) ? m_l_rf + 1 : 0); - for (unsigned i = upper_start; i < n; ++i) - upper_pos[rfs[i].ps_idx] = i; + // For sector: lower side is [0, m_l_rf], upper side is [m_u_rf, n) - // Collect both-side polynomial ps_idx's (must have valid positions on both sides) - std_vector both; + // 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 && lower_pos[i] != UINT_MAX && upper_pos[i] != UINT_MAX) - both.push_back(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]}); + } } - // Need at least threshold polynomials on both sides - if (both.size() < m_spanning_tree_threshold) - return false; + SASSERT(both.size() >= m_spanning_tree_threshold); - // Sort both by lower_pos (this is K in the lemma) + // Sort by lower_rf (root position on lower side) std::sort(both.begin(), both.end(), - [&](unsigned a, unsigned b) { return lower_pos[a] < lower_pos[b]; }); + [](both_info const& a, both_info const& b) { return a.lower_rf < b.lower_rf; }); - // Lemma construction: iteratively connect min(remaining) to element with min f-value in tree + // Build spanning tree using Reaching Orders Lemma: + // Start with element that has minimum upper_rf, root of out-arborescence on the right std_vector in_tree(both.size(), false); - - // Start with element that has minimum upper_pos (root of out-arborescence on f(K)) unsigned root_idx = 0; for (unsigned i = 1; i < both.size(); ++i) { - if (upper_pos[both[i]] < upper_pos[both[root_idx]]) + if (both[i].upper_rf < both[root_idx].upper_rf) root_idx = i; } in_tree[root_idx] = true; - // Add remaining elements + // Iteratively add remaining elements for (unsigned added = 1; added < both.size(); ++added) { - // Find minimum lower_pos element not yet in tree (element 'a' in lemma) + // Find min lower_rf element not in tree (first in sorted order) unsigned a_idx = UINT_MAX; - for (unsigned i = 0; i < both.size(); ++i) { - if (!in_tree[i]) { - a_idx = i; - break; // both is sorted by lower_pos, so first not-in-tree is min - } - } + for (unsigned i = 0; i < both.size(); ++i) + if (!in_tree[i]) { a_idx = i; break; } + SASSERT(a_idx != UINT_MAX); - // Find element in tree with minimum upper_pos (element 'c' in lemma) + // Find element in tree with min upper_rf unsigned c_idx = UINT_MAX; - unsigned min_upper = UINT_MAX; for (unsigned i = 0; i < both.size(); ++i) { - if (in_tree[i] && upper_pos[both[i]] < min_upper) { - min_upper = upper_pos[both[i]]; + if (in_tree[i] && (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} as ps_idx pair directly to m_rel.m_pairs - unsigned ps_a = both[a_idx]; - unsigned ps_c = both[c_idx]; + // Add edge {a, c} as ps_idx pair + 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)}); in_tree[a_idx] = true; } - return true; } // Sector spanning tree heuristic: // 1. Build spanning tree on both-side polynomials // 2. Extend with lowest_degree for single-side polynomials - void fill_relation_with_spanning_tree_heuristic() { + // 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(); - if (n == 0) return; - - // Need side_mask computed - compute_side_mask(); - - // Build spanning tree on both-side polynomials (adds to m_rel.m_pairs) - if (!build_both_side_spanning_tree()) { - // Not enough valid both-side polys, fall back - fill_relation_with_biggest_cell_heuristic(); - return; - } - - // Both-side polynomials exist, so both boundaries must be set - SASSERT(is_set(m_l_rf) && is_set(m_u_rf)); - - // Extend for single-side polynomials using biggest_cell style - // (connect all unconnected to the boundary) - - // Lower side: connect all unconnected polynomials to the boundary m_l_rf + + // 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) // single-side poly + if (m_side_mask[rfs[j].ps_idx] != 3) m_rel.add_pair(j, m_l_rf); } - // If boundary is not a both-side poly, connect spanning tree to boundary - // using the both-side poly with maximum lower position - if (m_side_mask[rfs[m_l_rf].ps_idx] != 3) { - unsigned max_lower_idx = UINT_MAX; - for (unsigned j = 0; j < m_l_rf; ++j) { - if (m_side_mask[rfs[j].ps_idx] == 3) - max_lower_idx = j; // last one found is maximum - } - if (max_lower_idx != UINT_MAX) - m_rel.add_pair(max_lower_idx, m_l_rf); - } - - // Upper side: connect all unconnected polynomials to the boundary m_u_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) // single-side poly + if (m_side_mask[rfs[j].ps_idx] != 3) m_rel.add_pair(m_u_rf, j); } - // If boundary is not a both-side poly, connect spanning tree to boundary - // using the both-side poly with minimum upper position - if (m_side_mask[rfs[m_u_rf].ps_idx] != 3) { - unsigned min_upper_idx = UINT_MAX; - for (unsigned j = m_u_rf + 1; j < n; ++j) { - if (m_side_mask[rfs[j].ps_idx] == 3) { - min_upper_idx = j; // first one found is minimum - break; - } - } - if (min_upper_idx != UINT_MAX) - m_rel.add_pair(m_u_rf, min_upper_idx); + } + + // 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); } - // Section spanning tree heuristic - void fill_relation_pairs_for_section_spanning_tree() { - 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); - - // Need side_mask computed - compute_side_mask(); - - // Build spanning tree on both-side polynomials (adds to m_rel.m_pairs) - if (!build_both_side_spanning_tree()) { - // Not enough valid both-side polys, fall back - fill_relation_pairs_for_section_biggest_cell(); - return; - } - - // Extend for single-side polynomials using biggest_cell style - // (connect all unconnected to the section m_l_rf) - - // Below section: connect all unconnected to section - for (unsigned j = 0; j < m_l_rf; ++j) { - if (m_side_mask[rfs[j].ps_idx] != 3) // single-side poly - m_rel.add_pair(j, m_l_rf); - } - // If section poly is not a both-side poly, connect spanning tree to section - // using the both-side poly with maximum lower position - if (m_side_mask[rfs[m_l_rf].ps_idx] != 3) { - unsigned max_lower_idx = UINT_MAX; - for (unsigned j = 0; j < m_l_rf; ++j) { - if (m_side_mask[rfs[j].ps_idx] == 3) - max_lower_idx = j; // last one found is maximum - } - if (max_lower_idx != UINT_MAX) - m_rel.add_pair(max_lower_idx, m_l_rf); - } - - // Above section: connect all unconnected to section - for (unsigned j = m_l_rf + 1; j < n; ++j) { - if (m_side_mask[rfs[j].ps_idx] != 3) // single-side poly - m_rel.add_pair(m_l_rf, j); - } - // If section poly is not a both-side poly, connect spanning tree to section - // using the both-side poly with minimum upper position - if (m_side_mask[rfs[m_l_rf].ps_idx] != 3) { - unsigned min_upper_idx = UINT_MAX; - for (unsigned j = m_l_rf + 1; j < n; ++j) { - if (m_side_mask[rfs[j].ps_idx] == 3) { - min_upper_idx = j; // first one found is minimum - break; - } - } - if (min_upper_idx != UINT_MAX) - m_rel.add_pair(m_l_rf, min_upper_idx); - } - } - - void fill_relation_pairs_for_section() { - SASSERT(m_I[m_level].section); - if (m_rel_mode == spanning_tree) - fill_relation_pairs_for_section_spanning_tree(); - else - fill_relation_pairs_for_section_biggest_cell(); - } - - void fill_relation_pairs_for_sector() { - SASSERT(!m_I[m_level].section); - if (m_rel_mode == spanning_tree) - fill_relation_with_spanning_tree_heuristic(); - else - fill_relation_with_biggest_cell_heuristic(); - } - // 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) { + anum const& v, std_vector& lhalf, + std_vector& uhalf) { scoped_anum_vector roots(m_am); @@ -910,14 +821,15 @@ namespace nlsat { } bool collect_partitioned_root_functions_around_sample(anum const& v, - std_vector& lhalf, std_vector& uhalf) { + 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) { + 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(); @@ -945,9 +857,9 @@ namespace nlsat { 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); }); + [&](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); }); + [&](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. @@ -981,7 +893,7 @@ namespace nlsat { 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 = UINT_MAX; + m_u_rf = m_l_rf; } else { m_I[m_level].l = r.ire.p; @@ -1018,35 +930,34 @@ namespace nlsat { 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"; - ); + 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, - tout << " resultant(" << pr.first << ", " << pr.second << "): "; - m_pm.display(tout, p1); - tout << " and "; - m_pm.display(tout, p2); - tout << "\n"; - ); + tout << " resultant(" << pr.first << ", " << pr.second << "): "; + m_pm.display(tout, p1); + tout << " and "; + m_pm.display(tout, p2); + tout << "\n"; + ); polynomial_ref res = psc_resultant(p1, p2, m_level); TRACE(lws, - tout << " resultant poly: "; - if (res) { - m_pm.display(tout, res); - tout << "\n resultant sign at sample: " << m_am.eval_sign_at(res, sample()); - } else { - tout << "(null)"; - } - tout << "\n"; - ); + tout << " resultant poly: "; + if (res) { + m_pm.display(tout, res); + tout << "\n resultant sign at sample: " << m_am.eval_sign_at(res, sample()); + } else { + tout << "(null)"; + } + tout << "\n"; + ); request_factorized(res, inv_req::ord); } } @@ -1063,6 +974,14 @@ namespace nlsat { 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]); @@ -1075,6 +994,18 @@ namespace nlsat { 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) { + tout << " [" << j << "] val="; + m_am.display_decimal(tout, root_vals[j].first, 5); + tout << " poly="; + m_pm.display(tout, root_vals[j].second); + tout << "\n"; + } + ); + std::set> added_pairs; for (unsigned j = 0; j + 1 < root_vals.size(); ++j) { poly* p1 = root_vals[j].second; @@ -1084,6 +1015,13 @@ namespace nlsat { if (p1 > p2) std::swap(p1, p2); if (!added_pairs.insert({p1, p2}).second) continue; + TRACE(lws, + tout << " Adjacent resultant pair: "; + m_pm.display(tout, p1); + tout << " and "; + m_pm.display(tout, p2); + tout << "\n"; + ); request_factorized(psc_resultant(p1, p2, m_n), inv_req::ord); } } @@ -1100,19 +1038,19 @@ namespace nlsat { void add_level_projections_sector() { TRACE(lws, - tout << "\n add_level_projections_sector at level " << m_level << "\n"; - tout << " Lower bound rf=" << m_l_rf << ", Upper bound rf=" << m_u_rf << "\n"; - if (is_set(m_l_rf)) { - tout << " lower poly idx=" << m_rel.m_rfunc[m_l_rf].ps_idx << ": "; - m_pm.display(tout, m_level_ps.get(m_rel.m_rfunc[m_l_rf].ps_idx)); - tout << "\n"; - } - if (is_set(m_u_rf)) { - tout << " upper poly idx=" << m_rel.m_rfunc[m_u_rf].ps_idx << ": "; - m_pm.display(tout, m_level_ps.get(m_rel.m_rfunc[m_u_rf].ps_idx)); - tout << "\n"; - } - ); + tout << "\n add_level_projections_sector at level " << m_level << "\n"; + tout << " Lower bound rf=" << m_l_rf << ", Upper bound rf=" << m_u_rf << "\n"; + if (is_set(m_l_rf)) { + tout << " lower poly idx=" << m_rel.m_rfunc[m_l_rf].ps_idx << ": "; + m_pm.display(tout, m_level_ps.get(m_rel.m_rfunc[m_l_rf].ps_idx)); + tout << "\n"; + } + if (is_set(m_u_rf)) { + tout << " upper poly idx=" << m_rel.m_rfunc[m_u_rf].ps_idx << ": "; + m_pm.display(tout, m_level_ps.get(m_rel.m_rfunc[m_u_rf].ps_idx)); + tout << "\n"; + } + ); // Lines 11-12 (Algorithm 1): add projections for each p // Note: Algorithm 1 adds disc + ldcf for ALL polynomials (classical delineability) // We additionally omit leading coefficients for rootless polynomials when possible @@ -1135,12 +1073,12 @@ namespace nlsat { bool is_upper_bound = is_set(m_u_rf) && i == m_rel.m_rfunc[m_u_rf].ps_idx; TRACE(lws, - tout << " poly[" << i << "] is_lower=" << is_lower_bound << " is_upper=" << is_upper_bound; - tout << " omit_lc[i]=" << (i < m_omit_lc.size() ? (m_omit_lc[i] ? "true" : "false") : "N/A"); - tout << " add_lc=" << add_lc << "\n"; - ); + tout << " poly[" << i << "] is_lower=" << is_lower_bound << " is_upper=" << is_upper_bound; + tout << " omit_lc[i]=" << (i < m_omit_lc.size() ? (m_omit_lc[i] ? "true" : "false") : "N/A"); + tout << " add_lc=" << add_lc << "\n"; + ); - if (add_lc && i < usize(m_poly_has_roots) && !m_poly_has_roots[i]) + if (add_lc && !poly_has_roots(i)) if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0) add_lc = false; @@ -1159,22 +1097,21 @@ namespace nlsat { if (add_disc && get_req(i) != inv_req::ord && i < m_omit_disc.size() && m_omit_disc[i]) add_disc = false; - add_projections_for(p, m_level, witness, add_lc, add_disc); + add_projection_for_poly(p, m_level, witness, add_lc, add_disc); } } + bool is_section() { return m_I[m_level].is_section();} + void add_level_projections_section() { - SASSERT(m_I[m_level].section); - poly* section_p = m_I[m_level].l.get(); - + SASSERT(is_section()); compute_omit_lc_section(); // m_omit_disc is computed by compute_omit_disc_for_spanning_tree() in process_level_section() - + poly* section_p = m_I[m_level].l.get(); + SASSERT(section_p); for (unsigned i = 0; i < m_level_ps.size(); ++i) { polynomial_ref p(m_level_ps.get(i), m_pm); - bool is_section_poly = section_p && p.get() == section_p; - bool has_roots = i < usize(m_poly_has_roots) && m_poly_has_roots[i]; - + bool is_section_poly = p.get() == section_p; polynomial_ref lc(m_pm); unsigned deg = m_pm.degree(p, m_level); lc = m_pm.coeff(p, m_level, deg); @@ -1184,7 +1121,7 @@ namespace nlsat { if (!is_section_poly && m_rel_mode == spanning_tree) { if (i < m_omit_lc.size() && m_omit_lc[i]) add_lc = false; - if (add_lc && !has_roots) + if (add_lc && !poly_has_roots(i)) if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0) add_lc = false; } @@ -1205,7 +1142,7 @@ namespace nlsat { if (lc && !is_zero(lc) && m_am.eval_sign_at(lc, sample()) != 0) witness = polynomial_ref(m_pm); - add_projections_for(p, m_level, witness, add_lc, add_disc); + add_projection_for_poly(p, m_level, witness, add_lc, add_disc); } } @@ -1214,15 +1151,22 @@ namespace nlsat { // 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() < 2) + if (m_rel.m_rfunc.size() < 4) // 2 different bounds + at least 2 same out of bounds return false; - compute_side_mask(); + 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 does not make sense: all same pairs will be obtained with the biggest cell + unsigned both_count = 0; - for (unsigned i = 0; i < m_level_ps.size(); ++i) { + 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; } @@ -1239,18 +1183,9 @@ namespace nlsat { void process_level_section(bool have_interval) { SASSERT(m_I[m_level].section); - clear_level_state(); // Clear stale state from previous level m_rel_mode = biggest_cell; // default if (have_interval) { - // Check spanning tree threshold first - if (should_use_spanning_tree()) { - fill_relation_pairs_for_section_spanning_tree(); - compute_omit_lc_both_sides(true); - compute_omit_disc_for_spanning_tree(); - m_rel_mode = spanning_tree; - } else { - fill_relation_pairs_for_section(); - } + fill_relation_pairs_for_section_biggest_cell(); SASSERT(relation_invariant()); } upgrade_bounds_to_ord(); @@ -1260,18 +1195,17 @@ namespace nlsat { void process_level_sector(bool have_interval) { SASSERT(!m_I[m_level].section); - clear_level_state(); // Clear stale state from previous level m_rel_mode = biggest_cell; // default if (have_interval) { // Check spanning tree threshold first if (should_use_spanning_tree()) { - fill_relation_with_spanning_tree_heuristic(); + fill_relation_sector_spanning_tree(); compute_omit_lc_both_sides(true); compute_omit_disc_for_spanning_tree(); m_rel_mode = spanning_tree; - } else { - fill_relation_pairs_for_sector(); - } + } else + fill_relation_sector_biggest_cell(); + SASSERT(relation_invariant()); } upgrade_bounds_to_ord(); @@ -1279,17 +1213,7 @@ namespace nlsat { add_relation_resultants(); } - void process_level() { - TRACE(lws, - tout << "\n--- process_level: level=" << m_level << " ---\n"; - tout << "Polynomials at this level (" << m_level_ps.size() << "):\n"; - for (unsigned i = 0; i < m_level_ps.size(); ++i) { - tout << " ps[" << i << "]: "; - m_pm.display(tout, m_level_ps.get(i)); - tout << "\n"; - } - ); - + 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()); @@ -1300,22 +1224,33 @@ namespace nlsat { 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, - tout << "Interval: "; - display(tout, m_solver, m_I[m_level]); - tout << "\n"; - tout << "Section? " << (m_I[m_level].section ? "yes" : "no") << "\n"; - tout << "have_interval=" << have_interval << ", rfunc.size=" << m_rel.m_rfunc.size() << "\n"; - for (unsigned i = 0; i < m_rel.m_rfunc.size(); ++i) { - tout << " rfunc[" << i << "]: ps_idx=" << m_rel.m_rfunc[i].ps_idx << ", val="; - m_am.display(tout, m_rel.m_rfunc[i].val); - tout << "\n"; - } - ); + tout << "Interval: "; + display(tout, m_solver, m_I[m_level]); + tout << "\n"; + tout << "Section? " << (m_I[m_level].section ? "yes" : "no") << "\n"; + tout << "have_interval=" << have_interval << ", rfunc.size=" << m_rel.m_rfunc.size() << "\n"; + for (unsigned i = 0; i < m_rel.m_rfunc.size(); ++i) { + tout << " rfunc[" << i << "]: ps_idx=" << m_rel.m_rfunc[i].ps_idx << ", val="; + m_am.display(tout, m_rel.m_rfunc[i].val); + tout << "\n"; + } + ); if (m_I[m_level].section) process_level_section(have_interval); @@ -1323,74 +1258,50 @@ namespace nlsat { process_level_sector(have_interval); } + bool poly_has_roots(unsigned i) { return i < usize(m_poly_has_roots) && m_poly_has_roots[i]; } + void process_top_level() { - TRACE(lws, - tout << "\n--- process_top_level: level=" << m_n << " ---\n"; - tout << "Polynomials at top level (" << m_level_ps.size() << "):\n"; - for (unsigned i = 0; i < m_level_ps.size(); ++i) { - tout << " ps[" << i << "]: "; - m_pm.display(tout, m_level_ps.get(i)); - tout << "\n"; - } - ); - - 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_n); - if (!w) - fail(); - m_witnesses.push_back(w); - } - - // Resultants between adjacent root functions (a lightweight ordering for the top level). + TRACE(lws, display_polys_at_level(tout);); + collect_non_null_witnesses(); add_adjacent_root_resultants(); // Projections (coeff witness, disc, leading coeff). - // Note: SMT-RAT's levelwise implementation additionally has dedicated heuristics for - // selecting resultants and selectively adding leading coefficients (see OneCellCAD.h, - // sectionHeuristic/sectorHeuristic). Z3's levelwise currently uses the relation_mode - // ordering heuristics instead of these specialized cases. for (unsigned i = 0; i < 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; // Let todo_set handle duplicates if witness == lc - if (i < usize(m_poly_has_roots) && !m_poly_has_roots[i]) { + 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; - } - // SMT-RAT-style coeffNonNull: if the leading coefficient is already non-zero at the sample + + // 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) && add_lc) + 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_projections_for(p, m_n, witness, add_lc, true); + 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) { - tout << " p[" << i << "]: "; - m_pm.display(tout, m_P.get(i)); - tout << "\n"; - } - tout << "Sample values:\n"; - for (unsigned j = 0; j < m_n; ++j) { - tout << " x" << j << " = "; - m_am.display(tout, sample().value(j)); - tout << "\n"; - } - ); + tout << "Input polynomials (" << m_P.size() << "):\n"; + for (unsigned i = 0; i < m_P.size(); ++i) { + tout << " p[" << i << "]: "; + m_pm.display(tout, m_P.get(i)) << "\n"; + } + tout << "Sample values:\n"; + for (unsigned j = 0; j < m_n; ++j) { + tout << " x" << j << " = "; + m_am.display(tout, sample().value(j)) << "\n"; + } + ); - if (m_n == 0) - return m_I; + if (m_n == 0) return m_I; m_todo.reset(); @@ -1402,8 +1313,7 @@ namespace nlsat { }); } - if (m_todo.empty()) - return m_I; + 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) { @@ -1417,19 +1327,19 @@ namespace nlsat { 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"; - ); + 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) { - tout << "I[" << i << "]: "; - display(tout, m_solver, m_I[i]); - tout << "\n"; - } - ); + for (unsigned i = 0; i < m_I.size(); ++i) { + tout << "I[" << i << "]: "; + display(tout, m_solver, m_I[i]); + tout << "\n"; + } + ); return m_I; } diff --git a/src/nlsat/nlsat_explain.cpp b/src/nlsat/nlsat_explain.cpp index 6c3a99d60..c7971c117 100644 --- a/src/nlsat/nlsat_explain.cpp +++ b/src/nlsat/nlsat_explain.cpp @@ -63,6 +63,9 @@ namespace nlsat { // 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; @@ -88,6 +91,7 @@ namespace nlsat { 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; @@ -1012,6 +1016,11 @@ namespace nlsat { */ 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()) @@ -1062,6 +1071,7 @@ namespace nlsat { 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); @@ -1849,6 +1859,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); } diff --git a/src/nlsat/nlsat_explain.h b/src/nlsat/nlsat_explain.h index 2fe967efb..3ff4fb982 100644 --- a/src/nlsat/nlsat_explain.h +++ b/src/nlsat/nlsat_explain.h @@ -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_params.pyg b/src/nlsat/nlsat_params.pyg index 7485905d0..881392d77 100644 --- a/src/nlsat/nlsat_params.pyg +++ b/src/nlsat/nlsat_params.pyg @@ -23,6 +23,6 @@ def_module_params('nlsat', ('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"), ('lws', BOOL, True, "apply levelwise."), - ('lws_spt_threshold', UINT, 3, "minimum both-side polynomial count to apply spanning tree optimization; < 2 disables spanning tree"), + ('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 97cf337fd..593514e4c 100644 --- a/src/nlsat/nlsat_solver.cpp +++ b/src/nlsat/nlsat_solver.cpp @@ -54,7 +54,6 @@ Revision History: namespace nlsat { - typedef chashtable ineq_atom_table; typedef chashtable root_atom_table; @@ -250,6 +249,7 @@ namespace nlsat { 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), @@ -1022,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); }; @@ -1058,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 { @@ -1085,12 +1088,20 @@ 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"; @@ -1130,7 +1141,7 @@ namespace nlsat { // Create a separate reslimit for the checker with 10 second timeout reslimit checker_rlimit; cancel_eh eh(checker_rlimit); - scoped_timer timer(1000, &eh); + scoped_timer timer(1000, &eh); // one second ctx c(checker_rlimit, m_ctx.m_params, m_ctx.m_incremental); solver solver2(c); @@ -1139,14 +1150,20 @@ 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)) + 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_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); } @@ -2125,9 +2142,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); @@ -2364,6 +2380,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)); } @@ -2609,9 +2627,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); @@ -4645,6 +4662,7 @@ namespace nlsat { 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 diff --git a/src/test/nlsat.cpp b/src/test/nlsat.cpp index 2684abb49..6d71468ed 100644 --- a/src/test/nlsat.cpp +++ b/src/test/nlsat.cpp @@ -1164,7 +1164,288 @@ static void tst_unsound_lws2380() { 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"; +} + void tst_nlsat() { + tst_unsound_lws_x3(); + std::cout << "------------------\n"; tst_unsound_lws2380(); std::cout << "------------------\n"; tst_polynomial_cache_mk_unique();