millerrabin.c 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289
  1. /*
  2. * millerrabin.c: Miller-Rabin probabilistic primality testing, as
  3. * declared in sshkeygen.h.
  4. */
  5. #include <assert.h>
  6. #include "ssh.h"
  7. #include "sshkeygen.h"
  8. #include "mpint.h"
  9. #include "mpunsafe.h"
  10. /*
  11. * The Miller-Rabin primality test is an extension to the Fermat
  12. * test. The Fermat test just checks that a^(p-1) == 1 mod p; this
  13. * is vulnerable to Carmichael numbers. Miller-Rabin considers how
  14. * that 1 is derived as well.
  15. *
  16. * Lemma: if a^2 == 1 (mod p), and p is prime, then either a == 1
  17. * or a == -1 (mod p).
  18. *
  19. * Proof: p divides a^2-1, i.e. p divides (a+1)(a-1). Hence,
  20. * since p is prime, either p divides (a+1) or p divides (a-1).
  21. * But this is the same as saying that either a is congruent to
  22. * -1 mod p or a is congruent to +1 mod p. []
  23. *
  24. * Comment: This fails when p is not prime. Consider p=mn, so
  25. * that mn divides (a+1)(a-1). Now we could have m dividing (a+1)
  26. * and n dividing (a-1), without the whole of mn dividing either.
  27. * For example, consider a=10 and p=99. 99 = 9 * 11; 9 divides
  28. * 10-1 and 11 divides 10+1, so a^2 is congruent to 1 mod p
  29. * without a having to be congruent to either 1 or -1.
  30. *
  31. * So the Miller-Rabin test, as well as considering a^(p-1),
  32. * considers a^((p-1)/2), a^((p-1)/4), and so on as far as it can
  33. * go. In other words. we write p-1 as q * 2^k, with k as large as
  34. * possible (i.e. q must be odd), and we consider the powers
  35. *
  36. * a^(q*2^0) a^(q*2^1) ... a^(q*2^(k-1)) a^(q*2^k)
  37. * i.e. a^((n-1)/2^k) a^((n-1)/2^(k-1)) ... a^((n-1)/2) a^(n-1)
  38. *
  39. * If p is to be prime, the last of these must be 1. Therefore, by
  40. * the above lemma, the one before it must be either 1 or -1. And
  41. * _if_ it's 1, then the one before that must be either 1 or -1,
  42. * and so on ... In other words, we expect to see a trailing chain
  43. * of 1s preceded by a -1. (If we're unlucky, our trailing chain of
  44. * 1s will be as long as the list so we'll never get to see what
  45. * lies before it. This doesn't count as a test failure because it
  46. * hasn't _proved_ that p is not prime.)
  47. *
  48. * For example, consider a=2 and p=1729. 1729 is a Carmichael
  49. * number: although it's not prime, it satisfies a^(p-1) == 1 mod p
  50. * for any a coprime to it. So the Fermat test wouldn't have a
  51. * problem with it at all, unless we happened to stumble on an a
  52. * which had a common factor.
  53. *
  54. * So. 1729 - 1 equals 27 * 2^6. So we look at
  55. *
  56. * 2^27 mod 1729 == 645
  57. * 2^108 mod 1729 == 1065
  58. * 2^216 mod 1729 == 1
  59. * 2^432 mod 1729 == 1
  60. * 2^864 mod 1729 == 1
  61. * 2^1728 mod 1729 == 1
  62. *
  63. * We do have a trailing string of 1s, so the Fermat test would
  64. * have been happy. But this trailing string of 1s is preceded by
  65. * 1065; whereas if 1729 were prime, we'd expect to see it preceded
  66. * by -1 (i.e. 1728.). Guards! Seize this impostor.
  67. *
  68. * (If we were unlucky, we might have tried a=16 instead of a=2;
  69. * now 16^27 mod 1729 == 1, so we would have seen a long string of
  70. * 1s and wouldn't have seen the thing _before_ the 1s. So, just
  71. * like the Fermat test, for a given p there may well exist values
  72. * of a which fail to show up its compositeness. So we try several,
  73. * just like the Fermat test. The difference is that Miller-Rabin
  74. * is not _in general_ fooled by Carmichael numbers.)
  75. *
  76. * Put simply, then, the Miller-Rabin test requires us to:
  77. *
  78. * 1. write p-1 as q * 2^k, with q odd
  79. * 2. compute z = (a^q) mod p.
  80. * 3. report success if z == 1 or z == -1.
  81. * 4. square z at most k-1 times, and report success if it becomes
  82. * -1 at any point.
  83. * 5. report failure otherwise.
  84. *
  85. * (We expect z to become -1 after at most k-1 squarings, because
  86. * if it became -1 after k squarings then a^(p-1) would fail to be
  87. * 1. And we don't need to investigate what happens after we see a
  88. * -1, because we _know_ that -1 squared is 1 modulo anything at
  89. * all, so after we've seen a -1 we can be sure of seeing nothing
  90. * but 1s.)
  91. */
  92. struct MillerRabin {
  93. MontyContext *mc;
  94. mp_int *pm1, *m_pm1;
  95. mp_int *lowbit, *two;
  96. };
  97. MillerRabin *miller_rabin_new(mp_int *p)
  98. {
  99. MillerRabin *mr = snew(MillerRabin);
  100. assert(mp_hs_integer(p, 2));
  101. assert(mp_get_bit(p, 0) == 1);
  102. mr->pm1 = mp_copy(p);
  103. mp_sub_integer_into(mr->pm1, mr->pm1, 1);
  104. /*
  105. * Standard bit-twiddling trick for isolating the lowest set bit
  106. * of a number: x & (-x)
  107. */
  108. mr->lowbit = mp_new(mp_max_bits(mr->pm1));
  109. mp_sub_into(mr->lowbit, mr->lowbit, mr->pm1);
  110. mp_and_into(mr->lowbit, mr->lowbit, mr->pm1);
  111. mr->two = mp_from_integer(2);
  112. mr->mc = monty_new(p);
  113. mr->m_pm1 = monty_import(mr->mc, mr->pm1);
  114. return mr;
  115. }
  116. void miller_rabin_free(MillerRabin *mr)
  117. {
  118. mp_free(mr->pm1);
  119. mp_free(mr->m_pm1);
  120. mp_free(mr->lowbit);
  121. mp_free(mr->two);
  122. monty_free(mr->mc);
  123. smemclr(mr, sizeof(*mr));
  124. sfree(mr);
  125. }
  126. /*
  127. * The main internal function that implements a single M-R test.
  128. *
  129. * Expects the witness integer to be in Montgomery representation.
  130. * (Since in live use witnesses are invented at random, this imposes
  131. * no extra cost on the callers, and saves effort in here.)
  132. */
  133. static struct mr_result miller_rabin_test_inner(MillerRabin *mr, mp_int *mw)
  134. {
  135. mp_int *acc = mp_copy(monty_identity(mr->mc));
  136. mp_int *spare = mp_new(mp_max_bits(mr->pm1));
  137. size_t bit = mp_max_bits(mr->pm1);
  138. /*
  139. * The obvious approach to Miller-Rabin would be to start by
  140. * calling monty_pow to raise w to the power q, and then square it
  141. * k times ourselves. But that introduces a timing leak that gives
  142. * away the value of k, i.e., how many factors of 2 there are in
  143. * p-1.
  144. *
  145. * Instead, we don't call monty_pow at all. We do a modular
  146. * exponentiation ourselves to compute w^((p-1)/2), using the
  147. * technique that works from the top bit of the exponent
  148. * downwards. That is, in each iteration we compute
  149. * w^floor(exponent/2^i) for i one less than the previous
  150. * iteration, by squaring the value we previously had and then
  151. * optionally multiplying in w if the next exponent bit is 1.
  152. *
  153. * At the end of that process, once i <= k, the division
  154. * (exponent/2^i) yields an integer, so the values we're computing
  155. * are not just w^(floor of that), but w^(exactly that). In other
  156. * words, the last k intermediate values of this modexp are
  157. * precisely the values M-R wants to check against +1 or -1.
  158. *
  159. * So we interleave those checks with the modexp loop itself, and
  160. * to avoid a timing leak, we check _every_ intermediate result
  161. * against (the Montgomery representations of) both +1 and -1. And
  162. * then we do bitwise masking to arrange that only the sensible
  163. * ones of those checks find their way into our final answer.
  164. */
  165. unsigned active = 0;
  166. struct mr_result result;
  167. result.passed = result.potential_primitive_root = 0;
  168. while (bit-- > 1) {
  169. /*
  170. * In this iteration, we're computing w^(2e) or w^(2e+1),
  171. * where we have w^e from the previous iteration. So we square
  172. * the value we had already, and then optionally multiply in
  173. * another copy of w depending on the next bit of the exponent.
  174. */
  175. monty_mul_into(mr->mc, acc, acc, acc);
  176. monty_mul_into(mr->mc, spare, acc, mw);
  177. mp_select_into(acc, acc, spare, mp_get_bit(mr->pm1, bit));
  178. /*
  179. * mr->lowbit is a number with only one bit set, corresponding
  180. * to the lowest set bit in p-1. So when that's the bit of the
  181. * exponent we've just processed, we'll detect it by setting
  182. * first_iter to true. That's our indication that we're now
  183. * generating intermediate results useful to M-R, so we also
  184. * set 'active', which stays set from then on.
  185. */
  186. unsigned first_iter = mp_get_bit(mr->lowbit, bit);
  187. active |= first_iter;
  188. /*
  189. * Check the intermediate result against both +1 and -1.
  190. */
  191. unsigned is_plus_1 = mp_cmp_eq(acc, monty_identity(mr->mc));
  192. unsigned is_minus_1 = mp_cmp_eq(acc, mr->m_pm1);
  193. /*
  194. * M-R must report success iff either: the first of the useful
  195. * intermediate results (which is w^q) is 1, or _any_ of them
  196. * (from w^q all the way up to w^((p-1)/2)) is -1.
  197. *
  198. * So we want to pass the test if is_plus_1 is set on the
  199. * first iteration, or if is_minus_1 is set on any iteration.
  200. */
  201. result.passed |= (first_iter & is_plus_1);
  202. result.passed |= (active & is_minus_1);
  203. /*
  204. * In the final iteration, is_minus_1 is also used to set the
  205. * 'potential primitive root' flag, because we haven't found
  206. * any exponent smaller than p-1 for which w^(that) == 1.
  207. */
  208. if (bit == 1)
  209. result.potential_primitive_root = is_minus_1;
  210. }
  211. mp_free(acc);
  212. mp_free(spare);
  213. return result;
  214. }
  215. /*
  216. * Wrapper on miller_rabin_test_inner for the convenience of
  217. * testcrypt. Expects the witness integer to be literal, so we
  218. * monty_import it before running the real test.
  219. */
  220. struct mr_result miller_rabin_test(MillerRabin *mr, mp_int *w)
  221. {
  222. mp_int *mw = monty_import(mr->mc, w);
  223. struct mr_result result = miller_rabin_test_inner(mr, mw);
  224. mp_free(mw);
  225. return result;
  226. }
  227. bool miller_rabin_test_random(MillerRabin *mr)
  228. {
  229. mp_int *mw = mp_random_in_range(mr->two, mr->pm1);
  230. struct mr_result result = miller_rabin_test_inner(mr, mw);
  231. mp_free(mw);
  232. return result.passed;
  233. }
  234. mp_int *miller_rabin_find_potential_primitive_root(MillerRabin *mr)
  235. {
  236. while (true) {
  237. mp_int *mw = mp_unsafe_shrink(mp_random_in_range(mr->two, mr->pm1));
  238. struct mr_result result = miller_rabin_test_inner(mr, mw);
  239. if (result.passed && result.potential_primitive_root) {
  240. mp_int *pr = monty_export(mr->mc, mw);
  241. mp_free(mw);
  242. return pr;
  243. }
  244. mp_free(mw);
  245. if (!result.passed) {
  246. return NULL;
  247. }
  248. }
  249. }
  250. unsigned miller_rabin_checks_needed(unsigned bits)
  251. {
  252. /* Table 4.4 from Handbook of Applied Cryptography */
  253. return (bits >= 1300 ? 2 : bits >= 850 ? 3 : bits >= 650 ? 4 :
  254. bits >= 550 ? 5 : bits >= 450 ? 6 : bits >= 400 ? 7 :
  255. bits >= 350 ? 8 : bits >= 300 ? 9 : bits >= 250 ? 12 :
  256. bits >= 200 ? 15 : bits >= 150 ? 18 : 27);
  257. }