32bitCPUで2^32未満の素数を列挙する

PentiumM-1.1GHzで148秒。
タスクマネージャで見たメモリ使用量は144MB程度。

//2,3,5の倍数以外の整数
const unsigned long CPrime::p[] = {1, 7, 11, 13, 17, 19, 23, 29};
//逆引き
const int CPrime::q[] = {
    0, 1, 0, 0, 0, 0, 0, 2, 0, 0,
    0, 3, 0, 4, 0, 0, 0, 5, 0, 6,
    0, 0, 0, 7, 0, 0, 0, 0, 0, 8,
};

CPrime::CPrime()//2^32未満の素数を求める
{
    unsigned long h, i, j, k, m, u, t;
    unsigned long _j;

    m = 0xffffffff / 30 + 1;
    bb = new char [m];//2,3,5の倍数を除いたビットマップ、136.5MB
    for (i = 0; i < m; i++){
        bb[i] = 0;
    }
    bb[0] = 1;//1は素数でない（素数でないビットを立てる）
    
    k = (0x10000 - 1) / 30 + 1;//2^16未満の素数の倍数を除去する
    for (u = 0; u < k; u++){
        for (h = 0; h < 8; h++){
            if (bb[u] >> h & 1) continue;//素数でないことが既知なら次へ
            i = u * 30 + p[h];//これが新たに素数と判明した数
            if (i >= 0x10000) break;//i * iのオーバーフロー回避
            for (_j = 0, j = i * i; _j < j; _j = j, j += i * 2){
                (t = q[j % 30]) && (bb[j / 30] |= 1 << (t - 1));
            }
        }
    }
}

エラトステネスのふるいを使う。
i * i未満の合成数は、必ずi未満の素数の倍数なので、j = i * i とする。
（ていうかエラトステネスのふるい自体がそういう手法）
また、iは奇素数なので、i * iも奇数、i * 2は偶数。
ビットマップから偶数は除いてあるので、j += i * 2 とする。
この2点は、簡単に書けるし、速度が何割か違ってくる重要ポイント。
更に、unsigned longで表せる最大数までの判定を行うため、少し処理を追加する。
整数30個単位で処理する都合上、iが2^16を超える危険があるため、
if (i >= 0x10000) break; としている。
jはiの倍数だが、これが2^32以上になったことを普通には検知できないので、
「i * 2を足したのに前回のjより小さくなった」ことで判定している。
（符号なし整数加算でオーバーフローすると0に戻るので）
C++レベルだと面倒だが、原理的には、
アセンブラレベルならCF（キャリーフラグ）を使って簡単高速に分岐できる。
こうして作ったビットマップで素数を判定するには次のようにする。

bool isPrime(unsigned long i)
{
    int t;
    bool b = (t = q[i % 30]) && !(bb[i / 30] >> (t - 1) & 1);
    return b || i == 2 || i == 3 || i == 5;
}

対象のiを30で割った余りでテーブルを引いて、
0ならビットマップを見るまでもなく非素数。
それ以外なら示されたビット位置を見て、ビットが0なら素数。
また、この方法では2,3,5の倍数の素数を判定できないため、
2,3,5については個別の条件を追加している。
これは多くの場合に無駄な処理となるだろうから、
ケースバイケースでこの処理を省略したバージョンを使うようにする。
ここで重要なのが、30で割る処理。
定数除算は乗算のみで可能なので、コンパイラにこの最適化をしてもらう。
（VC++6.0ではなぜか最適化してくれなかった。2010なら大丈夫）
この最適化をしないと、かなり（倍近く）遅くなってしまい、
速度的には、わざわざ3,5の倍数を除いた意味がなくなってしまう。
それでも、判定にかかる時間は、2の倍数のみを除いたビットマップと比べて
2〜3倍遅くなるので、メモリに余裕があって大量の判定をする場合は、
どの方法を使うのがベストか考える必要がある。
絶対的な時間は短いので、通常は問題ないと思う。
これを使って求めた2^32未満の素数の個数は203280221個（コードは下）。
素数・個数・2^32などで検索しても確認することができなかったが、
逆に203280221で検索すればよかった。
http://oeis.org/A007053

unsigned long i, s = 0;
CPrime p;
for (i = 1; i; i++){
    if (p.isPrime(i)) s++;
}

このコード自体は高速化の余地があるけど、isPrime()の性能を見るためなので。
これは42億回（2^32-1）のループで89秒だった。
1.1GHzは毎秒11億クロックなので、1回の判定は22〜23クロック程度。
ただし、これはキャッシュのヒット率が高いケースと思われる。
2,3,5の判定を省けば、もう少し速くなる。