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

5000兆回考えられたネタ。

# エラトステネスの篩とは

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

細かいことはググって。

# ここでの条件

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

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

# 基本のコード

```c
#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 の倍数のフラグを立てなくて済むが、実はそこまで性能差はでない。

```diff
- 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、… と飛ばしていけばよい。

```diff
- 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 とかについてはフラグを操作することがないので、これらを無視する。これで随分はやくなる。

```c
#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、…の倍数はフラグを起こさなくていいので、これも飛ばす。

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

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

### 後ろから見る

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

```diff
- 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 しかないんだから、ビット演算にしろ。かなり高速になる。

```c
#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 を素数扱いしてしまっている分の埋め合わせ。

```c

#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](https://twitter.com/dsk_saito)、[@D_Plius](https://twitter.com/D_Plius)、[@KRiver1](https://twitter.com/KRiver1)、[@lpha_z](https://twitter.com/lpha_z)、[@eukaryo](https://twitter.com/eukaryo)