prime.py 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201
  1. # -*- coding: utf-8 -*-
  2. #
  3. # Copyright 2011 Sybren A. Stüvel <sybren@stuvel.eu>
  4. #
  5. # Licensed under the Apache License, Version 2.0 (the "License");
  6. # you may not use this file except in compliance with the License.
  7. # You may obtain a copy of the License at
  8. #
  9. # https://www.apache.org/licenses/LICENSE-2.0
  10. #
  11. # Unless required by applicable law or agreed to in writing, software
  12. # distributed under the License is distributed on an "AS IS" BASIS,
  13. # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14. # See the License for the specific language governing permissions and
  15. # limitations under the License.
  16. """Numerical functions related to primes.
  17. Implementation based on the book Algorithm Design by Michael T. Goodrich and
  18. Roberto Tamassia, 2002.
  19. """
  20. from rsa._compat import range
  21. import rsa.common
  22. import rsa.randnum
  23. __all__ = ['getprime', 'are_relatively_prime']
  24. def gcd(p, q):
  25. """Returns the greatest common divisor of p and q
  26. >>> gcd(48, 180)
  27. 12
  28. """
  29. while q != 0:
  30. (p, q) = (q, p % q)
  31. return p
  32. def get_primality_testing_rounds(number):
  33. """Returns minimum number of rounds for Miller-Rabing primality testing,
  34. based on number bitsize.
  35. According to NIST FIPS 186-4, Appendix C, Table C.3, minimum number of
  36. rounds of M-R testing, using an error probability of 2 ** (-100), for
  37. different p, q bitsizes are:
  38. * p, q bitsize: 512; rounds: 7
  39. * p, q bitsize: 1024; rounds: 4
  40. * p, q bitsize: 1536; rounds: 3
  41. See: http://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf
  42. """
  43. # Calculate number bitsize.
  44. bitsize = rsa.common.bit_size(number)
  45. # Set number of rounds.
  46. if bitsize >= 1536:
  47. return 3
  48. if bitsize >= 1024:
  49. return 4
  50. if bitsize >= 512:
  51. return 7
  52. # For smaller bitsizes, set arbitrary number of rounds.
  53. return 10
  54. def miller_rabin_primality_testing(n, k):
  55. """Calculates whether n is composite (which is always correct) or prime
  56. (which theoretically is incorrect with error probability 4**-k), by
  57. applying Miller-Rabin primality testing.
  58. For reference and implementation example, see:
  59. https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
  60. :param n: Integer to be tested for primality.
  61. :type n: int
  62. :param k: Number of rounds (witnesses) of Miller-Rabin testing.
  63. :type k: int
  64. :return: False if the number is composite, True if it's probably prime.
  65. :rtype: bool
  66. """
  67. # prevent potential infinite loop when d = 0
  68. if n < 2:
  69. return False
  70. # Decompose (n - 1) to write it as (2 ** r) * d
  71. # While d is even, divide it by 2 and increase the exponent.
  72. d = n - 1
  73. r = 0
  74. while not (d & 1):
  75. r += 1
  76. d >>= 1
  77. # Test k witnesses.
  78. for _ in range(k):
  79. # Generate random integer a, where 2 <= a <= (n - 2)
  80. a = rsa.randnum.randint(n - 3) + 1
  81. x = pow(a, d, n)
  82. if x == 1 or x == n - 1:
  83. continue
  84. for _ in range(r - 1):
  85. x = pow(x, 2, n)
  86. if x == 1:
  87. # n is composite.
  88. return False
  89. if x == n - 1:
  90. # Exit inner loop and continue with next witness.
  91. break
  92. else:
  93. # If loop doesn't break, n is composite.
  94. return False
  95. return True
  96. def is_prime(number):
  97. """Returns True if the number is prime, and False otherwise.
  98. >>> is_prime(2)
  99. True
  100. >>> is_prime(42)
  101. False
  102. >>> is_prime(41)
  103. True
  104. """
  105. # Check for small numbers.
  106. if number < 10:
  107. return number in {2, 3, 5, 7}
  108. # Check for even numbers.
  109. if not (number & 1):
  110. return False
  111. # Calculate minimum number of rounds.
  112. k = get_primality_testing_rounds(number)
  113. # Run primality testing with (minimum + 1) rounds.
  114. return miller_rabin_primality_testing(number, k + 1)
  115. def getprime(nbits):
  116. """Returns a prime number that can be stored in 'nbits' bits.
  117. >>> p = getprime(128)
  118. >>> is_prime(p-1)
  119. False
  120. >>> is_prime(p)
  121. True
  122. >>> is_prime(p+1)
  123. False
  124. >>> from rsa import common
  125. >>> common.bit_size(p) == 128
  126. True
  127. """
  128. assert nbits > 3 # the loop wil hang on too small numbers
  129. while True:
  130. integer = rsa.randnum.read_random_odd_int(nbits)
  131. # Test for primeness
  132. if is_prime(integer):
  133. return integer
  134. # Retry if not prime
  135. def are_relatively_prime(a, b):
  136. """Returns True if a and b are relatively prime, and False if they
  137. are not.
  138. >>> are_relatively_prime(2, 3)
  139. True
  140. >>> are_relatively_prime(2, 4)
  141. False
  142. """
  143. d = gcd(a, b)
  144. return d == 1
  145. if __name__ == '__main__':
  146. print('Running doctests 1000x or until failure')
  147. import doctest
  148. for count in range(1000):
  149. (failures, tests) = doctest.testmod()
  150. if failures:
  151. break
  152. if count % 100 == 0 and count:
  153. print('%i times' % count)
  154. print('Doctests done')