Effective implementation of algorithms (Master Thesis)
Effective and error-free implementation of algorithms
|
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