本当は遅い「似非エラトステネスの篩」の罠
この記事は？
インターネット上でエラトステネスの篩の実装を検索すると、かなりの割合で、エラトステネスの篩とよんでいいのか怪しい「似非エラトステネスの篩」とでも称すべきものがみられます。
この記事では、「試し割り」「似非エラトステネスの篩」「エラトステネスの篩」の3つのアルゴリズムを比較して、その違いを解説します。
3つのアルゴリズムの比較
1. 試し割り
実装
まず他のアルゴリズムに対する評価基準として、試し割りのアルゴリズムを示します。
試し割りは、
import sys
# limit 以下の全ての素数を返す
def list_primes(limit):
primes = []
for i in range(2, limit + 1):
is_prime = True
for p in primes:
if p * p > i:
break
elif i % p == 0:
is_prime = False
break
if is_prime:
primes.append(i)
return primes
if __name__ == '__main__':
n = int(sys.argv[1])
print(len(list_primes(n)))
ベンチマーク
10000000以下の素数の数を求める時間を測ります。
$ python3 division.py 10000000
664579
python3 division.py 10000000 17.35s user 0.06s system 99% cpu 17.423 total
所要時間は17.35秒でした。
2. 似非エラトステネスの篩
実装
エラトステネスの篩は、数のリストの中から合成数を取り除いていくアルゴリズムです。
この原理を素朴にpythonで実装すると以下のようになるでしょう。
import bisect
import sys
# limit 以下の全ての素数を返す
def list_primes(limit):
primes = []
nums = [i for i in range (2, limit + 1)]
while len(nums) != 0:
p = nums.pop(0)
primes.append(p)
if p * p > limit:
break
split_index = bisect.bisect_left(nums, p*p)
small_nums = nums[:split_index]
large_nums = nums[split_index:]
nums = small_nums + [x for x in large_nums if x % p != 0]
return primes + nums
if __name__ == '__main__':
n = int(sys.argv[1])
print(len(list_primes(n)))
ベンチマーク
試し割りと同じ条件で時間を測ります。
$ python3 bad_sieve.py 10000000
664579
python3 bad_sieve.py 10000000 37.85s user 0.21s system 99% cpu 38.071 total
なんと試し割りの倍の時間がかかってしまいました。
何がいけないのか？
このアルゴリズムは、エラトステネスの篩の アイデアの実装としては 間違っていないのですが、 試し割りよりも高速なアルゴリズムであるとするのは 明確な間違い です。
実のところ、このアルゴリズムと試し割りは「
具体例として、「35が合成数と判明するまでの流れ」を考えてみましょう。
試し割りであれば、i が 35 になった時点におけるprimesの値は[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]になります。
そして、for p in primes: のループの中で、35 % 2 35 % 3 35 % 5 が順に評価され、35 % 5が0となることで35が合成数と判明します。
似非篩の場合はどうでしょうか。
まず、ループの1周目ではpが2であり、[x for x in large_nums if x % p != 0] の中で 35 % 2 が評価され、これは0でないので 35 はnumsに残ります。
次に、ループの2周目ではpが3であり、[x for x in large_nums if x % p != 0] の中で 35 % 3 が評価され、これは0でないので 35 はnumsに残ります。
そして、ループの3周目ではpが5であり、[x for x in large_nums if x % p != 0] の中で 35 % 5 が評価され、これは0なので35はnumsから消えます。
重要なのは、似非エラトステネスの篩の場合においても、試し割りの場合と同様に35 % 2 35 % 3 35 % 5 が全て評価されるという点です。
そしてこれは35に限らず全ての数について成立します。
結果、ここに挙げた二つの実装に同じ入力を与えたときに行われる剰余計算の回数は、全く同じになります。
実際剰余計算をするときに出力をしてみれば確認できますが、両者の違いは剰余の評価を「左辺が小さい順に実行する」か、「右辺が小さい順に実行する」か、でしかないのです。
そして割り算の回数が同じであれば、わざわざ全ての数を一旦リストに入れた上にそのリストを何度も加工するというオーバーヘッドがある「似非エラトステネスの篩」が、試し割りよりも遅いのは当然といえるでしょう。
3. エラトステネスの篩（本物）
実装
「正しい」エラトステネスの篩では、値の削除はリストからの削除ではなく、フラグの変更によって表現します。
import sys
# limit 以下の全ての素数を返す
def list_primes(limit):
primes = []
is_prime = [True] * (limit + 1)
is_prime[0] = False
is_prime[1] = False
for p in range (0, limit + 1):
if not is_prime[p]:
continue
primes.append(p)
for i in range(p*p, limit + 1, p):
is_prime[i] = False
return primes
if __name__ == '__main__':
n = int(sys.argv[1])
print(len(list_primes(n)))
ベンチマーク
試し割りと同じ条件で時間を測ります。
$ python3 good_sieve.py 10000000
664579
python3 good_sieve.py 10000000 1.37s user 0.04s system 99% cpu 1.411 total
試し割りに対して10倍以上高速になりました。
なぜ速い？
まず、リストのサイズ変更を行っていないので、領域の再確保やコピーの手間が発生しません。
そしてより重要なのは、リストのサイズが変わらないことで、pの倍数を除くための作業が「pごとにマークする」だけで済むということです。
試し割りや似非篩においては35 % 2のような割り切れないペアも評価する必要がありますが、正しい篩では割り切れないペアに対する処理は そもそも発生しない のです。
したがって ループの実行回数が少なく なり、試し割りよりも高速に素数を求めることができるのです。
Discussion