Coverage for /usr/lib/python3/dist-packages/sympy/polys/galoistools.py: 11%
721 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Dense univariate polynomials with coefficients in Galois fields. """
3from math import ceil as _ceil, sqrt as _sqrt, prod
5from sympy.core.random import uniform
6from sympy.external.gmpy import SYMPY_INTS
7from sympy.polys.polyconfig import query
8from sympy.polys.polyerrors import ExactQuotientFailed
9from sympy.polys.polyutils import _sort_factors
12def gf_crt(U, M, K=None):
13 """
14 Chinese Remainder Theorem.
16 Given a set of integer residues ``u_0,...,u_n`` and a set of
17 co-prime integer moduli ``m_0,...,m_n``, returns an integer
18 ``u``, such that ``u = u_i mod m_i`` for ``i = ``0,...,n``.
20 Examples
21 ========
23 Consider a set of residues ``U = [49, 76, 65]``
24 and a set of moduli ``M = [99, 97, 95]``. Then we have::
26 >>> from sympy.polys.domains import ZZ
27 >>> from sympy.polys.galoistools import gf_crt
29 >>> gf_crt([49, 76, 65], [99, 97, 95], ZZ)
30 639985
32 This is the correct result because::
34 >>> [639985 % m for m in [99, 97, 95]]
35 [49, 76, 65]
37 Note: this is a low-level routine with no error checking.
39 See Also
40 ========
42 sympy.ntheory.modular.crt : a higher level crt routine
43 sympy.ntheory.modular.solve_congruence
45 """
46 p = prod(M, start=K.one)
47 v = K.zero
49 for u, m in zip(U, M):
50 e = p // m
51 s, _, _ = K.gcdex(e, m)
52 v += e*(u*s % m)
54 return v % p
57def gf_crt1(M, K):
58 """
59 First part of the Chinese Remainder Theorem.
61 Examples
62 ========
64 >>> from sympy.polys.domains import ZZ
65 >>> from sympy.polys.galoistools import gf_crt1
67 >>> gf_crt1([99, 97, 95], ZZ)
68 (912285, [9215, 9405, 9603], [62, 24, 12])
70 """
71 E, S = [], []
72 p = prod(M, start=K.one)
74 for m in M:
75 E.append(p // m)
76 S.append(K.gcdex(E[-1], m)[0] % m)
78 return p, E, S
81def gf_crt2(U, M, p, E, S, K):
82 """
83 Second part of the Chinese Remainder Theorem.
85 Examples
86 ========
88 >>> from sympy.polys.domains import ZZ
89 >>> from sympy.polys.galoistools import gf_crt2
91 >>> U = [49, 76, 65]
92 >>> M = [99, 97, 95]
93 >>> p = 912285
94 >>> E = [9215, 9405, 9603]
95 >>> S = [62, 24, 12]
97 >>> gf_crt2(U, M, p, E, S, ZZ)
98 639985
100 """
101 v = K.zero
103 for u, m, e, s in zip(U, M, E, S):
104 v += e*(u*s % m)
106 return v % p
109def gf_int(a, p):
110 """
111 Coerce ``a mod p`` to an integer in the range ``[-p/2, p/2]``.
113 Examples
114 ========
116 >>> from sympy.polys.galoistools import gf_int
118 >>> gf_int(2, 7)
119 2
120 >>> gf_int(5, 7)
121 -2
123 """
124 if a <= p // 2:
125 return a
126 else:
127 return a - p
130def gf_degree(f):
131 """
132 Return the leading degree of ``f``.
134 Examples
135 ========
137 >>> from sympy.polys.galoistools import gf_degree
139 >>> gf_degree([1, 1, 2, 0])
140 3
141 >>> gf_degree([])
142 -1
144 """
145 return len(f) - 1
148def gf_LC(f, K):
149 """
150 Return the leading coefficient of ``f``.
152 Examples
153 ========
155 >>> from sympy.polys.domains import ZZ
156 >>> from sympy.polys.galoistools import gf_LC
158 >>> gf_LC([3, 0, 1], ZZ)
159 3
161 """
162 if not f:
163 return K.zero
164 else:
165 return f[0]
168def gf_TC(f, K):
169 """
170 Return the trailing coefficient of ``f``.
172 Examples
173 ========
175 >>> from sympy.polys.domains import ZZ
176 >>> from sympy.polys.galoistools import gf_TC
178 >>> gf_TC([3, 0, 1], ZZ)
179 1
181 """
182 if not f:
183 return K.zero
184 else:
185 return f[-1]
188def gf_strip(f):
189 """
190 Remove leading zeros from ``f``.
193 Examples
194 ========
196 >>> from sympy.polys.galoistools import gf_strip
198 >>> gf_strip([0, 0, 0, 3, 0, 1])
199 [3, 0, 1]
201 """
202 if not f or f[0]:
203 return f
205 k = 0
207 for coeff in f:
208 if coeff:
209 break
210 else:
211 k += 1
213 return f[k:]
216def gf_trunc(f, p):
217 """
218 Reduce all coefficients modulo ``p``.
220 Examples
221 ========
223 >>> from sympy.polys.galoistools import gf_trunc
225 >>> gf_trunc([7, -2, 3], 5)
226 [2, 3, 3]
228 """
229 return gf_strip([ a % p for a in f ])
232def gf_normal(f, p, K):
233 """
234 Normalize all coefficients in ``K``.
236 Examples
237 ========
239 >>> from sympy.polys.domains import ZZ
240 >>> from sympy.polys.galoistools import gf_normal
242 >>> gf_normal([5, 10, 21, -3], 5, ZZ)
243 [1, 2]
245 """
246 return gf_trunc(list(map(K, f)), p)
249def gf_from_dict(f, p, K):
250 """
251 Create a ``GF(p)[x]`` polynomial from a dict.
253 Examples
254 ========
256 >>> from sympy.polys.domains import ZZ
257 >>> from sympy.polys.galoistools import gf_from_dict
259 >>> gf_from_dict({10: ZZ(4), 4: ZZ(33), 0: ZZ(-1)}, 5, ZZ)
260 [4, 0, 0, 0, 0, 0, 3, 0, 0, 0, 4]
262 """
263 n, h = max(f.keys()), []
265 if isinstance(n, SYMPY_INTS):
266 for k in range(n, -1, -1):
267 h.append(f.get(k, K.zero) % p)
268 else:
269 (n,) = n
271 for k in range(n, -1, -1):
272 h.append(f.get((k,), K.zero) % p)
274 return gf_trunc(h, p)
277def gf_to_dict(f, p, symmetric=True):
278 """
279 Convert a ``GF(p)[x]`` polynomial to a dict.
281 Examples
282 ========
284 >>> from sympy.polys.galoistools import gf_to_dict
286 >>> gf_to_dict([4, 0, 0, 0, 0, 0, 3, 0, 0, 0, 4], 5)
287 {0: -1, 4: -2, 10: -1}
288 >>> gf_to_dict([4, 0, 0, 0, 0, 0, 3, 0, 0, 0, 4], 5, symmetric=False)
289 {0: 4, 4: 3, 10: 4}
291 """
292 n, result = gf_degree(f), {}
294 for k in range(0, n + 1):
295 if symmetric:
296 a = gf_int(f[n - k], p)
297 else:
298 a = f[n - k]
300 if a:
301 result[k] = a
303 return result
306def gf_from_int_poly(f, p):
307 """
308 Create a ``GF(p)[x]`` polynomial from ``Z[x]``.
310 Examples
311 ========
313 >>> from sympy.polys.galoistools import gf_from_int_poly
315 >>> gf_from_int_poly([7, -2, 3], 5)
316 [2, 3, 3]
318 """
319 return gf_trunc(f, p)
322def gf_to_int_poly(f, p, symmetric=True):
323 """
324 Convert a ``GF(p)[x]`` polynomial to ``Z[x]``.
327 Examples
328 ========
330 >>> from sympy.polys.galoistools import gf_to_int_poly
332 >>> gf_to_int_poly([2, 3, 3], 5)
333 [2, -2, -2]
334 >>> gf_to_int_poly([2, 3, 3], 5, symmetric=False)
335 [2, 3, 3]
337 """
338 if symmetric:
339 return [ gf_int(c, p) for c in f ]
340 else:
341 return f
344def gf_neg(f, p, K):
345 """
346 Negate a polynomial in ``GF(p)[x]``.
348 Examples
349 ========
351 >>> from sympy.polys.domains import ZZ
352 >>> from sympy.polys.galoistools import gf_neg
354 >>> gf_neg([3, 2, 1, 0], 5, ZZ)
355 [2, 3, 4, 0]
357 """
358 return [ -coeff % p for coeff in f ]
361def gf_add_ground(f, a, p, K):
362 """
363 Compute ``f + a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``.
365 Examples
366 ========
368 >>> from sympy.polys.domains import ZZ
369 >>> from sympy.polys.galoistools import gf_add_ground
371 >>> gf_add_ground([3, 2, 4], 2, 5, ZZ)
372 [3, 2, 1]
374 """
375 if not f:
376 a = a % p
377 else:
378 a = (f[-1] + a) % p
380 if len(f) > 1:
381 return f[:-1] + [a]
383 if not a:
384 return []
385 else:
386 return [a]
389def gf_sub_ground(f, a, p, K):
390 """
391 Compute ``f - a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``.
393 Examples
394 ========
396 >>> from sympy.polys.domains import ZZ
397 >>> from sympy.polys.galoistools import gf_sub_ground
399 >>> gf_sub_ground([3, 2, 4], 2, 5, ZZ)
400 [3, 2, 2]
402 """
403 if not f:
404 a = -a % p
405 else:
406 a = (f[-1] - a) % p
408 if len(f) > 1:
409 return f[:-1] + [a]
411 if not a:
412 return []
413 else:
414 return [a]
417def gf_mul_ground(f, a, p, K):
418 """
419 Compute ``f * a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``.
421 Examples
422 ========
424 >>> from sympy.polys.domains import ZZ
425 >>> from sympy.polys.galoistools import gf_mul_ground
427 >>> gf_mul_ground([3, 2, 4], 2, 5, ZZ)
428 [1, 4, 3]
430 """
431 if not a:
432 return []
433 else:
434 return [ (a*b) % p for b in f ]
437def gf_quo_ground(f, a, p, K):
438 """
439 Compute ``f/a`` where ``f`` in ``GF(p)[x]`` and ``a`` in ``GF(p)``.
441 Examples
442 ========
444 >>> from sympy.polys.domains import ZZ
445 >>> from sympy.polys.galoistools import gf_quo_ground
447 >>> gf_quo_ground(ZZ.map([3, 2, 4]), ZZ(2), 5, ZZ)
448 [4, 1, 2]
450 """
451 return gf_mul_ground(f, K.invert(a, p), p, K)
454def gf_add(f, g, p, K):
455 """
456 Add polynomials in ``GF(p)[x]``.
458 Examples
459 ========
461 >>> from sympy.polys.domains import ZZ
462 >>> from sympy.polys.galoistools import gf_add
464 >>> gf_add([3, 2, 4], [2, 2, 2], 5, ZZ)
465 [4, 1]
467 """
468 if not f:
469 return g
470 if not g:
471 return f
473 df = gf_degree(f)
474 dg = gf_degree(g)
476 if df == dg:
477 return gf_strip([ (a + b) % p for a, b in zip(f, g) ])
478 else:
479 k = abs(df - dg)
481 if df > dg:
482 h, f = f[:k], f[k:]
483 else:
484 h, g = g[:k], g[k:]
486 return h + [ (a + b) % p for a, b in zip(f, g) ]
489def gf_sub(f, g, p, K):
490 """
491 Subtract polynomials in ``GF(p)[x]``.
493 Examples
494 ========
496 >>> from sympy.polys.domains import ZZ
497 >>> from sympy.polys.galoistools import gf_sub
499 >>> gf_sub([3, 2, 4], [2, 2, 2], 5, ZZ)
500 [1, 0, 2]
502 """
503 if not g:
504 return f
505 if not f:
506 return gf_neg(g, p, K)
508 df = gf_degree(f)
509 dg = gf_degree(g)
511 if df == dg:
512 return gf_strip([ (a - b) % p for a, b in zip(f, g) ])
513 else:
514 k = abs(df - dg)
516 if df > dg:
517 h, f = f[:k], f[k:]
518 else:
519 h, g = gf_neg(g[:k], p, K), g[k:]
521 return h + [ (a - b) % p for a, b in zip(f, g) ]
524def gf_mul(f, g, p, K):
525 """
526 Multiply polynomials in ``GF(p)[x]``.
528 Examples
529 ========
531 >>> from sympy.polys.domains import ZZ
532 >>> from sympy.polys.galoistools import gf_mul
534 >>> gf_mul([3, 2, 4], [2, 2, 2], 5, ZZ)
535 [1, 0, 3, 2, 3]
537 """
538 df = gf_degree(f)
539 dg = gf_degree(g)
541 dh = df + dg
542 h = [0]*(dh + 1)
544 for i in range(0, dh + 1):
545 coeff = K.zero
547 for j in range(max(0, i - dg), min(i, df) + 1):
548 coeff += f[j]*g[i - j]
550 h[i] = coeff % p
552 return gf_strip(h)
555def gf_sqr(f, p, K):
556 """
557 Square polynomials in ``GF(p)[x]``.
559 Examples
560 ========
562 >>> from sympy.polys.domains import ZZ
563 >>> from sympy.polys.galoistools import gf_sqr
565 >>> gf_sqr([3, 2, 4], 5, ZZ)
566 [4, 2, 3, 1, 1]
568 """
569 df = gf_degree(f)
571 dh = 2*df
572 h = [0]*(dh + 1)
574 for i in range(0, dh + 1):
575 coeff = K.zero
577 jmin = max(0, i - df)
578 jmax = min(i, df)
580 n = jmax - jmin + 1
582 jmax = jmin + n // 2 - 1
584 for j in range(jmin, jmax + 1):
585 coeff += f[j]*f[i - j]
587 coeff += coeff
589 if n & 1:
590 elem = f[jmax + 1]
591 coeff += elem**2
593 h[i] = coeff % p
595 return gf_strip(h)
598def gf_add_mul(f, g, h, p, K):
599 """
600 Returns ``f + g*h`` where ``f``, ``g``, ``h`` in ``GF(p)[x]``.
602 Examples
603 ========
605 >>> from sympy.polys.domains import ZZ
606 >>> from sympy.polys.galoistools import gf_add_mul
607 >>> gf_add_mul([3, 2, 4], [2, 2, 2], [1, 4], 5, ZZ)
608 [2, 3, 2, 2]
609 """
610 return gf_add(f, gf_mul(g, h, p, K), p, K)
613def gf_sub_mul(f, g, h, p, K):
614 """
615 Compute ``f - g*h`` where ``f``, ``g``, ``h`` in ``GF(p)[x]``.
617 Examples
618 ========
620 >>> from sympy.polys.domains import ZZ
621 >>> from sympy.polys.galoistools import gf_sub_mul
623 >>> gf_sub_mul([3, 2, 4], [2, 2, 2], [1, 4], 5, ZZ)
624 [3, 3, 2, 1]
626 """
627 return gf_sub(f, gf_mul(g, h, p, K), p, K)
630def gf_expand(F, p, K):
631 """
632 Expand results of :func:`~.factor` in ``GF(p)[x]``.
634 Examples
635 ========
637 >>> from sympy.polys.domains import ZZ
638 >>> from sympy.polys.galoistools import gf_expand
640 >>> gf_expand([([3, 2, 4], 1), ([2, 2], 2), ([3, 1], 3)], 5, ZZ)
641 [4, 3, 0, 3, 0, 1, 4, 1]
643 """
644 if isinstance(F, tuple):
645 lc, F = F
646 else:
647 lc = K.one
649 g = [lc]
651 for f, k in F:
652 f = gf_pow(f, k, p, K)
653 g = gf_mul(g, f, p, K)
655 return g
658def gf_div(f, g, p, K):
659 """
660 Division with remainder in ``GF(p)[x]``.
662 Given univariate polynomials ``f`` and ``g`` with coefficients in a
663 finite field with ``p`` elements, returns polynomials ``q`` and ``r``
664 (quotient and remainder) such that ``f = q*g + r``.
666 Consider polynomials ``x**3 + x + 1`` and ``x**2 + x`` in GF(2)::
668 >>> from sympy.polys.domains import ZZ
669 >>> from sympy.polys.galoistools import gf_div, gf_add_mul
671 >>> gf_div(ZZ.map([1, 0, 1, 1]), ZZ.map([1, 1, 0]), 2, ZZ)
672 ([1, 1], [1])
674 As result we obtained quotient ``x + 1`` and remainder ``1``, thus::
676 >>> gf_add_mul(ZZ.map([1]), ZZ.map([1, 1]), ZZ.map([1, 1, 0]), 2, ZZ)
677 [1, 0, 1, 1]
679 References
680 ==========
682 .. [1] [Monagan93]_
683 .. [2] [Gathen99]_
685 """
686 df = gf_degree(f)
687 dg = gf_degree(g)
689 if not g:
690 raise ZeroDivisionError("polynomial division")
691 elif df < dg:
692 return [], f
694 inv = K.invert(g[0], p)
696 h, dq, dr = list(f), df - dg, dg - 1
698 for i in range(0, df + 1):
699 coeff = h[i]
701 for j in range(max(0, dg - i), min(df - i, dr) + 1):
702 coeff -= h[i + j - dg] * g[dg - j]
704 if i <= dq:
705 coeff *= inv
707 h[i] = coeff % p
709 return h[:dq + 1], gf_strip(h[dq + 1:])
712def gf_rem(f, g, p, K):
713 """
714 Compute polynomial remainder in ``GF(p)[x]``.
716 Examples
717 ========
719 >>> from sympy.polys.domains import ZZ
720 >>> from sympy.polys.galoistools import gf_rem
722 >>> gf_rem(ZZ.map([1, 0, 1, 1]), ZZ.map([1, 1, 0]), 2, ZZ)
723 [1]
725 """
726 return gf_div(f, g, p, K)[1]
729def gf_quo(f, g, p, K):
730 """
731 Compute exact quotient in ``GF(p)[x]``.
733 Examples
734 ========
736 >>> from sympy.polys.domains import ZZ
737 >>> from sympy.polys.galoistools import gf_quo
739 >>> gf_quo(ZZ.map([1, 0, 1, 1]), ZZ.map([1, 1, 0]), 2, ZZ)
740 [1, 1]
741 >>> gf_quo(ZZ.map([1, 0, 3, 2, 3]), ZZ.map([2, 2, 2]), 5, ZZ)
742 [3, 2, 4]
744 """
745 df = gf_degree(f)
746 dg = gf_degree(g)
748 if not g:
749 raise ZeroDivisionError("polynomial division")
750 elif df < dg:
751 return []
753 inv = K.invert(g[0], p)
755 h, dq, dr = f[:], df - dg, dg - 1
757 for i in range(0, dq + 1):
758 coeff = h[i]
760 for j in range(max(0, dg - i), min(df - i, dr) + 1):
761 coeff -= h[i + j - dg] * g[dg - j]
763 h[i] = (coeff * inv) % p
765 return h[:dq + 1]
768def gf_exquo(f, g, p, K):
769 """
770 Compute polynomial quotient in ``GF(p)[x]``.
772 Examples
773 ========
775 >>> from sympy.polys.domains import ZZ
776 >>> from sympy.polys.galoistools import gf_exquo
778 >>> gf_exquo(ZZ.map([1, 0, 3, 2, 3]), ZZ.map([2, 2, 2]), 5, ZZ)
779 [3, 2, 4]
781 >>> gf_exquo(ZZ.map([1, 0, 1, 1]), ZZ.map([1, 1, 0]), 2, ZZ)
782 Traceback (most recent call last):
783 ...
784 ExactQuotientFailed: [1, 1, 0] does not divide [1, 0, 1, 1]
786 """
787 q, r = gf_div(f, g, p, K)
789 if not r:
790 return q
791 else:
792 raise ExactQuotientFailed(f, g)
795def gf_lshift(f, n, K):
796 """
797 Efficiently multiply ``f`` by ``x**n``.
799 Examples
800 ========
802 >>> from sympy.polys.domains import ZZ
803 >>> from sympy.polys.galoistools import gf_lshift
805 >>> gf_lshift([3, 2, 4], 4, ZZ)
806 [3, 2, 4, 0, 0, 0, 0]
808 """
809 if not f:
810 return f
811 else:
812 return f + [K.zero]*n
815def gf_rshift(f, n, K):
816 """
817 Efficiently divide ``f`` by ``x**n``.
819 Examples
820 ========
822 >>> from sympy.polys.domains import ZZ
823 >>> from sympy.polys.galoistools import gf_rshift
825 >>> gf_rshift([1, 2, 3, 4, 0], 3, ZZ)
826 ([1, 2], [3, 4, 0])
828 """
829 if not n:
830 return f, []
831 else:
832 return f[:-n], f[-n:]
835def gf_pow(f, n, p, K):
836 """
837 Compute ``f**n`` in ``GF(p)[x]`` using repeated squaring.
839 Examples
840 ========
842 >>> from sympy.polys.domains import ZZ
843 >>> from sympy.polys.galoistools import gf_pow
845 >>> gf_pow([3, 2, 4], 3, 5, ZZ)
846 [2, 4, 4, 2, 2, 1, 4]
848 """
849 if not n:
850 return [K.one]
851 elif n == 1:
852 return f
853 elif n == 2:
854 return gf_sqr(f, p, K)
856 h = [K.one]
858 while True:
859 if n & 1:
860 h = gf_mul(h, f, p, K)
861 n -= 1
863 n >>= 1
865 if not n:
866 break
868 f = gf_sqr(f, p, K)
870 return h
872def gf_frobenius_monomial_base(g, p, K):
873 """
874 return the list of ``x**(i*p) mod g in Z_p`` for ``i = 0, .., n - 1``
875 where ``n = gf_degree(g)``
877 Examples
878 ========
880 >>> from sympy.polys.domains import ZZ
881 >>> from sympy.polys.galoistools import gf_frobenius_monomial_base
882 >>> g = ZZ.map([1, 0, 2, 1])
883 >>> gf_frobenius_monomial_base(g, 5, ZZ)
884 [[1], [4, 4, 2], [1, 2]]
886 """
887 n = gf_degree(g)
888 if n == 0:
889 return []
890 b = [0]*n
891 b[0] = [1]
892 if p < n:
893 for i in range(1, n):
894 mon = gf_lshift(b[i - 1], p, K)
895 b[i] = gf_rem(mon, g, p, K)
896 elif n > 1:
897 b[1] = gf_pow_mod([K.one, K.zero], p, g, p, K)
898 for i in range(2, n):
899 b[i] = gf_mul(b[i - 1], b[1], p, K)
900 b[i] = gf_rem(b[i], g, p, K)
902 return b
904def gf_frobenius_map(f, g, b, p, K):
905 """
906 compute gf_pow_mod(f, p, g, p, K) using the Frobenius map
908 Parameters
909 ==========
911 f, g : polynomials in ``GF(p)[x]``
912 b : frobenius monomial base
913 p : prime number
914 K : domain
916 Examples
917 ========
919 >>> from sympy.polys.domains import ZZ
920 >>> from sympy.polys.galoistools import gf_frobenius_monomial_base, gf_frobenius_map
921 >>> f = ZZ.map([2, 1, 0, 1])
922 >>> g = ZZ.map([1, 0, 2, 1])
923 >>> p = 5
924 >>> b = gf_frobenius_monomial_base(g, p, ZZ)
925 >>> r = gf_frobenius_map(f, g, b, p, ZZ)
926 >>> gf_frobenius_map(f, g, b, p, ZZ)
927 [4, 0, 3]
928 """
929 m = gf_degree(g)
930 if gf_degree(f) >= m:
931 f = gf_rem(f, g, p, K)
932 if not f:
933 return []
934 n = gf_degree(f)
935 sf = [f[-1]]
936 for i in range(1, n + 1):
937 v = gf_mul_ground(b[i], f[n - i], p, K)
938 sf = gf_add(sf, v, p, K)
939 return sf
941def _gf_pow_pnm1d2(f, n, g, b, p, K):
942 """
943 utility function for ``gf_edf_zassenhaus``
944 Compute ``f**((p**n - 1) // 2)`` in ``GF(p)[x]/(g)``
945 ``f**((p**n - 1) // 2) = (f*f**p*...*f**(p**n - 1))**((p - 1) // 2)``
946 """
947 f = gf_rem(f, g, p, K)
948 h = f
949 r = f
950 for i in range(1, n):
951 h = gf_frobenius_map(h, g, b, p, K)
952 r = gf_mul(r, h, p, K)
953 r = gf_rem(r, g, p, K)
955 res = gf_pow_mod(r, (p - 1)//2, g, p, K)
956 return res
958def gf_pow_mod(f, n, g, p, K):
959 """
960 Compute ``f**n`` in ``GF(p)[x]/(g)`` using repeated squaring.
962 Given polynomials ``f`` and ``g`` in ``GF(p)[x]`` and a non-negative
963 integer ``n``, efficiently computes ``f**n (mod g)`` i.e. the remainder
964 of ``f**n`` from division by ``g``, using the repeated squaring algorithm.
966 Examples
967 ========
969 >>> from sympy.polys.domains import ZZ
970 >>> from sympy.polys.galoistools import gf_pow_mod
972 >>> gf_pow_mod(ZZ.map([3, 2, 4]), 3, ZZ.map([1, 1]), 5, ZZ)
973 []
975 References
976 ==========
978 .. [1] [Gathen99]_
980 """
981 if not n:
982 return [K.one]
983 elif n == 1:
984 return gf_rem(f, g, p, K)
985 elif n == 2:
986 return gf_rem(gf_sqr(f, p, K), g, p, K)
988 h = [K.one]
990 while True:
991 if n & 1:
992 h = gf_mul(h, f, p, K)
993 h = gf_rem(h, g, p, K)
994 n -= 1
996 n >>= 1
998 if not n:
999 break
1001 f = gf_sqr(f, p, K)
1002 f = gf_rem(f, g, p, K)
1004 return h
1007def gf_gcd(f, g, p, K):
1008 """
1009 Euclidean Algorithm in ``GF(p)[x]``.
1011 Examples
1012 ========
1014 >>> from sympy.polys.domains import ZZ
1015 >>> from sympy.polys.galoistools import gf_gcd
1017 >>> gf_gcd(ZZ.map([3, 2, 4]), ZZ.map([2, 2, 3]), 5, ZZ)
1018 [1, 3]
1020 """
1021 while g:
1022 f, g = g, gf_rem(f, g, p, K)
1024 return gf_monic(f, p, K)[1]
1027def gf_lcm(f, g, p, K):
1028 """
1029 Compute polynomial LCM in ``GF(p)[x]``.
1031 Examples
1032 ========
1034 >>> from sympy.polys.domains import ZZ
1035 >>> from sympy.polys.galoistools import gf_lcm
1037 >>> gf_lcm(ZZ.map([3, 2, 4]), ZZ.map([2, 2, 3]), 5, ZZ)
1038 [1, 2, 0, 4]
1040 """
1041 if not f or not g:
1042 return []
1044 h = gf_quo(gf_mul(f, g, p, K),
1045 gf_gcd(f, g, p, K), p, K)
1047 return gf_monic(h, p, K)[1]
1050def gf_cofactors(f, g, p, K):
1051 """
1052 Compute polynomial GCD and cofactors in ``GF(p)[x]``.
1054 Examples
1055 ========
1057 >>> from sympy.polys.domains import ZZ
1058 >>> from sympy.polys.galoistools import gf_cofactors
1060 >>> gf_cofactors(ZZ.map([3, 2, 4]), ZZ.map([2, 2, 3]), 5, ZZ)
1061 ([1, 3], [3, 3], [2, 1])
1063 """
1064 if not f and not g:
1065 return ([], [], [])
1067 h = gf_gcd(f, g, p, K)
1069 return (h, gf_quo(f, h, p, K),
1070 gf_quo(g, h, p, K))
1073def gf_gcdex(f, g, p, K):
1074 """
1075 Extended Euclidean Algorithm in ``GF(p)[x]``.
1077 Given polynomials ``f`` and ``g`` in ``GF(p)[x]``, computes polynomials
1078 ``s``, ``t`` and ``h``, such that ``h = gcd(f, g)`` and ``s*f + t*g = h``.
1079 The typical application of EEA is solving polynomial diophantine equations.
1081 Consider polynomials ``f = (x + 7) (x + 1)``, ``g = (x + 7) (x**2 + 1)``
1082 in ``GF(11)[x]``. Application of Extended Euclidean Algorithm gives::
1084 >>> from sympy.polys.domains import ZZ
1085 >>> from sympy.polys.galoistools import gf_gcdex, gf_mul, gf_add
1087 >>> s, t, g = gf_gcdex(ZZ.map([1, 8, 7]), ZZ.map([1, 7, 1, 7]), 11, ZZ)
1088 >>> s, t, g
1089 ([5, 6], [6], [1, 7])
1091 As result we obtained polynomials ``s = 5*x + 6`` and ``t = 6``, and
1092 additionally ``gcd(f, g) = x + 7``. This is correct because::
1094 >>> S = gf_mul(s, ZZ.map([1, 8, 7]), 11, ZZ)
1095 >>> T = gf_mul(t, ZZ.map([1, 7, 1, 7]), 11, ZZ)
1097 >>> gf_add(S, T, 11, ZZ) == [1, 7]
1098 True
1100 References
1101 ==========
1103 .. [1] [Gathen99]_
1105 """
1106 if not (f or g):
1107 return [K.one], [], []
1109 p0, r0 = gf_monic(f, p, K)
1110 p1, r1 = gf_monic(g, p, K)
1112 if not f:
1113 return [], [K.invert(p1, p)], r1
1114 if not g:
1115 return [K.invert(p0, p)], [], r0
1117 s0, s1 = [K.invert(p0, p)], []
1118 t0, t1 = [], [K.invert(p1, p)]
1120 while True:
1121 Q, R = gf_div(r0, r1, p, K)
1123 if not R:
1124 break
1126 (lc, r1), r0 = gf_monic(R, p, K), r1
1128 inv = K.invert(lc, p)
1130 s = gf_sub_mul(s0, s1, Q, p, K)
1131 t = gf_sub_mul(t0, t1, Q, p, K)
1133 s1, s0 = gf_mul_ground(s, inv, p, K), s1
1134 t1, t0 = gf_mul_ground(t, inv, p, K), t1
1136 return s1, t1, r1
1139def gf_monic(f, p, K):
1140 """
1141 Compute LC and a monic polynomial in ``GF(p)[x]``.
1143 Examples
1144 ========
1146 >>> from sympy.polys.domains import ZZ
1147 >>> from sympy.polys.galoistools import gf_monic
1149 >>> gf_monic(ZZ.map([3, 2, 4]), 5, ZZ)
1150 (3, [1, 4, 3])
1152 """
1153 if not f:
1154 return K.zero, []
1155 else:
1156 lc = f[0]
1158 if K.is_one(lc):
1159 return lc, list(f)
1160 else:
1161 return lc, gf_quo_ground(f, lc, p, K)
1164def gf_diff(f, p, K):
1165 """
1166 Differentiate polynomial in ``GF(p)[x]``.
1168 Examples
1169 ========
1171 >>> from sympy.polys.domains import ZZ
1172 >>> from sympy.polys.galoistools import gf_diff
1174 >>> gf_diff([3, 2, 4], 5, ZZ)
1175 [1, 2]
1177 """
1178 df = gf_degree(f)
1180 h, n = [K.zero]*df, df
1182 for coeff in f[:-1]:
1183 coeff *= K(n)
1184 coeff %= p
1186 if coeff:
1187 h[df - n] = coeff
1189 n -= 1
1191 return gf_strip(h)
1194def gf_eval(f, a, p, K):
1195 """
1196 Evaluate ``f(a)`` in ``GF(p)`` using Horner scheme.
1198 Examples
1199 ========
1201 >>> from sympy.polys.domains import ZZ
1202 >>> from sympy.polys.galoistools import gf_eval
1204 >>> gf_eval([3, 2, 4], 2, 5, ZZ)
1205 0
1207 """
1208 result = K.zero
1210 for c in f:
1211 result *= a
1212 result += c
1213 result %= p
1215 return result
1218def gf_multi_eval(f, A, p, K):
1219 """
1220 Evaluate ``f(a)`` for ``a`` in ``[a_1, ..., a_n]``.
1222 Examples
1223 ========
1225 >>> from sympy.polys.domains import ZZ
1226 >>> from sympy.polys.galoistools import gf_multi_eval
1228 >>> gf_multi_eval([3, 2, 4], [0, 1, 2, 3, 4], 5, ZZ)
1229 [4, 4, 0, 2, 0]
1231 """
1232 return [ gf_eval(f, a, p, K) for a in A ]
1235def gf_compose(f, g, p, K):
1236 """
1237 Compute polynomial composition ``f(g)`` in ``GF(p)[x]``.
1239 Examples
1240 ========
1242 >>> from sympy.polys.domains import ZZ
1243 >>> from sympy.polys.galoistools import gf_compose
1245 >>> gf_compose([3, 2, 4], [2, 2, 2], 5, ZZ)
1246 [2, 4, 0, 3, 0]
1248 """
1249 if len(g) <= 1:
1250 return gf_strip([gf_eval(f, gf_LC(g, K), p, K)])
1252 if not f:
1253 return []
1255 h = [f[0]]
1257 for c in f[1:]:
1258 h = gf_mul(h, g, p, K)
1259 h = gf_add_ground(h, c, p, K)
1261 return h
1264def gf_compose_mod(g, h, f, p, K):
1265 """
1266 Compute polynomial composition ``g(h)`` in ``GF(p)[x]/(f)``.
1268 Examples
1269 ========
1271 >>> from sympy.polys.domains import ZZ
1272 >>> from sympy.polys.galoistools import gf_compose_mod
1274 >>> gf_compose_mod(ZZ.map([3, 2, 4]), ZZ.map([2, 2, 2]), ZZ.map([4, 3]), 5, ZZ)
1275 [4]
1277 """
1278 if not g:
1279 return []
1281 comp = [g[0]]
1283 for a in g[1:]:
1284 comp = gf_mul(comp, h, p, K)
1285 comp = gf_add_ground(comp, a, p, K)
1286 comp = gf_rem(comp, f, p, K)
1288 return comp
1291def gf_trace_map(a, b, c, n, f, p, K):
1292 """
1293 Compute polynomial trace map in ``GF(p)[x]/(f)``.
1295 Given a polynomial ``f`` in ``GF(p)[x]``, polynomials ``a``, ``b``,
1296 ``c`` in the quotient ring ``GF(p)[x]/(f)`` such that ``b = c**t
1297 (mod f)`` for some positive power ``t`` of ``p``, and a positive
1298 integer ``n``, returns a mapping::
1300 a -> a**t**n, a + a**t + a**t**2 + ... + a**t**n (mod f)
1302 In factorization context, ``b = x**p mod f`` and ``c = x mod f``.
1303 This way we can efficiently compute trace polynomials in equal
1304 degree factorization routine, much faster than with other methods,
1305 like iterated Frobenius algorithm, for large degrees.
1307 Examples
1308 ========
1310 >>> from sympy.polys.domains import ZZ
1311 >>> from sympy.polys.galoistools import gf_trace_map
1313 >>> gf_trace_map([1, 2], [4, 4], [1, 1], 4, [3, 2, 4], 5, ZZ)
1314 ([1, 3], [1, 3])
1316 References
1317 ==========
1319 .. [1] [Gathen92]_
1321 """
1322 u = gf_compose_mod(a, b, f, p, K)
1323 v = b
1325 if n & 1:
1326 U = gf_add(a, u, p, K)
1327 V = b
1328 else:
1329 U = a
1330 V = c
1332 n >>= 1
1334 while n:
1335 u = gf_add(u, gf_compose_mod(u, v, f, p, K), p, K)
1336 v = gf_compose_mod(v, v, f, p, K)
1338 if n & 1:
1339 U = gf_add(U, gf_compose_mod(u, V, f, p, K), p, K)
1340 V = gf_compose_mod(v, V, f, p, K)
1342 n >>= 1
1344 return gf_compose_mod(a, V, f, p, K), U
1346def _gf_trace_map(f, n, g, b, p, K):
1347 """
1348 utility for ``gf_edf_shoup``
1349 """
1350 f = gf_rem(f, g, p, K)
1351 h = f
1352 r = f
1353 for i in range(1, n):
1354 h = gf_frobenius_map(h, g, b, p, K)
1355 r = gf_add(r, h, p, K)
1356 r = gf_rem(r, g, p, K)
1357 return r
1360def gf_random(n, p, K):
1361 """
1362 Generate a random polynomial in ``GF(p)[x]`` of degree ``n``.
1364 Examples
1365 ========
1367 >>> from sympy.polys.domains import ZZ
1368 >>> from sympy.polys.galoistools import gf_random
1369 >>> gf_random(10, 5, ZZ) #doctest: +SKIP
1370 [1, 2, 3, 2, 1, 1, 1, 2, 0, 4, 2]
1372 """
1373 return [K.one] + [ K(int(uniform(0, p))) for i in range(0, n) ]
1376def gf_irreducible(n, p, K):
1377 """
1378 Generate random irreducible polynomial of degree ``n`` in ``GF(p)[x]``.
1380 Examples
1381 ========
1383 >>> from sympy.polys.domains import ZZ
1384 >>> from sympy.polys.galoistools import gf_irreducible
1385 >>> gf_irreducible(10, 5, ZZ) #doctest: +SKIP
1386 [1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4]
1388 """
1389 while True:
1390 f = gf_random(n, p, K)
1391 if gf_irreducible_p(f, p, K):
1392 return f
1395def gf_irred_p_ben_or(f, p, K):
1396 """
1397 Ben-Or's polynomial irreducibility test over finite fields.
1399 Examples
1400 ========
1402 >>> from sympy.polys.domains import ZZ
1403 >>> from sympy.polys.galoistools import gf_irred_p_ben_or
1405 >>> gf_irred_p_ben_or(ZZ.map([1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4]), 5, ZZ)
1406 True
1407 >>> gf_irred_p_ben_or(ZZ.map([3, 2, 4]), 5, ZZ)
1408 False
1410 """
1411 n = gf_degree(f)
1413 if n <= 1:
1414 return True
1416 _, f = gf_monic(f, p, K)
1417 if n < 5:
1418 H = h = gf_pow_mod([K.one, K.zero], p, f, p, K)
1420 for i in range(0, n//2):
1421 g = gf_sub(h, [K.one, K.zero], p, K)
1423 if gf_gcd(f, g, p, K) == [K.one]:
1424 h = gf_compose_mod(h, H, f, p, K)
1425 else:
1426 return False
1427 else:
1428 b = gf_frobenius_monomial_base(f, p, K)
1429 H = h = gf_frobenius_map([K.one, K.zero], f, b, p, K)
1430 for i in range(0, n//2):
1431 g = gf_sub(h, [K.one, K.zero], p, K)
1432 if gf_gcd(f, g, p, K) == [K.one]:
1433 h = gf_frobenius_map(h, f, b, p, K)
1434 else:
1435 return False
1437 return True
1440def gf_irred_p_rabin(f, p, K):
1441 """
1442 Rabin's polynomial irreducibility test over finite fields.
1444 Examples
1445 ========
1447 >>> from sympy.polys.domains import ZZ
1448 >>> from sympy.polys.galoistools import gf_irred_p_rabin
1450 >>> gf_irred_p_rabin(ZZ.map([1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4]), 5, ZZ)
1451 True
1452 >>> gf_irred_p_rabin(ZZ.map([3, 2, 4]), 5, ZZ)
1453 False
1455 """
1456 n = gf_degree(f)
1458 if n <= 1:
1459 return True
1461 _, f = gf_monic(f, p, K)
1463 x = [K.one, K.zero]
1465 from sympy.ntheory import factorint
1467 indices = { n//d for d in factorint(n) }
1469 b = gf_frobenius_monomial_base(f, p, K)
1470 h = b[1]
1472 for i in range(1, n):
1473 if i in indices:
1474 g = gf_sub(h, x, p, K)
1476 if gf_gcd(f, g, p, K) != [K.one]:
1477 return False
1479 h = gf_frobenius_map(h, f, b, p, K)
1481 return h == x
1483_irred_methods = {
1484 'ben-or': gf_irred_p_ben_or,
1485 'rabin': gf_irred_p_rabin,
1486}
1489def gf_irreducible_p(f, p, K):
1490 """
1491 Test irreducibility of a polynomial ``f`` in ``GF(p)[x]``.
1493 Examples
1494 ========
1496 >>> from sympy.polys.domains import ZZ
1497 >>> from sympy.polys.galoistools import gf_irreducible_p
1499 >>> gf_irreducible_p(ZZ.map([1, 4, 2, 2, 3, 2, 4, 1, 4, 0, 4]), 5, ZZ)
1500 True
1501 >>> gf_irreducible_p(ZZ.map([3, 2, 4]), 5, ZZ)
1502 False
1504 """
1505 method = query('GF_IRRED_METHOD')
1507 if method is not None:
1508 irred = _irred_methods[method](f, p, K)
1509 else:
1510 irred = gf_irred_p_rabin(f, p, K)
1512 return irred
1515def gf_sqf_p(f, p, K):
1516 """
1517 Return ``True`` if ``f`` is square-free in ``GF(p)[x]``.
1519 Examples
1520 ========
1522 >>> from sympy.polys.domains import ZZ
1523 >>> from sympy.polys.galoistools import gf_sqf_p
1525 >>> gf_sqf_p(ZZ.map([3, 2, 4]), 5, ZZ)
1526 True
1527 >>> gf_sqf_p(ZZ.map([2, 4, 4, 2, 2, 1, 4]), 5, ZZ)
1528 False
1530 """
1531 _, f = gf_monic(f, p, K)
1533 if not f:
1534 return True
1535 else:
1536 return gf_gcd(f, gf_diff(f, p, K), p, K) == [K.one]
1539def gf_sqf_part(f, p, K):
1540 """
1541 Return square-free part of a ``GF(p)[x]`` polynomial.
1543 Examples
1544 ========
1546 >>> from sympy.polys.domains import ZZ
1547 >>> from sympy.polys.galoistools import gf_sqf_part
1549 >>> gf_sqf_part(ZZ.map([1, 1, 3, 0, 1, 0, 2, 2, 1]), 5, ZZ)
1550 [1, 4, 3]
1552 """
1553 _, sqf = gf_sqf_list(f, p, K)
1555 g = [K.one]
1557 for f, _ in sqf:
1558 g = gf_mul(g, f, p, K)
1560 return g
1563def gf_sqf_list(f, p, K, all=False):
1564 """
1565 Return the square-free decomposition of a ``GF(p)[x]`` polynomial.
1567 Given a polynomial ``f`` in ``GF(p)[x]``, returns the leading coefficient
1568 of ``f`` and a square-free decomposition ``f_1**e_1 f_2**e_2 ... f_k**e_k``
1569 such that all ``f_i`` are monic polynomials and ``(f_i, f_j)`` for ``i != j``
1570 are co-prime and ``e_1 ... e_k`` are given in increasing order. All trivial
1571 terms (i.e. ``f_i = 1``) are not included in the output.
1573 Consider polynomial ``f = x**11 + 1`` over ``GF(11)[x]``::
1575 >>> from sympy.polys.domains import ZZ
1577 >>> from sympy.polys.galoistools import (
1578 ... gf_from_dict, gf_diff, gf_sqf_list, gf_pow,
1579 ... )
1580 ... # doctest: +NORMALIZE_WHITESPACE
1582 >>> f = gf_from_dict({11: ZZ(1), 0: ZZ(1)}, 11, ZZ)
1584 Note that ``f'(x) = 0``::
1586 >>> gf_diff(f, 11, ZZ)
1587 []
1589 This phenomenon does not happen in characteristic zero. However we can
1590 still compute square-free decomposition of ``f`` using ``gf_sqf()``::
1592 >>> gf_sqf_list(f, 11, ZZ)
1593 (1, [([1, 1], 11)])
1595 We obtained factorization ``f = (x + 1)**11``. This is correct because::
1597 >>> gf_pow([1, 1], 11, 11, ZZ) == f
1598 True
1600 References
1601 ==========
1603 .. [1] [Geddes92]_
1605 """
1606 n, sqf, factors, r = 1, False, [], int(p)
1608 lc, f = gf_monic(f, p, K)
1610 if gf_degree(f) < 1:
1611 return lc, []
1613 while True:
1614 F = gf_diff(f, p, K)
1616 if F != []:
1617 g = gf_gcd(f, F, p, K)
1618 h = gf_quo(f, g, p, K)
1620 i = 1
1622 while h != [K.one]:
1623 G = gf_gcd(g, h, p, K)
1624 H = gf_quo(h, G, p, K)
1626 if gf_degree(H) > 0:
1627 factors.append((H, i*n))
1629 g, h, i = gf_quo(g, G, p, K), G, i + 1
1631 if g == [K.one]:
1632 sqf = True
1633 else:
1634 f = g
1636 if not sqf:
1637 d = gf_degree(f) // r
1639 for i in range(0, d + 1):
1640 f[i] = f[i*r]
1642 f, n = f[:d + 1], n*r
1643 else:
1644 break
1646 if all:
1647 raise ValueError("'all=True' is not supported yet")
1649 return lc, factors
1652def gf_Qmatrix(f, p, K):
1653 """
1654 Calculate Berlekamp's ``Q`` matrix.
1656 Examples
1657 ========
1659 >>> from sympy.polys.domains import ZZ
1660 >>> from sympy.polys.galoistools import gf_Qmatrix
1662 >>> gf_Qmatrix([3, 2, 4], 5, ZZ)
1663 [[1, 0],
1664 [3, 4]]
1666 >>> gf_Qmatrix([1, 0, 0, 0, 1], 5, ZZ)
1667 [[1, 0, 0, 0],
1668 [0, 4, 0, 0],
1669 [0, 0, 1, 0],
1670 [0, 0, 0, 4]]
1672 """
1673 n, r = gf_degree(f), int(p)
1675 q = [K.one] + [K.zero]*(n - 1)
1676 Q = [list(q)] + [[]]*(n - 1)
1678 for i in range(1, (n - 1)*r + 1):
1679 qq, c = [(-q[-1]*f[-1]) % p], q[-1]
1681 for j in range(1, n):
1682 qq.append((q[j - 1] - c*f[-j - 1]) % p)
1684 if not (i % r):
1685 Q[i//r] = list(qq)
1687 q = qq
1689 return Q
1692def gf_Qbasis(Q, p, K):
1693 """
1694 Compute a basis of the kernel of ``Q``.
1696 Examples
1697 ========
1699 >>> from sympy.polys.domains import ZZ
1700 >>> from sympy.polys.galoistools import gf_Qmatrix, gf_Qbasis
1702 >>> gf_Qbasis(gf_Qmatrix([1, 0, 0, 0, 1], 5, ZZ), 5, ZZ)
1703 [[1, 0, 0, 0], [0, 0, 1, 0]]
1705 >>> gf_Qbasis(gf_Qmatrix([3, 2, 4], 5, ZZ), 5, ZZ)
1706 [[1, 0]]
1708 """
1709 Q, n = [ list(q) for q in Q ], len(Q)
1711 for k in range(0, n):
1712 Q[k][k] = (Q[k][k] - K.one) % p
1714 for k in range(0, n):
1715 for i in range(k, n):
1716 if Q[k][i]:
1717 break
1718 else:
1719 continue
1721 inv = K.invert(Q[k][i], p)
1723 for j in range(0, n):
1724 Q[j][i] = (Q[j][i]*inv) % p
1726 for j in range(0, n):
1727 t = Q[j][k]
1728 Q[j][k] = Q[j][i]
1729 Q[j][i] = t
1731 for i in range(0, n):
1732 if i != k:
1733 q = Q[k][i]
1735 for j in range(0, n):
1736 Q[j][i] = (Q[j][i] - Q[j][k]*q) % p
1738 for i in range(0, n):
1739 for j in range(0, n):
1740 if i == j:
1741 Q[i][j] = (K.one - Q[i][j]) % p
1742 else:
1743 Q[i][j] = (-Q[i][j]) % p
1745 basis = []
1747 for q in Q:
1748 if any(q):
1749 basis.append(q)
1751 return basis
1754def gf_berlekamp(f, p, K):
1755 """
1756 Factor a square-free ``f`` in ``GF(p)[x]`` for small ``p``.
1758 Examples
1759 ========
1761 >>> from sympy.polys.domains import ZZ
1762 >>> from sympy.polys.galoistools import gf_berlekamp
1764 >>> gf_berlekamp([1, 0, 0, 0, 1], 5, ZZ)
1765 [[1, 0, 2], [1, 0, 3]]
1767 """
1768 Q = gf_Qmatrix(f, p, K)
1769 V = gf_Qbasis(Q, p, K)
1771 for i, v in enumerate(V):
1772 V[i] = gf_strip(list(reversed(v)))
1774 factors = [f]
1776 for k in range(1, len(V)):
1777 for f in list(factors):
1778 s = K.zero
1780 while s < p:
1781 g = gf_sub_ground(V[k], s, p, K)
1782 h = gf_gcd(f, g, p, K)
1784 if h != [K.one] and h != f:
1785 factors.remove(f)
1787 f = gf_quo(f, h, p, K)
1788 factors.extend([f, h])
1790 if len(factors) == len(V):
1791 return _sort_factors(factors, multiple=False)
1793 s += K.one
1795 return _sort_factors(factors, multiple=False)
1798def gf_ddf_zassenhaus(f, p, K):
1799 """
1800 Cantor-Zassenhaus: Deterministic Distinct Degree Factorization
1802 Given a monic square-free polynomial ``f`` in ``GF(p)[x]``, computes
1803 partial distinct degree factorization ``f_1 ... f_d`` of ``f`` where
1804 ``deg(f_i) != deg(f_j)`` for ``i != j``. The result is returned as a
1805 list of pairs ``(f_i, e_i)`` where ``deg(f_i) > 0`` and ``e_i > 0``
1806 is an argument to the equal degree factorization routine.
1808 Consider the polynomial ``x**15 - 1`` in ``GF(11)[x]``::
1810 >>> from sympy.polys.domains import ZZ
1811 >>> from sympy.polys.galoistools import gf_from_dict
1813 >>> f = gf_from_dict({15: ZZ(1), 0: ZZ(-1)}, 11, ZZ)
1815 Distinct degree factorization gives::
1817 >>> from sympy.polys.galoistools import gf_ddf_zassenhaus
1819 >>> gf_ddf_zassenhaus(f, 11, ZZ)
1820 [([1, 0, 0, 0, 0, 10], 1), ([1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1], 2)]
1822 which means ``x**15 - 1 = (x**5 - 1) (x**10 + x**5 + 1)``. To obtain
1823 factorization into irreducibles, use equal degree factorization
1824 procedure (EDF) with each of the factors.
1826 References
1827 ==========
1829 .. [1] [Gathen99]_
1830 .. [2] [Geddes92]_
1832 """
1833 i, g, factors = 1, [K.one, K.zero], []
1835 b = gf_frobenius_monomial_base(f, p, K)
1836 while 2*i <= gf_degree(f):
1837 g = gf_frobenius_map(g, f, b, p, K)
1838 h = gf_gcd(f, gf_sub(g, [K.one, K.zero], p, K), p, K)
1840 if h != [K.one]:
1841 factors.append((h, i))
1843 f = gf_quo(f, h, p, K)
1844 g = gf_rem(g, f, p, K)
1845 b = gf_frobenius_monomial_base(f, p, K)
1847 i += 1
1849 if f != [K.one]:
1850 return factors + [(f, gf_degree(f))]
1851 else:
1852 return factors
1855def gf_edf_zassenhaus(f, n, p, K):
1856 """
1857 Cantor-Zassenhaus: Probabilistic Equal Degree Factorization
1859 Given a monic square-free polynomial ``f`` in ``GF(p)[x]`` and
1860 an integer ``n``, such that ``n`` divides ``deg(f)``, returns all
1861 irreducible factors ``f_1,...,f_d`` of ``f``, each of degree ``n``.
1862 EDF procedure gives complete factorization over Galois fields.
1864 Consider the square-free polynomial ``f = x**3 + x**2 + x + 1`` in
1865 ``GF(5)[x]``. Let's compute its irreducible factors of degree one::
1867 >>> from sympy.polys.domains import ZZ
1868 >>> from sympy.polys.galoistools import gf_edf_zassenhaus
1870 >>> gf_edf_zassenhaus([1,1,1,1], 1, 5, ZZ)
1871 [[1, 1], [1, 2], [1, 3]]
1873 Notes
1874 =====
1876 The case p == 2 is handled by Cohen's Algorithm 3.4.8. The case p odd is
1877 as in Geddes Algorithm 8.9 (or Cohen's Algorithm 3.4.6).
1879 References
1880 ==========
1882 .. [1] [Gathen99]_
1883 .. [2] [Geddes92]_ Algorithm 8.9
1884 .. [3] [Cohen93]_ Algorithm 3.4.8
1886 """
1887 factors = [f]
1889 if gf_degree(f) <= n:
1890 return factors
1892 N = gf_degree(f) // n
1893 if p != 2:
1894 b = gf_frobenius_monomial_base(f, p, K)
1896 t = [K.one, K.zero]
1897 while len(factors) < N:
1898 if p == 2:
1899 h = r = t
1901 for i in range(n - 1):
1902 r = gf_pow_mod(r, 2, f, p, K)
1903 h = gf_add(h, r, p, K)
1905 g = gf_gcd(f, h, p, K)
1906 t += [K.zero, K.zero]
1907 else:
1908 r = gf_random(2 * n - 1, p, K)
1909 h = _gf_pow_pnm1d2(r, n, f, b, p, K)
1910 g = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)
1912 if g != [K.one] and g != f:
1913 factors = gf_edf_zassenhaus(g, n, p, K) \
1914 + gf_edf_zassenhaus(gf_quo(f, g, p, K), n, p, K)
1916 return _sort_factors(factors, multiple=False)
1919def gf_ddf_shoup(f, p, K):
1920 """
1921 Kaltofen-Shoup: Deterministic Distinct Degree Factorization
1923 Given a monic square-free polynomial ``f`` in ``GF(p)[x]``, computes
1924 partial distinct degree factorization ``f_1,...,f_d`` of ``f`` where
1925 ``deg(f_i) != deg(f_j)`` for ``i != j``. The result is returned as a
1926 list of pairs ``(f_i, e_i)`` where ``deg(f_i) > 0`` and ``e_i > 0``
1927 is an argument to the equal degree factorization routine.
1929 This algorithm is an improved version of Zassenhaus algorithm for
1930 large ``deg(f)`` and modulus ``p`` (especially for ``deg(f) ~ lg(p)``).
1932 Examples
1933 ========
1935 >>> from sympy.polys.domains import ZZ
1936 >>> from sympy.polys.galoistools import gf_ddf_shoup, gf_from_dict
1938 >>> f = gf_from_dict({6: ZZ(1), 5: ZZ(-1), 4: ZZ(1), 3: ZZ(1), 1: ZZ(-1)}, 3, ZZ)
1940 >>> gf_ddf_shoup(f, 3, ZZ)
1941 [([1, 1, 0], 1), ([1, 1, 0, 1, 2], 2)]
1943 References
1944 ==========
1946 .. [1] [Kaltofen98]_
1947 .. [2] [Shoup95]_
1948 .. [3] [Gathen92]_
1950 """
1951 n = gf_degree(f)
1952 k = int(_ceil(_sqrt(n//2)))
1953 b = gf_frobenius_monomial_base(f, p, K)
1954 h = gf_frobenius_map([K.one, K.zero], f, b, p, K)
1955 # U[i] = x**(p**i)
1956 U = [[K.one, K.zero], h] + [K.zero]*(k - 1)
1958 for i in range(2, k + 1):
1959 U[i] = gf_frobenius_map(U[i-1], f, b, p, K)
1961 h, U = U[k], U[:k]
1962 # V[i] = x**(p**(k*(i+1)))
1963 V = [h] + [K.zero]*(k - 1)
1965 for i in range(1, k):
1966 V[i] = gf_compose_mod(V[i - 1], h, f, p, K)
1968 factors = []
1970 for i, v in enumerate(V):
1971 h, j = [K.one], k - 1
1973 for u in U:
1974 g = gf_sub(v, u, p, K)
1975 h = gf_mul(h, g, p, K)
1976 h = gf_rem(h, f, p, K)
1978 g = gf_gcd(f, h, p, K)
1979 f = gf_quo(f, g, p, K)
1981 for u in reversed(U):
1982 h = gf_sub(v, u, p, K)
1983 F = gf_gcd(g, h, p, K)
1985 if F != [K.one]:
1986 factors.append((F, k*(i + 1) - j))
1988 g, j = gf_quo(g, F, p, K), j - 1
1990 if f != [K.one]:
1991 factors.append((f, gf_degree(f)))
1993 return factors
1995def gf_edf_shoup(f, n, p, K):
1996 """
1997 Gathen-Shoup: Probabilistic Equal Degree Factorization
1999 Given a monic square-free polynomial ``f`` in ``GF(p)[x]`` and integer
2000 ``n`` such that ``n`` divides ``deg(f)``, returns all irreducible factors
2001 ``f_1,...,f_d`` of ``f``, each of degree ``n``. This is a complete
2002 factorization over Galois fields.
2004 This algorithm is an improved version of Zassenhaus algorithm for
2005 large ``deg(f)`` and modulus ``p`` (especially for ``deg(f) ~ lg(p)``).
2007 Examples
2008 ========
2010 >>> from sympy.polys.domains import ZZ
2011 >>> from sympy.polys.galoistools import gf_edf_shoup
2013 >>> gf_edf_shoup(ZZ.map([1, 2837, 2277]), 1, 2917, ZZ)
2014 [[1, 852], [1, 1985]]
2016 References
2017 ==========
2019 .. [1] [Shoup91]_
2020 .. [2] [Gathen92]_
2022 """
2023 N, q = gf_degree(f), int(p)
2025 if not N:
2026 return []
2027 if N <= n:
2028 return [f]
2030 factors, x = [f], [K.one, K.zero]
2032 r = gf_random(N - 1, p, K)
2034 if p == 2:
2035 h = gf_pow_mod(x, q, f, p, K)
2036 H = gf_trace_map(r, h, x, n - 1, f, p, K)[1]
2037 h1 = gf_gcd(f, H, p, K)
2038 h2 = gf_quo(f, h1, p, K)
2040 factors = gf_edf_shoup(h1, n, p, K) \
2041 + gf_edf_shoup(h2, n, p, K)
2042 else:
2043 b = gf_frobenius_monomial_base(f, p, K)
2044 H = _gf_trace_map(r, n, f, b, p, K)
2045 h = gf_pow_mod(H, (q - 1)//2, f, p, K)
2047 h1 = gf_gcd(f, h, p, K)
2048 h2 = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)
2049 h3 = gf_quo(f, gf_mul(h1, h2, p, K), p, K)
2051 factors = gf_edf_shoup(h1, n, p, K) \
2052 + gf_edf_shoup(h2, n, p, K) \
2053 + gf_edf_shoup(h3, n, p, K)
2055 return _sort_factors(factors, multiple=False)
2058def gf_zassenhaus(f, p, K):
2059 """
2060 Factor a square-free ``f`` in ``GF(p)[x]`` for medium ``p``.
2062 Examples
2063 ========
2065 >>> from sympy.polys.domains import ZZ
2066 >>> from sympy.polys.galoistools import gf_zassenhaus
2068 >>> gf_zassenhaus(ZZ.map([1, 4, 3]), 5, ZZ)
2069 [[1, 1], [1, 3]]
2071 """
2072 factors = []
2074 for factor, n in gf_ddf_zassenhaus(f, p, K):
2075 factors += gf_edf_zassenhaus(factor, n, p, K)
2077 return _sort_factors(factors, multiple=False)
2080def gf_shoup(f, p, K):
2081 """
2082 Factor a square-free ``f`` in ``GF(p)[x]`` for large ``p``.
2084 Examples
2085 ========
2087 >>> from sympy.polys.domains import ZZ
2088 >>> from sympy.polys.galoistools import gf_shoup
2090 >>> gf_shoup(ZZ.map([1, 4, 3]), 5, ZZ)
2091 [[1, 1], [1, 3]]
2093 """
2094 factors = []
2096 for factor, n in gf_ddf_shoup(f, p, K):
2097 factors += gf_edf_shoup(factor, n, p, K)
2099 return _sort_factors(factors, multiple=False)
2101_factor_methods = {
2102 'berlekamp': gf_berlekamp, # ``p`` : small
2103 'zassenhaus': gf_zassenhaus, # ``p`` : medium
2104 'shoup': gf_shoup, # ``p`` : large
2105}
2108def gf_factor_sqf(f, p, K, method=None):
2109 """
2110 Factor a square-free polynomial ``f`` in ``GF(p)[x]``.
2112 Examples
2113 ========
2115 >>> from sympy.polys.domains import ZZ
2116 >>> from sympy.polys.galoistools import gf_factor_sqf
2118 >>> gf_factor_sqf(ZZ.map([3, 2, 4]), 5, ZZ)
2119 (3, [[1, 1], [1, 3]])
2121 """
2122 lc, f = gf_monic(f, p, K)
2124 if gf_degree(f) < 1:
2125 return lc, []
2127 method = method or query('GF_FACTOR_METHOD')
2129 if method is not None:
2130 factors = _factor_methods[method](f, p, K)
2131 else:
2132 factors = gf_zassenhaus(f, p, K)
2134 return lc, factors
2137def gf_factor(f, p, K):
2138 """
2139 Factor (non square-free) polynomials in ``GF(p)[x]``.
2141 Given a possibly non square-free polynomial ``f`` in ``GF(p)[x]``,
2142 returns its complete factorization into irreducibles::
2144 f_1(x)**e_1 f_2(x)**e_2 ... f_d(x)**e_d
2146 where each ``f_i`` is a monic polynomial and ``gcd(f_i, f_j) == 1``,
2147 for ``i != j``. The result is given as a tuple consisting of the
2148 leading coefficient of ``f`` and a list of factors of ``f`` with
2149 their multiplicities.
2151 The algorithm proceeds by first computing square-free decomposition
2152 of ``f`` and then iteratively factoring each of square-free factors.
2154 Consider a non square-free polynomial ``f = (7*x + 1) (x + 2)**2`` in
2155 ``GF(11)[x]``. We obtain its factorization into irreducibles as follows::
2157 >>> from sympy.polys.domains import ZZ
2158 >>> from sympy.polys.galoistools import gf_factor
2160 >>> gf_factor(ZZ.map([5, 2, 7, 2]), 11, ZZ)
2161 (5, [([1, 2], 1), ([1, 8], 2)])
2163 We arrived with factorization ``f = 5 (x + 2) (x + 8)**2``. We did not
2164 recover the exact form of the input polynomial because we requested to
2165 get monic factors of ``f`` and its leading coefficient separately.
2167 Square-free factors of ``f`` can be factored into irreducibles over
2168 ``GF(p)`` using three very different methods:
2170 Berlekamp
2171 efficient for very small values of ``p`` (usually ``p < 25``)
2172 Cantor-Zassenhaus
2173 efficient on average input and with "typical" ``p``
2174 Shoup-Kaltofen-Gathen
2175 efficient with very large inputs and modulus
2177 If you want to use a specific factorization method, instead of the default
2178 one, set ``GF_FACTOR_METHOD`` with one of ``berlekamp``, ``zassenhaus`` or
2179 ``shoup`` values.
2181 References
2182 ==========
2184 .. [1] [Gathen99]_
2186 """
2187 lc, f = gf_monic(f, p, K)
2189 if gf_degree(f) < 1:
2190 return lc, []
2192 factors = []
2194 for g, n in gf_sqf_list(f, p, K)[1]:
2195 for h in gf_factor_sqf(g, p, K)[1]:
2196 factors.append((h, n))
2198 return lc, _sort_factors(factors)
2201def gf_value(f, a):
2202 """
2203 Value of polynomial 'f' at 'a' in field R.
2205 Examples
2206 ========
2208 >>> from sympy.polys.galoistools import gf_value
2210 >>> gf_value([1, 7, 2, 4], 11)
2211 2204
2213 """
2214 result = 0
2215 for c in f:
2216 result *= a
2217 result += c
2218 return result
2221def linear_congruence(a, b, m):
2222 """
2223 Returns the values of x satisfying a*x congruent b mod(m)
2225 Here m is positive integer and a, b are natural numbers.
2226 This function returns only those values of x which are distinct mod(m).
2228 Examples
2229 ========
2231 >>> from sympy.polys.galoistools import linear_congruence
2233 >>> linear_congruence(3, 12, 15)
2234 [4, 9, 14]
2236 There are 3 solutions distinct mod(15) since gcd(a, m) = gcd(3, 15) = 3.
2238 References
2239 ==========
2241 .. [1] https://en.wikipedia.org/wiki/Linear_congruence_theorem
2243 """
2244 from sympy.polys.polytools import gcdex
2245 if a % m == 0:
2246 if b % m == 0:
2247 return list(range(m))
2248 else:
2249 return []
2250 r, _, g = gcdex(a, m)
2251 if b % g != 0:
2252 return []
2253 return [(r * b // g + t * m // g) % m for t in range(g)]
2256def _raise_mod_power(x, s, p, f):
2257 """
2258 Used in gf_csolve to generate solutions of f(x) cong 0 mod(p**(s + 1))
2259 from the solutions of f(x) cong 0 mod(p**s).
2261 Examples
2262 ========
2264 >>> from sympy.polys.galoistools import _raise_mod_power
2265 >>> from sympy.polys.galoistools import csolve_prime
2267 These is the solutions of f(x) = x**2 + x + 7 cong 0 mod(3)
2269 >>> f = [1, 1, 7]
2270 >>> csolve_prime(f, 3)
2271 [1]
2272 >>> [ i for i in range(3) if not (i**2 + i + 7) % 3]
2273 [1]
2275 The solutions of f(x) cong 0 mod(9) are constructed from the
2276 values returned from _raise_mod_power:
2278 >>> x, s, p = 1, 1, 3
2279 >>> V = _raise_mod_power(x, s, p, f)
2280 >>> [x + v * p**s for v in V]
2281 [1, 4, 7]
2283 And these are confirmed with the following:
2285 >>> [ i for i in range(3**2) if not (i**2 + i + 7) % 3**2]
2286 [1, 4, 7]
2288 """
2289 from sympy.polys.domains import ZZ
2290 f_f = gf_diff(f, p, ZZ)
2291 alpha = gf_value(f_f, x)
2292 beta = - gf_value(f, x) // p**s
2293 return linear_congruence(alpha, beta, p)
2296def csolve_prime(f, p, e=1):
2297 """
2298 Solutions of f(x) congruent 0 mod(p**e).
2300 Examples
2301 ========
2303 >>> from sympy.polys.galoistools import csolve_prime
2305 >>> csolve_prime([1, 1, 7], 3, 1)
2306 [1]
2307 >>> csolve_prime([1, 1, 7], 3, 2)
2308 [1, 4, 7]
2310 Solutions [7, 4, 1] (mod 3**2) are generated by ``_raise_mod_power()``
2311 from solution [1] (mod 3).
2312 """
2313 from sympy.polys.domains import ZZ
2314 X1 = [i for i in range(p) if gf_eval(f, i, p, ZZ) == 0]
2315 if e == 1:
2316 return X1
2317 X = []
2318 S = list(zip(X1, [1]*len(X1)))
2319 while S:
2320 x, s = S.pop()
2321 if s == e:
2322 X.append(x)
2323 else:
2324 s1 = s + 1
2325 ps = p**s
2326 S.extend([(x + v*ps, s1) for v in _raise_mod_power(x, s, p, f)])
2327 return sorted(X)
2330def gf_csolve(f, n):
2331 """
2332 To solve f(x) congruent 0 mod(n).
2334 n is divided into canonical factors and f(x) cong 0 mod(p**e) will be
2335 solved for each factor. Applying the Chinese Remainder Theorem to the
2336 results returns the final answers.
2338 Examples
2339 ========
2341 Solve [1, 1, 7] congruent 0 mod(189):
2343 >>> from sympy.polys.galoistools import gf_csolve
2344 >>> gf_csolve([1, 1, 7], 189)
2345 [13, 49, 76, 112, 139, 175]
2347 References
2348 ==========
2350 .. [1] 'An introduction to the Theory of Numbers' 5th Edition by Ivan Niven,
2351 Zuckerman and Montgomery.
2353 """
2354 from sympy.polys.domains import ZZ
2355 from sympy.ntheory import factorint
2356 P = factorint(n)
2357 X = [csolve_prime(f, p, e) for p, e in P.items()]
2358 pools = list(map(tuple, X))
2359 perms = [[]]
2360 for pool in pools:
2361 perms = [x + [y] for x in perms for y in pool]
2362 dist_factors = [pow(p, e) for p, e in P.items()]
2363 return sorted([gf_crt(per, dist_factors, ZZ) for per in perms])