Implementing Sieve of Eratosthenes

エラトステネスの櫛はN以下の素数を列挙するアルゴリズムとして知られています。まずは素数を列挙するナイーブな実装を行い、計算量を落としていきながらエラトステネスの櫛を考えてみることにします。

次の実装の計算量はO(N²)です。

#include <iostream>
#include <vector>


#define N 100000


int main(int argc, char** argv)
{
  std::vector<bool> is_prime(N + 1, true);

  is_prime[0] = is_prime[1] = false;

  for (int i = 2; i <= N; i ++)
    for (int j = 2; j < i; j ++)
        if (i % j == 0) {
          is_prime[i] = false;

          break;
        }

  return 0;
}

この実装では、

2〜Nまでの全ての数iに関して、
2〜i-1までの数で割り切れるかどうか

を調べています。とても遅そうです。手元の端末ではNが10⁵のような比較的小さな値でも2秒かかりました。試しにNを10⁶としたら終わりませんでした。

次にO(N√N)の実装をしましょう。

                  :

  for (int i = 2; i <= N; i ++)
    for (int j = 2; j * j <= i; j ++)
        if (i % j == 0) {
          is_prime[i] = false;

          break;
        }

                  :

変更したのは2行目だけです。この実装では先日のエントリで紹介した知見を利用して、

2〜Nまでの全ての数iに関して、
iが2〜√iまでの整数で割り切れるかどうか

を調べています。

先ほどのナイーブな実装に比べるととても速く、Nが10⁵のとき、0.15秒で終わりました。しかし、Nが10⁶のときはまだ2.9秒ほどかかるようです。

次にエラトステネスの篩の要素を導入してみましょう。

素数の整数倍の数は素数ではありません。そのため、素数を見つける度にその整数倍の数が「素数ではない」と記録しておけば将来的にその数の素数判定をする必要がなくなるため無駄な計算が省けそうです。文字通り、探索の候補から「篩い落とす」のです。

                  :

  for (int i = 2; i <= N; i ++)
    if (is_prime[i])
      for (int j = 2; j * j <= i; j ++)
        if (i % j == 0) {
          for (int k = i; k <= N; k += i)
            is_prime[k] = false;

          break;
        }

                  :

この方法にすることによって、Nが10⁶のときでも0.3秒で終わるようになりました。

考えてみるとこの実装にも無駄があることが分かります。

素数の倍数は予め探索から篩い落とされています。数iに訪れたときにその数がまだ「素数ではない」と判定されていなかったとすると、数iは素数である、と言えます。その数を篩うような約数がなかったからです。そのため数iが素数であるかを改めてチェックする必要もないのです。

エラトステネスの櫛の動作原理をもう少し深く理解するために作図しました(図を考えるにあたって、Sieve of Eratosthenes - Wikipedia, the free encyclopediaを参考にしました。アニメーションを使った良質な資料です。合わせてご参照ください)。

初めに2からNまでの状態を記録する領域を確保します。

まず2を処理します。2には何も記録されていません。そのため、2を素数とします。また、将来的に処理するであろう2の整数倍の数は全て素数ではありません。これを記録します。

次に3を処理します。3は2で割り切れる数ではありませんでした。そのため、3には何も記録されていません。そこで、3を素数として記録します。将来的に処理するであろう3の整数倍の数も全て素数ではありません。先ほどと同様に、これも記録しましょう。

次に4を処理しますが、4は「素数ではない」と記録されています。そのため、単純にスキップします。

次は5です。5には何も記録されていないため、素数です。3と同様の処理をします。

6からもこれを繰り返します。

説明では「記録されていない」、「素数である」、「素数ではない」の3つの状態を考えていますが、実際には「素数である」、「素数ではない」の2つの状態のみを扱えばアルゴリズムとして成立します。その場合、初めに記憶領域は全て「素数である」の状態にしておきます。

ここまでを実装に反映したものが以下です。

                  :

  for (int i = 2; i <= N; i ++)
    if (is_prime[i])
      for (int k = i + i; k <= N; k += i)
        is_prime[k] = false;

                  :

その効果は絶大で、Nが10⁶で0.11秒で終わるようになりました。

さて、続けて7を処理した状態を考えてみましょう。

先ほどの実装では7に続き、11、13、...と処理していますが、実はその必要はありません。Nまでの素数を調べるのに2〜√Nまでを調べればよいのはエラトステネスの櫛も同様なのです。ここではNを120としていますので、√Nは10.95ほどです。そのため、10まで考えればいいことになります。

例えば120に近い113について考えてみましょう。√113は10.63です。先ほど計算した10.95よりも小さいですね。同様に、109、107、...、11までの平方根はそれぞれ11よりも小さいので、やはり高々10までを調べれば十分なのです。

今ひとつピンと来ない人には以下のように形式的に記述したもののほうが分かりやすいかもしれません。

ここまでの考察を元に、エラトステネスの篩の最終版を実装します。

                  :

  for (int i = 2; i * i <= N; i ++)
    if (is_prime[i])
      for (int k = i + i; k <= N; k += i)
        is_prime[k] = false;

                  :

WikipediaのSieve of Eratosthenesによれば計算量はO(NloglogN)だそうです。実装も随分簡素になりました。速度も速く、Nが10⁶のときは0.006秒で終わるようになりました。Nが10⁷であっても0.05秒ほどで終わり、10⁸であっても0.74秒で終わります。