mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-08-16 07:45:27 +00:00
Code Activity
z3/src/math/polynomial/polynomial_factorization.cpp
Christoph M. Wintersteiger f3441c6a9b tabs and indentation
2015-09-17 13:25:22 +01:00

1143 lines
43 KiB
C++
Raw Blame History

/*++
Copyright (c) 2011 Microsoft Corporation
Module Name:
polynomial_factorization.cpp
Abstract:
Implementation of polynomial factorization.
Author:
Dejan (t-dejanj) 2011-11-15
Notes:
[1] Elwyn Ralph Berlekamp. Factoring Polynomials over Finite Fields. Bell System Technical Journal,
46(8-10):1853-1859, 1967.
[2] Donald Ervin Knuth. The Art of Computer Programming, volume 2: Seminumerical Algorithms. Addison Wesley, third
edition, 1997.
[3] Henri Cohen. A Course in Computational Algebraic Number Theory. Springer Verlag, 1993.
--*/
#if 0
// disabled for reorg
#include"trace.h"
#include"util.h"
#include"polynomial_factorization.h"
#include"upolynomial_factorization_int.h"
#include"prime_generator.h"
using namespace std;
namespace polynomial {
typedef upolynomial::manager::scoped_numeral scoped_numeral;
/**
Generates a substitution of values for f -> f_univariate in order to reduce the factorization to the
univariate case.
@param f multivariate polynomial (square-free, primitive, vars(f) > 1)
@param x the variable we want to keep as the univarate one
@param f_lc the leading coefficient of f in x
@param vars the vector of all variables from f witouth x
@param size the bound to use for selecting the values, i.e. |a_i| <= size
@param a the output values corresponding to vairables in (place place for x should be ignored)
@param f_univariate the output substitution
*/
void generate_substitution_values(
polynomial_ref const & f, var x, polynomial_ref const & f_lc, var_vector const & vars, unsigned & size,
upolynomial::numeral_vector & a, upolynomial::manager upm, upolynomial::numeral_vector & f_u) {
SASSERT(a.size() == vars.size());
TRACE("polynomial::factorization", tout << "polynomial::generate_substitution_values(f = " << f << ", f_lc = " << f_lc << ")";);
// f = f_n x^n + ... + f_0, square-free and primitive
// this means
// f_lc = f_n
// since f is primitive,
// we are looking for numbers a_i such that
// (1) f_lc doesn't vanish
// (2) f_u = f(a_0, ..., a_n, x) is square-free
manager & pm = f.m();
numeral_manager & nm = pm.m();
// polynomial to use for subtituting into the lc(f)
polynomial_ref f_lc_subst(pm);
// random generator
random_gen generator;
// increase the size every once in a while (RETHINK THIS)
unsigned inc_size_c = 0;
unsigned inc_size_c_max = size;
while (true) {
// see if we should increase the size of the substitution
if ((++ inc_size_c) % inc_size_c_max == 0) {
size ++;
inc_size_c = 0;
inc_size_c_max *= 2;
}
// the head coefficient we'll substitute in
f_lc_subst = f_lc;
bool vanished = false;
for (unsigned i = 0; i < vars.size() && !vanished; ++ i) {
SASSERT(vars[i] != x);
// the value for x_i
nm.set(a[i], (int)generator(2*size+1) - (int)size);
// substitute
f_lc_subst = pm.substitute(f_lc_subst, x, a[i]);
// did it vanish
vanished = pm.is_zero(f_lc_subst);
}
if (vanished) {
// leading coefficient vanished, try again
continue;
}
// substitute into f and get the univariate one
polynomial_ref f_subst(pm);
f_subst = pm.substitute(f, vars.size(), vars.c_ptr(), a.c_ptr());
upm.to_numeral_vector(f_subst, f_u);
// if the result is not square-free we try again
if (!upm.is_square_free(f_u))
continue;
// found it, break
break;
}
}
/**
\brief Bound for the coefficients of the factorst of the multivariate polynomial f. R
Returns power of p -> p^e that covers the bound
We use the gelfond bound here:
d_i: degree of x_i in f(x1, ..., x_n)
bound = |f| * 2^(\sum d_i - (n-1)/2))
*/
void multivariate_factor_coefficient_bound(polynomial_ref const & f, var x, numeral const & p, unsigned & e, numeral & p_e, var2degree & d) {
manager & pm = f.m();
numeral_manager & nm = pm.m();
// compute g = lc(f)*f
polynomial_ref f_lc(pm), g(pm);
f_lc = pm.coeff(f, x, pm.degree(f, x));
g = pm.mul(f, f_lc);
// get the norm
scoped_numeral g_norm(nm);
pm.abs_norm(g, g_norm);
// get the variable degrees
var_vector vars;
pm.vars(f, vars);
unsigned power = 0;
for (unsigned i = 0; i < vars.size(); ++ i) {
unsigned c_d = pm.degree(g, vars[i]);
d.set_degree(vars[i], c_d + 1);
power += c_d;
}
power = power - (vars.size()-1)/2;
// compute the bound
scoped_numeral bound(nm);
nm.set(bound, 2);
nm.power(bound, power, bound);
nm.mul(g_norm, bound, bound);
// return the first power of two that is bigger than the norm
e = 1;
nm.set(p_e, p);
while (nm.lt(p_e, bound)) {
nm.mul(p_e, p_e, p_e);
e *= 2;
}
}
// check that A*S+B*T=C in zp mod ideal
bool check_solve(zp_manager & zp_pm, var2degree const & ideal,
zp_polynomial_ref const & A, zp_polynomial_ref const & B, zp_polynomial_ref const & C,
zp_polynomial_ref const & S, zp_polynomial_ref const & T) {
zp_polynomial_ref AS(zp_pm), BT(zp_pm), sum(zp_pm);
AS = zp_pm.mul(A, S); AS = zp_pm.mod_d(AS, ideal);
BT = zp_pm.mul(B, T); BT = zp_pm.mod_d(BT, ideal);
sum = zp_pm.add(AS, BT);
TRACE("polynomial::factorization::multivariate",
tout << "zp_pm = Z_" << zp_pm.m().m().to_string(zp_pm.m().p()) << endl;
tout << "ideal = " << ideal << endl;
tout << "A = " << A << endl;
tout << "B = " << B << endl;
tout << "S = " << S << endl;
tout << "T = " << T << endl;
tout << "C = " << C << endl;
tout << "sum = " << sum << endl;
);
bool result = zp_pm.eq(sum, C);
return result;
}
/**
Solve the equation A*S + B*T = C, given, AU + BV = 1, with deg(T) < deg(A)
S = U*C + tB
T = V*C - tA
we divide VC with A to get (T, t)
*/
template<typename output_manager>
void solve(zp_manager & zp_pm, var x, var2degree const & ideal,
zp_polynomial_ref const & A, zp_polynomial_ref const & U,
zp_polynomial_ref const & B, zp_polynomial_ref const & V,
zp_polynomial_ref const & C,
output_manager & out_pm,
typename output_manager::polynomial_ref & S_out,
typename output_manager::polynomial_ref & T_out) {
TRACE("polynomial::factorization::multivariate",
tout << "polynomial::solve(" << endl;
tout << "zp_pm = Z_" << zp_pm.m().m().to_string(zp_pm.m().p()) << endl;
tout << "ideal = " << ideal << endl;
tout << "A = " << A << endl;
tout << "B = " << B << endl;
tout << "U = " << U << endl;
tout << "V = " << V << endl;
tout << "C = " << C << endl;
);
// solution is S = C*U + tB, T = C*V - tA
zp_polynomial_ref CV(zp_pm);
CV = zp_pm.mul(C, V); CV = zp_pm.mod_d(CV, ideal);
zp_polynomial_ref CU(zp_pm);
CU = zp_pm.mul(C, U); CU = zp_pm.mod_d(CU, ideal);
zp_polynomial_ref t(zp_pm), T(zp_pm);
zp_pm.exact_pseudo_division_mod_d(CV, A, x, ideal, t, T);
zp_polynomial_ref tB(zp_pm);
tB = zp_pm.mul(t, B); tB = zp_pm.mod_d(tB, ideal);
zp_polynomial_ref S(zp_pm);
S = zp_pm.add(CU, tB);
SASSERT(check_solve(zp_pm, ideal, A, B, C, S, T));
// convert to the other manager
S_out = convert(zp_pm, S, out_pm);
T_out = convert(zp_pm, T, out_pm);
TRACE("polynomial::factorization::multivariate",
tout << "CU = " << CU << endl;
tout << "CV = " << CV << endl;
tout << "t = " << t << endl;
tout << "--- solution ---" << endl;
tout << "S = " << S_out << endl;
tout << "T = " << T_out << endl;
);
}
/**
A, B, U, V: multivariate polynomials in Z_p[x, ...], mod ..., also the output polynomials, y is not there
C: C = A*B in Z_p[x, ...] mod p, ...
ideal: all the vars we care about in the ideal
y, the variable we are lifting is not in ideal_vars, we will add it
A monic, A*U+B*V = 1
p is not necessary prime, it's a power of a prime
we're doing quadratic lifting here
output: added y, i.e.
* all polynomials in Z_p[x, ..., y] mod (..., y^d)
*/
void multivariate_hansel_lift_ideal(
zp_manager & zp_pm, var x,
zp_polynomial_ref const & C,
zp_polynomial_ref & A, zp_polynomial_ref & U, zp_polynomial_ref & B, zp_polynomial_ref & V,
var2degree & ideal, var y, unsigned d) {
numeral_manager & nm = zp_pm.m().m();
TRACE("polynomial::factorization::multivariate",
tout << "polynomial::multiratiate_hensel_lift_ideal" << endl;
tout << "zp_pm is Z_" << zp_pm.m().m().to_string(zp_pm.m().p()) << endl;
tout << "x = x" << x << endl;
tout << "y = x" << y << endl;
tout << "C = " << C << endl;
tout << "A = " << A << endl;
tout << "B = " << B << endl;
tout << "U = " << U << endl;
tout << "V = " << V << endl;
tout << "ideal = " << ideal << endl;
);
// constant 1
scoped_numeral one(nm);
nm.set(one, 1);
zp_polynomial_ref one_p(zp_pm);
one_p = zp_pm.mk_const(one);
SASSERT(zp_pm.degree(A, y) == 0 && zp_pm.degree(B, y) == 0 && zp_pm.degree(U, y) == 0 && zp_pm.degree(V, y) == 0);
// update the ideal, and start with y
ideal.set_degree(y, 1);
unsigned current_d = 1;
zp_polynomial_ref current_power(zp_pm);
current_power = zp_pm.mk_polynomial(y);
// lift quadratic until we are over the asked for
while (current_d < d) {
TRACE("polynomial::factorization::multivariate",
tout << "zp_pm = Z_" << nm.to_string(zp_pm.m().p()) << endl;
tout << "ideal = " << ideal << endl;
tout << "C = " << C << endl;
tout << "A = " << A << endl;
tout << "B = " << B << endl;
);
// again, classic hensel:
// since C = A*B mod (p, ideal, y^k) we know that (C - A*B) = 0 mod (p, ideal, y^k)
// f = (C - A*B) mod (p, ideal, y^k) and thus divisible by y^current_d = current_power
zp_polynomial_ref f(zp_pm);
f = zp_pm.mul(A, B);
TRACE("polynomial::factorization::multivariate",
tout << "zp_pm = Z_" << nm.to_string(zp_pm.m().p()) << endl;
tout << "ideal = " << ideal << endl;
tout << "C = " << C << endl;
tout << "A = " << A << endl;
tout << "B = " << B << endl;
tout << "f = " << f << endl;
);
f = zp_pm.sub(C, f);
f = zp_pm.exact_div(f, current_power);
f = zp_pm.mod_d(f, ideal);
TRACE("polynomial::factorization::multivariate",
tout << "A = " << A << endl;
tout << "B = " << B << endl;
);
// get the S, T, solution to A*S+B*T = f, with d(T, x) < d(A, x)
// but we know that S = U*f + Bt, T = V*f - At, so we do division
zp_polynomial_ref S(zp_pm), T(zp_pm);
solve<zp_manager>(zp_pm, x, ideal, A, U, B, V, f, zp_pm, S, T);
// now, lets lift A and B
// we want A1 = A + T*y^d, B1 = B + S*y^d with A1*B1 = C mod (ideal, y^(2*k))
// hence A*B + y^d*(A*S+B*T) = C
S = zp_pm.mul(S, current_power);
T = zp_pm.mul(T, current_power);
A = zp_pm.add(A, T);
B = zp_pm.add(B, S);
TRACE("polynomial::factorization::multivariate",
tout << "A = " << A << endl;
tout << "B = " << B << endl;
);
// again, classic quadratic hensel
// we need A*U1 + B*V1 = 1 mod (p, ideal, y^2), from above
// U1 = U + S*y^d, V1 = V + T*y^d, hence A*U + B*V + y^d(S + T) = 1 mod new ideal^2
// we know that y^d divides (1-UA-BV) so we compute f = (1-UA-BV)/y^d
// UA + VB + y^d(SA + TB) = 1 (mod ideal^2)
// SA + TB = f (mod ideal)
// we solve for S, T again, and do as above
zp_polynomial_ref UA(zp_pm), BV(zp_pm);
f = zp_pm.mk_const(one);
UA = zp_pm.mul(U, A);
BV = zp_pm.mul(V, B);
f = zp_pm.sub(f, UA); f = zp_pm.sub(f, BV);
TRACE("polynomial::factorization::multivariate",
tout << "ideal = " << ideal << endl;
tout << "current_power = " << current_power << endl;
tout << "UA = " << UA << endl;
tout << "BV = " << BV << endl;
tout << "f = " << f << endl;
tout << "x = x" << x << endl;
tout << "y = x" << y << endl;
);
f = zp_pm.exact_div(f, current_power);
f = zp_pm.mod_d(f, ideal);
// get the S, T, solution to A*S+B*T = f, with d(T, x) < d(A, x)
solve<zp_manager>(zp_pm, x, ideal, A, U, B, V, f, zp_pm, S, T);
// now, lets lift U and V
S = zp_pm.mul(S, current_power);
U = zp_pm.add(U, S);
T = zp_pm.mul(T, current_power);
V = zp_pm.add(V, T);
// lift the ideal
current_d *= 2;
current_power = zp_pm.mul(current_power, current_power);
ideal.set_degree(y, current_d);
// move, A, B, C, D into the ideal
A = zp_pm.mod_d(A, ideal);
B = zp_pm.mod_d(B, ideal);
S = zp_pm.mod_d(S, ideal);
T = zp_pm.mod_d(T, ideal);
TRACE("polynomial::factorization::multivariate",
tout << "current_d = " << d << endl;
tout << "zp_pm is Z_" << zp_pm.m().m().to_string(zp_pm.m().p()) << endl;
tout << "x = x" << x << endl;
tout << "y = x" << y << endl;
tout << "C = " << C << endl;
tout << "A = " << A << endl;
tout << "B = " << B << endl;
tout << "U = " << U << endl;
tout << "V = " << V << endl;
tout << "ideal = " << ideal << endl;
);
SASSERT(check_solve(zp_pm, ideal, A, B, one_p, U, V));
}
}
template<typename manager_to_check, typename manager_1, typename manager_2>
bool are_equal_in(
manager_to_check pm,
typename manager_1::polynomial_ref const & A,
typename manager_2::polynomial_ref const & B) {
typename manager_to_check::polynomial_ref A_pm(pm), B_pm(pm);
A_pm = convert(A.m(), A, pm);
B_pm = convert(B.m(), B, pm);
bool equal = pm.eq(A_pm, B_pm);
return equal;
}
/**
C: target multivariate polynomial mod ideal, p^e, the manager is in p^e
x: main variable
A, B, U, V: univariate polynomials in Z_p[x] such that U*A+B*V=1 mod ideal, these guys managers are in Z_p
output: A_lifted, B_lifted, A = A_lifted mod ideal
A_lifted*B_lifted = f mod x_i^d_i, p^e
*/
void multivariate_hansel_lift_zp(
manager & pm, zp_manager & zp_pm, zp_manager & zpe_pm,
zp_polynomial_ref const & C_pe, var x, unsigned e,
zp_polynomial_ref const & A_p, zp_polynomial_ref const & U_p,
zp_polynomial_ref const & B_p, zp_polynomial_ref const & V_p,
var2degree const & ideal,
zp_polynomial_ref & A_lifted, zp_polynomial_ref & B_lifted) {
TRACE("polynomial::factorization::multivariate",
tout << "polynomial::multiratiate_hensel_lift_zp:" << endl;
tout << "zp_pm = Z_" << zp_pm.m().m().to_string(zp_pm.m().p()) << endl;
tout << "zpe_pm = Z_" << zpe_pm.m().m().to_string(zpe_pm.m().p()) << endl;
tout << "x = x" << x << endl;
tout << "ideal = " << ideal << endl;
tout << "C_pe = " << C_pe << "," << endl;
tout << "A_p = " << A_p << "," << endl;
tout << "B_p = " << B_p << "," << endl;
);
// fixed zpe_pm
// upolynomial::zp_numeral_manager & zpe_nm = zpe_pm.m();
// numeral const & pe = zpe_nm.p();
// fixed zp_pm
upolynomial::zp_numeral_manager & zp_nm = zp_pm.m();
numeral const & p = zp_nm.p();
// regular numeral manager and mangager
numeral_manager & nm = zp_nm.m();
// sliding zpk_pm mod p^k
upolynomial::zp_numeral_manager zpk_nm(nm, p);
zp_manager zpk_pm(zpk_nm, &zp_pm.mm()); // in the end we copy the result over to zpe
unsigned k = 1;
upolynomial::scoped_numeral pk(nm);
nm.set(pk, zpk_nm.p());
// constant 1
scoped_numeral one(nm);
nm.set(one, 1);
zp_polynomial_ref one_p(zpk_pm);
one_p = zpk_pm.mk_const(one);
// lift until you get over the requested power of e
zp_polynomial_ref A_pk(zpk_pm), B_pk(zpk_pm), U_pk(zpk_pm), V_pk(zpk_pm);
A_pk = convert(zp_pm, A_p, zpk_pm);
B_pk = convert(zp_pm, B_p, zpk_pm);
U_pk = convert(zp_pm, U_p, zpk_pm);
V_pk = convert(zp_pm, V_p, zpk_pm);
TRACE("polynomial::factorization::multivariate",
tout << "zpk_pm = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
tout << "A_pk = " << A_pk << endl;
tout << "B_pk = " << B_pk << endl;
tout << "U_pk = " << U_pk << endl;
tout << "V_pk = " << V_pk << endl;
);
SASSERT(check_solve(zpk_pm, ideal, A_pk, B_pk, one_p, U_pk, V_pk));
while (k < e) {
// standard hensel:
// (C - AB) and is divisible by p^k, so we compute f = (C - AB)/p^k mod ideal in Z[...]
zp_polynomial_ref f_pk(zpk_pm);
polynomial_ref A_pk_in_Z(pm), B_pk_in_Z(pm), AB_in_Z(pm), f_in_Z(pm);
f_in_Z = convert(zpe_pm, C_pe, pm);
A_pk_in_Z = convert(zpk_pm, A_pk, pm);
B_pk_in_Z = convert(zpk_pm, B_pk, pm);
AB_in_Z = pm.mul(A_pk_in_Z, B_pk_in_Z);
AB_in_Z = pm.mod_d(AB_in_Z, ideal);
f_in_Z = pm.sub(f_in_Z, AB_in_Z);
f_in_Z = pm.exact_div(f_in_Z, pk);
f_in_Z = pm.mod_d(f_in_Z, ideal);
f_pk = convert(pm, f_in_Z, zpk_pm);
TRACE("polynomial::factorization::multivariate",
tout << "zpk_pm = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
tout << "f_pk = " << f_pk << endl;
);
// standard hensel we need to lift to p^(2k)
// we have U*A+V*B = 1, C = A*B, so p^k divides C - AB
// we want A1 = A + p^k*S, B1 = B + p^k*T, and also
// C - (A + p^k*S)*(B + p^k*T) = 0 mod (p^2k)
// C - A*B = p^k (T*A + S*B), i.e.
// f = (C - A*B)/p^k = (T*A + S*B), so we solve this equation in Z_p^k
polynomial_ref S_in_Z(pm), T_in_Z(pm);
solve<manager>(zpk_pm, x, ideal, A_pk, U_pk, B_pk, V_pk, f_pk, pm, T_in_Z, S_in_Z);
TRACE("polynomial::factorization::multivariate",
tout << "zpk_pm = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
tout << "S_in_Z = " << S_in_Z << endl;
tout << "T_in_Z = " << T_in_Z << endl;
);
// lift A and B to, A = A + p^k*S, B = B + p^k*T
polynomial_ref A_next_in_Z(pm), B_next_in_Z(pm);
S_in_Z = pm.mul(pk, S_in_Z);
S_in_Z = pm.mod_d(S_in_Z, ideal);
A_next_in_Z = convert(zpk_pm, A_pk, pm);
A_next_in_Z = pm.add(A_next_in_Z, S_in_Z);
T_in_Z = pm.mul(pk, T_in_Z);
T_in_Z = pm.mod_d(T_in_Z, ideal);
B_next_in_Z = convert(zpk_pm, B_pk, pm);
B_next_in_Z = pm.add(B_next_in_Z, T_in_Z);
TRACE("polynomial::factorization::multivariate",
tout << "pk = " << nm.to_string(pk) << endl;
tout << "zpk_pm = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
tout << "S_in_Z = " << S_in_Z << endl;
tout << "T_in_Z = " << T_in_Z << endl;
tout << "A_pk = " << A_pk << endl;
tout << "B_pk = " << B_pk << endl;
tout << "A_next_in_Z = " << A_next_in_Z << endl;
tout << "B_next_in_Z = " << B_next_in_Z << endl;
);
bool eq1 = are_equal_in<zp_manager, manager, zp_manager>(zpk_pm, A_next_in_Z, A_pk);
SASSERT(eq1);
bool eq2 = are_equal_in<zp_manager, manager, zp_manager>(zpk_pm, B_next_in_Z, B_pk);
SASSERT(eq2);
// again, classic quadratic hensel
// we need A*U1 + B*V1 = 1 mod p^2k, from above
// U1 = U + p^k*S, V1 = V + p^k*T, hence A*U + B*V + p^k*(S + T) = 1 mod (p^2k)
// we know that p^k divides (1-UA-BV) so we compute f = (1-UA-BV)/p^k
// UA + VB + p^k(SA + TB) = 1 (mod p^k)
// SA + TB = f (mod ideal)
// we solve for S, T again, and do as above
polynomial_ref U_pk_in_Z(pm), V_pk_in_Z(pm), UA_in_Z(pm), BV_in_Z(pm);
U_pk_in_Z = convert(zpk_pm, U_pk, pm);
V_pk_in_Z = convert(zpk_pm, V_pk, pm);
f_in_Z = pm.mk_const(one);
UA_in_Z = pm.mul(U_pk_in_Z, A_next_in_Z);
UA_in_Z = pm.mod_d(UA_in_Z, ideal);
BV_in_Z = pm.mul(V_pk_in_Z, B_next_in_Z);
BV_in_Z = pm.mod_d(BV_in_Z, ideal);
f_in_Z = pm.sub(f_in_Z, UA_in_Z);
f_in_Z = pm.sub(f_in_Z, BV_in_Z);
TRACE("polynomial::factorization::multivariate",
tout << "pk = " << nm.to_string(pk) << endl;
tout << "zpk_pm = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
tout << "U_pk_in_Z = " << U_pk_in_Z << endl;
tout << "V_pk_in_Z = " << V_pk_in_Z << endl;
tout << "UA_in_Z = " << UA_in_Z << endl;
tout << "BV_in_Z = " << BV_in_Z << endl;
tout << "f_in_Z = " << f_in_Z << endl;
);
f_in_Z = pm.exact_div(f_in_Z, pk);
f_pk = convert(pm, f_in_Z, zpk_pm);
// get the S, T, solution to A*S+B*T = f, with d(T, x) < d(A, x)
solve<manager>(zpk_pm, x, ideal, A_pk, U_pk, B_pk, V_pk, f_pk, pm, S_in_Z, T_in_Z);
TRACE("polynomial::factorization::multivariate",
tout << "pk = " << nm.to_string(pk) << endl;
tout << "zpk_pm = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
tout << "S_in_Z = " << S_in_Z << endl;
tout << "T_in_Z = " << T_in_Z << endl;
);
// go to the next zpk
scoped_numeral next_pk(nm);
nm.mul(pk, pk, next_pk);
zpk_nm.set_p(next_pk);
TRACE("polynomial::factorization::multivariate",
tout << "zp_pk = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
);
// lift U
zp_polynomial_ref S_pk(zpk_pm);
S_in_Z = pm.mul(pk, S_in_Z);
S_pk = convert(pm, S_in_Z, zpk_pm);
TRACE("polynomial::factorization::multivariate",
tout << "S_pk = " << S_pk << endl;
);
U_pk = zpk_pm.add(U_pk, S_pk);
// lift V
zp_polynomial_ref T_pk(zpk_pm);
T_in_Z = pm.mul(pk, T_in_Z);
T_pk = convert(pm, T_in_Z, zpk_pm);
TRACE("polynomial::factorization::multivariate",
tout << "T_pk = " << T_pk << endl;
);
V_pk = zpk_pm.add(V_pk, T_pk);
// lift A and B
TRACE("polynomial::factorization::multivariate",
tout << "A_pk_in_Z = " << A_pk_in_Z << endl;
tout << "B_pk_in_Z = " << B_pk_in_Z << endl;
);
A_pk = convert(pm, A_pk_in_Z, zpk_pm);
B_pk = convert(pm, B_pk_in_Z, zpk_pm);
// move to the next pk
k *= 2;
nm.set(pk, next_pk);
TRACE("polynomial::factorization::multivariate",
tout << "zp_pk = Z_" << zpk_pm.m().m().to_string(zpk_pm.m().p()) << endl;
tout << "A_pk = " << A_pk << endl;
tout << "B_pk = " << B_pk << endl;
tout << "U_pk = " << U_pk << endl;
tout << "V_pk = " << V_pk << endl;
tout << "C_pe = " << C_pe << endl;
);
SASSERT(check_solve(zpk_pm, ideal, A_pk, B_pk, one_p, U_pk, V_pk));
}
// now convert to the non-sliding zpe_manager
SASSERT(k == e);
A_lifted = convert(zpk_pm, A_pk, zpe_pm);
B_lifted = convert(zpk_pm, B_pk, zpe_pm);
}
/**
f: target multivariate polynomial mod x_i^d_i, p^e
x: main variable
all_vars: all variables (including x)
A, B, U, V: univariate polynomials in Z_p[x] such that U*A+B*V=1 from ext gcd
output: A_lifted, B_lifted d(A) = d(A_lifted), A = A_lifted mod x_i^d_i, p
A_lifted*B_lifted = f mod x_i^d_i, p^e
*/
void multivariate_hensel_lift(
manager & pm, zp_manager & zp_pm, zp_manager & zpe_pm,
zp_polynomial_ref const & f, var x, unsigned e, var_vector const & all_vars,
upolynomial::zp_manager & zp_upm,
upolynomial::numeral_vector const & U, upolynomial::numeral_vector const & A,
upolynomial::numeral_vector const & V, upolynomial::numeral_vector const & B,
var2degree & target_ideal, zp_polynomial_ref & A_lifted, zp_polynomial_ref & B_lifted) {
upolynomial::zp_numeral_manager & zp_nm = zp_upm.m();
upolynomial::numeral_manager & nm = zp_nm.m();
TRACE("polynomial::factorization::multivariate",
tout << "polynomial::multiratiate_hensel_lift(" << endl;
tout << "f = " << f << "," << endl;
tout << "x = x" << x << "," << endl;
tout << "e = " << e << "," << endl;
tout << "U = "; zp_upm.display(tout, U); tout << "," << endl;
tout << "A = "; zp_upm.display(tout, A); tout << "," << endl;
tout << "V = "; zp_upm.display(tout, V); tout << "," << endl;
tout << "B = "; zp_upm.display(tout, B); tout << "," << endl;
tout << "target_ideal = " << target_ideal << "," << endl;
tout << "p = " << nm.to_string(zp_pm.m().p()) << endl;
tout << "pe = " << nm.to_string(zpe_pm.m().p()) << endl;
);
// multivariate versions of A, B, U, V that we keep lifting over ideal x_i^d_i
zp_polynomial_ref A_m_p(zp_pm), B_m_p(zp_pm), U_m_p(zp_pm), V_m_p(zp_pm);
A_m_p = zp_pm.to_polynomial(A, x);
B_m_p = zp_pm.to_polynomial(B, x);
U_m_p = zp_pm.to_polynomial(U, x);
V_m_p = zp_pm.to_polynomial(V, x);
TRACE("polynomial::factorization::multivariate",
tout << "A_m_p = " << A_m_p << endl;
tout << "B_m_p = " << B_m_p << endl;
tout << "U_m_p = " << U_m_p << endl;
tout << "V_m_p = " << V_m_p << endl;
);
// the the target in Z_p[...]
zp_polynomial_ref C_m_p(zp_pm);
C_m_p = convert(zpe_pm, f, zp_pm);
// lift each variable individually
var2degree lifted_ideal;
unsigned_vector lifted_degs;
for (unsigned i = 0; i < all_vars.size(); ++ i) {
if (all_vars[i] == x) {
// skip the main variable
continue;
}
// current variable and degree we are lifting to, y^(d_y), at least
var y = all_vars[i];
// lift to y^(d_y)
multivariate_hansel_lift_ideal(zp_pm, x, C_m_p, A_m_p, U_m_p, B_m_p, V_m_p, lifted_ideal, y, target_ideal.degree(y));
}
TRACE("polynomial::factorization::multivariate",
tout << "A_m_p = " << A_m_p << endl;
tout << "B_m_p = " << B_m_p << endl;
tout << "U_m_p = " << U_m_p << endl;
tout << "V_m_p = " << V_m_p << endl;
tout << "lifted_ideal = " << lifted_ideal << endl;
);
// now lift it all to p^e
multivariate_hansel_lift_zp(pm, zp_pm, zpe_pm, f, x, e, A_m_p, U_m_p, B_m_p, V_m_p, lifted_ideal, A_lifted, B_lifted);
}
/**
f: multivariate polynomial
x: main variable
all_vars: all variables (including x)
f_u: f mod x_1, ..., x_n (excluding mod x), i.e. this is f(0, x), f_u is square_free
f_u_zp_factors: monic factors of f_u (mod p), pairvise gcd = 1
we're lifting the factors to mod x_1^d_1, ..., x_n&d_n (excliding x), mod p^e
i.e. such that f congruent to the new factors. output goes to f_zpe factors.
*/
void multivariate_hensel_lift(
manager & pm, zp_manager & zp_pm, zp_manager & zpe_pm,
polynomial_ref const & f, var x, unsigned e, var_vector const & all_vars,
upolynomial::manager & upm, upolynomial::numeral_vector const & f_u,
upolynomial::zp_factors const & f_u_zp_factors,
var2degree & target_ideal,
zp_factors & f_zpe_factors) {
SASSERT(f_u_zp_factors.distinct_factors() > 1);
TRACE("polynomial::factorization::multivariate",
tout << "polynomial::multivariate_hensel_lift(" << endl;
tout << "f = " << f << "," << endl;
tout << "x = x" << x << "," << endl;
tout << "e = " << e << "," << endl;
tout << "f_u = "; upm.display(tout, f_u); tout << "," << endl;
tout << "f_u_zp_factors" << f_u_zp_factors << "," << endl;
tout << "target_ideal = " << target_ideal << "," << endl;
tout << "f_zpe_factors = " << f_zpe_factors << ")" << endl;
);
// managers and all
numeral_manager & nm = pm.m();
upolynomial::zp_manager & zp_upm = f_u_zp_factors.upm();
// upolynomial::zp_numeral_manager & zp_nm = zp_upm.m();
upolynomial::zp_numeral_manager & zpe_nm = zpe_pm.m();
upolynomial::zp_manager zpe_upm(zpe_nm);
// the targed product we want (mod x_i^d_i, mod p^e)
zp_polynomial_ref f_target_zpe(zpe_pm);
f_target_zpe = convert(pm, f, zpe_pm);
f_target_zpe = zpe_pm.mod_d(f_target_zpe, target_ideal);
TRACE("polynomial::factorization::multivariate",
tout << "target_ideal = " << target_ideal << endl;
tout << "f_target_zpe = " << f_target_zpe << endl;
);
// we do the product by doing individual lifting like in the univarate case
zp_polynomial_ref B(zp_pm), C_p(zp_pm);
zp_polynomial_ref A_lifted(zpe_pm), B_lifted(zpe_pm);
upolynomial::scoped_numeral_vector B_u(nm), C_u(nm), tmp_u(nm);
upolynomial::scoped_numeral_vector U(nm), V(nm);
for (int i = 0, i_end = f_u_zp_factors.distinct_factors() - 1; i < i_end; ++ i)
{
// get the univarate ones to lift now
upolynomial::numeral_vector const & A_u = f_u_zp_factors[i];
// current starting product is f_target_zpe(0, x) in *Z_p*
zp_upm.to_numeral_vector(f_target_zpe, x, C_u);
// we get the rest into B (mod p)
zp_upm.exact_div(C_u, A_u, B_u);
TRACE("polynomial::factorization::multivariate",
tout << "p = " << nm.to_string(zp_upm.m().p()) << endl;
tout << "f_target_zpe = " << f_target_zpe << endl;
tout << "A_u = "; upm.display(tout, A_u); tout << endl;
tout << "B_u = "; upm.display(tout, B_u); tout << endl;
tout << "C_u = "; upm.display(tout, C_u); tout << endl;
);
// and get the U, V, such that A*U+B*V = 1
zp_upm.ext_gcd(A_u, B_u, U, V, tmp_u);
TRACE("polynomial::factorization::multivariate",
tout << "U = "; upm.display(tout, U); tout << endl;
tout << "V = "; upm.display(tout, V); tout << endl;
tout << "gcd = "; upm.display(tout, tmp_u); tout << endl;
);
// do the lifting for this pair
multivariate_hensel_lift(pm, zp_pm, zpe_pm, f_target_zpe, x, e, all_vars, zp_upm, U, A_u, V, B_u, target_ideal, A_lifted, B_lifted);
// add the lifted A to the output
f_zpe_factors.push_back(A_lifted, 1);
// move to the new target by dividing with the lifted A
f_target_zpe = zpe_pm.exact_div(f_target_zpe, A_lifted);
}
// add the last f_target
f_zpe_factors.push_back(f_target_zpe, 1);
}
class mfactorization_combination_iterator : public upolynomial::factorization_combination_iterator_base<zp_factors> {
/** main variable */
var m_x;
public:
mfactorization_combination_iterator(zp_factors const & factors, var x)
: upolynomial::factorization_combination_iterator_base<zp_factors>(factors)
{}
/**
\brief Filter the ones not in the degree set.
*/
bool filter_current() const {
return false;
}
/**
\brief Returns the degree of the current selection.
*/
unsigned current_degree() const {
unsigned degree = 0;
zp_manager & pm = m_factors.pm();
for (unsigned i = 0; i < left_size(); ++ i) {
degree += pm.degree(m_factors[m_current[i]], m_x);
}
return degree;
}
void left(zp_polynomial_ref & out) const {
SASSERT(m_current_size > 0);
zp_manager & zp_pm = m_factors.pm();
out = m_factors[m_current[0]];
for (int i = 1; i < m_current_size; ++ i) {
out = zp_pm.mul(out, m_factors[m_current[i]]);
}
}
void get_left_tail_coeff(numeral const & m, numeral & out) {
zp_manager & zp_pm = m_factors.pm();
upolynomial::zp_numeral_manager & zp_nm = zp_pm.m();
zp_nm.set(out, m);
for (int i = 0; i < m_current_size; ++ i) {
zp_nm.mul(out, zp_pm.numeral_tc(m_factors[m_current[i]]), out);
}
}
void get_right_tail_coeff(numeral const & m, numeral & out) {
zp_manager & zp_pm = m_factors.pm();
upolynomial::zp_numeral_manager & zp_nm = zp_pm.m();
zp_nm.set(out, m);
unsigned current = 0;
unsigned selection_i = 0;
// selection is ordered, so we just take the ones in between that are not disable
while (current < m_factors.distinct_factors()) {
if (!m_enabled[current]) {
// by skipping the disabled we never skip a selected one
current ++;
} else {
if (selection_i >= m_current.size() || (int) current < m_current[selection_i]) {
SASSERT(m_factors.get_degree(current) == 1);
zp_nm.mul(out, zp_pm.numeral_tc(m_factors[current]), out);
current ++;
} else {
current ++;
selection_i ++;
}
}
}
}
void right(zp_polynomial_ref & out) const {
SASSERT(m_current_size > 0);
zp_manager & zp_pm = m_factors.pm();
upolynomial::zp_numeral_manager & zp_nm = zp_pm.m();
unsigned current = 0;
unsigned selection_i = 0;
numeral one;
zp_nm.set(one, 1);
out = zp_pm.mk_const(one);
// selection is ordered, so we just take the ones in between that are not disable
while (current < m_factors.distinct_factors()) {
if (!m_enabled[current]) {
// by skipping the disabled we never skip a selected one
current ++;
} else {
if (selection_i >= m_current.size() || (int) current < m_current[selection_i]) {
SASSERT(m_factors.get_degree(current) == 1);
out = zp_pm.mul(out, m_factors[current]);
current ++;
} else {
current ++;
selection_i ++;
}
}
}
}
};
// the multivariate factorization
bool factor_square_free_primitive(polynomial_ref const & f, factors & f_factors) {
TRACE("polynomial::factorization", tout << "polynomial::factor_square_free_primitive(f = " << f << ", factors = " << f_factors << ")" << endl;);
manager & pm = f.m();
numeral_manager & nm = pm.m();
// to start with, maybe this should be part of input
var x = pm.max_var(f);
// get all the variables
var_vector vars, vars_no_x;
pm.vars(f, vars);
for(unsigned i = 0; i < vars.size(); ++ i) {
if (vars[i] != x) {
vars_no_x.push_back(vars[i]);
}
}
SASSERT(vars.size() > 1);
// degree of the main variable
unsigned x_degree = pm.degree(f, x);
// the leading coefficient
polynomial_ref f_lc(pm);
f_lc = pm.coeff(f, x, x_degree);
// the vector of values we substitute
upolynomial::scoped_numeral_vector a(nm);
// the univariate polynomial
upolynomial::manager upm(nm);
upolynomial::scoped_numeral_vector f_u(upm);
// generate the values to substitute and substitute them to get f_u(x) = f(a, x), the univariate version of f
unsigned size = 1;
a.resize(vars_no_x.size());
for (unsigned i = 0; i < a.size(); ++ i) { nm.reset(a[i]); }
generate_substitution_values(f, x, f_lc, vars_no_x, size, a, upm, f_u);
TRACE("polynomial::factorization::multivariate",
tout << "f_u = "; upm.display(tout, f_u); tout << endl;
tout << "substitution:" << endl;
for (unsigned i = 0; i < vars_no_x.size(); ++ i) {
tout << "x" << vars[i] << " -> " << nm.to_string(a[i]) << endl;
});
// the primitive part of f_u
scoped_numeral f_u_lc(nm);
upolynomial::scoped_numeral_vector f_u_pp(nm);
upm.get_primitive_and_content(f_u, f_u_pp, f_u_lc);
TRACE("polynomial::factorization::multivariate",
tout << "f_u_lc" << nm.to_string(f_u_lc) << endl;
tout << "f_u_pp = "; upm.display(tout, f_u_pp); tout << endl;
);
// factor the univariate one
upolynomial::factors factors_u(upm);
upolynomial::factor_square_free(upm, f_u, factors_u);
TRACE("polynomial::factorization::multivariate",
tout << "factors_u = " << factors_u << endl;
);
// if there is no univariate factors, it's irreducible
if (factors_u.distinct_factors() == 1) {
f_factors.push_back(f, 1);
return false;
}
// translate f with a, so that we work modulo x_i^k and not (x_i - a_i)^k
polynomial_ref f_t(pm);
// Do the translation, we must have that a[x] = 0
pm.translate(f, vars_no_x, a, f_t);
TRACE("polynomial::factorization::multivariate",
tout << "f_t = " << f_t << endl;
);
// the zp manager stuff, we'll be changing the base
upolynomial::zp_numeral_manager zp_nm(nm, 2);
upolynomial::zp_manager zp_upm(zp_nm);
// get the first prime number p that keeps f square-free
scoped_numeral p(nm);
prime_iterator prime_it;
upolynomial::scoped_numeral_vector f_u_pp_zp(nm);
do {
// create Z_p with the next prime
nm.set(p, prime_it.next());
// translate to Zp[x]
zp_nm.set_p(p);
upolynomial::to_zp_manager(zp_upm, f_u_pp, f_u_pp_zp);
} while (!zp_upm.is_square_free(f_u_pp_zp));
TRACE("polynomial::factorization::multivariate",
tout << "p = " << p << endl;
tout << "f_t = " << f_t << endl;
);
// convert factors of f to factors modulo p (monic ones)
upolynomial::zp_factors factors_u_zp(zp_upm);
upolynomial::scoped_numeral_vector current_factor(nm);
for (unsigned i = 0; i < factors_u.distinct_factors(); ++ i) {
upolynomial::to_zp_manager(zp_upm, factors_u[i], current_factor);
zp_upm.mk_monic(current_factor);
factors_u_zp.push_back_swap(current_factor, 1);
}
TRACE("polynomial::factorization::multivariate",
tout << "factors_u_zp = " << factors_u_zp << endl;
);
// compute factor coefficient bound (pover p^e) of the translated f with the leading coefficient added, i.e.
// f_t*lc(f_t, x) = f_t*lc(f, x)
unsigned e;
scoped_numeral p_e(nm);
var2degree target_ideal;
upolynomial::scoped_numeral f_t_lc(nm);
nm.set(f_t_lc, pm.numeral_lc(f_t, x));
polynomial_ref f_t_with_lc(pm);
f_t_with_lc = pm.mul(f_t_lc, f_t);
multivariate_factor_coefficient_bound(f_t_with_lc, x, p, e, p_e, target_ideal);
TRACE("polynomial::factorization::multivariate",
tout << "target_ideal = " << target_ideal << endl;
);
// do the multivariate lifting of the translated one f_t
upolynomial::zp_numeral_manager zpe_nm(nm, p_e);
zp_manager zpe_pm(zpe_nm, &pm.mm());
zp_manager zp_pm(zp_nm, &pm.mm());
zp_factors factors_zpe(zpe_pm);
multivariate_hensel_lift(pm, zp_pm, zpe_pm, f_t, x, e, vars, upm, f_u_pp_zp, factors_u_zp, target_ideal, factors_zpe);
TRACE("polynomial::factorization::multivariate",
tout << "factors_zpe = " << factors_zpe << endl;
);
// try the factors from the lifted combinations
factors f_t_factors(pm);
bool remove = false;
mfactorization_combination_iterator it(factors_zpe, x);
unsigned max_degre = pm.degree(f_t, x) / 2;
zp_polynomial_ref zpe_trial_factor(zpe_pm);
while (it.next(remove)) {
//
// our bound ensures we can extract the right factors of degree at most 1/2 of the original
// so, if out trial factor has degree bigger than 1/2, we need to take the rest of the factors
// but, if we take the rest and it works, it doesn't meen that the rest is factorized, so we still take out
// the original factor
//
bool using_left = it.current_degree() <= max_degre;
if (using_left) {
// do a quick check first
it.left(zpe_trial_factor);
} else {
// do a quick check first
it.right(zpe_trial_factor);
}
// add the lc(f_pp) to the trial divisor
zpe_trial_factor = zpe_pm.mul(f_t_lc, zpe_trial_factor);
polynomial_ref trial_factor(pm), trial_factor_quo(pm);
trial_factor = convert(zpe_pm, zpe_trial_factor, pm);
bool true_factor = true;
trial_factor_quo = pm.exact_div(f_t_with_lc, trial_factor);
// if division is precise we have a factor
if (true_factor) {
if (!using_left) {
// as noted above, we still use the original factor
trial_factor.swap(trial_factor_quo);
}
// We need to get the content out of the factor
upolynomial::scoped_numeral trial_factor_cont(nm);
pm.int_content(f_t, trial_factor_cont);
trial_factor = pm.exact_div(trial_factor, trial_factor_cont);
// add the factor
f_t_factors.push_back(trial_factor, 1);
// we continue with the int-primitive quotient (with the lc added back)
// but we also have to keep lc(f_t)*f_t
pm.int_content(trial_factor_quo, f_t_lc); // content
trial_factor_quo = pm.exact_div(trial_factor_quo, f_t_lc);
nm.set(f_t_lc, pm.numeral_lc(trial_factor_quo, x));
f_t = pm.mul(f_t_lc, trial_factor_quo);
// but we also remove it from the iterator
remove = true;
} else {
// don't remove this combination
remove = false;
}
}
// translate the factor back
for (unsigned i = 0; i < a.size(); ++ i) {
nm.neg(a[i]);
}
for (unsigned i = 0; i < f_t_factors.distinct_factors(); ++ i) {
polynomial_ref factor(pm);
pm.translate(f_t_factors[i], vars_no_x, a, factor);
f_factors.push_back(factor, 1);
}
TRACE("polynomial::factorization", tout << "polynomial::factor_square_free_primitive(f = " << f << ") => " << f_factors << endl;);
return true;
}
}; // end polynomial namespace
#endif
View git blame Copy permalink