2014年3月22日土曜日

ビットごとの秘法と技法 から

最も左のビットの位置を知る法

二進法で表したxの最も左の1のビットの位置λ(x)は, 2を底とする対数lgがあれば簡単だ.

λ(x)=⌊lg (x)⌋

たしかに
 
(lg 1) => 0
(lg 2) => 1.
(lg 3) => 1.5849625007211563
(lg 4) => 2.
(lg 5) => 2.321928094887362
(lg 6) => 2.584962500721156
(lg 7) => 2.807354922057604
(lg 8) => 3.
になっている.

TAOCPでの最初の方法は浮動小数点演算命令による.
 FLOTU y,ROUND_DOWN,x; SUB y,y,fone; SR lam,y,52
ここで fone=#3ff0000000000000.

これでうまく行く理由だが, MMIXの浮動小数点はIEEE/ANSI standard 754で, その形は



Sは1ビット, Eは11ビット, Fは52ビットの整数で, S=0なら正, =1なら負. 0<E<2047なら
±2E-1023(1+F/252)
を表す.

だからx=1ならFの部分は0で, Eの部分は1023; x=2ならFの部分 は0で, Eの部分は1024; のようになる. 従って上の命令のように, FOUTU(convert fixed to floating unsigned) で変換し, Eの部分から1023を引いて52ビット右シフトしたのがλになるわけだ.

一方, 浮動小数点を使わないMMIX流の方法は次の通り.
 SRU y,x,32; ZSNZ lam,y,32;
 ADD t,lam,16; SRU y,x,t; CSNZ lam,y,t;
 ADD t,lam,8; SRU y,x,t; CSNZ lam,y,t;
 SRU y,x,lam; LDB t,lamtab,y; ADD lam,lam,t;
まずxを32ビット右シフトしてyに置く.

ZSNZ lam,y,32はyが≠0なら(つまりxの左半分に1があれば)lamに32を, そうでないなら0を入れる. 次はそのlamに16を足してtに置いておく. lamは16か48である. xを16か48ビット右シフトしてyに置く. そのyが≠0ならlamにtを入れ, そうでないならlamはそのまま. 次は8ビットの範囲で同様なことをするので, lamは最も左の1を含む8ビットの範囲の右の境界の値になっている.

SRU y,x,lamでその範囲をレジスタの右端に移動し, 256バイトの表lamtabを引いてlamに足すのである. lamtabは0の場所は使わない. 1のところは0, 2,3のところは1, 4〜7は2, 8〜15は3, ..., 2k〜2k+1-1はk, 128〜255は7になっている.

最も右のビットを取り出すには x & -x とやったが, 左に対してはうまい方法はない. 「ハッカーのたのしみ」の著者 Warrenの考えた方法というのはこうだ.

yxとする. レジスタの長さを2dとして0 ≤ k <dについて yy | (y » 2k)とする. y - (y » 1) がxの最も左ビットである.

yにあるxの最も左のビットは次々とシフトされて右方向に2ビット, 4ビット, 8ビットと広がり, レジスタの右端まで1で埋める(塗りつぶす). 他の右の方にある1はこの塗りつぶしに埋没する. 左端の0の並びはそのままである. そこでyからyを1ビット右にシフトしたものを引くと, xの最も左のビットがとれるわけだ.

TAOCPには

λx = λy   if and only if    xyx & y

という式がある. これも証明は簡単だ.

x, y の最も左のビットの位置が同じとすると, xyのその位置は0になり, x & yのその位置は1になるから, ⊕の方が小さい. 最も左の1の位置が異なると, より左にある1の位置では, ⊕は1, &は0なので⊕の方が大きいことになる.

気になるのは等号の場合だ. x, yが共に0の時, λxもλyも不定になる. 不定同士を=にしたいとすると, ⊕も&も0だからこの場合は等号になるのであろう.

2014年3月21日金曜日

ビットごとの秘法と技法 から

最も右のビットの位置を知る法

TAOCPにある第二の方法はde Bruijnサイクルを利用するものだ.

4ビットのde Bruijnサイクルは例えば次のようなものである.

0000111101001011
サイクルというから円状に描いたのがこの図で, 11時の方向にある0から時計回りに上の数字が並んでいる. 各数字から右方向へ4文字で0から15までのいずれかの二進数になっている.


内側の色の付いた円弧は, その上のビットの二進数が下の同じ色の十進数であることを示す.

この作り方は次の通り.

まず0000と書く. つまり0だ. この4ビットの左端を削除し, 右端に0か1を挿入すると0000と0001が出来る. つまり0の次は0と1だ.

1000からも0と1が出来る. 8の次も0と1だ.

要するに0から15までの数について, 8以上なら8を引き, 8未満ならそのままで, それを2倍したものと2倍して1を足したものが次として出来るわけだ.

0 → 0,18 → 0,1
1 → 2,39 → 2,3
2 → 4,510 → 4,5
3 → 6,711 → 6,7
4 → 8,912 → 8,9
5 → 10,1113 → 10,11
6 → 12,1314 → 12,13
7 → 14,1515 → 14,15
この遷移を有向グラフにしたのが次の図である. それぞれのノードは入次数も出次数も2の正則グラフである.


このグラフですべてのノードを経由して元へ戻るHamilton経路の一例が赤い矢印で, それが上にあったde Bruijnサイクルである. この例は 0000 の右に 1111 を置き, 0100 の後に反対の 1011 を置いた形になっている.

4ビットのde Bruijnサイクルがどのくらいあるか, 手元に計算機があるから全解探索をしてみた. その結果が下の図で16通りあった. 上の目の子で探した解は最下段の左から2つ目である.


そのそれぞれに対するサイクルは次のようだ.

0000100110101111
0000100111101011
0000101001101111
0000101001111011
0000101100111101
0000101101001111
0000101111001101
0000101111010011
0000110010111101
0000110100101111
0000110101111001
0000110111100101
0000111100101101
0000111101001011
0000111101011001
0000111101100101

そこでρ関数の話になる. 最も右の1のビットにするのは常套手段x & -xで, それにde Bruijnサイクルの定数を掛け, 64ビットの左端から6ビットを取ると, 最初の1の位置により, 0から63のいずれかが得られる.

それを先程の4ビットの例でいうと, サイクルは0000111101001011だったから, 1を掛けた時は0000が, 2を掛けた時は1ビット左にシフトしたのの左端4ビットだから0001, 4を掛けた時は0011, ... のように得られる. その様子を次の図で示す.



横長の枠がde Bruijnサイクルを何ビットか左シフトした位置である. 16ビットのレジスタの左端の4ビットに相当する場所にだけ二進数が書いてある. 下の方, 213, 214, 215を掛ける辺りはサイクルの右端の4ビットがレジスタからはみ出しているが, サイクルの左端が0000だったのが幸いしてちょうど輪になっている.

さてこれから分るように

×1は20を掛けたのだから, 表の0番には0を置く.
×2は21を掛けたのだから, 表の1番には1を置く.
×4は22を掛けたのだから, 表の3番には2を置く.
×8は23を掛けたのだから, 表の7番には3を置く.
... とやって出来上がった表が

0, 1, 10, 2, 8, 11, 13, 3, 15, 9, 7, 12, 14, 6, 5, 4

であり, ρ関数は最も右端の1にde Bruijnサイクルの定数を掛け, 左端の4ビットの位置を表で見ると得られるのである.

TAOCPの記述はこうだ.



TAOCPのレジスタは64ビットだから6ビットのde Bruijnサイクルを使う. それは#03f79d71b4ca8b09である. やはり左端が000000で始まっていることに注意しよう. まず6ビットずつ取ると0から63が得られることを確認しよう.
(define as
(map (lambda (n) 
 (quotient (modulo (* #x03f79d71b4ca8b09 (expt 2 n)) 
  (expt 2 64)) (expt 2 58)))
 (a2b 0 64)))

as =>
(0 1 3 7 15 31 63 62 61 59 55 47 30 60 57 51 39 14 29 58 53
43 23 46 28 56 49 35 6 13 27 54 45 26 52 41 19 38 12 25 50 37
10 21 42 20 40 17 34 5 11 22 44 24 48 33 2 4 9 18 36 8 16 32)
ソートして0から63が1回ずつなことを確かめる.
(sort as <)
=>
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
63)
decodeで使う表を作る.
(define bs 
(map (lambda (n) (- 64 (length (member n as)))) (a2b 0 64)))

bs
=>
(0 1 56 2 57 49 28 3 61 58 42 50 38 29 17 4 62 47 59 36 45 43
51 22 53 39 33 30 24 18 12 5 63 55 48 27 60 41 37 16 46 35 44
21 52 32 23 11 54 26 40 15 34 20 31 10 25 14 19 9 13 8 7 6)
もちろんTAOCPにある表と合っている.

2014年3月19日水曜日

ビットごとの秘法と技法 から

最も右のビットの位置を知る法

TAOCP 7.1.3はbitwise tricks and techniques(ビットごとの秘法と技法)で楽しい話題が満載だ. このブログで以前取り上げたビットスワップもそこにある.

今回の話は計算機にある語xの最も右にある1のビットの位置を計算するもので, x = 0なら1のビットはないから, エラーか不定とするか無限大とするかだが, その辺はどうでもいいのでx≠0 の場合を考える.

xに対してこの位置をρ(x)で表すので, 位置を求めるアルゴリズムをρ関数という. ρは右(right)のRに対応するギリシア語のアルファベット(小文字)である. これに対し最も左のビットはTAOCPではleftのLのギリシア語のλという. λはLisp屋にはちょっと困る命名だ.

二進法での計算に慣れている人は, 最も右の1のビットを取り出すのが簡単なことは知っている. x & -xでそのビットが得られる. x = 0ならどのビットも立たず0である.

例えばxが十進法の10なら, 二進法では1010, -10は二進法では ...0110(左の方は1が何桁も並んでいる)だからandをとると...0010となる. 二進法のビットを右からa0, a1,...と番号をつけると, ...0010の1はa1だから位置としては1とする. ρ(10)=1だ.

奇数は右端が1なのでρ(奇数)=0である.

ρ関数は原始的には右端のビットが1になるまでxを右シフトすればよい.

(define (rho x)
 (do ((i 0 (+ i 1))) ((odd? x) i)
  (set! x (quotient x 2))))

(map rho (a2b 1 17))
=> (0 1 0 2 0 1 0 3 0 1 0 2 0 1 0 4)
x & -xで最も右のビットだけが得られたら, その2を底とする対数をとるという方法も考えられる. TAOCPでは2を底とする対数をlgと書くから,
(define (lg x)
 (/ (log x) (log 2)))
2の63乗までのlgを取ってみると
(do ((i 0 (+ i 1))) ((= i 64))
 (display (lg (expt 2 i))) (newline))
0
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.000000000000004
30.
31.000000000000004
32.
33.
34.
35.
36.
37.
38.
39.00000000000001
40.
41.
42.
43.
44.
45.
46.
47.00000000000001
48.
49.
50.
51.00000000000001
52.
53.
54.
55.00000000000001
56.
57.
58.00000000000001
59.00000000000001
60.
61.
62.00000000000001
63.
ところどころに誤差が出るからroundしinexact->exactすると
(define (rho x)
 (inexact->exact (round (/ (log x) (log 2)))))

(map (lambda (i) (rho (expt 2 i))) (a2b 0 64))
=>
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
 63)
しかしこれらはアルゴリズムとしてはお粗末で, いろいろな工夫がされている. 例えば語の中の1を数える命令(sideway addition SADD)を計算機が持つているなら, x & -xで右端の1だけが取れた後, 1を引くと右端の1は0になり, その右の0は1になるから1の数を数えればよい.

MMIXのSADDはこのためにあるような命令で,
 SUBU t,x,1; SADD rho,t,x
SUBU でxから1を引き, SADDはtのビットが1で, xのビットが0のものの和をとるからrhoに答えが得られる. x = 0の時は, 1を引くとtがオール1になり, xが0だから64ビットの1が足されて64が得られる. TAOCPにはこれを無限大とみてよいと書いてあった.

SADDのような奇妙な命令がない場合のTAOCPの方法はこうだ.

MMIXというアセンブリ言語で書いてあるがそれを写すと
m0 GREG #5555555555555555 ;m1 GREG #3333333333333333;
m2 GREG #0f0f0f0f0f0f0f0f ;m3 GREG #00ff00ff00ff00ff;
m4 GREG #0000ffff0000ffff ;m5 GREG #0000ffff0000ffff;
 NEGU y,x; AND y,x,y; AND q,y,m5; ZSZ rho,q 32;
 AND q,y,m4; ADD t,rho,16; CSZ rho,q,t;
 AND q,y,m3; ADD t,rho,8; CSZ rho,q,t;
 AND q,y,m2; ADD t,rho,4; CSZ rho,q,t;
 AND q,y,m1; ADD t,rho,2; CSZ rho,q,t;
 AND q,y,m0; ADD t,rho,1; CSZ rho,q,t;
最初の3行はマスク用の定数である. GREGは大域レジスタを設定する命令でm0というラベルの場所を十六進(#で示す)の5555555555555555にするというような命令が並ぶ.

1行目の中程「 ;m1 」のセミコロンの次にすぐm1が続いているのは, m1をラベルにしたいからである. 4行目から命令の前に空白があるのは, ラベルがないことを示す. 行頭やセミコロンの直後に文字があるとラベルに扱われる.

さて命令の部分に入って, NEGU y,x; はxのunsigned negativeをとってyに入れる. AND y,x,y;はx & -xを作ってそれをyとした. yには最も右の1のビットだけが立っている.

AND q,y,m5; はyのビットが64ビットの左半分にあればqが0になり, 右半分にあれば, yのままの場所に1が残る.

ZSZ rho,q,32; はzero or set if zeroで, qが0ならrhoに32を入れ, そうでないなら0を入れる. 従ってyのビットが左半分にあればrhoは32になる.


次の行はm4のマスクで調べ, 以前のrhoに一応16を足したものをtに用意した上で, 上の図の16の範囲にあれば前のrhoにtを入れる. CSZは条件付きでqが0ならtをrhoに入れるが, そうでないなら何もしない.

そういう次第で, 上の図の上の2段目の区画のどこにyのビットがあるかでrhoに0, 16, 32, 48 のいずれかが入る.

同様にして8ビットごとの区間のどこにあるかでrhoに8を足す足さないとしてyのある区画を決める.

後も同様で, 変数rhoには答えの位置の番号が入るのである.

TAOCPのこのアルゴリズムの後にはまだ何か書いてある. 「最後の3行を
 SRU y,y,rho; LDB t,rhotab,y; ADD rho,rho,t;
としてもよい. rhotabは129バイトの表(その8個だけを使う)の原点である.」 これは何か.

最後の3行に来たときにはyのビットは上の図の最下段のどれかの8ビットの区画にあり, その右端のビット位置がrhoに入っている. 従ってyを右にrhoビットシフトすると, yのビットは右端の8ビットの区画に移動する. つまり

のどれかになっている. この一番上なら(1なら)0, その次なら(2なら)1, 4なら2, ..., 128なら7をrhoに足せばよい. 従って次のような表を用意すればよいわけである. たしかに129の場所があり, そのうち8カ所だけを使っている. 

MIT Schemeのbit-stringを使って実装してみた.
(define (bmake x) (signed-integer->bit-string 64 x))
(define (band x y) (bit-string-and x y))
(define (bzero? x) (bit-string-zero? x))

(define m0 (bmake #x5555555555555555))
(define m1 (bmake #x3333333333333333))
(define m2 (bmake #x0f0f0f0f0f0f0f0f))
(define m3 (bmake #x00ff00ff00ff00ff))
(define m4 (bmake #x0000ffff0000ffff))
(define m5 (bmake #x00000000ffffffff))

(define (rho x)
 (let ((bx (bmake x)) (by (bmake (- x))) (bq 0) (rho 0) (t 0))
  (set! by (band bx by))
  (set! bq (band by m5))
  (set! rho (if (bzero? bq) 32 0))
  (set! bq (band by m4))
  (set! t (+ rho 16))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m3))
  (set! t (+ rho 8))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m2))
  (set! t (+ rho 4))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m1))
  (set! t (+ rho 2))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m0))
  (set! t (+ rho 1))
  (if (bzero? bq) (set! rho t))
  rho))

(rho 1000)=>3
(rho 10000)=>4
(rho (+ (expt 2 62) (expt 2 32)))=>32
≥263ではbit-stringを作るところでエラーになったが, 他はうまくいっている.

2014年3月13日木曜日

計算機による音楽演奏

3月11日から13日まで, 東京電機大学東京千住キャンパスで情報処理学会第76回全国大会が開催された. 私は「〜コンピュータパイオニアが語る〜「私の詩と真実」」のセッションで昔話をさせられた.

その中で思い出深いパラメトロン計算機PC-1による音楽演奏にも触れた. この事は約半世紀前には一部の計算機屋によく知られていたが, 最近は知る人もほとんどいないであろうから, ブログに記録していおこうと思う.

話は56年前, つまり1958年にさかのぼる. 私は秋の半ばからニューヨークへ出張した. ニューヨーク滞在中にあちこちに出掛けたが, 12月だったかボストンへ行きMITを訪問した. ちょうどMITにいた高橋研究室の先輩がMITの名所を案内してくれた. 記憶に鮮明に残っているもののひとつは最先端のトランジスタ計算機TX-0であった. いろいろなデモを見せて貰ったなかに, コンソールのタイプライタのキーを押すとキー毎にいろいろな音が出て音楽が演奏できるのがあった. (TX-0は私が1973年から74年にMITでファカルティメンバーとして在籍した頃も学内のどこかにあったらしく, TX-0で出力したというレポートをもってきた学生がいた.)

それはとても楽しそうだったから, 東京に戻ったらなんとかしてPC-1でもやってみたいと思いつつ帰国した.

という次第で以下に述べる音楽演奏の実験をしたのは1959年の2月か3月頃であったろうか.

PC-1は理学部の研究者の科学計算にも使われていたが, もともとは高橋研の実験機だったから, 何かに使えるだろうということでフリップフロップが2個組み込まれていた. 30と31という名前である.

PC-1の命令はアルファベット1文字か, それにバリエーションを示すl(エル)がついているかである. 通常の演算には使わないアルファベットyがあったので, y30とy31がそれぞれ30と31のフリップフロップをセットする命令, yl30とyl31がリセットする命令である.

私はこのフリップフロップに目をつけた. PC-1に研究室に転がっていたマグネチックスピーカを固定し, フリップフロップのセット, リセットでスピーカのコーンを押したり引いたりしたら 音が出せるのではないか.

その話を同僚の相馬君に伝えると早速回路を作ってくれた. それが下の図である.



この辺でパラメトロン回路の見方を説明しなければならない. 図の左下がパラメトロン回路で丸がパラメトロンを表す. パラメトロンは左側が入力, 右側が出力だ. 丸の中のプラスやマイナスも入力だが, これらは定数で, プラスは常に1, マイナスは常に0である.

もうひとつ大事なことはパラメトロンは多数決素子であって, 入力は定数も含めて奇数になっている. 真を1, 偽を0で表すと, 定数がマイナスのパラメトロンはAndであり, 定数がプラスのはOrである.

入力の直前の線と直交する短かい線分はNotを示す.

またパラメトロンは3拍励振といって, それぞれのパラメトロンはI, II, IIIの3つの相のどれかであり, I相の出力がII相の入力, II相の出力がIII相の入力, III相の出力がI相の入力になる. 従って上の図でff31の下の3個のパラメトロンは, 左のプラスのあるパラメトロンをI相とすると, その右のマイナスのパラメトロンがII相, そのさらに右の白丸がIII相になる. また入出力の関係から, 左端のyやlのパラメトロンはII相になる.

他の命令を解読実行している時は左端のyのパラメトロンは0の状態で, その右隣りのAndの出力も0である. この0がフリップフロップの最左の素子に入るが, 定数の1と打ち消してIII相からの情報がそのまま出力になる, つまりフリップフロップの状態は保たれる.

フリップフロップのII相のマイナスのある素子は, 左下からの入力が0だとNotによって1になり, 定数と打ち消される.

yが1, lが0だと上のAndの出力は1になり, フリップフロップは1にセットされる.

yもlも1の時は下のAndが1になり, フリップフロップのII相の素子は2つの入力が0になるから, フリップフロップはリセットされる.

パラメトロン素子の0と1は, 実際には発振の位相が基準のパラメトロンの位相と合っているとき1, 反対位相のとき0という.

そこでIII相の素子の下のプラスだけのパラメトロンは, 1の位相で発振しており, フリップフロップも1の位相のとき, その右のトランスから出力波が得られる. 0だと上下2つの波が消し合い, 出力はない. こういう出力を整流増幅してスピーカに入れている.

プログラムとしては, 一定の時間ごとにy 31とyl 31の命令を繰り返えせばよいから簡単である.

それが下の図だ. 0番地と3番地にy 31とyl 31がある. それらの命令の次に時間を調整するための左シフト命令l nがある. この命令の実行時間はシフトする量nによって変るのでこういう時に便利である. 各命令の右端の{,}内は実行時間で, 基本的な命令は4τである. 但しτはクロック時間, 各相が出力を出す間隔であり, 100マイクロ秒くらいである.



プログラムの5番地には0番地に戻るためにo rという命令が置いてあり, これと時間を合せるために2番地に3番地へジャンプするjl 3rという命令を置く.

ここで今度は命令oを説明しなければならない. oはoutputのoで, 本来はアキュムレータの最上位6ビットをテレタイプに送る命令である. 情報の元はアキュムレータだからこの命令にはアドレス部はいらいない筈である. しかしテレタイプは計算に比べて遲いから, 次のビットを送ろうとしても, テレタイプはまだ前の文字の処理中である公算は大きい. そこでo命令は出力しようとしたとき, テレタイプが準備出来ていなければ, その番地部の示す命令にジャンプするという方式設計になっていた. (入出力と計算を同時に出来るということから, 十人の訴えを一度に聞いたといわれる「聖徳太子の機能」といった.)

仮にテレタイプが1秒に10文字印字するなら, o命令は一旦ビットを送信した後は100ミリ秒は番地の示す先にジャンプする. したがって先のプログラムでは100ミリ秒のループが作れる.

100ミリ秒経つとダミーの情報をテレタイプに送り, ループから下へ抜ける. そこでは次の音高に従い, nの値を設定してまたループに戻ることになる. (ビットが送られると テレタイプはがちゃがちゃいってうるさいから電源を切っておく.)

最後に決めるのはnの値である.



この表の左frequencyの欄はdoに対する各音の周波数を示す. しかし我々は上の図のように周期で考える方が分り易いので, 右の欄を使う. いろいろ試行したが, 下のdoの周期を120にし, reをその9/10倍, miをreの8/9倍, ...のようにすると, reのところだけ胡麻化しだが, なんとか収まるようである. それに対応するnの値を一番右の欄に示した.

doを例にすると, 1周期は120τ, 半分の波の時間は60τ, 命令の固定の時間は4τ + 6τ + 4τ=14τだから, 60から14を引いてn=46になる. 先程のプログラムのシフト命令にこれらのnの値を順々にいれて走らせるとスケールの音が聞こえるのである.

計算機で音を出すことには成功したが, 四分音符しか演奏できないから, 実用にはならなかった. でも手元に計算機があると面白いことが出来るという見本である.

2014年1月14日火曜日

開平法

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,...と引くのは電動計算機の方法というべきであったかも しれない.

2013年12月11日水曜日

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が返るのである.

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

2013年11月8日金曜日

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, 通常の十六進ではe716=23110から順に格納される.

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はどうせそうなるからというので, 読み込む先はアキュムレータであった.

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

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