# Source code for diofant.polys.galoistools

"""Dense univariate polynomials with coefficients in Galois fields. """

import math
import random

from ..ntheory import factorint
from .densearith import dup_lshift
from .densebasic import (dmp_convert, dmp_degree_in, dmp_from_dict, dmp_normal,
dmp_strip)
from .polyconfig import query
from .polyerrors import ExactQuotientFailed
from .polyutils import _sort_factors

[docs]def gf_from_dict(f, p, K):
"""
Create a GF(p)[x] polynomial from a dict.

Examples
========

>>> gf_from_dict({10: 4, 4: 33, 0: -1}, 5, ZZ)
[4, 0, 0, 0, 0, 0, 3, 0, 0, 0, 4]

"""
f = dmp_from_dict(f, 0, K.finite_field(p))
f = dmp_normal(f, 0, K.finite_field(p))

return dmp_convert(f, 0, K.finite_field(p), K)

[docs]def gf_neg(f, p, K):
"""
Negate a polynomial in GF(p)[x].

Examples
========

>>> gf_neg([3, 2, 1, 0], 5, ZZ)
[2, 3, 4, 0]

"""
return [-coeff % p for coeff in f]

[docs]def gf_add_ground(f, a, p, K):
"""
Compute f + a where f in GF(p)[x] and a in GF(p).

Examples
========

>>> gf_add_ground([3, 2, 4], 2, 5, ZZ)
[3, 2, 1]

"""
if not f:
a = a % p
else:
a = (f[-1] + a) % p

if len(f) > 1:
return f[:-1] + [a]

if not a:
return []
else:
return [a]

[docs]def gf_sub_ground(f, a, p, K):
"""
Compute f - a where f in GF(p)[x] and a in GF(p).

Examples
========

>>> gf_sub_ground([3, 2, 4], 2, 5, ZZ)
[3, 2, 2]

"""
if not f:
a = -a % p
else:
a = (f[-1] - a) % p

if len(f) > 1:
return f[:-1] + [a]

if not a:
return []
else:
return [a]

[docs]def gf_mul_ground(f, a, p, K):
"""
Compute f * a where f in GF(p)[x] and a in GF(p).

Examples
========

>>> gf_mul_ground([3, 2, 4], 2, 5, ZZ)
[1, 4, 3]

"""
if not a:
return []
else:
return [(a*b) % p for b in f]

[docs]def gf_quo_ground(f, a, p, K):
"""
Compute f/a where f in GF(p)[x] and a in GF(p).

Examples
========

>>> gf_quo_ground([3, 2, 4], 2, 5, ZZ)
[4, 1, 2]

"""
return gf_mul_ground(f, K.invert(a, p), p, K)

[docs]def gf_add(f, g, p, K):
"""
Add polynomials in GF(p)[x].

Examples
========

>>> gf_add([3, 2, 4], [2, 2, 2], 5, ZZ)
[4, 1]

"""
if not f:
return g
if not g:
return f

df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

if df == dg:
return dmp_strip([(a + b) % p for a, b in zip(f, g)], 0)
else:
k = abs(df - dg)

if df > dg:
h, f = f[:k], f[k:]
else:
h, g = g[:k], g[k:]

return h + [(a + b) % p for a, b in zip(f, g)]

[docs]def gf_sub(f, g, p, K):
"""
Subtract polynomials in GF(p)[x].

Examples
========

>>> gf_sub([3, 2, 4], [2, 2, 2], 5, ZZ)
[1, 0, 2]

"""
if not g:
return f
if not f:
return gf_neg(g, p, K)

df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

if df == dg:
return dmp_strip([(a - b) % p for a, b in zip(f, g)], 0)
else:
k = abs(df - dg)

if df > dg:
h, f = f[:k], f[k:]
else:
h, g = gf_neg(g[:k], p, K), g[k:]

return h + [(a - b) % p for a, b in zip(f, g)]

[docs]def gf_mul(f, g, p, K):
"""
Multiply polynomials in GF(p)[x].

Examples
========

>>> gf_mul([3, 2, 4], [2, 2, 2], 5, ZZ)
[1, 0, 3, 2, 3]

"""
df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

dh = df + dg
if dh < 0:
return []

h = [0]*(dh + 1)

for i in range(dh + 1):
coeff = K.zero

for j in range(max(0, i - dg), min(i, df) + 1):
coeff += f[j]*g[i - j]

h[i] = coeff % p

return dmp_strip(h, 0)

[docs]def gf_sqr(f, p, K):
"""
Square polynomials in GF(p)[x].

Examples
========

>>> gf_sqr([3, 2, 4], 5, ZZ)
[4, 2, 3, 1, 1]

"""
df = dmp_degree_in(f, 0, 0)
if df < 0:
return []

dh = 2*df
h = [0]*(dh + 1)

for i in range(dh + 1):
coeff = K.zero

jmin = max(0, i - df)
jmax = min(i, df)

n = jmax - jmin + 1

jmax = jmin + n // 2 - 1

for j in range(jmin, jmax + 1):
coeff += f[j]*f[i - j]

coeff += coeff

if n & 1:
elem = f[jmax + 1]
coeff += elem**2

h[i] = coeff % p

return dmp_strip(h, 0)

[docs]def gf_div(f, g, p, K):
"""
Division with remainder in GF(p)[x].

Given univariate polynomials f and g with coefficients in a
finite field with p elements, returns polynomials q and r
(quotient and remainder) such that f = q*g + r.

Examples
========

>>> gf_div([1, 0, 1, 1], [1, 1, 0], 2, ZZ)
([1, 1], [1])

References
==========

* :cite:Monagan1993inplace
* :cite:Gathen1999modern

"""
df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

if not g:
raise ZeroDivisionError("polynomial division")
elif df < dg:
return [], f

inv = K.invert(g[0], p)

h, dq, dr = list(f), df - dg, dg - 1

for i in range(df + 1):
coeff = h[i]

for j in range(max(0, dg - i), min(df - i, dr) + 1):
coeff -= h[i + j - dg] * g[dg - j]

if i <= dq:
coeff *= inv

h[i] = coeff % p

return h[:dq + 1], dmp_strip(h[dq + 1:], 0)

[docs]def gf_rem(f, g, p, K):
"""
Compute polynomial remainder in GF(p)[x].

Examples
========

>>> gf_rem([1, 0, 1, 1], [1, 1, 0], 2, ZZ)
[1]

"""
return gf_div(f, g, p, K)[1]

[docs]def gf_quo(f, g, p, K):
"""
Compute exact quotient in GF(p)[x].

Examples
========

>>> gf_quo([1, 0, 1, 1], [1, 1, 0], 2, ZZ)
[1, 1]
>>> gf_quo([1, 0, 3, 2, 3], [2, 2, 2], 5, ZZ)
[3, 2, 4]

"""
df = dmp_degree_in(f, 0, 0)
dg = dmp_degree_in(g, 0, 0)

if not g:
raise ZeroDivisionError("polynomial division")
elif df < dg:
return []

inv = K.invert(g[0], p)

h, dq, dr = f[:], df - dg, dg - 1

for i in range(dq + 1):
coeff = h[i]

for j in range(max(0, dg - i), min(df - i, dr) + 1):
coeff -= h[i + j - dg] * g[dg - j]

h[i] = (coeff * inv) % p

return h[:dq + 1]

[docs]def gf_exquo(f, g, p, K):
"""
Compute polynomial quotient in GF(p)[x].

Examples
========

>>> gf_exquo([1, 0, 3, 2, 3], [2, 2, 2], 5, ZZ)
[3, 2, 4]

>>> gf_exquo([1, 0, 1, 1], [1, 1, 0], 2, ZZ)
Traceback (most recent call last):
...
ExactQuotientFailed: [1, 1, 0] does not divide [1, 0, 1, 1]

"""
q, r = gf_div(f, g, p, K)

if not r:
return q
else:
raise ExactQuotientFailed(f, g)

[docs]def gf_frobenius_monomial_base(g, p, K):
"""
return the list of x**(i*p) mod g in Z_p for i = 0, .., n - 1
where n = dmp_degree_in(g, 0, 0)

Examples
========

>>> gf_frobenius_monomial_base([1, 0, 2, 1], 5, ZZ)
[[1], [4, 4, 2], [1, 2]]

"""
n = dmp_degree_in(g, 0, 0)
if n == 0:
return []
b = [0]*n
b[0] = [1]
if p < n:
for i in range(1, n):
mon = dup_lshift(b[i - 1], p, K)
b[i] = gf_rem(mon, g, p, K)
elif n > 1:
b[1] = gf_pow_mod([K.one, K.zero], p, g, p, K)
for i in range(2, n):
b[i] = gf_mul(b[i - 1], b[1], p, K)
b[i] = gf_rem(b[i], g, p, K)

return b

[docs]def gf_frobenius_map(f, g, b, p, K):
"""
compute gf_pow_mod(f, p, g, p, K) using the Frobenius map

Parameters
==========

f, g : polynomials in GF(p)[x]
b : frobenius monomial base
p : prime number
K : domain

Examples
========

>>> f = [2, 1, 0, 1]
>>> g = [1, 0, 2, 1]
>>> p = 5
>>> b = gf_frobenius_monomial_base(g, p, ZZ)
>>> r = gf_frobenius_map(f, g, b, p, ZZ)
>>> gf_frobenius_map(f, g, b, p, ZZ)
[4, 0, 3]

"""
m = dmp_degree_in(g, 0, 0)
if dmp_degree_in(f, 0, 0) >= m:
f = gf_rem(f, g, p, K)
if not f:
return []
n = dmp_degree_in(f, 0, 0)
sf = [f[-1]]
for i in range(1, n + 1):
v = gf_mul_ground(b[i], f[n - i], p, K)
sf = gf_add(sf, v, p, K)
return sf

def _gf_pow_pnm1d2(f, n, g, b, p, K):
"""
utility function for gf_edf_zassenhaus
Compute f**((p**n - 1) // 2) in GF(p)[x]/(g)
f**((p**n - 1) // 2) = (f*f**p*...*f**(p**n - 1))**((p - 1) // 2)

"""
f = gf_rem(f, g, p, K)
h = f
r = f
for i in range(1, n):
h = gf_frobenius_map(h, g, b, p, K)
r = gf_mul(r, h, p, K)
r = gf_rem(r, g, p, K)

res = gf_pow_mod(r, (p - 1)//2, g, p, K)
return res

[docs]def gf_pow_mod(f, n, g, p, K):
"""
Compute f**n in GF(p)[x]/(g) using repeated squaring.

Given polynomials f and g in GF(p)[x] and a non-negative
integer n, efficiently computes f**n (mod g) i.e. the remainder
of f**n from division by g, using the repeated squaring algorithm.

Examples
========

>>> gf_pow_mod([3, 2, 4], 3, [1, 1], 5, ZZ)
[]

References
==========

* :cite:Gathen1999modern

"""
if not n:
return [K.one]
elif n == 1:
return gf_rem(f, g, p, K)
elif n == 2:
return gf_rem(gf_sqr(f, p, K), g, p, K)

h = [K.one]

while True:
if n & 1:
h = gf_mul(h, f, p, K)
h = gf_rem(h, g, p, K)
n -= 1

n >>= 1

if not n:
break

f = gf_sqr(f, p, K)
f = gf_rem(f, g, p, K)

return h

[docs]def gf_gcd(f, g, p, K):
"""
Euclidean Algorithm in GF(p)[x].

Examples
========

>>> gf_gcd([3, 2, 4], [2, 2, 3], 5, ZZ)
[1, 3]

"""
while g:
f, g = g, gf_rem(f, g, p, K)

return gf_monic(f, p, K)[1]

[docs]def gf_monic(f, p, K):
"""
Compute LC and a monic polynomial in GF(p)[x].

Examples
========

>>> gf_monic([3, 2, 4], 5, ZZ)
(3, [1, 4, 3])

"""
if not f:
return K.zero, []
else:
lc = f[0]

if lc == K.one:
return lc, list(f)
else:
return lc, gf_quo_ground(f, lc, p, K)

[docs]def gf_compose_mod(g, h, f, p, K):
"""
Compute polynomial composition g(h) in GF(p)[x]/(f).

Examples
========

>>> gf_compose_mod([3, 2, 4], [2, 2, 2], [4, 3], 5, ZZ)
[4]

"""
if not g:
return []

comp = [g[0]]

for a in g[1:]:
comp = gf_mul(comp, h, p, K)
comp = gf_add_ground(comp, a, p, K)
comp = gf_rem(comp, f, p, K)

return comp

[docs]def gf_trace_map(a, b, c, n, f, p, K):
"""
Compute polynomial trace map in GF(p)[x]/(f).

Given a polynomial f in GF(p)[x], polynomials a, b,
c in the quotient ring GF(p)[x]/(f) such that b = c**t
(mod f) for some positive power t of p, and a positive
integer n, returns a mapping::

a -> a**t**n, a + a**t + a**t**2 + ... + a**t**n (mod f)

In factorization context, b = x**p mod f and c = x mod f.
This way we can efficiently compute trace polynomials in equal
degree factorization routine, much faster than with other methods,
like iterated Frobenius algorithm, for large degrees.

Examples
========

>>> gf_trace_map([1, 2], [4, 4], [1, 1], 4, [3, 2, 4], 5, ZZ)
([1, 3], [1, 3])

References
==========

* :cite:Gathen1992frobenious

"""
u = gf_compose_mod(a, b, f, p, K)
v = b

if n & 1:
U = gf_add(a, u, p, K)
V = b
else:
U = a
V = c

n >>= 1

while n:
u = gf_add(u, gf_compose_mod(u, v, f, p, K), p, K)
v = gf_compose_mod(v, v, f, p, K)

if n & 1:
U = gf_add(U, gf_compose_mod(u, V, f, p, K), p, K)
V = gf_compose_mod(v, V, f, p, K)

n >>= 1

return gf_compose_mod(a, V, f, p, K), U

def _gf_trace_map(f, n, g, b, p, K):
"""
utility for gf_edf_shoup

"""
f = gf_rem(f, g, p, K)
h = f
r = f
for i in range(1, n):
h = gf_frobenius_map(h, g, b, p, K)
r = gf_add(r, h, p, K)
r = gf_rem(r, g, p, K)
return r

[docs]def gf_random(n, p, K):
"""
Generate a random polynomial in GF(p)[x] of degree n.

Examples
========

>>> gf_random(10, 5, ZZ) #doctest: +SKIP
[1, 2, 3, 2, 1, 1, 1, 2, 0, 4, 2]

"""
return [K.one] + [K(int(random.uniform(0, p))) for i in range(n)]

[docs]def gf_irreducible(n, p, K):
"""
Generate random irreducible polynomial of degree n in GF(p)[x].

Examples
========

>>> gf_irreducible(10, 5, ZZ) #doctest: +SKIP
[1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4]

"""
while True:
f = gf_random(n, p, K)
if gf_irreducible_p(f, p, K):
return f

[docs]def gf_irred_p_ben_or(f, p, K):
"""
Ben-Or's polynomial irreducibility test over finite fields.

Examples
========

>>> gf_irred_p_ben_or([1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4], 5, ZZ)
True
>>> gf_irred_p_ben_or([3, 2, 4], 5, ZZ)
False

References
==========

* :cite:Ben-Or1981ff

"""
n = dmp_degree_in(f, 0, 0)

if n <= 1:
return True

_, f = gf_monic(f, p, K)
if n < 5:
H = h = gf_pow_mod([K.one, K.zero], p, f, p, K)

for i in range(n//2):
g = gf_sub(h, [K.one, K.zero], p, K)

if gf_gcd(f, g, p, K) == [K.one]:
h = gf_compose_mod(h, H, f, p, K)
else:
return False
else:
b = gf_frobenius_monomial_base(f, p, K)
H = h = gf_frobenius_map([K.one, K.zero], f, b, p, K)
for i in range(n//2):
g = gf_sub(h, [K.one, K.zero], p, K)
if gf_gcd(f, g, p, K) == [K.one]:
h = gf_frobenius_map(h, f, b, p, K)
else:
return False

return True

[docs]def gf_irred_p_rabin(f, p, K):
"""
Rabin's polynomial irreducibility test over finite fields.

Examples
========

>>> gf_irred_p_rabin([1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4], 5, ZZ)
True
>>> gf_irred_p_rabin([3, 2, 4], 5, ZZ)
False

"""
n = dmp_degree_in(f, 0, 0)

if n <= 1:
return True

_, f = gf_monic(f, p, K)

x = [K.one, K.zero]

indices = {n//d for d in factorint(n)}

b = gf_frobenius_monomial_base(f, p, K)
h = b[1]

for i in range(1, n):
if i in indices:
g = gf_sub(h, x, p, K)

if gf_gcd(f, g, p, K) != [K.one]:
return False

h = gf_frobenius_map(h, f, b, p, K)

return h == x

_irred_methods = {
'ben-or': gf_irred_p_ben_or,
'rabin': gf_irred_p_rabin,
}

[docs]def gf_irreducible_p(f, p, K):
"""
Test irreducibility of a polynomial f in GF(p)[x].

Examples
========

>>> gf_irreducible_p([1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4], 5, ZZ)
True
>>> gf_irreducible_p([3, 2, 4], 5, ZZ)
False

"""
method = query('GF_IRRED_METHOD')

return _irred_methods[method](f, p, K)

[docs]def gf_Qmatrix(f, p, K):
"""
Calculate Berlekamp's Q matrix.

Examples
========

>>> gf_Qmatrix([3, 2, 4], 5, ZZ)
[[1, 0],
[3, 4]]

>>> gf_Qmatrix([1, 0, 0, 0, 1], 5, ZZ)
[[1, 0, 0, 0],
[0, 4, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 4]]

"""
n, r = dmp_degree_in(f, 0, 0), int(p)

q = [K.one] + [K.zero]*(n - 1)
Q = [list(q)] + [[]]*(n - 1)

for i in range(1, (n - 1)*r + 1):
qq, c = [(-q[-1]*f[-1]) % p], q[-1]

for j in range(1, n):
qq.append((q[j - 1] - c*f[-j - 1]) % p)

if not (i % r):
Q[i//r] = list(qq)

q = qq

return Q

[docs]def gf_Qbasis(Q, p, K):
"""
Compute a basis of the kernel of Q.

Examples
========

>>> gf_Qbasis(gf_Qmatrix([1, 0, 0, 0, 1], 5, ZZ), 5, ZZ)
[[1, 0, 0, 0], [0, 0, 1, 0]]

>>> gf_Qbasis(gf_Qmatrix([3, 2, 4], 5, ZZ), 5, ZZ)
[[1, 0]]

"""
Q, n = [list(q) for q in Q], len(Q)

for k in range(n):
Q[k][k] = (Q[k][k] - K.one) % p

for k in range(n):
for i in range(k, n):
if Q[k][i]:
break
else:
continue

inv = K.invert(Q[k][i], p)

for j in range(n):
Q[j][i] = (Q[j][i]*inv) % p

for j in range(n):
t = Q[j][k]
Q[j][k] = Q[j][i]
Q[j][i] = t

for i in range(n):
if i != k:
q = Q[k][i]

for j in range(n):
Q[j][i] = (Q[j][i] - Q[j][k]*q) % p

for i in range(n):
for j in range(n):
if i == j:
Q[i][j] = (K.one - Q[i][j]) % p
else:
Q[i][j] = (-Q[i][j]) % p

basis = []

for q in Q:
if any(q):
basis.append(q)

return basis

[docs]def gf_berlekamp(f, p, K):
"""
Factor a square-free f in GF(p)[x] for small p.

Examples
========

>>> gf_berlekamp([1, 0, 0, 0, 1], 5, ZZ)
[[1, 0, 2], [1, 0, 3]]

"""
Q = gf_Qmatrix(f, p, K)
V = gf_Qbasis(Q, p, K)

for i, v in enumerate(V):
V[i] = dmp_strip(list(reversed(v)), 0)

factors = [f]

for k in range(1, len(V)):
for f in list(factors):
s = K.zero

while s < p:
g = gf_sub_ground(V[k], s, p, K)
h = gf_gcd(f, g, p, K)

if h != [K.one] and h != f:
factors.remove(f)

f = gf_quo(f, h, p, K)
factors.extend([f, h])

if len(factors) == len(V):
return _sort_factors(factors, multiple=False)

s += K.one

return _sort_factors(factors, multiple=False)

[docs]def gf_ddf_zassenhaus(f, p, K):
"""
Cantor-Zassenhaus: Deterministic Distinct Degree Factorization

Given a monic square-free polynomial f in GF(p)[x], computes
partial distinct degree factorization f_1 ... f_d of f where
deg(f_i) != deg(f_j) for i != j. The result is returned as a
list of pairs (f_i, e_i) where deg(f_i) > 0 and e_i > 0
is an argument to the equal degree factorization routine.

Examples
========

>>> f = gf_from_dict({15: ZZ(1), 0: ZZ(-1)}, 11, ZZ)

>>> gf_ddf_zassenhaus(f, 11, ZZ)
[([1, 0, 0, 0, 0, 10], 1), ([1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1], 2)]

To obtain factorization into irreducibles, use equal degree factorization
procedure (EDF) with each of the factors.

References
==========

* :cite:Gathen1999modern
* :cite:Geddes1992algorithms

"""
i, g, factors = 1, [K.one, K.zero], []

b = gf_frobenius_monomial_base(f, p, K)
while 2*i <= dmp_degree_in(f, 0, 0):
g = gf_frobenius_map(g, f, b, p, K)
h = gf_gcd(f, gf_sub(g, [K.one, K.zero], p, K), p, K)

if h != [K.one]:
factors.append((h, i))

f = gf_quo(f, h, p, K)
g = gf_rem(g, f, p, K)
b = gf_frobenius_monomial_base(f, p, K)

i += 1

if f != [K.one]:
return factors + [(f, dmp_degree_in(f, 0, 0))]
else:
return factors

[docs]def gf_edf_zassenhaus(f, n, p, K):
"""
Cantor-Zassenhaus: Probabilistic Equal Degree Factorization

Given a monic square-free polynomial f in GF(p)[x] and
an integer n, such that n divides deg(f), returns all
irreducible factors f_1,...,f_d of f, each of degree n.
EDF procedure gives complete factorization over Galois fields.

Examples
========

>>> gf_edf_zassenhaus([1, 1, 1, 1], 1, 5, ZZ)
[[1, 1], [1, 2], [1, 3]]

References
==========

* :cite:Gathen1999modern
* :cite:Geddes1992algorithms

"""
factors = [f]

if dmp_degree_in(f, 0, 0) <= n:
return factors

N = dmp_degree_in(f, 0, 0) // n
if p != 2:
b = gf_frobenius_monomial_base(f, p, K)

while len(factors) < N:
r = gf_random(2*n - 1, p, K)

if p == 2:
h = r

for i in range(2**(n*N - 1)):
r = gf_pow_mod(r, 2, f, p, K)
h = gf_add(h, r, p, K)

g = gf_gcd(f, h, p, K)
else:
h = _gf_pow_pnm1d2(r, n, f, b, p, K)
g = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)

if g != [K.one] and g != f:
factors = gf_edf_zassenhaus(g, n, p, K) \
+ gf_edf_zassenhaus(gf_quo(f, g, p, K), n, p, K)

return _sort_factors(factors, multiple=False)

[docs]def gf_ddf_shoup(f, p, K):
"""
Kaltofen-Shoup: Deterministic Distinct Degree Factorization

Given a monic square-free polynomial f in GF(p)[x], computes
partial distinct degree factorization f_1,...,f_d of f where
deg(f_i) != deg(f_j) for i != j. The result is returned as a
list of pairs (f_i, e_i) where deg(f_i) > 0 and e_i > 0
is an argument to the equal degree factorization routine.

This algorithm is an improved version of Zassenhaus algorithm for
large deg(f) and modulus p (especially for deg(f) ~ lg(p)).

Examples
========

>>> f = gf_from_dict({6: ZZ(1), 5: ZZ(-1), 4: ZZ(1), 3: ZZ(1), 1: ZZ(-1)}, 3, ZZ)

>>> gf_ddf_shoup(f, 3, ZZ)
[([1, 1, 0], 1), ([1, 1, 0, 1, 2], 2)]

References
==========

* :cite:Kaltofen1998subquadratic
* :cite:Shoup1995factor
* :cite:Gathen1992frobenious

"""
n = dmp_degree_in(f, 0, 0)
k = math.ceil(math.sqrt(n//2))
b = gf_frobenius_monomial_base(f, p, K)
h = gf_frobenius_map([K.one, K.zero], f, b, p, K)
# U[i] = x**(p**i)
U = [[K.one, K.zero], h] + [K.zero]*(k - 1)

for i in range(2, k + 1):
U[i] = gf_frobenius_map(U[i-1], f, b, p, K)

h, U = U[k], U[:k]
# V[i] = x**(p**(k*(i+1)))
V = [h] + [K.zero]*(k - 1)

for i in range(1, k):
V[i] = gf_compose_mod(V[i - 1], h, f, p, K)

factors = []

for i, v in enumerate(V):
h, j = [K.one], k - 1

for u in U:
g = gf_sub(v, u, p, K)
h = gf_mul(h, g, p, K)
h = gf_rem(h, f, p, K)

g = gf_gcd(f, h, p, K)
f = gf_quo(f, g, p, K)

for u in reversed(U):
h = gf_sub(v, u, p, K)
F = gf_gcd(g, h, p, K)

if F != [K.one]:
factors.append((F, k*(i + 1) - j))

g, j = gf_quo(g, F, p, K), j - 1

if f != [K.one]:
factors.append((f, dmp_degree_in(f, 0, 0)))

return factors

[docs]def gf_edf_shoup(f, n, p, K):
"""
Gathen-Shoup: Probabilistic Equal Degree Factorization

Given a monic square-free polynomial f in GF(p)[x] and integer
n such that n divides deg(f), returns all irreducible factors
f_1,...,f_d of f, each of degree n. This is a complete
factorization over Galois fields.

This algorithm is an improved version of Zassenhaus algorithm for
large deg(f) and modulus p (especially for deg(f) ~ lg(p)).

Examples
========

>>> gf_edf_shoup([1, 2837, 2277], 1, 2917, ZZ)
[[1, 852], [1, 1985]]

References
==========

* :cite:Shoup1991ffactor
* :cite:Gathen1992frobenious

"""
N, q = dmp_degree_in(f, 0, 0), int(p)

if not N:
return []
if N <= n:
return [f]

factors, x = [f], [K.one, K.zero]

r = gf_random(N - 1, p, K)

if p == 2:
h = gf_pow_mod(x, q, f, p, K)
H = gf_trace_map(r, h, x, n - 1, f, p, K)[1]
h1 = gf_gcd(f, H, p, K)
h2 = gf_quo(f, h1, p, K)

factors = gf_edf_shoup(h1, n, p, K) \
+ gf_edf_shoup(h2, n, p, K)
else:
b = gf_frobenius_monomial_base(f, p, K)
H = _gf_trace_map(r, n, f, b, p, K)
h = gf_pow_mod(H, (q - 1)//2, f, p, K)

h1 = gf_gcd(f, h, p, K)
h2 = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)
h3 = gf_quo(f, gf_mul(h1, h2, p, K), p, K)

factors = gf_edf_shoup(h1, n, p, K) \
+ gf_edf_shoup(h2, n, p, K) \
+ gf_edf_shoup(h3, n, p, K)

return _sort_factors(factors, multiple=False)

[docs]def gf_zassenhaus(f, p, K):
"""
Factor a square-free f in GF(p)[x] for medium p.

Examples
========

>>> gf_zassenhaus([1, 4, 3], 5, ZZ)
[[1, 1], [1, 3]]

"""
factors = []

for factor, n in gf_ddf_zassenhaus(f, p, K):
factors += gf_edf_zassenhaus(factor, n, p, K)

return _sort_factors(factors, multiple=False)

[docs]def gf_shoup(f, p, K):
"""
Factor a square-free f in GF(p)[x] for large p.

Examples
========

>>> gf_shoup([1, 4, 3], 5, ZZ)
[[1, 1], [1, 3]]

"""
factors = []

for factor, n in gf_ddf_shoup(f, p, K):
factors += gf_edf_shoup(factor, n, p, K)

return _sort_factors(factors, multiple=False)

_factor_methods = {
'berlekamp': gf_berlekamp,  # p : small
'zassenhaus': gf_zassenhaus,  # p : medium
'shoup': gf_shoup,      # p : large
}

[docs]def gf_factor_sqf(f, p, K):
"""
Factor a square-free polynomial f in GF(p)[x].

Returns its complete factorization into irreducibles::

f_1(x) f_2(x) ... f_d(x)

where each f_i is a monic polynomial and gcd(f_i, f_j) == 1,
for i != j.  The result is given as a tuple consisting of the
leading coefficient of f and a list of factors of f.

Square-free factors of f can be factored into irreducibles over
GF(p) using three very different methods:

Berlekamp
efficient for very small values of p (usually p < 25)
Cantor-Zassenhaus
efficient on average input and with "typical" p
Shoup-Kaltofen-Gathen
efficient with very large inputs and modulus

If you want to use a specific factorization method, instead of the default
one, set GF_FACTOR_METHOD with one of berlekamp, zassenhaus or
shoup values.

Examples
========

>>> gf_factor_sqf([3, 2, 4], 5, ZZ)
(3, [[1, 1], [1, 3]])

References
==========

* :cite:Gathen1999modern

"""
lc, f = gf_monic(f, p, K)

if dmp_degree_in(f, 0, 0) < 1:
return lc, []

method = query('GF_FACTOR_METHOD')

return lc, _factor_methods[method](f, p, K)