Coverage for sympy/polys/factortools.py : 67%

"""Polynomial factorization routines in characteristic zero. """ 

from __future__ import print_function, division 

from sympy.polys.galoistools import ( 

gf_from_int_poly, gf_to_int_poly, 

gf_lshift, gf_add_mul, gf_mul, 

gf_div, gf_rem, 

gf_gcdex, 

gf_sqf_p, 

gf_factor_sqf, gf_factor) 

from sympy.polys.densebasic import ( 

dup_LC, dmp_LC, dmp_ground_LC, 

dup_TC, 

dup_convert, dmp_convert, 

dup_degree, dmp_degree, 

dmp_degree_in, dmp_degree_list, 

dmp_from_dict, 

dmp_zero_p, 

dmp_one, 

dmp_nest, dmp_raise, 

dup_strip, 

dmp_ground, 

dup_inflate, 

dmp_exclude, dmp_include, 

dmp_inject, dmp_eject, 

dup_terms_gcd, dmp_terms_gcd) 

from sympy.polys.densearith import ( 

dup_neg, dmp_neg, 

dup_add, dmp_add, 

dup_sub, dmp_sub, 

dup_mul, dmp_mul, 

dup_sqr, 

dmp_pow, 

dup_div, dmp_div, 

dup_quo, dmp_quo, 

dmp_expand, 

dmp_add_mul, 

dup_sub_mul, dmp_sub_mul, 

dup_lshift, 

dup_max_norm, dmp_max_norm, 

dup_l1_norm, 

dup_mul_ground, dmp_mul_ground, 

dup_quo_ground, dmp_quo_ground) 

from sympy.polys.densetools import ( 

dup_clear_denoms, dmp_clear_denoms, 

dup_trunc, dmp_ground_trunc, 

dup_content, 

dup_monic, dmp_ground_monic, 

dup_primitive, dmp_ground_primitive, 

dmp_eval_tail, 

dmp_eval_in, dmp_diff_eval_in, 

dmp_compose, 

dup_shift, dup_mirror) 

from sympy.polys.euclidtools import ( 

dmp_primitive, 

dup_inner_gcd, dmp_inner_gcd) 

from sympy.polys.sqfreetools import ( 

dup_sqf_p, 

dup_sqf_norm, dmp_sqf_norm, 

dup_sqf_part, dmp_sqf_part) 

from sympy.polys.polyutils import _sort_factors 

from sympy.polys.polyconfig import query 

from sympy.polys.polyerrors import ( 

ExtraneousFactors, DomainError, CoercionFailed, EvaluationFailed) 

from sympy.ntheory import nextprime, isprime, factorint 

from sympy.utilities import subsets 

from math import ceil as _ceil, log as _log 

from sympy.core.compatibility import range 

def dup_trial_division(f, factors, K): 

"""Determine multiplicities of factors using trial division. """ 

result = [] 

for factor in factors: 

k = 0 

while True: 

q, r = dup_div(f, factor, K) 

if not r: 

f, k = q, k + 1 

else: 

break 

result.append((factor, k)) 

return _sort_factors(result) 

def dmp_trial_division(f, factors, u, K): 

"""Determine multiplicities of factors using trial division. """ 

result = [] 

for factor in factors: 

k = 0 

while True: 

q, r = dmp_div(f, factor, u, K) 

if dmp_zero_p(r, u): 

f, k = q, k + 1 

else: 

break 

result.append((factor, k)) 

return _sort_factors(result) 

def dup_zz_mignotte_bound(f, K): 

"""Mignotte bound for univariate polynomials in `K[x]`. """ 

a = dup_max_norm(f, K) 

b = abs(dup_LC(f, K)) 

n = dup_degree(f) 

return K.sqrt(K(n + 1))*2**n*a*b 

def dmp_zz_mignotte_bound(f, u, K): 

"""Mignotte bound for multivariate polynomials in `K[X]`. """ 

a = dmp_max_norm(f, u, K) 

b = abs(dmp_ground_LC(f, u, K)) 

n = sum(dmp_degree_list(f, u)) 

return K.sqrt(K(n + 1))*2**n*a*b 

def dup_zz_hensel_step(m, f, g, h, s, t, K): 

""" 

One step in Hensel lifting in `Z[x]`. 

Given positive integer `m` and `Z[x]` polynomials `f`, `g`, `h`, `s` 

and `t` such that:: 

f == g*h (mod m) 

s*g + t*h == 1 (mod m) 

lc(f) is not a zero divisor (mod m) 

lc(h) == 1 

deg(f) == deg(g) + deg(h) 

deg(s) < deg(h) 

deg(t) < deg(g) 

returns polynomials `G`, `H`, `S` and `T`, such that:: 

f == G*H (mod m**2) 

S*G + T**H == 1 (mod m**2) 

References 

========== 

1. [Gathen99]_ 

""" 

M = m**2 

e = dup_sub_mul(f, g, h, K) 

e = dup_trunc(e, M, K) 

q, r = dup_div(dup_mul(s, e, K), h, K) 

q = dup_trunc(q, M, K) 

r = dup_trunc(r, M, K) 

u = dup_add(dup_mul(t, e, K), dup_mul(q, g, K), K) 

G = dup_trunc(dup_add(g, u, K), M, K) 

H = dup_trunc(dup_add(h, r, K), M, K) 

u = dup_add(dup_mul(s, G, K), dup_mul(t, H, K), K) 

b = dup_trunc(dup_sub(u, [K.one], K), M, K) 

c, d = dup_div(dup_mul(s, b, K), H, K) 

c = dup_trunc(c, M, K) 

d = dup_trunc(d, M, K) 

u = dup_add(dup_mul(t, b, K), dup_mul(c, G, K), K) 

S = dup_trunc(dup_sub(s, d, K), M, K) 

T = dup_trunc(dup_sub(t, u, K), M, K) 

return G, H, S, T 

def dup_zz_hensel_lift(p, f, f_list, l, K): 

""" 

Multifactor Hensel lifting in `Z[x]`. 

Given a prime `p`, polynomial `f` over `Z[x]` such that `lc(f)` 

is a unit modulo `p`, monic pair-wise coprime polynomials `f_i` 

over `Z[x]` satisfying:: 

f = lc(f) f_1 ... f_r (mod p) 

and a positive integer `l`, returns a list of monic polynomials 

`F_1`, `F_2`, ..., `F_r` satisfying:: 

f = lc(f) F_1 ... F_r (mod p**l) 

F_i = f_i (mod p), i = 1..r 

References 

========== 

1. [Gathen99]_ 

""" 

r = len(f_list) 

lc = dup_LC(f, K) 

if r == 1: 

F = dup_mul_ground(f, K.gcdex(lc, p**l)[0], K) 

return [ dup_trunc(F, p**l, K) ] 

m = p 

k = r // 2 

d = int(_ceil(_log(l, 2))) 

g = gf_from_int_poly([lc], p) 

for f_i in f_list[:k]: 

g = gf_mul(g, gf_from_int_poly(f_i, p), p, K) 

h = gf_from_int_poly(f_list[k], p) 

for f_i in f_list[k + 1:]: 

h = gf_mul(h, gf_from_int_poly(f_i, p), p, K) 

s, t, _ = gf_gcdex(g, h, p, K) 

g = gf_to_int_poly(g, p) 

h = gf_to_int_poly(h, p) 

s = gf_to_int_poly(s, p) 

t = gf_to_int_poly(t, p) 

for _ in range(1, d + 1): 

(g, h, s, t), m = dup_zz_hensel_step(m, f, g, h, s, t, K), m**2 

return dup_zz_hensel_lift(p, g, f_list[:k], l, K) \ 

+ dup_zz_hensel_lift(p, h, f_list[k:], l, K) 

def _test_pl(fc, q, pl): 

if q > pl // 2: 

q = q - pl 

if not q: 

return True 

return fc % q == 0 

def dup_zz_zassenhaus(f, K): 

"""Factor primitive square-free polynomials in `Z[x]`. """ 

n = dup_degree(f) 

if n == 1: 

return [f] 

fc = f[-1] 

A = dup_max_norm(f, K) 

b = dup_LC(f, K) 

B = int(abs(K.sqrt(K(n + 1))*2**n*A*b)) 

C = int((n + 1)**(2*n)*A**(2*n - 1)) 

gamma = int(_ceil(2*_log(C, 2))) 

bound = int(2*gamma*_log(gamma)) 

a = [] 

# choose a prime number `p` such that `f` be square free in Z_p 

# if there are many factors in Z_p, choose among a few different `p` 

# the one with fewer factors 

for px in range(3, bound + 1): 

if not isprime(px) or b % px == 0: 

continue 

px = K.convert(px) 

F = gf_from_int_poly(f, px) 

if not gf_sqf_p(F, px, K): 

continue 

fsqfx = gf_factor_sqf(F, px, K)[1] 

a.append((px, fsqfx)) 

if len(fsqfx) < 15 or len(a) > 4: 

break 

p, fsqf = min(a, key=lambda x: len(x[1])) 

l = int(_ceil(_log(2*B + 1, p))) 

modular = [gf_to_int_poly(ff, p) for ff in fsqf] 

g = dup_zz_hensel_lift(p, f, modular, l, K) 

sorted_T = range(len(g)) 

T = set(sorted_T) 

factors, s = [], 1 

pl = p**l 

while 2*s <= len(T): 

for S in subsets(sorted_T, s): 

# lift the constant coefficient of the product `G` of the factors 

# in the subset `S`; if it is does not divide `fc`, `G` does 

# not divide the input polynomial 

if b == 1: 

q = 1 

for i in S: 

q = q*g[i][-1] 

q = q % pl 

if not _test_pl(fc, q, pl): 

continue 

else: 

G = [b] 

for i in S: 

G = dup_mul(G, g[i], K) 

G = dup_trunc(G, pl, K) 

G = dup_primitive(G, K)[1] 

q = G[-1] 

if q and fc % q != 0: 

continue 

H = [b] 

S = set(S) 

T_S = T - S 

if b == 1: 

G = [b] 

for i in S: 

G = dup_mul(G, g[i], K) 

G = dup_trunc(G, pl, K) 

for i in T_S: 

H = dup_mul(H, g[i], K) 

H = dup_trunc(H, pl, K) 

G_norm = dup_l1_norm(G, K) 

H_norm = dup_l1_norm(H, K) 

if G_norm*H_norm <= B: 

T = T_S 

sorted_T = [i for i in sorted_T if i not in S] 

G = dup_primitive(G, K)[1] 

f = dup_primitive(H, K)[1] 

factors.append(G) 

b = dup_LC(f, K) 

break 

else: 

s += 1 

return factors + [f] 

def dup_zz_irreducible_p(f, K): 

"""Test irreducibility using Eisenstein's criterion. """ 

lc = dup_LC(f, K) 

tc = dup_TC(f, K) 

e_fc = dup_content(f[1:], K) 

if e_fc: 

e_ff = factorint(int(e_fc)) 

for p in e_ff.keys(): 

if (lc % p) and (tc % p**2): 

return True 

def dup_cyclotomic_p(f, K, irreducible=False): 

""" 

Efficiently test if ``f`` is a cyclotomic polnomial. 

Examples 

======== 

>>> from sympy.polys import ring, ZZ 

>>> R, x = ring("x", ZZ) 

>>> f = x**16 + x**14 - x**10 + x**8 - x**6 + x**2 + 1 

>>> R.dup_cyclotomic_p(f) 

False 

>>> g = x**16 + x**14 - x**10 - x**8 - x**6 + x**2 + 1 

>>> R.dup_cyclotomic_p(g) 

True 

""" 

if K.is_QQ: 

try: 

K0, K = K, K.get_ring() 

f = dup_convert(f, K0, K) 

except CoercionFailed: 

return False 

elif not K.is_ZZ: 

return False 

lc = dup_LC(f, K) 

tc = dup_TC(f, K) 

if lc != 1 or (tc != -1 and tc != 1): 

return False 

if not irreducible: 

coeff, factors = dup_factor_list(f, K) 

if coeff != K.one or factors != [(f, 1)]: 

return False 

n = dup_degree(f) 

g, h = [], [] 

for i in range(n, -1, -2): 

g.insert(0, f[i]) 

for i in range(n - 1, -1, -2): 

h.insert(0, f[i]) 

g = dup_sqr(dup_strip(g), K) 

h = dup_sqr(dup_strip(h), K) 

F = dup_sub(g, dup_lshift(h, 1, K), K) 

if K.is_negative(dup_LC(F, K)): 

F = dup_neg(F, K) 

if F == f: 

return True 

g = dup_mirror(f, K) 

if K.is_negative(dup_LC(g, K)): 

g = dup_neg(g, K) 

if F == g and dup_cyclotomic_p(g, K): 

return True 

G = dup_sqf_part(F, K) 

if dup_sqr(G, K) == F and dup_cyclotomic_p(G, K): 

return True 

return False 

def dup_zz_cyclotomic_poly(n, K): 

"""Efficiently generate n-th cyclotomic polnomial. """ 

h = [K.one, -K.one] 

for p, k in factorint(n).items(): 

h = dup_quo(dup_inflate(h, p, K), h, K) 

h = dup_inflate(h, p**(k - 1), K) 

return h 

def _dup_cyclotomic_decompose(n, K): 

H = [[K.one, -K.one]] 

for p, k in factorint(n).items(): 

Q = [ dup_quo(dup_inflate(h, p, K), h, K) for h in H ] 

H.extend(Q) 

for i in range(1, k): 

Q = [ dup_inflate(q, p, K) for q in Q ] 

H.extend(Q) 

return H 

def dup_zz_cyclotomic_factor(f, K): 

""" 

Efficiently factor polynomials `x**n - 1` and `x**n + 1` in `Z[x]`. 

Given a univariate polynomial `f` in `Z[x]` returns a list of factors 

of `f`, provided that `f` is in the form `x**n - 1` or `x**n + 1` for 

`n >= 1`. Otherwise returns None. 

Factorization is performed using using cyclotomic decomposition of `f`, 

which makes this method much faster that any other direct factorization 

approach (e.g. Zassenhaus's). 

References 

========== 

1. [Weisstein09]_ 

""" 

lc_f, tc_f = dup_LC(f, K), dup_TC(f, K) 

if dup_degree(f) <= 0: 

return None 

if lc_f != 1 or tc_f not in [-1, 1]: 

return None 

if any(bool(cf) for cf in f[1:-1]): 

return None 

n = dup_degree(f) 

F = _dup_cyclotomic_decompose(n, K) 

if not K.is_one(tc_f): 

return F 

else: 

H = [] 

for h in _dup_cyclotomic_decompose(2*n, K): 

if h not in F: 

H.append(h) 

return H 

def dup_zz_factor_sqf(f, K): 

"""Factor square-free (non-primitive) polyomials in `Z[x]`. """ 

cont, g = dup_primitive(f, K) 

n = dup_degree(g) 

if dup_LC(g, K) < 0: 

cont, g = -cont, dup_neg(g, K) 

if n <= 0: 

return cont, [] 

elif n == 1: 

return cont, [g] 

if query('USE_IRREDUCIBLE_IN_FACTOR'): 

if dup_zz_irreducible_p(g, K): 

return cont, [g] 

factors = None 

if query('USE_CYCLOTOMIC_FACTOR'): 

factors = dup_zz_cyclotomic_factor(g, K) 

if factors is None: 

factors = dup_zz_zassenhaus(g, K) 

return cont, _sort_factors(factors, multiple=False) 

def dup_zz_factor(f, K): 

""" 

Factor (non square-free) polynomials in `Z[x]`. 

Given a univariate polynomial `f` in `Z[x]` computes its complete 

factorization `f_1, ..., f_n` into irreducibles over integers:: 

f = content(f) f_1**k_1 ... f_n**k_n 

The factorization is computed by reducing the input polynomial 

into a primitive square-free polynomial and factoring it using 

Zassenhaus algorithm. Trial division is used to recover the 

multiplicities of factors. 

The result is returned as a tuple consisting of:: 

(content(f), [(f_1, k_1), ..., (f_n, k_n)) 

Consider polynomial `f = 2*x**4 - 2`:: 

>>> from sympy.polys import ring, ZZ 

>>> R, x = ring("x", ZZ) 

>>> R.dup_zz_factor(2*x**4 - 2) 

(2, [(x - 1, 1), (x + 1, 1), (x**2 + 1, 1)]) 

In result we got the following factorization:: 

f = 2 (x - 1) (x + 1) (x**2 + 1) 

Note that this is a complete factorization over integers, 

however over Gaussian integers we can factor the last term. 

By default, polynomials `x**n - 1` and `x**n + 1` are factored 

using cyclotomic decomposition to speedup computations. To 

disable this behaviour set cyclotomic=False. 

References 

========== 

1. [Gathen99]_ 

""" 

cont, g = dup_primitive(f, K) 

n = dup_degree(g) 

if dup_LC(g, K) < 0: 

cont, g = -cont, dup_neg(g, K) 

if n <= 0: 

return cont, [] 

elif n == 1: 

return cont, [(g, 1)] 

if query('USE_IRREDUCIBLE_IN_FACTOR'): 

if dup_zz_irreducible_p(g, K): 

return cont, [(g, 1)] 

g = dup_sqf_part(g, K) 

H = None 

if query('USE_CYCLOTOMIC_FACTOR'): 

H = dup_zz_cyclotomic_factor(g, K) 

if H is None: 

H = dup_zz_zassenhaus(g, K) 

factors = dup_trial_division(f, H, K) 

return cont, factors 

def dmp_zz_wang_non_divisors(E, cs, ct, K): 

"""Wang/EEZ: Compute a set of valid divisors. """ 

result = [ cs*ct ] 

for q in E: 

q = abs(q) 

for r in reversed(result): 

while r != 1: 

r = K.gcd(r, q) 

q = q // r 

if K.is_one(q): 

return None 

result.append(q) 

return result[1:] 

def dmp_zz_wang_test_points(f, T, ct, A, u, K): 

"""Wang/EEZ: Test evaluation points for suitability. """ 

if not dmp_eval_tail(dmp_LC(f, K), A, u - 1, K): 

raise EvaluationFailed('no luck') 

g = dmp_eval_tail(f, A, u, K) 

if not dup_sqf_p(g, K): 

raise EvaluationFailed('no luck') 

c, h = dup_primitive(g, K) 

if K.is_negative(dup_LC(h, K)): 

c, h = -c, dup_neg(h, K) 

v = u - 1 

E = [ dmp_eval_tail(t, A, v, K) for t, _ in T ] 

D = dmp_zz_wang_non_divisors(E, c, ct, K) 

if D is not None: 

return c, h, E 

else: 

raise EvaluationFailed('no luck') 

def dmp_zz_wang_lead_coeffs(f, T, cs, E, H, A, u, K): 

"""Wang/EEZ: Compute correct leading coefficients. """ 

C, J, v = [], [0]*len(E), u - 1 

for h in H: 

c = dmp_one(v, K) 

d = dup_LC(h, K)*cs 

for i in reversed(range(len(E))): 

k, e, (t, _) = 0, E[i], T[i] 

while not (d % e): 

d, k = d//e, k + 1 

if k != 0: 

c, J[i] = dmp_mul(c, dmp_pow(t, k, v, K), v, K), 1 

C.append(c) 

if any(not j for j in J): 

raise ExtraneousFactors # pragma: no cover 

CC, HH = [], [] 

for c, h in zip(C, H): 

d = dmp_eval_tail(c, A, v, K) 

lc = dup_LC(h, K) 

if K.is_one(cs): 

cc = lc//d 

else: 

g = K.gcd(lc, d) 

d, cc = d//g, lc//g 

h, cs = dup_mul_ground(h, d, K), cs//d 

c = dmp_mul_ground(c, cc, v, K) 

CC.append(c) 

HH.append(h) 

if K.is_one(cs): 

return f, HH, CC 

CCC, HHH = [], [] 

for c, h in zip(CC, HH): 

CCC.append(dmp_mul_ground(c, cs, v, K)) 

HHH.append(dmp_mul_ground(h, cs, 0, K)) 

f = dmp_mul_ground(f, cs**(len(H) - 1), u, K) 

return f, HHH, CCC 

def dup_zz_diophantine(F, m, p, K): 

"""Wang/EEZ: Solve univariate Diophantine equations. """ 

if len(F) == 2: 

a, b = F 

f = gf_from_int_poly(a, p) 

g = gf_from_int_poly(b, p) 

s, t, G = gf_gcdex(g, f, p, K) 

s = gf_lshift(s, m, K) 

t = gf_lshift(t, m, K) 

q, s = gf_div(s, f, p, K) 

t = gf_add_mul(t, q, g, p, K) 

s = gf_to_int_poly(s, p) 

t = gf_to_int_poly(t, p) 

result = [s, t] 

else: 

G = [F[-1]] 

for f in reversed(F[1:-1]): 

G.insert(0, dup_mul(f, G[0], K)) 

S, T = [], [[1]] 

for f, g in zip(F, G): 

t, s = dmp_zz_diophantine([g, f], T[-1], [], 0, p, 1, K) 

T.append(t) 

S.append(s) 

result, S = [], S + [T[-1]] 

for s, f in zip(S, F):