Source code for diofant.polys.galoistools

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

import math
import random

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

[docs]def gf_crt(U, M, K=None):
"""
Chinese Remainder Theorem.

Given a set of integer residues u_0,...,u_n and a set of
co-prime integer moduli m_0,...,m_n, returns an integer
u, such that u = u_i mod m_i for i = 0,...,n.

As an example consider a set of residues U = [49, 76, 65]
and a set of moduli M = [99, 97, 95]. Then we have::

>>> from diofant.ntheory.modular import solve_congruence

>>> gf_crt([49, 76, 65], [99, 97, 95], ZZ)
639985

This is the correct result because::

>>> [639985 % m for m in [99, 97, 95]]
[49, 76, 65]

Note: this is a low-level routine with no error checking.

========

diofant.ntheory.modular.crt : a higher level crt routine
diofant.ntheory.modular.solve_congruence
"""
p = prod(M, start=K.one)
v = K.zero

for u, m in zip(U, M):
e = p // m
s, _, _ = K.gcdex(e, m)
v += e*(u*s % m)

return v % p

[docs]def gf_crt1(M, K):
"""
First part of the Chinese Remainder Theorem.

Examples
========

>>> gf_crt1([99, 97, 95], ZZ)
(912285, [9215, 9405, 9603], [62, 24, 12])
"""
E, S = [], []
p = prod(M, start=K.one)

for m in M:
E.append(p // m)
S.append(K.gcdex(E[-1], m)[0] % m)

return p, E, S

[docs]def gf_crt2(U, M, p, E, S, K):
"""
Second part of the Chinese Remainder Theorem.

Examples
========

>>> U = [49, 76, 65]
>>> M = [99, 97, 95]
>>> p = 912285
>>> E = [9215, 9405, 9603]
>>> S = [62, 24, 12]

>>> gf_crt2(U, M, p, E, S, ZZ)
639985
"""
v = K.zero

for u, m, e, s in zip(U, M, E, S):
v += e*(u*s % m)

return v % p

[docs]def gf_int(a, p):
"""
Coerce a mod p to an integer in the range [-p/2, p/2].

Examples
========

>>> gf_int(2, 7)
2
>>> gf_int(5, 7)
-2

"""
if a <= p // 2:
return a
else:
return a - p

[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_from_int_poly(f, p):
"""
Create a GF(p)[x] polynomial from Z[x].

Examples
========

>>> gf_from_int_poly([7, -2, 3], 5)
[2, 3, 3]
"""
return dmp_strip([a % p for a in f], 0)

[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 ]

"""
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)

"""
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_add_mul(f, g, h, p, K):
"""
Returns f + g*h where f, g, h in GF(p)[x].

Examples
========

>>> gf_add_mul([3, 2, 4], [2, 2, 2], [1, 4], 5, ZZ)
[2, 3, 2, 2]
"""
return gf_add(f, gf_mul(g, h, p, K), p, K)

[docs]def gf_sub_mul(f, g, h, p, K):
"""
Compute f - g*h where f, g, h in GF(p)[x].

Examples
========

>>> gf_sub_mul([3, 2, 4], [2, 2, 2], [1, 4], 5, ZZ)
[3, 3, 2, 1]
"""
return gf_sub(f, gf_mul(g, h, p, K), p, K)

[docs]def gf_expand(F, p, K):
"""
Expand results of :func:~diofant.polys.polytools.factor in GF(p)[x].

Examples
========

>>> gf_expand([([3, 2, 4], 1), ([2, 2], 2), ([3, 1], 3)], 5, ZZ)
[4, 3, 0, 3, 0, 1, 4, 1]
"""
if type(F) is tuple:
lc, F = F
else:
lc = K.one

g = [lc]

for f, k in F:
f = gf_pow(f, k, p, K)
g = gf_mul(g, f, p, K)

return g

[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.

Consider polynomials x**3 + x + 1 and x**2 + x in GF(2)::

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

As result we obtained quotient x + 1 and remainder 1, thus::

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

References
==========

* [Monagan93]_
* [Gathen99]_
"""
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_pow(f, n, p, K):
"""
Compute f**n in GF(p)[x] using repeated squaring.

Examples
========

>>> gf_pow([3, 2, 4], 3, 5, ZZ)
[2, 4, 4, 2, 2, 1, 4]
"""
if not n:
return [K.one]
elif n == 1:
return f
elif n == 2:
return gf_sqr(f, p, K)

h = [K.one]

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

n >>= 1

if not n:
break

f = gf_sqr(f, p, K)

return h

[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
========

>>> g = [1, 0, 2, 1]
>>> gf_frobenius_monomial_base(g, 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
==========

* [Gathen99]_
"""
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_lcm(f, g, p, K):
"""
Compute polynomial LCM in GF(p)[x].

Examples
========

>>> gf_lcm([3, 2, 4], [2, 2, 3], 5, ZZ)
[1, 2, 0, 4]
"""
if not f or not g:
return []

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

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

[docs]def gf_cofactors(f, g, p, K):
"""
Compute polynomial GCD and cofactors in GF(p)[x].

Examples
========

>>> gf_cofactors([3, 2, 4], [2, 2, 3], 5, ZZ)
([1, 3], [3, 3], [2, 1])
"""
if not f and not g:
return [], [], []

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

return (h, gf_quo(f, h, p, K),
gf_quo(g, h, p, K))

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

Given polynomials f and g in GF(p)[x], computes polynomials
s, t and h, such that h = gcd(f, g) and s*f + t*g = h.
The typical application of EEA is solving polynomial diophantine equations.

Consider polynomials f = (x + 7) (x + 1), g = (x + 7) (x**2 + 1)
in GF(11)[x]. Application of Extended Euclidean Algorithm gives::

>>> s, t, g = gf_gcdex([1, 8, 7], [1, 7, 1, 7], 11, ZZ)
>>> (s, t, g)
([5, 6], [6], [1, 7])

As result we obtained polynomials s = 5*x + 6 and t = 6, and
additionally gcd(f, g) = x + 7. This is correct because::

>>> S = gf_mul(s, [1, 8, 7], 11, ZZ)
>>> T = gf_mul(t, [1, 7, 1, 7], 11, ZZ)

[1, 7]

References
==========

* [Gathen99]_
"""
if not (f or g):
return [K.one], [], []

p0, r0 = gf_monic(f, p, K)
p1, r1 = gf_monic(g, p, K)

if not f:
return [], [K.invert(p1, p)], r1
if not g:
return [K.invert(p0, p)], [], r0

s0, s1 = [K.invert(p0, p)], []
t0, t1 = [], [K.invert(p1, p)]

while True:
Q, R = gf_div(r0, r1, p, K)

if not R:
break

(lc, r1), r0 = gf_monic(R, p, K), r1

inv = K.invert(lc, p)

s = gf_sub_mul(s0, s1, Q, p, K)
t = gf_sub_mul(t0, t1, Q, p, K)

s1, s0 = gf_mul_ground(s, inv, p, K), s1
t1, t0 = gf_mul_ground(t, inv, p, K), t1

return s1, t1, r1

[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_diff(f, p, K):
"""
Differentiate polynomial in GF(p)[x].

Examples
========

>>> gf_diff([3, 2, 4], 5, ZZ)
[1, 2]
"""
df = dmp_degree_in(f, 0, 0)
if df < 0:
return []

h, n = [K.zero]*df, df

for coeff in f[:-1]:
coeff *= K(n)
coeff %= p

if coeff:
h[df - n] = coeff

n -= 1

return dmp_strip(h, 0)

[docs]def gf_eval(f, a, p, K):
"""
Evaluate f(a) in GF(p) using Horner scheme.

Examples
========

>>> gf_eval([3, 2, 4], 2, 5, ZZ)
0
"""
result = K.zero

for c in f:
result *= a
result += c
result %= p

return result

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

Examples
========

>>> gf_compose([3, 2, 4], [2, 2, 2], 5, ZZ)
[2, 4, 0, 3, 0]
"""
if len(g) <= 1:
return dmp_strip([gf_eval(f, dmp_LC(g, K), p, K)], 0)

if not f:
return []

h = [f[0]]

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

return h

[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
==========

* [Gathen92]_
"""
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
==========

* [BenOr81]_
"""
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_sqf_list(f, p, K):
"""
Return the square-free decomposition of a GF(p)[x] polynomial.

Given a polynomial f in GF(p)[x], returns the leading coefficient
of f and a square-free decomposition f_1**e_1 f_2**e_2 ... f_k**e_k
such that all f_i are monic polynomials and (f_i, f_j) for i != j
are co-prime and e_1 ... e_k are given in increasing order. All trivial
terms (i.e. f_i = 1) aren't included in the output.

Consider polynomial f = x**11 + 1 over GF(11)[x]::

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

Note that f'(x) = 0::

>>> gf_diff(f, 11, ZZ)
[]

This phenomenon doesn't happen in characteristic zero. However we can
still compute square-free decomposition of f using gf_sqf()::

>>> gf_sqf_list(f, 11, ZZ)
(1, [([1, 1], 11)])

We obtained factorization f = (x + 1)**11. This is correct because::

>>> gf_pow([1, 1], 11, 11, ZZ) == f
True

References
==========

* [Geddes92]_
"""
n, sqf, factors, r = 1, False, [], int(p)

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

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

while True:
F = gf_diff(f, p, K)

if F != []:
g = gf_gcd(f, F, p, K)
h = gf_quo(f, g, p, K)

i = 1

while h != [K.one]:
G = gf_gcd(g, h, p, K)
H = gf_quo(h, G, p, K)

if dmp_degree_in(H, 0, 0) > 0:
factors.append((H, i*n))

g, h, i = gf_quo(g, G, p, K), G, i + 1

if g == [K.one]:
sqf = True
else:
f = g

if not sqf:
d = dmp_degree_in(f, 0, 0) // r

for i in range(d + 1):
f[i] = f[i*r]

f, n = f[:d + 1], n*r
else:
break

return lc, factors

[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.

Consider the polynomial x**15 - 1 in GF(11)[x]::

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

Distinct degree factorization gives::

>>> 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)]

which means x**15 - 1 = (x**5 - 1) (x**10 + x**5 + 1). To obtain
factorization into irreducibles, use equal degree factorization
procedure (EDF) with each of the factors.

References
==========

* [Gathen99]_
* [Geddes92]_
"""
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.

Consider the square-free polynomial f = x**3 + x**2 + x + 1 in
GF(5)[x]. Let's compute its irreducible factors of degree one::

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

References
==========

* [Gathen99]_
* [Geddes92]_
"""
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
==========

* [Kaltofen98]_
* [Shoup95]_
* [Gathen92]_
"""
n = dmp_degree_in(f, 0, 0)
k = int(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
==========

* [Shoup91]_
* [Gathen92]_
"""
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].

Examples
========

>>> gf_factor_sqf([3, 2, 4], 5, ZZ)
(3, [[1, 1], [1, 3]])
"""
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)

[docs]def gf_factor(f, p, K):
"""
Factor (non square-free) polynomials in GF(p)[x].

Given a possibly non square-free polynomial f in GF(p)[x],
returns its complete factorization into irreducibles::

f_1(x)**e_1 f_2(x)**e_2 ... f_d(x)**e_d

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 with
their multiplicities.

The algorithm proceeds by first computing square-free decomposition
of f and then iteratively factoring each of square-free factors.

Consider a non square-free polynomial f = (7*x + 1) (x + 2)**2 in
GF(11)[x]. We obtain its factorization into irreducibles as follows::

>>> gf_factor([5, 2, 7, 2], 11, ZZ)
(5, [([1, 2], 1), ([1, 8], 2)])

We arrived with factorization f = 5 (x + 2) (x + 8)**2. We didn't
recover the exact form of the input polynomial because we requested to
get monic factors of f and its leading coefficient separately.

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.

References
==========

* [Gathen99]_
"""
lc, f = gf_monic(f, p, K)

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

factors = []

for g, n in gf_sqf_list(f, p, K)[1]:
for h in gf_factor_sqf(g, p, K)[1]:
factors.append((h, n))

return lc, _sort_factors(factors)

[docs]def gf_value(f, a):
"""
Value of polynomial 'f' at 'a' in field R.

Examples
========

>>> gf_value([1, 7, 2, 4], 11)
2204

"""
result = 0
for c in f:
result *= a
result += c
return result

[docs]def linear_congruence(a, b, m):
"""
Returns the values of x satisfying a*x congruent b mod(m)

Here m is positive integer and a, b are natural numbers.
This function returns only those values of x which are distinct mod(m).

Examples
========

>>> linear_congruence(3, 12, 15)
[4, 9, 14]

There are 3 solutions distinct mod(15) since gcd(a, m) = gcd(3, 15) = 3.

References
==========

* https//en.wikipedia.org/wiki/Linear_congruence_theorem
"""
from .polytools import gcdex
if a % m == 0:
if b % m == 0:
return list(range(m))
else:
return []
r, _, g = gcdex(a, m)
if b % g != 0:
return []
return [(r * b // g + t * m // g) % m for t in range(g)]

def _raise_mod_power(x, s, p, f):
"""
Used in gf_csolve to generate solutions of f(x) cong 0 mod(p**(s + 1))
from the solutions of f(x) cong 0 mod(p**s).

Examples
========

These is the solutions of f(x) = x**2 + x + 7 cong 0 mod(3)

>>> f = [1, 1, 7]
>>> csolve_prime(f, 3)
[1]
>>> [i for i in range(3) if not (i**2 + i + 7) % 3]
[1]

The solutions of f(x) cong 0 mod(9) are constructed from the
values returned from _raise_mod_power:

>>> x, s, p = 1, 1, 3
>>> V = _raise_mod_power(x, s, p, f)
>>> [x + v * p**s for v in V]
[1, 4, 7]

And these are confirmed with the following:

>>> [i for i in range(3**2) if not (i**2 + i + 7) % 3**2]
[1, 4, 7]

"""
from ..domains import ZZ
f_f = gf_diff(f, p, ZZ)
alpha = gf_value(f_f, x)
beta = - gf_value(f, x) // p**s
return linear_congruence(alpha, beta, p)

[docs]def csolve_prime(f, p, e=1):
"""
Solutions of f(x) congruent 0 mod(p**e).

Examples
========

>>> csolve_prime([1, 1, 7], 3, 1)
[1]
>>> csolve_prime([1, 1, 7], 3, 2)
[1, 4, 7]

Solutions [7, 4, 1] (mod 3**2) are generated by _raise_mod_power()
from solution [1] (mod 3).
"""
from ..domains import ZZ
X1 = [i for i in range(p) if gf_eval(f, i, p, ZZ) == 0]
if e == 1:
return X1
X = []
S = list(zip(X1, [1]*len(X1)))
while S:
x, s = S.pop()
if s == e:
X.append(x)
else:
s1 = s + 1
ps = p**s
S.extend([(x + v*ps, s1) for v in _raise_mod_power(x, s, p, f)])
return sorted(X)

[docs]def gf_csolve(f, n):
"""
To solve f(x) congruent 0 mod(n).

n is divided into canonical factors and f(x) cong 0 mod(p**e) will be
solved for each factor. Applying the Chinese Remainder Theorem to the

Examples
========

Solve [1, 1, 7] congruent 0 mod(189):

>>> gf_csolve([1, 1, 7], 189)
[13, 49, 76, 112, 139, 175]

References
==========

* [Niven91]_
"""
from ..domains import ZZ
P = factorint(n)
X = [csolve_prime(f, p, e) for p, e in P.items()]
pools = list(map(tuple, X))
perms = [[]]
for pool in pools:
perms = [x + [y] for x in perms for y in pool]
dist_factors = [pow(p, e) for p, e in P.items()]
return sorted(gf_crt(per, dist_factors, ZZ) for per in perms)