miller_rabin.hpp 6.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215
  1. ///////////////////////////////////////////////////////////////
  2. // Copyright 2012 John Maddock. Distributed under the Boost
  3. // Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
  5. #ifndef BOOST_MP_MR_HPP
  6. #define BOOST_MP_MR_HPP
  7. #include <boost/random.hpp>
  8. #include <boost/multiprecision/integer.hpp>
  9. namespace boost {
  10. namespace multiprecision {
  11. namespace detail {
  12. template <class I>
  13. bool check_small_factors(const I& n)
  14. {
  15. constexpr const std::uint32_t small_factors1[] = {
  16. 3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u};
  17. constexpr const std::uint32_t pp1 = 223092870u;
  18. std::uint32_t m1 = integer_modulus(n, pp1);
  19. for (unsigned i = 0; i < sizeof(small_factors1) / sizeof(small_factors1[0]); ++i)
  20. {
  21. BOOST_ASSERT(pp1 % small_factors1[i] == 0);
  22. if (m1 % small_factors1[i] == 0)
  23. return false;
  24. }
  25. constexpr const std::uint32_t small_factors2[] = {
  26. 29u, 31u, 37u, 41u, 43u, 47u};
  27. constexpr const std::uint32_t pp2 = 2756205443u;
  28. m1 = integer_modulus(n, pp2);
  29. for (unsigned i = 0; i < sizeof(small_factors2) / sizeof(small_factors2[0]); ++i)
  30. {
  31. BOOST_ASSERT(pp2 % small_factors2[i] == 0);
  32. if (m1 % small_factors2[i] == 0)
  33. return false;
  34. }
  35. constexpr const std::uint32_t small_factors3[] = {
  36. 53u, 59u, 61u, 67u, 71u};
  37. constexpr const std::uint32_t pp3 = 907383479u;
  38. m1 = integer_modulus(n, pp3);
  39. for (unsigned i = 0; i < sizeof(small_factors3) / sizeof(small_factors3[0]); ++i)
  40. {
  41. BOOST_ASSERT(pp3 % small_factors3[i] == 0);
  42. if (m1 % small_factors3[i] == 0)
  43. return false;
  44. }
  45. constexpr const std::uint32_t small_factors4[] = {
  46. 73u, 79u, 83u, 89u, 97u};
  47. constexpr const std::uint32_t pp4 = 4132280413u;
  48. m1 = integer_modulus(n, pp4);
  49. for (unsigned i = 0; i < sizeof(small_factors4) / sizeof(small_factors4[0]); ++i)
  50. {
  51. BOOST_ASSERT(pp4 % small_factors4[i] == 0);
  52. if (m1 % small_factors4[i] == 0)
  53. return false;
  54. }
  55. constexpr const std::uint32_t small_factors5[6][4] = {
  56. {101u, 103u, 107u, 109u},
  57. {113u, 127u, 131u, 137u},
  58. {139u, 149u, 151u, 157u},
  59. {163u, 167u, 173u, 179u},
  60. {181u, 191u, 193u, 197u},
  61. {199u, 211u, 223u, 227u}};
  62. constexpr const std::uint32_t pp5[6] =
  63. {
  64. 121330189u,
  65. 113u * 127u * 131u * 137u,
  66. 139u * 149u * 151u * 157u,
  67. 163u * 167u * 173u * 179u,
  68. 181u * 191u * 193u * 197u,
  69. 199u * 211u * 223u * 227u};
  70. for (unsigned k = 0; k < sizeof(pp5) / sizeof(*pp5); ++k)
  71. {
  72. m1 = integer_modulus(n, pp5[k]);
  73. for (unsigned i = 0; i < 4; ++i)
  74. {
  75. BOOST_ASSERT(pp5[k] % small_factors5[k][i] == 0);
  76. if (m1 % small_factors5[k][i] == 0)
  77. return false;
  78. }
  79. }
  80. return true;
  81. }
  82. inline bool is_small_prime(unsigned n)
  83. {
  84. constexpr const unsigned char p[] =
  85. {
  86. 3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u, 29u, 31u,
  87. 37u, 41u, 43u, 47u, 53u, 59u, 61u, 67u, 71u, 73u,
  88. 79u, 83u, 89u, 97u, 101u, 103u, 107u, 109u, 113u,
  89. 127u, 131u, 137u, 139u, 149u, 151u, 157u, 163u,
  90. 167u, 173u, 179u, 181u, 191u, 193u, 197u, 199u,
  91. 211u, 223u, 227u};
  92. for (unsigned i = 0; i < sizeof(p) / sizeof(*p); ++i)
  93. {
  94. if (n == p[i])
  95. return true;
  96. }
  97. return false;
  98. }
  99. template <class I>
  100. typename std::enable_if<std::is_convertible<I, unsigned>::value, unsigned>::type
  101. cast_to_unsigned(const I& val)
  102. {
  103. return static_cast<unsigned>(val);
  104. }
  105. template <class I>
  106. typename std::enable_if<!std::is_convertible<I, unsigned>::value, unsigned>::type
  107. cast_to_unsigned(const I& val)
  108. {
  109. return val.template convert_to<unsigned>();
  110. }
  111. } // namespace detail
  112. template <class I, class Engine>
  113. typename std::enable_if<number_category<I>::value == number_kind_integer, bool>::type
  114. miller_rabin_test(const I& n, unsigned trials, Engine& gen)
  115. {
  116. using number_type = I;
  117. if (n == 2)
  118. return true; // Trivial special case.
  119. if (bit_test(n, 0) == 0)
  120. return false; // n is even
  121. if (n <= 227)
  122. return detail::is_small_prime(detail::cast_to_unsigned(n));
  123. if (!detail::check_small_factors(n))
  124. return false;
  125. number_type nm1 = n - 1;
  126. //
  127. // Begin with a single Fermat test - it excludes a lot of candidates:
  128. //
  129. number_type q(228), x, y; // We know n is greater than this, as we've excluded small factors
  130. x = powm(q, nm1, n);
  131. if (x != 1u)
  132. return false;
  133. q = n - 1;
  134. unsigned k = lsb(q);
  135. q >>= k;
  136. // Declare our random number generator:
  137. boost::random::uniform_int_distribution<number_type> dist(2, n - 2);
  138. //
  139. // Execute the trials:
  140. //
  141. for (unsigned i = 0; i < trials; ++i)
  142. {
  143. x = dist(gen);
  144. y = powm(x, q, n);
  145. unsigned j = 0;
  146. while (true)
  147. {
  148. if (y == nm1)
  149. break;
  150. if (y == 1)
  151. {
  152. if (j == 0)
  153. break;
  154. return false; // test failed
  155. }
  156. if (++j == k)
  157. return false; // failed
  158. y = powm(y, 2, n);
  159. }
  160. }
  161. return true; // Yeheh! probably prime.
  162. }
  163. template <class I>
  164. typename std::enable_if<number_category<I>::value == number_kind_integer, bool>::type
  165. miller_rabin_test(const I& x, unsigned trials)
  166. {
  167. static mt19937 gen;
  168. return miller_rabin_test(x, trials, gen);
  169. }
  170. template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class Engine>
  171. bool miller_rabin_test(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& n, unsigned trials, Engine& gen)
  172. {
  173. using number_type = typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type;
  174. return miller_rabin_test(number_type(n), trials, gen);
  175. }
  176. template <class tag, class Arg1, class Arg2, class Arg3, class Arg4>
  177. bool miller_rabin_test(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& n, unsigned trials)
  178. {
  179. using number_type = typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type;
  180. return miller_rabin_test(number_type(n), trials);
  181. }
  182. }} // namespace boost::multiprecision
  183. #endif