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の値を順々にいれて走らせるとスケールの音が聞こえるのである.

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