素数生成プログラム「エラトステネスの篩」を実装｜Pythonで数学を学ぼう！第7回

今回はPythonでエラトステネスの篩（ふるい）をプログラミングしてみましょう。

素数生成プログラム

エラトステネスの篩は素数を見つける方法の1つです。

次がエラトステネスの篩によって素数を生成するプログラムです。

>>> n =int(input("n以下の素数と個数を表示。n="))
>>> p = [i for i in range(n + 1)]

>>> for i in p[3:]:
>>> if p[i] % 2 == 0:
>>> p[i] = 0

>>> root_n = n ** 0.5
>>> for i in range(3, n):
>>>      if i > root_n:
>>>           break
>>>      if p[i] != 0:
>>>           for j in range(i, n + 1, 2):
>>>                if i * j >= n + 1 :
>>>                     break
>>>                p[i * j] = 0

>>> plist = sorted(list(set(p)))[2:]
>>> plistn = str(len(plist))+'個'

>>> print(plist)
>>> print(plistn)

実行するとnの値をたずねられるので100と代入してみると、次のように表示されます。

n以下の素数と個数を表示。n=100
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
25個

100以下の素数が小さい順に表示されます。最後にその個数が25個と表示されます。

エラトステネスの篩（ふるい）

エラトステネスの篩のアルゴリズムをみていきましょう。

100以下の素数を探すことにします。

ステップ1

1番目の素数2に対して、2より大きな2の倍数4,6,8,10,12,…,100を消します。

ステップ2

1番目の素数2より次に大きな数を残った数の中から探します。それが2番目の素数3です。

残った数の中から3より大きな3の倍数である9,15,21,27,33,39,45,51,57,63,69,75,81,87,93,99を消します。

ステップ3

2番目の素数3より次に大きな数を残った数の中から探します。それが3番目の素数5です。

残った数の中から5より大きな5の倍数である25,35,55,65,85,95を消します。

ステップ4

3番目の素数5より次に大きな数を残った数の中から探します。それが4番目の素数7です。

残った数の中から7より大きな7の倍数である49,77,91を消します。

ステップ5

4番目の素数7より次に大きな数を残った数の中から探します。それが5番目の素数11です。

ここで、11は√100=10よりも大きくなったのでSTOPします。ここで残った数が素数となります。

なぜ100の平方根(>0)よりも大きな素数が見つかったらSTOPするのか、わかりやすく小さな数36で説明してみます。

36をi×jに因数分解してみます。
i≦j
2×18
3×12
4×9
6×6

すると、最後の6が36の平方根だと分かります。

iを6より大きな数にしてしまうと、それは既に計算されたjとなります。

ですから36までの素数をエラトステネスの篩で求める場合には、素数7が見つかった場合にSTOPしていいのです。

プログラムの解説

>>> p = [i for i in range(n + 1)]
により、
p[0]=0
p[1]=1
p[2]=2
…
p[100]=100
を生成します。

>>> for i in p[3:]:
>>> if p[i] % 2 == 0:
>>> p[i] = 0
p[3],p[4],…,p[100]の中で2の倍数のものを0に置き換えます。

これが先のステップ1に相当します。ステップ1での「消す」代わりに「0に置き換える」ことにします。

>>> root_n = n ** 0.5
>>> for i in range(3, n):
>>>      if i > root_n:
>>>           break
>>>       if p[i] != 0:
>>>           for j in range(i, n + 1, 2):
>>>                if i * j >= n + 1 :
>>>                     break
>>>                p[i * j] = 0

先に説明したようにi×jのiは√nより大きくなったらSTOPします。

それまでの間、0でないp[i]=iの倍数を「消す」代わりに「0に置き換える」ことを繰り返します。

>>> plist = sorted(list(set(p)))[2:]
これがこのプログラムの肝になるコードです。

ここまでに、
p[0]=0, p[1]=1, p[2]=2, p[3]=3, p[4]=0, p[5]=5, p[6]=0, p[7]=7, p[8]=0, p[9]=0, p[10]=0, p[11]=11,…, p[99]=0, p[100]=0
となっています。

そこで、set(p)とすることでp[i]を要素とする集合（set）を生成します。

このとき、前回説明したように重複する要素は自動的に削除され1つになります。

p[0]=0, p[4]=0, p[5]=5, p[6]=0, p[9]=0, p[10]=0, …, p[99]=0, p[100]=0
といった0が重複しているので、これらはまとめて1つの要素0となります。

はたして、
set(p)={11,29,2,61,7,0,…,3}
のように、集合の中身は大小バラバラの順に整理されます。

そこで、まずこれをlist()によりリスト化します。
list(set(p)) =[11,29,2,61,7,0,…,3]

その上で、sorted()することで大きさの順に並べ変えます。すると、
sorted(list(set(p)))=[0,1,2,3,5,7,11,…,97]
となります。

最後に、最初の2つの0（0番目）と1（1番目）が不必要なので、2番目の2から表示させます。
sorted(list(set(p)))[2:] =[2,3,5,7,11,…,97]

この最後の処理過程をわかるようにしたのが次のプログラムです。