comparison DEPENDENCIES/generic/include/boost/multiprecision/miller_rabin.hpp @ 16:2665513ce2d3

Add boost headers

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