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)

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

0 件のコメント: