/*++ 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 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_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_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 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(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(zpk_pm, A_next_in_Z, A_pk); SASSERT(eq1); bool eq2 = are_equal_in(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(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 { /** main variable */ var m_x; public: mfactorization_combination_iterator(zp_factors const & factors, var x) : upolynomial::factorization_combination_iterator_base(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