横方向の和の和

昨日のブログではさぼったnusum (2^r-1) から nusum 2^rを計算してみた.

前回のような表を書くと

n	nusum	νn	es
2r-1	nusum 2r-1	r	(r-1 r-2 r-3 ... 0)
2r	nusum 2r	1	(r)

羃の列が(r-1 r-2 r-3 ... 0)のnusum (2^r-1) の式は

nusum (2^r-1) = (r-1)2^r-2 + (r-2+2)2^r-3 + (r-3+4)2^r-4 + ... + (0+2r-2)2^-1

で, これにrを足したものが

nusum 2^r = r.2^r-1

になればよい.

nusum (2^r-1)
={各項の係数の中を計算する}
(r-1)2^r-2+r2^r-3+(r+1)2^r-4+...+(r+(r-2))2^-1
={降順の羃の列で同じ係数を持つものをならべる}
(r-1)2^r-2+(r-1)2^r-3+(r-1)2^r-4+...+(r-1)2^-1  {1行目}
             +2^r-3+2^r-4+...+2^-1                    {2行目}
                    +2^r-4+...+2^-1                    {3行目}
                          +2^r-5+...+2^-1              {4行目}
                          +...
                          +2^-1                          {r行目}
={各行を2^p+2^p-1+... +2^q=2^p+1-2^q で置き換える}
(r-1)(2^r-1-2^-1)
+2^r-2-2^-1
+2^r-3-2^-1
+...
+2⁰-2^-1
={1行目を展開し2行目以下の第1項を足し合せる 2行目以下のr-1個の2^-1を足す}
(r-1)2^r-1-(r-1)2^-1+(2^r-1-2⁰) -(r-1)2^-1
={いろいろ消えて}
r2^r-1-2^r-1-(r-1)+2^r-1-1
=r2^r-1-r.

rを足せばnusum 2^rになる.

存外簡単であった.

横方向の和の和

TAOCPの演習問題713-42にこういうのがある.

「問題42. [M21] e₁ > ... > e_r ≥ 0について, n = 2^e₁ + ... + 2^e_r なら, 和 ∑_k=0^n-1νkを, 冪e₁, ..., e_rを使って表せ.」

ここでνkはkを二進法で表わした時の1の数であり, 横方向の和といったりする.

ν0=0; ν1=1; ν2=1; ν3=2; ν4=1; ν5=2; ν6=2; ν7=3; ν8=1; ν9=2; ...

従って問題はnが与えれた時, n-1までの各整数にあった1の総和nusumを計算するものである.

n=10とすれば, 上のνの値を使いnusum 10 = 0+1+1+2+1+2+2+3+1+2=15

しかし途中のνを全部計算して求めるのはしんどいから, 10を二進法で1010と表わし, 10 = 8 + 2 = 2³ + 2¹. そのe₁=3, e₂=1からnusumを得たいという問題である.

私はこの解答にある J.-P. Allouche, J. Shallitの本 Automatic Sequences (2003)にあるというフラクタルの図が見たかったが, その本が簡単に見られそうもなかったので, PostScriptで描いてみることにした. そこで解答は次のようだ.

「解答. rについての帰納法により,

e₁2^e₁-1+ (e₂+2)2^e₂-1 + ... + (e_r+2r-2)2^e_r-1.

D. E. Knuth, Proc. IFIP Congress (1971), 1, 19--27. この和のフラクタル模様はAlloucheとShallitの本の図3.1と3.2に示してある.] またS_n'(1)も考察せよ. ただし S_n(z) = ∑_k=0^n-1z^νk = (1+z)^e₁+ z(1+z)^e₂+... +z^r-1(1+z)^e_r.」

つまりe₁=3, e₂=1だったから 3*2^3-1 + (1+2)*2^1-1 = 3*2² + 3*2⁰ = 12 + 3 = 15なわけだ.

もう少し計算してみる. 下の表は左端がn. 次がnusum. 各行のnusumは 1行上のnusumに1行上のνnを足したものになっている. その次のかっこ 内はe₁...e_r. さらにその右は上の 式による計算とその和を示す.

n	nusum	νn	es
1	0	1	(0)	0.2-1	0
2	1	1	(1)	1.20	1
3	2	2	(1 0)	1.20 + 2.2-1	2
4	4	1	(2)	2.21	4
5	5	2	(2 0)	2.21 + 2.2-1	5
6	7	2	(2 1)	2.21 + 3.20	7
7	9	3	(2 1 0)	2.21 + 3.20 + 4.2-1	9
8	12	1	(3)	3.22	12
9	13	2	(3 0)	3.22 + 2.2-1	13
10	15		(3 1)	3.22 + 3.20	15

さてこれでrによる帰納法が出来るであろうか.

偶数のnからそれより1大きいn'に進むとき, eの列は 最後にr+1項目として0が追加される. だから

nusum n+1 = e₁2^e₁-1+ (e₂+2)2^e₂-1 + ... + (e_r+2r-2)2^e_r-1 +(e_r+1+2(r+1)-2)2^{e_{r +1}-1} .

e_r+1=0とすると

nusum n+1 = e₁2^e₁-1+ (e₂+2)2^e₂-1 + ... + (e_r+2r-2)2^e_r-1 +2r.2^-1
=nusum n+r = nusum n+νn = nusum n'.

一方n=7から8へのように繰上がりがあるときはどうするか.

31から32の計算をしてみると, eは(4 3 2 1 0) から(5)になる.

このときのnusumは

4.2³+5.2²+6.2¹+7.2⁰+8.2^-1

これを次のように計算してみる.

3(23+22+21+20+2-1)
+(23+22+21+20+2-1)
+(22+21+20+2-1)
+(21+20+2-1)
+(20+2-1)
+(2-1)
=
3(24-2-1)
+(24-2-1)
+(23-2-1)
+(22-2-1)
+(21-2-1)
+(20-2-1)
=
3(24-2-1)
+(25-20)
-5(2-1)
=
3*15.5+31-2.5
=75

nusum 32は(5)だから, 5.2⁴ = 80. さっきの75にν31=5を足して 80になる.

一般的には上の(4 3 2 1 0)が(r-1 r-2 ... 0)なのでr について同様の計算をしなければならない. しかしもう出来そうになったから フラクタル図の方へ急ごう.

nusumの式をPostScriptで書くことになる.

/es {4 dict begin /s exch def /m exch def /n exch def
/l s length def
n 0 eq{s}{n 2 mod 0 eq{n 2 idiv m 1 add s es}
{n 2 idiv m 1 add m s aload pop l 1 add array astore es}
ifelse}ifelse end} def

/nusum {4 dict begin /n exch def
/ee n 0 [] es def /sum 0 def
0 1 ee length 1 sub{/i exch def /e ee i get def
/sum e i 2 mul add 2 e 1 sub exp mul sum add def} for
sum end} def

0.33 setlinewidth
0 1 560{/x exch def
x 0 moveto 0 x nusum 4 div rlineto stroke} for

最初の手続きesはnの冪の列を返す. n 0 [] esのように呼ぶ. 10 0 [] es => [3 1]

手続きnusumが総和の計算である. 32 nusum => 80.0, 560 nusum => 2480.0

n=0から560までを描いたのが下の図だ. フラクタルと言われればそう見えなくもないが.



横軸の目盛は100ごと, 縦軸の目盛は500ごとである.

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

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を作るところでエラーになったが, 他はうまくいっている.