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_
6 #ifndef BOOST_MP_MR_HPP
7 #define BOOST_MP_MR_HPP
9 #include <boost/random.hpp>
10 #include <boost/multiprecision/integer.hpp>
13 namespace multiprecision{
17 bool check_small_factors(const I& n)
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;
23 boost::uint32_t m1 = integer_modulus(n, pp1);
25 for(unsigned i = 0; i < sizeof(small_factors1) / sizeof(small_factors1[0]); ++i)
27 BOOST_ASSERT(pp1 % small_factors1[i] == 0);
28 if(m1 % small_factors1[i] == 0)
32 static const boost::uint32_t small_factors2[] = {
33 29u, 31u, 37u, 41u, 43u, 47u };
34 static const boost::uint32_t pp2 = 2756205443u;
36 m1 = integer_modulus(n, pp2);
38 for(unsigned i = 0; i < sizeof(small_factors2) / sizeof(small_factors2[0]); ++i)
40 BOOST_ASSERT(pp2 % small_factors2[i] == 0);
41 if(m1 % small_factors2[i] == 0)
45 static const boost::uint32_t small_factors3[] = {
46 53u, 59u, 61u, 67u, 71u };
47 static const boost::uint32_t pp3 = 907383479u;
49 m1 = integer_modulus(n, pp3);
51 for(unsigned i = 0; i < sizeof(small_factors3) / sizeof(small_factors3[0]); ++i)
53 BOOST_ASSERT(pp3 % small_factors3[i] == 0);
54 if(m1 % small_factors3[i] == 0)
58 static const boost::uint32_t small_factors4[] = {
59 73u, 79u, 83u, 89u, 97u };
60 static const boost::uint32_t pp4 = 4132280413u;
62 m1 = integer_modulus(n, pp4);
64 for(unsigned i = 0; i < sizeof(small_factors4) / sizeof(small_factors4[0]); ++i)
66 BOOST_ASSERT(pp4 % small_factors4[i] == 0);
67 if(m1 % small_factors4[i] == 0)
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 }
79 static const boost::uint32_t pp5[6] =
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
89 for(unsigned k = 0; k < sizeof(pp5) / sizeof(*pp5); ++k)
91 m1 = integer_modulus(n, pp5[k]);
93 for(unsigned i = 0; i < 4; ++i)
95 BOOST_ASSERT(pp5[k] % small_factors5[k][i] == 0);
96 if(m1 % small_factors5[k][i] == 0)
103 inline bool is_small_prime(unsigned n)
105 static const unsigned char p[] =
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,
114 for(unsigned i = 0; i < sizeof(p) / sizeof(*p); ++i)
123 typename enable_if_c<is_convertible<I, unsigned>::value, unsigned>::type
124 cast_to_unsigned(const I& val)
126 return static_cast<unsigned>(val);
129 typename disable_if_c<is_convertible<I, unsigned>::value, unsigned>::type
130 cast_to_unsigned(const I& val)
132 return val.template convert_to<unsigned>();
135 } // namespace detail
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)
142 #pragma warning(push)
143 #pragma warning(disable:4127)
145 typedef I number_type;
148 return true; // Trivial special case.
149 if(bit_test(n, 0) == 0)
150 return false; // n is even
152 return detail::is_small_prime(detail::cast_to_unsigned(n));
154 if(!detail::check_small_factors(n))
157 number_type nm1 = n - 1;
159 // Begin with a single Fermat test - it excludes a lot of candidates:
161 number_type q(228), x, y; // We know n is greater than this, as we've excluded small factors
170 // Declare our random number generator:
171 boost::random::uniform_int_distribution<number_type> dist(2, n - 2);
173 // Execute the trials:
175 for(unsigned i = 0; i < trials; ++i)
188 return false; // test failed
191 return false; // failed
195 return true; // Yeheh! probably prime.
202 typename enable_if_c<number_category<I>::value == number_kind_integer, bool>::type
203 miller_rabin_test(const I& x, unsigned trials)
206 return miller_rabin_test(x, trials, gen);
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)
212 typedef typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type number_type;
213 return miller_rabin_test(number_type(n), trials, gen);
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)
219 typedef typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type number_type;
220 return miller_rabin_test(number_type(n), trials);