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/multiprecision/random.hpp>
10#include <boost/multiprecision/integer.hpp>
11
12namespace boost{
13namespace multiprecision{
14namespace detail{
15
16template <class I>
17bool 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
103inline 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
122template <class I>
123typename 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}
128template <class I>
129typename 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
137template <class I, class Engine>
138typename 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(n: 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
201template <class I>
202typename 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
209template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class Engine>
210bool 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
216template <class tag, class Arg1, class Arg2, class Arg3, class Arg4>
217bool 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
227

source code of boost/boost/multiprecision/miller_rabin.hpp