mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2026-03-26 06:18:38 +00:00
Code Activity

refactor try_bivar_hensel_lift and outline the algorithm

Browse source 
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2026-03-23 14:18:42 -10:00 • committed by Lev Nachmanson
parent 117da362f0
commit 44e84dc5d0
3 changed files with 262 additions and 171 deletions
Download patch file Download diff file Expand all files Collapse all files

12
src/math/lp/hnf.h
View file

@ -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
}

419
src/math/polynomial/polynomial.cpp
View file

@ -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<upolynomial::scoped_numeral_vector*> const & q_coeffs,
vector<upolynomial::scoped_numeral_vector*> const & F1_coeffs,
vector<upolynomial::scoped_numeral_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<int>(prime));
scoped_numeral half_p(m_manager);
m_manager.set(half_p, static_cast<int>(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<upolynomial::scoped_numeral_vector*> const & q_coeffs,
vector<upolynomial::scoped_numeral_vector*> & F1_coeffs,
vector<upolynomial::scoped_numeral_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<upolynomial::scoped_numeral_vector*> & q_coeffs,
vector<upolynomial::scoped_numeral_vector*> & F1_coeffs,
vector<upolynomial::scoped_numeral_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<a, 0<b} F1_a * F2_b
// so that F1_j * F2_0 + F2_j * F1_0 = e_j,
// i.e. F1_j * (lc_val * f2) + F2_j * f1 = e_j.
// From step 1 we have the Bezout identity s*f1 + t*f2 = 1 in Zp[x].
// Multiplying by e_j: (e_j*s)*f1 + (e_j*t)*f2 = e_j.
// Reducing to match degree constraints gives:
// F2_j = (e_j * s) mod f2
// F1_j = ((e_j * t) mod f1) / lc_val
// all in Zp[x].
// 4. Map Zp coefficients to centered representatives in (-p/2, p/2]
// and verify q = F1 * F2 over Z[x, y].
//
// Preconditions:
// - q is a bivariate polynomial in x and y, i.e. every coeff(q, y, j)
// is either zero, a constant, or univariate in x.
// - uf1_monic, uf2_monic are monic univariate polynomials in Z[x],
// stored as coefficient vectors, such that
// q(x, 0) = lc_val * uf1_monic(x) * uf2_monic(x).
// - gcd(uf1_monic, uf2_monic) = 1 over Q[x].
// - prime does not divide lc_val.
//
// Returns true if the lift succeeds, i.e. q = F1 * F2 holds over Z.
// Returns false if coprimality fails in Zp, if q is not bivariate,
// or if p is too small for the centered representative trick to work.
bool try_bivar_hensel_lift(
polynomial const * q,
var x, var y,
unsigned deg_y,
upolynomial::numeral_vector const & uf1_monic, // monic factor of q(x,0)/lc_val in Z[x], as coefficient vector
upolynomial::numeral_vector const & uf2_monic, // monic factor of q(x,0)/lc_val in Z[x], coprime to uf1_monic
numeral const & lc_val, // lc(q, x) evaluated at y=0, so q(x,0) = lc_val * uf1_monic * uf2_monic
uint64_t prime,
polynomial_ref & F1_out,
polynomial_ref & F2_out) {
@ -7019,165 +7246,27 @@ namespace polynomial {
if (d_vec.size() != 1 || !znm.is_one(d_vec[0]))
return false;
// Extract coefficients of q w.r.t. y: q_coeffs[j] = coeff(q, y, j) as upolynomials in Zp[x]
vector<upolynomial::scoped_numeral_vector*> 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<upolynomial::scoped_numeral_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<upolynomial::scoped_numeral_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, 0<a<j, 0<b<j} 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; // 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<int>(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

2
src/test/lp/lp.cpp
View file

@ -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
}