第五回: Pythonで素数を探す----エラトステネスのふるい

概要
素数を調べるエラトステネスの篩と呼ばれるアルゴリズムをPythonで実装する。それを利用して、双子素数を求めるプログラムを作る。

# Function of prime numbers, i.e., Sieve of Eratosthenes
def prime(n):
    l = list(range(2, n+1))
    i = 0
    while i < len(l):
        l = [j for j in l if not(j%l[i] == 0) or (j/l[i] == 1)]
        i+= 1
    else:
        return l

################
# Sieve of Eratosthenes: version 2
def prime_2(n):
    l = list(range(2, n + 1))
    l_empty = []
    while len(l) >= 1:
        l_empty.append(l[0])
        s = set(l) - {j for j in l if j%l[0] == 0}
        list_s = list(s)
        l = sorted(list_s)
    return l_empty, len(l_empty) 
   
################
# Example of prime numbers
print(prime(20))
[2, 3, 5, 7, 11, 13, 17, 19]
################

# Twin prime numbers
def twin(n):
    l = prime(n)
    l2 = []
    for i in l:
        j = i + 2
        if j in l:
            l2.append((i, j))
    return l2

################
# Example of twin prime numbers
print(twin(20))
[(3, 5), (5, 7), (11, 13), (17, 19)]        

エラトステネスの篩の数学的基礎

どのような整数が素数であるかどうかを調べたいとき、エラトステネスの篩(ふるい)(Sieve of Eratosthenes)と呼ばれるアルゴリズムが最も有名です。それは次のようなものです。例えば、10以下の素数を考えるとき、2から10までの整数の集まりを考えます。
(0)

(1)
2で割り切れる整数を から除く。つまり、
2, 3, 4, 5, 6, 7, 8, 9, 10
(2)
2の次に小さい整数(つまり3)で割り切れる整数を から除く。つまり、
2, 3, 4, 5, 6, 7, 8, 9, 10
(3)
これを最後まで繰り返し、残った数が素数である。つまり、
3の次に小さい数(つまり5)で割り切れる整数を から除く。
2, 3, 4, 5, 6, 7, 8, 9, 10
5の次に小さい数(つまり7)で割り切れる整数を から除く。
2, 3, 4, 5, 6, 7, 8, 9, 10
(4)
10以下の素数は となる。

エラトステネスの篩のプログラムの概要

もう少しPythonのプログラムに適した形に表現してみます。 とします。
(0)
リスト [ ] を生成する。
空リスト [ ] を定義する。
(1)
リスト の最初の数 [ ] を取る。つまり、リストの0番目の数を取る。ここでは、2を取る。その数を に入れる。つまり、 [ ]となる。
(2)
リスト の中で、リストの0番目の数(つまり )で割り切れる数を(集合として)集める。つまり、
%
(3)
リスト を集合にする:
集合 から集合 を引いた集合 を考える。ここで、 より、 は の補集合である。つまり、
よって、 である。
(4)
集合 をリスト化したものを新たに とする。
ここでは、 [ ]
(5)
これを の数 が0個になるまで繰り返す。つまり、リストの要素の数が1個まで上の操作を繰り返し、もし0個になったら を表示する。
[ ] [ ]
[ ] [ ]
[ ] [ ]
[ ] [ ]
[ ]

あとは、これをプログラムとして書けばいいです。その前にPythonのいくつかの知識をまとめます。

Pythonの基礎知識

このエラトステネスの篩のプログラムを作るために、以下では特にrangeとwhileとsetを説明します。

range()

と入れたとき、 から までのリスト [ ] を作ります。このためにrangeを使います。これは次のように使います。

>>> list(range(1, 10))
[1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> list(range(1, 10, 2))
[1, 3, 5, 7, 9]
>>> list(range(2, 11))
[2, 3, 4, 5, 6, 7, 8, 9, 10]

最初の例は1から9までのリストの生成であり、二つ目の例は1から9までの奇数のリストの生成です。要は1から始まり、+2の数をリストの中に入れます。今回はそのような機能を使いませんので、2番目の例はわからなくても構いません。
これらの例から容易にわかるように、 で を含むリストを作るためには としなければなりません。

while

次に繰り返しの操作をする を説明します。 の後に書かれる条件を満たすまで繰り返すという意味です。だいたいこんな感じです。

>>> i = 0                  # define i = 0
>>> while i < 4:           # i < 4
...     print(i)           # print i
...     i = i + 1          # add 1 to i. Instead it' s OK to write i+= 1.
... 
0           # print(0)
1           # print(0 + 1 == 1)
2           # print(1 + 1 == 2)
3           # print(2 + 1 == 3)
       # i = 3 + 1 == 4:  i is not smaller than 4. 

を使うときは無限ループに陥らないように注意しなければなりません。無限ループを回避するためのテクニックがありますが、今回は省略します。

# The infinite loop
>>> l = [2, 3, 4, 5, 6]
>>> while len(l) >= 1:
...     print('Yes')
... 
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
......
...... (infinite) # control + C
# loop
>>> l = [2, 3, 4, 5, 6]
>>> while len(l) >= 1:
...     print('Yes')
...     l.pop(0) # delete the zeroth number from the list l, i.e., l == [3, 4, 5, 6]
... 
Yes
2
Yes
3
Yes
4
Yes
5
Yes
6
# pop
>>> l = [2, 3, 4, 5, 6]
>>> l.pop(0)
2
>>> print(l)
[3, 4, 5, 6]

集合 set()

集合 は集合の形に変えます。集合は直接 で書いてもいいです。

# set
>>> set(range(1, 3))
{1, 2}
>>> {1, 2}
{1, 2}
# converse lists to sets
>>> l = [1, 2]
>>> set(l)
{1, 2}
# converse sets to lists
>>> s = {1, 2}
>>> list(s)
[1, 2]

数学で集合論を勉強したのと同じように、集合間には演算が定義できます。
はじめに共通部分 、和集合 、差集合と排他的OR(どちらか一方に含まれているが、両方には含まれていない要素の集合)です。

>>> s1 = {1, 2, 3, 4, 5}
>>> s2 = {2, 4, 6}
# intersection
>>> s1 & s2
{2, 4}
>>> s1.intersection(s2)
{2, 4}
# union
>>> s1 | s2
{1, 2, 3, 4, 5, 6}
>>> s1.union(s2)
{1, 2, 3, 4, 5, 6}
# difference
>>> s1 - s2
{1, 3, 5}
>>> s1.difference(s2)
{1, 3, 5}
# exclusive OR, i.e., symmetric difference
>>> s1 ^ s2
{1, 3, 5, 6}
>>> s1.symmetric_difference(s2)
{1, 3, 5, 6}

2つの集合の関係が部分集合かどうか、真部分集合か判定することができます。

>>> s1 = {1, 2, 3, 4, 5}
>>> s2 = {2, 4, 6}
>>> s3 = {1, 3, 5}
# whether a set is a subset of another set or not.
>>> s1 <= s1
True
>>> s1.issubset(s1)
True
>>> s2 <= s1
False
>>> s3 <= s1
True
# whether a set is a proper subset of another set or not.
>>> s1 < s1
False
>>> s2 < s1
False
>>> s3 < s1
True
# complement
>>> s3 <= s1
True
>>> s1 - s3
{2, 4}

今回は差集合しか使いません。だいたいこんな感じです。

>>> l = list(range(1, 11)) # l = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> {j for j in l if j%2 == 0 } # the set of numbers whose reminder divided by 2 is zero, i.e., the set of even numbers.
{2, 4, 6, 8, 10}
>>> {j for j in l if j%2 == 0 } <= set(l) # subset
True
>>> set(l) - {j for j in l if j%2 == 0} # difference 
{1, 3, 5, 7, 9}

最後に注意事項を書きます。どういうわけか、差集合を作ると、要素の順番が変わってしまうときがあります。集合自体は要素間の順序は関係ありませんが、リストでは重要ですので、集合をリスト化するときに、順番も で変えます。

>>> {3, 5, 7, 9, 11, 13} - {3, 9}
{13, 11, 5, 7}
>>> s = {3, 5, 7, 9, 11, 13} - {3, 9}
>>> list(s)
[13, 11, 5, 7]
>>> sorted(list(s))
[5, 7, 11, 13]

以上でエラトステネスの篩のプログラムを書くための準備ができました。それでは次にそれを書いてみます。

エラトステネスの篩のプログラム

それではエラトステネスの篩のプログラムを書きます(こちらは、Version 2です。最初のやつは友達に手伝ってもらいました。Ver. 2は自力で作りました)。

>>> # Sieve of Eratosthenes: Version.2
... n = 10
>>> l = list(range(2, n + 1)) # l == [2, 3, 4, 5, 6, 7, 8, 9, 10]
>>> l_empty = []
>>> while len(l) >= 1:
...     l_empty.append(l[0])
...     s = set(l) - {j for j in l if j%l[0] == 0}
...     list_s = list(s)
...     l = sorted(list_s)
... 
... 
>>> print(l_empty)
[2, 3, 5, 7]
>>> # Sieve of Eratosthenes
>>> n = 10
>>> l = list(range(2, n + 1))
>>> i = 0
>>> while i < len(l):
...     l = [j for j in l if not(j%l[i] == 0) or (j/l[i] == 1)]
...     i+= 1
... else:
...     print(l)
... 
[2, 3, 5, 7]

今回は集合の差集合を使い、リスト の要素を削除しましたが、リストの や などを使っても同様のものが作れるかもしれません。リストのみで作ろうかと試しましたが、よくわかりませんでした。

あとは、これを関数として定義すれば、素数のプログラムは完成です。リストで表示される代わりにタプル( )で表示したいのならば、 とすればいいです。

>>> # Sieve of Eratosthenes: version 2
>>> def prime_2(n):
...     l = list(range(2, n + 1))
...     l_empty = []
...     while len(l) >= 1:
...         l_empty.append(l[0])
...         s = set(l) - {j for j in l if j%l[0] == 0}
...         list_s = list(s)
...         l = sorted(list_s)
...     return l_empty 
... 
>>> prime_2(23)
[2, 3, 5, 7, 11, 13, 17, 19, 23] 
##########
# Function of prime numbers, i.e., Sieve of Eratosthenes
>>> def prime(n):
...     l = list(range(2, n+1))
...     i = 0
...     while i < len(l):
...         l = [j for j in l if not(j%l[i] == 0) or (j/l[i] == 1)]
...         i+= 1
...     else:
...         return l
... 
>>> prime(23)
[2, 3, 5, 7, 11, 13, 17, 19, 23]

どちらのプログラムの方が効率的なのかどうか気になりますが、どう判断すればいいのかわかりませんので、割愛します。
最後に、エラトステネスの篩のプログラムの応用として、双子素数を検出するプログラムを書いて終わりにします。

エラトステネスの篩のプログラムの応用: 双子素数の検出

素数を数えるプログラムができたので、それを応用してみます。双子素数(twin prime numbers) とは、(3, 5)のように二つの素数の差が2であるような素数の組のことを言います。例えば、20以下の整数ならば、双子素数は (3, 5), (5, 7), (11, 13), (17, 19)の4組です。これを求めるプログラムを作ってみます。
(0)
まず、 で20以下の素数のリストを求める。つまり、 [ ]
空リスト [ ] を用意する。
(1)
次にリスト の各要素 に対して、 とおく。
(2)
もし、 がリスト の中に入っているのならば、つまり、 がリストの要素であるならば、 を に入れる。
(3)
をプリントする。

これをプログラムに書くと次のようになります。

>>> # twin prime numbers up to 20
>>> # (0)
>>> l = prime(20) # l == [2, 3, 5, 7, 11, 13, 17, 19]
>>> l2 = []
>>> # (1)
>>> for i in l:
...     j = i + 2
...  # (2)
...     if j in l:
...         l2.append((i, j))
... 
>>> # (3)
>>> print(l2)
[(3, 5), (5, 7), (11, 13), (17, 19)]

したがって、双子素数を求める関数は次のようになります。

>>> #twin prime numbers
... def twin(n):
...     l = prime(n)
...     l2 = []
...     for i in l:
...         j = i + 2
...         if j in l:
...             l2.append((i, j))
...     return l2    
... 
>>> twin(20)
[(3, 5), (5, 7), (11, 13), (17, 19)]
>>> twin(1000)
[(3, 5), (5, 7), (11, 13), (17, 19), (29, 31), (41, 43), (59, 61), (71, 73), (101, 103), (107, 109), (137, 139), (149, 151), (179, 181), (191, 193), (197, 199), (227, 229), (239, 241), (269, 271), (281, 283), (311, 313), (347, 349), (419, 421), (431, 433), (461, 463), (521, 523), (569, 571), (599, 601), (617, 619), (641, 643), (659, 661), (809, 811), (821, 823), (827, 829), (857, 859), (881, 883)]

素数が無限に存在することはユークリッドの頃から知られていますが、双子素数が無限に存在するかどうかはまだわかっていないそうです。しかし、素数が存在する間隔についての理論はかなり展開されているそうです。いつか勉強してみたいです。

僕から以上