Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/ecdsa/numbertheory.py: 33%

1#! /usr/bin/env python

3# Provide some simple capabilities from number theory.

5# Version of 2008.11.14.

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

25 GMPY2 = True

26 GMPY = False

27except ImportError:

28 GMPY2 = False

29 try:

30 from gmpy import mpz

32 GMPY = True

33 except ImportError:

34 GMPY = False

36import math

37import warnings

40class Error(Exception):

41 """Base class for exceptions in this module."""

43 pass

46class JacobiError(Error):

47 pass

50class SquareRootError(Error):

51 pass

54class NegativeExponentError(Error):

55 pass

58def modular_exp(base, exponent, modulus): # pragma: no cover

59 """Raise base to exponent, reducing by modulus"""

60 # deprecated in 0.14

61 warnings.warn(

62 "Function is unused in library code. If you use this code, "

63 "change to pow() builtin.",

64 DeprecationWarning,

65 )

66 if exponent < 0:

67 raise NegativeExponentError(

68 "Negative exponents (%d) not allowed" % exponent

69 )

70 return pow(base, exponent, modulus)

73def polynomial_reduce_mod(poly, polymod, p):

74 """Reduce poly by polymod, integer arithmetic modulo p.

76 Polynomials are represented as lists of coefficients

77 of increasing powers of x."""

79 # This module has been tested only by extensive use

80 # in calculating modular square roots.

82 # Just to make this easy, require a monic polynomial:

83 assert polymod[-1] == 1

85 assert len(polymod) > 1

87 while len(poly) >= len(polymod):

88 if poly[-1] != 0:

89 for i in xrange(2, len(polymod) + 1):

90 poly[-i] = (poly[-i] - poly[-1] * polymod[-i]) % p

91 poly = poly[0:-1]

93 return poly

96def polynomial_multiply_mod(m1, m2, polymod, p):

97 """Polynomial multiplication modulo a polynomial over ints mod p.

99 Polynomials are represented as lists of coefficients

100 of increasing powers of x."""

101

102 # This is just a seat-of-the-pants implementation.

103

104 # This module has been tested only by extensive use

105 # in calculating modular square roots.

106

107 # Initialize the product to zero:

108

109 prod = (len(m1) + len(m2) - 1) * [0]

110

111 # Add together all the cross-terms:

112

113 for i in xrange(len(m1)):

114 for j in xrange(len(m2)):

115 prod[i + j] = (prod[i + j] + m1[i] * m2[j]) % p

116

117 return polynomial_reduce_mod(prod, polymod, p)

118

119

120def polynomial_exp_mod(base, exponent, polymod, p):

121 """Polynomial exponentiation modulo a polynomial over ints mod p.

122

123 Polynomials are represented as lists of coefficients

124 of increasing powers of x."""

125

126 # Based on the Handbook of Applied Cryptography, algorithm 2.227.

127

128 # This module has been tested only by extensive use

129 # in calculating modular square roots.

130

131 assert exponent < p

132

133 if exponent == 0:

134 return [1]

135

136 G = base

137 k = exponent

138 if k % 2 == 1:

139 s = G

140 else:

141 s = [1]

142

143 while k > 1:

144 k = k // 2

145 G = polynomial_multiply_mod(G, G, polymod, p)

146 if k % 2 == 1:

147 s = polynomial_multiply_mod(G, s, polymod, p)

148

149 return s

150

151

152def jacobi(a, n):

153 """Jacobi symbol"""

154

155 # Based on the Handbook of Applied Cryptography (HAC), algorithm 2.149.

156

157 # This function has been tested by comparison with a small

158 # table printed in HAC, and by extensive use in calculating

159 # modular square roots.

160

161 if not n >= 3:

162 raise JacobiError("n must be larger than 2")

163 if not n % 2 == 1:

164 raise JacobiError("n must be odd")

165 a = a % n

166 if a == 0:

167 return 0

168 if a == 1:

169 return 1

170 a1, e = a, 0

171 while a1 % 2 == 0:

172 a1, e = a1 // 2, e + 1

173 if e % 2 == 0 or n % 8 == 1 or n % 8 == 7:

174 s = 1

175 else:

176 s = -1

177 if a1 == 1:

178 return s

179 if n % 4 == 3 and a1 % 4 == 3:

180 s = -s

181 return s * jacobi(n % a1, a1)

182

183

184def square_root_mod_prime(a, p):

185 """Modular square root of a, mod p, p prime."""

186

187 # Based on the Handbook of Applied Cryptography, algorithms 3.34 to 3.39.

188

189 # This module has been tested for all values in [0,p-1] for

190 # every prime p from 3 to 1229.

191

192 assert 0 <= a < p

193 assert 1 < p

194

195 if a == 0:

196 return 0

197 if p == 2:

198 return a

199

200 jac = jacobi(a, p)

201 if jac == -1:

202 raise SquareRootError("%d has no square root modulo %d" % (a, p))

203

204 if p % 4 == 3:

205 return pow(a, (p + 1) // 4, p)

206

207 if p % 8 == 5:

208 d = pow(a, (p - 1) // 4, p)

209 if d == 1:

210 return pow(a, (p + 3) // 8, p)

211 assert d == p - 1

212 return (2 * a * pow(4 * a, (p - 5) // 8, p)) % p

213

214 if PY2:

215 # xrange on python2 can take integers representable as C long only

216 range_top = min(0x7FFFFFFF, p)

217 else:

218 range_top = p

219 for b in xrange(2, range_top):

220 if jacobi(b * b - 4 * a, p) == -1:

221 f = (a, -b, 1)

222 ff = polynomial_exp_mod((0, 1), (p + 1) // 2, f, p)

223 if ff[1]:

224 raise SquareRootError("p is not prime")

225 return ff[0]

226 raise RuntimeError("No b found.")

227

228

229# because all the inverse_mod code is arch/environment specific, and coveralls

230# expects it to execute equal number of times, we need to waive it by

231# adding the "no branch" pragma to all branches

232if GMPY2: # pragma: no branch

233

234 def inverse_mod(a, m):

235 """Inverse of a mod m."""

236 if a == 0: # pragma: no branch

237 return 0

238 return powmod(a, -1, m)

239

240elif GMPY: # pragma: no branch

241

242 def inverse_mod(a, m):

243 """Inverse of a mod m."""

244 # while libgmp does support inverses modulo, it is accessible

245 # only using the native `pow()` function, and `pow()` in gmpy sanity

246 # checks the parameters before passing them on to underlying

247 # implementation

248 if a == 0: # pragma: no branch

249 return 0

250 a = mpz(a)

251 m = mpz(m)

252

253 lm, hm = mpz(1), mpz(0)

254 low, high = a % m, m

255 while low > 1: # pragma: no branch

256 r = high // low

257 lm, low, hm, high = hm - lm * r, high - low * r, lm, low

258

259 return lm % m

260

261elif sys.version_info >= (3, 8): # pragma: no branch

262

263 def inverse_mod(a, m):

264 """Inverse of a mod m."""

265 if a == 0: # pragma: no branch

266 return 0

267 return pow(a, -1, m)

268

269else: # pragma: no branch

270

271 def inverse_mod(a, m):

272 """Inverse of a mod m."""

273

274 if a == 0: # pragma: no branch

275 return 0

276

277 lm, hm = 1, 0

278 low, high = a % m, m

279 while low > 1: # pragma: no branch

280 r = high // low

281 lm, low, hm, high = hm - lm * r, high - low * r, lm, low

282

283 return lm % m

284

285

286try:

287 gcd2 = math.gcd

288except AttributeError:

289

290 def gcd2(a, b):

291 """Greatest common divisor using Euclid's algorithm."""

292 while a:

293 a, b = b % a, a

294 return b

295

296

297def gcd(*a):

298 """Greatest common divisor.

299

300 Usage: gcd([ 2, 4, 6 ])

301 or: gcd(2, 4, 6)

302 """

303

304 if len(a) > 1:

305 return reduce(gcd2, a)

306 if hasattr(a[0], "__iter__"):

307 return reduce(gcd2, a[0])

308 return a[0]

309

310

311def lcm2(a, b):

312 """Least common multiple of two integers."""

313

314 return (a * b) // gcd(a, b)

315

316

317def lcm(*a):

318 """Least common multiple.

319

320 Usage: lcm([ 3, 4, 5 ])

321 or: lcm(3, 4, 5)

322 """

323

324 if len(a) > 1:

325 return reduce(lcm2, a)

326 if hasattr(a[0], "__iter__"):

327 return reduce(lcm2, a[0])

328 return a[0]

329

330

331def factorization(n):

332 """Decompose n into a list of (prime,exponent) pairs."""

333

334 assert isinstance(n, integer_types)

335

336 if n < 2:

337 return []

338

339 result = []

340

341 # Test the small primes:

342

343 for d in smallprimes:

344 if d > n:

345 break

346 q, r = divmod(n, d)

347 if r == 0:

348 count = 1

349 while d <= n:

350 n = q

351 q, r = divmod(n, d)

352 if r != 0:

353 break

354 count = count + 1

355 result.append((d, count))

356

357 # If n is still greater than the last of our small primes,

358 # it may require further work:

359

360 if n > smallprimes[-1]:

361 if is_prime(n): # If what's left is prime, it's easy:

362 result.append((n, 1))

363 else: # Ugh. Search stupidly for a divisor:

364 d = smallprimes[-1]

365 while 1:

366 d = d + 2 # Try the next divisor.

367 q, r = divmod(n, d)

368 if q < d: # n < d*d means we're done, n = 1 or prime.

369 break

370 if r == 0: # d divides n. How many times?

371 count = 1

372 n = q

373 while d <= n: # As long as d might still divide n,

374 q, r = divmod(n, d) # see if it does.

375 if r != 0:

376 break

377 n = q # It does. Reduce n, increase count.

378 count = count + 1

379 result.append((d, count))

380 if n > 1:

381 result.append((n, 1))

382

383 return result

384

385

386def phi(n): # pragma: no cover

387 """Return the Euler totient function of n."""

388 # deprecated in 0.14

389 warnings.warn(

390 "Function is unused by library code. If you use this code, "

391 "please open an issue in "

392 "https://github.com/tlsfuzzer/python-ecdsa",

393 DeprecationWarning,

394 )

395

396 assert isinstance(n, integer_types)

397

398 if n < 3:

399 return 1

400

401 result = 1

402 ff = factorization(n)

403 for f in ff:

404 e = f[1]

405 if e > 1:

406 result = result * f[0] ** (e - 1) * (f[0] - 1)

407 else:

408 result = result * (f[0] - 1)

409 return result

410

411

412def carmichael(n): # pragma: no cover

413 """Return Carmichael function of n.

414

415 Carmichael(n) is the smallest integer x such that

416 m**x = 1 mod n for all m relatively prime to n.

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/tlsfuzzer/python-ecdsa",

423 DeprecationWarning,

424 )

425

426 return carmichael_of_factorized(factorization(n))

427

428

429def carmichael_of_factorized(f_list): # pragma: no cover

430 """Return the Carmichael function of a number that is

431 represented as a list of (prime,exponent) pairs.

432 """

433 # deprecated in 0.14

434 warnings.warn(

435 "Function is unused by library code. If you use this code, "

436 "please open an issue in "

437 "https://github.com/tlsfuzzer/python-ecdsa",

438 DeprecationWarning,

439 )

440

441 if len(f_list) < 1:

442 return 1

443

444 result = carmichael_of_ppower(f_list[0])

445 for i in xrange(1, len(f_list)):

446 result = lcm(result, carmichael_of_ppower(f_list[i]))

447

448 return result

449

450

451def carmichael_of_ppower(pp): # pragma: no cover

452 """Carmichael function of the given power of the given prime."""

453 # deprecated in 0.14

454 warnings.warn(

455 "Function is unused by library code. If you use this code, "

456 "please open an issue in "

457 "https://github.com/tlsfuzzer/python-ecdsa",

458 DeprecationWarning,

459 )

460

461 p, a = pp

462 if p == 2 and a > 2:

463 return 2 ** (a - 2)

464 else:

465 return (p - 1) * p ** (a - 1)

466

467

468def order_mod(x, m): # pragma: no cover

469 """Return the order of x in the multiplicative group mod m."""

470 # deprecated in 0.14

471 warnings.warn(

472 "Function is unused by library code. If you use this code, "

473 "please open an issue in "

474 "https://github.com/tlsfuzzer/python-ecdsa",

475 DeprecationWarning,

476 )

477

478 # Warning: this implementation is not very clever, and will

479 # take a long time if m is very large.

480

481 if m <= 1:

482 return 0

483

484 assert gcd(x, m) == 1

485

486 z = x

487 result = 1

488 while z != 1:

489 z = (z * x) % m

490 result = result + 1

491 return result

492

493

494def largest_factor_relatively_prime(a, b): # pragma: no cover

495 """Return the largest factor of a relatively prime to b."""

496 # deprecated in 0.14

497 warnings.warn(

498 "Function is unused by library code. If you use this code, "

499 "please open an issue in "

500 "https://github.com/tlsfuzzer/python-ecdsa",

501 DeprecationWarning,

502 )

503

504 while 1:

505 d = gcd(a, b)

506 if d <= 1:

507 break

508 b = d

509 while 1:

510 q, r = divmod(a, d)

511 if r > 0:

512 break

513 a = q

514 return a

515

516

517def kinda_order_mod(x, m): # pragma: no cover

518 """Return the order of x in the multiplicative group mod m',

519 where m' is the largest factor of m relatively prime to x.

520 """

521 # deprecated in 0.14

522 warnings.warn(

523 "Function is unused by library code. If you use this code, "

524 "please open an issue in "

525 "https://github.com/tlsfuzzer/python-ecdsa",

526 DeprecationWarning,

527 )

528

529 return order_mod(x, largest_factor_relatively_prime(m, x))

530

531

532def is_prime(n):

533 """Return True if x is prime, False otherwise.

534

535 We use the Miller-Rabin test, as given in Menezes et al. p. 138.

536 This test is not exact: there are composite values n for which

537 it returns True.

538

539 In testing the odd numbers from 10000001 to 19999999,

540 about 66 composites got past the first test,

541 5 got past the second test, and none got past the third.

542 Since factors of 2, 3, 5, 7, and 11 were detected during

543 preliminary screening, the number of numbers tested by

544 Miller-Rabin was (19999999 - 10000001)*(2/3)*(4/5)*(6/7)

545 = 4.57 million.

546 """

547

548 # (This is used to study the risk of false positives:)

549 global miller_rabin_test_count

550

551 miller_rabin_test_count = 0

552

553 if n <= smallprimes[-1]:

554 if n in smallprimes:

555 return True

556 else:

557 return False

558

559 if gcd(n, 2 * 3 * 5 * 7 * 11) != 1:

560 return False

561

562 # Choose a number of iterations sufficient to reduce the

563 # probability of accepting a composite below 2**-80

564 # (from Menezes et al. Table 4.4):

565

566 t = 40

567 n_bits = 1 + int(math.log(n, 2))

568 for k, tt in (

569 (100, 27),

570 (150, 18),

571 (200, 15),

572 (250, 12),

573 (300, 9),

574 (350, 8),

575 (400, 7),

576 (450, 6),

577 (550, 5),

578 (650, 4),

579 (850, 3),

580 (1300, 2),

581 ):

582 if n_bits < k:

583 break

584 t = tt

585

586 # Run the test t times:

587

588 s = 0

589 r = n - 1

590 while (r % 2) == 0:

591 s = s + 1

592 r = r // 2

593 for i in xrange(t):

594 a = smallprimes[i]

595 y = pow(a, r, n)

596 if y != 1 and y != n - 1:

597 j = 1

598 while j <= s - 1 and y != n - 1:

599 y = pow(y, 2, n)

600 if y == 1:

601 miller_rabin_test_count = i + 1

602 return False

603 j = j + 1

604 if y != n - 1:

605 miller_rabin_test_count = i + 1

606 return False

607 return True

608

609

610def next_prime(starting_value):

611 """Return the smallest prime larger than the starting value."""

612

613 if starting_value < 2:

614 return 2

615 result = (starting_value + 1) | 1

616 while not is_prime(result):

617 result = result + 2

618 return result

619

620

621smallprimes = [

622 2,

623 3,

624 5,

625 7,

626 11,

627 13,

628 17,

629 19,

630 23,

631 29,

632 31,

633 37,

634 41,

635 43,

636 47,

637 53,

638 59,

639 61,

640 67,

641 71,

642 73,

643 79,

644 83,

645 89,

646 97,

647 101,

648 103,

649 107,

650 109,

651 113,

652 127,

653 131,

654 137,

655 139,

656 149,

657 151,

658 157,

659 163,

660 167,

661 173,

662 179,

663 181,

664 191,

665 193,

666 197,

667 199,

668 211,

669 223,

670 227,

671 229,

672 233,

673 239,

674 241,

675 251,

676 257,

677 263,

678 269,

679 271,

680 277,

681 281,

682 283,

683 293,

684 307,

685 311,

686 313,

687 317,

688 331,

689 337,

690 347,

691 349,

692 353,

693 359,

694 367,

695 373,

696 379,

697 383,

698 389,

699 397,

700 401,

701 409,

702 419,

703 421,

704 431,

705 433,

706 439,

707 443,

708 449,

709 457,

710 461,

711 463,

712 467,

713 479,

714 487,

715 491,

716 499,

717 503,

718 509,

719 521,

720 523,

721 541,

722 547,

723 557,

724 563,

725 569,

726 571,

727 577,

728 587,

729 593,

730 599,

731 601,

732 607,

733 613,

734 617,

735 619,

736 631,

737 641,

738 643,

739 647,

740 653,

741 659,

742 661,

743 673,

744 677,

745 683,

746 691,

747 701,

748 709,

749 719,

750 727,

751 733,

752 739,

753 743,

754 751,

755 757,

756 761,

757 769,

758 773,

759 787,

760 797,

761 809,

762 811,

763 821,

764 823,

765 827,

766 829,

767 839,

768 853,

769 857,

770 859,

771 863,

772 877,

773 881,

774 883,

775 887,

776 907,

777 911,

778 919,

779 929,

780 937,

781 941,

782 947,

783 953,

784 967,

785 971,

786 977,

787 983,

788 991,

789 997,

790 1009,

791 1013,

792 1019,

793 1021,

794 1031,

795 1033,

796 1039,

797 1049,

798 1051,

799 1061,

800 1063,

801 1069,

802 1087,

803 1091,

804 1093,

805 1097,

806 1103,

807 1109,

808 1117,

809 1123,

810 1129,

811 1151,

812 1153,

813 1163,

814 1171,

815 1181,

816 1187,

817 1193,

818 1201,

819 1213,

820 1217,

821 1223,

822 1229,

823]

824

825miller_rabin_test_count = 0