# Source code for sympy.polys.galoistools

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

from __future__ import print_function, division

from random import uniform
from math import ceil as _ceil, sqrt as _sqrt

from sympy.core.compatibility import SYMPY_INTS, range
from sympy.core.mul import prod
from sympy.polys.polyutils import _sort_factors
from sympy.polys.polyconfig import query
from sympy.polys.polyerrors import ExactQuotientFailed

from sympy.ntheory import factorint

[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 sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_crt
>>> from sympy.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.

========

sympy.ntheory.modular.crt : a higher level crt routine
sympy.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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_crt1

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_crt2

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

>>> from sympy.polys.galoistools import gf_int

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

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

[docs]def gf_degree(f):
"""
Return the leading degree of f.

Examples
========

>>> from sympy.polys.galoistools import gf_degree

>>> gf_degree([1, 1, 2, 0])
3
>>> gf_degree([])
-1

"""
return len(f) - 1

[docs]def gf_LC(f, K):
"""
Return the leading coefficient of f.

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_LC

>>> gf_LC([3, 0, 1], ZZ)
3

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

[docs]def gf_TC(f, K):
"""
Return the trailing coefficient of f.

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_TC

>>> gf_TC([3, 0, 1], ZZ)
1

"""
if not f:
return K.zero
else:
return f[-1]

[docs]def gf_strip(f):
"""
Remove leading zeros from f.

Examples
========

>>> from sympy.polys.galoistools import gf_strip

>>> gf_strip([0, 0, 0, 3, 0, 1])
[3, 0, 1]

"""
if not f or f[0]:
return f

k = 0

for coeff in f:
if coeff:
break
else:
k += 1

return f[k:]

[docs]def gf_trunc(f, p):
"""
Reduce all coefficients modulo p.

Examples
========

>>> from sympy.polys.galoistools import gf_trunc

>>> gf_trunc([7, -2, 3], 5)
[2, 3, 3]

"""
return gf_strip([ a % p for a in f ])

[docs]def gf_normal(f, p, K):
"""
Normalize all coefficients in K.

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_normal

>>> gf_normal([5, 10, 21, -3], 5, ZZ)
[1, 2]

"""
return gf_trunc(list(map(K, f)), p)

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_from_dict

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

"""
n, h = max(f.keys()), []

if isinstance(n, SYMPY_INTS):
for k in range(n, -1, -1):
h.append(f.get(k, K.zero) % p)
else:
(n,) = n

for k in range(n, -1, -1):
h.append(f.get((k,), K.zero) % p)

return gf_trunc(h, p)

[docs]def gf_to_dict(f, p, symmetric=True):
"""
Convert a GF(p)[x] polynomial to a dict.

Examples
========

>>> from sympy.polys.galoistools import gf_to_dict

>>> gf_to_dict([4, 0, 0, 0, 0, 0, 3, 0, 0, 0, 4], 5)
{0: -1, 4: -2, 10: -1}
>>> gf_to_dict([4, 0, 0, 0, 0, 0, 3, 0, 0, 0, 4], 5, symmetric=False)
{0: 4, 4: 3, 10: 4}

"""
n, result = gf_degree(f), {}

for k in range(0, n + 1):
if symmetric:
a = gf_int(f[n - k], p)
else:
a = f[n - k]

if a:
result[k] = a

return result

[docs]def gf_from_int_poly(f, p):
"""
Create a GF(p)[x] polynomial from Z[x].

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_from_int_poly

>>> gf_from_int_poly([7, -2, 3], 5)
[2, 3, 3]

"""
return gf_trunc(f, p)

[docs]def gf_to_int_poly(f, p, symmetric=True):
"""
Convert a GF(p)[x] polynomial to Z[x].

Examples
========

>>> from sympy.polys.galoistools import gf_to_int_poly

>>> gf_to_int_poly([2, 3, 3], 5)
[2, -2, -2]
>>> gf_to_int_poly([2, 3, 3], 5, symmetric=False)
[2, 3, 3]

"""
if symmetric:
return [ gf_int(c, p) for c in f ]
else:
return f

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_neg

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

>>> from sympy.polys.domains import ZZ

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_sub_ground

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_mul_ground

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_quo_ground

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

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

"""
Add polynomials in GF(p)[x].

Examples
========

>>> from sympy.polys.domains import ZZ

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

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

df = gf_degree(f)
dg = gf_degree(g)

if df == dg:
return gf_strip([ (a + b) % p for a, b in zip(f, g) ])
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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_sub

>>> 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 = gf_degree(f)
dg = gf_degree(g)

if df == dg:
return gf_strip([ (a - b) % p for a, b in zip(f, g) ])
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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_mul

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

"""
df = gf_degree(f)
dg = gf_degree(g)

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

for i in range(0, 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 gf_strip(h)

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_sqr

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

"""
df = gf_degree(f)

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

for i in range(0, 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 gf_strip(h)

[docs]def gf_add_mul(f, g, h, p, K):
"""
Returns f + g*h where f, g, h in GF(p)[x].

Examples
========

>>> from sympy.polys.domains import ZZ
>>> 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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_sub_mul

>>> 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:factor in GF(p)[x].

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_expand

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_div, gf_add_mul

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

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

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

References
==========

1. [Monagan93]_
2. [Gathen99]_

"""
df = gf_degree(f)
dg = gf_degree(g)

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(0, 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], gf_strip(h[dq + 1:])

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_rem

>>> gf_rem(ZZ.map([1, 0, 1, 1]), ZZ.map([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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_quo

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

"""
df = gf_degree(f)
dg = gf_degree(g)

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(0, 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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_exquo

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

>>> gf_exquo(ZZ.map([1, 0, 1, 1]), ZZ.map([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_lshift(f, n, K):
"""
Efficiently multiply f by x**n.

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_lshift

>>> gf_lshift([3, 2, 4], 4, ZZ)
[3, 2, 4, 0, 0, 0, 0]

"""
if not f:
return f
else:
return f + [K.zero]*n

[docs]def gf_rshift(f, n, K):
"""
Efficiently divide f by x**n.

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_rshift

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

"""
if not n:
return f, []
else:
return f[:-n], f[-n:]

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_pow

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

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 = gf_degree(g)

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_frobenius_monomial_base
>>> g = ZZ.map([1, 0, 2, 1])
>>> gf_frobenius_monomial_base(g, 5, ZZ)
[[1], [4, 4, 2], [1, 2]]

"""
n = gf_degree(g)
if n == 0:
return []
b = [0]*n
b[0] = [1]
if p < n:
for i in range(1, n):
mon = gf_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

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_frobenius_monomial_base, gf_frobenius_map
>>> f = ZZ.map([2, 1 , 0, 1])
>>> g = ZZ.map([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 = gf_degree(g)
if gf_degree(f) >= m:
f = gf_rem(f, g, p, K)
if not f:
return []
n = gf_degree(f)
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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_pow_mod

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

References
==========

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_gcd

>>> gf_gcd(ZZ.map([3, 2, 4]), ZZ.map([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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_lcm

>>> gf_lcm(ZZ.map([3, 2, 4]), ZZ.map([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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_cofactors

>>> gf_cofactors(ZZ.map([3, 2, 4]), ZZ.map([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::

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_gcdex, gf_mul, gf_add

>>> s, t, g = gf_gcdex(ZZ.map([1, 8, 7]), ZZ.map([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, ZZ.map([1, 8, 7]), 11, ZZ)
>>> T = gf_mul(t, ZZ.map([1, 7, 1, 7]), 11, ZZ)

>>> gf_add(S, T, 11, ZZ) == [1, 7]
True

References
==========

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_monic

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

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

if K.is_one(lc):
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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_diff

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

"""
df = gf_degree(f)

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 gf_strip(h)

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_eval

>>> 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_multi_eval(f, A, p, K):
"""
Evaluate f(a) for a in [a_1, ..., a_n].

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_multi_eval

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

"""
return [ gf_eval(f, a, p, K) for a in A ]

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_compose

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

"""
if len(g) <= 1:
return gf_strip([gf_eval(f, gf_LC(g, K), p, K)])

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_compose_mod

>>> gf_compose_mod(ZZ.map([3, 2, 4]), ZZ.map([2, 2, 2]), ZZ.map([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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_trace_map

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

References
==========

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_random
>>> gf_random(10, 5, ZZ) #doctest: +SKIP
[1, 2, 3, 2, 1, 1, 1, 2, 0, 4, 2]

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

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_irreducible
>>> 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

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_irred_p_ben_or

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

"""
n = gf_degree(f)

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(0, 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(0, 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

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_irred_p_rabin

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

"""
n = gf_degree(f)

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_irreducible_p

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

"""
method = query('GF_IRRED_METHOD')

if method is not None:
irred = _irred_methods[method](f, p, K)
else:
irred = gf_irred_p_rabin(f, p, K)

return irred

[docs]def gf_sqf_p(f, p, K):
"""
Return True if f is square-free in GF(p)[x].

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_sqf_p

>>> gf_sqf_p(ZZ.map([3, 2, 4]), 5, ZZ)
True
>>> gf_sqf_p(ZZ.map([2, 4, 4, 2, 2, 1, 4]), 5, ZZ)
False

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

if not f:
return True
else:
return gf_gcd(f, gf_diff(f, p, K), p, K) == [K.one]

[docs]def gf_sqf_part(f, p, K):
"""
Return square-free part of a GF(p)[x] polynomial.

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_sqf_part

>>> gf_sqf_part(ZZ.map([1, 1, 3, 0, 1, 0, 2, 2, 1]), 5, ZZ)
[1, 4, 3]

"""
_, sqf = gf_sqf_list(f, p, K)

g = [K.one]

for f, _ in sqf:
g = gf_mul(g, f, p, K)

return g

[docs]def gf_sqf_list(f, p, K, all=False):
"""
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]::

>>> from sympy.polys.domains import ZZ

>>> from sympy.polys.galoistools import (
...     gf_from_dict, gf_diff, gf_sqf_list, gf_pow,
... )
... # doctest: +NORMALIZE_WHITESPACE

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

1. [Geddes92]_

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

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

if gf_degree(f) < 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 gf_degree(H) > 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 = gf_degree(f) // r

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

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

if all:
raise ValueError("'all=True' is not supported yet")

return lc, factors

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

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_Qmatrix

>>> 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 = gf_degree(f), 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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_Qmatrix, gf_Qbasis

>>> 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(0, n):
Q[k][k] = (Q[k][k] - K.one) % p

for k in range(0, 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(0, n):
Q[j][i] = (Q[j][i]*inv) % p

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

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

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

for i in range(0, n):
for j in range(0, 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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_berlekamp

>>> 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] = gf_strip(list(reversed(v)))

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)

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_from_dict

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

Distinct degree factorization gives::

>>> from sympy.polys.galoistools import gf_ddf_zassenhaus

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

1. [Gathen99]_
2. [Geddes92]_

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

b = gf_frobenius_monomial_base(f, p, K)
while 2*i <= gf_degree(f):
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, gf_degree(f))]
else:
return factors

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::

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_edf_zassenhaus

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

References
==========

1. [Gathen99]_
2. [Geddes92]_

"""
factors, q = [f], int(p)

if gf_degree(f) <= n:
return factors

N = gf_degree(f) // 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(0, 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)

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_ddf_shoup, gf_from_dict

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

1. [Kaltofen98]_
2. [Shoup95]_
3. [Gathen92]_

"""
n = gf_degree(f)
k = int(_ceil(_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, gf_degree(f)))

return factors

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_edf_shoup

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

References
==========

1. [Shoup91]_
2. [Gathen92]_

"""
N, q = gf_degree(f), 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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_zassenhaus

>>> gf_zassenhaus(ZZ.map([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
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_shoup

>>> gf_shoup(ZZ.map([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, method=None):
"""
Factor a square-free polynomial f in GF(p)[x].

Examples
========

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_factor_sqf

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

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

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

method = method or query('GF_FACTOR_METHOD')

if method is not None:
factors = _factor_methods[method](f, p, K)
else:
factors = gf_zassenhaus(f, p, K)

return lc, factors

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

>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_factor

>>> gf_factor(ZZ.map([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
==========

1. [Gathen99]_

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

if gf_degree(f) < 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
========

>>> from sympy.polys.galoistools import gf_value

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

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

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

>>> from sympy.polys.galoistools import linear_congruence

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

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

**Reference**
1) Wikipedia http://en.wikipedia.org/wiki/Linear_congruence_theorem

"""
from sympy.polys.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
========

>>> from sympy.polys.galoistools import _raise_mod_power
>>> from sympy.polys.galoistools import csolve_prime

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 sympy.polys.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)

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

Examples
========

>>> from sympy.polys.galoistools import csolve_prime

>>> 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 sympy.polys.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):

>>> from sympy.polys.galoistools import gf_csolve
>>> gf_csolve([1, 1, 7], 189)
[13, 49, 76, 112, 139, 175]

References
==========

[1] 'An introduction to the Theory of Numbers' 5th Edition by Ivan Niven,
Zuckerman and Montgomery.

"""
from sympy.polys.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])