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