diff --git a/src/math/lp/nla_stellensatz.cpp b/src/math/lp/nla_stellensatz.cpp index 44728f35c..8fa128a4e 100644 --- a/src/math/lp/nla_stellensatz.cpp +++ b/src/math/lp/nla_stellensatz.cpp @@ -27,19 +27,44 @@ necessarily derivable with the old constraints. Strategy 3: The lemma uses the new constraints. - --*/ + +Add also factoring: + +(16) - 16384*j6 - j8 + j10 >= 0 +(65) j32 - j36 >= 1 +(18) j8 >= 0 +(66) j34 - 16384*j36 - j38 <= 0 +(52) - j18 + j30 <= 0 +(72) j38 <= 16383 +(58) j4 - j12 - j32 >= 0 +(139) - 32769*j4 + j4*j18 <= -1 +(12) j4 - j6 <= 0 +(108) 2*j10 - j10*j12 <= 0 +(60) j10 - j30 - j34 <= 0 +(111) j4 >= 6 + ==> false + + Note that: +j10 > 0 => 2 - j12 <= 0 +j10 < 0 => 2 - j12 >= 0 +Then add assumption j10 > 0 or j10 < 0 + +p <= 0, factor into p*q <= 0, then assume p <= 0, q >= 0 + + */ #include "math/lp/nla_core.h" #include "math/lp/nla_coi.h" #include "math/lp/nla_stellensatz.h" - namespace nla { stellensatz::stellensatz(core* core) : common(core), m_solver(*this), - m_coi(*core) {} + m_coi(*core), + m_pm(c().reslim(), m_qm, &m_allocator) + {} lbool stellensatz::saturate() { lp::explanation ex; @@ -62,6 +87,7 @@ namespace nla { m_old_constraints.reset(); m_new_bounds.reset(); m_to_refine.reset(); + m_factored_constraints.reset(); m_coi.init(); init_vars(); } @@ -427,15 +453,13 @@ namespace nla { continue; svector _vars; _vars.push_back(v); - auto cs = var2cs[v]; - for (auto ci : cs) { + for (auto ci : var2cs[v]) { for (auto [coeff, u] : m_solver.lra().constraints()[ci].coeffs()) { if (u == v) saturate_constraint(ci, j, _vars); } } } - #if 0 for (auto v : vars) { if (v >= vars2cs.size()) continue; @@ -445,7 +469,6 @@ namespace nla { saturate_constraint(ci, j, m_mon2vars[u]); } } - #endif } IF_VERBOSE(5, c().lra_solver().display(verbose_stream() << "original\n"); @@ -519,6 +542,85 @@ namespace nla { insert_monomials_from_constraint(new_ci); m_ci2dep.setx(new_ci, nullptr, nullptr); m_old_constraints.insert(new_ci, old_ci); + } + + // call polynomial factorization using the polynomial manager + // return true if there is a non-trivial factorization + // need the following: + // term -> polynomial translation + // polynomial -> term translation + // the call to factorization + bool stellensatz::get_factors(term_offset& t, vector& factors) { + auto p = to_poly(t); + polynomial::factors fs(m_pm); + m_pm.factor(p, fs); + if (fs.distinct_factors() <= 1) + return false; + unsigned df = fs.distinct_factors(); + for (unsigned i = 0; i < df; ++i) { + auto f = fs[i]; + auto degree = fs.get_degree(i); + term_offset tf = to_term(*f); + for (unsigned j = 0; j < degree; ++j) + factors.push_back(tf); + } + return true; + } + + polynomial::polynomial_ref stellensatz::to_poly(term_offset const& t) { + + ptr_vector ms; + vector coeffs; + rational den(1); + for (lp::lar_term::ival kv : t.first) { + // TODO add to ms + den = lcm(den, denominator(kv.coeff())); + } + + for (auto kv : t.first) + coeffs.push_back(den * kv.coeff()); + coeffs.push_back(den * t.second); + + polynomial::polynomial_ref p(m_pm); + p = m_pm.mk_polynomial(ms.size(), coeffs.data(), ms.data()); + return p; + } + + stellensatz::term_offset stellensatz::to_term(polynomial::polynomial const& p) { + throw nullptr; + } + + void stellensatz::saturate_factors(lp::constraint_index ci) { + if (m_factored_constraints.contains(ci)) + return; + m_factored_constraints.insert(ci); + auto const &con = m_solver.lra().constraints()[ci]; + if (all_of(con.coeffs(), [&](auto const& cv) { return !m_mon2vars.contains(cv.second); })) + return; + term_offset t; + // ignore nested terms and nested polynomials.. + for (auto const& [coeff, v] : con.coeffs()) + t.first.add_monomial(coeff, v); + t.second = -con.rhs(); + + vector factors; + if (get_factors(t, factors)) + return; + + vector values; + for (auto const &t : factors) + values.push_back(value(t.first) + t.second); + + // p * q >= 0 => (p >= 0 & q >= 0) or (p <= 0 & q <= 0) + // assert both factors with bound assumptions, and reference to original constraint + // p >= 0 <- q >= 0 and + // q >= 0 <- p >= 0 + // the values of p, q in the current interpretation may not satisfy p*q >= 0 + // therefore, assume all but one polynomial satisfies the bound (ignoring multiplicities in this argument). + // Then derive the bound on the remaining polynomial from the others. + + + } // @@ -791,15 +893,14 @@ namespace nla { for (auto v : s.m_coi.vars()) s.c().lra_solver().set_column_value(v, lp::impq(s.m_values[v], rational(0))); } - #if 1 for (auto j : s.m_to_refine) { auto val_j = values[j]; auto const &vars = s.m_mon2vars[j]; auto val_vars = s.m_values[j]; + TRACE(arith, tout << "refine j" << j << " " << vars << "\n"); if (val_j != val_vars) s.saturate_basic_linearize(j, val_j, vars, val_vars); } - #endif return satisfies_products; } } diff --git a/src/math/lp/nla_stellensatz.h b/src/math/lp/nla_stellensatz.h index b36da638d..6a5208adc 100644 --- a/src/math/lp/nla_stellensatz.h +++ b/src/math/lp/nla_stellensatz.h @@ -6,6 +6,7 @@ #include "math/lp/nla_coi.h" #include "math/lp/int_solver.h" +#include "math/polynomial/polynomial.h" #include namespace nla { @@ -55,6 +56,10 @@ namespace nla { map, eq> m_vars2mon; u_map m_mon2vars; + small_object_allocator m_allocator; + unsynch_mpq_manager m_qm; + polynomial::manager m_pm; + struct resolvent { lp::constraint_index old_ci; lpvar mi; @@ -101,9 +106,15 @@ namespace nla { void saturate_monotonicity(lpvar j, rational const &val_j, svector const &vars, rational const &val_vars); void saturate_monotonicity(lpvar j, rational const & val_j, lpvar x, int sign_x, lpvar y, int sign_y); void saturate_monotonicity(lpvar j, rational const &val_j, lpvar x, lpvar y); - void saturate_squares(lpvar j, rational const &val_j, svector const &vars); + // polynomial tricks + using term_offset = std::pair; // term and its offset + uint_set m_factored_constraints; + bool get_factors(term_offset &t, vector &factors); + polynomial::polynomial_ref to_poly(term_offset const &t); + term_offset to_term(polynomial::polynomial const &p); + void saturate_factors(lp::constraint_index ci); // lemmas void add_lemma(lp::explanation const& ex); diff --git a/src/math/lp/nra_solver.cpp b/src/math/lp/nra_solver.cpp index 3a7f40093..0d75eb1e1 100644 --- a/src/math/lp/nra_solver.cpp +++ b/src/math/lp/nra_solver.cpp @@ -120,7 +120,7 @@ struct solver::imp { } } m_nlsat->collect_statistics(st); - TRACE(nra, + CTRACE(nra, false, m_nlsat->display(tout << r << "\n"); display(tout); for (auto [j, x] : m_lp2nl) tout << "j" << j << " := x" << x << "\n";); @@ -150,11 +150,11 @@ struct solver::imp { for (auto c : core) { unsigned idx = static_cast(static_cast(c) - this); ex.push_back(idx); - TRACE(nra, lra.display_constraint(tout << "ex: " << idx << ": ", idx) << "\n";); } nla::lemma_builder lemma(m_nla_core, __FUNCTION__); lemma &= ex; m_nla_core.set_use_nra_model(true); + TRACE(nra, tout << lemma << "\n"); break; } case l_undef: @@ -342,6 +342,7 @@ struct solver::imp { ex.push_back(ci); nla::lemma_builder lemma(m_nla_core, __FUNCTION__); lemma &= ex; + TRACE(nra, tout << lemma << "\n"); break; } case l_undef: