Created
November 14, 2018 12:11
-
-
Save jakobkogler/7b251da1a509e7c5eac69c34b169b765 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <bits/stdc++.h> | |
using namespace std; | |
const int N = 1e8; | |
bool comp[N]; | |
vector<int> primes; | |
void normal_sieve() { | |
memset(comp, 0, sizeof comp); | |
primes.clear(); | |
bool done = false; | |
for (int i = 2; i < N; i++) { | |
if (!comp[i]) { | |
primes.push_back(i); | |
if (i * i > N) { | |
done = true; | |
} | |
if (!done) | |
for (int j = i * i; j < N; j += i) | |
comp[j] = 1; | |
} | |
} | |
} | |
void normal_optimized_sieve() { | |
memset(comp, 0, (sizeof comp) / 2); | |
primes.clear(); | |
primes.push_back(2); | |
int n = N / 2; | |
int sr = round(sqrt(N)); | |
for (int i = 1; i < n; i++) { | |
if (!comp[i]) { | |
int p = 2 * i + 1; | |
primes.push_back(p); | |
if (p > sr) | |
continue; | |
for (int j = (p + 1) * i; j < n; j += p) | |
comp[j] = 1; | |
} | |
} | |
} | |
void linear_sieve() { | |
memset(comp, 0, sizeof comp); | |
primes.clear(); | |
for (int i = 2; i < N; i++) { | |
if (!comp[i]) | |
primes.push_back(i); | |
for (int j = 0, si = primes.size(); j < si && i * primes[j] < N; j++) { | |
comp[primes[j] * i] = 1; | |
if (i % primes[j] == 0) | |
break; | |
} | |
} | |
} | |
void linear_sieve_optimized() { | |
memset(comp, 0, (sizeof comp) / 2); | |
primes.clear(); | |
primes.push_back(2); | |
for (int i = 1, p = 3; p < N; i++, p+=2) { | |
if (!comp[i]) | |
primes.push_back(p); | |
for (int j = 1, si = primes.size(); j < si && p * primes[j] < N; j++) { | |
comp[(primes[j] * p) >> 1] = 1; | |
if (p % primes[j] == 0) | |
break; | |
} | |
} | |
} | |
void segmented_sieve() { | |
int n = N; | |
const int S = 50000; | |
primes.clear(); | |
int nsqrt = round(sqrt(n)); | |
vector<char> is_prime(nsqrt + 1, true); | |
::primes.push_back(2); | |
vector<char> block(S); | |
vector<int> primes, start_indices; | |
for (int i = 3; i <= nsqrt; i += 2) { | |
if (is_prime[i]) { | |
primes.push_back(i); | |
start_indices.push_back(i * i); | |
for (int j = i * i; j <= nsqrt; j += 2 * i) | |
is_prime[j] = false; | |
} | |
} | |
for (int k = 0; k * S <= n; k++) { | |
fill(block.begin(), block.end(), true); | |
int start = k * S; | |
for (auto i = 0u; i < primes.size(); i++) { | |
int p = primes[i]; | |
int idx = start_indices[i]; | |
for (; idx < S; idx += 2 * p) | |
block[idx] = false; | |
start_indices[i] = idx - S; | |
} | |
if (k == 0) | |
block[1] = false; | |
for (int i = 1; i < S && start + i <= n; i += 2) { | |
if (block[i]) | |
::primes.push_back(start + i); | |
} | |
}; | |
} | |
void segmented_sieve_optimized() { | |
int n = N; | |
const int S = 30000; | |
primes.clear(); | |
int nsqrt = round(sqrt(n)); | |
vector<char> is_prime(nsqrt + 1, true); | |
vector<int> primes, start_indices; | |
for (int i = 3; i <= nsqrt; i += 2) { | |
if (is_prime[i]) { | |
primes.push_back(i); | |
start_indices.push_back((i * i - 1) / 2); | |
for (int j = i * i; j <= nsqrt; j += 2 * i) | |
is_prime[j] = false; | |
} | |
} | |
::primes.push_back(2); | |
vector<char> block(S); | |
int high = (n - 1) / 2; | |
for (int low = 0; low <= high; low += S) { | |
fill(block.begin(), block.end(), true); | |
for (auto i = 0u; i < primes.size(); i++) { | |
int p = primes[i]; | |
int idx = start_indices[i]; | |
for (; idx < S; idx += p) | |
block[idx] = false; | |
start_indices[i] = idx - S; | |
} | |
if (low == 0) | |
block[0] = false; | |
for (int i = 0; i < S && low + i <= high; i++) { | |
if (block[i]) | |
::primes.push_back((low + i) * 2 + 1); | |
} | |
}; | |
} | |
void aitken() { | |
int n = N; | |
std::vector<bool> m(n + 1, false); | |
int root = std::sqrt(n) + 1; | |
for (int i = 1; i <= root; i++) { | |
for (int j = 1; j <= root; j++) { | |
int a = 4 * i * i + j * j; | |
int b = 3 * i * i + j * j; | |
int c = 3 * i * i - j * j; | |
if (a <= n && (a % 12 == 1 || a % 12 == 5)) | |
m[a].flip(); | |
if (b <= n && b % 12 == 7) | |
m[b].flip(); | |
if (i > j && c <= n && c % 12 == 11) | |
m[c].flip(); | |
} | |
} | |
for (int i = 5; i * i <= n; i++) { | |
if (m[i]) { | |
for (int j = 1; j * i * i <= n; j++) | |
m[j * i * i] = false; | |
} | |
} | |
primes.clear(); | |
primes.push_back(2); | |
primes.push_back(3); | |
for (int i = 5; i < n; i++) { | |
if (m[i]) | |
primes.push_back(i); | |
} | |
} | |
template <typename T> | |
void time_compare(string name, T func, vector<int> comparison) { | |
primes.clear(); | |
auto t0 = clock(); | |
(*func)(); | |
auto t1 = clock(); | |
if (!comparison.empty()) | |
assert(comparison == primes); | |
cout << setw(20) << name << ": " << (1.0 * (t1 - t0)) / CLOCKS_PER_SEC << " seconds" | |
<< endl; | |
} | |
int main() { | |
::primes.reserve(6'000'00); | |
normal_sieve(); | |
auto comparison = primes; | |
time_compare("normal", normal_sieve, comparison); | |
time_compare("normal_optimized", normal_optimized_sieve, comparison); | |
time_compare("linear", linear_sieve, comparison); | |
time_compare("linear_optimized", linear_sieve_optimized, comparison); | |
time_compare("segmented", segmented_sieve, comparison); | |
time_compare("segmented_optimized", segmented_sieve_optimized, {}); | |
time_compare("aitken", aitken, comparison); | |
} |
Author
jakobkogler
commented
Nov 14, 2018
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment