エラトステネスのふるいとは素数を求めるアルゴリズムです。素数とは 1 とその数自身以外で割り切る事ができない 1 より大きい整数のことです。
エラトステネスのふるいでは次のような手順で素数を求めます。
たとえば 50 未満の整数に対して、このアルゴリズムを視覚的に観察するプログラムを書いてみましょう。どのような感じで、整数がふるいにかけられるのかを眺めてみる事にします。
50 個の要素(その値は True か False)をもつリストを作成して、その要素が True であれば素数である。False であれば素数ではないという風にします。たとえば、リストの n 番目の要素がもし True であれば n は素数です。
以下の printn 関数は引数のリストの要素が True であれば、そのインデックスを表示します。
def printn(lst):
size = len(lst)
for i in range(0, size, 10):
for j in range(10):
n = i+j
if n < size and lst[n]:
print "%3d" % n,
else:
print " ",
print>>> printn([True] * 50) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
リストに * 演算子を適用すると、左オペランドのリストが右オペランドの値(整数)の回数だけ連結した新しいリストを返します。
>>> [10] * 5 [10, 10, 10, 10, 10] >>> [1, 2, 3] * 3 [1, 2, 3, 1, 2, 3, 1, 2, 3]
range 関数を引数 1 つで呼び出すと、0 から始まって引数から 1 引いた数までの範囲からなるリストを返します。
>>> range(3) [0, 1, 2] >>> range(10) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
range 関数に引数を 2 つ渡すと、最初の引数がリストの開始値となります。
>>> range(1, 3) [1, 2] >>> range(5, 10) [5, 6, 7, 8, 9]
range 関数に引数を 3 つ渡すと、3 番目の引数が刻み幅になります。
>>> range(1, 10, 2) [1, 3, 5, 7, 9] >>> range(3, 10, 5) [3, 8] >>> range(0, 50, 10) [0, 10, 20, 30, 40] >>> range(10, 0, -1) [10, 9, 8, 7, 6, 5, 4, 3, 2, 1]
print 文はカンマで終了すると、最後に改行を行いません。
>>> def foo():
print "Hello",
print "World"
>>> foo()
Hello World
>>> def bar():
... print "Hello", "World",
... print "Python"
>>> bar()
Hello World Python
文字列に対して % 演算子を適用すると、文字列の整形が行われます。文字列の中の変換指令(たとえば %d や %sなど)の所に、右オペランドのオブジェクトが文字列として埋め込まれ、新たな文字列が返されます。
>>> "%s" % 100 '100' >>> "%5d" % 100 ' 100' >>> "%-5d" % 100 '100 ' >>> "total = %d" % 120 'total = 120' >>> "%d" % "foo" TypeError: int argument required
%5d とすると、整数が文字列に変換されたとき 5 文字に満たなければ、残りを空白で埋めてくれます(%-5d とすると左詰になります)。
%d はオブジェクトとして int(整数) を要求します。ので、文字列を %d に適用する事は出来ません。文字列に対する変換指令は %s です。
>>> "Hello, %s" % "World" 'Hello, World'
しかし %s は文字列オブジェクトに限らず任意のオブジェクトを受け取る汎用変換指令として機能します。
>>> "%s" % 100 '100' >>> "%s" % [1, 2] '[1, 2]' >>> "%s" % True 'True' >>> "%s" % None 'None'
文字列の中に複数の変換指令がある場合は、右オペランドをそれと同数の要素をもつタプルにします。
>>> "%d %d" % (1, 2)
'1 2'
>>> "%s %s" % ("Hello", "World")
'Hello World'
また、文字列を囲むクオートですが、これはシングルクオート、ダブルクオートのどちらで囲ってもかまいません。2 つの書き方があることで、文字列中のクオートをエスケープする必要がなくなります。
>>> "He said to me, \"Are you hungry?\"" 'He said to me, "Are you hungry?"' >>> 'He said to me, "Are you hungry?"' 'He said to me, "Are you hungry?"' >>> "I'm very happy to meet you." "I'm very happy to meet you." >>> 'I\'m very happy to meet you.' "I'm very happy to meet you."
printn 関数ができたので、あとはエラトステネスのふるいで、整数がふるいにかけられる様子を観察してみましょう。
def sieve(size):
primes = [True] * size
primes[0] = False # 0 は素数でない
primes[1] = False # 1 は素数でない
printn(primes)
for i in range(2, size):
if primes[i]:
print "%d の倍数を削除します(%d を除く)" % (i, i)
for j in range(i+i, size, i):
primes[j] = False
printn(primes)>>> sieve(50)
2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18 19
20 21 22 23 24 25 26 27 28 29
30 31 32 33 34 35 36 37 38 39
40 41 42 43 44 45 46 47 48 49
2 の倍数を削除します(2 を除く)
2 3 5 7 9
11 13 15 17 19
21 23 25 27 29
31 33 35 37 39
41 43 45 47 49
3 の倍数を削除します(3 を除く)
2 3 5 7
11 13 17 19
23 25 29
31 35 37
41 43 47 49
5 の倍数を削除します(5 を除く)
2 3 5 7
11 13 17 19
23 29
31 37
41 43 47 49
……途中省略……
47 の倍数を削除します(47 を除く)
2 3 5 7
11 13 17 19
23 29
31 37
41 43 47エラトステネスのふるいでは、
というごく当たり前のことを述べています。2 以上の整数に対して、そこから(2 より大きい) 2 の倍数をふるいにかけ、(3 より大きい) 3 の倍数をふるいにかけ…、そして最終的に残ったものが素数である、というようにして素数を求めています。
エラトステネスのふるいは、ある数 n が素数なら、n の倍数は素数でない、ということで素数を求めていました。このアルゴリズムとは別に、ある数 n を実際にそれよりも小さな値で割ってみるという方法もあります。もし n を割り切る数がないなら、n は素数です。実際には n 以下の素数で割り切れるか調べればよいですね。
def isprime(n, lst):
for i in lst:
if n % i == 0:
return False # 割り切れたら素数でない
return True # 素数表のどの素数でも割り切れないなら素数である
def primes(n):
prms = [2] # 素数表に 2 を入れておく
for i in range(3, n+1, 2):
if isprime(i, prms):
prms.append(i) # 素数を素数表に追加
return prms>>> primes(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]
2 より大きい偶数は素数ではないことは予め分かっているので、3 から値を 2 ずつ増やし、奇数に対してのみ素数かどうかを調べるようにしています。
isprime 関数は、引数の(素数かどうかを調べたい) n に対して、それまでに見つかった全ての素数で割り切れるかどうかを調べています。しかし、実際にはルート n 以下の素数で割り切れるかを調べるだけで構いません。
たとえば、n が素数でないなら、n は a × b の形で表すことが出来ます。このとき、 a または b のどちらかはルート n 以下でなくてはなりません。もしそうでなければ、a × b は n よりも大きい値ということになっていまいます。
何故そうなるかを図形的に考えてみましょう。n が a と b の積で表されるとき、a と b を四角形の縦と横の長さだと考えます。そうすると n は四角形の面積となります。
縦の横の長さで短い方(図では a)が、もしルート n よりも大きいとします。n はルート n の二乗ですので、
a^2 > n (なぜなら、 a > ルートn)
ということになります(a^2 は a の二乗)。これはおかしいですね。
したがって、n が素数であるかどうかを調べるには、ルートn までの値を調べれば良いことになります。ルート n 以下の素数で割り切ることができなければ、n は素数です。
先ほどのプログラムを改良します。
import math
def isprime2(n, primes):
for i in range(int(math.sqrt(n))):
if n % primes[i] == 0:
return False
return True
def primes2(n):
prms = [2]
for i in range(3, n+1, 2):
if isprime2(i, prms):
prms.append(i)
return prms>>> print primes(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]
sqrt は平方根を求める関数です(import についてはモジュールの所で説明を行う予定です)。
ある数 n が素数かどうかを調べるときに、ルート n までの素数を調べればよいというだけの変更ですが、2 つのプログラムでどの程度差が出るのか 10000 までの素数を求めて比べて見ます。
5002 function calls in 0.649 CPU seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.649 0.649 <string>:1(?)
1 0.000 0.000 0.649 0.649 profile:0(primes(10000))
0 0.000 0.000 profile:0(profiler)
1 0.053 0.053 0.649 0.649 sieve.py:10(primes)
4999 0.596 0.000 0.596 0.000 sieve.py:4(isprime)
5002 function calls in 0.246 CPU seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.246 0.246 <string>:1(?)
1 0.001 0.001 0.246 0.246 profile:0(primes2(10000))
0 0.000 0.000 profile:0(profiler)
4999 0.191 0.000 0.191 0.000 sieve.py:17(isprime2)
1 0.055 0.055 0.246 0.246 sieve.py:23(primes2)輪郭分析 で作成した行実行回数プロファイラでも試してみます。
import math
4999 def isprime(n, lst):
771632 for i in lst:
3771 if n % i == 0:
return False
return True
1 def primes(n):
prms = [2]
4999 for i in range(3, n+1, 2):
1228 if isprime(i, prms):
prms.append(i)
return prms
4999 def isprime2(n, primes):
94496 for i in range(int(math.sqrt(n))):
3771 if n % primes[i] == 0:
return False
return True
1 def primes2(n):
prms = [2]
4999 for i in range(3, n+1, 2):
1228 if isprime2(i, prms):
prms.append(i)
return prms
primes(10000)
primes2(10000)大きな数を素因数分解することは大変難しい問題のようです。ですが、10 桁程度の数であれば、上記の「小さい素数で順に割っていく」というアルゴリズムでも十分間に合います。次の数は、2 つの素数の積になっています。それぞれ素因数分解してください(解答)。