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 |
このグラフですべてのノードを経由して元へ戻る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にある表と合っている.
0 件のコメント:
コメントを投稿