世の中に多くあるエラトステネスの篩の実装、多くあるくせにちょっとしか高速化してないのが悲しいので高速化エラトステネスの篩を書いてみることにします。7 本目(利用者が限られる)を書くかどうか迷い中。

問題設定

とりあえず C++11 くらいで動く物を作りたいと思います。インラインアセンブラや SIMD、並列化は明示的には入れない方針です。あと、エラトステネスの篩自体の高速化を目指すので、結果は一部の具体的な素数の確認と $\pi(x)$ の値で確認することにしますし、その時間は考慮しないことにします。
また、ライブラリ的に使えるよう、将来的には区間篩 (大きな $x$ に対して $[x, x+y)$ の区間での素数を求める) も計測することにします。

アルゴリズム

タイトルにあるエラトステネスの篩のアルゴリズムをざっくり説明すると

2 以上 $x$ 未満の整数に対して $True$ を用意する
$\sqrt{x}$ 以下の数 $p$ を小さい方から見ていく
1. $p$ の対応が $True$ なら、$p^2$ 以上の $p$ の倍数の対応を $False$ とする
2 で選ばれた $p$、並びに残った $True$ に対応する整数が素数である

という形になります。基礎的なアルゴリズム説明だと、$x$ まですべての整数で篩う説明だったり、篩い始める箇所を $2p$ にしたりしますが、今回はそういう定番高速化は最初から導入しています。プログラムをざっくり書くと

void Sieve(int64 x) {
  std::vector<bool> flags(x, true);
  flags[0] = false;
  flags[1] = false;
  const int64 sqrt_x = std::ceil(std::sqrt(x) + 0.1);  // 少し誤差を考えてる
  for (int64 p = 2; p < sqrt_x; ++p) {
    if (!flags[p])
      continue;
    for (int64 mult = p * p; mult < x; mult += p)
      flags[mult] = false;
  }
}

という形になります。これをベースに今後高速化していきます。

偶数の省力化

この高速化はよく知られた方法で、素直に実装すればだいたい 2 倍速になり、素直に実装していれば使うメモリの量も半減します。そもそも最初の 2 で払い落としてる ($False$ に設定する) 数は対象とする数の半分あるので、それらを計算対象から外すわけなので考察もそうむずかしくありません。
新たな実装で flags[i] は整数 $2i+1$ に対応しています。なので　$p=2i+1$ に対して $p^2=4i^2+4i+1=2(2i(i+1))+1$ となるので、篩い落とす最初のインデックスは 2*i*(i+1) となります。

void Sieve(int64 x) {
  std::vector<bool> flags(x / 2, true);
  flags[0] = false;
  const uint64 sqrt_x = std::ceil(std::sqrt(x) + 0.1);
  const uint64 sqrt_xi = sqrt_x / 2;
  for (uint64 i = 1; i < sqrt_xi; ++i) {
    if (!flags[i])
      continue;
    const uint64 p = 2 * i + 1;
    for (uint64 mult = 2 * i * (i + 1); mult < flags_.size(); mult += p)
      flags_[mult] = false;
  }
}

応用として 3 や 5 の倍数を除く処理も知られていますが、その処理については次回。
ちなみに今回(を含む一連)のプログラムは github に纏めています。今回のは version 0 と 1 に該当します。

計算時間

目安として、計算時間を記録しておきます。単位は秒。

$x$	$\pi(x)$	Version 0	Version 1
$10^7$	664579	0.047	0.017
$10^8$	5761455	1.039	0.352
$10^9$	50847534	-	6.180

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

問題設定

アルゴリズム

偶数の省力化

計算時間