Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.11/site-packages/ecdsa/numbertheory.py: 19%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1#! /usr/bin/env python
2#
3# Provide some simple capabilities from number theory.
4#
5# Version of 2008.11.14.
6#
7# Written in 2005 and 2006 by Peter Pearson and placed in the public domain.
8# Revision history:
9# 2008.11.14: Use pow(base, exponent, modulus) for modular_exp.
10# Make gcd and lcm accept arbitrarily many arguments.
12from __future__ import division
14import sys
15from six import integer_types, PY2
16from six.moves import reduce
18try:
19 xrange
20except NameError:
21 xrange = range
22try:
23 from gmpy2 import powmod, mpz
25 GMPY2 = True
26 GMPY = False
27except ImportError: # pragma: no branch
28 GMPY2 = False
29 try:
30 from gmpy import mpz
32 GMPY = True
33 except ImportError:
34 GMPY = False
37if GMPY2 or GMPY: # pragma: no branch
38 integer_types = tuple(integer_types + (type(mpz(1)),))
41import math
42import warnings
43import random
44from .util import bit_length
47class Error(Exception):
48 """Base class for exceptions in this module."""
50 pass
53class JacobiError(Error):
54 pass
57class SquareRootError(Error):
58 pass
61class NegativeExponentError(Error):
62 pass
65def modular_exp(base, exponent, modulus): # pragma: no cover
66 """Raise base to exponent, reducing by modulus"""
67 # deprecated in 0.14
68 warnings.warn(
69 "Function is unused in library code. If you use this code, "
70 "change to pow() builtin.",
71 DeprecationWarning,
72 )
73 if exponent < 0:
74 raise NegativeExponentError(
75 "Negative exponents (%d) not allowed" % exponent
76 )
77 return pow(base, exponent, modulus)
80def polynomial_reduce_mod(poly, polymod, p):
81 """Reduce poly by polymod, integer arithmetic modulo p.
83 Polynomials are represented as lists of coefficients
84 of increasing powers of x."""
86 # This module has been tested only by extensive use
87 # in calculating modular square roots.
89 # Just to make this easy, require a monic polynomial:
90 assert polymod[-1] == 1
92 assert len(polymod) > 1
94 while len(poly) >= len(polymod):
95 if poly[-1] != 0:
96 for i in xrange(2, len(polymod) + 1):
97 poly[-i] = (poly[-i] - poly[-1] * polymod[-i]) % p
98 poly = poly[0:-1]
100 return poly
103def polynomial_multiply_mod(m1, m2, polymod, p):
104 """Polynomial multiplication modulo a polynomial over ints mod p.
106 Polynomials are represented as lists of coefficients
107 of increasing powers of x."""
109 # This is just a seat-of-the-pants implementation.
111 # This module has been tested only by extensive use
112 # in calculating modular square roots.
114 # Initialize the product to zero:
116 prod = (len(m1) + len(m2) - 1) * [0]
118 # Add together all the cross-terms:
120 for i in xrange(len(m1)):
121 for j in xrange(len(m2)):
122 prod[i + j] = (prod[i + j] + m1[i] * m2[j]) % p
124 return polynomial_reduce_mod(prod, polymod, p)
127def polynomial_exp_mod(base, exponent, polymod, p):
128 """Polynomial exponentiation modulo a polynomial over ints mod p.
130 Polynomials are represented as lists of coefficients
131 of increasing powers of x."""
133 # Based on the Handbook of Applied Cryptography, algorithm 2.227.
135 # This module has been tested only by extensive use
136 # in calculating modular square roots.
138 assert exponent < p
140 if exponent == 0:
141 return [1]
143 G = base
144 k = exponent
145 if k % 2 == 1:
146 s = G
147 else:
148 s = [1]
150 while k > 1:
151 k = k // 2
152 G = polynomial_multiply_mod(G, G, polymod, p)
153 if k % 2 == 1:
154 s = polynomial_multiply_mod(G, s, polymod, p)
156 return s
159def jacobi(a, n):
160 """Jacobi symbol"""
162 # Based on the Handbook of Applied Cryptography (HAC), algorithm 2.149.
164 # This function has been tested by comparison with a small
165 # table printed in HAC, and by extensive use in calculating
166 # modular square roots.
168 if not n >= 3:
169 raise JacobiError("n must be larger than 2")
170 if not n % 2 == 1:
171 raise JacobiError("n must be odd")
172 a = a % n
173 if a == 0:
174 return 0
175 if a == 1:
176 return 1
177 a1, e = a, 0
178 while a1 % 2 == 0:
179 a1, e = a1 // 2, e + 1
180 if e % 2 == 0 or n % 8 == 1 or n % 8 == 7:
181 s = 1
182 else:
183 s = -1
184 if a1 == 1:
185 return s
186 if n % 4 == 3 and a1 % 4 == 3:
187 s = -s
188 return s * jacobi(n % a1, a1)
191def square_root_mod_prime(a, p):
192 """Modular square root of a, mod p, p prime."""
194 # Based on the Handbook of Applied Cryptography, algorithms 3.34 to 3.39.
196 # This module has been tested for all values in [0,p-1] for
197 # every prime p from 3 to 1229.
199 assert 0 <= a < p
200 assert 1 < p
202 if a == 0:
203 return 0
204 if p == 2:
205 return a
207 jac = jacobi(a, p)
208 if jac == -1:
209 raise SquareRootError("%d has no square root modulo %d" % (a, p))
211 if p % 4 == 3:
212 return pow(a, (p + 1) // 4, p)
214 if p % 8 == 5:
215 d = pow(a, (p - 1) // 4, p)
216 if d == 1:
217 return pow(a, (p + 3) // 8, p)
218 assert d == p - 1
219 return (2 * a * pow(4 * a, (p - 5) // 8, p)) % p
221 if PY2:
222 # xrange on python2 can take integers representable as C long only
223 range_top = min(0x7FFFFFFF, p)
224 else:
225 range_top = p
226 for b in xrange(2, range_top): # pragma: no branch
227 if jacobi(b * b - 4 * a, p) == -1:
228 f = (a, -b, 1)
229 ff = polynomial_exp_mod((0, 1), (p + 1) // 2, f, p)
230 if ff[1]:
231 raise SquareRootError("p is not prime")
232 return ff[0]
233 # just an assertion
234 raise RuntimeError("No b found.") # pragma: no cover
237# because all the inverse_mod code is arch/environment specific, and coveralls
238# expects it to execute equal number of times, we need to waive it by
239# adding the "no branch" pragma to all branches
240if GMPY2: # pragma: no branch
242 def inverse_mod(a, m):
243 """Inverse of a mod m."""
244 if a == 0: # pragma: no branch
245 return 0
246 return powmod(a, -1, m)
248elif GMPY: # pragma: no branch
250 def inverse_mod(a, m):
251 """Inverse of a mod m."""
252 # while libgmp does support inverses modulo, it is accessible
253 # only using the native `pow()` function, and `pow()` in gmpy sanity
254 # checks the parameters before passing them on to underlying
255 # implementation
256 if a == 0: # pragma: no branch
257 return 0
258 a = mpz(a)
259 m = mpz(m)
261 lm, hm = mpz(1), mpz(0)
262 low, high = a % m, m
263 while low > 1: # pragma: no branch
264 r = high // low
265 lm, low, hm, high = hm - lm * r, high - low * r, lm, low
267 return lm % m
269elif sys.version_info >= (3, 8): # pragma: no branch
271 def inverse_mod(a, m):
272 """Inverse of a mod m."""
273 if a == 0: # pragma: no branch
274 return 0
275 return pow(a, -1, m)
277else: # pragma: no branch
279 def inverse_mod(a, m):
280 """Inverse of a mod m."""
282 if a == 0: # pragma: no branch
283 return 0
285 lm, hm = 1, 0
286 low, high = a % m, m
287 while low > 1: # pragma: no branch
288 r = high // low
289 lm, low, hm, high = hm - lm * r, high - low * r, lm, low
291 return lm % m
294try:
295 gcd2 = math.gcd
296except AttributeError:
298 def gcd2(a, b):
299 """Greatest common divisor using Euclid's algorithm."""
300 while a:
301 a, b = b % a, a
302 return b
305def gcd(*a):
306 """Greatest common divisor.
308 Usage: gcd([ 2, 4, 6 ])
309 or: gcd(2, 4, 6)
310 """
312 if len(a) > 1:
313 return reduce(gcd2, a)
314 if hasattr(a[0], "__iter__"):
315 return reduce(gcd2, a[0])
316 return a[0]
319def lcm2(a, b):
320 """Least common multiple of two integers."""
322 return (a * b) // gcd(a, b)
325def lcm(*a):
326 """Least common multiple.
328 Usage: lcm([ 3, 4, 5 ])
329 or: lcm(3, 4, 5)
330 """
332 if len(a) > 1:
333 return reduce(lcm2, a)
334 if hasattr(a[0], "__iter__"):
335 return reduce(lcm2, a[0])
336 return a[0]
339def factorization(n):
340 """Decompose n into a list of (prime,exponent) pairs."""
342 assert isinstance(n, integer_types)
344 if n < 2:
345 return []
347 result = []
349 # Test the small primes:
351 for d in smallprimes:
352 if d > n:
353 break
354 q, r = divmod(n, d)
355 if r == 0:
356 count = 1
357 while d <= n: # pragma: no branch
358 n = q
359 q, r = divmod(n, d)
360 if r != 0:
361 break
362 count = count + 1
363 result.append((d, count))
365 # If n is still greater than the last of our small primes,
366 # it may require further work:
368 if n > smallprimes[-1]:
369 if is_prime(n): # If what's left is prime, it's easy:
370 result.append((n, 1))
371 else: # Ugh. Search stupidly for a divisor:
372 d = smallprimes[-1]
373 while 1:
374 d = d + 2 # Try the next divisor.
375 q, r = divmod(n, d)
376 if q < d: # n < d*d means we're done, n = 1 or prime.
377 break
378 if r == 0: # d divides n. How many times?
379 count = 1
380 n = q
381 # As long as d might still divide n,
382 while d <= n: # pragma: no branch
383 q, r = divmod(n, d) # see if it does.
384 if r != 0:
385 break
386 n = q # It does. Reduce n, increase count.
387 count = count + 1
388 result.append((d, count))
389 if n > 1:
390 result.append((n, 1))
392 return result
395def phi(n): # pragma: no cover
396 """Return the Euler totient function of n."""
397 # deprecated in 0.14
398 warnings.warn(
399 "Function is unused by library code. If you use this code, "
400 "please open an issue in "
401 "https://github.com/tlsfuzzer/python-ecdsa",
402 DeprecationWarning,
403 )
405 assert isinstance(n, integer_types)
407 if n < 3:
408 return 1
410 result = 1
411 ff = factorization(n)
412 for f in ff:
413 e = f[1]
414 if e > 1:
415 result = result * f[0] ** (e - 1) * (f[0] - 1)
416 else:
417 result = result * (f[0] - 1)
418 return result
421def carmichael(n): # pragma: no cover
422 """Return Carmichael function of n.
424 Carmichael(n) is the smallest integer x such that
425 m**x = 1 mod n for all m relatively prime to n.
426 """
427 # deprecated in 0.14
428 warnings.warn(
429 "Function is unused by library code. If you use this code, "
430 "please open an issue in "
431 "https://github.com/tlsfuzzer/python-ecdsa",
432 DeprecationWarning,
433 )
435 return carmichael_of_factorized(factorization(n))
438def carmichael_of_factorized(f_list): # pragma: no cover
439 """Return the Carmichael function of a number that is
440 represented as a list of (prime,exponent) pairs.
441 """
442 # deprecated in 0.14
443 warnings.warn(
444 "Function is unused by library code. If you use this code, "
445 "please open an issue in "
446 "https://github.com/tlsfuzzer/python-ecdsa",
447 DeprecationWarning,
448 )
450 if len(f_list) < 1:
451 return 1
453 result = carmichael_of_ppower(f_list[0])
454 for i in xrange(1, len(f_list)):
455 result = lcm(result, carmichael_of_ppower(f_list[i]))
457 return result
460def carmichael_of_ppower(pp): # pragma: no cover
461 """Carmichael function of the given power of the given prime."""
462 # deprecated in 0.14
463 warnings.warn(
464 "Function is unused by library code. If you use this code, "
465 "please open an issue in "
466 "https://github.com/tlsfuzzer/python-ecdsa",
467 DeprecationWarning,
468 )
470 p, a = pp
471 if p == 2 and a > 2:
472 return 2 ** (a - 2)
473 else:
474 return (p - 1) * p ** (a - 1)
477def order_mod(x, m): # pragma: no cover
478 """Return the order of x in the multiplicative group mod m."""
479 # deprecated in 0.14
480 warnings.warn(
481 "Function is unused by library code. If you use this code, "
482 "please open an issue in "
483 "https://github.com/tlsfuzzer/python-ecdsa",
484 DeprecationWarning,
485 )
487 # Warning: this implementation is not very clever, and will
488 # take a long time if m is very large.
490 if m <= 1:
491 return 0
493 assert gcd(x, m) == 1
495 z = x
496 result = 1
497 while z != 1:
498 z = (z * x) % m
499 result = result + 1
500 return result
503def largest_factor_relatively_prime(a, b): # pragma: no cover
504 """Return the largest factor of a relatively prime to b."""
505 # deprecated in 0.14
506 warnings.warn(
507 "Function is unused by library code. If you use this code, "
508 "please open an issue in "
509 "https://github.com/tlsfuzzer/python-ecdsa",
510 DeprecationWarning,
511 )
513 while 1:
514 d = gcd(a, b)
515 if d <= 1:
516 break
517 b = d
518 while 1:
519 q, r = divmod(a, d)
520 if r > 0:
521 break
522 a = q
523 return a
526def kinda_order_mod(x, m): # pragma: no cover
527 """Return the order of x in the multiplicative group mod m',
528 where m' is the largest factor of m relatively prime to x.
529 """
530 # deprecated in 0.14
531 warnings.warn(
532 "Function is unused by library code. If you use this code, "
533 "please open an issue in "
534 "https://github.com/tlsfuzzer/python-ecdsa",
535 DeprecationWarning,
536 )
538 return order_mod(x, largest_factor_relatively_prime(m, x))
541def is_prime(n):
542 """Return True if x is prime, False otherwise.
544 We use the Miller-Rabin test, as given in Menezes et al. p. 138.
545 This test is not exact: there are composite values n for which
546 it returns True.
548 In testing the odd numbers from 10000001 to 19999999,
549 about 66 composites got past the first test,
550 5 got past the second test, and none got past the third.
551 Since factors of 2, 3, 5, 7, and 11 were detected during
552 preliminary screening, the number of numbers tested by
553 Miller-Rabin was (19999999 - 10000001)*(2/3)*(4/5)*(6/7)
554 = 4.57 million.
555 """
557 # (This is used to study the risk of false positives:)
558 global miller_rabin_test_count
560 miller_rabin_test_count = 0
562 if n <= smallprimes[-1]:
563 if n in smallprimes:
564 return True
565 else:
566 return False
567 # 2310 = 2 * 3 * 5 * 7 * 11
568 if gcd(n, 2310) != 1:
569 return False
571 # Choose a number of iterations sufficient to reduce the
572 # probability of accepting a composite below 2**-80
573 # (from Menezes et al. Table 4.4):
575 t = 40
576 n_bits = 1 + bit_length(n)
577 assert 11 <= n_bits <= 16384
578 for k, tt in (
579 (100, 27),
580 (150, 18),
581 (200, 15),
582 (250, 12),
583 (300, 9),
584 (350, 8),
585 (400, 7),
586 (450, 6),
587 (550, 5),
588 (650, 4),
589 (850, 3),
590 (1300, 2),
591 ):
592 if n_bits < k:
593 break
594 t = tt
596 # Run the test t times:
598 s = 0
599 r = n - 1
600 while (r % 2) == 0:
601 s = s + 1
602 r = r // 2
603 for i in xrange(t):
604 a = random.choice(smallprimes)
605 y = pow(a, r, n)
606 if y != 1 and y != n - 1:
607 j = 1
608 while j <= s - 1 and y != n - 1:
609 y = pow(y, 2, n)
610 if y == 1:
611 miller_rabin_test_count = i + 1
612 return False
613 j = j + 1
614 if y != n - 1:
615 miller_rabin_test_count = i + 1
616 return False
617 return True
620def next_prime(starting_value):
621 """Return the smallest prime larger than the starting value."""
623 if starting_value < 2:
624 return 2
625 result = (starting_value + 1) | 1
626 while not is_prime(result):
627 result = result + 2
628 return result
631smallprimes = [
632 2,
633 3,
634 5,
635 7,
636 11,
637 13,
638 17,
639 19,
640 23,
641 29,
642 31,
643 37,
644 41,
645 43,
646 47,
647 53,
648 59,
649 61,
650 67,
651 71,
652 73,
653 79,
654 83,
655 89,
656 97,
657 101,
658 103,
659 107,
660 109,
661 113,
662 127,
663 131,
664 137,
665 139,
666 149,
667 151,
668 157,
669 163,
670 167,
671 173,
672 179,
673 181,
674 191,
675 193,
676 197,
677 199,
678 211,
679 223,
680 227,
681 229,
682 233,
683 239,
684 241,
685 251,
686 257,
687 263,
688 269,
689 271,
690 277,
691 281,
692 283,
693 293,
694 307,
695 311,
696 313,
697 317,
698 331,
699 337,
700 347,
701 349,
702 353,
703 359,
704 367,
705 373,
706 379,
707 383,
708 389,
709 397,
710 401,
711 409,
712 419,
713 421,
714 431,
715 433,
716 439,
717 443,
718 449,
719 457,
720 461,
721 463,
722 467,
723 479,
724 487,
725 491,
726 499,
727 503,
728 509,
729 521,
730 523,
731 541,
732 547,
733 557,
734 563,
735 569,
736 571,
737 577,
738 587,
739 593,
740 599,
741 601,
742 607,
743 613,
744 617,
745 619,
746 631,
747 641,
748 643,
749 647,
750 653,
751 659,
752 661,
753 673,
754 677,
755 683,
756 691,
757 701,
758 709,
759 719,
760 727,
761 733,
762 739,
763 743,
764 751,
765 757,
766 761,
767 769,
768 773,
769 787,
770 797,
771 809,
772 811,
773 821,
774 823,
775 827,
776 829,
777 839,
778 853,
779 857,
780 859,
781 863,
782 877,
783 881,
784 883,
785 887,
786 907,
787 911,
788 919,
789 929,
790 937,
791 941,
792 947,
793 953,
794 967,
795 971,
796 977,
797 983,
798 991,
799 997,
800 1009,
801 1013,
802 1019,
803 1021,
804 1031,
805 1033,
806 1039,
807 1049,
808 1051,
809 1061,
810 1063,
811 1069,
812 1087,
813 1091,
814 1093,
815 1097,
816 1103,
817 1109,
818 1117,
819 1123,
820 1129,
821 1151,
822 1153,
823 1163,
824 1171,
825 1181,
826 1187,
827 1193,
828 1201,
829 1213,
830 1217,
831 1223,
832 1229,
833]
835miller_rabin_test_count = 0