ceph/src/boost/boost/multiprecision/miller_rabin.hpp

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