Chris@16: ///////////////////////////////////////////////////////////////
Chris@16: //  Copyright 2012 John Maddock. Distributed under the Boost
Chris@16: //  Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
Chris@16: 
Chris@16: #ifndef BOOST_MP_MR_HPP
Chris@16: #define BOOST_MP_MR_HPP
Chris@16: 
Chris@16: #include <boost/multiprecision/random.hpp>
Chris@16: #include <boost/multiprecision/integer.hpp>
Chris@16: 
Chris@16: namespace boost{
Chris@16: namespace multiprecision{
Chris@16: namespace detail{
Chris@16: 
Chris@16: template <class I>
Chris@16: bool check_small_factors(const I& n)
Chris@16: {
Chris@16:    static const boost::uint32_t small_factors1[] = {
Chris@16:       3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u };
Chris@16:    static const boost::uint32_t pp1 = 223092870u;
Chris@16: 
Chris@16:    boost::uint32_t m1 = integer_modulus(n, pp1);
Chris@16: 
Chris@16:    for(unsigned i = 0; i < sizeof(small_factors1) / sizeof(small_factors1[0]); ++i)
Chris@16:    {
Chris@16:       BOOST_ASSERT(pp1 % small_factors1[i] == 0);
Chris@16:       if(m1 % small_factors1[i] == 0)
Chris@16:          return false;
Chris@16:    }
Chris@16: 
Chris@16:    static const boost::uint32_t small_factors2[] = {
Chris@16:       29u, 31u, 37u, 41u, 43u, 47u };
Chris@16:    static const boost::uint32_t pp2 = 2756205443u;
Chris@16: 
Chris@16:    m1 = integer_modulus(n, pp2);
Chris@16: 
Chris@16:    for(unsigned i = 0; i < sizeof(small_factors2) / sizeof(small_factors2[0]); ++i)
Chris@16:    {
Chris@16:       BOOST_ASSERT(pp2 % small_factors2[i] == 0);
Chris@16:       if(m1 % small_factors2[i] == 0)
Chris@16:          return false;
Chris@16:    }
Chris@16: 
Chris@16:    static const boost::uint32_t small_factors3[] = {
Chris@16:       53u, 59u, 61u, 67u, 71u };
Chris@16:    static const boost::uint32_t pp3 = 907383479u;
Chris@16: 
Chris@16:    m1 = integer_modulus(n, pp3);
Chris@16: 
Chris@16:    for(unsigned i = 0; i < sizeof(small_factors3) / sizeof(small_factors3[0]); ++i)
Chris@16:    {
Chris@16:       BOOST_ASSERT(pp3 % small_factors3[i] == 0);
Chris@16:       if(m1 % small_factors3[i] == 0)
Chris@16:          return false;
Chris@16:    }
Chris@16: 
Chris@16:    static const boost::uint32_t small_factors4[] = {
Chris@16:       73u, 79u, 83u, 89u, 97u };
Chris@16:    static const boost::uint32_t pp4 = 4132280413u;
Chris@16: 
Chris@16:    m1 = integer_modulus(n, pp4);
Chris@16: 
Chris@16:    for(unsigned i = 0; i < sizeof(small_factors4) / sizeof(small_factors4[0]); ++i)
Chris@16:    {
Chris@16:       BOOST_ASSERT(pp4 % small_factors4[i] == 0);
Chris@16:       if(m1 % small_factors4[i] == 0)
Chris@16:          return false;
Chris@16:    }
Chris@16: 
Chris@16:    static const boost::uint32_t small_factors5[6][4] = {
Chris@16:       { 101u, 103u, 107u, 109u },
Chris@16:       { 113u, 127u, 131u, 137u },
Chris@16:       { 139u, 149u, 151u, 157u },
Chris@16:       { 163u, 167u, 173u, 179u },
Chris@16:       { 181u, 191u, 193u, 197u },
Chris@16:       { 199u, 211u, 223u, 227u }
Chris@16:    };
Chris@16:    static const boost::uint32_t pp5[6] = 
Chris@16:    { 
Chris@16:       121330189u, 
Chris@16:       113u * 127u * 131u * 137u, 
Chris@16:       139u * 149u * 151u * 157u,
Chris@16:       163u * 167u * 173u * 179u,
Chris@16:       181u * 191u * 193u * 197u,
Chris@16:       199u * 211u * 223u * 227u
Chris@16:    };
Chris@16: 
Chris@16:    for(unsigned k = 0; k < sizeof(pp5) / sizeof(*pp5); ++k)
Chris@16:    {
Chris@16:       m1 = integer_modulus(n, pp5[k]);
Chris@16: 
Chris@16:       for(unsigned i = 0; i < 4; ++i)
Chris@16:       {
Chris@16:          BOOST_ASSERT(pp5[k] % small_factors5[k][i] == 0);
Chris@16:          if(m1 % small_factors5[k][i] == 0)
Chris@16:             return false;
Chris@16:       }
Chris@16:    }
Chris@16:    return true;
Chris@16: }
Chris@16: 
Chris@16: inline bool is_small_prime(unsigned n)
Chris@16: {
Chris@16:    static const unsigned char p[] = 
Chris@16:    {
Chris@16:       3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u, 29u, 31u, 
Chris@16:       37u, 41u, 43u, 47u, 53u, 59u, 61u, 67u, 71u, 73u, 
Chris@16:       79u, 83u, 89u, 97u, 101u, 103u, 107u, 109u, 113u, 
Chris@16:       127u, 131u, 137u, 139u, 149u, 151u, 157u, 163u, 
Chris@16:       167u, 173u, 179u, 181u, 191u, 193u, 197u, 199u, 
Chris@16:       211u, 223u, 227u
Chris@16:    };
Chris@16:    for(unsigned i = 0; i < sizeof(p) / sizeof(*p); ++i)
Chris@16:    {
Chris@16:       if(n == p[i])
Chris@16:          return true;
Chris@16:    }
Chris@16:    return false;
Chris@16: }
Chris@16: 
Chris@16: template <class I>
Chris@16: typename enable_if_c<is_convertible<I, unsigned>::value, unsigned>::type
Chris@16:    cast_to_unsigned(const I& val)
Chris@16: {
Chris@16:    return static_cast<unsigned>(val);
Chris@16: }
Chris@16: template <class I>
Chris@16: typename disable_if_c<is_convertible<I, unsigned>::value, unsigned>::type
Chris@16:    cast_to_unsigned(const I& val)
Chris@16: {
Chris@16:    return val.template convert_to<unsigned>();
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class I, class Engine>
Chris@16: typename enable_if_c<number_category<I>::value == number_kind_integer, bool>::type 
Chris@16:    miller_rabin_test(const I& n, unsigned trials, Engine& gen)
Chris@16: {
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(push)
Chris@16: #pragma warning(disable:4127)
Chris@16: #endif
Chris@16:    typedef I number_type;
Chris@16: 
Chris@16:    if(bit_test(n, 0) == 0)
Chris@16:       return false;  // n is even
Chris@16:    if(n <= 227)
Chris@16:       return detail::is_small_prime(detail::cast_to_unsigned(n));
Chris@16: 
Chris@16:    if(!detail::check_small_factors(n))
Chris@16:       return false;
Chris@16: 
Chris@16:    number_type nm1 = n - 1;
Chris@16:    //
Chris@16:    // Begin with a single Fermat test - it excludes a lot of candidates:
Chris@16:    //
Chris@16:    number_type q(228), x, y; // We know n is greater than this, as we've excluded small factors
Chris@16:    x = powm(q, nm1, n);
Chris@16:    if(x != 1u)
Chris@16:       return false;
Chris@16: 
Chris@16:    q = n - 1;
Chris@16:    unsigned k = lsb(q);
Chris@16:    q >>= k;
Chris@16: 
Chris@16:    // Declare our random number generator:
Chris@16:    boost::random::uniform_int_distribution<number_type> dist(2, n - 2);
Chris@16:    //
Chris@16:    // Execute the trials:
Chris@16:    //
Chris@16:    for(unsigned i = 0; i < trials; ++i)
Chris@16:    {
Chris@16:       x = dist(gen);
Chris@16:       y = powm(x, q, n);
Chris@16:       unsigned j = 0;
Chris@16:       while(true)
Chris@16:       {
Chris@16:          if(y == nm1)
Chris@16:             break;
Chris@16:          if(y == 1)
Chris@16:          {
Chris@16:             if(j == 0)
Chris@16:                break;
Chris@16:             return false; // test failed
Chris@16:          }
Chris@16:          if(++j == k)
Chris@16:             return false; // failed
Chris@16:          y = powm(y, 2, n);
Chris@16:       }
Chris@16:    }
Chris@16:    return true;  // Yeheh! probably prime.
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(pop)
Chris@16: #endif
Chris@16: }
Chris@16: 
Chris@16: template <class I>
Chris@16: typename enable_if_c<number_category<I>::value == number_kind_integer, bool>::type 
Chris@16:    miller_rabin_test(const I& x, unsigned trials)
Chris@16: {
Chris@16:    static mt19937 gen;
Chris@16:    return miller_rabin_test(x, trials, gen);
Chris@16: }
Chris@16: 
Chris@16: template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class Engine>
Chris@16: bool miller_rabin_test(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4> & n, unsigned trials, Engine& gen)
Chris@16: {
Chris@16:    typedef typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type number_type;
Chris@16:    return miller_rabin_test(number_type(n), trials, gen);
Chris@16: }
Chris@16: 
Chris@16: template <class tag, class Arg1, class Arg2, class Arg3, class Arg4>
Chris@16: bool miller_rabin_test(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4> & n, unsigned trials)
Chris@16: {
Chris@16:    typedef typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type number_type;
Chris@16:    return miller_rabin_test(number_type(n), trials);
Chris@16: }
Chris@16: 
Chris@16: }} // namespaces
Chris@16: 
Chris@16: #endif
Chris@16: