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)

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

Dijkstraのプログラム

オランダの計算機科学者 Edsger W. Dijkstra は, ものを書くのが好きで, BachのBWVにならってEWD何番というドキュメントが沢山ある. (http://www.cs.utexas.edu/~EWD/ )

それを眺めていたら, EWD221に Raw code for computing DeBruijn-sequence というのがあった.

例えば00011101がB(2,3)で, 最初のパラメータ2は, 記号が0と1の2種類ということ, 後のパラメータ3は, 3個ずつとると, 全てのパターンが得られるということである.

つまり先頭から1ビットずつずらしながら, 3ビットとると

000
001
011
111
110
101
010
100

(最後の2つは, 3ビットとれないので, 始めに循環している)つまり上から順に0, 1, 3, 7, 6, 5, 2, 4とすべて現れる.

Dijkstraのプログラムは, B(2,5)を作るものであった. まず元のプログラムを見ると

array d[0:36];
d[0] := 0; d[1] := 0; d[2] := 0; d[3] := 0; d[4] := 0; 
d[5] := 0;
k := 5;
4 < k and k < 36 herhaal
begin j := k-1; geluk := true;
  j >= 4 and geluk herhaal
  begin gelijk := true; h := 0;
    h <= 4 and gelijk herhaal
    begin gelijk := (d[k-h] = d[j-h]); h := h+1 end;
    gelijk zoja geluk := false; j := j-1
  end;
  geluk dan
  begin k := k+1; 
    k <= 31 dan d[k] := 0 anders d[k] := d[k-32] end
    anders
  begin k > 31 zoja k := 31;
    d[k] = 1 zoja k := k-1;
    d[k] := 1
  end
end

見慣れないキーワードがあるが, オランダ語らしい. それでも何となく意味が分かるので, Algol 60風に書き直してみた. ただし型宣言は省略した.

array d[0:36];
d[0] := 0; d[1] := 0; d[2] := 0; d[3] := 0; d[4] := 0; 
d[5] := 0;
k := 5;
while 4 < k and k < 36 repeat
begin j := k-1; geluk := true;
  while j >= 4 and geluk repeat
  begin gelijk := true; h := 0;
    while h <= 4 and gelijk repeat
    begin gelijk := (d[k-h] = d[j-h]); h := h+1 end;
      if gelijk then geluk := false; j := j-1
  end;
  if geluk then
  begin k := k+1; if k <= 31 then d[k] := 0 
    else d[k] := d[k-32] end
    else
  begin if k > 31 then k := 31;
    if d[k] = 1 then k := k-1;
    d[k] := 1
  end
end

最初の37語の配列 d は, 32ビットのde Bruijn-sequenceに循環して使う部分の5ビットを追加したものらしい. 最初の6個にとりあえず0を入れた. またkを5にした. kは配列 d で作業中の場所を指しているらしい.

その次 while 4 < k and k < 36 repeat は, kが増えたり減ったりしなが進行し, この範囲を脱出したら終りということであろう.

次は, k, k-1, k-2, k-3, k-4の5ビットと同じパターンがそれ以前のj, j-1, j-2, j-3, j-4ににあるかを見るループである. したがって, まずjをk-1にし, kの5つ組とjの5つ組を比較する.
途中の比較の結果を覚えておくのが, Boolean変数 glukで, まず真にしておく. hとglijkが5つ組の比較を制御する.

hを0から4まで変えながら, 5つ組に相違があればglijkを偽にし, ただちにループから脱出する. 5つが一致していれば, glijkは真のままだ. 一致のものが1組でもあれば, kからの5つ組は使えない. それを伝えるのがglukで, glukが真なら, 過去に同じものはなかった; 偽なら同じものがあったことになる.

それを使うのが, 中ほどの if glukである. 同じものがなければ, kの5つ組は合格で, kを増やし, d[k]に0を入れ, テストに戻る. 失敗なら, d[k]が1の間, kを減らし, 最初の0が見つけ, それを1にする.

とまぁこういうプログラムであった. kが31を越えてからは, はじめの方の値をコピーする. (d[k] := d[k-32])

大体そういうことが分かったので, 実行してみた. それにはSchemeに書き換える.

(define (dijkdebr)
 (define d (make-vector 37))
 (define k 5)
 (define (db0)
  (if (and (< 4 k) (< k 36))
   (let ((j (- k 1)) (geluk #t))
    (define (db1)
     (if (and (>= j 4) geluk)
      (let ((gelijk #t) (h 0))
       (define (db2)
        (if (and (<= h 4) gelijk)
         (begin (set! gelijk (= (vector-ref d (- k h))
                             (vector-ref d (- j h))))
          (set! h (+ h 1)) (db2))))
       (db2)
       (if gelijk (set! geluk #f))
       (set! j (- j 1)) (db1))))
    (db1)
    (if geluk
     (begin (set! k (+ k 1))
       (vector-set! d k (if (<= k 31) 0 
                            (vector-ref d (- k 32)))))
     (begin (if (> k 31) (set! k 31))
            (if (= (vector-ref d k) 1) (set! k (- k 1)))
            (vector-set! d k 1)))
    (db0))))
 (for-each (lambda (i) (vector-set! d i 0))
   '(0 1 2 3 4 5))
 (db0)
 (display d) 'ok)
(dijkdebr)

一番そとのループに入った時点で, 配列 d を出力したのが, 次である.

左は k の値. その右がdのその時確定している部分である. kを減じている時は, 短くなっている.

 5 000000
 5 000001
 6 0000010
 7 00000100
 8 000001000
 9 0000010000
10 00000100000
10 00000100001
 9 0000010001
10 00000100010
10 00000100011
11 000001000110
12 0000010001100
13 00000100011000
14 000001000110000
15 0000010001100000
15 0000010001100001
14 000001000110001
13 00000100011001
14 000001000110010
15 0000010001100100
15 0000010001100101
16 00000100011001010
17 000001000110010100
18 0000010001100101000
18 0000010001100101001
19 00000100011001010010
19 00000100011001010011
20 000001000110010100110
20 000001000110010100111
21 0000010001100101001110
22 00000100011001010011100
23 000001000110010100111000
24 0000010001100101001110000
25 00000100011001010011100000
25 00000100011001010011100001
24 0000010001100101001110001
25 00000100011001010011101011
26 000001000110010100111010110
27 0000010001100101001110101100
27 0000010001100101001110101101
28 00000100011001010011101011010
28 00000100011001010011101011011
29 000001000110010100111010110110
29 000001000110010100111010110111
30 0000010001100101001110101101110
30 0000010001100101001110101101111
31 00000100011001010011101011011110
32 000001000110010100111010110111100
33 0000010001100101001110101101111000
34 00000100011001010011101011011110000
35 000001000110010100111010110111100000
31 000001000110010100111010110111110000
32 000001000110010100111010110111110000
33 000001000110010100111010110111110000
34 000001000110010100111010110111110000
35 000001000110010100111010110111110000
36 0000010001100101001110101101111100000

Dijkstraには, Stepwise refinementでプログラムを作る論文がいくつもあり, それらは, 彼がどう考えてプログラムを考えたかが, 克明に記されている.

このプログラムのように, 結果だけ書いてあるのは, 非常に珍しいが, これでも, 何を考えたかは, 多少推測される.

P.S. プログラムに現れる2つのBoolean 変数のgelukとgelijkだが, 手元のオランダ語の辞書によると, gelukは幸運, gelijkは等価という意味らしい.