💡

de Bruijn列でBSR32をO(1)で実装できる理由

2021/04/07に公開

背景

mimalloc が公開された
その中で BSR32(bit scan reverse) がde Bruijn sequence を利用してO(1)で実装されていると話題に
なんでde Bruijn sequence でBSRが実装できるのか知りたくなって調べたがなかなかわかんなかった
唐突にわかった(←ｲﾏｺｺ)

TL;DR

これ考えたやつはすげぇ

用語

mimalloc - Microsoftにより公開された高パフォーマンスなメモリアロケーター。jemallocとかtcmallocとかの競合
bit scan reverse - 「1になっている最も高位のビットのインデックスを0-indexedで計算」(引用元) 0b0001 なら 0、 0b0010 なら 1、 0b0100 なら 2 という感じ。 0b0000 については未定義となる。BSR32なら32ビットなので、結果は0～31になる。
de Bruijn sequence - 詳しくはWikipedia参照。簡潔に言うと「k 種類のアルファベットの n 文字の全組み合わせを含む最小の文字列のこと」となる。

コード

static inline uint8_t mi_bsr32(uint32_t x) {
  // de Bruijn multiplication, see <http://supertech.csail.mit.edu/papers/debruijn.pdf>
  static const uint8_t debruijn[32] = {
     31,  0, 22,  1, 28, 23, 18,  2, 29, 26, 24, 10, 19,  7,  3, 12,
     30, 21, 27, 17, 25,  9,  6, 11, 20, 16,  8,  5, 15,  4, 14, 13,
  };
  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;
  x++;
  return debruijn[(x*0x076be629) >> 27];
}

これだけ見てなんのことだかわかった人はめちゃくちゃすごい。そしてこれを思いついた人はアホほど天才だろう。

解説

このコードは前半部分と後半部分にわかれている。また前半部分は2つにわかれているので、合計3つに分かれていることになる。1つずつ見ていこう。

前半

前半では1である最上位ビットの1つ上(左)を1にし、残りをすべて0にしている。詳しく見ていく。なお最初の定数配列は忘れて良い。

前半Aパート

まずは以下の部分。この部分では1である最上位のビット以下全部を1に立てている。

  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;

つまり3ビットで考えるならば以下のようになる。

before	after A
0b0000	0b0000
0b0001	0b0001
0b0010	0b0011
0b0011	0b0011
0b0100	0b0111
0b0101	0b0111
0b0110	0b0111
0b0111	0b0111

前半Bパート

Bパートははシンプルに1行だけだが、Aパートによりある箇所より下のビットがすべて立っているため、繰り上がって1つ左だけが1となり残りはすべて0となる。

x++;

以下に前半パートのA,B合わせた結果を示す。

before	after A	after B
0b0000	0b0000	0b0001
0b0001	0b0001	0b0010
0b0010	0b0011	0b0100
0b0011	0b0011	0b0100
0b0100	0b0111	0b1000
0b0101	0b0111	0b1000
0b0110	0b0111	0b1000
0b0111	0b0111	0b1000

前半の出力は 2, 4, 8, 16... という 2^m (m > 0) になることを抑えておけば良い。before が 0b0000 のケースは未定義だから考慮する必要はない。ちなみに m は BSR+1 になっている。

後半

後半のコードのキモは (x*0x076be629) >> 27 に集約される。この部分は2つに分けて考えると良い。

外側の (...) >> 27 は左辺の先頭5ビットだけを取り出す意図がある。

内側の (x * 0x076be629) はなんのこっちゃだが、このとき x が 2^m (m > 0) になっていることを思い出して欲しい。つまりこれはビットシフトを意味してて、 0x076be629 << m ということになる。

つまり後半は定数 0x76be629 を左ビットシフトして先頭5ビットを取り出しているわけだ。となると 0x76be629 とはなんなのかが問題だ。

`0x76be629` とは

結論を言うと 0x076be629 が k=2, n=5 の de Bruijn sequence だった。つまりアルファベットは 0 と 1 の2種類、長さが5の単語すべてを含んでる。もうちょと噛み砕いて言うと 00000, 00001, 00010, ..., 11101, 11110, 11111 の32種類の単語=2進数で考えれば数字、を含むのが 0x076be629 である。

ためしに2進数で展開してみよう。

00000111011010111110011000101001

ここから5ビットずつ順番に取り出してみると次のようになる。見事に0~31の全単射を表現できている。28以降は右側に0が補われることに注意。

`m`	bits	output
0	`00000`	0
1	`00001`	1
2	`00011`	3
3	`00111`	7
4	`01110`	14
5	`11101`	29
6	`11011`	27
7	`10110`	22
8	`01101`	13
9	`11010`	26
10	`10101`	21
11	`01011`	11
12	`10111`	23
13	`01111`	15
14	`11111`	31
15	`11110`	30
16	`11100`	28
17	`11001`	25
18	`10011`	19
19	`00110`	6
20	`01100`	12
21	`11000`	24
22	`10001`	17
23	`00010`	2
24	`00101`	5
25	`01010`	10
26	`10100`	20
27	`01001`	9
28	`10010`	18
29	`00100`	4
30	`01000`	8
31	`10000`	16

これにより0~31の数値が、0~31の別の数値に全単射できている。あとはこの数値をさらに配列を用いた全単射で32種類の任意の数値=実際のBSRの値に変換している。

まとめ

mimallocのBSRは、入力を 2^(m+1) へ算術変換してから、de Bruijn sequenceにビットシフトとマスク演算を利用して m をBSRとして算出していた。