細かすぎて伝わらないエラトステネスの篩の高速化

5000兆回考えられたネタ。

エラトステネスの篩とは

指定された数以下の素数の一覧を求める単純なアルゴリズム。
計算量は $O(n\log\log n)$ らしい。

細かいことはググって。

ここでの条件

• C言語で実装
• エラトステネスの篩を使う
• 素数の数を求める（100 なら 25）
• コンパイラ最適化とかパイプラインとかには適度に期待する
• 最大値はそれなりに小さいとする

計算時間は手元の Mac で測っている（コンパイラは clang）。コンパイラオプションは -O3 -march=native。-march=native なのは popcnt と比較するため。

基本のコード

#include <stdio.h>
#define MAXVAL 200000000

int main() {
    int i, j, count = 0;
    static int flags[MAXVAL + 1] = {};
    for (i = 2; i <= MAXVAL; i++) {
        if (flags[i] == 0) {
            ++count;
            for (j = i * 2; j <= MAXVAL; j += i) {
                flags[j] = 1;
            }
        }
    }
    printf("%d\n", count);
    return 0;
}

5.54秒くらいかかる。おっせ〜。

高速化手法

外周の高速化

2 の倍数をとばす

あるあるの手法。2 の倍数を飛ばす場合、2 の倍数のフラグを立てなくて済むが、実はそこまで性能差はでない。

- int i, j, count = 0;
+ int i, j, count = 1;

- for (i = 2; i <= MAXVAL; i++) {
+ for (i = 3; i <= MAXVAL; i += 2) {

3 の倍数もとばす

2、4、2、4、… と飛ばしていけばよい。

- int i, j, count = 1;
+ int i, j, f, count = 2;

- for (i = 3; i <= MAXVAL; i += 2) {
+ for (i = 5, f = 4; i <= MAXVAL; i += (f = 6 - f)) {

最大値のルートまでしか回さない

100 までの場合、11 とか 13 とかについてはフラグを操作することがないので、これらを無視する。これで随分はやくなる。

#include <stdio.h>
#include <math.h>
#define MAXVAL 200000000

int main() {
    int i, j, f, count = 2;
    int max = (int)sqrt(MAXVAL) + 1;
    static int flags[MAXVAL + 1] = {};
    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!flags[i]) {
            ++count;
            for (j = i * 2; j <= MAXVAL; j += i) {
                flags[j] = 1;
            }
        }
    }
    for (; i <= MAXVAL; i += (f = 6 - f))
        if (!flags[i])
            ++count;
    printf("%d\n", count);
    return 0;
}

内周の高速化

開始インデックスと 2 の倍数

冷静に考えると、i * i から計算すれば十分である。なぜなら、たとえば 7 の倍数を消すとき、14 とか 35 とかは、それぞれ 2 と 5 を見ているときに消しているはずだからだ。

あと、2、4、…の倍数はフラグを起こさなくていいので、これも飛ばす。

これでもかなり速くなる。

-             for (j = i * 2; j <= MAXVAL; j += i) {
+             for (j = i * (i | 1); j <= MAXVAL; j += i * 2) {

後ろから見る

知っている合成数 × 今見ている素数 のパターンは、既にフラグが立っているはずなので、これはフラグを起こさなくてよい。
ただ、下から見てしまうと、このフラグを勘違いしてしまうことがある。
たとえば、11 を見ているとき、121 の合成数フラグを立ててしまうと、1331 のフラグを立てなくていいと勘違いしてしまう。
これを避けるために、後ろから for 文を回す。

- int i, j, f, count = 2;
+ int i, j, f, s, count = 2;

-             for (j = i * 2; j <= MAXVAL; j += i) {
-                 flags[j] = 1;
+             s = MAXVAL / i;
+             for (j = s - !(s & 1); j >= i; j -= 2) {
+                 if (!flags[j]) flags[j * i] = 1;

ちなみにこの高速化は、フラグをビット演算すると無駄になる（キャッシュの方が勝るので、どんどん書いてしまう方が速い）。

フラグをビット演算化する

フラグをビット演算化する

フラグは 0 と 1 しかないんだから、ビット演算にしろ。かなり高速になる。

#include <stdio.h>
#include <math.h>
#define MAXVAL 200000000
#define BITS (sizeof(int) * 8)
#define FLAGS_NUM (MAXVAL / BITS + 1)

static int flags[FLAGS_NUM] = {};

inline static void upflag(int i) {
    flags[i / BITS] |= 1 << (i % BITS);
}

inline static int getflag(int i) {
    return (flags[i / BITS] >> (i % BITS)) & 1;
}

int main() {
    int i, j, f, s, count = 2;
    int max = (int)sqrt(MAXVAL) + 1;
    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!getflag(i)) {
            ++count;
            s = MAXVAL / i;
            for (j = s - !(s & 1); j >= i; j -= 2) {
                if (!getflag(i))
                    upflag(j * i);
            }
        }
    }
    for (; i <= MAXVAL; i += (f = 6 - f))
        if (!getflag(i))
            ++count;
    printf("%d\n", count);
    return 0;
}

popcnt を使う

1 になっているビットの数を数える関数。2 と 3 の倍数のフラグを立てる必要があるので、これはテーブルを作ってしまう。

平方根までは count をインクリメントしてもいいかな、と思ったけどどっちが速いんだろう。

count = -2 は、0・2・3 を合成数扱いしてしまっているのと、1 を素数扱いしてしまっている分の埋め合わせ。

#include <stdio.h>
#include <math.h>
#include <x86intrin.h>
#define MAXVAL 200000000
#define BITS (sizeof(int) * 8)
#define FLAGS_NUM (MAXVAL / BITS + 1)

static int defaults[] = { 0x5D75D75D, 0xD75D75D7, 0x75D75D75 };
static int flags[FLAGS_NUM];

inline static void upflag(int i) {
    flags[i / BITS] |= 1 << (i % BITS);
}

inline static int getflag(int i) {
    return (flags[i / BITS] >> (i % BITS)) & 1;
}

int main() {
    int i, j, f, count = -2;
    int max = (int)sqrt(MAXVAL) + 1;

    for (i = 0; i < FLAGS_NUM; i++)
        flags[i] = defaults[i % 3];
    flags[FLAGS_NUM - 1] &= (1 << (MAXVAL % BITS + 1)) - 1;

    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!getflag(i)) {
            for (j = i * (i | 1); j <= MAXVAL; j += i * 2)
                upflag(j);
        }
    }

    for (i = 0; i < FLAGS_NUM; i++)
        count += _popcnt32(flags[i]);

    printf("%d\n", MAXVAL - count);
    return 0;
}

手元の環境では、これで 0.41 秒までになった。

そのほか

あまり意味のなかった手法

• 数値型を変える
  • unsigned にしたり 64bit にしたりしたけど、ほとんど変わらなかった
  • 組み合わせによっては若干速くなることがある
  • コンパイラ最適化の問題?
• flags の初期化
  • 実は {} はいらない
  • とるのは未初期化っぽい雰囲気が出るのであまり好きじゃない（わがまま）
  • ちなみにあってもなくてもほとんど変わらない（ないほうがほんの少し速い）

ダメだった手法

• while 文にする
  • for 文より while 文の方が速いことがあるが、スクリプト言語じゃないしまあ無意味
  • むしろ遅い
• 平方根の改良
  • math.h が意外と速い
  • 開平法を実装しても無意味だった

まだやってない手法

• 5 の倍数も飛ばす
  • 計算が複雑化して遅くなりそう（？）
  • 飛ばし方が 4, 2, 4, 2, 4, 6, 2, 6 となる（煩雑！）
  • このパターンを高速に作れれば速くなるかも
• 2 の倍数のフラグを作らない
  • 直接インデックスでアクセスできなくなって遅くなりそう（？）

まとめ

環境によっては選び方によって違うかも。なんか他にアイデアがあったらください。

Special Thanks: @dsk_saito、@D_Plius、@KRiver1、@lpha_z、@eukaryo