From 44e84dc5d06ef3a317a96164a8f9f686f99b275d Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Mon, 23 Mar 2026 14:18:42 -1000 Subject: [PATCH] refactor try_bivar_hensel_lift and outline the algorithm Signed-off-by: Lev Nachmanson --- src/math/lp/hnf.h | 12 +- src/math/polynomial/polynomial.cpp | 419 +++++++++++++++++------------ src/test/lp/lp.cpp | 2 +- 3 files changed, 262 insertions(+), 171 deletions(-) diff --git a/src/math/lp/hnf.h b/src/math/lp/hnf.h index 527d3f681..5ee6301ab 100644 --- a/src/math/lp/hnf.h +++ b/src/math/lp/hnf.h @@ -613,11 +613,13 @@ public: #endif calculate_by_modulo(); #ifdef Z3DEBUG - CTRACE(hnf_calc, m_H != m_W, - tout << "A = "; m_A_orig.print(tout, 4); tout << std::endl; - tout << "H = "; m_H.print(tout, 4); tout << std::endl; - tout << "W = "; m_W.print(tout, 4); tout << std::endl;); - SASSERT (m_H == m_W); + if (!m_cancelled) { + CTRACE(hnf_calc, m_H != m_W, + tout << "A = "; m_A_orig.print(tout, 4); tout << std::endl; + tout << "H = "; m_H.print(tout, 4); tout << std::endl; + tout << "W = "; m_W.print(tout, 4); tout << std::endl;); + SASSERT (m_H == m_W); + } #endif } diff --git a/src/math/polynomial/polynomial.cpp b/src/math/polynomial/polynomial.cpp index 4cf5417f6..1b06fb2b0 100644 --- a/src/math/polynomial/polynomial.cpp +++ b/src/math/polynomial/polynomial.cpp @@ -6964,22 +6964,249 @@ namespace polynomial { } } - // Bivariate Hensel lifting for multivariate factorization. - // Given q(x, y) with q(x, 0) = lc_val * h1(x) * h2(x) where h1, h2 are coprime monic, - // lift to q(x, y) = F1(x, y) * F2(x, y). - // Works modulo a single prime p: all arithmetic is in Zp[x], then coefficients are - // converted to centered representatives in (-p/2, p/2). This succeeds when p is large - // enough that all true integer coefficients of the factors lie in that range. - // A Mignotte-style bound could be used to choose p more precisely, but for now we - // rely on the caller trying increasingly large primes from a hardcoded list. - // Returns true on success, false if this split doesn't yield a factorization. - bool try_bivar_hensel_lift( - polynomial const * q, // bivariate poly, q(x, 0) has known factorization + // Convert Zp-lifted coefficient arrays to centered Z representatives, + // build multivariate polynomials F1(x,y) and F2(x,y), and verify q == F1*F2. + // Returns true on success. Does NOT deallocate the coefficient vectors. + bool reconstruct_lifted_factors( + polynomial const * q, var x, var y, unsigned deg_y, - upolynomial::numeral_vector const & uf1_monic, // first monic univariate factor - upolynomial::numeral_vector const & uf2_monic, // second monic univariate factor - numeral const & lc_val, // leading coefficient value: lc(q, x)(y=0) + uint64_t prime, + vector const & q_coeffs, + vector const & F1_coeffs, + vector const & F2_coeffs, + polynomial_ref & F1_out, + polynomial_ref & F2_out) { + + scoped_numeral p_num(m_manager); + m_manager.set(p_num, static_cast(prime)); + scoped_numeral half_p(m_manager); + m_manager.set(half_p, static_cast(prime)); + m_manager.div(half_p, mpz(2), half_p); + + polynomial_ref F1_poly(pm()), F2_poly(pm()); + F1_poly = mk_zero(); + F2_poly = mk_zero(); + + for (unsigned j = 0; j <= deg_y; j++) { + // Center the coefficients: if coeff > p/2, subtract p + for (unsigned i = 0; i < F1_coeffs[j]->size(); i++) + if (m_manager.gt((*F1_coeffs[j])[i], half_p)) + m_manager.sub((*F1_coeffs[j])[i], p_num, (*F1_coeffs[j])[i]); + for (unsigned i = 0; i < F2_coeffs[j]->size(); i++) + if (m_manager.gt((*F2_coeffs[j])[i], half_p)) + m_manager.sub((*F2_coeffs[j])[i], p_num, (*F2_coeffs[j])[i]); + + // Build y^j * F_coeffs[j](x) + if (!F1_coeffs[j]->empty()) { + polynomial_ref univ(pm()); + univ = to_polynomial(F1_coeffs[j]->size(), F1_coeffs[j]->data(), x); + if (!is_zero(univ)) { + monomial_ref yj(pm()); + yj = mk_monomial(y, j); + polynomial_ref term(pm()); + term = mul(yj, univ); + F1_poly = add(F1_poly, term); + } + } + if (!F2_coeffs[j]->empty()) { + polynomial_ref univ(pm()); + univ = to_polynomial(F2_coeffs[j]->size(), F2_coeffs[j]->data(), x); + if (!is_zero(univ)) { + monomial_ref yj(pm()); + yj = mk_monomial(y, j); + polynomial_ref term(pm()); + term = mul(yj, univ); + F2_poly = add(F2_poly, term); + } + } + } + + // Verify: q == F1 * F2 over Z[x, y] + polynomial_ref product(pm()); + product = mul(F1_poly, F2_poly); + if (!eq(product, q)) + return false; + + F1_out = F1_poly; + F2_out = F2_poly; + return true; + } + + // Hensel lifting loop: for j = 1..deg_y, compute F1_coeffs[j] and F2_coeffs[j] + // using Bezout coefficients s, t such that s*f1 + t*f2 = 1 in Zp[x]. + // F1_coeffs[0] and F2_coeffs[0] must already be initialized. + void hensel_lift_coefficients( + upolynomial::zp_manager & zp_upm, + unsigned deg_y, + upolynomial::scoped_numeral_vector const & f1_p, + upolynomial::scoped_numeral_vector const & f2_p, + upolynomial::scoped_numeral_vector const & s_vec, + upolynomial::scoped_numeral_vector const & t_vec, + numeral const & lc_val, + vector const & q_coeffs, + vector & F1_coeffs, + vector & F2_coeffs) { + + auto & nm = upm().m(); + auto & znm = zp_upm.m(); + + scoped_numeral lc_inv(m_manager); + znm.set(lc_inv, lc_val); + znm.inv(lc_inv); + + for (unsigned j = 1; j <= deg_y; j++) { + checkpoint(); + + // Compute e_j = q_coeffs[j] - sum_{a+b=j, a>0, b>0} F1_coeffs[a] * F2_coeffs[b] + upolynomial::scoped_numeral_vector e_j(nm); + if (j < q_coeffs.size() && !q_coeffs[j]->empty()) + zp_upm.set(q_coeffs[j]->size(), q_coeffs[j]->data(), e_j); + + for (unsigned a = 0; a <= j; a++) { + unsigned b = j - a; + if (b > deg_y) continue; + if (a == j && b == 0) continue; + if (a == 0 && b == j) continue; + if (F1_coeffs[a]->empty() || F2_coeffs[b]->empty()) continue; + + upolynomial::scoped_numeral_vector prod(nm); + zp_upm.mul(F1_coeffs[a]->size(), F1_coeffs[a]->data(), + F2_coeffs[b]->size(), F2_coeffs[b]->data(), prod); + upolynomial::scoped_numeral_vector new_e(nm); + zp_upm.sub(e_j.size(), e_j.data(), prod.size(), prod.data(), new_e); + e_j.swap(new_e); + } + + if (e_j.empty()) continue; + + // Solve using Bezout coefficients: A = (e_j * t) mod f1 + upolynomial::scoped_numeral_vector e_t(nm), A_val(nm), Q_tmp(nm); + zp_upm.mul(e_j.size(), e_j.data(), t_vec.size(), t_vec.data(), e_t); + zp_upm.div_rem(e_t.size(), e_t.data(), f1_p.size(), f1_p.data(), Q_tmp, A_val); + + // B = (e_j * s) mod f2 + upolynomial::scoped_numeral_vector e_s(nm), B_val(nm); + zp_upm.mul(e_j.size(), e_j.data(), s_vec.size(), s_vec.data(), e_s); + zp_upm.div_rem(e_s.size(), e_s.data(), f2_p.size(), f2_p.data(), Q_tmp, B_val); + + // F1[j] = A * lc_inv, F2[j] = B + zp_upm.mul(A_val, lc_inv); + zp_upm.set(A_val.size(), A_val.data(), *F1_coeffs[j]); + zp_upm.set(B_val.size(), B_val.data(), *F2_coeffs[j]); + } + } + + // Extract coefficients of q w.r.t. y as upolynomials in Zp[x], and initialize + // the lifted factor coefficient arrays with F1[0] = f1, F2[0] = lc_val * f2. + // Returns false if q is not truly bivariate in x and y. + bool extract_and_init_lift_coefficients( + upolynomial::zp_manager & zp_upm, + polynomial const * q, + var y, + unsigned deg_y, + upolynomial::scoped_numeral_vector const & f1_p, + upolynomial::scoped_numeral_vector const & f2_p, + numeral const & lc_val, + vector & q_coeffs, + vector & F1_coeffs, + vector & F2_coeffs) { + + auto & nm = upm().m(); + auto & znm = zp_upm.m(); + + for (unsigned j = 0; j <= deg_y; j++) { + polynomial_ref cj(pm()); + cj = coeff(q, y, j); + auto * vec = alloc(upolynomial::scoped_numeral_vector, nm); + if (!is_zero(cj) && is_univariate(cj)) + upm().to_numeral_vector(cj, *vec); + else if (!is_zero(cj) && is_const(cj)) { + vec->push_back(numeral()); + nm.set(vec->back(), cj->a(0)); + } + else if (!is_zero(cj)) { + dealloc(vec); + return false; + } + for (unsigned i = 0; i < vec->size(); i++) + znm.p_normalize((*vec)[i]); + zp_upm.trim(*vec); + q_coeffs.push_back(vec); + } + + // Initialize lifted factor coefficient arrays + for (unsigned j = 0; j <= deg_y; j++) { + F1_coeffs.push_back(alloc(upolynomial::scoped_numeral_vector, nm)); + F2_coeffs.push_back(alloc(upolynomial::scoped_numeral_vector, nm)); + } + + // F1[0] = f1, F2[0] = lc_val * f2 + zp_upm.set(f1_p.size(), f1_p.data(), *F1_coeffs[0]); + scoped_numeral lc_p(m_manager); + znm.set(lc_p, lc_val); + upolynomial::scoped_numeral_vector lc_f2(nm); + zp_upm.set(f2_p.size(), f2_p.data(), lc_f2); + zp_upm.mul(lc_f2, lc_p); + zp_upm.set(lc_f2.size(), lc_f2.data(), *F2_coeffs[0]); + return true; + } + + // Bivariate Hensel lifting for multivariate factorization. + // + // Mathematical setup: + // We have q(x, y) in Z[x, y] with degree deg_y in y. + // Evaluating at y = 0 gives a univariate factorization + // q(x, 0) = lc_val * uf1_monic(x) * uf2_monic(x) + // where uf1_monic and uf2_monic are monic, coprime polynomials in Z[x], + // and lc_val = lc(q, x)(y=0) is an integer. + // + // Goal: lift to q(x, y) = F1(x, y) * F2(x, y) over Z[x, y]. + // + // Method (linear Hensel lifting in Zp[x]): + // 1. Reduce uf1_monic, uf2_monic to f1, f2 in Zp[x] and compute + // Bezout coefficients s, t with s*f1 + t*f2 = 1 in Zp[x]. + // This requires gcd(f1, f2) = 1 in Zp[x], i.e. the prime p + // must not divide the resultant of f1, f2. + // 2. Expand q, F1, F2 as polynomials in y with coefficients in Zp[x]: + // q = q_0(x) + q_1(x)*y + ... + q_{deg_y}(x)*y^{deg_y} + // F1 = F1_0(x) + F1_1(x)*y + ... + // F2 = F2_0(x) + F2_1(x)*y + ... + // The y^0 coefficients are known: F1_0 = f1, F2_0 = lc_val * f2. + // 3. For j = 1, ..., deg_y, matching the y^j coefficient of q = F1 * F2: + // q_j = sum_{a+b=j} F1_a * F2_b + // The unknowns are F1_j and F2_j. Set + // e_j = q_j - sum_{a+b=j, 0 q_coeffs; - for (unsigned j = 0; j <= deg_y; j++) { - polynomial_ref cj(pm()); - cj = coeff(q, y, j); - auto * vec = alloc(upolynomial::scoped_numeral_vector, nm); - if (!is_zero(cj) && is_univariate(cj)) { - upm().to_numeral_vector(cj, *vec); - } - else if (!is_zero(cj) && is_const(cj)) { - vec->push_back(numeral()); - nm.set(vec->back(), cj->a(0)); - } - else if (!is_zero(cj)) { - // q is not bivariate, abort - for (auto * v : q_coeffs) dealloc(v); - dealloc(vec); - return false; - } - // Convert to Zp - for (unsigned i = 0; i < vec->size(); i++) - znm.p_normalize((*vec)[i]); - zp_upm.trim(*vec); - q_coeffs.push_back(vec); + vector q_coeffs, F1_coeffs, F2_coeffs; + if (!extract_and_init_lift_coefficients(zp_upm, q, y, deg_y, + f1_p, f2_p, lc_val, + q_coeffs, F1_coeffs, F2_coeffs)) { + for (auto * v : q_coeffs) dealloc(v); + return false; } - // Initialize lifted factor coefficient arrays - // F1[j], F2[j] are the coefficient of y^j in the lifted factors, as upolynomials in Zp[x] - vector F1_coeffs, F2_coeffs; - for (unsigned j = 0; j <= deg_y; j++) { - F1_coeffs.push_back(alloc(upolynomial::scoped_numeral_vector, nm)); - F2_coeffs.push_back(alloc(upolynomial::scoped_numeral_vector, nm)); - } + hensel_lift_coefficients(zp_upm, deg_y, f1_p, f2_p, + s_vec, t_vec, lc_val, + q_coeffs, F1_coeffs, F2_coeffs); - // F1[0] = f1, F2[0] = lc_val * f2 (absorb leading coefficient) - zp_upm.set(f1_p.size(), f1_p.data(), *F1_coeffs[0]); - scoped_numeral lc_p(m_manager); - znm.set(lc_p, lc_val); - upolynomial::scoped_numeral_vector lc_f2(nm); - zp_upm.set(f2_p.size(), f2_p.data(), lc_f2); - zp_upm.mul(lc_f2, lc_p); - zp_upm.set(lc_f2.size(), lc_f2.data(), *F2_coeffs[0]); + bool ok = reconstruct_lifted_factors(q, x, y, deg_y, prime, + q_coeffs, F1_coeffs, F2_coeffs, + F1_out, F2_out); - // Hensel lifting: for j = 1, ..., deg_y - // At each step, solve F1[j]*F2[0] + F2[j]*F1[0] = e_j for the new coefficients. - // Since F1[0]=f1, F2[0]=lc*f2, and s*f1 + t*f2 = 1 (Bézout), we get: - // A = (e_j * t) mod f1, B = (e_j * s) mod f2 - // F1[j] = A * lc_inv, F2[j] = B - - // Precompute lc_inv = lc_val^{-1} in Zp, used in every lifting step - scoped_numeral lc_inv(m_manager); - znm.set(lc_inv, lc_val); - znm.inv(lc_inv); - - for (unsigned j = 1; j <= deg_y; j++) { - checkpoint(); - - // Compute e_j = q_coeffs[j] - sum_{a+b=j, 0empty()) - zp_upm.set(q_coeffs[j]->size(), q_coeffs[j]->data(), e_j); - - for (unsigned a = 0; a <= j; a++) { - unsigned b = j - a; - if (b > deg_y) continue; - if (a == j && b == 0) continue; // skip F1[j]*F2[0], that's what we're computing - if (a == 0 && b == j) continue; // skip F1[0]*F2[j] - if (F1_coeffs[a]->empty() || F2_coeffs[b]->empty()) continue; - - upolynomial::scoped_numeral_vector prod(nm); - zp_upm.mul(F1_coeffs[a]->size(), F1_coeffs[a]->data(), - F2_coeffs[b]->size(), F2_coeffs[b]->data(), prod); - upolynomial::scoped_numeral_vector new_e(nm); - zp_upm.sub(e_j.size(), e_j.data(), prod.size(), prod.data(), new_e); - e_j.swap(new_e); - } - - if (e_j.empty()) continue; - - // Solve A*f2 + B*f1 = e_j using Bézout coefficients - // A = (e_j * t) mod f1 - upolynomial::scoped_numeral_vector e_t(nm), A_val(nm), Q_tmp(nm); - zp_upm.mul(e_j.size(), e_j.data(), t_vec.size(), t_vec.data(), e_t); - zp_upm.div_rem(e_t.size(), e_t.data(), f1_p.size(), f1_p.data(), Q_tmp, A_val); - - // B = (e_j * s) mod f2 - upolynomial::scoped_numeral_vector e_s(nm), B_val(nm); - zp_upm.mul(e_j.size(), e_j.data(), s_vec.size(), s_vec.data(), e_s); - zp_upm.div_rem(e_s.size(), e_s.data(), f2_p.size(), f2_p.data(), Q_tmp, B_val); - - // F1[j] = A * lc_inv in Zp - zp_upm.mul(A_val, lc_inv); - zp_upm.set(A_val.size(), A_val.data(), *F1_coeffs[j]); - - // F2[j] = B - zp_upm.set(B_val.size(), B_val.data(), *F2_coeffs[j]); - } - - // Convert Zp coefficients to centered Z representatives and build multivariate polynomials - // F1(x, y) = sum_j F1_coeffs[j](x) * y^j - // For centered rep: if coeff > p/2, subtract p - scoped_numeral half_p(m_manager); - m_manager.set(half_p, static_cast(prime)); - m_manager.div(half_p, mpz(2), half_p); - - polynomial_ref F1_poly(pm()), F2_poly(pm()); - F1_poly = mk_zero(); - F2_poly = mk_zero(); - - for (unsigned j = 0; j <= deg_y; j++) { - // Center the coefficients - for (unsigned i = 0; i < F1_coeffs[j]->size(); i++) { - if (m_manager.gt((*F1_coeffs[j])[i], half_p)) - m_manager.sub((*F1_coeffs[j])[i], p_num, (*F1_coeffs[j])[i]); - } - for (unsigned i = 0; i < F2_coeffs[j]->size(); i++) { - if (m_manager.gt((*F2_coeffs[j])[i], half_p)) - m_manager.sub((*F2_coeffs[j])[i], p_num, (*F2_coeffs[j])[i]); - } - - // Build y^j * F_coeffs[j](x) - if (!F1_coeffs[j]->empty()) { - polynomial_ref univ(pm()); - univ = to_polynomial(F1_coeffs[j]->size(), F1_coeffs[j]->data(), x); - if (!is_zero(univ)) { - monomial_ref yj(pm()); - yj = mk_monomial(y, j); - polynomial_ref term(pm()); - term = mul(yj, univ); - F1_poly = add(F1_poly, term); - } - } - if (!F2_coeffs[j]->empty()) { - polynomial_ref univ(pm()); - univ = to_polynomial(F2_coeffs[j]->size(), F2_coeffs[j]->data(), x); - if (!is_zero(univ)) { - monomial_ref yj(pm()); - yj = mk_monomial(y, j); - polynomial_ref term(pm()); - term = mul(yj, univ); - F2_poly = add(F2_poly, term); - } - } - } - - // Cleanup allocated vectors for (auto * v : q_coeffs) dealloc(v); for (auto * v : F1_coeffs) dealloc(v); for (auto * v : F2_coeffs) dealloc(v); - // Verify: q == F1 * F2 over Z[x, y] - polynomial_ref product(pm()); - product = mul(F1_poly, F2_poly); - if (!eq(product, q)) - return false; - - F1_out = F1_poly; - F2_out = F2_poly; - return true; + return ok; } // Evaluation points used for multivariate factorization diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 160caaa46..ee1989c20 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -1711,7 +1711,7 @@ void test_dio() { enable_trace("dioph_eq"); enable_trace("dioph_eq_fresh"); #ifdef Z3DEBUG - auto r = i_solver.dio_test(); + i_solver.dio_test(); #endif }