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. See Also ======== 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 ]
[docs]def gf_add_ground(f, a, p, K): """ Compute ``f + a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``. Examples ======== >>> gf_add_ground([3, 2, 4], 2, 5, ZZ) [3, 2, 1] """ if not f: a = a % p else: a = (f[-1] + a) % p if len(f) > 1: return f[:-1] + [a] if not a: return [] else: return [a]
[docs]def gf_sub_ground(f, a, p, K): """ Compute ``f - a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``. Examples ======== >>> gf_sub_ground([3, 2, 4], 2, 5, ZZ) [3, 2, 2] """ if not f: a = -a % p else: a = (f[-1] - a) % p if len(f) > 1: return f[:-1] + [a] if not a: return [] else: return [a]
[docs]def gf_mul_ground(f, a, p, K): """ Compute ``f * a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``. Examples ======== >>> gf_mul_ground([3, 2, 4], 2, 5, ZZ) [1, 4, 3] """ if not a: return [] else: return [ (a*b) % p for b in f ]
[docs]def gf_quo_ground(f, a, p, K): """ Compute ``f/a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``. Examples ======== >>> gf_quo_ground([3, 2, 4], 2, 5, ZZ) [4, 1, 2] """ return gf_mul_ground(f, K.invert(a, p), p, K)
[docs]def gf_add(f, g, p, K): """ Add polynomials in ``GF(p)[x]``. Examples ======== >>> gf_add([3, 2, 4], [2, 2, 2], 5, ZZ) [4, 1] """ if not f: return g if not g: return f df = dmp_degree_in(f, 0, 0) dg = dmp_degree_in(g, 0, 0) if df == dg: return dmp_strip([(a + b) % p for a, b in zip(f, g)], 0) else: k = abs(df - dg) if df > dg: h, f = f[:k], f[k:] else: h, g = g[:k], g[k:] return h + [ (a + b) % p for a, b in zip(f, g) ]
[docs]def gf_sub(f, g, p, K): """ Subtract polynomials in ``GF(p)[x]``. Examples ======== >>> gf_sub([3, 2, 4], [2, 2, 2], 5, ZZ) [1, 0, 2] """ if not g: return f if not f: return gf_neg(g, p, K) df = dmp_degree_in(f, 0, 0) dg = dmp_degree_in(g, 0, 0) if df == dg: return dmp_strip([(a - b) % p for a, b in zip(f, g)], 0) else: k = abs(df - dg) if df > dg: h, f = f[:k], f[k:] else: h, g = gf_neg(g[:k], p, K), g[k:] return h + [ (a - b) % p for a, b in zip(f, g) ]
[docs]def gf_mul(f, g, p, K): """ Multiply polynomials in ``GF(p)[x]``. Examples ======== >>> gf_mul([3, 2, 4], [2, 2, 2], 5, ZZ) [1, 0, 3, 2, 3] """ df = dmp_degree_in(f, 0, 0) dg = dmp_degree_in(g, 0, 0) dh = df + dg if dh < 0: return [] h = [0]*(dh + 1) for i in range(dh + 1): coeff = K.zero for j in range(max(0, i - dg), min(i, df) + 1): coeff += f[j]*g[i - j] h[i] = coeff % p return dmp_strip(h, 0)
[docs]def gf_sqr(f, p, K): """ Square polynomials in ``GF(p)[x]``. Examples ======== >>> gf_sqr([3, 2, 4], 5, ZZ) [4, 2, 3, 1, 1] """ df = dmp_degree_in(f, 0, 0) if df < 0: return [] dh = 2*df h = [0]*(dh + 1) for i in range(dh + 1): coeff = K.zero jmin = max(0, i - df) jmax = min(i, df) n = jmax - jmin + 1 jmax = jmin + n // 2 - 1 for j in range(jmin, jmax + 1): coeff += f[j]*f[i - j] coeff += coeff if n & 1: elem = f[jmax + 1] coeff += elem**2 h[i] = coeff % p return dmp_strip(h, 0)
[docs]def gf_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) >>> gf_add(S, T, 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 results returns the final answers. 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)