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