Primality.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356
  1. # ===================================================================
  2. #
  3. # Copyright (c) 2014, Legrandin <helderijs@gmail.com>
  4. # All rights reserved.
  5. #
  6. # Redistribution and use in source and binary forms, with or without
  7. # modification, are permitted provided that the following conditions
  8. # are met:
  9. #
  10. # 1. Redistributions of source code must retain the above copyright
  11. # notice, this list of conditions and the following disclaimer.
  12. # 2. Redistributions in binary form must reproduce the above copyright
  13. # notice, this list of conditions and the following disclaimer in
  14. # the documentation and/or other materials provided with the
  15. # distribution.
  16. #
  17. # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  18. # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  19. # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  20. # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  21. # COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  22. # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  23. # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  24. # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  25. # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  26. # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  27. # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28. # POSSIBILITY OF SUCH DAMAGE.
  29. # ===================================================================
  30. """Functions to create and test prime numbers.
  31. :undocumented: __package__
  32. """
  33. from Crypto.Math.Numbers import Integer
  34. from Crypto import Random
  35. COMPOSITE = 0
  36. PROBABLY_PRIME = 1
  37. def miller_rabin_test(candidate, iterations, randfunc=None):
  38. """Perform a Miller-Rabin primality test on an integer.
  39. The test is specified in Section C.3.1 of `FIPS PUB 186-4`__.
  40. :Parameters:
  41. candidate : integer
  42. The number to test for primality.
  43. iterations : integer
  44. The maximum number of iterations to perform before
  45. declaring a candidate a probable prime.
  46. randfunc : callable
  47. An RNG function where bases are taken from.
  48. :Returns:
  49. ``Primality.COMPOSITE`` or ``Primality.PROBABLY_PRIME``.
  50. .. __: http://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf
  51. """
  52. if not isinstance(candidate, Integer):
  53. candidate = Integer(candidate)
  54. if candidate.is_even():
  55. return COMPOSITE
  56. one = Integer(1)
  57. minus_one = Integer(candidate - 1)
  58. if randfunc is None:
  59. randfunc = Random.new().read
  60. # Step 1 and 2
  61. m = Integer(minus_one)
  62. a = 0
  63. while m.is_even():
  64. m >>= 1
  65. a += 1
  66. # Skip step 3
  67. # Step 4
  68. for i in xrange(iterations):
  69. # Step 4.1-2
  70. base = 1
  71. while base in (one, minus_one):
  72. base = Integer.random_range(min_inclusive=2,
  73. max_inclusive=candidate - 2)
  74. assert(2 <= base <= candidate - 2)
  75. # Step 4.3-4.4
  76. z = pow(base, m, candidate)
  77. if z in (one, minus_one):
  78. continue
  79. # Step 4.5
  80. for j in xrange(1, a):
  81. z = pow(z, 2, candidate)
  82. if z == minus_one:
  83. break
  84. if z == one:
  85. return COMPOSITE
  86. else:
  87. return COMPOSITE
  88. # Step 5
  89. return PROBABLY_PRIME
  90. def lucas_test(candidate):
  91. """Perform a Lucas primality test on an integer.
  92. The test is specified in Section C.3.3 of `FIPS PUB 186-4`__.
  93. :Parameters:
  94. candidate : integer
  95. The number to test for primality.
  96. :Returns:
  97. ``Primality.COMPOSITE`` or ``Primality.PROBABLY_PRIME``.
  98. .. __: http://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf
  99. """
  100. if not isinstance(candidate, Integer):
  101. candidate = Integer(candidate)
  102. # Step 1
  103. if candidate.is_even() or candidate.is_perfect_square():
  104. return COMPOSITE
  105. # Step 2
  106. def alternate():
  107. sgn = 1
  108. value = 5
  109. for x in xrange(10):
  110. yield sgn * value
  111. sgn, value = -sgn, value + 2
  112. for D in alternate():
  113. js = Integer.jacobi_symbol(D, candidate)
  114. if js == 0:
  115. return COMPOSITE
  116. if js == -1:
  117. break
  118. else:
  119. return COMPOSITE
  120. # Found D. P=1 and Q=(1-D)/4 (note that Q is guaranteed to be an integer)
  121. # Step 3
  122. # This is \delta(n) = n - jacobi(D/n)
  123. K = candidate + 1
  124. # Step 4
  125. r = K.size_in_bits() - 1
  126. # Step 5
  127. # U_1=1 and V_1=P
  128. U_i = Integer(1)
  129. V_i = Integer(1)
  130. U_temp = Integer(0)
  131. V_temp = Integer(0)
  132. # Step 6
  133. for i in xrange(r - 1, -1, -1):
  134. # Square
  135. # U_temp = U_i * V_i % candidate
  136. U_temp.set(U_i)
  137. U_temp *= V_i
  138. U_temp %= candidate
  139. # V_temp = (((V_i ** 2 + (U_i ** 2 * D)) * K) >> 1) % candidate
  140. V_temp.set(U_i)
  141. V_temp *= U_i
  142. V_temp *= D
  143. V_temp.multiply_accumulate(V_i, V_i)
  144. if V_temp.is_odd():
  145. V_temp += candidate
  146. V_temp >>= 1
  147. V_temp %= candidate
  148. # Multiply
  149. if K.get_bit(i):
  150. # U_i = (((U_temp + V_temp) * K) >> 1) % candidate
  151. U_i.set(U_temp)
  152. U_i += V_temp
  153. if U_i.is_odd():
  154. U_i += candidate
  155. U_i >>= 1
  156. U_i %= candidate
  157. # V_i = (((V_temp + U_temp * D) * K) >> 1) % candidate
  158. V_i.set(V_temp)
  159. V_i.multiply_accumulate(U_temp, D)
  160. if V_i.is_odd():
  161. V_i += candidate
  162. V_i >>= 1
  163. V_i %= candidate
  164. else:
  165. U_i.set(U_temp)
  166. V_i.set(V_temp)
  167. # Step 7
  168. if U_i == 0:
  169. return PROBABLY_PRIME
  170. return COMPOSITE
  171. from Crypto.Util.number import sieve_base as _sieve_base
  172. ## The optimal number of small primes to use for the sieve
  173. ## is probably dependent on the platform and the candidate size
  174. _sieve_base = _sieve_base[:100]
  175. def test_probable_prime(candidate, randfunc=None):
  176. """Test if a number is prime.
  177. A number is qualified as prime if it passes a certain
  178. number of Miller-Rabin tests (dependent on the size
  179. of the number, but such that probability of a false
  180. positive is less than 10^-30) and a single Lucas test.
  181. For instance, a 1024-bit candidate will need to pass
  182. 4 Miller-Rabin tests.
  183. :Parameters:
  184. candidate : integer
  185. The number to test for primality.
  186. randfunc : callable
  187. The routine to draw random bytes from to select Miller-Rabin bases.
  188. :Returns:
  189. ``PROBABLE_PRIME`` if the number if prime with very high probability.
  190. ``COMPOSITE`` if the number is a composite.
  191. For efficiency reasons, ``COMPOSITE`` is also returned for small primes.
  192. """
  193. if randfunc is None:
  194. randfunc = Random.new().read
  195. if not isinstance(candidate, Integer):
  196. candidate = Integer(candidate)
  197. # First, check trial division by the smallest primes
  198. try:
  199. map(candidate.fail_if_divisible_by, _sieve_base)
  200. except ValueError:
  201. return False
  202. # These are the number of Miller-Rabin iterations s.t. p(k, t) < 1E-30,
  203. # with p(k, t) being the probability that a randomly chosen k-bit number
  204. # is composite but still survives t MR iterations.
  205. mr_ranges = ((220, 30), (280, 20), (390, 15), (512, 10),
  206. (620, 7), (740, 6), (890, 5), (1200, 4),
  207. (1700, 3), (3700, 2))
  208. bit_size = candidate.size_in_bits()
  209. try:
  210. mr_iterations = list(filter(lambda x: bit_size < x[0],
  211. mr_ranges))[0][1]
  212. except IndexError:
  213. mr_iterations = 1
  214. if miller_rabin_test(candidate, mr_iterations,
  215. randfunc=randfunc) == COMPOSITE:
  216. return COMPOSITE
  217. if lucas_test(candidate) == COMPOSITE:
  218. return COMPOSITE
  219. return PROBABLY_PRIME
  220. def generate_probable_prime(**kwargs):
  221. """Generate a random probable prime.
  222. The prime will not have any specific properties
  223. (e.g. it will not be a *strong* prime).
  224. Random numbers are evaluated for primality until one
  225. passes all tests, consisting of a certain number of
  226. Miller-Rabin tests with random bases followed by
  227. a single Lucas test.
  228. The number of Miller-Rabin iterations is chosen such that
  229. the probability that the output number is a non-prime is
  230. less than 1E-30 (roughly 2^{-100}).
  231. This approach is compliant to `FIPS PUB 186-4`__.
  232. :Keywords:
  233. exact_bits : integer
  234. The desired size in bits of the probable prime.
  235. It must be at least 160.
  236. randfunc : callable
  237. An RNG function where candidate primes are taken from.
  238. prime_filter : callable
  239. A function that takes an Integer as parameter and returns
  240. True if the number can be passed to further primality tests,
  241. False if it should be immediately discarded.
  242. :Return:
  243. A probable prime in the range 2^exact_bits > p > 2^(exact_bits-1).
  244. .. __: http://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf
  245. """
  246. exact_bits = kwargs.pop("exact_bits", None)
  247. randfunc = kwargs.pop("randfunc", None)
  248. prime_filter = kwargs.pop("prime_filter", lambda x: True)
  249. if kwargs:
  250. print "Unknown parameters:", kwargs.keys()
  251. if exact_bits is None:
  252. raise ValueError("Missing exact_bits parameter")
  253. if exact_bits < 160:
  254. raise ValueError("Prime number is not big enough.")
  255. if randfunc is None:
  256. randfunc = Random.new().read
  257. result = COMPOSITE
  258. while result == COMPOSITE:
  259. candidate = Integer.random(exact_bits=exact_bits,
  260. randfunc=randfunc) | 1
  261. if not prime_filter(candidate):
  262. continue
  263. result = test_probable_prime(candidate, randfunc)
  264. return candidate
  265. def generate_probable_safe_prime(**kwargs):
  266. """Generate a random, probable safe prime.
  267. Note this operation is much slower than generating a simple prime.
  268. :Keywords:
  269. exact_bits : integer
  270. The desired size in bits of the probable safe prime.
  271. randfunc : callable
  272. An RNG function where candidate primes are taken from.
  273. :Return:
  274. A probable safe prime in the range
  275. 2^exact_bits > p > 2^(exact_bits-1).
  276. """
  277. exact_bits = kwargs.pop("exact_bits", None)
  278. randfunc = kwargs.pop("randfunc", None)
  279. if kwargs:
  280. print "Unknown parameters:", kwargs.keys()
  281. if randfunc is None:
  282. randfunc = Random.new().read
  283. result = COMPOSITE
  284. while result == COMPOSITE:
  285. q = generate_probable_prime(exact_bits=exact_bits - 1, randfunc=randfunc)
  286. candidate = q * 2 + 1
  287. if candidate.size_in_bits() != exact_bits:
  288. continue
  289. result = test_probable_prime(candidate, randfunc=randfunc)
  290. return candidate