イニシアルオーダーは十進二進変換もするが, その際の10倍に乗算を使わず, シフトしたり足したりしているのは, 当時, 乗算命令が未実装だったからである.
従って, 始めの頃は, 乗算や除算のサブルーチンが存在した. 私は数値計算には興味がないので, 除算は必要なく, 高橋先生が作られたはずの除算のサブルーチンがどうなっていたかは知らない. 先生のライブラリのファイルが手元に来たので, 除算サブルーチンを探したが, 見つからなかった.
ところで, PC-1がお手本にした, Cambridge大学のEDSACには, 最後まで除算はなく, サブルーチンを使っていた. そういえば, EDSACでは除算はどうやったか知りたくなり, 例のEDSACの本を見ると, 思いの他簡単であった.
a0←被除数,
c0←除数-1
と初期化し, 漸化式
ak+1←-akck+ak,
ck+1←-ck2.
で繰り返すと, ck=0の時, akが商になる.
これをそのままSchemeでプログラムにする.
(define (div n d) ;n 被除数 d 除数
(define (dv a c)
(display (list a c)) ;途中の様子を見る
(if (< (abs c) 0.0001) a
(dv (- a (* a c)) (- (* c c)))))
(dv n (- d 1)))
PC-1もEDSACも当時の計算機は, -1≤n<1の固定小数点数しか扱わない. 利用者は, 計算がこの範囲になるよう, 絶えず位取りに注意する必要があった.
そこで, 上のサブルーチンを(div 0.4 0.8), (div 0.4 0.5)でやってみると,
> (div 0.4 0.8)
(.4 -.19999999999999996)
(.48 -.03999999999999998)
(.4992 -1.5999999999999983e-3)
(.49999871999999995 -2.5599999999999945e-6)
=> .49999871999999995
> (div 0.4 0.5)
(.4 -.5)
(.6000000000000001 -.25)
(.7500000000000001 -.0625)
(.7968750000000001 -.00390625)
(.7999877929687501 -.0000152587890625)
=> .7999877929687501
正しい商は, 0.5と0.8である. 目がちらちらするから, 途中のak, ckを整数でシミュレートしてから同様な表にすると,
0.4/0.8
k a c delta
0 0.4 -0.2 0.1
1 0.48 -0.04 0.02
2 0.4992 -0.0016 0.0008
3 0.4999872 -0.00000256 0.0000128
4 0.49999999 -0.00000000 0.00000000
99967232 00065536 00032768
0.4/0.5
k a c delta
0 0.4 -0.5 0.4
1 0.60 -0.25 0.2
2 0.7500 -0.0625 0.05
3 0.796875 -0.00390625 0.003125
4 0.79998779 -0.00001525 0.00001220
296875 87890625 703125
左端のkはaとcの添字, 右端のdeltaは正確なn/dとakとの差である. 最後のk=4の行は長いので, 少数点以下8桁目で折り返してある.
こういうのの性質をしらべるには, 各回におけるinvariantを発見するのが常道である. ためつすがめつするまでもなく, この例ではcとdeltaには比例関係があることが分かる. つまり, 上の計算ではdelta/cは-1/2, 下のではdelta/cは-4/5, つまり除算の商の符号を変えたものであった.
すなわち, n/d=ak-(n/d)*ckということだ.
ここまで分かると, 次は帰納法で証明するだけである.
k=0の時
a0-(n/d)*c0=n-(n/d)*(d-1)=n/d.
k=kで
ak-(n/d)*ck=n/d
が成立しているとする. その時k+1の式は
ak+1-(a/d)*ck+1
=-akck+ak-(n/d)*(-ck2)
=-akck+ak+(n/d)*ck2
=ck(ak-(n/d)*ck)+ak
=ak-(n/d)*ck)=n/d.
やはり, n/dになる. 一方 ckはどんどん小さくなっていくから, ck=0の時, ak→n/d.
そういう風にEDSACの除算サブルーチンは出来ていたのであった.