From d94b2ceac54a234c626f887aad4f0a29607ec69b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 7 Mar 2026 01:48:49 +0000 Subject: [PATCH] Revert multivariate polynomial GCDheu changes, keep integer Lehmer GCD improvement Co-authored-by: levnach <5377127+levnach@users.noreply.github.com> --- src/math/polynomial/polynomial.cpp | 321 ----------------------------- src/test/polynomial.cpp | 160 -------------- src/util/trace_tags.def | 1 - 3 files changed, 482 deletions(-) diff --git a/src/math/polynomial/polynomial.cpp b/src/math/polynomial/polynomial.cpp index daaf2fca2..7bd0fa2d6 100644 --- a/src/math/polynomial/polynomial.cpp +++ b/src/math/polynomial/polynomial.cpp @@ -4691,321 +4691,6 @@ namespace polynomial { gcd_prs(u, v, max_var(u), r); } - // ----------------------------------------------------------------------- - // Heuristic GCD (GCDheu) for multivariate polynomials over Z. - // - // Algorithm by Char, Geddes, Gonnet (1989): - // 1. Choose evaluation base xi > 2 * max|coeff(u,v)|. - // 2. Apply the Kronecker substitution: map variable x_j -> xi^{E_j} - // where E_0=1, E_{j+1} = E_j * (deg_bound_j + 1). - // 3. Evaluate u and v at these integer points to get integers U, V. - // 4. Compute G = gcd(U, V) using integer arithmetic. - // 5. Reconstruct a polynomial g from G using balanced base-xi representation - // (inverse Kronecker substitution). - // 6. Check if g divides both u and v. If yes, return it. - // - // The procedure is tried up to HGCD_MAX_ITERS times with increasing xi. - // Returns true when the GCD was successfully computed, false otherwise. - // ----------------------------------------------------------------------- - - static constexpr unsigned HGCD_MAX_ITERS = 6; - static constexpr unsigned HGCD_E_TOTAL_LIMIT = 1000; // max Kronecker "size" - // Maximum estimated bit-length of xi^E_total. Keeps integer GCD fast. - static constexpr unsigned HGCD_MAX_BITS = 4096; - - // Return max absolute coefficient of p (Lāˆž norm). - void max_abs_coeff(polynomial const * p, numeral & result) { - unsynch_mpz_manager & zm = m().m(); - zm.reset(result); - scoped_numeral tmp(m_manager); - unsigned sz = p->size(); - for (unsigned i = 0; i < sz; ++i) { - zm.set(tmp, p->a(i)); - zm.abs(tmp); - if (zm.gt(tmp, result)) - zm.set(result, tmp); - } - } - - // Evaluate polynomial p at integer substitution xs[j] -> vs[j] for all j. - // Stores the resulting integer in 'val'. - // xs must list every variable that appears in p. - void hgcd_eval(polynomial const * p, - unsigned num_vars, - var const * xs, - numeral const * vs, - numeral & val) { - unsynch_mpz_manager & zm = m().m(); - scoped_numeral term(m_manager), pow_val(m_manager); - zm.set(val, 0); - unsigned sz = p->size(); - // Build a position look-up table from variable id -> xs-index. - // Sized to the largest variable id that appears plus one. - unsigned_vector var2pos; - for (unsigned k = 0; k < num_vars; ++k) { - if (xs[k] >= var2pos.size()) - var2pos.resize(xs[k] + 1, UINT_MAX); - var2pos[xs[k]] = k; - } - for (unsigned i = 0; i < sz; ++i) { - zm.set(term, p->a(i)); - monomial * mon = p->m(i); - unsigned msz = mon->size(); - for (unsigned j = 0; j < msz; ++j) { - var x = mon->get_var(j); - unsigned deg = mon->degree(j); - SASSERT(x < var2pos.size() && var2pos[x] != UINT_MAX); - zm.power(vs[var2pos[x]], deg, pow_val); - zm.mul(term, pow_val, term); - } - zm.add(val, term, val); - } - } - - // Reconstruct a multivariate polynomial from an integer G obtained by the - // Kronecker substitution x_j -> xi^{E_j} with E[j] = prod_{i counter(num_vars, 0u); - - // Temporary for building monomials. - sbuffer pws; - - // Accumulate the polynomial in m_cheap_som_buffer. - SASSERT(m_cheap_som_buffer.empty()); - scoped_numeral digit(m_manager); - - while (!zm.is_zero(H)) { - // Extract balanced digit: digit = H mod xi, centered in (-xi/2, xi/2]. - zm.mod(H, xi, digit); // digit in [0, xi) - if (zm.gt(digit, half_xi)) - zm.sub(digit, xi, digit); // bring into (-xi/2, xi/2] - - // Advance H: H = (H - digit) / xi (exact division). - zm.sub(H, digit, H); - zm.div(H, xi, H); - - // If the digit is non-zero, add it to the result polynomial. - if (!zm.is_zero(digit)) { - // Build monomial x_{vars[0]}^{counter[0]} * ... * x_{vars[k-1]}^{counter[k-1]}. - pws.reset(); - for (unsigned j = 0; j < num_vars; ++j) { - if (counter[j] > 0) - pws.push_back(power(vars[j], counter[j])); - } - monomial * mon = mk_monomial(pws.size(), pws.data()); - m_cheap_som_buffer.add(digit, mon); - } - - // Advance the mixed-radix counter. - for (unsigned j = 0; j < num_vars; ++j) { - counter[j]++; - if (counter[j] <= degree_bounds[j]) - break; // no carry - counter[j] = 0; - if (j + 1 == num_vars) { - // Counter wrapped around: all positions exhausted. - // If H is still non-zero the degree bound was violated. - if (!zm.is_zero(H)) { - m_cheap_som_buffer.reset(); - return false; - } - } - } - } - - result = m_cheap_som_buffer.mk(); - return true; - } - - // Main Heuristic GCD entry point. - // pre: u and v share the same set of variables (non-modular mode). - // post: returns true and stores the GCD in r on success; false otherwise. - bool heuristic_gcd(polynomial const * u, - polynomial const * v, - power_buffer const & u_var_degrees, - power_buffer const & v_var_degrees, - polynomial_ref & r) { - SASSERT(!m().modular()); - SASSERT(u_var_degrees.size() == v_var_degrees.size()); - - unsigned num_vars = u_var_degrees.size(); - - // Only apply heuristic GCD to genuinely multivariate polynomials. - // For univariate inputs, uni_mod_gcd is already efficient. - if (num_vars < 2) return false; - - // ---- Step 1: Compute per-variable degree bounds and Kronecker exponents. - sbuffer vars(num_vars, null_var); - sbuffer deg_bounds(num_vars, 0u); - for (unsigned j = 0; j < num_vars; ++j) { - SASSERT(u_var_degrees[j].get_var() == v_var_degrees[j].get_var()); - vars[j] = u_var_degrees[j].get_var(); - deg_bounds[j] = std::min(u_var_degrees[j].degree(), - v_var_degrees[j].degree()); - } - - // E[j] = product of (deg_bounds[i]+1) for i < j. - sbuffer E(num_vars, 0u); - E[0] = 1u; - for (unsigned j = 1; j < num_vars; ++j) { - // Guard against overflow / excessively large evaluation points. - if (E[j-1] > HGCD_E_TOTAL_LIMIT / (deg_bounds[j-1] + 1u)) - return false; - E[j] = E[j-1] * (deg_bounds[j-1] + 1u); - } - // Total Kronecker size = E[num_vars-1] * (deg_bounds[num_vars-1]+1). - unsigned E_total; - { - unsigned last = num_vars - 1; - if (E[last] > HGCD_E_TOTAL_LIMIT / (deg_bounds[last] + 1u)) - return false; - E_total = E[last] * (deg_bounds[last] + 1u); - } - - // ---- Step 2: Compute xi = 2 * max(max_abs_coeff(u), max_abs_coeff(v)) + 2. - // - // We use the Lāˆž norm (max absolute coefficient) rather than the L1 norm, - // so that xi is minimal. Smaller xi means smaller integers throughout, - // which is the dominant cost driver. - unsynch_mpz_manager & zm = m().m(); - scoped_numeral norm_u(m_manager), norm_v(m_manager), xi(m_manager); - max_abs_coeff(u, norm_u); - max_abs_coeff(v, norm_v); - if (zm.lt(norm_u, norm_v)) - zm.set(xi, norm_v); - else - zm.set(xi, norm_u); - zm.mul2k(xi, 1); // xi = 2 * max_coeff - zm.inc(xi); - zm.inc(xi); // xi = 2 * max_coeff + 2 (ensures xi >= 2) - - // Safety check: make sure the integers we'll deal with are not astronomically - // large. The evaluation image |U_val| is bounded by roughly xi^E_total. - // Estimate the bit length of xi^E_total and bail if it would exceed the limit. - { - // Approximate bit length of xi: log2(xi) ā‰ˆ number of significant bits. - // Use zm.log2 if available, otherwise fall back to a simpler check. - // If xi fits in 64 bits, compute precisely; otherwise it's already big. - if (zm.is_uint64(xi)) { - uint64_t xi_val = zm.get_uint64(xi); - unsigned xi_bits = 1; - while (xi_val >>= 1) ++xi_bits; - if ((uint64_t)xi_bits * E_total > HGCD_MAX_BITS) - return false; - } - // If xi doesn't fit in 64 bits it's definitely too large. - else { return false; } - } - - // ---- Step 3: Extract integer contents of u and v so we work with - // primitive parts during the reconstruction. - scoped_numeral c_u(m_manager), c_v(m_manager), d_a(m_manager); - polynomial_ref pp_u(m_wrapper), pp_v(m_wrapper); - ic(u, c_u, pp_u); - ic(v, c_v, pp_v); - zm.gcd(c_u, c_v, d_a); - - // Pre-compute xi^{E[j]} for each variable substitution. - scoped_numeral_vector xi_pows(m_manager); - xi_pows.resize(num_vars); - for (unsigned j = 0; j < num_vars; ++j) - zm.power(xi, E[j], xi_pows[j]); - - // ---- Step 4: Iterate with increasing xi. - // The integer GCD gcd(eval(pp_u), eval(pp_v)) may equal k * eval(gcd_pp) - // for some integer k that is a power of xi (arising from common factors in - // the Kronecker images). We therefore try stripping leading powers of xi - // from the integer GCD before reconstruction, which peels off such spurious - // factors and recovers the true encoding of the polynomial GCD. - scoped_numeral U_val(m_manager), V_val(m_manager), G_val(m_manager); - scoped_numeral G_try(m_manager), rem_tmp(m_manager); - polynomial_ref candidate(m_wrapper); - - for (unsigned iter = 0; iter < HGCD_MAX_ITERS; ++iter) { - checkpoint(); - - // Re-check bit-size safety after xi is doubled. - if (!zm.is_uint64(xi)) return false; - { - uint64_t xi_val = zm.get_uint64(xi); - unsigned xi_bits = 1; - while (xi_val >>= 1) ++xi_bits; - if ((uint64_t)xi_bits * E_total > HGCD_MAX_BITS) - return false; - } - - // Evaluate the primitive parts at the Kronecker substitution. - hgcd_eval(pp_u, num_vars, vars.data(), xi_pows.data(), U_val); - hgcd_eval(pp_v, num_vars, vars.data(), xi_pows.data(), V_val); - - // Integer GCD. - zm.gcd(U_val, V_val, G_val); - zm.abs(G_val); - - // Try reconstruction with G_try = G_val, G_val/xi, G_val/xi^2, ... - // Stripping factors of xi exposes the "true" Kronecker encoding when - // the integer GCD has absorbed common xi-powers from both evaluations. - zm.set(G_try, G_val); - for (unsigned strip = 0; strip <= E_total; ++strip) { - if (zm.is_zero(G_try)) break; - - if (hgcd_reconstruct(G_try, xi, num_vars, vars.data(), - deg_bounds.data(), candidate)) { - // Make the candidate primitive (remove integer content). - scoped_numeral cont(m_manager); - polynomial_ref pp_cand(m_wrapper); - ic(candidate, cont, pp_cand); - - // Verify: pp_cand must divide both pp_u and pp_v. - if (!is_zero(pp_cand) && - divides(pp_cand, pp_u) && divides(pp_cand, pp_v)) { - r = mul(d_a, pp_cand); - flip_sign_if_lm_neg(r); - TRACE(hgcd, tout << "heuristic_gcd succeeded, iter=" << iter - << " strip=" << strip << "\nr: " << r << "\n";); - return true; - } - } - - // Divide G_try by xi for the next inner attempt. - zm.rem(G_try, xi, rem_tmp); - if (!zm.is_zero(rem_tmp)) break; // xi no longer divides G_try - zm.div(G_try, xi, G_try); - } - - // Double xi for the next outer attempt. - zm.mul2k(xi, 1); - // Recompute xi_pows. - for (unsigned j = 0; j < num_vars; ++j) - zm.power(xi, E[j], xi_pows[j]); - } - - TRACE(hgcd, tout << "heuristic_gcd failed\n";); - return false; - } - void gcd(polynomial const * u, polynomial const * v, polynomial_ref & r) { power_buffer u_var_degrees; power_buffer v_var_degrees; @@ -5123,12 +4808,6 @@ namespace polynomial { #endif if (!m().modular() && !m_use_prs_gcd) { SASSERT(max_var(u) == max_var(v)); - - // Fast path: try heuristic GCD first. - // It succeeds quickly for polynomials with small coefficients / low degree. - if (heuristic_gcd(u, v, u_var_degrees, v_var_degrees, r)) - return; - if (is_univariate(u)) { SASSERT(is_univariate(v)); uni_mod_gcd(u, v, r); diff --git a/src/test/polynomial.cpp b/src/test/polynomial.cpp index 1eaa5c343..adea1521a 100644 --- a/src/test/polynomial.cpp +++ b/src/test/polynomial.cpp @@ -23,7 +23,6 @@ Notes: #include "math/polynomial/linear_eq_solver.h" #include "util/rlimit.h" #include -#include static void tst1() { std::cout << "\n----- Basic testing -------\n"; @@ -1813,162 +1812,6 @@ static void tst_divides() { std::cout << "divides(q, p): " << m.divides(q, p) << "\n"; } -// ----------------------------------------------------------------------- -// tst_hgcd: dedicated tests for the Heuristic GCD (GCDheu) algorithm. -// These cover correctness for univariate and multivariate cases and provide -// a simple timing comparison against the existing modular-GCD baseline. -// ----------------------------------------------------------------------- - -static void tst_hgcd() { - std::cout << "\n----- tst_hgcd -------\n"; - polynomial::numeral_manager nm; - reslimit rl; - polynomial::manager m(rl, nm); - - polynomial_ref x0(m), x1(m), x2(m), x3(m); - x0 = m.mk_polynomial(m.mk_var()); - x1 = m.mk_polynomial(m.mk_var()); - x2 = m.mk_polynomial(m.mk_var()); - x3 = m.mk_polynomial(m.mk_var()); - - polynomial_ref one(m); - one = m.mk_const(mpz(1)); - - // ---- univariate tests ---- - - // gcd(x^2-1, x-1) = x-1 - tst_gcd((x0^2) - 1, x0 - 1, x0 - 1); - - // gcd(x^2 - 4x + 4, x^2 - 2x) = x-2 (factor (x-2)^2 vs x(x-2)) - tst_gcd((x0^2) - 4*x0 + 4, (x0^2) - 2*x0, x0 - 2); - - // gcd(x^3 - x, x^2 - 1) = x-1 ... wait, gcd(x(x-1)(x+1), (x-1)(x+1)) = (x-1)(x+1) - tst_gcd((x0^3) - x0, (x0^2) - 1, (x0^2) - 1); - - // gcd with content: gcd(6x^2, 9x) = 3x - { - polynomial_ref three(m), three_x(m); - three = m.mk_const(mpz(3)); - three_x = three * x0; - tst_gcd(6*(x0^2), 9*x0, three_x); - } - - // gcd(f, f) = f - { - polynomial_ref f(m); - f = (x0^3) - 2*(x0^2) + x0 - 5; - tst_gcd(f, f, f); - } - - // gcd(f, 1) = 1 - tst_gcd((x0^5) + x0 + 1, one, one); - - // gcd(f, 0) = f - { - polynomial_ref zero(m), f(m); - zero = m.mk_const(mpz(0)); - f = (x0^3) + x0 + 7; - tst_gcd(f, zero, f); - } - - // ---- bivariate tests ---- - - // gcd((x+y)*(x-y), (x+y)^2) = x+y - tst_gcd((x0 + x1)*(x0 - x1), (x0 + x1)*(x0 + x1), x0 + x1); - - // gcd((x+y)(2x+3y), (x+y)(x-y)) = x+y - tst_gcd((x0 + x1)*(2*x0 + 3*x1), (x0 + x1)*(x0 - x1), x0 + x1); - - // gcd(x^2 + 2xy + y^2, x^2 - y^2) = x+y ((x+y)^2, (x+y)(x-y)) - tst_gcd((x0^2) + 2*x0*x1 + (x1^2), (x0^2) - (x1^2), x0 + x1); - - // gcd with non-trivial content: gcd(6(x+y)(x+1), 9(x+y)(y-1)) = 3(x+y) - { - polynomial_ref three(m), g_expected(m); - three = m.mk_const(mpz(3)); - g_expected = three * (x0 + x1); - tst_gcd(6*(x0 + x1)*(x0 + 1), 9*(x0 + x1)*(x1 - 1), g_expected); - } - - // gcd where result = 1 - tst_gcd((x0^2) + x1 + 1, x0 + (x1^2) + 3, one); - - // gcd of same polynomial (bivariate, leading monomial has positive coefficient) - { - polynomial_ref f(m); - f = x0*(x1^2) - (x0^2)*x1 + 3*x0 + 5; - tst_gcd(f, f, f); - } - - // ---- trivariate tests ---- - - // gcd((x+y+z)(x-y), (x+y+z)(x+z)) = x+y+z - tst_gcd((x0+x1+x2)*(x0-x1), (x0+x1+x2)*(x0+x2), x0+x1+x2); - - // gcd((x0+1)(x1+1)(x2+1), (x0+1)(x1-1)(x2+1)) = (x0+1)(x2+1) - tst_gcd((x0+1)*(x1+1)*(x2+1), (x0+1)*(x1-1)*(x2+1), (x0+1)*(x2+1)); - - // ---- 4-variable tests ---- - - // gcd((x0^2 + x1*x2 + 1)*(x3*x1 + 2), (x0^2 + x1*x2 + 1)*(x3*x1 - 3)) - // = x0^2 + x1*x2 + 1 - tst_gcd(((x0^2) + x1*x2 + 1)*(x3*x1 + 2), - ((x0^2) + x1*x2 + 1)*(x3*x1 - 3), - (x0^2) + x1*x2 + 1); - - // gcd with scalar content - tst_gcd(15*((x0^2) + x1*x2 + 1)*(x3*x1 + 2), - 10*((x0^2) + x1*x2 + 1)*(x3*x1 - 3), - 5*((x0^2) + x1*x2 + 1)); - - std::cout << "tst_hgcd: all correctness checks passed\n"; - - // ---- simple performance comparison ---- - // Compute the same GCD 1000 times and report wall time. - { - std::cout << "\n--- Performance: GCD repeated 1000 times ---\n"; - polynomial_ref u(m), v(m), g(m); - // Use a moderately complex 4-variable GCD: - // u = (x0^2 + x1*x2 + 1)*(x2*x2 + x3 + 2)*(x3*x1 + 2)*(x3*x1^2 + x1*x2 + 1) - // v = (x0^2 + x1*x2 + 1)*(x3*x1^2 + x1*x2 + 1)*(x3*x1 + x1*x2 + 17) - u = ((x0^2) + x1*x2 + 1)*(x2*x2 + x3 + 2)*(x3*x1 + 2)*(x3*(x1^2) + x1*x2 + 1); - v = ((x0^2) + x1*x2 + 1)*(x3*(x1^2) + x1*x2 + 1)*(x3*x1 + x1*x2 + 17); - unsigned N = 1000; - auto t0 = std::chrono::high_resolution_clock::now(); - for (unsigned i = 0; i < N; ++i) - g = gcd(u, v); - auto t1 = std::chrono::high_resolution_clock::now(); - double ms = std::chrono::duration(t1 - t0).count(); - std::cout << "GCD result: " << g << "\n"; - std::cout << N << " GCD calls: " << ms << " ms (" - << ms/N << " ms/call)\n"; - } - - // ---- dense bivariate performance ---- - { - std::cout << "\n--- Performance: dense bivariate GCD ---\n"; - // Build u = A*B, v = A*C where A, B, C have moderate degree. - polynomial_ref A(m), B(m), C(m), u2(m), v2(m), g2(m); - // A = x0^4 + x1^3 + x0*x1^2 + 2 - A = (x0^4) + (x1^3) + x0*(x1^2) + 2; - // B = x0^3 - x1^2 + x0*x1 + 1 - B = (x0^3) - (x1^2) + x0*x1 + 1; - // C = x0^2 + x1^4 - x0*x1 + 3 - C = (x0^2) + (x1^4) - x0*x1 + 3; - u2 = A * B; - v2 = A * C; - unsigned N2 = 500; - auto t0 = std::chrono::high_resolution_clock::now(); - for (unsigned i = 0; i < N2; ++i) - g2 = gcd(u2, v2); - auto t1 = std::chrono::high_resolution_clock::now(); - double ms = std::chrono::duration(t1 - t0).count(); - std::cout << "GCD result: " << g2 << "\n"; - std::cout << N2 << " dense-bivariate GCD calls: " << ms << " ms (" - << ms/N2 << " ms/call)\n"; - } -} - void tst_polynomial() { set_verbosity_level(1000); // enable_trace("factor"); @@ -1979,9 +1822,6 @@ void tst_polynomial() { // enable_trace("Lazard"); // enable_trace("eval_bug"); // enable_trace("mgcd"); - tst_gcd2(); - tst_gcd(); - tst_hgcd(); tst_psc(); return; tst_eval(); diff --git a/src/util/trace_tags.def b/src/util/trace_tags.def index 45dfc9314..1ad305c2d 100644 --- a/src/util/trace_tags.def +++ b/src/util/trace_tags.def @@ -582,7 +582,6 @@ X(Global, mgcd, "mgcd") X(Global, mgcd_bug, "mgcd bug") X(Global, mgcd_call, "mgcd call") X(Global, mgcd_detail, "mgcd detail") -X(Global, hgcd, "heuristic gcd") X(Global, missing_instance_detail, "missing instance detail") X(Global, missing_propagation, "missing propagation") X(Global, mk_and_bug, "mk and bug")