Effective implementation of algorithms (Master Thesis)
Effective and error-free implementation of algorithms
src/math/prime_sieve/segmented_sieve.h
Go to the documentation of this file.
00001 
00013 #include "utils/preconditions/preconditions.h"
00014 #include "debug/ppdebug.h"
00015 #include "template/template.h"
00016 #include "math/prime_sieve/eratosthenes_basic.h"
00017 
00018 namespace math {
00019 namespace prime_sieve {
00020 
00026 class SieveCallback {
00027  public:
00028   virtual void foundPrime(long long int p) = 0;
00029 };
00030 
00041 class SegmentedSieve {
00042   private:
00046     static void sieve(
00047         const vector<int> &primes,
00048         long long int segment_start,
00049         long long int delta,
00050         long long int n,
00051         SieveCallback* callback) {
00052       Preconditions::check(segment_start % 2 == 1); 
00053       Preconditions::check(segment_start > 0.1 + sqrt(delta));
00054       // segment size is 2 * delta as we do not bother to sieve even values
00055 
00056       vector<long long int> L(primes.size());
00057       // Using vector<char> over vector<bool> gains 20% performance
00058       // and we do not need so tight memory optimization
00059       vector<char> sieve(delta, true);
00060 
00061       // prepare values of L
00062       FOR(i, primes.size()) {
00063         long long int p = primes[i];
00064         L[i] = (p - ((segment_start - p)/2) % p) % p;
00065       }
00066 
00067 
00068       while (segment_start < n) {
00069         // sieve all values in range [segment_start, segment_start + 2 * delta)
00070         FOR(i, primes.size()) {
00071           while (L[i] < delta) {
00072             sieve[L[i]] = false;
00073             L[i] += primes[i];
00074           }
00075           L[i] -= delta; // prepare L for next segment
00076         }
00077 
00078         for (int i = 0; i < delta; i++) {
00079           if (sieve[i]) {
00080             long long int p = segment_start + 2 * i;
00081             if (p >= n) break; // we should return only primes < n
00082             callback->foundPrime(segment_start + 2 * i);
00083           } else {
00084             sieve[i] = true; // reset sieve array to true
00085           }
00086         }
00087         segment_start += 2 * delta;
00088       }
00089 
00090     }
00091 
00092   public:
00100     static void findPrimes(long long int n, SieveCallback* callback) {
00101       Preconditions::check(n > 2);
00102       callback->foundPrime(2);
00103 
00104       // find delta >= sqrt(n)
00105 
00106       long long int start = 1 + (int) sqrt((double) n);
00107       if (start % 2 == 0) start++; // we need odd number
00108 
00109       // find all odd primes < start
00110       // (i.e. all primes which may factor numbers up to n)
00111       EratosthenesBasic e;
00112       // FIXME(ppershing): nepretecie?
00113       e.initialize((unsigned int) start);
00114 
00115       vector<int> primes;
00116       for (int i = 3; i < start; i += 2 ) {
00117         if (e.isPrime((unsigned) i)) {
00118           callback->foundPrime(i);
00119           primes.push_back(i);
00120         }
00121       }
00122 
00123       long long int delta = start;
00124       sieve(primes, start, delta, n, callback);
00125     };
00126 };
00127 
00128 } // namespace prime_sieve
00129 } // namespace math
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines