Coverage for /usr/lib/python3/dist-packages/sympy/ntheory/residue_ntheory.py: 6%
786 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import annotations
3from sympy.core.function import Function
4from sympy.core.numbers import igcd, igcdex, mod_inverse
5from sympy.core.power import isqrt
6from sympy.core.singleton import S
7from sympy.polys import Poly
8from sympy.polys.domains import ZZ
9from sympy.polys.galoistools import gf_crt1, gf_crt2, linear_congruence
10from .primetest import isprime
11from .factor_ import factorint, trailing, totient, multiplicity, perfect_power
12from sympy.utilities.misc import as_int
13from sympy.core.random import _randint, randint
15from itertools import cycle, product
18def n_order(a, n):
19 """Returns the order of ``a`` modulo ``n``.
21 The order of ``a`` modulo ``n`` is the smallest integer
22 ``k`` such that ``a**k`` leaves a remainder of 1 with ``n``.
24 Parameters
25 ==========
27 a : integer
28 n : integer, n > 1. a and n should be relatively prime
30 Examples
31 ========
33 >>> from sympy.ntheory import n_order
34 >>> n_order(3, 7)
35 6
36 >>> n_order(4, 7)
37 3
38 """
39 from collections import defaultdict
40 a, n = as_int(a), as_int(n)
41 if n <= 1:
42 raise ValueError("n should be an integer greater than 1")
43 a = a % n
44 # Trivial
45 if a == 1:
46 return 1
47 if igcd(a, n) != 1:
48 raise ValueError("The two numbers should be relatively prime")
49 # We want to calculate
50 # order = totient(n), factors = factorint(order)
51 factors = defaultdict(int)
52 for px, kx in factorint(n).items():
53 if kx > 1:
54 factors[px] += kx - 1
55 for py, ky in factorint(px - 1).items():
56 factors[py] += ky
57 order = 1
58 for px, kx in factors.items():
59 order *= px**kx
60 # Now the `order` is the order of the group.
61 # The order of `a` divides the order of the group.
62 for p, e in factors.items():
63 for _ in range(e):
64 if pow(a, order // p, n) == 1:
65 order //= p
66 else:
67 break
68 return order
71def _primitive_root_prime_iter(p):
72 """
73 Generates the primitive roots for a prime ``p``
75 Examples
76 ========
78 >>> from sympy.ntheory.residue_ntheory import _primitive_root_prime_iter
79 >>> list(_primitive_root_prime_iter(19))
80 [2, 3, 10, 13, 14, 15]
82 References
83 ==========
85 .. [1] W. Stein "Elementary Number Theory" (2011), page 44
87 """
88 # it is assumed that p is an int
89 v = [(p - 1) // i for i in factorint(p - 1).keys()]
90 a = 2
91 while a < p:
92 for pw in v:
93 # a TypeError below may indicate that p was not an int
94 if pow(a, pw, p) == 1:
95 break
96 else:
97 yield a
98 a += 1
101def primitive_root(p):
102 """
103 Returns the smallest primitive root or None.
105 Parameters
106 ==========
108 p : positive integer
110 Examples
111 ========
113 >>> from sympy.ntheory.residue_ntheory import primitive_root
114 >>> primitive_root(19)
115 2
117 References
118 ==========
120 .. [1] W. Stein "Elementary Number Theory" (2011), page 44
121 .. [2] P. Hackman "Elementary Number Theory" (2009), Chapter C
123 """
124 p = as_int(p)
125 if p < 1:
126 raise ValueError('p is required to be positive')
127 if p <= 2:
128 return 1
129 f = factorint(p)
130 if len(f) > 2:
131 return None
132 if len(f) == 2:
133 if 2 not in f or f[2] > 1:
134 return None
136 # case p = 2*p1**k, p1 prime
137 for p1, e1 in f.items():
138 if p1 != 2:
139 break
140 i = 1
141 while i < p:
142 i += 2
143 if i % p1 == 0:
144 continue
145 if is_primitive_root(i, p):
146 return i
148 else:
149 if 2 in f:
150 if p == 4:
151 return 3
152 return None
153 p1, n = list(f.items())[0]
154 if n > 1:
155 # see Ref [2], page 81
156 g = primitive_root(p1)
157 if is_primitive_root(g, p1**2):
158 return g
159 else:
160 for i in range(2, g + p1 + 1):
161 if igcd(i, p) == 1 and is_primitive_root(i, p):
162 return i
164 return next(_primitive_root_prime_iter(p))
167def is_primitive_root(a, p):
168 """
169 Returns True if ``a`` is a primitive root of ``p``.
171 ``a`` is said to be the primitive root of ``p`` if gcd(a, p) == 1 and
172 totient(p) is the smallest positive number s.t.
174 a**totient(p) cong 1 mod(p)
176 Parameters
177 ==========
179 a : integer
180 p : integer, p > 1. a and p should be relatively prime
182 Examples
183 ========
185 >>> from sympy.ntheory import is_primitive_root, n_order, totient
186 >>> is_primitive_root(3, 10)
187 True
188 >>> is_primitive_root(9, 10)
189 False
190 >>> n_order(3, 10) == totient(10)
191 True
192 >>> n_order(9, 10) == totient(10)
193 False
195 """
196 a, p = as_int(a), as_int(p)
197 if p <= 1:
198 raise ValueError("p should be an integer greater than 1")
199 a = a % p
200 if igcd(a, p) != 1:
201 raise ValueError("The two numbers should be relatively prime")
202 # Primitive root of p exist only for
203 # p = 2, 4, q**e, 2*q**e (q is odd prime)
204 if p <= 4:
205 # The primitive root is only p-1.
206 return a == p - 1
207 t = trailing(p)
208 if t > 1:
209 return False
210 q = p >> t
211 if isprime(q):
212 group_order = q - 1
213 factors = set(factorint(q - 1).keys())
214 else:
215 m = perfect_power(q)
216 if not m:
217 return False
218 q, e = m
219 if not isprime(q):
220 return False
221 group_order = q**(e - 1)*(q - 1)
222 factors = set(factorint(q - 1).keys())
223 factors.add(q)
224 return all(pow(a, group_order // prime, p) != 1 for prime in factors)
227def _sqrt_mod_tonelli_shanks(a, p):
228 """
229 Returns the square root in the case of ``p`` prime with ``p == 1 (mod 8)``
231 References
232 ==========
234 .. [1] R. Crandall and C. Pomerance "Prime Numbers", 2nt Ed., page 101
236 """
237 s = trailing(p - 1)
238 t = p >> s
239 # find a non-quadratic residue
240 while 1:
241 d = randint(2, p - 1)
242 r = legendre_symbol(d, p)
243 if r == -1:
244 break
245 #assert legendre_symbol(d, p) == -1
246 A = pow(a, t, p)
247 D = pow(d, t, p)
248 m = 0
249 for i in range(s):
250 adm = A*pow(D, m, p) % p
251 adm = pow(adm, 2**(s - 1 - i), p)
252 if adm % p == p - 1:
253 m += 2**i
254 #assert A*pow(D, m, p) % p == 1
255 x = pow(a, (t + 1)//2, p)*pow(D, m//2, p) % p
256 return x
259def sqrt_mod(a, p, all_roots=False):
260 """
261 Find a root of ``x**2 = a mod p``.
263 Parameters
264 ==========
266 a : integer
267 p : positive integer
268 all_roots : if True the list of roots is returned or None
270 Notes
271 =====
273 If there is no root it is returned None; else the returned root
274 is less or equal to ``p // 2``; in general is not the smallest one.
275 It is returned ``p // 2`` only if it is the only root.
277 Use ``all_roots`` only when it is expected that all the roots fit
278 in memory; otherwise use ``sqrt_mod_iter``.
280 Examples
281 ========
283 >>> from sympy.ntheory import sqrt_mod
284 >>> sqrt_mod(11, 43)
285 21
286 >>> sqrt_mod(17, 32, True)
287 [7, 9, 23, 25]
288 """
289 if all_roots:
290 return sorted(sqrt_mod_iter(a, p))
291 try:
292 p = abs(as_int(p))
293 it = sqrt_mod_iter(a, p)
294 r = next(it)
295 if r > p // 2:
296 return p - r
297 elif r < p // 2:
298 return r
299 else:
300 try:
301 r = next(it)
302 if r > p // 2:
303 return p - r
304 except StopIteration:
305 pass
306 return r
307 except StopIteration:
308 return None
311def _product(*iters):
312 """
313 Cartesian product generator
315 Notes
316 =====
318 Unlike itertools.product, it works also with iterables which do not fit
319 in memory. See https://bugs.python.org/issue10109
321 Author: Fernando Sumudu
322 with small changes
323 """
324 inf_iters = tuple(cycle(enumerate(it)) for it in iters)
325 num_iters = len(inf_iters)
326 cur_val = [None]*num_iters
328 first_v = True
329 while True:
330 i, p = 0, num_iters
331 while p and not i:
332 p -= 1
333 i, cur_val[p] = next(inf_iters[p])
335 if not p and not i:
336 if first_v:
337 first_v = False
338 else:
339 break
341 yield cur_val
344def sqrt_mod_iter(a, p, domain=int):
345 """
346 Iterate over solutions to ``x**2 = a mod p``.
348 Parameters
349 ==========
351 a : integer
352 p : positive integer
353 domain : integer domain, ``int``, ``ZZ`` or ``Integer``
355 Examples
356 ========
358 >>> from sympy.ntheory.residue_ntheory import sqrt_mod_iter
359 >>> list(sqrt_mod_iter(11, 43))
360 [21, 22]
361 """
362 a, p = as_int(a), abs(as_int(p))
363 if isprime(p):
364 a = a % p
365 if a == 0:
366 res = _sqrt_mod1(a, p, 1)
367 else:
368 res = _sqrt_mod_prime_power(a, p, 1)
369 if res:
370 if domain is ZZ:
371 yield from res
372 else:
373 for x in res:
374 yield domain(x)
375 else:
376 f = factorint(p)
377 v = []
378 pv = []
379 for px, ex in f.items():
380 if a % px == 0:
381 rx = _sqrt_mod1(a, px, ex)
382 if not rx:
383 return
384 else:
385 rx = _sqrt_mod_prime_power(a, px, ex)
386 if not rx:
387 return
388 v.append(rx)
389 pv.append(px**ex)
390 mm, e, s = gf_crt1(pv, ZZ)
391 if domain is ZZ:
392 for vx in _product(*v):
393 r = gf_crt2(vx, pv, mm, e, s, ZZ)
394 yield r
395 else:
396 for vx in _product(*v):
397 r = gf_crt2(vx, pv, mm, e, s, ZZ)
398 yield domain(r)
401def _sqrt_mod_prime_power(a, p, k):
402 """
403 Find the solutions to ``x**2 = a mod p**k`` when ``a % p != 0``
405 Parameters
406 ==========
408 a : integer
409 p : prime number
410 k : positive integer
412 Examples
413 ========
415 >>> from sympy.ntheory.residue_ntheory import _sqrt_mod_prime_power
416 >>> _sqrt_mod_prime_power(11, 43, 1)
417 [21, 22]
419 References
420 ==========
422 .. [1] P. Hackman "Elementary Number Theory" (2009), page 160
423 .. [2] http://www.numbertheory.org/php/squareroot.html
424 .. [3] [Gathen99]_
425 """
426 pk = p**k
427 a = a % pk
429 if k == 1:
430 if p == 2:
431 return [ZZ(a)]
432 if not (a % p < 2 or pow(a, (p - 1) // 2, p) == 1):
433 return None
435 if p % 4 == 3:
436 res = pow(a, (p + 1) // 4, p)
437 elif p % 8 == 5:
438 sign = pow(a, (p - 1) // 4, p)
439 if sign == 1:
440 res = pow(a, (p + 3) // 8, p)
441 else:
442 b = pow(4*a, (p - 5) // 8, p)
443 x = (2*a*b) % p
444 if pow(x, 2, p) == a:
445 res = x
446 else:
447 res = _sqrt_mod_tonelli_shanks(a, p)
449 # ``_sqrt_mod_tonelli_shanks(a, p)`` is not deterministic;
450 # sort to get always the same result
451 return sorted([ZZ(res), ZZ(p - res)])
453 if k > 1:
454 # see Ref.[2]
455 if p == 2:
456 if a % 8 != 1:
457 return None
458 if k <= 3:
459 s = set()
460 for i in range(0, pk, 4):
461 s.add(1 + i)
462 s.add(-1 + i)
463 return list(s)
464 # according to Ref.[2] for k > 2 there are two solutions
465 # (mod 2**k-1), that is four solutions (mod 2**k), which can be
466 # obtained from the roots of x**2 = 0 (mod 8)
467 rv = [ZZ(1), ZZ(3), ZZ(5), ZZ(7)]
468 # hensel lift them to solutions of x**2 = 0 (mod 2**k)
469 # if r**2 - a = 0 mod 2**nx but not mod 2**(nx+1)
470 # then r + 2**(nx - 1) is a root mod 2**(nx+1)
471 n = 3
472 res = []
473 for r in rv:
474 nx = n
475 while nx < k:
476 r1 = (r**2 - a) >> nx
477 if r1 % 2:
478 r = r + (1 << (nx - 1))
479 #assert (r**2 - a)% (1 << (nx + 1)) == 0
480 nx += 1
481 if r not in res:
482 res.append(r)
483 x = r + (1 << (k - 1))
484 #assert (x**2 - a) % pk == 0
485 if x < (1 << nx) and x not in res:
486 if (x**2 - a) % pk == 0:
487 res.append(x)
488 return res
489 rv = _sqrt_mod_prime_power(a, p, 1)
490 if not rv:
491 return None
492 r = rv[0]
493 fr = r**2 - a
494 # hensel lifting with Newton iteration, see Ref.[3] chapter 9
495 # with f(x) = x**2 - a; one has f'(a) != 0 (mod p) for p != 2
496 n = 1
497 px = p
498 while 1:
499 n1 = n
500 n1 *= 2
501 if n1 > k:
502 break
503 n = n1
504 px = px**2
505 frinv = igcdex(2*r, px)[0]
506 r = (r - fr*frinv) % px
507 fr = r**2 - a
508 if n < k:
509 px = p**k
510 frinv = igcdex(2*r, px)[0]
511 r = (r - fr*frinv) % px
512 return [r, px - r]
515def _sqrt_mod1(a, p, n):
516 """
517 Find solution to ``x**2 == a mod p**n`` when ``a % p == 0``
519 see http://www.numbertheory.org/php/squareroot.html
520 """
521 pn = p**n
522 a = a % pn
523 if a == 0:
524 # case gcd(a, p**k) = p**n
525 m = n // 2
526 if n % 2 == 1:
527 pm1 = p**(m + 1)
528 def _iter0a():
529 i = 0
530 while i < pn:
531 yield i
532 i += pm1
533 return _iter0a()
534 else:
535 pm = p**m
536 def _iter0b():
537 i = 0
538 while i < pn:
539 yield i
540 i += pm
541 return _iter0b()
543 # case gcd(a, p**k) = p**r, r < n
544 f = factorint(a)
545 r = f[p]
546 if r % 2 == 1:
547 return None
548 m = r // 2
549 a1 = a >> r
550 if p == 2:
551 if n - r == 1:
552 pnm1 = 1 << (n - m + 1)
553 pm1 = 1 << (m + 1)
554 def _iter1():
555 k = 1 << (m + 2)
556 i = 1 << m
557 while i < pnm1:
558 j = i
559 while j < pn:
560 yield j
561 j += k
562 i += pm1
563 return _iter1()
564 if n - r == 2:
565 res = _sqrt_mod_prime_power(a1, p, n - r)
566 if res is None:
567 return None
568 pnm = 1 << (n - m)
569 def _iter2():
570 s = set()
571 for r in res:
572 i = 0
573 while i < pn:
574 x = (r << m) + i
575 if x not in s:
576 s.add(x)
577 yield x
578 i += pnm
579 return _iter2()
580 if n - r > 2:
581 res = _sqrt_mod_prime_power(a1, p, n - r)
582 if res is None:
583 return None
584 pnm1 = 1 << (n - m - 1)
585 def _iter3():
586 s = set()
587 for r in res:
588 i = 0
589 while i < pn:
590 x = ((r << m) + i) % pn
591 if x not in s:
592 s.add(x)
593 yield x
594 i += pnm1
595 return _iter3()
596 else:
597 m = r // 2
598 a1 = a // p**r
599 res1 = _sqrt_mod_prime_power(a1, p, n - r)
600 if res1 is None:
601 return None
602 pm = p**m
603 pnr = p**(n-r)
604 pnm = p**(n-m)
606 def _iter4():
607 s = set()
608 pm = p**m
609 for rx in res1:
610 i = 0
611 while i < pnm:
612 x = ((rx + i) % pn)
613 if x not in s:
614 s.add(x)
615 yield x*pm
616 i += pnr
617 return _iter4()
620def is_quad_residue(a, p):
621 """
622 Returns True if ``a`` (mod ``p``) is in the set of squares mod ``p``,
623 i.e a % p in set([i**2 % p for i in range(p)]).
625 Examples
626 ========
628 If ``p`` is an odd
629 prime, an iterative method is used to make the determination:
631 >>> from sympy.ntheory import is_quad_residue
632 >>> sorted(set([i**2 % 7 for i in range(7)]))
633 [0, 1, 2, 4]
634 >>> [j for j in range(7) if is_quad_residue(j, 7)]
635 [0, 1, 2, 4]
637 See Also
638 ========
640 legendre_symbol, jacobi_symbol
641 """
642 a, p = as_int(a), as_int(p)
643 if p < 1:
644 raise ValueError('p must be > 0')
645 if a >= p or a < 0:
646 a = a % p
647 if a < 2 or p < 3:
648 return True
649 if not isprime(p):
650 if p % 2 and jacobi_symbol(a, p) == -1:
651 return False
652 r = sqrt_mod(a, p)
653 if r is None:
654 return False
655 else:
656 return True
658 return pow(a, (p - 1) // 2, p) == 1
661def is_nthpow_residue(a, n, m):
662 """
663 Returns True if ``x**n == a (mod m)`` has solutions.
665 References
666 ==========
668 .. [1] P. Hackman "Elementary Number Theory" (2009), page 76
670 """
671 a = a % m
672 a, n, m = as_int(a), as_int(n), as_int(m)
673 if m <= 0:
674 raise ValueError('m must be > 0')
675 if n < 0:
676 raise ValueError('n must be >= 0')
677 if n == 0:
678 if m == 1:
679 return False
680 return a == 1
681 if a == 0:
682 return True
683 if n == 1:
684 return True
685 if n == 2:
686 return is_quad_residue(a, m)
687 return _is_nthpow_residue_bign(a, n, m)
690def _is_nthpow_residue_bign(a, n, m):
691 r"""Returns True if `x^n = a \pmod{n}` has solutions for `n > 2`."""
692 # assert n > 2
693 # assert a > 0 and m > 0
694 if primitive_root(m) is None or igcd(a, m) != 1:
695 # assert m >= 8
696 for prime, power in factorint(m).items():
697 if not _is_nthpow_residue_bign_prime_power(a, n, prime, power):
698 return False
699 return True
700 f = totient(m)
701 k = int(f // igcd(f, n))
702 return pow(a, k, int(m)) == 1
705def _is_nthpow_residue_bign_prime_power(a, n, p, k):
706 r"""Returns True/False if a solution for `x^n = a \pmod{p^k}`
707 does/does not exist."""
708 # assert a > 0
709 # assert n > 2
710 # assert p is prime
711 # assert k > 0
712 if a % p:
713 if p != 2:
714 return _is_nthpow_residue_bign(a, n, pow(p, k))
715 if n & 1:
716 return True
717 c = trailing(n)
718 return a % pow(2, min(c + 2, k)) == 1
719 else:
720 a %= pow(p, k)
721 if not a:
722 return True
723 mu = multiplicity(p, a)
724 if mu % n:
725 return False
726 pm = pow(p, mu)
727 return _is_nthpow_residue_bign_prime_power(a//pm, n, p, k - mu)
730def _nthroot_mod2(s, q, p):
731 f = factorint(q)
732 v = []
733 for b, e in f.items():
734 v.extend([b]*e)
735 for qx in v:
736 s = _nthroot_mod1(s, qx, p, False)
737 return s
740def _nthroot_mod1(s, q, p, all_roots):
741 """
742 Root of ``x**q = s mod p``, ``p`` prime and ``q`` divides ``p - 1``
744 References
745 ==========
747 .. [1] A. M. Johnston "A Generalized qth Root Algorithm"
749 """
750 g = primitive_root(p)
751 if not isprime(q):
752 r = _nthroot_mod2(s, q, p)
753 else:
754 f = p - 1
755 assert (p - 1) % q == 0
756 # determine k
757 k = 0
758 while f % q == 0:
759 k += 1
760 f = f // q
761 # find z, x, r1
762 f1 = igcdex(-f, q)[0] % q
763 z = f*f1
764 x = (1 + z) // q
765 r1 = pow(s, x, p)
766 s1 = pow(s, f, p)
767 h = pow(g, f*q, p)
768 t = discrete_log(p, s1, h)
769 g2 = pow(g, z*t, p)
770 g3 = igcdex(g2, p)[0]
771 r = r1*g3 % p
772 #assert pow(r, q, p) == s
773 res = [r]
774 h = pow(g, (p - 1) // q, p)
775 #assert pow(h, q, p) == 1
776 hx = r
777 for i in range(q - 1):
778 hx = (hx*h) % p
779 res.append(hx)
780 if all_roots:
781 res.sort()
782 return res
783 return min(res)
787def _help(m, prime_modulo_method, diff_method, expr_val):
788 """
789 Helper function for _nthroot_mod_composite and polynomial_congruence.
791 Parameters
792 ==========
794 m : positive integer
795 prime_modulo_method : function to calculate the root of the congruence
796 equation for the prime divisors of m
797 diff_method : function to calculate derivative of expression at any
798 given point
799 expr_val : function to calculate value of the expression at any
800 given point
801 """
802 from sympy.ntheory.modular import crt
803 f = factorint(m)
804 dd = {}
805 for p, e in f.items():
806 tot_roots = set()
807 if e == 1:
808 tot_roots.update(prime_modulo_method(p))
809 else:
810 for root in prime_modulo_method(p):
811 diff = diff_method(root, p)
812 if diff != 0:
813 ppow = p
814 m_inv = mod_inverse(diff, p)
815 for j in range(1, e):
816 ppow *= p
817 root = (root - expr_val(root, ppow) * m_inv) % ppow
818 tot_roots.add(root)
819 else:
820 new_base = p
821 roots_in_base = {root}
822 while new_base < pow(p, e):
823 new_base *= p
824 new_roots = set()
825 for k in roots_in_base:
826 if expr_val(k, new_base)!= 0:
827 continue
828 while k not in new_roots:
829 new_roots.add(k)
830 k = (k + (new_base // p)) % new_base
831 roots_in_base = new_roots
832 tot_roots = tot_roots | roots_in_base
833 if tot_roots == set():
834 return []
835 dd[pow(p, e)] = tot_roots
836 a = []
837 m = []
838 for x, y in dd.items():
839 m.append(x)
840 a.append(list(y))
841 return sorted({crt(m, list(i))[0] for i in product(*a)})
844def _nthroot_mod_composite(a, n, m):
845 """
846 Find the solutions to ``x**n = a mod m`` when m is not prime.
847 """
848 return _help(m,
849 lambda p: nthroot_mod(a, n, p, True),
850 lambda root, p: (pow(root, n - 1, p) * (n % p)) % p,
851 lambda root, p: (pow(root, n, p) - a) % p)
854def nthroot_mod(a, n, p, all_roots=False):
855 """
856 Find the solutions to ``x**n = a mod p``.
858 Parameters
859 ==========
861 a : integer
862 n : positive integer
863 p : positive integer
864 all_roots : if False returns the smallest root, else the list of roots
866 Examples
867 ========
869 >>> from sympy.ntheory.residue_ntheory import nthroot_mod
870 >>> nthroot_mod(11, 4, 19)
871 8
872 >>> nthroot_mod(11, 4, 19, True)
873 [8, 11]
874 >>> nthroot_mod(68, 3, 109)
875 23
876 """
877 a = a % p
878 a, n, p = as_int(a), as_int(n), as_int(p)
880 if n == 2:
881 return sqrt_mod(a, p, all_roots)
882 # see Hackman "Elementary Number Theory" (2009), page 76
883 if not isprime(p):
884 return _nthroot_mod_composite(a, n, p)
885 if a % p == 0:
886 return [0]
887 if not is_nthpow_residue(a, n, p):
888 return [] if all_roots else None
889 if (p - 1) % n == 0:
890 return _nthroot_mod1(a, n, p, all_roots)
891 # The roots of ``x**n - a = 0 (mod p)`` are roots of
892 # ``gcd(x**n - a, x**(p - 1) - 1) = 0 (mod p)``
893 pa = n
894 pb = p - 1
895 b = 1
896 if pa < pb:
897 a, pa, b, pb = b, pb, a, pa
898 while pb:
899 # x**pa - a = 0; x**pb - b = 0
900 # x**pa - a = x**(q*pb + r) - a = (x**pb)**q * x**r - a =
901 # b**q * x**r - a; x**r - c = 0; c = b**-q * a mod p
902 q, r = divmod(pa, pb)
903 c = pow(b, q, p)
904 c = igcdex(c, p)[0]
905 c = (c * a) % p
906 pa, pb = pb, r
907 a, b = b, c
908 if pa == 1:
909 if all_roots:
910 res = [a]
911 else:
912 res = a
913 elif pa == 2:
914 return sqrt_mod(a, p, all_roots)
915 else:
916 res = _nthroot_mod1(a, pa, p, all_roots)
917 return res
920def quadratic_residues(p) -> list[int]:
921 """
922 Returns the list of quadratic residues.
924 Examples
925 ========
927 >>> from sympy.ntheory.residue_ntheory import quadratic_residues
928 >>> quadratic_residues(7)
929 [0, 1, 2, 4]
930 """
931 p = as_int(p)
932 r = {pow(i, 2, p) for i in range(p // 2 + 1)}
933 return sorted(r)
936def legendre_symbol(a, p):
937 r"""
938 Returns the Legendre symbol `(a / p)`.
940 For an integer ``a`` and an odd prime ``p``, the Legendre symbol is
941 defined as
943 .. math ::
944 \genfrac(){}{}{a}{p} = \begin{cases}
945 0 & \text{if } p \text{ divides } a\\
946 1 & \text{if } a \text{ is a quadratic residue modulo } p\\
947 -1 & \text{if } a \text{ is a quadratic nonresidue modulo } p
948 \end{cases}
950 Parameters
951 ==========
953 a : integer
954 p : odd prime
956 Examples
957 ========
959 >>> from sympy.ntheory import legendre_symbol
960 >>> [legendre_symbol(i, 7) for i in range(7)]
961 [0, 1, 1, -1, 1, -1, -1]
962 >>> sorted(set([i**2 % 7 for i in range(7)]))
963 [0, 1, 2, 4]
965 See Also
966 ========
968 is_quad_residue, jacobi_symbol
970 """
971 a, p = as_int(a), as_int(p)
972 if not isprime(p) or p == 2:
973 raise ValueError("p should be an odd prime")
974 a = a % p
975 if not a:
976 return 0
977 if pow(a, (p - 1) // 2, p) == 1:
978 return 1
979 return -1
982def jacobi_symbol(m, n):
983 r"""
984 Returns the Jacobi symbol `(m / n)`.
986 For any integer ``m`` and any positive odd integer ``n`` the Jacobi symbol
987 is defined as the product of the Legendre symbols corresponding to the
988 prime factors of ``n``:
990 .. math ::
991 \genfrac(){}{}{m}{n} =
992 \genfrac(){}{}{m}{p^{1}}^{\alpha_1}
993 \genfrac(){}{}{m}{p^{2}}^{\alpha_2}
994 ...
995 \genfrac(){}{}{m}{p^{k}}^{\alpha_k}
996 \text{ where } n =
997 p_1^{\alpha_1}
998 p_2^{\alpha_2}
999 ...
1000 p_k^{\alpha_k}
1002 Like the Legendre symbol, if the Jacobi symbol `\genfrac(){}{}{m}{n} = -1`
1003 then ``m`` is a quadratic nonresidue modulo ``n``.
1005 But, unlike the Legendre symbol, if the Jacobi symbol
1006 `\genfrac(){}{}{m}{n} = 1` then ``m`` may or may not be a quadratic residue
1007 modulo ``n``.
1009 Parameters
1010 ==========
1012 m : integer
1013 n : odd positive integer
1015 Examples
1016 ========
1018 >>> from sympy.ntheory import jacobi_symbol, legendre_symbol
1019 >>> from sympy import S
1020 >>> jacobi_symbol(45, 77)
1021 -1
1022 >>> jacobi_symbol(60, 121)
1023 1
1025 The relationship between the ``jacobi_symbol`` and ``legendre_symbol`` can
1026 be demonstrated as follows:
1028 >>> L = legendre_symbol
1029 >>> S(45).factors()
1030 {3: 2, 5: 1}
1031 >>> jacobi_symbol(7, 45) == L(7, 3)**2 * L(7, 5)**1
1032 True
1034 See Also
1035 ========
1037 is_quad_residue, legendre_symbol
1038 """
1039 m, n = as_int(m), as_int(n)
1040 if n < 0 or not n % 2:
1041 raise ValueError("n should be an odd positive integer")
1042 if m < 0 or m > n:
1043 m %= n
1044 if not m:
1045 return int(n == 1)
1046 if n == 1 or m == 1:
1047 return 1
1048 if igcd(m, n) != 1:
1049 return 0
1051 j = 1
1052 while m != 0:
1053 while m % 2 == 0 and m > 0:
1054 m >>= 1
1055 if n % 8 in [3, 5]:
1056 j = -j
1057 m, n = n, m
1058 if m % 4 == n % 4 == 3:
1059 j = -j
1060 m %= n
1061 return j
1064class mobius(Function):
1065 """
1066 Mobius function maps natural number to {-1, 0, 1}
1068 It is defined as follows:
1069 1) `1` if `n = 1`.
1070 2) `0` if `n` has a squared prime factor.
1071 3) `(-1)^k` if `n` is a square-free positive integer with `k`
1072 number of prime factors.
1074 It is an important multiplicative function in number theory
1075 and combinatorics. It has applications in mathematical series,
1076 algebraic number theory and also physics (Fermion operator has very
1077 concrete realization with Mobius Function model).
1079 Parameters
1080 ==========
1082 n : positive integer
1084 Examples
1085 ========
1087 >>> from sympy.ntheory import mobius
1088 >>> mobius(13*7)
1089 1
1090 >>> mobius(1)
1091 1
1092 >>> mobius(13*7*5)
1093 -1
1094 >>> mobius(13**2)
1095 0
1097 References
1098 ==========
1100 .. [1] https://en.wikipedia.org/wiki/M%C3%B6bius_function
1101 .. [2] Thomas Koshy "Elementary Number Theory with Applications"
1103 """
1104 @classmethod
1105 def eval(cls, n):
1106 if n.is_integer:
1107 if n.is_positive is not True:
1108 raise ValueError("n should be a positive integer")
1109 else:
1110 raise TypeError("n should be an integer")
1111 if n.is_prime:
1112 return S.NegativeOne
1113 elif n is S.One:
1114 return S.One
1115 elif n.is_Integer:
1116 a = factorint(n)
1117 if any(i > 1 for i in a.values()):
1118 return S.Zero
1119 return S.NegativeOne**len(a)
1122def _discrete_log_trial_mul(n, a, b, order=None):
1123 """
1124 Trial multiplication algorithm for computing the discrete logarithm of
1125 ``a`` to the base ``b`` modulo ``n``.
1127 The algorithm finds the discrete logarithm using exhaustive search. This
1128 naive method is used as fallback algorithm of ``discrete_log`` when the
1129 group order is very small.
1131 Examples
1132 ========
1134 >>> from sympy.ntheory.residue_ntheory import _discrete_log_trial_mul
1135 >>> _discrete_log_trial_mul(41, 15, 7)
1136 3
1138 See Also
1139 ========
1141 discrete_log
1143 References
1144 ==========
1146 .. [1] "Handbook of applied cryptography", Menezes, A. J., Van, O. P. C., &
1147 Vanstone, S. A. (1997).
1148 """
1149 a %= n
1150 b %= n
1151 if order is None:
1152 order = n
1153 x = 1
1154 for i in range(order):
1155 if x == a:
1156 return i
1157 x = x * b % n
1158 raise ValueError("Log does not exist")
1161def _discrete_log_shanks_steps(n, a, b, order=None):
1162 """
1163 Baby-step giant-step algorithm for computing the discrete logarithm of
1164 ``a`` to the base ``b`` modulo ``n``.
1166 The algorithm is a time-memory trade-off of the method of exhaustive
1167 search. It uses `O(sqrt(m))` memory, where `m` is the group order.
1169 Examples
1170 ========
1172 >>> from sympy.ntheory.residue_ntheory import _discrete_log_shanks_steps
1173 >>> _discrete_log_shanks_steps(41, 15, 7)
1174 3
1176 See Also
1177 ========
1179 discrete_log
1181 References
1182 ==========
1184 .. [1] "Handbook of applied cryptography", Menezes, A. J., Van, O. P. C., &
1185 Vanstone, S. A. (1997).
1186 """
1187 a %= n
1188 b %= n
1189 if order is None:
1190 order = n_order(b, n)
1191 m = isqrt(order) + 1
1192 T = {}
1193 x = 1
1194 for i in range(m):
1195 T[x] = i
1196 x = x * b % n
1197 z = mod_inverse(b, n)
1198 z = pow(z, m, n)
1199 x = a
1200 for i in range(m):
1201 if x in T:
1202 return i * m + T[x]
1203 x = x * z % n
1204 raise ValueError("Log does not exist")
1207def _discrete_log_pollard_rho(n, a, b, order=None, retries=10, rseed=None):
1208 """
1209 Pollard's Rho algorithm for computing the discrete logarithm of ``a`` to
1210 the base ``b`` modulo ``n``.
1212 It is a randomized algorithm with the same expected running time as
1213 ``_discrete_log_shanks_steps``, but requires a negligible amount of memory.
1215 Examples
1216 ========
1218 >>> from sympy.ntheory.residue_ntheory import _discrete_log_pollard_rho
1219 >>> _discrete_log_pollard_rho(227, 3**7, 3)
1220 7
1222 See Also
1223 ========
1225 discrete_log
1227 References
1228 ==========
1230 .. [1] "Handbook of applied cryptography", Menezes, A. J., Van, O. P. C., &
1231 Vanstone, S. A. (1997).
1232 """
1233 a %= n
1234 b %= n
1236 if order is None:
1237 order = n_order(b, n)
1238 randint = _randint(rseed)
1240 for i in range(retries):
1241 aa = randint(1, order - 1)
1242 ba = randint(1, order - 1)
1243 xa = pow(b, aa, n) * pow(a, ba, n) % n
1245 c = xa % 3
1246 if c == 0:
1247 xb = a * xa % n
1248 ab = aa
1249 bb = (ba + 1) % order
1250 elif c == 1:
1251 xb = xa * xa % n
1252 ab = (aa + aa) % order
1253 bb = (ba + ba) % order
1254 else:
1255 xb = b * xa % n
1256 ab = (aa + 1) % order
1257 bb = ba
1259 for j in range(order):
1260 c = xa % 3
1261 if c == 0:
1262 xa = a * xa % n
1263 ba = (ba + 1) % order
1264 elif c == 1:
1265 xa = xa * xa % n
1266 aa = (aa + aa) % order
1267 ba = (ba + ba) % order
1268 else:
1269 xa = b * xa % n
1270 aa = (aa + 1) % order
1272 c = xb % 3
1273 if c == 0:
1274 xb = a * xb % n
1275 bb = (bb + 1) % order
1276 elif c == 1:
1277 xb = xb * xb % n
1278 ab = (ab + ab) % order
1279 bb = (bb + bb) % order
1280 else:
1281 xb = b * xb % n
1282 ab = (ab + 1) % order
1284 c = xb % 3
1285 if c == 0:
1286 xb = a * xb % n
1287 bb = (bb + 1) % order
1288 elif c == 1:
1289 xb = xb * xb % n
1290 ab = (ab + ab) % order
1291 bb = (bb + bb) % order
1292 else:
1293 xb = b * xb % n
1294 ab = (ab + 1) % order
1296 if xa == xb:
1297 r = (ba - bb) % order
1298 try:
1299 e = mod_inverse(r, order) * (ab - aa) % order
1300 if (pow(b, e, n) - a) % n == 0:
1301 return e
1302 except ValueError:
1303 pass
1304 break
1305 raise ValueError("Pollard's Rho failed to find logarithm")
1308def _discrete_log_pohlig_hellman(n, a, b, order=None):
1309 """
1310 Pohlig-Hellman algorithm for computing the discrete logarithm of ``a`` to
1311 the base ``b`` modulo ``n``.
1313 In order to compute the discrete logarithm, the algorithm takes advantage
1314 of the factorization of the group order. It is more efficient when the
1315 group order factors into many small primes.
1317 Examples
1318 ========
1320 >>> from sympy.ntheory.residue_ntheory import _discrete_log_pohlig_hellman
1321 >>> _discrete_log_pohlig_hellman(251, 210, 71)
1322 197
1324 See Also
1325 ========
1327 discrete_log
1329 References
1330 ==========
1332 .. [1] "Handbook of applied cryptography", Menezes, A. J., Van, O. P. C., &
1333 Vanstone, S. A. (1997).
1334 """
1335 from .modular import crt
1336 a %= n
1337 b %= n
1339 if order is None:
1340 order = n_order(b, n)
1342 f = factorint(order)
1343 l = [0] * len(f)
1345 for i, (pi, ri) in enumerate(f.items()):
1346 for j in range(ri):
1347 gj = pow(b, l[i], n)
1348 aj = pow(a * mod_inverse(gj, n), order // pi**(j + 1), n)
1349 bj = pow(b, order // pi, n)
1350 cj = discrete_log(n, aj, bj, pi, True)
1351 l[i] += cj * pi**j
1353 d, _ = crt([pi**ri for pi, ri in f.items()], l)
1354 return d
1357def discrete_log(n, a, b, order=None, prime_order=None):
1358 """
1359 Compute the discrete logarithm of ``a`` to the base ``b`` modulo ``n``.
1361 This is a recursive function to reduce the discrete logarithm problem in
1362 cyclic groups of composite order to the problem in cyclic groups of prime
1363 order.
1365 It employs different algorithms depending on the problem (subgroup order
1366 size, prime order or not):
1368 * Trial multiplication
1369 * Baby-step giant-step
1370 * Pollard's Rho
1371 * Pohlig-Hellman
1373 Examples
1374 ========
1376 >>> from sympy.ntheory import discrete_log
1377 >>> discrete_log(41, 15, 7)
1378 3
1380 References
1381 ==========
1383 .. [1] https://mathworld.wolfram.com/DiscreteLogarithm.html
1384 .. [2] "Handbook of applied cryptography", Menezes, A. J., Van, O. P. C., &
1385 Vanstone, S. A. (1997).
1387 """
1388 n, a, b = as_int(n), as_int(a), as_int(b)
1389 if order is None:
1390 order = n_order(b, n)
1392 if prime_order is None:
1393 prime_order = isprime(order)
1395 if order < 1000:
1396 return _discrete_log_trial_mul(n, a, b, order)
1397 elif prime_order:
1398 if order < 1000000000000:
1399 return _discrete_log_shanks_steps(n, a, b, order)
1400 return _discrete_log_pollard_rho(n, a, b, order)
1402 return _discrete_log_pohlig_hellman(n, a, b, order)
1406def quadratic_congruence(a, b, c, p):
1407 """
1408 Find the solutions to ``a x**2 + b x + c = 0 mod p.
1410 Parameters
1411 ==========
1413 a : int
1414 b : int
1415 c : int
1416 p : int
1417 A positive integer.
1418 """
1419 a = as_int(a)
1420 b = as_int(b)
1421 c = as_int(c)
1422 p = as_int(p)
1423 a = a % p
1424 b = b % p
1425 c = c % p
1427 if a == 0:
1428 return linear_congruence(b, -c, p)
1429 if p == 2:
1430 roots = []
1431 if c % 2 == 0:
1432 roots.append(0)
1433 if (a + b + c) % 2 == 0:
1434 roots.append(1)
1435 return roots
1436 if isprime(p):
1437 inv_a = mod_inverse(a, p)
1438 b *= inv_a
1439 c *= inv_a
1440 if b % 2 == 1:
1441 b = b + p
1442 d = ((b * b) // 4 - c) % p
1443 y = sqrt_mod(d, p, all_roots=True)
1444 res = set()
1445 for i in y:
1446 res.add((i - b // 2) % p)
1447 return sorted(res)
1448 y = sqrt_mod(b * b - 4 * a * c, 4 * a * p, all_roots=True)
1449 res = set()
1450 for i in y:
1451 root = linear_congruence(2 * a, i - b, 4 * a * p)
1452 for j in root:
1453 res.add(j % p)
1454 return sorted(res)
1457def _polynomial_congruence_prime(coefficients, p):
1458 """A helper function used by polynomial_congruence.
1459 It returns the root of a polynomial modulo prime number
1460 by naive search from [0, p).
1462 Parameters
1463 ==========
1465 coefficients : list of integers
1466 p : prime number
1467 """
1469 roots = []
1470 rank = len(coefficients)
1471 for i in range(0, p):
1472 f_val = 0
1473 for coeff in range(0,rank - 1):
1474 f_val = (f_val + pow(i, int(rank - coeff - 1), p) * coefficients[coeff]) % p
1475 f_val = f_val + coefficients[-1]
1476 if f_val % p == 0:
1477 roots.append(i)
1478 return roots
1481def _diff_poly(root, coefficients, p):
1482 """A helper function used by polynomial_congruence.
1483 It returns the derivative of the polynomial evaluated at the
1484 root (mod p).
1486 Parameters
1487 ==========
1489 coefficients : list of integers
1490 p : prime number
1491 root : integer
1492 """
1494 diff = 0
1495 rank = len(coefficients)
1496 for coeff in range(0, rank - 1):
1497 if not coefficients[coeff]:
1498 continue
1499 diff = (diff + pow(root, rank - coeff - 2, p)*(rank - coeff - 1)*
1500 coefficients[coeff]) % p
1501 return diff % p
1504def _val_poly(root, coefficients, p):
1505 """A helper function used by polynomial_congruence.
1506 It returns value of the polynomial at root (mod p).
1508 Parameters
1509 ==========
1511 coefficients : list of integers
1512 p : prime number
1513 root : integer
1514 """
1515 rank = len(coefficients)
1516 f_val = 0
1517 for coeff in range(0, rank - 1):
1518 f_val = (f_val + pow(root, rank - coeff - 1, p)*
1519 coefficients[coeff]) % p
1520 f_val = f_val + coefficients[-1]
1521 return f_val % p
1524def _valid_expr(expr):
1525 """
1526 return coefficients of expr if it is a univariate polynomial
1527 with integer coefficients else raise a ValueError.
1528 """
1530 if not expr.is_polynomial():
1531 raise ValueError("The expression should be a polynomial")
1532 polynomial = Poly(expr)
1533 if not polynomial.is_univariate:
1534 raise ValueError("The expression should be univariate")
1535 if not polynomial.domain == ZZ:
1536 raise ValueError("The expression should should have integer coefficients")
1537 return polynomial.all_coeffs()
1540def polynomial_congruence(expr, m):
1541 """
1542 Find the solutions to a polynomial congruence equation modulo m.
1544 Parameters
1545 ==========
1547 coefficients : Coefficients of the Polynomial
1548 m : positive integer
1550 Examples
1551 ========
1553 >>> from sympy.ntheory import polynomial_congruence
1554 >>> from sympy.abc import x
1555 >>> expr = x**6 - 2*x**5 -35
1556 >>> polynomial_congruence(expr, 6125)
1557 [3257]
1558 """
1559 coefficients = _valid_expr(expr)
1560 coefficients = [num % m for num in coefficients]
1561 rank = len(coefficients)
1562 if rank == 3:
1563 return quadratic_congruence(*coefficients, m)
1564 if rank == 2:
1565 return quadratic_congruence(0, *coefficients, m)
1566 if coefficients[0] == 1 and 1 + coefficients[-1] == sum(coefficients):
1567 return nthroot_mod(-coefficients[-1], rank - 1, m, True)
1568 if isprime(m):
1569 return _polynomial_congruence_prime(coefficients, m)
1570 return _help(m,
1571 lambda p: _polynomial_congruence_prime(coefficients, p),
1572 lambda root, p: _diff_poly(root, coefficients, p),
1573 lambda root, p: _val_poly(root, coefficients, p))