エラトステネスの定数倍高速化術

　Project Eulerのとある問題でエラトステネスの篩を高速化している人がいた。
なかなか興味深かったので、コードを読みといて自分なりにアレンジしてみた。

ローカル環境で約2倍の高速化に成功した。1e8までのふるいが0.7s程度で実行出来た。

普通のふるい

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 1e8;
bool pr[N];

int main(int argc, char **argv) {
    
    fill(pr + 2, pr + N, true);

    for (int i = 0; i * i < N; i++) {
        if (!pr[i])
            continue;

        for (int j = i * i; j < N; j += i)
            pr[j] = false;
    }

    int x = 0;
    int y = 0;
    for (int i = 0; i < N; i++)
        if (pr[i])
            x ^= i, y = 23 * y + i;
    
    cout << x << endl;
    cout << y << endl;

    return 0;
}

定数倍高速版ふるい

アイディアとしては、

1. 2以外の偶数は素数ではないため、ふるいの配列は奇数に対してのみ用意する。
2. 素数か合成数かのフラグをboolではなくintのビットで管理する。

というもの。

1.で約2倍高速になった。これは直感的にも理解しやすい。一番内側の処理の実行回数が半分程度になるため2倍程度速くなりそうである。

1.の後に2.を施すことで、さらに1.4倍ほど高速になった。32ビット毎に何らかの処理を纏めたというわけではなく、単純にフラグの管理をビットで行うようにしただけであるが意外と速くなった。使用するメモリがコンパクトになったためキャッシュに乗りやすくなったということだろうか。

#include <iostream>
#include <algorithm>

using namespace std;

const int N = 1e8;
const int BITS = 32;

#define BACKET(n) (n / 2 / BITS)
#define MASK(n) (1 << n / 2 % BITS)

int pr[BACKET(N) + 1];    // sieve arrray containing only odd numbers  [1,3,5,7,,9,....]

int main(int argc, char **argv) {
    
    fill(pr, pr + BACKET(N) + 1, -1);
    pr[0] = -2;   // 1 is not prime

    for (int i = 1; i * i < N; i += 2) {
        if (!(pr[BACKET(i)] & MASK(i)))
            continue;
        for (int j = i * i; j < N; j += 2 * i)
            pr[BACKET(j)] &= ~MASK(j);
    }

    int x = 0;
    int y = 0;
    for (int i = 0; i < N; i++)
        if (i == 2 || i % 2 && (pr[BACKET(i)] & MASK(i)))
            x ^= i, y = 23 * y + i;

    cout << x << endl;
    cout << y << endl;

    return 0;
}