2009年4月22日水曜日

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は等価という意味らしい.

2009年4月13日月曜日

九元連立方程式求解機

情報処理学会歴史特別委員会では, 本年3月2日に情報処理技術の貴重な物件に対して「情報処理技術遺産」を認定した. その詳細は学会のページを見て欲しい.

その遺産の中に, 九元連立方程式求解機というものがある. 要するに連立方程式を解くアナログ計算機である. アナログ計算機と聞くと, なとなく眉に唾をつけたくなるのが, 現代の計算機屋の人情だが, なかなか面白い面もある.

その1つの課題が, 精度の高め方だ. われわれはアナコンと聞くや, 精度は悪いと思い勝ちである. まずこの機械の目盛リは, マイクロメーターで設定し, 滑車にはボールベアリングを使うという気の使いようではあるが, それでも程度は知れていると思う.

この計算機の使い方を読むと, 逐次的に機械を使うことで, 精度を高めていくと書いてあったので, その理由を考えてみた.

2元連立方程式
a00x + a01y = c0
a10x + a11y = c1 (0)
を解きたいとする. そして上のアナコンで, x',y'が得られたとする. その時
a00x' + a01y' = c0'
a10x' + a11y' = c1' (1)
と書くと, x'=x, y'=y (つまり正確な値が得られた) なら c0'=c0, c1'=c1 である.
違うなら, (0)-(1)
a00(x-x') + a01(y-y') = c0 - c0'
a10(x-x') + a11(y-y') = c1 - c1' (2)
を作る.

右辺の定数を c0 - c0', c1 - c1'として, もう一度この連立方程式を解くと, 前のx'、y'の代りに, x-x', y-y'が得られる. これをそれぞれ前に得られた値, x', y'に加えると, x, y が得られるということで, これを繰り返すらしい.

手元に計算機があるから, さっそくやってみることにする.

添字が面倒なので, ここでは
ax+by=c
dx+ey=f
と書くと
x=(ce-bf)/(ae-bd), y=(af-cd)/(ae-bd)
と解ける. これにわざと誤差を入れてみる.

Schemeのプログラムはこう書いた.

(define (solv a b c d e f)
(let* ((det (- (* a e) (* b d))) ;分母の行列式
(x (* (/ (- (* c e) (* b f)) det) 1.1)) ;誤差を入れる
(y (* (/ (- (* a f) (* c d)) det) 0.9)) ;誤差を入れる
(c0 (+ (* a x) (* b y)))
(f0 (+ (* d x) (* e y))))
(list (list x y) a b (- c c0) d e (- f f0))))

c0とf0は上のc0'とc1'のこと.
結果は最初に新しいxとyのリスト. 続いて次に解くための引数をリストにして出す.

x=2, y=3と決め. a00=2, a01=1, a10=1, a11=2とすると,

c0=a00x+a01y=2*2+1*3=7
c1=a10x+a11y=1*2+2*3=8

だから

2x+y=7
x+2y=8

を解くことになる
(solv 2 1 7 1 2 8)=>
((2.2 2.7) 2 1 -.1 1 2 .4)

最初の近似解は x=2.2, y=2.7である. また
c0'-c0=-.1,
c1'-c1=.4.

以後次のように進行する.
(solv 2 1 -.1 1 2 .4)=>
((-.22 .27) 2 1 .07 1 2 .08)

(+ 2.2 -.22)=>1.98
(+ 2.7 .27)=>2.97

(solv 2 1 .07 1 2 .08)=>
((2.2e-2 2.7e-2) 2 1 -1.e-3 1 2 4.e-3)

(+ 1.98 2.2e-2)=>2.002
(+ 2.97 2.7e-2)=>2.997

(solv 2 1 -1.e-3 1 2 4.e-3)=>
((-2.2e-3 2.7e-3) 2 1 7.e-4 1 2 8.e-4)

(+ 2.002 -2.2e-3)=>1.9998
(+ 2.997 2.7e-3)=>2.9997

(solv 2 1 7.e-4 1 2 8.e-4)=>
((2.2e-4 2.7e-4) 2 1 -1.e-5 1 2 4.e-5)

(+ 1.9998 2.2e-4)=>2.00002
(+ 2.9997 2.7e-4)=>2.99997

従って
xは2.2 1.98 2.002 1.9998 2.00002
yは2.7 2.97 2.997 2.9997 2.99997
のように2と3に近づく.

我々は計算機パワーを大活用して, 一気に計算しようとするが, アナコンのような低精度の計算機を使っても, 工夫次第では有効な計算が出来ると改めて思う.

2009年4月9日木曜日

moebiusの帯

今回もコンピュータ作図の話題である.

以前の小田急の社章は「小」の文字の縦の線を線路の断面の「エ」のような形にし, 両側の点を円弧にして囲み, 円に縦棒が通っているような形と記憶している. 最近の車体にはMoebiusの帯のような図案の後に, Odakyuと書いてある. それを見ているうちに, またMoebiusの帯をプログラムで描いてみたくなった.

moebiusの帯では, Escherの版画にアリが歩いている有名な絵がある. 例えばこれ.
アリが動いている動画もあった.

それほど芸術的には出来ないが, とりあえず

のような書き方を考えた. まず上の図は, 平行四辺形ABCDが, くっきりと折り目をつけたmoebiusの帯である. 外周が六角形になる. 帯の長さ方向の辺ABの長さをaとし, 短い方の辺BCの長さをbとする. この平行四辺形を120度ずつ回転して3回書くわけだ. 回転の中心をOとし, 辺ABの中心HからOまでの長さをcとする. またOHとCDの交点をG. DGの長さをd, OGの長さをfとする.

d=(a-b)/2 c=(a/2+b)tan 30 f=d tan 30

である. 従ってOを原点とすれば, dだけ左, fだけ下のDへ行き, b/2だけ左, c-fだけ下へ移動してA, aだけ右へ移動してB, b/2だけ右, c-fだけ上へ移動してC, そして線を閉じる.

d neg f neg moveto %D
b -2 div f c sub rlineto %A
a 0 rlineto %B
b 2 div c f sub rlineto %D
closepath stroke

で描ける. このサブルーチンを例えばpと定義すると, p 120 rotate p 120 rotate p 120 rotate でくっきり線のmoebiusの帯が出来る.

でもやはりゆるやかに曲がった帯が描きたい. しかしそう大変でもないことが分かる. 上の図で, A,B,Cの角を円弧にすればよい. PostScriptの関数arcは円弧の中心のx,y座標, 円の半径r, 反時計まわりに中心から見た出発の角度と到着の角度をパラメータにし, 前に描き終った点から, この出発点までを直線で結び, その後に到着点までの円弧を描いてくれる. 時計回りにしたければarcnを使う.

上のプログラムでは rlineto で相対的な位置までの直線を描いたが, 今度は, translateを使い, 座標をA,B,Cに移動しながら描くことにする.

d neg f neg translate 0 0 moveto
b -2 div f c sub translate e r r 150 270 arc
a 0 translate e neg r r 270 330 arc
b 2 div c f sub translate e neg r r 330 270 arcn
closepath stroke a neg 0 translate d f translate} def

最初の行は, Dへ原点を移し, 0,0へペンを置いた.
次の行は, Aへ原点を移し, e,rを中心, 半径rで150度から270度まで円弧を描いた.
次の行は, Bへ原点を移し, -e,rを中心, 半径rで270度から330度まで円弧を描いた.
次の行は, Cへ原点を移し, -e,rを中心, 半径rで330度から270度まで時計まわりに円弧を描いた.
次の行のclosepathで0,0へ線を引き, 続く2つのtranslateで原点をOへ戻している.

最後にtranslateで戻す代りに, 描画環境をスタックに保存してくれる, gsaveとgrestoreを使う方が簡単であろう.