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