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

TAOCPの演習問題7.1.3-20にGosperのハックというのがある.

正の整数xについて

u ← x&-x;
v ← x + u;
y ← v + (((x⊕v)) »2);
としてyを返す.

とにかくやってみよう.

(define (gosper x)
 (let* ((u (band x (- (expt 2 16) x)))
        (v (+ x u))
        (y (+ v (>> (/ (bxor v x) u) 2))))
 y))

(gosper 1) => 2
(gosper 2) => 4
(gosper 4) => 8
(gosper 8) => 16
(gosper 3) => 5
(gosper 5) => 6
(gosper 6) => 9
(gosper 7) => 11
(gosper 11) => 13

という次第で, 1のビットの数が同じの次に大きい整数が得られる.

その理由を考えてみた.

x=α01^a0^b

とする. つまりxは左に適当な0と1の列αがあり, その右に0が1個, その右に1がa >0, その右に0がb ≥0個あるとする.

uの式は有名な最も右の1を取り出すものだから

u=10^b

である.

v=x + uは1の並びの右端に1を足すから

v=α10^a0^b.

x ⊕ v は左のαがキャンセルされるから

x ⊕ v=11^a0^b

になり, これをuで割るから右の0の列が消えて

(x ⊕ v)/u=11^a

つまり1が右端にa +1個ならぶ. これを2ビット右シフトするから右端に1がa -1個ならぶ. これをvに足すから

y = α10^b01^a-1     (a >0に注意)

で次のyが得られたのであった. たしかにハックだね. これはMITのHAKMEMの175番にある.

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

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

TACPにはもっと凄いλ関数の計算法がある. 今回はそれを説明しよう.

xのサイズは知らなければならないが, それがどんなに大きくてもbignumのような2adic integer(2個進数)での処理だ.

まず準備としてあるアルゴリズムを考える. 具体的に4ビットとしよう. 4ビットだから対象は0≤x<16 である. c≤8と最も左のビットだけが1のh=8が登場する.

h|x-cは引かれる方が8≤ <16だから, 0≤x<16について, 差は8-c≤ <16-cの2回繰り返しになる.

たとえばc=3とすれば, 5,6,...,12,5,6,...,12になる. つまりそれぞれの最初のc個だけが<8になる. これにxをビットごとorすれば, 右半分はxの最も左のビットが1なので≥8になる. 従って全体では左端のc個だけが<8, それ以外は≥8になって, hでビットごとandをとると左端のc個が0, それ以外は8になる. c=3の例では0,0,0,8,8,...,8だ. このアルゴリズムで重要なのは, 4ビットの範囲の左から繰下げがないことである.

以下のアルゴリズムにはこれを多用する.

xが如何に大きくても, gを2の羃, nをg²とし, x<2ⁿの⌊lg (x)⌋を計算する.



以下の説明ではg=4, n=16とする. そうすると最初の行にあるlは2¹⁶-1, つまり1111111111111111を2⁴-1, つまり1111で割るから0001000100010001になる. 従ってhはそれを3ビット左へシフトした1000100010001000になる. 最初の行は16ビットを4ビットの区間4つのそれぞれについて, 4ビットがオール0かそうでないかを上のアルゴリズムで調べる. c=1の場合だ. 各区間はその区間がオール0か否かにより, 0か8になる. その結果がt₁だ. だからt₁は各々の区間の結果をp,q,r,sのビットで表すと, p000q000r000s000になる. p,q,r,sの少なくても1つは1である. 次の行はまずaを計算する. 2¹²-1つまり111111111111を2³-1つまり111で割るからa=001001001001である. これを先ほどのt₁に掛けるから, 積は

         p000q000r000s000
      p000q000r000s000
   p000q000r000s000
p000q000r000s000

の和,
p00pq0pqrpqrsqrs0rs00s000
222221111111111
4321098765432109876543210←ビット位置

になり, この下から12,13,14,15ビットをとるからyはpqrsにlを掛けたもの, pqrspqrspqrspqrsになる. 言い換えれば区間の情報を4ビットに詰め込んで4(=g)回コピーしたものが出来た.

今度はt₂の計算だ. bは2²⁰-1を2⁵-1で割るからb=1000010000100001となる. だから式中のbは左端の区間から順に1000, 0100, 0010, 0001で, これらをcとしてまたあのアルゴリズムを使う.

pqrs≥1000>100>10>1ならすべての区間で8になり,
1000>pqrs≥100>10>1なら0888になり,
1000>100>pqrs≥10>1なら0088になり,
1000>100>10>pqrs≥1なら0008になる.
従ってt₂は二進法では
1000100010001000
0000100010001000
0000000010001000
0000000000001000
のいずれかになり, このそれぞれからmは
1111111111111111
0000111111111111
0000000011111111
0000000000001111
になり, 4ビット右シフトしたものとxorするとmとして
1111000000000000
0000111100000000
0000000011110000
0000000000001111
が出来る.

次の行は最も左の1を持つ区間をmでマスクしてとりだし, lを掛けて左へ4個ならべ16ビットに区切ってから12ビット右シフトして右端の4ビットに置き, さらにlを掛けて4区間全体に展開する. 先ほどのpqrsから最も左の区間を選んだように, 今度は最も左のビットの位置がt₃に出来る.

整理すると最も左の1のある区間が右端にあればt₂は0008, その1つ左にあれば0088, 左端の区間なら8888である.

その区間で最も左の1が右端にあればt₃は0008, 左端にあれば8888である.

これからビットの位置λを計算するにはこうする.

t₂を5ビット右シフトするから, 0,4,44,444のいずれかになる.
t₃を7ビット右シフトするから, 0,1,11,111のいずれかになる.

これらを上下から取って足し, lを掛けて12,13,14,15ビット目をとるとそれがλである. たとえば444と1だとすると和は445になり,

   445
  445
 445
445

この4+4+5で13が得られる. 下からの繰り上げはないか心配だが, この計算は十六進であり, 最大の444と111でも5+5+5になって繰り上げはない.

このアルゴリズムは繰り返しがないのでO(1)と書いてあるが, 乗算が5回もあるから実用にはならないとTAOCPはいう.

次はDoループはあるがステップ間のジャンプはないアルゴリズムである.



d=4とするとn=d・2^d=64になる. この値を使ってアルゴリズムを調べてみよう.

B1. λ=0とした後, ⌈lg n⌉=6だから, k=5,4についてxを32ビット, 16ビットの幅に1があるかどうかでxを右シフトし, λを増やす. このステップが終わるとxは最も左のビットを含む16ビットの区間になり,　その区間がどこにあるかに従ってλは0, 16, 32, 48のいずれかになっている.

B2. 16ビットの区間に縮められたxを4つの区間にコピーする.

次にμ_d,_kというのがでてくるが, 右から2^kごとに1と0を繰り返す長さ2^dのビット列である.

μ_d,k=(2^2d-1)/(2^2k+1)

今はd=4だから
μ_4,3=0000000011111111
μ_4,2=0000111100001111
μ_4,1=0011001100110011
μ_4,0=0101010101010101

μ_4,3の0の部分を左半分, 1の部分を右半分といおう. 同様に
μ_4,2の0の部分を左四半分, 1の部分を右四半分といい,
μ_4,1の0の部分を左八半分, 1の部分を右八半分といい,
μ_4,0の0の部分を左十六半分, 1の部分を右十六半分といおう.

すると下の図で灰色の部分が左なんとか半分, 白い部分が右なんとか半分である.


B3. これらのμを繋げてビットごとにnotをとったものとxのandをyとする. つまりxの左なんとか半分の1をとりだす.

B4. x⊕yは右なんとか半分の1になるわけだ. そしてh=8000800080008000としておなじみのアルゴリズムを使う.

そうすると左部分をアルゴリズムのx, 右部分をcとすることになるので, x<cの時に0, x≥cの時に8だったように, 左≥右なら0, そうでなければ8000になる.

上の図で最上段の0〜fはビット番号で, そのビット6に最も左の1があるとすると上から順に
左<右
左>右
左>右
左<右
だから, B4のtはそれぞれの区画で0, 8000, 8000, 0になる.

B5. ⌈lgd⌉=2なので, 0,1のkについて
2⁴-2⁰=15ビット,
2⁵-2¹=30ビット
ずつ左シフトして加える. 8の立つビットを左からp,q,r,sとするとそれらのビット位置はそれぞれp=63, q=47, r=31, s=15である.

まず15ビット左へ移るから, qは62, sは30になりpq, rsは並ぶ. それをさらに30ビット左シフトすると31,30にあったrsは61,60になり, pqrsは並ぶ.

2⁶⁴の剰余をとるから64ビットレジスタの左端に出来る.

先ほどの例ではその4ビットに0110, つまり最も左のビット位置の6ができる. それをn-d=60ビット右シフトし, すでに得られていたλと足すと最も左の1のビット位置が得られるのである.

ところで理解して頂けたであろうか. 私はまずSchemeでプログラムを書いて実行したりして, 様子を見ながらアルゴリズムがどうなっているかを考えた. 分かってしまえば簡単なのだが, 最初に述べたアルゴリズムが何をしているか分からないうちは, なんとも不思議に思うアルゴリズムであった.

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

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

二進法で表した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なら
±2^E-1023(1+F/2⁵²)
を表す.

だから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, ..., 2^k〜2^k+1-1はk, 128〜255は7になっている.

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

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

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

TAOCPには

λx = λy   if and only if    x ⊕ y ≤ x & y

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

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

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

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

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

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,1	8 → 0,1
1 → 2,3	9 → 2,3
2 → 4,5	10 → 4,5
3 → 6,7	11 → 6,7
4 → 8,9	12 → 8,9
5 → 10,11	13 → 10,11
6 → 12,13	14 → 12,13
7 → 14,15	15 → 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ビットに相当する場所にだけ二進数が書いてある. 下の方, 2¹³, 2¹⁴, 2¹⁵を掛ける辺りはサイクルの右端の4ビットがレジスタからはみ出しているが, サイクルの左端が0000だったのが幸いしてちょうど輪になっている.

さてこれから分るように

×1は2⁰を掛けたのだから, 表の0番には0を置く.
×2は2¹を掛けたのだから, 表の1番には1を置く.
×4は2²を掛けたのだから, 表の3番には2を置く.
×8は2³を掛けたのだから, 表の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にある表と合っている.

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

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

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となる. 二進法のビットを右からa₀, a₁,...と番号をつけると, ...0010の1はa₁だから位置としては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

≥2⁶³ではbit-stringを作るところでエラーになったが, 他はうまくいっている.

計算機による音楽演奏

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

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

開平法

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