Effective implementation of algorithms (Master Thesis)
Effective and error-free implementation of algorithms
src/math/primes/primes_fast.h
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines