Illiacのブートストラップ

アメリカイリノイ州のIllinois大学で1950年頃に作られた計算機がIlliacである.

この計算機は最初のプログラムの入れ方が面白かった. 最近 同僚が神田の古書店で, 電気通信研究所から放出された図書の山から「M1プログラム作製法」を見つけてきたというので, 見せてもらった.

早速プログラムを最初に入れるブートストラップのところを読んだ. 今回はそれをちょっと書いてみたい.

1959年頃の計算機は, 内部は二進法であっても, 機械語のプログラムは十進数で書く. 相対番地も使える. それを読み込み, 二進に直し, 相対番地を絶対番地に直す短いプログラムが用意してあった. イギリスCambridge大学のEDSAC計算機ではそれをイニシアルオーダーといった. アメリカIllinois大学のIlliacではD.O.I.(Decimal Order Input)といった. 東大のPC-1はEDSACがお手本だったのでやはりイニシアルオーダーといった. これらのプログラムは始めから二進法で作ってある.

これらの二進十進変換プログラムの読込み法がそれぞれの計算機でいろいろ工夫された. EDSACではイニシアルオーダーが電話交換機用のロータリースイッチに配線され, スタートボタンを押すとそれが回転して読み込まれる.

PC-1では二進法でパンチした紙テープの読込み方が決められ, それに従ってパンチされた.

Illiacの方法はこれらにくらべてユニークであった. ブートストラップという. つまり靴の紐のことで, 編み上げの靴紐をしばるように順々に完成していくのである.

Illiacは二進法40ビットの計算機である. 命令は1語の左と右に20ビットずつ2つ収める.

20ビットの命令は左8ビットが機能部, 右12ビットが番地部で, 機械語ではそれを十六進法で表す. ただしIlliacの十六進法は10,11,...がA,B,..ではなく, K,S,N,J,F,Lという不思議な文字であった. (king sized number just for loveと覚える.)

ブートストラップの命令は次のようだ.

0 80 028 40 001
1 80 028 40 002
2 19 026 26 000
1 80 028 40 000
0 L4 001 40 001
1 80 028 40 0F6

左端の0,1,2はその右の語が格納される番地である.

ブートストラップが始まるとシーケンスコントロールが0になり, 最初の80028 40 001が0番地と命令対レジスタに入る.

その次の行からがテープにパンチされている. 途中の空白や改行は見易いように入れてあるので, 実際にはない.

命令80はアドレスで示す数の1/4桁をアキュムレータに読めであり, 右命令の40はアキュムレータを番地へ格納せよである. 従って十六進の28は40だから10桁, つまり次の命令80 028 40 002が読み込まれ, 1番地へ格納される.

すると直ぐこの命令が実行されて, 19 026 26 000 が2番地に入る.

命令19はアキュムレータに1/2を置き, それを番地部の数だけ右シフトするので, 1/2は左端の符号ビットの次の位置が1で, 38ビット右シフトするから, アキュムレータの右端に1を置くことになる. 26は番地部の左命令にジャンプ.

そこで0番地の命令対を再び実行し, 80 028 40 000 が1番地に入り, 続いてこれを実行するからL4 001 40 001が0番地に入る.

L4は番地部の数をアキュムレータに加算する命令である.

従って1番地80 028 40 000は80 028 40 001になり, この命令により最後の80 028 40 0F6が1番地に入る. 現時点では

0 L4 001 40 001
1 80 028 40 0F6
2 19 026 26 000

となって2番地を実行しようとしている.

アキュムレータに1を置き0に戻り1番地の右の番地を1殖やして次の10桁を読み格納するというループが廻り始める. だから続く命令はF7, 通常の十六進ではe7₁₆=231₁₀から順に格納される.

M1の記憶装置は256語だったから, 231から255まで25語に50命令が入ることになる. Illiacは1024語だったかも知れないから, そうなら1番地右の定数が違っている.

IlliacのD.O.I.はCambridge大学から來ていたDavid Wheelerが作ったといわれる. このブートストラップもWheelerの考案かも知れない.

こういうブートストラップが作れたのは, Illiacに10桁読むという命令があったからで, EDSACやPC-1では1文字読む命令しかないのでこういうことはできない. EDSACそれも記憶場所に読み込むから, それをまたアキュムレータに読み出しで処理することになった. PC-1はどうせそうなるからというので, 読み込む先はアキュムレータであった.

最初のプログラムの読込み方は, その後ミニコン時代にも話題になった. 大方は正面のパネルにスイッチを使って始めの何語かを手でいれる方式であった.

ミニコンのメーカーからその入れ方が指示されるが, それにあきたりないユーザーは別の入力シーケンスを考え, どれだけ短くなったかを競い合った. 最近はそういう楽しみはなくなった.

数学体験館

この10月にオープンした東京理科大学の数学体験館を見た. パラメトロン計算機や微分解析機や沢山の機械式計算器のある近代科学資料館の地下に ある. 直接地下に行ける入口も作られた. 入場無料. 館長は秋山仁さん. 開館時間は12時(土曜は10時)〜16時. 日曜・月曜・祝日は休み.

実はそれを見るのではなく, 近代科学資料館の所用があっていったら, 体験館はまだ30分は開いていると聞き, いそいで見て回ったので, そのうち改めてゆっくり見学したい.

写真を撮影してはいけないというので, ここに示す図はちらしをスキャンしたものである.





ちらしの絵から面白そうなものを拾うと


これは(1/4)+(1/4)²+(1/4)³+...=1/3の証明である.

正三角形が4つの相似な正三角形に分割され, 上の1つが更に4つに分割され,...た図である. 例えば右下の青の正三角形は全体の1/4であり, 上の同じ大きさの正三角形は全体の1/4であって, その右下の青の正三角形はその1/4だから全体の(1/4)². このようにして青の正三角形は上の式の左辺になり, 青は赤や黄とともに3色で全体になるから, 青の面積は1/3というのである.

しかし展示は図だけなので, 証明は自分で考えなければならない.



これは最小公倍数・最大公約数算出器. 名前だけみると素晴しい計算機だが, 素因数分解は先に出来ていて, 例になっている90は2の球が1個, 3の球が2個, 5の球が1個で並んでおり, 24の方は2の球3個と3の球1個が並んでいる.

左右に分類器があって, 球をいれるとそれが分類されて出てくるから, それぞれの球の最大個数, 最小個数がわかれば公倍数や公約数が得られる.

これはどうもなぁという装置である.



球の体積の公式を導く話だ. 右にあるスイカをサッカーボールのように筋を入れて, 球の中心との母線で錐体に切り出す. 左にあるようにスイカの錐体が沢山できる. 底面は球の一部だがどんどん細い錐体にすると底面はほとんど平面になり, その錐体の体積は底面積に高さ(つまり球の半径r)を掛けた1/3である. 錐体全部の底面積はスイカの表面積だから4πr². 従ってスイカの体積は4πr³/3である.

表面積を知っている人なら体積も知っていそうなものではあるが.



三平方スライド. 円のなかに直角三角形と各辺を1辺とする正方形が出来ている. 左の図は小と中の正方形に青と黄の板が入っているところ. 円板を回転すると青と黄の板が分割されスライドして大の正方形にうまく収まることをしめす.

昔ボストンのサイエンスミュージアムに同じ趣向だが水が入っているのを見たことがあった.



二項分布パチンコ. あまりよくは見えぬが右上が坂の上になっていて, パチンコの玉が沢山入っている. 出口を開けると一斉に左下に流れ落ちるが, 途中に通路をランダムにする釘が打ってあり, それによってパチンコ玉は下の区画のどれかに落ちる. その形が二項分布に見えるという実験である.

なにかの本でも見たような気がする.

まぁこういう実験道具が沢山おいてあって, 数学がこういうふうに体験できると分かって面白い.

最初のちらしのフロアーマップにあるように, 立派な数学工房がある. 体験館の実験具の多くはどうもそこで試作されたらしい. お願いすると使わせても貰えるらしい. なにか作ってみるものはないだろうか.

計算機屋からみると, ユークリッドの互除法, エラトステネスのふるい, ハノイの塔のようなアルゴリズム関係の展示もあればいいのにと思うが, すでに展示場は一杯である.

微分解析機

前回の微分解析機のシミュレータの続きである.

実物の微分解析機には出力装置があり, 2つの変数軸をxとyに入れると解の図形が得られる.

MIT Schemeにあるscheme graphicsの最も簡単な使い方は

(define mydevice (make-graphics-device 'x))

で描画の装置を宣言する. 装置の描く場所の正方形は左右がx軸, 上下がy軸であり, 左端はx=-1, 右端はx=1, 下端はy=-1, 上端はy=1である. この座標系で(x₀,y₀)から(x₁,y₁) まで直線を引くには

(graphics-draw-line device x0 y0 x1 y1)

とする.

そこで前回のブログのsinz, coszを描くには

(define (draw device step x0 y0 x1 y1)
 (if (> step 0) (begin
  (graphics-draw-line device x0 y0
   (stream-car x1) (stream-car y1))
  (draw device (- step 1) (stream-car x1) (stream-car y1) 
    (stream-cdr x1) (stream-cdr y1)))))

を使う. (draw ..)が描画の関数で, deviceに対してstep回描画する. x0, y0は解の曲線の始点; x1, y1はstreamである.

一方 駆動する方は, 座標を描いた後

(for-each (lambda (y) (graphics-draw-line mydevice -1 y 1 y))
  '(-0.75 0 0.75))
(for-each (lambda (x) (graphics-draw-line mydevice x -1 x 1))
  '(-0.75 0 0.75))

(draw mydevice 6282 0 0.75
 (scale-stream sinz 0.75) (scale-stream cosz 0.75))

のようになっている. stepは積分がzの1/1000で進行し, 2πで一周するから6282ステップ. 次の0と0.75はxとyの初期値(つまり鉛筆を置く位置)で, 1だと画面すれすれになるので0.75倍にしてみた.

次のsinzとcoszが出力軸の変数で, こちらも0.75にギアダウンしている. そういうわけで, この呼び方も出力装置をシミュレートしているといえるであろう.

こうして描いたサークルテストの図が次である.



円の始点と終点の真上でちょっと食い違っているようにも見えるが....

次はe^zを描く. 描きたい範囲を-1と1の間に変換しなければならない. 下の図を見てほしい. 左と上が実際の座標で右と下がGraphicsの座標である. これによるとx軸は0.4倍すればよい. y軸の方は1が-0.2になっているから0.4倍してから0.6を引く.



従ってプログラムは次のようだ. 始めのfor-eachは座標を描く.

(for-each (lambda (y) (graphics-draw-line mydevice -1 y 1 y))
  '(-0.6 -0.2 0.2 0.6))
(for-each (lambda (x) (graphics-draw-line mydevice x -1 x 1))
  '(-0.8 -0.4 0 0.4 0.8))
 (scale-stream z 0.4) 
  (add-streams (scale-stream ez 0.4)
    (scale-stream ones -0.6)))

回数の1386は丁度上の端に相当するxの値がlog 4=1.386であることによる.



こういう図を描くと, e^-zで左半分も描きたくなる.

(define e-z (integrator (delay (scale-stream e-z -1))
  (delay z)
 1.0))

を用意し,

(draw mydevice 2500 0 -0.2
 (scale-stream z -0.4) 
  (add-streams (scale-stream e-z 0.4) 
   (scale-stream ones -0.6)))

を追加すれば, 左も描ける.

微分解析機

8月30日のブログに「これらのシミュレーションにはSchemeのstream処理が適していると思うが, まだ手が付かずにいる.」と書いた.

やっとどうやらシミュレーションが出来るようになったので, 今回はその話題である.

SICPの206ページのintegral関数のように披積分関数を遅延にするのが味噌である. SICPの例題では独立変数がdtだけだが, 微分解析機では「回す」値になにが入るかわからないので, こちらもstreamに対応させなければならない.

独立変数に方は回転角の差が必要になるので, stream-cadrからstream-carを引いたものに, 披積分関数のstream-carを掛けることになる.

まずintegratorは次のようだ.

(define (integrator delayed-integrand delayed-variable 
 initial-value)
 (define (mult-streams a b)
  (cons-stream (* (stream-car a)
            (- (stream-car (stream-cdr b)) (stream-car b)))
    (mult-streams (stream-cdr a) (stream-cdr b))))
 (define int
  (cons-stream initial-value
   (let ((integrand (force delayed-integrand))
         (variable (force delayed-variable)))
   (add-stream (mult-streams integrand variable) int))))
int)

integratorが貰う引数は遅延になっており, intの計算のなかでforceを使って強制する.

次はシミュレーションに必要な補助関数である. これらはSICPにあるものだ.

(define (scale-stream stream factor)
  (stream-map (lambda (x) (* x factor)) stream))
(define (add-streams s1 s2)
 (stream-map + s1 s2))
(define (integers-starting-from n)
 (cons-stream n (integers-starting-from (+ n 1))))
(define integers (integers-starting-from 1))

さてこれらが用意できたら, 8月14日のブログにあった関数からやってみる. ブログの図を見ながらプログラムを眺めてほしい.

最初はsin zとcos z. 独立変数はzである. zは整数のstreamで, 1/1000にして使う. sin zのπ/4, つまり45度での値をstream-refで見ることにする.

(define z (scale-stream integers 0.001))
(define sinz (integrator (delay cosz) (delay z) 0))
(define cosz (integrator (delay (scale-stream sinz -1))
  (delay z) 1.0))
(display (stream-ref sinz 785)) ;=>.7071024791302126

値としてはよさそうである.

次はe^z.

(define z (scale-stream integers 0.001))
(define ez (integrator (delay ez) (delay z) 1.0))
(stream-ref ez 1000) ;=>2.716923932235898

z²とz³は次のようにする.

(define z (scale-stream integers 0.001))
(define izdz (integrator (delay z) (delay z) 0))
(define z2 (scale-stream izdz 2))
(define iz2dz (integrator (delay z2) (delay z) 0))
(define z3 (scale-stream iz2dz 3))
(stream-ref z3 500) ;=>.12499950000000003
(stream-ref z2 500) ;=>.25049999999999994

0.5の3乗が0.125, 2乗が0.25ならまぁ問題はない.

tan zは多少手強いが,

(define ones (cons-stream 1 ones))
(define z (scale-stream integers 0.001))
(define tanz (integrator (delay 1+tanzsq) (delay z) 0))
(define itanzdtanz (integrator (delay tanz) (delay tanz) 0))
(define tanzsq (scale-stream itanzdtanz 2))
(define 1+tanzsq (add-streams tanzsq ones))
(stream-ref tanz 785) ;=>.9979528633948679
(stream-ref tanz (quotient 3142 6)) ;=>.5761873674775257

テストしたのはtan π/4とtan π/6で, 1/√ 3=.5773502691896258だからこれも合格である.

1/zとlog_ezはzを1から積分する. z=1での初期値はどちらも0である.

(define z (scale-stream integers 0.001))
(define logez (integrator (delay 1/z) (delay z) 0))
(define 1/z (integrator (delay (scale-stream 1/z -1)) 
   (delay logez) 1))
(stream-ref 1/z 500) ;=>.6664863143078377
(stream-ref logez 1718) ;=>.999948126092027

1/zの500での値は1/1.5なので0.6666.... log_ezの1718はlog_e2.718だから1でいいわけである. zの値は1からのはずだが, 回転角の差を使うから, integersで構わない.

微分解析機には結果を描画する出力装置もついているから, 次はscheme-graphicsを駆動するようなプログラムにしたい.

微分解析機

微分解析機でバックラッシュを減らす仕掛けには, フロントラッシュ以外にも, ラッシュロック(lashlock)なるものがあった.

微分解析機では独立変数軸で回転されるディスクは, 被積分関数軸の値によって前後に動かされる. つまりディスクの載っている台座が, 被積分関数軸に繋がる送りネジ(lead screw)に嵌めたナットと一緒になっていて, 移動するわけだ.

被積分関数の値は結果の積分値に大きく影響するので, この送りネジの工作精度は重要で, ケンブリッジの微分解析機では, 送りネジの作製に微分解析機全体の1/10のコストがかかったという.

ここでもバックラッシュを減らす方法が検討されている. 東京理科大学に保存されている微分解析機にもラッシュロックがちゃんとあった.

下の図でハッチのあるのが送りネジである. それに左の大きいナットと右の小さいナットが嵌めてあり, 大きいナットの上(破線の上)に台座が固定されている. だから送りネジはこのナットを経て台座を動かす.



2つのナットの間には, 太い線で示すバネが挟んであり, 兩ナットを離そうとしている.

小ナットには縁に何か所が切れ目があり, 大ナットから突出した棒の先の爪が差し込めて, 小ナットの回転を防いでいる.

つまり送りネジのネジ山を両端に押してバックラッシュをなくそうとしているのである. なるほどすごい仕掛けだ.

Crank本によると, この仕掛けもトルクアンプと同じく, ベツレヘムスチール社のNiemanが発明したそうだ.

微分解析機

微分解析機ではバスシャフト, クロスシャフトに積分機, トルクアンプ, 入力卓, 出力卓, 加算器などが接続してある. しかし佐々木本によると, 「ガタ補正装置」というものもあったらしい.

これは一体なにかと思う一方, 理科大の微分解析機になにやら不思議な機器がついていた. Crank本で判明したのだが, その不思議なものこそガタ補正装置であった. 英語では「フロントラッシュ」という.

微分解析機の最大の泣き所はトルクアンプでのバックラッシュであったらしい.

バックラッシュとは, 通常の平歯車ではこう起きる. 次の図のように2つの平歯車MとNがかみ合っていて, Mが駆動する側, Nが駆動される側とする.



Mが右回転していると, MとNのかみ合う部分では, Mの歯の下がNの歯を上から押し下げてNを左回転させる. Mの歯の下とNの歯の上は常に接しているから連続的に駆動できる.

Mが停止し逆回転を始めると, 今度はMの歯の上がNの歯の下を押し上げることになるが, こちらには隙間があって, MがNに当たるまでNは逆回転を始めずに停止している. 従ってMとNの回転角の関係を図に描くと, 回転が逆になると違う関係になる. これがバックラッシュである.

トルクアンプでも逆回転が始まると, 逆側のドラムの回りの糸が締まるまでは出力が始まらず, しばらく不感時間があってバックラッシュ現象が生じる.

このバックラッシュを出来るだけ早くキャンセルしようとするのがフロントラッシュである. 名前からして結構ふざけているがなかなかの機構である.

次がフロントラッシュのメカの分解図である. Cambridge大学のWilkes先生のところで撮った写真だと断りがあった.



フロントラッシュはバスシャフトの間に入れ, バスシャフトの回転数を加速する機能を持つ.

右が入力シャフトで, その左に両隣のシャフトへの引掛りをもったベルトに挟まれたドラムがある. ベルトはドラムが勝手に回らないよう軽くブレーキをかけている. ドラムに入力シャフトが通っているが, 入力シャフトとドラムはシャフトの回りではフリーになっている.

ドラムの左は遊星歯車. そして左端に内歯車とそれと一体になった出力シャフトがある.

ドラムには右方にペグが2本. この挿入位置は円周上に点在し, 場所は変更し調整出来る. 一方入力シャフトにはペグにぶつかるピンがある.

通常は入力シャフトがピンでドラムのペグを押し, 同時回転し, 従って遊星歯車も同時回転し, 入力シャフトと出力シャフトは同時回転する. しかし入力シャフトが逆転すると, ピンがペグから離れ, もう一方のペグに当るまで, ドラムはバンドの摩擦で停止する.

ドラムと入力シャフトの回転速度に差があると, 遊星歯車の出力は加速され, 出力シャフトも少し早く回り, やがてピンがもう一方のペグに当たって通常の回転数に戻るという仕掛けである.

この辺の事情を図にしたのが下だ.



左の図で横軸が入力, 縦軸が出力である. DからAは入力と出力が一体となってある方向に回転しているところとする. Aで入力が逆転し始めるても, 出力はすぐには追従しないからBまでは出力は動かぬ. BからCで出力は入力にすこし遅れて逆転する.

右の図はフロントラッシュの場合で, 横軸は入力シャフトのまだ先にある真の回転角である. Bまでは逆転直後で入力軸が遅れているところ. やがて遊星歯車が回転して出力は真の回転角に追いつく様子を示す. 逆側のC, Dでも事情は同じである.

なんとしてもバックラッシュを防ごうとう努力の後が見えるではないか.

微分解析機

初等関数を積分機を使って作りだす方法で, 前回のブログの続きである.

tan z

tan zはπ/2で無限大になるから要注意

d tan z/dz =1+tan²だから
tanz=∫ 1+tan²z

tan²のつもりのtanzsqなる変数を用意し

tanz = ∫ (1 + tanzsq) dz
tanzsq = 2∫ tanz d tanz

tanzは1+tanzsqで「動かし」, zで「回す」.
tanzsqはtanzで「動かし」, tanzで「回し」2倍する.

図は次のようだ.



3本のバスシャフトにまたがる[・Σ・]のような記号は加算器で, ・のバスシャフトの値の和がΣのバスシャフトに出力される.

1/z, log_ez

log_ez = ∫ 1/z dz
1/z = -∫ 1/z d log_ez だから

log_ezは1/zで「動かし」, zで「回す」.
1/zは-1/zで「動かし」, log_ezで「回す」.

図は次に示す.



これらのシミュレーションにはSchemeのstream処理が適していると思うが, まだ手が付かずにいる.