boost/multiprecision/miller_rabin.hpp

0001 ///////////////////////////////////////////////////////////////
0002 //  Copyright 2012 John Maddock. Distributed under the Boost
0003 //  Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
0005
0006 #ifndef BOOST_MP_MR_HPP
0007 #define BOOST_MP_MR_HPP
0008
0009 #include <random>
0010 #include <cstdint>
0011 #include <type_traits>
0012 #include <boost/multiprecision/detail/standalone_config.hpp>
0013 #include <boost/multiprecision/integer.hpp>
0014 #include <boost/multiprecision/detail/uniform_int_distribution.hpp>
0015 #include <boost/multiprecision/detail/assert.hpp>
0016
0017 namespace boost {
0018 namespace multiprecision {
0019 namespace detail {
0020
0021 template <class I>
0022 bool check_small_factors(const I& n)
0023 {
0024    constexpr std::uint32_t small_factors1[] = {
0025        3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u};
0026    constexpr std::uint32_t pp1 = 223092870u;
0027
0028    std::uint32_t m1 = integer_modulus(n, pp1);
0029
0030    for (std::size_t i = 0; i < sizeof(small_factors1) / sizeof(small_factors1[0]); ++i)
0031    {
0032       BOOST_MP_ASSERT(pp1 % small_factors1[i] == 0);
0033       if (m1 % small_factors1[i] == 0)
0034          return false;
0035    }
0036
0037    constexpr std::uint32_t small_factors2[] = {
0038        29u, 31u, 37u, 41u, 43u, 47u};
0039    constexpr std::uint32_t pp2 = 2756205443u;
0040
0041    m1 = integer_modulus(n, pp2);
0042
0043    for (std::size_t i = 0; i < sizeof(small_factors2) / sizeof(small_factors2[0]); ++i)
0044    {
0045       BOOST_MP_ASSERT(pp2 % small_factors2[i] == 0);
0046       if (m1 % small_factors2[i] == 0)
0047          return false;
0048    }
0049
0050    constexpr std::uint32_t small_factors3[] = {
0051        53u, 59u, 61u, 67u, 71u};
0052    constexpr std::uint32_t pp3 = 907383479u;
0053
0054    m1 = integer_modulus(n, pp3);
0055
0056    for (std::size_t i = 0; i < sizeof(small_factors3) / sizeof(small_factors3[0]); ++i)
0057    {
0058       BOOST_MP_ASSERT(pp3 % small_factors3[i] == 0);
0059       if (m1 % small_factors3[i] == 0)
0060          return false;
0061    }
0062
0063    constexpr std::uint32_t small_factors4[] = {
0064        73u, 79u, 83u, 89u, 97u};
0065    constexpr std::uint32_t pp4 = 4132280413u;
0066
0067    m1 = integer_modulus(n, pp4);
0068
0069    for (std::size_t i = 0; i < sizeof(small_factors4) / sizeof(small_factors4[0]); ++i)
0070    {
0071       BOOST_MP_ASSERT(pp4 % small_factors4[i] == 0);
0072       if (m1 % small_factors4[i] == 0)
0073          return false;
0074    }
0075
0076    constexpr std::uint32_t small_factors5[6][4] = {
0077        {101u, 103u, 107u, 109u},
0078        {113u, 127u, 131u, 137u},
0079        {139u, 149u, 151u, 157u},
0080        {163u, 167u, 173u, 179u},
0081        {181u, 191u, 193u, 197u},
0082        {199u, 211u, 223u, 227u}};
0083    constexpr std::uint32_t pp5[6] =
0084        {
0085            121330189u,
0086            113u * 127u * 131u * 137u,
0087            139u * 149u * 151u * 157u,
0088            163u * 167u * 173u * 179u,
0089            181u * 191u * 193u * 197u,
0090            199u * 211u * 223u * 227u};
0091
0092    for (std::size_t k = 0; k < sizeof(pp5) / sizeof(*pp5); ++k)
0093    {
0094       m1 = integer_modulus(n, pp5[k]);
0095
0096       for (std::size_t i = 0; i < 4; ++i)
0097       {
0098          BOOST_MP_ASSERT(pp5[k] % small_factors5[k][i] == 0);
0099          if (m1 % small_factors5[k][i] == 0)
0100             return false;
0101       }
0102    }
0103    return true;
0104 }
0105
0106 inline bool is_small_prime(std::size_t n)
0107 {
0108    constexpr unsigned char p[] =
0109        {
0110            3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u, 29u, 31u,
0111            37u, 41u, 43u, 47u, 53u, 59u, 61u, 67u, 71u, 73u,
0112            79u, 83u, 89u, 97u, 101u, 103u, 107u, 109u, 113u,
0113            127u, 131u, 137u, 139u, 149u, 151u, 157u, 163u,
0114            167u, 173u, 179u, 181u, 191u, 193u, 197u, 199u,
0115            211u, 223u, 227u};
0116    for (std::size_t i = 0; i < sizeof(p) / sizeof(*p); ++i)
0117    {
0118       if (n == p[i])
0119          return true;
0120    }
0121    return false;
0122 }
0123
0124 template <class I>
0125 typename std::enable_if<std::is_convertible<I, unsigned>::value, unsigned>::type
0126 cast_to_unsigned(const I& val)
0127 {
0128    return static_cast<unsigned>(val);
0129 }
0130 template <class I>
0131 typename std::enable_if<!std::is_convertible<I, unsigned>::value, unsigned>::type
0132 cast_to_unsigned(const I& val)
0133 {
0134    return val.template convert_to<unsigned>();
0135 }
0136
0137 } // namespace detail
0138
0139 template <class I, class Engine>
0140 typename std::enable_if<number_category<I>::value == number_kind_integer, bool>::type
0141 miller_rabin_test(const I& n, std::size_t trials, Engine& gen)
0142 {
0143    using number_type = I;
0144
0145    if (n == 2)
0146       return true; // Trivial special case.
0147    if (bit_test(n, 0) == 0)
0148       return false; // n is even
0149    if (n <= 227)
0150       return detail::is_small_prime(detail::cast_to_unsigned(n));
0151
0152    if (!detail::check_small_factors(n))
0153       return false;
0154
0155    number_type nm1 = n - 1;
0156    //
0157    // Begin with a single Fermat test - it excludes a lot of candidates:
0158    //
0159    number_type q(228), x, y; // We know n is greater than this, as we've excluded small factors
0160    x = powm(q, nm1, n);
0161    if (x != 1u)
0162       return false;
0163
0164    q          = n - 1;
0165    std::size_t k = lsb(q);
0166    q >>= k;
0167
0168    // Declare our random number generator:
0169    boost::multiprecision::uniform_int_distribution<number_type> dist(2, n - 2);
0170
0171    //
0172    // Execute the trials:
0173    //
0174    for (std::size_t i = 0; i < trials; ++i)
0175    {
0176       x          = dist(gen);
0177       y          = powm(x, q, n);
0178       std::size_t j = 0;
0179       while (true)
0180       {
0181          if (y == nm1)
0182             break;
0183          if (y == 1)
0184          {
0185             if (j == 0)
0186                break;
0187             return false; // test failed
0188          }
0189          if (++j == k)
0190             return false; // failed
0191          y = powm(y, 2, n);
0192       }
0193    }
0194    return true; // Yeheh! probably prime.
0195 }
0196
0197 template <class I>
0198 typename std::enable_if<number_category<I>::value == number_kind_integer, bool>::type
0199 miller_rabin_test(const I& x, std::size_t trials)
0200 {
0201    static std::mt19937 gen;
0202    return miller_rabin_test(x, trials, gen);
0203 }
0204
0205 template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class Engine>
0206 bool miller_rabin_test(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& n, std::size_t trials, Engine& gen)
0207 {
0208    using number_type = typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type;
0209    return miller_rabin_test(number_type(n), trials, gen);
0210 }
0211
0212 template <class tag, class Arg1, class Arg2, class Arg3, class Arg4>
0213 bool miller_rabin_test(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& n, std::size_t trials)
0214 {
0215    using number_type = typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type;
0216    return miller_rabin_test(number_type(n), trials);
0217 }
0218
0219 }} // namespace boost::multiprecision
0220
0221 #endif