除算サブルーチン

パラメトロン計算機 PC-1 が完成し稼働し始めたのは, 1958年3月26日ということになっているが, この時は乗算も除算も出来なかった. その後, それらの回路を設計し, やがて乗算が付き, 夏頃には除算も出来るようになった.

イニシアルオーダーは十進二進変換もするが, その際の10倍に乗算を使わず, シフトしたり足したりしているのは, 当時, 乗算命令が未実装だったからである.

従って, 始めの頃は, 乗算や除算のサブルーチンが存在した. 私は数値計算には興味がないので, 除算は必要なく, 高橋先生が作られたはずの除算のサブルーチンがどうなっていたかは知らない. 先生のライブラリのファイルが手元に来たので, 除算サブルーチンを探したが, 見つからなかった.

ところで, PC-1がお手本にした, Cambridge大学のEDSACには, 最後まで除算はなく, サブルーチンを使っていた. そういえば, EDSACでは除算はどうやったか知りたくなり, 例のEDSACの本を見ると, 思いの他簡単であった.

a₀←被除数,
c₀←除数-1
と初期化し, 漸化式
a_k+1←-a_kc_k+a_k,
c_k+1←-c_k².
で繰り返すと, c_k=0の時, a_kが商になる.

これをそのまま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である. 目がちらちらするから, 途中のa_k, c_kを整数でシミュレートしてから同様な表にすると,

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とa_kとの差である. 最後のk=4の行は長いので, 少数点以下8桁目で折り返してある.

こういうのの性質をしらべるには, 各回におけるinvariantを発見するのが常道である. ためつすがめつするまでもなく, この例ではcとdeltaには比例関係があることが分かる. つまり, 上の計算ではdelta/cは-1/2, 下のではdelta/cは-4/5, つまり除算の商の符号を変えたものであった.

すなわち, n/d=a_k-(n/d)*c_kということだ.

ここまで分かると, 次は帰納法で証明するだけである.

k=0の時

a₀-(n/d)*c₀=n-(n/d)*(d-1)=n/d.

k=kで

a_k-(n/d)*c_k=n/d

が成立しているとする. その時k+1の式は

a_k+1-(a/d)*c_k+1
=-a_kc_k+a_k-(n/d)*(-c_k²)
=-a_kc_k+a_k+(n/d)*c_k²
=c_k(a_k-(n/d)*c_k)+a_k
=a_k-(n/d)*c_k)=n/d.

やはり, n/dになる. 一方 c_kはどんどん小さくなっていくから, c_k=0の時, a_k→n/d.

そういう風にEDSACの除算サブルーチンは出来ていたのであった.