prime_sieve_benchmark.cpp

#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);
}

              normal: 0.975406 seconds
    normal_optimized: 0.503574 seconds
              linear: 0.663002 seconds
    linear_optimized: 0.3865 seconds
           segmented: 0.192711 seconds
 segmented_optimized: 0.145394 seconds
              aitken: 0.892283 seconds