最近LINEの新卒採用試験ズバリ問題解説～アルゴリズム問題編～という記事が少し話題になっていた（というか競技プログラミング勢から突っ込まれしまい、後日お詫びが掲載された）のを見かけたので、そこで取り扱われていたエラトステネスの篩（ふるい）について、記事を書いてみようと思う。
エラトステネスの篩とは、ある上限値までの素数列を簡単に作成することができるアルゴリズムだ。
簡単に流れを説明すると、次のようなものである。
- 0と1は素数ではないので除外して、2～上限値までの数列を用意する。
- 残っている数の中で、まだ素数として確定していない一番小さい数（最初は2になる）を素数として確定する。
- さきほど素数として確定した数の倍数は、素数ではないものとして数列から除外する。
- 以下2.と3.を繰り返す。素数として確定させた数が上限の平方根以上になったところで、残っている数全てが素数となる。
具体的に、上限を20として手動で操作したときの様子を見てみよう。
数列を準備
[2, 3, 4, 5, …, 18, 19, 20]
2を素数として確定し、2の倍数を数列から除外
[2, 3, 5, 7, 9, 11, 13, 15, 17, 19]
3を素数として確定し、3の倍数を数列から除外
[2, 3, 5, 7, 11, 13, 17, 19]
次に5が素数として確定するが、これは20の平方根(4.47…)以上なので、（今回は存在しないが）5の倍数を取り除いて停止。残っている数列は、全て素数となる……ホントに？……ホントだ。
念のために書いておくと、なぜ上限の平方根までで大丈夫なのかというと、ある素数でない数nがp * q(ただしp <= qとする)として表される場合、pは必ずnの平方根以下になるため、ちゃんと除外されるからだ。
本質としては、まさしくふるいに掛けるように、素数以外の数を除外していけば残ったものが素数になるという、極めて分かりやすいアルゴリズムである。簡単な上、わりと高速に動くので、単純に素数列が欲しい時に利用しやすい。
実装する際には、単純に上記のやり方をまるごと再現してしまうと、リストの途中からの削除操作をかなりの回数行うことになるため、遅くなる。そこで、いちいちリストから数を削除するのではなく、各数について素数判定かどうかをTrue / Falseで持つような配列を用意してやり、それを操作するようにする。
具体的なコードは、以下のような感じになる。
# エラトステネスの篩を使い、上限limitまでの素数列を取得する
def get_prime_list(limit):
# nが素数なら、primep[n]==Trueとする配列を準備
primep = [True] * (limit + 1)
# 0, 1は素数から除外
primep[0], primep[1] = False, False
# 2～limitの平方根まで順番に見ていく
for n in range(2, int(limit ** 0.5) + 1):
# nが素数と確定したら、その倍数を全て素数から除外
if primep[n] == True:
for p in range(n * 2, limit + 1, n):
primep[p] = False
# 最後に素数だけを取り出す
return [p for p in range(limit + 1) if primep[p]==True]
僕の手持ちの環境で、limit=10**6の時に0.2秒、10**7で2.4秒くらいだろうか。Pythonのforループは遅いので有名だが、それでもそれなりに高速に動作する。
ちなみに、PyPyを使うと10**7でも0.4秒程度だった。PyPy偉い。
ところで、エラトステネスの篩よりも効率がいいアルゴリズムとして、アトキンの篩というものがある……らしいのだが、正直なところ僕の数学力（≒0）では原理がいまいち理解できなかった。というわけで、ここでは実装を紹介しない。分からないものは、実装もしないのだ。えっへん。
その代わりに、エラトステネスの篩をもう少しだけ（定数倍）改善してみたので、そのコードを紹介しておこうと思う。発想は簡単で、奇数だけを調べることにより、計算量を半分に削減しようというものだ。
# エラトステネスの篩を使うが、2の倍数をあらかじめ除いておく
def get_prime_list_ex(limit):
if limit < 2:
return []
# primep[n]==Trueのとき、(2 * n) + 1 が素数とする
primep = [True] * ((limit - 1) // 2 + 1)
primep[0] = False
# 少しややこしく見えるが、やっていることは同じ
# f.e. primep[3]==True(=7)の場合、
# primep[10](=21), primep[17](=35)...と
# 素数から除外されていくことになる
for n in range(1, int((limit ** 0.5) - 1) // 2 + 1):
if primep[n] == True:
p = 2 * n + 1
for i in range(n + p, len(primep), p):
primep[i] = False
# 2だけはしょうがないので最後に追加する
return [2] + [2 * p + 1 for p in range(len(primep)) if primep[p]==True]
これなら、limit=10**7でも1.1秒、PyPyを使えば0.25秒程度なので、もう少しだけ実行時間を詰めたいというときには使えるかもしれない。実装も、上記の通りそれほど難しくはない。
この調子でもっといけるんじゃね？と思い、2, 3の倍数を除くパターン……みたいなのも考えてみたのだが、6周期で余りが1, 5の整数をチェックすることになり、実装が難しそうな上に、そこまで効果が見込めなさそうなので、こっちはとりあえず構想だけに留めておいた。
今回は、以上。
追記：
2個めのバージョンについて、実はもっと簡潔に書けるのではないかと思い、こういうのも書いてみた。
# エラトステネスの篩を使うが、2の倍数をあらかじめ除いておく
def get_prime_list_ex2(limit):
if limit < 2:
return []
primep = [True] * (limit+1)
# 3以上の奇数のみを順番に見ていく
for n in range(3, int(limit ** 0.5) + 1 , 2):
if primep[n] == True:
# 素数の奇数倍のみ書き換えすればOK
for i in range(n * 3, limit + 1, n * 2):
primep[i] = False
return [2] + [p for p in range(3, limit+1, 2) if primep[p]==True]
こちらの方が見通しが良く簡単だが、ごくわずかに(10%程度）遅い。まあ、誤差程度なのでお好みで……といったところだろうか。