EDSACのプログラム技法

EDSACのプログラムを集積したウェブページがある. (詳しくはそこの "To browse the collection of software, click this link."をクリックする.) ここには私の書いたエラトステネスの篩も置いてあって嬉しい(21行目). そこにDavid J. Wheelerの素数のプログラムがあった(15行目). Martin Campbell-Kellyの書いたEDSACの文献(Programming the EDSAC: Early Programming Activity at the University of Cambridge, IEEE Annals of the History of Computing, Vol 1, No. 2.)によると, EDSACには「EDSACの本」で有名なイニシアルオーダーの前にもう一つのイニシアルオーダーがあって, 1949年5月6日の運転開始の日にはその方が使われていた. 同年6月22日にケンブリッジで開催された計算機会議でEDSACが公開され, その時, Wilkesの書いた平方の表とWheelerの書いた素数のプログラムのデモが行われた.

Wheelerは1927年生まれだから, 22歳のWheelerが当時どういうプログラムを書いていたかには非常に興味がある. これはもう読むしかない.

ところがその内容たるや

[Primes]
T107SO92SO93SO94SS5ST6SO95SO95ST7S
A96SR4SS97SE42SL4SA97SG45SS98SG100ST7S
A97SA4ST97SH97SN97SL64SL64SA96SE39ST7S
A96SU1SA4ST96SA99ST97SS88ST7SH91SA1S
E72SV91SS89SE71SA89STLOSH90SV1SL4S
TLA7SA98SG67SA6SA4SG36SE33SP2SP500S
JSP16S#S@S&S!SP2LP1LPLP1L
T7SA99ST97SA4SA96ST96SE39S

だけなのだ. これを読む手がかりはここにある.

読んでみると, やはり面白いプログラムであった. 特にそのプリントルーチンには感心した. 以下にその解説をしたい.





31番地 T107S. このプログラムは最初のイニシアルオーダーInitial Order1で読み込むが, 読み込まれるプログラムが106番地まであるなら, 最初をT107Sとするのであった. つまりイニシアルオーダーには格納命令があり, 最初はT31Sになっていて, 31番地から命令を格納する. つまり, T31S, T32S, ..と変わっていく. そして格納作業の終りは, 格納命令から31番地にある命令を引き, それが正であるかを見る. T107Sを引くと0(正)になるので31番地へ進んでくる. このT命令はアキュムレータのクリア用にも見えるが, すでに0になっている.

32番地は92番地にある#S(figure shift)を出力する. 以下のプログラムが数字を出力するので, その準備.

33,34番地では93, 94番地にある@S(cr), $S(lf)を出力. 次のS5Sはイニシアルオーダーの5番地にあるP5S(10*2^-16)を引き, -10を作り, 6番地(m[6])に入れる. 6番地は -10,-8,-6,..と変わり, 1行に5個の素数を印刷する.

37, 38番地. !S(space)を2個出力. (このプログラムのあちこちにある)T7Sは(67番地のものを除いて)アキュムレータをクリアする.

このプログラムでは素数の候補pから3,5,...と奇数(d)を引き(以後 dを減数という), 0になるかどうかを見ている. (前にも書いたようにEDSACには除算命令はない.) pは96番地, dは97番地にあって, それぞれ5と3に初期化されている.

40から46番地はそのpからdを引く. ただ初めは16*dを負になるまで引き, 次にdを正になるまで足す. これでpをdで割った剰余が得られる. 実際には16*dを引くのではなく, pを4ビット右シフトして引き, 4ビット左シフトして足す. (おおざっぱにいうと単にdを引くより, 16倍加速されている. こういう剰余の取り方はEratosthenesの篩でも使える.)

47番地. EDSACには零ジャンプがないので, 正になったところで1を引き, 負なら零であったということで, 合成数の判定をして100番地へ進む.

剰余が0でなければ, 50,51,52番地でdを2増やし, 53番地からH97Sで乗算レジスタにdを置き, N97S, L64S, L64Sで-d^2を作る. それにpを足し, 減数の二乗が素数候補より小さければ39番地へ戻り,　次のテストに掛かる.

減数が大きくなって素数と判定出来ると, 60番地から先で素数の出力にとり掛かる.

U1Sで素数を1番地に入れ, A4Sでpを2増やして63番地で96番地に入れておく.

素数は4桁で出力するつもりで, 66番地で-4を7番地のカウンタへ入れる.

68番地から76番地は1桁の出力で, 1番地に1234があるとして説明すると, 68番地H91Sで乗算レジスタに2^-11を入れる(32*2^-16=2^-11).

アキュムレータに1234を足す. 72番地へ飛んで89番地にある1000を引く. アキュムレータは234(*2^-16)になり, 正なので71番地へ行く. V91Sだからアキュムレータは234*2^-16+2^-22になる.

再び1000を引くと今度は負になり, 74番地で1000を足し, アキュムレータは234*2^-16+2^-22に戻る.

これを0番地の長語に入れる. つまり左半分の1番地は234, 右半分の0番地はファンクション部に1となる. そこでo命令は0番地を出力するから1が印字される(EDSACの文字コードは2012年9月28日のこのブログ参照).

77番地から80番地は234を10倍する. H90Sで乗算レジスタに10*2^-4を置き, V1Sでアキュムレータは234*10*2^-4になり, 次の左シフトで2340になり, 1番地に戻される.

81番地からは7番地のカウンタを1足し, まだ4桁出力していなければ, 67番地へ戻る.

次は2340だから, 71番地のループは2回実行され, 右半分の方は2*2^-11になって2が印字され, 1番地は3400になる.

4桁出力すると84番地へ来て, 6番地のカウンタにA4Sで2を足し, 5個未満なら36番地, 5個終わっていたら33番地へ戻る.

一方合成数の時は100番地へ進み, 101, 102番地でdを3に戻し, 103,104,105番地でpを2増やして39へ戻る.

この解釈が正しければ, 期待される出力は

  0005  0007  0011  0013  0019
  0019  0023  0029  0031  0037
...

なる筈であった. サイボウズ・ラボの西尾君がEDSACのシミュレータを持っているというので, 試してもらったらその通りであったとメイルが来た.

このプログラムの疑問は素数候補が5から始まることである.

減数を3から始めると素数候補の3は割り切れるから, 3は素数にならない. 3を素数にするには2から割り始めなければならないからなのだ.

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でプログラマー人生を始めた頃はまったくこの通りであった.

EDSACのプログラム技法

英国の計算機歴史学者Martin Campbell-KellyがProgramming the EDSAC: Early Programming Activity at the University of Cambridgeという長い長い論文をIEEE Annals of the History of Computing Vol.2 No.1 January 1980に書いている.

この最後の方にInfluence on programming for other machinesという章があり, 東大のTACのことも述べているが, E. Wada, who was responsible for programming, designed a new set of initial orders based on those of Wheeler, but incorporating floating addresses and some other small refinements.と書かれているのを知った.

今回のブログの目的はその論文にあるEDSACのプログラムの例を説明することである. MartinはEDSACの例にTPKプログラムを使った. なんの変哲もない問題だが, アセンブリ械語のプログラムとしては長さ的に適当かも知れない.

11個のデータを読み込み, 配列aの0から10に置き, t =a[t ]として,

y =√|t|+5×t ³ (*)

を計算, 印刷するものである.

t が<1なら, 最初の項は平方根だからtより小さくなるから無視し, 次の項は3乗なのでこれが結果の大きさに利く. t が±10位ならyは±5000程度のt³の曲線を考えればよい.

EDSACの頃の計算機は, 扱う値を±1以下としていたから, tが10程度までなら1/16にし, yも1/2¹³=1/8192として計算する.

つまり *の式にy=2¹³y', t=2⁴t'を入れ, y' =2^-11√t' +5×2^-1t'³ を計算し, 後で2¹³倍する. t=1とすると, t'=2^-4. その平方根は2^-2. 2^-11を掛けると2^-13. またtの3乗は2^-12. 5.2^-1を掛けて5.2^-13. 両方を足して2^13を掛けると6になる.

EDSACに用意してあるライブラリのサブルーチンもみなこの範囲と思って働く.
さてMartinの書いた主プログラムはこうだ.



上から数行目の 0 |A θ と書いてある行を0番地としたアドレスにはθを最後に付ける. 下の方 M 0|P 4 D の行を0番地とするときは, Mを付ける. そのMは38θであること, set M-parameterの行で設定する.

0と1行目はWheelerリンケージで, 56絶対番地のR1へサブルーチンジャンプする. 2行目 T20Dはサブルーチンに20Dからデータを格納せよとの指定である.

データはテープに15+8+6-95+23+99+21+21-6+0001+0002-Fとパンチしてあり, Fが來るまで0.15, 0.8, -0.6, ...のように20番地, 22番地, 24番地, ...,40番地に読み込まれる. そして3番地に戻ってくる.

O 10M, O 11Mは10MにCarriage return(復帰), 11MにLine feed(改行)のコードがあって, 出力装置を改行する. こういう機能鍵にはアルファベットが書いてないので, パンチの都合上ギリシア文字が割り当てられている.

4,5番地にT Dが2度あるが, 最初のはアキュムレータを0にするもの, 次はその0を0番地の長語に格納するものだ. そして7M番地にある10を0番地の短語に設定し, 9,10番地でP7へサブルーチンジャンプする. EDSACの短語は, EDSACのプログラム技法のブログの2012年9月22日の命令語の図にあるように, L/Sビットのところが整数の1の桁であり, アドレスはその上に1ビットずれているから, P5Fと書くと10になり, 8MのようにPD, つまりアドレス部が0でL/Sビットがあると, 整数の1になる. (ファンクション部のPは0である.)

11番地へ戻るとO 12Mでスペースを2文字出力する. 次はtを16で割りt'にするのだが, もともと読み込むデータが1/10になっているとしたので, (15は0.15と読まれたが, 1.5のつもりであった.) 10/16を掛けることになる.

EDSACの乗算は一方の数を乗算レジスタへ置き, ある番地にある数と掛けて積をアキュムレータに置くので, H 4Mで10/16を乗算レジスタに置く. 4Mを見ると J F である. F はL/Sビットが0のことだが, Jはなにかというと, 前のブログでEDSACの文字コードを見るとJは01010である. 左端の0の次に小数点があるから, これが10/16なのである.

14番地のV40Dで40番地にある長語のデータと乗算が出来た. それを8Dに置いてサブルーチンへ行き, y' を計算する. y' も8Dに置いてある.

まだまだ続くが, 今回はここまでにしよう.

EDSACのプログラム技法

私が2000年ころ書いたEDSAC用Eratosthenesの篩の説明 少し長いが興味と元気があればフォローして欲しい.

まずプログラムの先頭はこうなっている.



これはすべて定数か作業場所である.

以下のが本体.



本体のプログラムは先頭の T 96 K G K で以下のプログラムを96番地から格納する. 最後のE 96 K P F で96番地から実行を開始する. 以下相対番地を前回のように ' で示す.

0'の T F でアキュムレータをクリア. S 850 Dで850,851の P D, P F つまり長語の1を引き35ビットオール1を作る. T 992 Dで 992,993番地へ入れる. 2'番地をアキュムレータに置き, A 849 Fでアドレス部に2を足して2'へ戻す. 848の T 1024 D を引き, 負なら0'へ. このループで長語の992から1022まで16長語をオール1にして篩を用意する.

8'からのループは A 1022 Dでオール1にし, 850 D の1を引き, 111...10(負のストロブ)を作って長語852に入れる. アキュムレータはクリアされていてそれに 850 D の1を足し, 000...01 (正のストロブ)を作って長語922に入れる. これを左1ビットシフトし, 000...10 にして 850 D へ戻す.

16'から22'までは11'と13'の格納命令を2ずつ増やし, 結局長語852-920には負のストロブ, 922-990には正のストロブを用意する.

23', 24', 25'でFIGS, CR, LFを出力. 26', 27', 28'で2を作り, 29'で0番地へ格納して30', 31'のWheelerリンケージで出力ルーチン P6 へサブルーチンジャンプして最初の素数2を印字する. 33'からがいよいよ篩だ.

以下の表で左端の852の列は長語の番地 次の0の列はストロブや篩の長語の番号 その右の二進表示が長語の内容である.

 852  0 11...110 負ストロブ
 854  1 11...101
 856  2 11...011
...
 920 34 01...111

 922  0 00...001 正ストロブ
 924  1 00...010
 926  2 00...100
...
 990 34 10...000

            9753←対応する奇数の素数
 992  0 11...111 篩
 994  1 11...111
 996  2 11...111
...
1022 15 11...111 

0番の篩の長語の各ビットはその上に書いてあるように右から3,5,7,9,...に対応する.従って最後の15番の篩の左端のビットは1121に対応する. EDSACのメモリー容量からすれば篩はもっと大きく出来るが, 篩の様子を水銀タンクで眺められるように, ちょうど1つのタンクの大きさにした.

篩は正のストロブを順に使い, 篩の右から順のビットとandをとって, 素数か合成数かをみる. 合成数なら次のストロブで次の篩のビットを調べに進む. 素数ならその素数をp とすると, 左の図のようにpビットごとにpの倍数があるから, 負のストロブを使っていま調べたビットからpビットおきに篩のビットを0に変えていく. 手順としてはこれだけだが, ストロブが下まで来たときに, 篩を次の長語にする; ストロブを先頭に戻す作業が必要である.

33'で正のストロブを乗数レジスタに置き34'で篩とandをとる. 844Dの長語の1を引き負なら合成数だったので75'へ. そうでないなら37'から素数出力; 841に素数の4倍があるので42'で右2ビットシフトし, 44',45'でP6へいく.

46'からは素数の倍数の篩を0にする, つまり篩う仕事を開始.

以下抜き書きしたプログラムは左から命令のある相対番地, 命令, 番地の示す場所の内容, つまりオペランド, 演算結果のアキュムレータ, 格納命令の場合は格納番地と内容を示す.

      命令     オペランド  アキュムレータ
 47   A  34 @   C 992 D      C 992 D 
 48   U  71 @                              71 C 992 D
 49   S 843 F   L     F      T 992 D
 50   T  72 @                              72 T 992 D
 51   A  33 @   H 922 D      H 922 D
 52   S 842 F   P  70 F      H 852 D
 53   A 841 F   P   6 F      H 858 D       
 54   U  70 @                              70 H 998 D
 55   S 840 F   H 992 D      P   6 F 
 56   G  69 @

47'で今素数を調べた篩を取り出し71'へ入れ, ファンクション部をTにした命令を72'へ入れる. 51'は素数を調べた正ストロブを乗数レジスタに入れる命令で, 篩を0にするには対応する負ストロブが必要である. それは7番地若い方にあるからまず70を引き, 素数のビット数だけ先から篩うので841番地にある素数の4倍の値を足し, 消すストロブの 番地を作る. 55'は負ストロブの範囲を超えたかのチェック.

範囲を超えていなければ69'へ来る. 以下のプログラムの70',71',72'は先ほどのプログラムが設定したものだ.

 69   T     F
 70   H 998 D
 71   C 992 D
 72   T 992 D
 73   A  70 @   H 998 D      H 998 D
 74   G  53 @
 53   A 841 F   P   6 F      H1004 D
 54   U  70 @                        70 H1004 D

73'からは篩をクリアする命令を更新する. H 998 Dをもって53'へ. そこで再び素数の数だけストロボのアドレスを増やす. こうしている内に負ストロブの範囲を飛び出し57'へ進む. 58'で篩の番地を取り出し, アドレスを2増やし71'へ戻す. 72'のT命令も増やす. しかし篩の方が範囲を超える心配があるので, 63'でT1024 Dを引き, 篩の範囲が終わっていなければストロブの番地をもとへ戻す. それが66'からでストロブを取り出し, 番地から70を引いて1周戻し, 54'へいって次のサイクルに入る

篩の最後まで来たら, 合成数の時と同じく75'へ行く.

75'からは次の奇数の素数テストになる. 76'からは素数の4倍をもっていたカウンタ841を4増やす. 79'からは正ストロプを乗数レジスタに置く33'の命令を2増やし, 82'で正ストロブの範囲をチェックし, 範囲内なら32'へ行って左隣のビットをテストする. 正ストロブが範囲を超えているなら, 84'で正ストロブを先頭に戻し, 86'から篩の長語を次にすすめ, 89'で篩も範囲を超えたなら91'でプログラムを終止する.

初期の機械語命令の時代はこういうプログラムは普通であった.

このプログラムが出来た時, EDSACのシミュレータを公開している英国Warick大学のMartin Campbell-Kelly君に送ったらシミュレータのホームページにWadaSieveとして公開してくれた.

彼は"The program is beautifully written and very fast. A little master work."とコメントしている.

EDSACのプログラム技法

EDSACのプログラムの読み書きに重要なのが文字コードである. コードの知識なしでは仕事にならない. それ以前の計算機での入力法はよくわからないが, EDSACではプログラムを紙テープにパンチし, それをイニシアルオーダーで読み込むことが例の本で公開されたので, イニシアルオーダーを解読するためにも文字コードに関心をもつことになった.

EDSACは大学で作った計算機なので, 入出力は市販の機器を用いることになる. 当時はもちろんテレタイプは存在していたから, 当然それらを使うわけだ.

その頃 欧米で使われていたテレタイプは, 英大文字と数字と若干の記号のもので, 5単位テープが標準であった. このコードはInternational Telegraphic Alphabet No.2 (ITA2)という. (規格はここのFreely available itemsのEnglishから得られる) この種のテレタイプは昔は国際電電にいくと見られたがとおの昔に姿を消した.

テレタイプとしてはこのような形であった.



なかなか可愛らしくて1台ほしい. 数年前にアメリカのある空港で耳の不自由な人のための似たようなキーボードを見たことがあった.

コード表は次のとおり.



Figure shift(FIGS)を打つとその後は右の数字と記号を送受信することになる. Letter shift(LTRS)で英大文字になる.

このコードの特徴は, 使用頻度の高い文字は1が, つまりテープの孔が少ないことだ. Wikipediaでletter frequencyを見て, 頻度の順に孔の数をならべると

e t a i n o s r l d h c u m f p y g w v b k x j q z
1 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 3 4 4 3 4 2

孔が少ないとテープが丈夫なのか, ごみを減らそうとしたのか知らないが, Morseコードと同様な精神で出来ている. ITA2の特殊記号も欧文Morseコードと同じである.

さてEDSACではこういう機器を使おうとしたが, 数字のコードが

P 0 01101
Q 1 11101
W 2 11001
E 3 10000
R 4 01010
T 5 00001
Y 6 10101
U 7 11100
I 8 01100
O 9 00011

と数値に関係がない. そこで最上段の文字のコードを00000から01001に変更した. ただそうするとP 0のコードはまったく孔が開かず, ブランクテープと同じになってしまうので, テープでは一番左の位置はコードが0の時に孔が開くようにした.

変更は最小にしたらしいが, T→F→Y→T, O→D→W→O, R→X→U→S→R, I→A→I, P→B→E→P, Z→Q→Zと交換した. ITA2と同じなのはH, N, M, L, G, C, V, J, Kである.

これがEDSACの文字コードだ.



左の方がPerforator, テープ鑽孔機で, そのすぐ右がTeleprinter, 出力用タイプライターである. 図形文字でない機能鍵のFigure shiftもプログラムでは文字として使いたいというので, 鑽孔機にはπの文字がついているという具合である.

EDSACの最初の本のコードでは, 鑽孔機の記号の位置はタイプライターのそれとずいぶん違っていた. その辺は上の表では省略した.

EDSACのプログラムでは命令を A 3 F のように書くが, 数字の前にFigure shiftを打つかというとそうではない. 命令の最初は文字と思い, そのあとPQWERTYUIOJが続くあいだは数字として扱う. Jは10として使え, そう使うプログラムを存在する.

一方数値だけの疑命令でも, 先頭にはPを書かなければならない. 命令はコード表のFからVまでの文字で終わる. π, S, Z, Kには別の機能があった.

EDSACはこのようにコードを変更したが, 我々のパラメトロン計算機では, コードは市販の機器のままで, あとはプログラムで挑戦するという方法をとった. 日本国内では6単位が標準だったので文字数も多く, プログラムが見やすいシステムが作れた. その話はまたいつか.

EDSACのプログラム技法

英国ケンブリッジ大学のEDSAC(Electronic Delay Storage Automatic Calculator)が稼働し始めたのは1949年5月6日であった.

EDSACのメンバーはプログラムの解説書(Maurice V. Wilkes, David J. Wheeler, Stanley Gill: The Preparation of Programs for an Electronic Digital Computer)を早々に出版した. 世界で最初のプログラムの本であり, 多くの人がそれでプログラミングの楽しさを知った. 私もその1人である.

遠い昔のテクニックはどんどん忘れ去られる. 最近EDSACのプログラムを眺める機会があったので, その頃のプログラムを解説したくなった.

当時のメモリーは水銀遅延線かブラウン管であった. EDSACはそのDelay Storageから分かるように, 水銀遅延線の記憶装置を持つ. つまり水銀タンクの一方からPiezo効果を使い音波を送り, 反対側で受けた音波を逆Piezo効果で電 気信号にもどす. その遅れを記憶装置として利用する. (水銀遅延線の写真はここに)

EDSACのアーキテクチャは簡単である. レジスタとしては71ビットのアキュムレータと35ビットの乗算レジスタ. 記憶装置は17ビットの短語が1024語. 2n番地と2n+1番地の短語を2語つなぐと35ビットの長語になる. 不思議な計算だがその秘密は, 水銀遅延線での短語は18ビット分あり, 語と語の切替えに1ビットの時間がかかるので, 長語は36-1ビット, 短語は18-1ビット, アキュムレータは72-1ビットということだ.



命令語は左端の5ビットが英文字1字で表すファンクション部, 次の11ビットがアドレス部, 最後の1ビットがL/S(long or short)部で, 0だとアドレス部の短語, 1だと長語を示す.

数値は2の補数表示で, -1≤x<1の範囲の値を持つ. 左端のビットが1なら負数である. 具体的には左端の方だけを示すと, 0.01₂は0.25, 0.1₂は0.5, 0.11₂は0.75, 1.00₂は-1, 1.1₂は-0.5, 1.01₂は-0.75である.

EDSACの命令は A n F のように書く. 先頭の英字がファンクションで, つづいて十進でアドレスを書き(0なら何も書かない), 最後のコードレターという英字を書く. コードレターはアドレスの終りを示す他, FはL/Sビットが0, Dは1を示す.

A n はn番地の内容をアキュムレータに足す; S n は引く. T n はアキュムレータをn番地に格納してアキュムレータをクリア. U n は格納するがクリアせず. E n Fはアキュムレータが正ならn 番地へジャンプ, G n F は負ならジャンプで, 無条件ジャンプはなかった.

乗算は乗数を H 命令で乗算レジスタに置き, V n はn番地と乗算レジスタを掛けて積をアキュムレータに足す; N n は積を引く.

EDSACに除算命令はない.

基本的な説明は以上で終る. サブルーチンジャンプの説明も要ろう. n番地からm番地にあるサブルーチンを呼ぶには, 引数を指定された場所に入れ(T 命令で入れるからアキュムレータはクリアされている)n番地に A n F の命令を置く. n+1番地に E m F を置く. 英文字Aは左端のビットが1なので, アキュムレータは負, 従って E 命令でm番地へジャンプする.

m番地のサブルーチンに来たときは, アキュムレータに A n F があるから, これに U 2 F を足す. EDSACの文字コードでは A+U=E だから和は E n+2 F の命令になり, これをサブルーチンの最後に格納する. サブルーチンからはこの命令を実行してn+2番地へ戻る. これをWheelerリンケージという.

次は二進法の小数を十進に変換して出力する方法.

0.0001₂は十進では0.0625である. これを10倍する. 二進の方は0.1010₂, 十進の方は0.625 この時の整数部が元の小数の1桁目だが, 今は0だから0を出力.

また10倍する. 二進は110.0100₂, 十進は6.25. 整数部は6だから6を出力. 整数部を捨ててまた10倍する. 二進は10.1₂, 十進は2.5. そこで2を出力. さらに10倍して5を出力. という次第で0625が出来る.

EDSACのプリントルーチンP6は次のように書いてある.



左の0から31は相対行番号である. 欄外の24→9のようなのは, 24番地から9番地へジャンプして來るの意味. 25行目の(E F)のかっこは, この命令は変更されること. その下の横線は, 無条件ジャンプを示す. 29番地から31番地の左の2本の縦線は, 偽命令, つまり命令の形はしているが実行されず, 定数であることを示す.

最初の G K はKで終っているから, 制御指令で, 次の命令が格納される番地を相対番地の原点に設定する. つまり命令がθ(短語)やπ(長語)で終っている時は, A 3 F 命令の場所を0 とする. 以下では相対番地を'で示す.

0',1'番地はWheelerリンケージである. 帰り命令は25'番地に格納される. 2'番地は29'番地の偽命令を乗算レジスタに置く. Jは01010, 995は二進では1111100011₂で, 最後がFつまり0だから, 命令の形としては
01010 01111100011 0
であり, 数値的には
(/ #b01010011111000110 65536.0) => .655364990234375
だから, 2¹⁶/10⁵を上に丸めたものである.

0番地には短語の範囲に5桁の整数があり, それを5桁で出力するのだが, 最大の99999と最小の1とを見ると,

1 ]=> (* (/ 99999 65536) .655364990234375)

;Value: .9999976144172251

1 ]=> (* (/ 1 65536) .655364990234375)

;Value: 1.00000761449337e-5

だから, 純小数にして5桁出力する方針である.

4'番地の T 4 D でこうして出来た長語を4,5番地に置く. 5',6'番地は3'番地の命令V F を0番地に入れる. Vは11111なので, コメントのように-1/16を置くことになる. 0番地の負は, 出力の先頭の0をスペースにする目印である. 7'番地は乗算レジスタに10/16を置き, 10倍の準備をする. 8'番地は6'番地の命令を0から引き, Tが00101, 5なので1番地のカウンタを-5/16にする. 10'番地は4,5の長語に10/16を掛ける, すなわち10倍して小数点を4ビット右へ. つまり整数部が先頭の5ビットに収まる.

11'番地ではアキュムレータを残したまま積を戻し, 12'番地の A F で6'番地で入れた0番地の-1/16を足す, つまり整数部が0だったら-1/16になり, 26'番地へ進む. 26'番地で先程引いた1/16を足し, O 31 θ で空白を出力, 20'へ戻る.

整数部が0でなければ, 13'番地のジャンプは起きず, 14'番地の T F でアキュムレータをクリアし, 15'番地の T F で0番地を0にする.

16'番地の O 5 F は, 4番地の長語の先頭の5ビットが数字のコードなので, それを出力する.

17'番地で4番地の長語を取り出し, F 4 F では O 5 F で出力レジスタにいれた数字のコードを4番地の先頭の5ビットに取り出す. 残りのビットは0. それを19'番地で引き, 整数部をなくしてから, L 4 F で左へ4ビットシフト, 4番地へ戻す.

EDSACのシフトは難しい. 命令の下位から見て最初の1のビットが現れるまでシフトする. だから L D は1ビット, L 1 F は2ビット, L 2 F は3ビット, L 4 F は4ビットシフトになる. それより上に1のビットがあっても影響しない.

22'番地からは5桁出力したかのチェック. 最初に入れた1番地の-5/16から3'番地の-1/16を引く. 初めの4回は負なので9'番地へ戻る. 5回済むと25'から主ルーチンに戻る.

思わず長い説明になったが, EDSACの頃のプログラムはこういう感じであった.

このプログラムの問題点は0を出力すると, 5個の空白になることだ. 当時はメモリーが少く, プログラムを短かくするのが重要であった.

EDSACについては, シミュレータのTutorial Guideを参照されたし. その23ページにP6の記述がある.