2012年1月26日木曜日

除算サブルーチン

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

イニシアルオーダーは十進二進変換もするが, その際の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の除算サブルーチンは出来ていたのであった.

2012年1月25日水曜日

世界一周プログラム

1月15日は私の恩師, 高橋秀俊先生の誕生日で, OBが先生のお宅に集まった. 先生が他界されたのは, 1985年だから, それから27年が経った. 成人の日は1月の第2月曜に変ったから, 15日ではなくなったが, 成人の日には今でもOBが某所に集まっている.

今回, その会で思わぬものを渡された. 高橋先生が東大にいらした頃の, パラメトロン計算機 PC-1 のライブラリのファイルである.



先生は1975年に東大を定年退官され, 慶応義塾大学へ移られた. 東大の教授室を後片付けした筑波大学の長谷部君が, そのファイルを見つけ, 自分の部屋で保管していたが, その長谷部君も定年になり, 筑波大学の教授室を片付けた. そしてそのファイルを私に委ねてきたのであった.

PC-1のライブラリを整備していた頃, 紙の質はまだ悪く, 日焼けで茶色になり, 縁もぼろぼろになってきた. 私が東大工学部にいた時, このままでは持たないと思って全部コピーしなおし, 製本した. その折, もう1部を作って高橋先生に差し上げたら大変喜ばれ, 「これを見ると, PC-1をやっていた頃が昨日のように思い出されます」という, 葉書の礼状を下さった.

先生がPC-1のライブラリファイルを, 東大に残していかれたとは, この正月まで知らなかった.

こうして手元にめぐってきたライブラリを見ると, 私の書いたイニシアルオーダーの前に, 先生が書かれたイニシアルオーダーがあった. 確かに1958年3月26日稼働開始の直後は先生のイニシアルオーダーが使われていたが, 長いし遅いので, 連休中に私が書き直したのが, その後1964年まで, 長く使われた.

先生のファイルには, "self reproducing program"もあった. これは計算機のメモリーをテストするため, 短いプログラムを走らせると, そのプログラムが自分を, アドレス部を修正しながら, 少し先のメモリーにコピーし, コピーが終わると, 新しい部分が走るようなものである.

Illinois大学のIlliac Iで使っていたと聞いていて, 高橋先生も書いてみられたのだろう.

Illiacのは, アドレス部がそれぞれその場所になっている, 2つのプログラムがあり, その差から新しいコピーを作っていると聞いたが, PC-1のは違っていた. (Illiacではこのプログラムをleapfrog(カエル跳び)と呼んでいたらしい.)

先生のスケッチは途中らしく, 正しくなさそうなので, 実際に使われていたプログラムを紹介しよう.



やはり, PC-1の命令から説明しなければなるまい. 1語は18ビットで, そのうち左端の6ビットが命令コード, 上のプログラムで, 各命令の先頭の"o"や"p"や"j"が入る. その右の1ビットが, l/sビットで, オペランドが長語(36ビット)か, 短語(18ビット)を区別し, 長語なら, 命令に"l"をつけ, そのビットが1になる. 残りの11ビットがアドレス部だが, PC-1は512語しかなかったので, 9ビットだけが必要, 上の2ビットは無視する.



0番地は"o"命令, outputの意味で, アキュムレータの左端の6ビットをテレタイプコードとして送り出す. アドレス部は, 命令の置いてある番地とすることになっているから, "0"である.

1番地は"p 16". "p"はpositive loadで, 16番地の短語をアキュムレータに取り出す. 従ってアキュムレータは, 16番地の命令"jl 17"になる.

2番地"jl"は, jump命令で, "j"の場合, lは長短ではなく, この場合は無条件ジャンプである. 従って, 3番地へ飛ぶ.

3番地"x 8"は, アキュムレータのアドレス部だけ, 8番地のアドレス部へ入れる命令. したがって, 8番地は "p 17"になる. 8番地のアドレスにかっこがついているのは, このアドレスは変更されるという意味である.

4番地の"a"はadd, 5番地の"s"はsubtractで, アキュムレータの"p 17"に, 16番地の"jl 17"を足し, 0番地の"o 0"を引く. 命令コード部も変るが今は関係なく, 17に17を足し0を引くから, アドレス部は17増えて34になり, 6番地, 7番地の"x"命令で, 9番地と12番地のアドレス部に入れられる.

このように, 定数を足すのにも, アドレス部は絶えず変るので, 定数を持っているわけにはいかず, 差を使う必要があった.

8番地に来て, "p 17", つまり"jl 3"をアキュムレータに持ってきて, それを9番地の"t", transfer命令で34番地にいれる. 10, 11番地は4, 5番地と同様で, アキュムレータのアドレス部を17増やし, 20にしてそれを34番地のアドレス部に入れる. 最終的に34番地は"jl 20"になる.

13番地は"p 8"だから, アキュムレータは"p 17"になる. 14, 15番地は, 14を足し, 15を引くから, 1引くことになり, アドレス部は16になり, 17番地へ飛び, さらに3番地へ飛ぶと, 今度は8番地, 9番地, 12番地のアドレス部がそれぞれ1減り, 16番地の"jl 17"が"jl 34"になって33番地に入る.

これを繰り返していると, やがて1番地の"p 16"が18番地に"p 33"として入り, 次のサイクルでは, "o 0"が"o 17"になって7番地の"jl 3"に上書きされる. コピーサイクルが17番地に到達した時点からコピーされたプログラムが走り出し, 同じことを繰り返すのである.

制御が新しいコピーに移ると, そこに出力命令があるので, プリンターが1文字打出す. つまりプリンターが文字を打ち続けている限り, メモリーは正常であるわけだ. 512/17=30だから, 30回のコピーでメモリーの始めに戻る. PC-1では, これを世界とみて, このプログラムを世界一周と呼んだ. 80日世界一周という映画があったから, ○○秒世界一周ともいったが, それが何秒だったかは, 残念ながら思い出せない.

2012年1月23日月曜日

素因数探し

昔は計算機はなくても, 計算用の道具はいろいろあった. 12月3日のブログに書いたfactor stencilもそういうものの1つである. 竹内端三先生はfactor stencilを「因数型紙」と訳された.

Derrick Norman Lehmerがfactor stencilを作ったのは1925年頃のことである. それ以前, 1914年にLehmerは10006721までの素数表を完成し出版した. といってみるだけは簡単であるが, 665000個の素数があるわけで, (Legendreの素数定理で計算すると(define (pi x) (/ x (- (log x) 1.08366))) (pi 10006721) => 665556.99...) 1ページに5000個ずつ, 133ページに収められている.

その第1ページ目には48593までの素数が縦100行, 横50列に印刷してある. factor stencilはこの第1ページの素数の, それぞれのRに対応する位置に穴を開けて作られた. Lehmerの当時の寄書に, "The device for cutting the stencils is already constructed and..."とあるから, 何かの器具を工夫したのであろうが, 詳しいことは分からない.

かくして作られた初版のfactor stencilには, しかし誤りが多数あったようで, 1937年にMichigan大学のJ.D.Elderがそれらを指摘するとともに, そういうものを作るならHollerithカードで作るのはどうかと提案した. Hollerithカードは, 今から50年くらい前, FortranプログラムをパンチしたいわゆるIBMカードで, 20世紀の初めの頃のアメリカの国勢調査の集計用にHerman Hollerithが開発した. 統計機械で読むため, 穴の位置は正確であり, パンチ用の機械も沢山あったはずだ.

Lermerもその案に賛成し, Elderに多くのノウハウを提供した. そういう次第でHollerithカード版のfactor stencilが完成出版されたらしい. 1枚のHollerithカードの, 縦10行横80列の位置を使い, Lehmerの5000個の代りに, 7枚の組で5600個の素数を収めたという. これが完成する直前にLehmerは他界した.

Hollerithカードのfactor stencilがどのようなものであったかに私は興味を持った. Lehmer流に素数を1から始めると1枚目のカードの最後の800番目は6131であり, そこまでの素数について, R=-1とR=2(12月3日のブログの表の1行目と2行目

3 5 7111317192329313741434753596167717379838997
-1 X X X X X X X X X X X
2 X X X X X X X X X X X

)について, 多分このようであったろうと想像して描いてみたのが下の図だ.




IBMカードは昔のプログラマにはお馴染だが, 最近は見たこともない人が多かろう. この図のように, 横187ミリ, 縦83ミリ程度のカードで, 向きが分かるように, 隅に1ヶ所にコーナーカットがある. 下の拡大図で分かるように, 0から9の行番号は大きい数字で, 0から79の列番号は小さい数字で示す. 本来のIBMカードの列番号は1から80であり, 行は0の上にXとYとがある.





列0, 1の各行に対する素数は

0 1 2 3 4 5 6 7 8 9
列0 1 2 3 5 7 11 13 17 19 23
列1 29 31 37 41 43 47 53 59 61 67

である. 列2以降に対しても, 素数表があれば, 穴の位置と素数は対応づけられる. 上の表と図の穴の位置とを較べてみて欲しい.

こういう絵を描くと, 実際に穴のあちあカードを作ってみたくなる. 鎌倉のFabLabの田中君と相談し, craftroboを使って試作したR=2のカードが下だ. まだ調整の必要があるらしく, うまく切れていない穴もあるが, 感触は得られた.



プログラムのカードと違い, 穴の数が多くてバイナリカードのようであり, 丁寧に扱わないとすぐに切れそうなのが心配だ. Edlerの作ったカードのセットはどこかに残っていないかしら.