EDSACのプログラム技法

前回のブログではy' を計算するサブルーチンにジャンプするところまで述べた.

今回はそのサブルーチンの中身の説明から始めよう.



0行目と1行目はlinkageで主ルーチンにもどる命令を作り, 22行目に入れる.

2行目は8Dにあるt' をアキュムレータに取り出す. 3行目のEはアキュムレータが正ならジャンプする条件ジャンプ命令なので, t' が正なら6行目(6番地)へ行く. 負の場合は8Dを2回引いて正にし, 6行目で合流した時には絶対値になっている.

平方根をとるサブルーチンS2は, 被開平数を4Dで渡すのでT4D. そしてS2へジャンプ. 結果は4Dにある.

9行目からはt' の3乗を計算する. まずH8Dで8Dにあるt' を乗数レジスタに置く. V8Dで2乗をつくり, TDで0D番地へ置く. VDで3乗が出来る.

13行目のRDはアキュムレータの1ビット右シフトである. EDSACのシフト命令のシフト数は, 命令の一番右のビットの位置で決まる. だからRDなら1ビット右シフト. R1Fなら2ビット, R2Fなら3ビット, R4Fなら4ビットである. 番地部が2ⁿならn+2ビットシフトする. 命令語の上の方のビットは影響しないから, 1ビットシフトは番地部が奇数ならなんでもよい.

14行目からは5倍の作業だ. まずUDで0D番地に入れる. UはTと同様な格納命令だが, アキュムレータをクリアしない. L1Fで左へ2ビットシフトで4倍する. それに0Dにあった元の数を足して5倍を実現する. それを一旦0Dへ退避し, 4Dにあった平方根をとりだし, R512Fで(512は2^9なので)右へ11ビットシフトし, 2^-11倍する. それに5.2^-1t' ³を加え, y' を8Dに置いて主ルーチンに戻る. T8Dをやったので, アキュムレータは0になり, 0は正なので, 正ジャンプのE命令で戻ってきた.

以後の作業には前回のブログのプログラムを見てほしい.

18行目. y' を乗算レジスタに入れる. 19行目. アキュムレータに取り出す.

次はyが400を超えるときはtoo largeと出力する仕様なので, y' を400.2^-13と比較する.

アドレスのところが2^-15だから, P400*2^2Fと比較すればよく, 5Mのような定数になる. ただ5MにはP1600Dとあるから, 400.125*2^-13と比較している. ちょうど400の時はtoo largeにしたくないというつもりだろうか.

21行目のGは負ジャンプ. 400を超えた場合は22番地で999.2^-13を乗算レジスタに入れる. これも定数はアドレス部に999の4倍が書いてある. (too largeと出力する代わりに999で代用する.)

ところで, これは正の方向に400を超ているを見ているだけで, -400より小さいときは別になにもしないようだ. 不思議な仕様である. 他の例を見ても正の場合しか考慮していない.

続いて2¹³=8192を掛けるのだが, 全体を小数にする必要があるので, 十進の小数の4桁目に小数点が來るように, 2¹³/10000を掛ける.

乗算VnDは, 乗算レジスタにある数とnDにある数を掛けてアキュムレータに足すから, TDでアミュムレータをクリアしておく. そして2MDにある数を掛ける. しかし2MDのようにコードレターを2つ書くことは, EDSACのイニシアルオーダーでは出来なかったので, Dの代わりにπを書いていた.

2Mと3Mにあるのが, 2¹³/10000である. なぜそうかというと

0.8192₁₀= 0.11010001101101110001011101011000111000100001100101101₂

これをEDSACの短語に対応させてみると, EDSACでは偶数語が右, それ+1の奇数語が左なので



のようになる. この図で見ると2Mの語の最後はDだが, プログラムではFなのはこの辺で切り捨てたのであろうか.

掛け算が済んだので, 0Dに入れ, 小数の印刷ルーチンP14を呼ぶ.

このP14は28行目にあるプログラムパラメタで, レイアウトできる.

例えば小数0.123456789を
のように出力するには, スペースの位置の下の数と, 最後の桁の上の数を足したものがパラメタの値であり, 28行目は3072+32=3104となっている.

この後はt を取り出すアドレスと, カウンタi の変更が残っている.

29行目は14行目の命令を40D, 38D, ...と減らすもので, A14θ, S 9M, T 14θ で変更する. 32行目からは, A 7M, S 8M, U 7Mとカウンタを減らし, 35行目で, アキュムレータに正の数が残っていれば, 3番地へ戻る. 負になったらジャンプせず, 36行目のZ命令で停止する.

論文にあった出力はこうだ.

10  0000 04472
 9  0000 03162
 8  0999 00000
 7 -0044 85586
 6  0047 75414
 5  0999 00000
 4  0062 35157
 3  0999 00000
 2 -1077 55051
 1  0999 00000
    0018 09974

最下行 0 がプリントされていない. 論文には印刷ルーチンの虫だろうと書いてある. そうだろうか.

私がSchemeで計算したのは次だ.

(define (tpk t)
(+ (sqrt (abs t)) (* 5 (expt t 3))))

(define data 
'(1.5 8 -6 9.5 2.3 9.9 2.1 -2.1 6 0.001 -0.002))

(map tpk data)
=>
(18.09974487139159 2562.828427124746 -1077.5505102572167
 4289.957207001485 62.35157508881029 4854.641426544511
 47.75413767461895 -44.85586232538106 1082.4494897427833
 .03162278160168379 .044721319549995794)

という次第でこのプログラムは解明出来た.

いまから思うと大変な手間であったが, 私がパラメトロン計算機PC-1でプログラマー人生を始めた頃はまったくこの通りであった.