開平法

CurtaのウェブページにDibble-Dabble法という平方根の アルゴリズムがあった.

 2. Dibble-dabble method

2.A Shift the carriage full CCW and zeroize.

2.B Enter "1" in the most significant input slider (IS).

2.C Add ONCE. If the value in the total register exceeds the 
target number, then subtract out the last number added and go 
to step F.

2.D Increase the value in the currently active IS by mentally 
adding "2". If adding "2" to "9" change the "9" to "1" and 
add "1" the the next left input slider (not relevant first 
time through).

2.E Repeat steps C & D.

2.F Reduce the setting of the current IS by "1" and shift the 
carriage CW one step. If the limit of carriage travel is 
reached, go to step I.

2.G Enter "1" in the next most significant IS (the one to the 
right of the one previously used).

2.H Repeat steps C through G.

2.I The number in the "turns counter" register is the true 
square root.

やってみよう. その前に日本語にすると

2.A キャリージを反時計回り(counter clock wise)に一杯に 回し, 零にする(zeroize).

2.B 置数レジスタ(input slider)の最上位に1を置く. 2.C 一度加算する. 結果レジスタ(total register)の値が, 目標の数値を越えるなら, 最後に足した数を引き, Fへ進む.

2.D 置数レジスタを2増やす. 2を9に足す時は, 9を1にし, 置数レジスタの左の 桁に1足す.

2.E CとDのステップを繰り返す.

2.F 置数レジスタの値を1減らし, キャリージを時計回りに1ステップ回す. キャリージ が回らないならIへ進む.

2.G 置数レジスタの次の最上位(前回の桁の右の桁)に1を置く.

2.H CからGのステップを繰り返す.

2.I 回転レジスタの数が平方根になる.

ここに Curtaの精巧なシミュレータがあるからそれを使って2の平方根を求めてみよう.

Curta計算機の部分の名称は下の図を参照してほしい.



Curtaの操作の基本は次の通り.

加減する数は左の絵の下の置数レジスタにつまみを上下して入れる.

加減算の結果は右の絵の結果レジスタに出る. レジスタへの加減する位置は キャリージを回して決める. それにはキャリージの左右の矢印をクリックする.

減算の時は左の絵の上の操作ハンドルを上に引き上げ, 加算の時は 押し下げ, 右の絵の操作ハンドルを回す.

結果レジスタや回転レジスタは右の絵のクリアレバーを回す.

そこで2の平方根を小数点以下5桁求めるために20000000000の平方根を 計算する. (下の説明の行末の番号は後のトレースの行番号.)

まずキャリージの右の矢印を5回押してキャリージを回す.

置数レジスタの6の桁を1にして加算する. 
結果レジスタは10000000000, 回転レジスタは100000になる. (1)
置数レジスタの6の桁を3にして加算する. 
結果レジスタは40000000000, 回転レジスタは200000になる. (2)
40000000000は20000000000より大きいから減算する. (3)
6の桁を1減らして2にし, キャリージを1桁戻す.

置数レジスタの5の桁を1にして(置数レジスタは210000になる)加算
する.(4)
結果レジスタは12100000000, 回転レジスタは110000になる.
置数レジスタの5の桁を3にして加算する. 5にして, 7にして, 9にし
て加算する.
結果レジスタは22500000000, 回転レジスタは150000になる.(8)
20000000000より大きいから減算する.(9)
5の桁を1減らして8にし, キャリージを1桁戻す.

置数レジスタの4の桁を1にして(置数レジスタは281000になる)加算す
る.(10)
3にして加算する.(11)
結果レジスタは20164000000, 回転レジスタは142000になる.
20000000000より大きいから減算する.(12)
4の桁を1減らして2にし, キャリージを1桁戻す.

置数レジスタの3の桁を1にして(置数レジスタは282100になる)加算す
る.
9まで加算すると(17)
結果レジスタは20022250000, 回転レジスタは141500になる.

20000000000より大きいから減算する.(18)
3の桁を1減らして8にし, キャリージを1桁戻す.

この辺までで大体分かったからプログラムを書いてみる.

(define (dibbledabble radicand)
 (let ((s0 1) (s1 1) (s 0) (a 0) (c 0))
  (define (loopp) (display "+") (display s) (display " ")
   (set! a (+ a s)) (set! c (+ c s1)) (display a)
   (display " ") (display c) (newline)
   (if (> a radicand)
    (begin (display "-") (display s) (display " ") 
     (set! a (- a s)) (set! c (- c s1)) (display a)
     (display " ") (display c) (newline)
     (set! s (/ (- s s0) 10))
     (if (> s0 1)
      (begin (set! s0 (/ s0 100)) (set! s1 (/ s1 10)) 
       (set! s (+ s s0)) (loopp)) c))
    (begin (set! s (+ s s0 s0)) (loopp))))
  (do ((n 100 (* n 100))) ((> n radicand))
   (set! s0 (* s0 100)) (set! s1 (* s1 10)))
  (display radicand) (display " ") (display s0)
  (display " ") (display s1) (newline) (set! s s0)
  (loopp)))

(dibbledabble 20000000000)

上のプログラムでletで用意したs0は置数レジスタを増減する数, s1は回転 カウンタを増減する数, sは結果レジスタに加減する数, aは結果レジスタ, cは回転カウンタ, radicandは平方根をとる数である.

radicandを20000000000とした時, まずs0とs1を用意する. それが本体の do loopである. s0は10000000000, s1は100000になる.

上の方のloopが1桁分の計算で, aにsを, cにs1を足す. (if (> a radicand) は結果レジスタが目標を超えた場合で, 置数レジスタと回転レジスタを 戻し, sからも1引いてキャリージを回すように1/10にする. またs0を1/100, s1を1/10にする.

超えない時はsにs0を2回足して加算を繰り返す.

このプログラムを実行した結果が次だ. 左端の斜体の行番号は後から追加した. 0行目はtarget, s0, s1を示す. 後の行はsを足したか引いたかを+, -で示し, s, a, cを順に示している.

0 20000000000 10000000000 100000
1+10000000000 10000000000 100000
2+30000000000 40000000000 200000
3-30000000000 10000000000 100000
4+2100000000 12100000000 110000
5+2300000000 14400000000 120000
6+2500000000 16900000000 130000
7+2700000000 19600000000 140000
8+2900000000 22500000000 150000
9-2900000000 19600000000 140000
10+281000000 19881000000 141000
11+283000000 20164000000 142000
12-283000000 19881000000 141000
13+28210000 19909210000 141100
14+28230000 19937440000 141200
15+28250000 19965690000 141300
16+28270000 19993960000 141400
17+28290000 20022250000 141500
18-28290000 19993960000 141400
19+2828100 19996788100 141410
20+2828300 19999616400 141420
21+2828500 20002444900 141430
22-2828500 19999616400 141420
23+282841 19999899241 141421
24+282843 20000182084 141422
25-282843 19999899241 141421

この1, 3, 5,...と引くのは要するに手回しの機械式計算機の平方根の計算法で あって, それがdibble-dabbleという名前だとは知らなかった. 以前のブログの 5倍してから5, 15, 25,...と引くのは電動計算機の方法というべきであったかも しれない.

Christopher StracheyのGPM

Christopher Stracheyという名前を初めて見たのは1959年にパリで開催された第1回のIFIPの会議の報告書で, そのpp.336-341に彼の
Time sharing in large, fast computers
という論文があった.

要するに割込みの機能を使ってのマルチプログラミングの提案であるが, 当時我々東大物理の高橋研のパラメトロン計算機でも, すでに割込みと並列計算の実験を始めていたので, 同じようなことは, 洋の東西で同じ頃に気附くものだと思った.

そのStracheyに出逢ったのは, 1971年4月 英国Warwick大学で開かれたIFIP WG2.2の会議でであった. この会議には主査のMike Woodger(NPL)を始め, Edsger Dijkstra, Willem van der Poel, Simulaを開発したOle-Johan Dahlなど有名な計算機科学者が大勢いた. その会議中 パーティーが主催者のJohn Buxtonの超古い家で開催され, みなが車に分乗して行くことになった. その時昔の小さいMiniに乗せてくれたのが大柄なStracheyであった. 「Eiitiは後に入れ」と私を後部座席に押し込み, 助手席に乘ったBrian Randellと機関銃のような英語で話し合っていた.

さて最近のIEEE Annals of the History of Computing(July-September 2013)を眺めていたら, 何年か前に私が情報処理学会誌に「Wilkes先生を悼む」を執筆した時にWilkes先生の写真を 送ってくれたDavid HartleyさんがCPLのことを書いていて, その参考文献を辿っているうちに, CPLのアセンブラをGPM(General Purpose Macrogenerator)で作ったという話から, GPMを開発したStracheyがComputer Journalに寄せた論文
C. Strachey, A general purpose macrogenerator, Computer Journal Vol.8, No.3, 225-241
を探しあてた.

CPLはCambridge Plus Londonとも思えるが, Combined Programming Languageのアクロニムである. この言語はその後Martin RichardsがMITにいるころ, BCPLに単純化して開発し, さらにベル研究所に行ってC言語になった, ということを覚えている人はもう少ないに違いない.

ところでそのmacrogeneratorである.

macroというからにはマクロ定義とマクロ呼出しがあるわけで, aというマクロをb~1dと定義すると, aがマクロ名, b~1dがマクロ本体で, 本体のうちbとdは定数, ~1は第1引数である. 一方, $a,c;がマクロ呼出しで, (もとの論文では先頭は§であるがASCII文字にないので, 以下では$を使う.) マクロ名の本体に対して, cが第1実引数となって, 呼出しの結果はbcdとなる.

マクロ呼出しは $マクロ名,引数1,引数2,...; の形. このマクロ名は引数0としても使うことが出来る.

マクロ呼出しは $ と ; で, 定数は < と > で囲まれているから, マクロ名や引数の中に再帰的に書くことが出来る. その場合, 局所的なマクロ呼出しは呼出しの結果, 局所的な定数は, 一番外の < と > を外した定数本体でその部分を置き換える. 置き換えて出来た結果の文字列で外側のマクロを呼び出す.

GPMではマクロ定義もマクロ呼出しの形で, 上のマクロaの定義は

$def,a,<b~1d>;

のように書く. このマクロ呼出しは空の文字列を結果として返す.

マクロ本体は, この例のように, 本体を評価されたくない時は < と > で囲む.

マクロ呼出しの入れ子の例は次の通り:

引数にマクロ呼出しがある例

$def,a,<b~1d>;

と定義し,

$a,$a,c;;

と呼び出すと,

bbcdd

が返る.

マクロ名にマクロ呼出しがある例

$def,a,<b~1d>;$def,bcd,<b~1c~2d>;

と定義し

$$a,c;,e,f;

と呼び出すと

bdcfd

が返る.

Stracheyの論文にはsucとsuccessorというマクロの例がある.

$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$def,successor,
 <$~2,$def,~2,~1<,$suc,>~2<;>;$def,9,<$suc,>~1<;,0>;;>;

sucは$suc,4;のように呼び出す. すると本体の1というマクロを呼び出すが, 引数の評価中にマクロ1が定義される. 引数~1が4なのでこのマクロの本体は~4になり, 1のマクロから5がとられて4の次の5が返る.

successorは 2桁の数の次を返す. $successor,3,4;だと3,5が, $successor,2,9;だと3,0が返るものだ.

3,4の場合はsuccessorの本体が

$4,$def,4,3,$suc,4;;$def,9.$suc,3;,0;;

となりマクロ4を呼び出そうとする. その$に対応する;は一番右にあるので, その前にマクロ4と9を定義する.

$def,4,3,$suc,4;;

だからマクロ4は3,5になり

$def,9,$suc,3;,0:

だからマクロ9は4,0になる. この状態でマクロ4を呼び出すから3,5か返る.

1の桁が9だと, マクロ9が2回定義され, その状態でマクロ9を呼び出すから後で定義した繰り上げのあるマクロの結果が返るのである.

この例にあるように, 本体評価中に定義されたものは, 評価が終わると定義のリストから削除されることに注意しよう.

Stracheyの元の論文にはCPLで書いたGPMの実装が載っている. CPLには再帰呼出しなどないから, スタックの上にリンクをつなげたりしていて, あまり読みたいとも思わない.

以下は私がSchemeで実装したものである. StracheyのGPMには四則演算の機能もあるが, そういうのはGPMの本質的なものではないから, 割愛した.

(define ch '())
(define (gpm str env args)
 (let ((outs "") (index 0))
  (define (getch)
   (set! ch (string-ref str index))
   (set! index (+ index 1))
   (if (or (char=? ch #\space) (char=? ch #\newline)) (getch)))
  (define (readquote)
   (let ((outs ""))
    (define (rq)
      (getch)
      (cond ((char=? ch #\>) outs)
            ((char=? ch #\<)
             (set! outs
              (string-append outs "<" (readquote) ">"))
             (rq))
            (else
             (set! outs (string-append outs (string ch)))
             (rq))))
    (rq)))
  (define (readstring)
   (let ((outs ""))
    (define (rs)
     (if (= index (string-length str)) outs
      (begin (getch)
       (cond ((char=? ch #\,) (cons ch outs))
             ((char=? ch #\;) (cons ch outs))
             ((char=? ch #\<)
              (set! outs (string-append outs (readquote)))
              (rs))
             ((char=? ch #\~)
              (getch)
              (set! outs (string-append outs
               (list-ref args (- (char->integer ch) 48))))
              (rs))
             ((char=? ch #\$)
              (set! outs (string-append outs (readmacrocall)))
              (rs))
             (else
              (set! outs (string-append outs (string ch)))
              (rs))))))
    (rs)))
  (define (readmacrocall)
    (let ((actuals '()))
     (define (rm)
      (let* ((result (readstring))
             (ch (car result)) (s (cdr result)))
       (cond ((char=? ch #\,)
              (set! actuals (cons s actuals)) (rm))
             ((char=? ch #\;)
              (set! actuals (cons s actuals))
              (macrocall (reverse actuals))))))
     (rm)))
  (define (macrocall actuals)
   (cond ((string=? (car actuals) "def")
          (set! env (cons (cdr actuals) env)) "")
         (else
          (let ((def (assoc (car actuals) env)))
           (apply (cons (cadr def) actuals))))))
  (define (apply expr)
   (gpm (car expr) env (cdr expr)))
  (define (gloop)
   (let ((result (readstring)))
    (cond ((string? result) (string-append outs result))
          ((char=? (car result) #\,)
           (set! outs (string-append outs (cdr result) ","))
           (gloop)))))

  (gloop)))

上の例を実行してみよう.

(define str "$def,a,;$a,c;")
(gpm str '() '()) => bcd

(define str "$def,a,;$a,$a,c;;")
(gpm str '() '()) => bbcdd

(define str "$def,a,;$def,bcd,;$$a,c;,e,f;")
(gpm str '() '()) => becfd

(define str "$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$suc,4;")
(gpm str '() '()) => 5

(define str "$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$def,successor,
 <$~2,$def,~2,~1<,$suc,>~2<;>;$def,9,<$suc,>~1<;,0>;;>;
$successor,3,4;,$successor,2,9;")
(gpm str '() '()) => 3,5,3,0

ついでに, 論文にあったsumは

$sum,α,β,γ;

で2桁の十進数α,βに1桁の数γを足すもので, $3,4,2;は以下のように3,6になる.

(define str "$def,suc,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$def,successor,
 <$~2,$def,~2,~1<,$suc,>~2<;>;$def,9,<$suc,>~1<;,0>;;>;
$def,sum,<$s,~1,~2,0,
 $def,s,<$~3,$def,~3,<$s,>$successor,~1,~2;
  <,>$suc,~3;<;>;$def,>~3<,~1<,>~2;;>;;>;
$sum,3,4,2;")
(gpm str '() '()) => 3,6

これは難しいからマクロ呼出しのトレースをしてみる.

00 (def suc $1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;)
01 (def successor $~2,$def,~2,~1<,$suc,>~2<;>;
  $def,9,<$suc,>~1<;,0>;;)
02 (def sum $s,~1,~2,0,$def,s,<$~3,$def,~3,<$s,>
  $successor,~1,~2;<,>$suc,~3;<;>;$def,>~3<,~1<,>~2;;>;;)
03 (sum 3 4 2)
04 (def s $~3,$def,~3,<$s,>$successor,~1,~2;<,>
  $suc,~3;<;>;$def,2,~1<,>~2;;)
05 (s 3 4 0 )
06 (successor 3 4)
07 (def 4 3,$suc,4;)
08 (def 9 $suc,3;,0)
09 (4 )
10 (suc 4)
11 (def 1 ~4)
12 (1 2 3 4 5 6 7 8 9 10 )
13 (suc 0)
14 (def 1 ~0)
15 (1 2 3 4 5 6 7 8 9 10 )
16 (def 0 $s,3,5,1;)
17 (def 2 3,4)
18 (0 )
19 (s 3 5 1)
20 (successor 3 5)
21 (def 5 3,$suc,5;)
22 (def 9 $suc,3;,0)
23 (5 )
24 (suc 5)
25 (def 1 ~5)
26 (1 2 3 4 5 6 7 8 9 10 )
27 (suc 1)
28 (def 1 ~1)
29 (1 2 3 4 5 6 7 8 9 10 )
30 (def 1 $s,3,6,2;)
31 (def 2 3,5)
32 (1 )
33 (s 3 6 2)
34 (successor 3 6)
35 (def 6 3,$suc,6;)
36 (def 9 $suc,3;,0)
37 (6 )
38 (suc 6)
39 (def 1 ~6)
40 (1 2 3 4 5 6 7 8 9 10 )
41 (suc 2)
42 (def 1 ~2)
43 (1 2 3 4 5 6 7 8 9 10 )
44 (def 2 $s,3,7,3;)
45 (def 2 3,6)
46 (2 )

行0,1,2はマクロ定義, 3が呼出しだ. sumを呼び出すと$s,3,4,0とマクロを呼び出す準備をするが ;はまだないので$def,s,...をやってsを定義する.

この際殆どの引数は< >に囲まれているが, 右から15文字目くらいの~3はsumの最後の引数2になる. 4行目で見るとおり. ここでsを呼び出す($s,3,4,0;).

sの本体は$0の呼出しをしようとするが$def,0,があるので0を定義する(16行目). 3,5は$successor,~1,~2;で, 1は$suc,~3;で作られる. それを作るのがトレースの15行目までだ. 次にマクロ2を定義する(17行目).

やっと0が呼び出せ(18行目), Ss,3,5,1;を呼び出す(19行目).

今度はマクロ1と2を同じように定義し, $s,3,6,2;を呼び出す(33行目).

するとマクロ2と2が定義され(44, 45行目), 後から定義された2の本体3,6が返るのである.

万歳. うまくいっているようだ.

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が発明したそうだ.