［C++］de Bruijn sequenceを用いたLSB/MSB位置検出テクニック

de Bruijn sequence（ド・ブラウン列）

de Bruijn sequenceとは、いくつかの文字である長さの文字列を作ることを考えた時、その組み合わせの全てを重複なく含んだ文字列のことを言います。

例えば、文字a, bを使った長さ3の文字列はaaa, aab, aba, baa, abb, bab, bbb, bbaの8通りなので、このde Bruijn sequenceは例えばaaababbbになります（通常いくつか存在します）。
aaababbbを先頭から3文字、1文字づつ右にずらしながら見ていくと（途中で巡回します）、確かに8通りの組み合わせ全てを含んでおり、重複が無い事が分かります。

aaababbb
aaa|||||
 aab||||
  aba|||
   bab||
    abb|
     bbb
      bba (先頭へ戻る
       baa

文字a, bを0, 1に置き換えてやれば、2進数列のde Bruijn sequenceを考える事ができそうです。

de Bruijn sequenceによるLSB位置検出

この2進数列のde Bruijn sequenceを用いて、高速にLSB位置を検出するアルゴリズムがあります。

Wikipediaより

unsigned int v;   
int r;

static const int MultiplyDeBruijnBitPosition[32] = 
{
  0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, 
  31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
};
r = MultiplyDeBruijnBitPosition[((uint32_t)((v & -v) * 0x077CB531U)) >> 27];

これによって符号なし32bit整数値の最下位ビット位置を得ることができます。
何やってるのかさっぱりわかりません・・・

なにがおきているの？

少しコードを整理してみます（ついでにC++的になおします）。

int lsb_pos(std::uint32_t v) {

  //1
  static constexpr int MultiplyDeBruijnBitPosition[32] = 
  {
    0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, 
    31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
  };

  //2
  std::uint32_t pop_lsb = (v & -v);

  //3
  std::uint32_t hash = std::uint32_t(pop_lsb * 0x077CB531U) >> 27;

  //4
  int r = MultiplyDeBruijnBitPosition[hash];

  return r;
}

少し見やすくなったでしょうか、番号を振ってあるところを順に見て行きますと

1. 謎の配列、ここに入ってる値がビット位置になっている様子
2. 謎のビット演算テク、これは最下位ビットだけを残すポピュラーなテクニック
3. 一番意味分からない所、後続の処理を見るにどうやら[0, 31]の数値を出力している
4. 3の結果をインデックスとして1の配列を参照。参照先がビット位置に対応している

3番以外は何やってるのかわかるのではないかと思います。3番が分からないので全部わからないのですが・・・

しかし、整理したコードをよく見るとこのアルゴリズムはハッシュテーブルによって高速にビット位置を求めていることが見えてきます。
すると、意味分からない3番目の処理はハッシュ値を求めている事に相当しそうです。そして、その入力は最下位ビットのみが立っている状態になっています。

つまり3番目の処理は、最下位ビットだけが立った値に対して完全ハッシュ（ダブりが無いハッシュ）を求めていることになります。

de Bruijn sequenceによる完全ハッシュ

//3
std::uint32_t hash = std::uint32_t(pop_lsb * 0x077CB531U) >> 27;

3番目の処理に出てくる謎の数字0x077CB531ですが、これが実はde Bruijn sequenceになっています。2進数に直して5文字づつ見て行くと確かに[0, 31]の数字（2進数列）がすべて含まれ、なおかつ重複がないことが確認できるでしょう（右から4ビットの所は先頭へ循環して見る必要があります）。

5文字・・・2^5 = 32であることに気付くと、27ビット右シフトというのは最上位5ビット分を残す処理である事が分かります（32 - 5 = 27）。

すると残った所はde Bruijn sequenceとpop_lsbのかけ算です。普通の数値のかけ算なら結果がどうなるかを考えるのは少し難しいですが、pop_lsbはどこか1ビットだけが立った値です。
つまり、その値は必ず2^nの値になります。その2のべき乗数値との掛け算はすなわちnビット左シフトに相当します。

ここでの2進de Bruijn sequenceは左から5桁づつ重複なく5ビットの表現全てを含んでいます。
整数型のビット数（今は32）未満であれば、nビット左シフトして左から5桁分を数値として読み取ると、n毎に異なった値が得られます。
しかも、5ビット数値なので2^5 = 32未満 = [0, 31]の値が得られます。

今、入力は最下位ビットのみが立った値であり、32ビット符号なし整数型なら取りえる値は32個だけです。従って、そのビット位置に応じた[0, 31]の個別の値をこの一行は計算しています。
あとは、テーブルによってその値をそのビット位置を示す数字に対応させてやれば、処理は完了です。つまり、1番初めに宣言されている配列はこの対応を取っているものである事が分かります。

8ビットの例

文字だけだと分からないので例を見てみましょう、しかし32ビットは長いので8ビットで見てみます。
8ビットの値は8桁なので、3ビットでその位置を表現可能です。従って、0, 1を使った3文字を尽くすような長さ8のde Bruijn sequenceが必要になります（元論文から拾ってきます・・・）。

先ほどの3番目の処理は以下のようになります。

std::uint8_t hash = std::uint8_t(pop_lsb * 0x1DU) >> 5;

当然ですがpop_lsbはどこか1ビットだけが立った8ビット符号なし整数値であることを前提とします。
つまり、入力となるpop_lsbは8個の値しかとりえません。それぞれについて処理を見てみると以下のようになります。

pop_lsb	pop_lsb * 0x1D	hash	index
1	0001 1101	000	0
2	0011 1010	001	1
4	0111 0100	011	3
8	1110 1000	111	7
16	1101 0000	110	6
32	1010 0000	101	5
64	0100 0000	010	2
128	1000 0000	100	4

こうしてみるとどこか1ビットだけが立った値をうまいこと3ビット数値に押し込めた完全ハッシュになっている事が分かるでしょう。
後はこのindex位置に、対応する桁数を持つような配列を用意してあげるだけです。

constexpr int table[8] = { 1, 2, 7, 3, 8, 6, 5, 4 };

上の表で得られたindexをtable[index]として値を取得すると、元のpop_lsbの立っているビットの位置が得られる事が分かるでしょう。

Nビット整数への一般化

Nは2のべき乗である必要がありますが、上記アルゴリズムは以下のようにNビットへ一般化できます。

ここで、Nは符号なし整数型の幅、debruijnは0, 1を使ったlog_2(N)文字の組み合わせを尽くすような長さNの適切なde Bruijn sequenceです。
そのようなde Bruijn sequenceを求める方法はいくつかあるようです。しかし良く分からない・・・

なお、使用するde Bruijn sequenceは先頭log_2(N)桁が0で始まる必要があります。これはオーバーフロー対策と、掛け算（左シフト演算）時に全体が自然に循環するようにするためです。

64ビット数値のLSB/MSB位置を求める

では、64ビット符号なし整数型の最上位/最下位ビット位置を求める処理をC++で実装してみます。

最下位はここまでやってきたことの流れで実装できますが、最上位ビット位置は少し違った処理が必要です。
とはいっても、最上位ビット位置だけを残すビット演算テクニックが必要になるだけで、幸いそれはネットに落ちてました。

64bitで使用するde Bruijn sequenceは0x03F566ED27179461になります。これもネットに落ちてました・・・

最後の右シフト量は上記式から64 - log_2(64) = 58と求められます。

最終的に桁位置に写すテーブルは簡単な計算で求められます。

constexpr auto hash_64(std::uint64_t x) -> int {
  return std::uint64_t(x * 0x03F566ED27179461UL) >> 58;
}

inline constexpr char hash2pos[] = {1, 2, 60, 3, 61, 41, 55, 4, 62, 33, 50, 42, 56, 20, 36, 5, 63, 53, 31, 34, 51, 13, 15, 43, 57, 17, 28, 21, 37, 24, 45, 6, 64, 59, 40, 54, 32, 49, 19, 35, 52, 30, 12, 14, 16, 27, 23, 44, 58, 39, 48, 18, 29, 11, 26, 22, 38, 47, 10, 25, 46, 9, 8, 7};

constexpr auto msb_pos(std::uint64_t x) -> int {
  if (x == 0) return 0;

  //最上位ビットだけを残す
  x |= (x >> 1);
  x |= (x >> 2);
  x |= (x >> 4);
  x |= (x >> 8);
  x |= (x >> 16);
  x |= (x >> 32);
  x = x ^ (x >> 1);

  int h = hash_64(x);

  return hash2pos[h];
}

constexpr auto lsb_pos(std::uint64_t x) -> int {
  if (x == 0) return 0;

  std::uint64_t v = x & -x; //最下位ビットだけを残す

  int h = hash_64(v);

  return hash2pos[h];
}

[Wandbox]三へ( へ՞ਊ ՞)へ ﾊｯﾊｯ

なお、このテクニックを応用することで更なるビット演算黒魔術を行えるようです（1が2つ並んでいる最上位の位置を求めるとか・・・）。
また、de Bruijn sequenceは南京錠の様な対象への総当たり攻撃の効率化や、ハミルトン路を求める問題をオイラー路を求める問題へ変換するなど色々応用できるみたいです（良く分かってない）。

参考文献

• Using de Bruijn Sequences to Index a 1 in a Computer Word
• 一番右端の立っているビット位置を求める「ものすごい」コード - 当面C#と.NETな記録
• De Bruijn sequence - Wikipedia
• De Bruijn列 - Thoth Children!!
• 組合せとグラフの理論 ( 塩田 )
• de Bruijn Graph を使った de novo アセンブリの発想がすごい件 - ほくそ笑む

この記事のMarkdownソース