今回はそのサブルーチンの中身の説明から始めよう.
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ビットである. 番地部が2nなら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' 3を加え, 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より小さいときは別になにもしないようだ. 不思議な仕様である. 他の例を見ても正の場合しか考慮していない.
続いて213=8192を掛けるのだが, 全体を小数にする必要があるので, 十進の小数の4桁目に小数点が來るように, 213/10000を掛ける.
乗算VnDは, 乗算レジスタにある数とnDにある数を掛けてアキュムレータに足すから, TDでアミュムレータをクリアしておく. そして2MDにある数を掛ける. しかし2MDのようにコードレターを2つ書くことは, EDSACのイニシアルオーダーでは出来なかったので, Dの代わりにπを書いていた.
2Mと3Mにあるのが, 213/10000である. なぜそうかというと
0.819210= 0.110100011011011100010111010110001110001000011001011012
これを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でプログラマー人生を始めた頃はまったくこの通りであった.
0 件のコメント:
コメントを投稿