2009年5月28日木曜日

pi関数

TAOCPにpi関数というのが登場する. 真理値表が 1100 1001 0000 1111 となる4変数の関数である. これは, 種明かしがしてあるが, 円周率πを二進で表わしたものだ. ランダムな真理値表が欲しいとき, πや素数表を使う.

十進のπの値は良く知られている.

TAOCPの第1巻の終りの方に八進のπの値が書いてあり,

3.11037 55242 10264 30215 14230 63050 56006 70163 2112

これを二進にすれば

11.001 001 000 011 111 101 101 010 100 010 ...

だから, 始めの方はたしかに上の真理値表だ.

ところが512ビットの真理値表というのも現れる.

1100100100011...00101

だというのである. 本当かな?と思って計算してみた. ついでだから, 1024ビットまで求めて, 十六進で書いたのが下のものだ.


2行目の最後が512桁で, 45だから00101だった.

この十六進の文字の分布も計算した. 下のリストの左から0, 1, 2,... , fの出現数で, 合計は256になっている.

(apply + '(15 16 17 16 21 17 21 10 14 13 14 19 18 14 17 14))
=> 256

ところでこの計算法だが, Lispのbignumを活用する. 十進法のπの値はspigotアルゴリズムで計算した. もっともここの段階で十六進に出来るはずだが, 面倒だったので, 十進の値を出した. 1024ビットは十進では310桁くらいなので, 320桁を用意した.


小数点は下から319桁にあるから, 10の319乗を用意する.

(define pow10 (expt 10 319))

とする. また最初を1100にしたいから, piを4倍する.

(set! pi (* pi 4))

十六進の印字プログラムも必要

(define (hexprint n)
(define hex "0123456789abcdef")
(display (string-ref hex n)))

そして計算する.

(hexprint (quotient pi pow10))
(set! pi (modulo pi pow10))

ここで最初のcが出力される.

(do ((i 0 (+ i 1))) ((= i 256))
(hexprint (quotient (* pi 16) pow10))
(set! pi (modulo (* pi 16) pow10)))

ここから90fdaa...が得られた.

折角プログラムがあるのだから, 自然対数の底のeでも同じ計算をしよう.


TAOCPの八進の値は

2.55760 52130 50535 51246 52773 42542 00471 72363 6166

で, 十六進1024ビットは以下のとおりである.

2009年5月10日日曜日

Dijkstraのプログラム

2008年8月25日のこのブログに, Derrik Henry Lehmerの篩機械のことを書いた. 彼の父親の, Derrik Norman Lehmerも数学の教授であった. 1914年に10006721までの素数表を出版している. このように, どこそこまでの素数表を作るという場合には, Eratosthenesの篩は有用だが, 例えば最初の100個の素数表となると, どこまでの篩で始めればよいか, 問題である. もちろん x/log x=100からxを解けばよいが, 面倒でもある.

DijkstraのEWD227に, 最初の1000個の素数表を計算するプログラムがあった. 面白いので紹介したい. これはAlgol 60風の言語で書かれている. (a,b,c):=(1,2,3) という構文はAlgol 60にはないが, 気持はよく分かる. 説明用に行番号をつけた.


0 integer array p[1:1000];
1 begin integer k,j,inc,ord; p[1]:= 2; p[2]:= 3; k:= 3;
2 (j,inc,ord):= (1,4,1);
3 while k ≤ 1000 do
4 begin begin boolean jprime;
5 repeat (j,inc):= (j + inc, 6 - inc);
6 while p[ord]^2 ≤ j do ord:= ord + 1;
7 begin integer n;
8 n:= 3; jprime:= true;
9 while n < ord and jprime do
10 begin jprime:=(j != (j/p[n])*p[n]);
11 n:= n + 1
12 end
13 end
14 until jprime
15 end;
16 p[k]:= j; k:= k + 1
17 end
18 end


行0 素数の場所pを配列として1000個とる. 添字は1から始まる. C言語の配列は0始番だが, Algolは下限と上限を指定するので, Dijkstraは下限を1にしている. 私なら0にしたいところだ. この0で始めるか1で始めるかに関しては, TAOCPの演習問題2.2.2-19の解答に,
Even though it is almost always better for computer scientists to start counting at 0, the rest of the world will probably never change to 0-origin indexing. Even Edsger Dijkstra counts "1-2-3-4|1-2-3-4" when he plays the piano!
とある.

通常の素数表だと, 偶数は省略し, 奇数だけを順にテストするわけだが, 従って2だけ例外として, 特別に最初に表におくわけだが, このプログラムは一味違って, テストするのは6N±1である.(マイナスが先) N=1なら5と7, N=2なら11と13, N=3なら17と19,...従って2と3を例外として, 行1のp[1]:=2; p[2]:=3;と設定する. 6N±1を次々と作るには, incを交互に4,2,4,2とし, jの初期値1に順に足せばよい. それが行6で, j:=j+inc; inc:=6-inc; である. 値をa,b,a,b,..と変える常套手段はx:=aと初期化し, x=aを使ってから, x:=(a+b)-xとするのである.

次にjが素数であると分かるのは, jの平方根までの素数で順に割る. その上限の添字がordで, 行6でordを増やしている. それに続き, jを割ってみるわけだが, 2と3で割ってみる必要はないから, p[3]から割ってみれば良く, そのため行8でn:=3とする. jprimeというboolean変数は行4で宣言し, 行8でtrueにした. 行9のwhile文では, jprimeとnを更新する. jprimeは「jはp[n]で整除出来ない」と設定するので, 整除でき, jが合成数と分かれば, whileループから脱出する. n=ordでも脱出する. 行14のuntilはjprimeがfalseなら, つまりjが合成数なら, repeat文を繰り返し, 次の候補のjをテストする. 素数が見つかれば, untilの条件がtrueになり, repeat文から抜け, 素数表にそのjを登録するのである. kを増やし, while文の先頭の戻り, kが1000を越えていたら終る.

jprime(行4)やn(行7)のように, ローカル変数は必要最小限のところで宣言しているのが見所でもある.

ところで本年2月28日に, ShibuyaLispという会合があり, 参加してきた. 数理システムの黒田さんが, SchemeはS式で書いたAlgol 60であるといわれたのに, まったく同感であった.

さて, 上のAlgol 60もどきのプログラムを, Schemeで書いてみることにした. 1000個は遠慮して100個にした.

Goto Statements を harmful と喝破したDijkstraのプログラムだから, goto文は当然ないが, while文とrepeat until文が沢山ある. 通常ならこれを下請け関数にするのだが, 今回はwhileとrepeatのマクロ定義から始めた.

それぞれ以下のように定義し, 以下のように使う.

(define-syntax while
(syntax-rules ()
((while test command ...)
(letrec
((loop (lambda ()
(if test (begin command ... (loop))))))
(loop)))))

(define i 0)
(while (< i 10) (set! i (+ i 1)))
i=>10

(define-syntax repeat
(syntax-rules ()
((repeat command test)
(letrec
((loop (lambda ()
(begin command (if (not test) (loop))))))
(loop)))))

(define j 1)
(repeat (begin (set! j (* j 2)) (set! j (+ j 1))) (> j 20))
j=>31

この準備をした後,

(define (primes100)
(let ((p (make-vector 100)) (k 2) (j 1) (inc 4) (ord 0))
(vector-set! p 0 2) (vector-set! p 1 3)
(while (< k 100) (let ((jprime #t))
(repeat (begin (set! j (+ j inc)) (set! inc (- 6 inc))
(while (let ((pord (vector-ref p ord)))
(<= (* pord pord) j))
(set! ord (+ ord 1)))
(let ((n 2)) (set! jprime #t)
(while (and (< n ord) jprime)
(let ((pn (vector-ref p n)))
(set! jprime (not (= j (* (quotient j pn) pn))))
(set! n (+ n 1)))))) jprime)
(vector-set! p k j) (set! k (+ k 1))))
(display p)))

(primes100)

このプログラムは, pを0から取っているので, pの添字になるk, ord, nは1ずつ小さく初期化しなければならない.

結果は

#(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 101 103 107 109 113 127
131 137 139 149 151 157 163 167 173 179 181 191
193 197 199 211 223 227 229 233 239 241 251 257
263 269 271 277 281 283 293 307 311 313 317 331
337 347 349 353 359 367 373 379 383 389 397 401
409 419 421 431 433 439 443 449 457 461 463 467
479 487 491 499 503 509 521 523 541)

プログラムはやはり楽しい.

2009年5月9日土曜日

素数定理

2009年5月7日の素数定理のブログに, 青赤の線で素数分布の絵を描いた. それを描いた自分自身でも, その図で意外と思うのは, xが5万くらいのところで, どちらの線でも振幅が急に大きくなることだ. 素数の分布が変わるか? と一瞬思ったが, (まさか!) 青と赤の線で振幅の変わるxの位置が違うというのも怪しげなことである.

その原因を調べようと, 違う範囲の図を描いてみたら, またまた妙なことに, この振幅の変わる位置が, そのたびに違っていた. これはどうもPostscript(のシミュレータgs)に原因があるのではないか, というので理由を考えている. (まだ分かっていない.)

とりあえず, 描画の線幅を1から, 0.2に変えてみたら, お! 実にちゃんとした図が描けたのである. 0.5だと振幅の変化が相変わらず見られたが, うまくいった0.4で描いた図を下に示す. 縦の倍率は前回の5倍にしてある.

振幅の変化もなく, 赤い線が美しく1に接近しているのは感激だ.

さて, 傾向がよく見えるのはGaussの青線の方で, xが素数だとπ(x)が増えるので, 線は右上がりになり, 合成数だと増えないから, 右下がりになるのがよく分かる. つまり2から始まって, 3へは上昇, 4へは下降, 5で上昇, 6で下降, 7で上昇の後, 8,9,10で下降するのが見える. その後多少顕著なのは, x=100の先で, 青も赤も急上昇する. これは前回のQ0の左の方の816d, 二進では1000 0001 0110 1101.
このビットは右から97, 99, 101, 103, 105, 107, 109, 111, 113, 115, 117, 119, 121, 123, 125, 127に対応し, 赤字が素数で, 97から113で上昇が見られるのである.

素数分布で話題になるのは, 大きい方での双子の素数と, ギャップである. 偶数が邪魔するので, 2つの続く整数が素数のことはないが, 1つ置きになら存在でき, それがどの辺まで分かっているかとか, どの辺にずーっと素数のない区間があるかとか, 関心が持たれている. TAOCPのMMIXでの素数表の解答のところには, このようなプログラムを使い, 418032645936712127とその次の素数の間に, 1370のギャップを見つけたと書いてある.

2009年5月8日金曜日

素数定理

The Art of Computer Programming, Volume 4, Pre-Fascicle 1Aに奇数に対する素数表がある. (5ページ)



早くいえば, 前回の話題の最後の方の十六進のビットマップを二進で表示したものである.

Q0のビット列の右からの位置, i=0,1,...の各ビットは整数p=2i+1に対応し, pが素数なら1, 合成数なら0になっている. Q0の左端はi=63, p=127, これは27-1なので, 有名なMersenne素数であり, ビットも1になっている.

各Qは64ビット, 8語全体で512ビットなので, 最後(Q7の左端)はi=511, p=1023に相当する. 1023は, 210-1だが, 冪の10が合成数なので, 1023ももちろん合成数である. 1023=3*11*31

TAOCPは最近はMMIXというRISCマシンで記述されている. MMIXは1語最大64ビットである. Eratosthenesの篩を使った, 上の表の作り方を, 今日はC言語で書いて説明したい. C言語は, 私にとっては第一外国語みたいなもので, 大体はOKだが, 多少怪しいという言語だ. これはTAOCPの演習問題7.1.3-24であり, 同書にMMIXALによるプログラムも載っている.


#include <stdio.h>
main (){
unsigned long long int qs [8];
unsigned long long int
x33=0x9249249249249249llu|0x4008010020040080llu;
unsigned long long int
x35=0x8421084210842108llu|0x0408102040810204llu;
unsigned long long q;
int j,pp=1,p=13,m;
qs[0]=0x816d129a64b4cb6eull;
for(j=1;j<8;j++){
qs[j]=(x33|x35)^0xffffffffffffffffull;
x33=x33<<2|x33>>31;x35=x35<<6|x35>>29;}
for(j=0;j<8;j++)printf("%16llx\n",qs[j]);printf("\n");
q=qs[0]>>6;
do{
for(m=p*p/2;m<8*64;m=m+p){qs[m/64]=qs[m/64]&
(0x00000000000000001llu<<(m%64)
^0xffffffffffffffffull);}
q=q>>1;p=p+2;
while(q==0){q=qs[pp];pp=pp+1;p=(p|0x7f)+2;}
while((q&1)==0){;q=q>>1;p=p+2;}
}while(p<1024);
for(j=0;j<8;j++)printf("%16llx\n",qs[j]);}
}

14行目と24行目で, qsを十六進で出力している.

816d129a64b4cb6e
2996c20d865a4c32
b49961225a2534c9
4a6982d12d86134c
c36996132ca45b0
948a6894d325264b
4b48b6186d30d225
65a4c92916d128a6

816d129a64b4cb6e
2196820d864a4c32
a48961205a0434c9
4a2882d129861144
834992132424030
148a48844225064b
b40b4086c304205
65048928125108a0

この後半の十六進表記を二進にしたのが, 最初の表である.

qs[0]は127までの素数表で, これは手作業で作る. 100までの素数は25個しかなく, 127は13*13より小さいから, なんでもない. qs[1]からqs[7]までは, 篩だから, 通常はすべて1にするが, このプログラムでは, 3か5か7か11かで割れる数(に対応するビット)は, 最初から0にしておき, 13から篩い始める.

初期パターンを作るため, x33とx35を用意する. x33のビット毎ORの左は, 129から始めて3で割れる位置にビットをたててある. 129は3で割れるから, 右端は1, そのあと左へ向って001001と続ければ良い. 同じ語の右は11で割れるもののパターンで, 121が割りきれるから, 次に割りきれる奇数は143になる. これは129を0ビット目とすると, 7ビット目なので, そこから左へ1ビットめ毎に1をたてる. x35は同様なことを5と7についてやっている.

nで割れるパターンは長さnで繰り返す. 従って3と11で割れるパターンは33で繰り返し, 5と7で割れるパターンは35で繰り返す. MMIXは1語64ビットなので, 3と11, 5と7の繰り返しのパターンは1語に収まり, つまり合成できる. このことは私がKnuthに提案したので, TAOCPの解答には, We can combine 3 with 11 and 5 with 7, as suggested by E. Wada. と書いてある.

これらのパターンで1が立っているところは, 初期パターンでは0にするので, ビットの否定など, 多少注意が必要だ. またx33もx35も, 次の初期パターンに対しては, ローテートする必要がある.

とにかくこうしてqs[1]からqs[7]が設定出来た. 実行結果の上半分は, この設定直後のqsを示す.

次は篩だ. qにqs[0]を入れ, 11まで済んだので, 6ビット右シフとし, pも13に設定する. そしてqs[1]から篩始める. 最初の場所はp*p/2から. pの場所は素数を示すから, 消してはいけない. 次のpの倍数で奇数の場所はpだけ先である. しかし, 例えば5で篩う場合, 15はすでに3で篩われているから, 25から篩うので充分である. 25はi=25/2のビットに対応する. (小数部は棄てる.)

ある素数で篩ったら, qを右シフトし, pを2増やしながら, 次のビット位置を調べる. そうこうしている内に, qが0になったらどうするか. qを64ビットシフトしたか確認したくなるのが人情である. このプログラムではそういうことはせず, qsの次の語を取ってくる. 前のqsの語が終ったとき, ビット位置は64の整数倍引く1で, pは128の整数倍引く1の筈である. 従ってpは0x7fをビット毎ORして次のqsに対応する値にしている. この技もちょっとしたもである.

MMIXALのプログラムはもう少し凝っているが, 趣旨は殆んど同じである.

2009年5月7日木曜日

素数定理

高木貞治先生の近世数学史談に, Abelの手紙の日付が3√6064321219 とあり, それが1823年8月4日と解されている.

計算してみると, この数の立方根は
(expt 6064321219 1/3)=>1823.5908275197114
であり, この小数部分が本当に8月4日であろうかという気分になる. そこでユリウス日を計算する関数jdを使い,

(define (abel m d)
(exact->inexact
(/ (- (jd 1823 m d) (jd 1823 1 1))
(- (jd 1824 1 1) (jd 1823 1 1)))))

(abel 8 4)=> .589041095890411
(abel 8 5)=> .5917808219178082

となり, 8月4日を確認した.

あるいは

(* .5908275197114 365) => 215.65204469466102

1月から3月は90日, 4月から6月は91日, 7月は31日だから, そこまでで212日. つまり212.なにがしなら8月1日であり, 215.なにがしは8月4日である.

実は日付はどうでもよく, 問題はその手紙にAbelがLegendreの素数定理の式を書いていることである.

π(x)=x/(log x - 1.08366) (a)

昔読んだ素数定理は, Gaussが考えたもので,

π(x)=x/log x (b)

という気持のいいものであるが, それはあまりよく合わないといわれている.

早速絵を描いてみる.

x=2から1000000まで, 2つの式で左辺/右辺をプロットした. 青が(b), 赤が(a)である.999999点を描いているのだが, 最近の計算機ではあっという間である.

Legendreのは流石によくあっている. 高尚な理論あるようだが, 最小自乗法で合わせたかと思うほどである.

π(x)は, Eratosthenesの篩を使い, 1000000までの素数の奇数のビットマップを作った. MITのSchemeにはbitstringという便利なものがある.

(define x 500000) (define ps '(2)) ;xは上限 psは素数の表
(define sieve (make-bit-string x #t)) ;篩を作る
(bit-string-clear! sieve 0) ;素数でない1の場所をクリアする
(do ((i 0 (+ i 1)) (c 1)) ((>= i x))
(if (bit-string-ref sieve i) ;iの場所が1なら2i+1が素数
(let ((p (+ i i 1))) ;pの値
(set! c (+ c 1)) (set! ps (cons p ps)) ;c,psを更新
(do ((q (/ (- (* p p) 1) 2) (+ q p))) ((>= q x))
(bit-string-clear! sieve q)))));(p*p-1)/2からp置きにクリア

(define port (open-output-file "bittable"))
(do ((i 0 (+ i 128))) ((>= i x));128bitごとに符号なし整数に
(display (substring (number->string (+
(bit-string->unsigned-integer (bit-substring sieve i
(min x (+ i 128)))) (expt 2 128)) 16) 1 33) port)
(newline port)) ;十六進で出力
(close-output-port port)

これで3907行の出力が得られる. 最初の8行は

2196820D864A4C32816D129A64B4CB6E
4A2882D129861144A48961205A0434C9
148A48844225064B0834992132424030
65048928125108A00B40B4086C304205
C02104C94112422180124496804C3098
220825B0826896810804490000982D32
69009244304340069004265940A28948
086122D22400C06012410DA408088210

先頭の語の最後の 6E は 01101110で, 右から3,5,7,11,13が素数. 1,9,15は合成数を示す.