From e2235d81d345b01ef9c88c1fe9b85dba41e2a983 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Mon, 1 Sep 2025 16:37:21 -0700 Subject: [PATCH] add option for gcd-test to grobner --- src/math/dd/dd_pdd.cpp | 65 ++++++++++++++++++ src/math/dd/dd_pdd.h | 45 +++++++++++++ src/math/lp/nla_grobner.cpp | 111 ++++++++++++++++++++++++++----- src/math/lp/nla_grobner.h | 9 ++- src/math/lp/nla_pp.cpp | 25 ++++--- src/params/smt_params_helper.pyg | 1 + 6 files changed, 224 insertions(+), 32 deletions(-) diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp index 3ad64acfd..3c0eb9173 100644 --- a/src/math/dd/dd_pdd.cpp +++ b/src/math/dd/dd_pdd.cpp @@ -1208,6 +1208,35 @@ namespace dd { return is_binary(p.root); } + /* + * Retrieve variables that occur as powers. + */ + void pdd_manager::get_powers(pdd const& p, svector>& powers) { + powers.reset(); + init_mark(); + m_todo.push_back(p.root); + while (!m_todo.empty()) { + PDD r = m_todo.back(); + m_todo.pop_back(); + if (is_val(r) || is_marked(r)) + continue; + set_mark(r); + unsigned v = var(r); + unsigned d = 1; + while (!is_val(r) && var(hi(r)) == v) { + d++; + m_todo.push_back(lo(r)); + r = hi(r); + } + if (d > 1) + powers.push_back({ v, d }); + if (!is_val(r)) { + m_todo.push_back(hi(r)); + m_todo.push_back(lo(r)); + } + } + } + /** Determine if p is a monomial. */ @@ -1986,6 +2015,42 @@ namespace dd { pdd_iterator pdd::begin() const { return pdd_iterator(*this, true); } pdd_iterator pdd::end() const { return pdd_iterator(*this, false); } + pdd_coeff_iterator pdd::pdd_coeffients::begin() const { return m_pdd.begin_coeff(); } + pdd_coeff_iterator pdd::pdd_coeffients::end() const { return m_pdd.end_coeff(); } + + void pdd_coeff_iterator::next() { + auto& m = m_pdd.manager(); + while (!m_nodes.empty()) { + auto p = m_nodes.back(); + m_nodes.pop_back(); + SASSERT(!m.is_val(p)); + unsigned n = m.lo(p); + if (m.is_val(n) && m.val(n).is_zero()) + continue; + while (!m.is_val(n)) { + m_nodes.push_back(n); + n = m.hi(n); + } + m_value = { m.val(n), m_nodes.empty()}; + return; + } + at_end = true; + } + + void pdd_coeff_iterator::first() { + unsigned n = m_pdd.root; + auto& m = m_pdd.manager(); + while (!m.is_val(n)) { + m_nodes.push_back(n); + n = m.hi(n); + } + at_end = false; + m_value = { m.val(n), m_nodes.empty() }; + } + + pdd_coeff_iterator pdd::begin_coeff() const { return pdd_coeff_iterator(*this, true); } + pdd_coeff_iterator pdd::end_coeff() const { return pdd_coeff_iterator(*this, false); } + std::ostream& operator<<(std::ostream& out, pdd_monomial const& m) { if (!m.coeff.is_one()) { out << m.coeff; diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h index 4fe208429..b96e5f04b 100644 --- a/src/math/dd/dd_pdd.h +++ b/src/math/dd/dd_pdd.h @@ -46,6 +46,7 @@ namespace dd { class pdd_manager; class pdd_iterator; class pdd_linear_iterator; + class pdd_coeff_iterator; class pdd_manager { public: @@ -55,6 +56,7 @@ namespace dd { friend pdd; friend pdd_iterator; friend pdd_linear_iterator; + friend pdd_coeff_iterator; typedef unsigned PDD; typedef vector> monomials_t; @@ -366,6 +368,8 @@ namespace dd { bool is_binary(PDD p); bool is_binary(pdd const& p); + void get_powers(pdd const& p, svector>& powers); + bool is_monomial(PDD p); bool is_univariate(PDD p); @@ -406,6 +410,7 @@ namespace dd { friend class pdd_manager; friend class pdd_iterator; friend class pdd_linear_iterator; + friend class pdd_coeff_iterator; unsigned root; pdd_manager* m; pdd(unsigned root, pdd_manager& pm): root(root), m(&pm) { m->inc_ref(root); } @@ -440,6 +445,7 @@ namespace dd { bool is_unary() const { return !is_val() && lo().is_zero() && hi().is_val(); } bool is_offset() const { return !is_val() && lo().is_val() && hi().is_one(); } bool is_binary() const { return m->is_binary(root); } + void get_powers(svector>& powers) const { m->get_powers(*this, powers); } bool is_monomial() const { return m->is_monomial(root); } bool is_univariate() const { return m->is_univariate(root); } bool is_univariate_in(unsigned v) const { return m->is_univariate_in(root, v); } @@ -508,6 +514,20 @@ namespace dd { pdd_iterator begin() const; pdd_iterator end() const; + + pdd_coeff_iterator begin_coeff() const; + pdd_coeff_iterator end_coeff() const; + + class pdd_coeffients { + friend class pdd; + pdd const& m_pdd; + pdd_coeffients(pdd const& p) : m_pdd(p) {} + public: + pdd_coeff_iterator begin() const; + pdd_coeff_iterator end() const; + }; + pdd_coeffients coefficients() const { return pdd_coeffients(*this); } + class pdd_linear_monomials { friend class pdd; pdd const& m_pdd; @@ -593,6 +613,31 @@ namespace dd { bool operator!=(pdd_linear_iterator const& other) const { return m_next != other.m_next; } }; + // iterator class that visits each coefficient with information + // if it is a constant or is used as a coefficient to a monomial. + + struct coeff_value { + rational coeff; + bool is_constant; // true if it is a constant term. + }; + + class pdd_coeff_iterator { + friend class pdd; + pdd m_pdd; + unsigned_vector m_nodes; + coeff_value m_value; + bool at_end = true; + pdd_coeff_iterator(pdd const& p, bool at_start): m_pdd(p) { if (at_start) first(); } + void first(); + void next(); + public: + coeff_value operator*() const { return m_value; } + coeff_value const* operator->() const { return &m_value; } + pdd_coeff_iterator& operator++() { next(); return *this; } + pdd_coeff_iterator operator++(int) { auto tmp = *this; next(); return tmp; } + bool operator!=(pdd_coeff_iterator const& other) const { return at_end != other.at_end; } + }; + class val_pp { pdd_manager const& m; rational const& val; diff --git a/src/math/lp/nla_grobner.cpp b/src/math/lp/nla_grobner.cpp index 384d4a3c3..e039f8dc7 100644 --- a/src/math/lp/nla_grobner.cpp +++ b/src/math/lp/nla_grobner.cpp @@ -31,6 +31,7 @@ namespace nla { void grobner::updt_params(params_ref const& p) { smt_params_helper ph(p); m_config.m_propagate_quotients = ph.arith_nl_grobner_propagate_quotients(); + m_config.m_gcd_test = ph.arith_nl_grobner_gcd_test(); } lp::lp_settings& grobner::lp_settings() { @@ -81,6 +82,9 @@ namespace nla { if (propagate_quotients()) return; + + if (propagate_gcd_test()) + return; } catch (...) { @@ -99,14 +103,6 @@ namespace nla { IF_VERBOSE(5, diagnose_pdd_miss(verbose_stream())); } - dd::solver::equation_vector const& grobner::core_equations(bool all_eqs) { - flet _add_all(m_add_all_eqs, all_eqs); - find_nl_cluster(); - if (!configure()) - throw dd::pdd_manager::mem_out(); - return m_solver.equations(); - } - bool grobner::is_conflicting() { for (auto eq : m_solver.equations()) { if (is_conflicting(*eq)) { @@ -228,6 +224,76 @@ namespace nla { return {t, offset}; } + bool grobner::propagate_gcd_test() { + if (!m_config.m_gcd_test) + return false; + unsigned changed = 0; + for (auto eq : m_solver.equations()) + if (propagate_gcd_test(*eq) && ++changed >= m_solver.number_of_conflicts_to_report()) + return true; + return changed > 0; + } + + // x^2 - y = 0 => y mod 4 = { 0, 1 } + // x^k - y = 0 => y mod k^2 = {0, 1, .., (k-1)^k } + bool grobner::propagate_gcd_test(dd::solver::equation const& eq) { + dd::pdd p = eq.poly(); + if (p.is_linear()) + return false; + if (p.is_val()) + return false; + auto v = p.var(); + if (!c().var_is_int(v)) + return false; + for (auto v : p.free_vars()) + if (!c().var_is_int(v)) + return false; + dd::pdd_eval eval; + eval.var2val() = [&](unsigned j) { return val(j); }; + if (eval(p) == 0) + return false; + p.get_powers(m_powers); + if (m_powers.empty()) + return false; + rational d(1); + for (auto const& value : p.coefficients()) + d = lcm(d, denominator(value.coeff)); + if (d != 1) + p *= d; + auto& m = p.manager(); + auto fails_gcd_test = [&](dd::pdd const& q) { + rational g(0); + rational k(0); + for (auto const& value : q.coefficients()) { + if (value.is_constant) { + k += value.coeff; + continue; + } + g = gcd(g, value.coeff); + if (g == 1) + return false; + } + return k != 0 && !divides(g, k); + }; + for (auto [x, k] : m_powers) { + SASSERT(k > 1); + bool gcd_fail = true; + dd::pdd kx = m.mk_var(x) * m.mk_val(k); + for (unsigned r = 0; gcd_fail && r < k; r++) { + dd::pdd kx_plus_r = kx + m.mk_val(r); + auto q = p.subst_pdd(x, kx_plus_r); + if (!fails_gcd_test(q)) + gcd_fail = false; + } + if (gcd_fail) { + lemma_builder lemma(c(), "pdd-power"); + add_dependencies(lemma, eq); + return true; + } + } + return false; + } + bool grobner::propagate_quotients() { if (!m_config.m_propagate_quotients) return false; @@ -249,8 +315,6 @@ namespace nla { // z < 0 & x < 0 => -x <= -z bool grobner::propagate_quotients(dd::solver::equation const& eq) { dd::pdd const& p = eq.poly(); - dd::pdd_eval eval; - eval.var2val() = [&](unsigned j) { return val(j); }; if (p.is_linear()) return false; if (p.is_val()) @@ -261,6 +325,8 @@ namespace nla { for (auto v : p.free_vars()) if (!c().var_is_int(v)) return false; + dd::pdd_eval eval; + eval.var2val() = [&](unsigned j) { return val(j); }; if (eval(p) == 0) return false; tracked_uint_set nl_vars; @@ -577,8 +643,12 @@ namespace nla { for (const factor& fc: fcn) q.push_back(var(fc)); } - - if (c().var_is_fixed(j)) + else if (c().lra.column_has_term(j)) { + const lp::lar_term& t = c().lra.get_term(j); + for (auto k : t) + q.push_back(k.j()); + } + else if (c().var_is_fixed(j)) return; const auto& matrix = lra.A_r(); for (auto & s : matrix.m_columns[j]) { @@ -587,12 +657,18 @@ namespace nla { continue; m_rows.insert(row); unsigned k = lra.get_base_column_in_row(row); - // grobner bassis does not know about integer constraints - if (lra.column_is_free(k) && !m_add_all_eqs && k != j) - continue; // a free column over the reals can be assigned - if (lra.column_is_free(k) && k != j && !lra.var_is_int(k)) - continue; + if (lra.column_is_free(k) && k != j) { + if (!lra.var_is_int(k)) + continue; + // free integer columns are ignored unless m_add_all_eqs is set or we are doing gcd test. + if (!m_add_all_eqs && !m_config.m_gcd_test) + continue; + // a free integer column with integer coefficients can be assigned. + if (!m_add_all_eqs && all_of(c().lra.get_row(row), [&](auto& ri) { return ri.coeff().is_int();})) + continue; + } + CTRACE(grobner, matrix.m_rows[row].size() > c().params().arith_nl_grobner_row_length_limit(), tout << "ignore the row " << row << " with the size " << matrix.m_rows[row].size() << "\n";); // limits overhead of grobner equations, unless this is for extracting a complete COI of the non-satisfied subset. @@ -618,6 +694,7 @@ namespace nla { while (!vars.empty()) { j = vars.back(); vars.pop_back(); + if (c().params().arith_nl_grobner_subs_fixed() > 0 && c().var_is_fixed_to_zero(j)) { r = m_pdd_manager.mk_val(val_of_fixed_var_with_deps(j, zero_dep)); dep = zero_dep; diff --git a/src/math/lp/nla_grobner.h b/src/math/lp/nla_grobner.h index d5c018815..61f2f3234 100644 --- a/src/math/lp/nla_grobner.h +++ b/src/math/lp/nla_grobner.h @@ -21,6 +21,7 @@ namespace nla { class grobner : common { struct config { bool m_propagate_quotients = false; + bool m_gcd_test = false; }; dd::pdd_manager m_pdd_manager; dd::solver m_solver; @@ -51,6 +52,10 @@ namespace nla { bool propagate_quotients(); bool propagate_quotients(dd::solver::equation const& eq); + svector> m_powers; + bool propagate_gcd_test(); + bool propagate_gcd_test(dd::solver::equation const& eq); + std::pair linear_to_term(dd::pdd q); void add_dependencies(lemma_builder& lemma, dd::solver::equation const& eq); @@ -74,7 +79,7 @@ namespace nla { bool is_solved(dd::pdd const& p, unsigned& v, dd::pdd& r); void add_eq(dd::pdd& p, u_dependency* dep); const rational& val_of_fixed_var_with_deps(lpvar j, u_dependency*& dep); - dd::pdd pdd_expr(const rational& c, lpvar j, u_dependency*& dep); + dd::pdd pdd_expr(const rational& c, lpvar j, u_dependency*& dep); void display_matrix_of_m_rows(std::ostream& out) const; std::ostream& diagnose_pdd_miss(std::ostream& out); @@ -83,6 +88,6 @@ namespace nla { grobner(core *core); void operator()(); void updt_params(params_ref const& p); - dd::solver::equation_vector const& core_equations(bool all_eqs); + // dd::solver::equation_vector const& core_equations(bool all_eqs); }; } diff --git a/src/math/lp/nla_pp.cpp b/src/math/lp/nla_pp.cpp index 2143427d4..feb41faee 100644 --- a/src/math/lp/nla_pp.cpp +++ b/src/math/lp/nla_pp.cpp @@ -19,13 +19,13 @@ std::ostream& core::print_product(const T& m, std::ostream& out) const { bool first = true; for (lpvar v : m) { if (!first) - out << "*"; + out << " * "; else first = false; if (lp_settings().print_external_var_name()) - out << "(" << lra.get_variable_name(v) << "=" << val(v) << ")"; + out << lra.get_variable_name(v); else - out << "(j" << v << " = " << val(v) << ")"; + out << "j" << v; } return out; } @@ -69,13 +69,13 @@ std::ostream& core::print_factor_with_vars(const factor& f, std::ostream& out) c std::ostream& core::print_monic(const monic& m, std::ostream& out) const { if (lp_settings().print_external_var_name()) - out << "([" << m.var() << "] = " << lra.get_variable_name(m.var()) << " = " << val(m.var()) << " = "; + out << "[" << m.var() << "] := " << lra.get_variable_name(m.var()) << " := "; else - out << "(j" << m.var() << " = " << val(m.var()) << " = "; - print_product(m.vars(), out) << ")\n"; + out << "j" << m.var() << " := "; + print_product(m.vars(), out) << "\n"; return out; } - + std::ostream& core::print_bfc(const factorization& m, std::ostream& out) const { SASSERT(m.size() == 2); out << "( x = " << pp(m[0]) << "* y = " << pp(m[1]) << ")"; @@ -85,6 +85,7 @@ std::ostream& core::print_bfc(const factorization& m, std::ostream& out) const { std::ostream& core::print_monic_with_vars(lpvar v, std::ostream& out) const { return print_monic_with_vars(m_emons[v], out); } + template std::ostream& core::print_product_with_vars(const T& m, std::ostream& out) const { print_product(m, out) << "\n"; @@ -141,9 +142,8 @@ std::ostream& core::print_var(lpvar j, std::ostream& out) const { } std::ostream& core::print_monics(std::ostream& out) const { - for (auto& m : m_emons) { - print_monic_with_vars(m, out); - } + for (auto& m : m_emons) + print_monic_with_vars(m, out); return out; } @@ -312,8 +312,7 @@ std::ostream& core::display_row(std::ostream& out, lp::row_strip const& } std::ostream& core::display(std::ostream& out) { - print_monics(out); - for (unsigned i = 0; i < lra.row_count(); ++i) - display_row(out, lra.get_row(i)); + for (auto& m : m_emons) + print_monic(m, out); return out; } \ No newline at end of file diff --git a/src/params/smt_params_helper.pyg b/src/params/smt_params_helper.pyg index 59e1ab10b..448253354 100644 --- a/src/params/smt_params_helper.pyg +++ b/src/params/smt_params_helper.pyg @@ -79,6 +79,7 @@ def_module_params(module_name='smt', ('arith.nl.grobner_max_simplified', UINT, 10000, 'grobner\'s maximum number of simplifications'), ('arith.nl.grobner_cnfl_to_report', UINT, 1, 'grobner\'s maximum number of conflicts to report'), ('arith.nl.grobner_propagate_quotients', BOOL, False, 'detect conflicts x*y + z = 0 where x doesn\'t divide z'), + ('arith.nl.grobner_gcd_test', BOOL, False, 'detect gcd conflicts for polynomial powers x^k - y = 0'), ('arith.nl.gr_q', UINT, 10, 'grobner\'s quota'), ('arith.nl.grobner_subs_fixed', UINT, 1, '0 - no subs, 1 - substitute, 2 - substitute fixed zeros only'), ('arith.nl.delay', UINT, 10, 'number of calls to final check before invoking bounded nlsat check'),