Effective implementation of algorithms (Master Thesis)
Effective and error-free implementation of algorithms
|
00001 #ifndef H_MATH_PRIMES_PRIMES_FAST 00002 #define H_MATH_PRIMES_PRIMES_FAST 00003 00006 #include "utils/preconditions/preconditions.h" 00007 #include "math/powermod/powermod.h" 00008 #include "utils/static_assert/static_assert.h" 00009 #include <limits> 00010 00011 namespace math { 00012 namespace primes { 00013 00031 template <class PowerModImpl> 00032 class PrimesFast_ { 00034 typedef long long int BaseType; 00035 // check base type 00036 STATIC_ASSERT(std::numeric_limits<BaseType>::digits > 50, 00037 "We need at least 50-bit integer as a base tyep"); 00038 00039 public: 00048 static bool isPrime(BaseType p) { 00049 // theoretical limit of the test based on original article 00050 const BaseType TEST_MAX = 341550071728321LL; 00051 const BaseType witnesses[] = {2, 3, 5, 7, 11, 13, 17, 0}; 00052 00053 Preconditions::check(p > 0, "tested integer should be positive!"); 00054 Preconditions::check(p < TEST_MAX, 00055 "Sorry, but we support only calculations up to 3*10^14"); 00056 00057 // small trivial cases 00058 if (p == 2) { 00059 return true; 00060 } else if (p == 1 || p % 2 == 0) { 00061 return false; 00062 } 00063 00064 // p is odd, now we want to write p-1 = 2^s * d 00065 BaseType s = 0; 00066 BaseType d = p - 1; 00067 while (d % 2 == 0) { 00068 s++; 00069 d /= 2; 00070 } 00071 // Now apply Miller-Rabin algorithm 00072 // witnesses up to 4*10^9 00073 for (int i = 0; witnesses[i] != 0 && witnesses[i] < p; i++) { 00074 BaseType a = witnesses[i]; 00075 bool ok; 00076 /* According to the article 00077 When n = 1 + 2^h * d with d odd and h > 0 00078 and when n is a composite number, 00079 then n is called a "strong pseudoprime to base b" 00080 if either 00081 (2a) b^d = 1 mod n or 00082 (2b) b^{2^k*d} = -l mod n for 0 <= k < h 00083 */ 00084 BaseType t = PowerModImpl::powermod(a, d, p); 00085 ok = (t == 1); 00086 for (int r = 0; (r < s) && !ok; r++) { 00087 ok |= (t == p-1); 00088 // Optimization of computation (2b) by iterative squaring 00089 t = PowerModImpl::multmod(t, t, p); 00090 } 00091 00092 if (!ok) { 00093 return false; 00094 } 00095 } 00096 return true; 00097 } 00098 }; 00099 00100 // Can't use simple powermod as we are using 64bit integers. 00101 typedef PrimesFast_<math::powermod::PowermodExtended> PrimesFast; 00102 00103 } // namespace primes 00104 } // namespace math 00105 #endif