前回のあらすじ
前回のブログでは、The PrimePagesで公開されている素数のデータベースをダウンロードし、これを読み込む部分を(pythonを用いて）書きあげるところまで完了した。
今回は、いよいよエラトステネスの篩と呼ばれる、古来より知られる有名なアルゴリズムを用いて、最初の10000個(1万個）の素数を計算してみようと思う。(計算結果は、The PrimePagesからダウンロードしたデータベースと比較し、正しく計算できたか確認する計画。）
エラトステネスの篩
高校数学でやったという記憶はないが、最近の高校の教科書では紹介されているようである。 うろ覚えながら、その手順は次のような感じとなる。
- 0は除く。
- 1は除く。
- 素数を見つけたい範囲を決める（たとえば0から100までの中に含まれる素数、といった具合）
- 2が最初の素数。
- 最初の素数（2）の倍数を、素数を見つけたい範囲から消す（2,4,6,8,...100、つまり100以下の偶数を消す）。
- 4.のプロセスで消去されなかった最初の数が、次の素数となる（3）。
- 2番目の素数（3）の倍数を、素数を見つけたい範囲から消す（3,6,9,.....,99)。
- ここまでで生き延びた最小の整数を3番目の素数とする（5）。
- 3-5のプロセスを繰り返し、最大の素数が100を超えたらアルゴリズムを終了する(97の次の素数が101となるので、最大の素数は97として計算終了）。
3-5のプロセスがメインの繰り返しジョブとなり、8の終了判定が重要である。
このアルゴリズムをより効率よく書き換えて、数値的に計算が楽になる有名なコードもおそらく存在するであろう。しかし、今回はまず自分でプログラムを組んでみることにする。
pythonによる実装
pythonはどの程度まで大きな整数を扱うことができるのであろうか？FORTRANやCではその最大値が決まっていて、gccの場合はlimits.hで調べることができるという。
実際にコードを書いて確かめてみた。
#include <stdio.h> #include <limits.h> int main(){ printf("Min for signed int: %d\n", INT_MIN); printf("Max for signed int: %d\n", INT_MAX); printf("Max for unsigned int: %u\n", UINT_MAX); }
実行結果は次のようになった。
Min for signed int: -2147483648 Max for signed int: 2147483647 Max for unsigned int: 4294967295
素数の計算では負数は出てこないから、unsigned intを用いることができる。となると、扱える最大の整数は4,294,967,295、つまり43億程度まで扱えるということである（裏金の帳簿管理プログラムにはギリギリ使える値だが、国家予算の帳簿プログラムには使えないかも）。この値は32ビットの最大値に対応しており、bcで確かみるとその通りである。
$ echo '2^32' | bc -ql > 4294967296
符号付きだと1ビット持って行かれてしまうので、31ビットとなり、上の値の半分となる。
$ echo '2^31' | bc -ql > 2147483648
ここで実験をやってみる。最大の整数に対して、あえて1を足してみるのである。水をなみなみと注いだグラスに、一滴水を垂らすと何が起きるであろうか？たぶん、これはマシンによって結果は異なると思うが。
#include <stdio.h> #include <limits.h> int main(){ int n; unsigned int m; n = INT_MAX; m = UINT_MAX; printf("Out of bounds: %d %d\n",n,n+1); printf("Out of bounds: %u %u\n",m,m+1); }
ここでのポイントは、符号なし整数の識別子は%uと表す点である（今日まで知らなかった....）。結果は次のとおり。
Out of bounds: 2147483647 -2147483648 Out of bounds: 4294967295 0
どうやら、丸い地球を平らだと思って探検に出かけ、水平線を越えてみたら出発点に戻ってしまった、という感じである。
FORTRANでは整数を使った数値計算はあまりやらないのだが、たぶんC言語と同じような状況になっているはずである。特にgfortranはgccの付加体として定義されているからなおのことである。（もしかするとgfortranなら#include <limits.h>を受け付けるかもしれないと思って、やってみたがダメであった。）
いろいろとFortranについての資料を探してみたが、整数型についての説明はあまりなかった。これは当然で、科学技術計算では浮動点小数計算が主流であり、あまり整数を使った計算はやらないからだ。とはいえ、最近では素数を利用した暗号技術などが重要になってきたので、Fortranで色々計算している人も多少はいるのではないかと想像する。また、デジタル情報の処理は基本的に整数でできるから、工学系の人で整数を扱うFortranプログラミングをしている人もいるのかもしれない。とはいえ、調べた限りでは、手持ちのgfortranの整数がどの程度の範囲まで扱えるのか調べきることができなかった。
こういう時は、自分で試してみるしかない。
Fortran95で整数の最大値を確かめる
たとえば、整数の4ビット表現の場合、ビットが4つあるので \begin{equation} [a][b][c][d] \end{equation} となる。ビット表現は2進数なので、a,b,c,dは1か0の値しかとれない。たとえば、 $[1][0][0][1]$といったような感じになる。これを10進数に変換する方法は（高校数学でも解説があるように） \begin{equation} 1\cdot 2^3 + 0\cdot 2^2 + 0\cdot 2^1 + 1\cdot 2^0 \end{equation} つまり、 $8+0+0+1 =9$である。したがって、4ビットの一般的な表現を10進数に変換する公式は \begin{equation} [a][b][c][d] = a\cdot 2^3 + b\cdot 2^2 + c\cdot 2^1 + d\cdot 2^0 \end{equation} である。
4ビットの場合は、2の3乗、つまり（二進数の）3桁までが登場する。したがって、4ビット表現で表せる最大の整数は $8+4+2+1=15$となる。これは$2^4-1$とも書けるが、よくよく考えると2の等比級数(n次まで）になっているから当然である。 \begin{equation} S_n = 1+2+2^2+\cdots +2^n = \frac{1-2^{n+1}}{1-2} \end{equation} しかし、この「当然」はプログラミングにおいては「当然」ではない。なぜなら、もし4ビット表現が採用されていた場合、$2^4$を計算した瞬間にオーバーフローが発生し、そこから1を引いたとしても正しい値に戻ってきてくれるかどうかは保証されないからである。したがって、上の等比級数の公式を「鵜呑み状態」で利用するのは危険である。より安全なのは、等比級数の定義式、つまり$2^3$までの等比数列を文字通り足しあげるアプローチである。
数値実験のためのプログラムは次のようなものである。
program max_int_limit implicit none integer :: n integer :: r,sum !integer*4 :: r,sum !integer*8 :: r,sum !integer*16 :: r,sum sum = 0 r = 1 do n=1,35 sum = sum + r print *,n,r,sum, sum + 1 r = r*2 enddo end program max_int_limit
integerの宣言にはさまざまなタイプがあるようだ。デフォルトのもの、$*4, *8, *16$などがある。C言語にある「符号なし」の宣言は見つからなかった。
プログラム中の変数の意味は、$n$が次数を表し（ビット数に相当）、$r$が等比数列、そして$sum$が等比級数である。
プログラム後半のループでは、$n$を一つずつ増加させながら、整数が正しく計算され表現できているかを確かめていく。破綻した表現が出た場所が、Fortran95が扱う整数の限界（最大値）となる。最大値がどこかを確かめれば、何ビット表現になっているかわかる、という狙いである。
さて、gfortranでコンパイルし、実行してみる。デフォルト、4, 8, *16の順番でやってみる。
デフォルトのinteger :: r, sumの場合
まずは$r$と$sum$を、デフォルトのintegerとして宣言して計算した場合を見てみる。結果は次のようになった。 (出力は、$n, r_n, S_n, S_n+1$、つまり次数、等比数列、級数、級数に1を足したもの、の順番）
1 1 1 2
2 2 3 4
3 4 7 8
4 8 15 16
5 16 31 32
6 32 63 64
7 64 127 128
8 128 255 256
9 256 511 512
10 512 1023 1024
11 1024 2047 2048
12 2048 4095 4096
13 4096 8191 8192
14 8192 16383 16384
15 16384 32767 32768
16 32768 65535 65536
17 65536 131071 131072
18 131072 262143 262144
19 262144 524287 524288
20 524288 1048575 1048576
21 1048576 2097151 2097152
22 2097152 4194303 4194304
23 4194304 8388607 8388608
24 8388608 16777215 16777216
25 16777216 33554431 33554432
26 33554432 67108863 67108864
27 67108864 134217727 134217728
28 134217728 268435455 268435456
29 268435456 536870911 536870912
30 536870912 1073741823 1073741824
31 1073741824 2147483647 -2147483648
32 -2147483648 -1 0
33 0 -1 0
34 0 -1 0
35 0 -1 0
$n=30$までは順調に計算は進んでいるが、$n=31$のとき、sum+1の値がおかしくなっている。つまり、gfortranのintegerは、32ビット表現になっていて、扱える最大の整数値は$2147483647 = 1+2+4+\cdots 2^{31}$というわけである。
$n=32$ではすでにオーバーフローしてしまい、$r$自体の値が変になっていて、級数（和）を計算してももはや何の意味もないことが確認できる。
integer*4 の場合
次にinteger4を試してみると、上の計算結果と同じとなった。つまり、デフォルトのintegerはinteger4と同じであり、この２つは32ビット表現になっているのである。
integer*8の場合
ではinteger*8を試してみよう。おそらく64bit表現になっているはずだから、nのループの最大値を65にしてみる。
program max_int_limit implicit none integer :: n integer*8 :: r,sum .... do n=1,65 ... enddo end program max_int_limit
結果はこうなった（前半部分省略）
56 36028797018963968 72057594037927935 72057594037927936
57 72057594037927936 144115188075855871 144115188075855872
58 144115188075855872 288230376151711743 288230376151711744
59 288230376151711744 576460752303423487 576460752303423488
60 576460752303423488 1152921504606846975 1152921504606846976
61 1152921504606846976 2305843009213693951 2305843009213693952
62 2305843009213693952 4611686018427387903 4611686018427387904
63 4611686018427387904 9223372036854775807 -9223372036854775808
64 -9223372036854775808 -1 0
65 0 -1 0
予想通りである。$n=63$の$sum$までは計算できているが、そこに1を足すと変なこと（つまり負数になってしまう）が発生し、その後の計算は「もうダメ」状態である。つまり、integer*8は64bit表現である。「900京」をこえる整数が扱えるというのは、すごいことだ！
この上となると「垓」の世界である（アボガドロ数の世界へ近づくということでもある！）。integer*16でコンパイルエラーが出なければ、とんでもなく大きさな整数をgfortranは扱えるということになりそうである。果たしてどうだろうか？
integer*16にトライ！
予想ではinteger*16は128bit表現になるはずである。$2^{128}=10^x$とおいて$x$を計算すると \begin{equation} x = 128\log_{10}2 \sim 128\cdot 0.3010 = 38.53... \end{equation} つまり$10^{38}$の世界である！アボガドロ数のオーダーが$10^{23}$であったから、それを遥かに上回る「宇宙の果て」の世界に似た感じがする大きさの整数である。
実行結果は次のようになって、ちゃんと「宇宙の果て」の雰囲気を醸し出す巨大な整数が表示されたのであった！
まちがいなく、128bitである！
ということで、どうやらFortran95で整数の研究をやっている数学者／情報工学者がいそうな感じに見えてきた。しかし、暗号資産とか電子通貨といったものはこれを超える桁での素数を研究しているはずだ。
pythonではどうなのか？
大きな素数の研究をするならFortran95でやった方がよいのだろうか？といった気持ちが芽生えてきたわけだが、果たして我々が利用しているPythonはどうなっているのであろうか？
ググってみると、どうやら「限界なし」というおそろしいことになっているらしいのである。
この辺りは次回にもう少し詳しく見てみたいと思う。有用そうな資料としては次のブログも注目である。
エラトステネスの篩のコードまでいけなかった....
コードはできているのだが、最大整数の扱いのところで知らなかったことが山ほど出てきてしまい、今回はそこで一杯一杯となってしまった。次回でPythonでの最大整数の扱いについて議論した後、その次でエラトステネスへいけるのではないかと思う。
