高速化
素数
エラトステネスの篩

エラトステネスの篩の高速化手法まとめ


はじめに

エラトステネスの篩は「整数の集合から素数の倍数(合成数)をふるい落とし、素数の集合を得る」アルゴリズムです。

この記事ではエラトステネスの篩の高速化手法を実装し、どれくらい速くなるのか実験しました。

完全なコードはGitHubに置いてあります。


参考資料

エラトステネスの篩の高速化について書かれた他の記事です。

以下の記事を整理して自分用にまとめたいと思ったのが、この記事を書くきっかけになりました。

高速化の話は以下が分かりやすく、いろいろ参考にさせていただきました。


エラトステネスの篩の実装例

以下はエラトステネスの篩のナイーブな実装です。

$N$ = 200,000,000 までの整数を篩にかけ、素数の数を出力します。

flagsは値 0(素数)で初期化された配列で、素数iの倍数i * jに対応する要素に値 1(合成数)がセットされます。

sieveは合成数のふるい落とし、countは素数の数を数える関数です。


sieve1.c

#include <math.h>

#include <stdint.h>

#define N (200000000)

#define IS_PRIME(n) (flags[n] == 0)
#define SET_COMPOSITE(n) (flags[n] = 1)

static char flags[N + 1];

void sieve(void) {
for (uint32_t i = 2, h = sqrt(N); i <= h; ++i)
if (IS_PRIME(i))
for (uint32_t j = i, h = N / i; j <= h; ++j)
SET_COMPOSITE(i * j);
}

uint32_t count(void) {
uint32_t c = 0;

for (uint32_t i = 2; i <= N; ++i)
if (IS_PRIME(i))
++c;

return c;
}


以下はmain関数から上記を呼び出す部分です。

このコードは以降も同じものを使うので、今後は省略します。


main.c

#include <inttypes.h>

#include <stdint.h>
#include <stdio.h>

extern void sieve(void);
extern uint32_t count(void);

int main(void) {
sieve();
printf("%" PRIu32 "\n", count());

return 0;
}



Wheel factorization

以下は 2, 3 の倍数をふるい落とすことを省略して計算量を減らしたものです。

6 (= 2 * 3) の周期では 6n ± 1 以外には素数が現れないので、外側のループでiの増分を 2, 4, 2, 4, ... とすることで素数になり得ない数を読み飛ばしています。

このアルゴリズムを一般化したものは Wheel factorization と呼ばれています。

flags[2, 3 の倍数]は 0(素数)のままなので、この部分を誤って素数としてカウントしないようにcount関数も修正する必要があります(省略)。


sieve2.c

...

void sieve(void) {
for (uint32_t i = 5, g = 2, h = sqrt(N); i <= h; i += g, g ^= 6)
if (IS_PRIME(i))
for (uint32_t j = i, h = N / i; j <= h; j += 2)
SET_COMPOSITE(i * j);
}

...



省メモリ化

以下はflagsの型を符号なし 64 ビット整数にして、ある数が素数かどうかを 1 ビットで保持するよう変更したものです。

flagsから数値に対応するビットを取り出す、および合成数のフラグをたてる操作はビット演算になります。


sieve3.c

...

#define IS_PRIME(n) (((flags[(n) / 64] >> (n) % 64) & 1) == 0)
#define SET_COMPOSITE(n) (flags[(n) / 64] |= (UINT64_C(1) << ((n) % 64)))

static uint64_t flags[N / 64 + 1];

...



Segmented sieve(区間篩)

ナイーブな実装ではflags配列の先頭近くから末尾までを順番にアクセスし、それを素数の数だけ繰り返していました。

これはキャッシュミスが多発するので遅いです。

以下はメモリアクセスの順番を変えて上記の点を改善した「区間篩」のアルゴリズムを実装したものです。

$\sqrt N$まで普通の篩にかけた(最初のループ)後、二番目のループでは $B$ = 0x8000 ずつの範囲に区切って篩にかけることで、メモリアクセスを局所化しています。


sieve4.c

...

#define B (0x8000)

...

void sieve(void) {
uint32_t s = sqrt(N);

for (uint32_t i = 2, h = sqrt(s); i <= h; ++i)
if (IS_PRIME(i))
for (uint32_t j = i, h = s / i; j <= h; ++j)
SET_COMPOSITE(i * j);

for (uint32_t b = s + 1; b <= N; b += B)
for (uint32_t i = 2, h = sqrt(b + B); i <= h; ++i)
if (IS_PRIME(i))
for (uint32_t j = b / i, h = (b + B) / i; j < h; ++j)
SET_COMPOSITE(i * j);
}

...



実行時間の計測

以下の条件でコンパイルして実行時間を計測しました。


  • PC


    • MacBook(以下は「このMacについて」メニューで確認)

    • プロセッサ: 1.3 GHz Intel Core M

    • メモリ: 8 GB 1600 MHz DDR3



  • コンパイル


    • clang

    • -O3 オプション



  • 計測方法


    • time コマンド

    • 3 回実行して平均値をとる



結果は以下のようになりました。

sieve1.c
sieve2.c
sieve3.c
sieve4.c

real
0m2.601s
0m1.627s
0m1.973s
0m1.093s

user
0m2.509s
0m1.535s
0m1.954s
0m1.007s

sys
0m0.081s
0m0.084s
0m0.014s
0m0.079s


まとめ

局所的でないメモリアクセスを避けましょう。キャッシュミスのペナルティは高くつきます。

データを小さくするだけでもキャッシュミスが減るので結果的に速くなります。

今回紹介したアルゴリズムはすべて組み合わせることが可能であり、そうすることでより速くなるでしょう。