"""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.ntheory import factorint
from sympy.polys.polyconfig import query
from sympy.polys.polyerrors import ExactQuotientFailed
from sympy.polys.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``.
Examples
========
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.
See Also
========
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 ]
[docs]def gf_add_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_add_ground
>>> 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)
[docs]def gf_add(f, g, p, K):
"""
Add polynomials in ``GF(p)[x]``.
Examples
========
>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.galoistools import gf_add
>>> 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
>>> from sympy.polys.galoistools import gf_add_mul
>>> 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 = [f]
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.
References
==========
.. [1] https://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
results returns the final answers.
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])
