簡潔ビットベクトル(完備辞書)

この記事はISer Advent Calendar 2018の15日目に書かれた記事です。昨日はどーくんさんによる奇天烈音楽入門でした。

アドカレに向けて最近Wavelet Matrixを実装したのですが、そこに使われる「簡潔ビットベクトル」というデータ構造の解説で力尽きました。それでも簡潔ビットベクトル自体も中々奥が深いデータ構造なので、こちらだけでも上げておきます。

機能と計算量

「簡潔ビットベクトル」とは、長さ $N$ のbit列に対して以下の操作を実現するデータ構造である。

名称	内容	時間計算量
( $build$ )	前計算	$O(N)$
$access$	$K$ 番目の要素を求める。	$O(1)$
$rank$	特定区間に含まれる $1$ の数を求める。	$O(1)$
$select$	$K$ 番目に現れる $1$ の位置を求める。	$O(\log N)$

$select$ を $O(1)$ に落とすことも可能らしいが、実装が重いらしいためここは一般に妥協されている。

そしてこれが重要なのだが、このデータ構造に必要な空間計算量は $N + o(N)$ bitとなる。

それ累積和でよくね？

実のところ、上の機能を実現したいだけなら「累積和¹」を使えば事足りる。これも上の機能を同様の時間計算量で実現できる上に、実装が非常に楽。

しかし簡潔ビットベクトルには累積和を上回るポイントがある。それが「空間計算量」だ。

「bit列の累積和」の空間計算量は $O(N)$ と捉えられがちだが、厳密に評価すると $N \lg N$ bit²となる。この $\lg N$ は何かと言うと、累積和の値を保持するのに必要なbit数である。

例えばbit列の長さが $255 = 2^8 - 1$ だったとしよう。このとき累積和の各値は $[0, 2^8)$ の範囲に収まる。したがってこれらは $8$ bit符号なし整数で保持すればよい。

これと同様の理屈で、長さ $N$ のbit列の累積和は、 $\lg N$ bit整数で値を保持する必要がある。よって全体で $N \lg N$ bitが必要となるわけだ³。

かたや簡潔ビットベクトルはどうか。後に解説する通り、それを考慮した上でも空間計算量は $N + o(N)$ bitとなる。

具体的に必要なbit数を概算すると以下の通り。

$N$	$2^{16}$	$2^{32}$
累積和	$2^{20}$	$2^{37}$
簡潔ビットベクトル	$2^{17} + 2^{12}$	$2^{32} + 2^{31} + 2^{29}$

所詮 $\log$ が落ちた程度の差なので劇的に小さくなるわけではないが、上のように確実に小さくはなる。

そういうわけで、「bit列の累積和を大量に用意したいけど、累積和だとメモリを食いすぎてしまう...」という場合⁴に簡潔ビットベクトルが役に立つ。

rankの仕組み

簡潔ビットベクトルの実装でメインとなるのは $rank$ である。逆にこれさえ実装できれば完成したも同然。というのも、

$access$ は元のデータを参照するだけ
$select$ は $rank$ によって二分探索すればいい

からだ。

そういうわけで、以降は $rank$ の仕組みに重点を置いて話すことにする。

まず $rank$ の仕組みは、「平方分割」という技法の応用のようなものになっている。これを念頭に置けばより理解が進みやすいと思われるので、先にそちらを軽く説明する。

平方分割

「平方分割」とは、長さ $N$ の数列を長さ $\sqrt N$ の数列 $\sqrt N$ 個に分割することで、区間に対する計算を高速化する技法である。

以下は長さ $16$ の数列の区間和を求める際に、平方分割を適用した例である。このように区間を「スキップする」感覚で、計算量を $O(\sqrt N)$ に落とすことができる。

f:id:misteer:20181215112156j:plain

chunkとblock

これと同様にビットベクトルもbit列をブロックに分割するのだが、平方分割と違って2段階に分割する。

まずは「chunk⁵」という大ブロックに分割する。そして各chunkの境界に、そこまでに現れた1の個数の合計を記録しておく。先程の平方分割では1つのブロック内だけを見ていたが、今回はそこまでの全ブロックを見ていることに注意。

そして次に各chunkを「block⁶」という小ブロックに分割する。こちらも保持するデータはchunkと同様だが、各chunk毎に独立となる。

具体例は以下の通り。

f:id:misteer:20181215112222j:plain

このデータによって、 $rank$ は以下のように実現できる。要領は平方分割と同じ。

f:id:misteer:20181215112238j:plain

chunk、block幅

まだchunkとblockの幅(以降それぞれ $(c, b)$ とおく)について触れていなかったが、これは $((\lg N)^2, (\lg N) / 2)$ がオススメ。この値にすることで、このデータ構造が簡潔性を持つことになる(簡潔性は後の項で証明する)。

例えば $N = 2^{16}$ なら、 $(c, b) = (256, 8)$ となる。最後に載せている実装例は、このケースに合わせている。

ただし注意点として、この幅を $N$ に応じて変える必要はない。これは単に実装が非常に重くなるからだ。想定している最大ケースに合わせて、固定してしまうのがよい。

popcount

上図にて黒で囲まれたbit列の処理(通称popcount)だが、1つずつbitを見ていると $O(b)$ かかってしまう。

ここの計算量を減らす方法として、以下の3つが挙げられる。

「bitパターン」と「1の数」を対応させるtableを作っておく。
- 理論的にはこれが最善で $O(1)$ 。
- ただし実装は若干重くなる。
- 加えて厳密に $O(1)$ にする必要が大抵ないため、実用上はあまり使われていない。
__builtin_popcountを使う。
- GCCのビルトイン関数の一種。
- 一番楽だが環境依存が大きい。
シフト演算を上手く使う。
- シフト演算で上手いことやると時間計算量 $O(\log w)$ でのpopcountを自力で実装することができる。
- __builtin_popcountがない環境ではこれがオススメ。
- 詳しい実装は、記事後半の参考記事を参考のこと。

計算量をわざわざ $O(\log w)$ から $O(1)$ に落とすメリットが小さいため、実用的には実装が楽で空間計算量の小さい2か3が主に使われる。

しかし今回は理論を重視し、以降1による実装で話を進めることにする。

簡潔性の証明

それでは厳密な空間計算量を計算することで、簡潔ビットベクトルの簡潔性を示す。理論的な話になるので、特に興味がなければ「実装」の項に移っていただいて構わない。

まず上で述べた通り、chunkとblockの幅 $(c, b)$ は $((\lg N)^2, (\lg N) / 2)$ とする。このとき、各変数に必要な空間計算量は以下のようになる。

元データ

これは普通に $N$ bit必要となる。

chunk

chunkの数は全部で $N / c$ 個。各値は $N$ を上回らないので、 $\lg N$ bit整数で保持できる。

よって必要な空間計算量は $\displaystyle \dfrac{N}{(\lg N)^2} \times \lg N = \dfrac{N}{\lg N}$ bitとなる。

block

blockの数は全部で $N / b$ 個。各値は $c$ を上回らないので、 $\lg c = 2 \lg (\lg N)$ bit整数で保持できる。

よって必要な空間計算量は $\displaystyle \dfrac{N}{(\lg N) / 2} \times 2 \lg (\lg N) = \dfrac{4 N \lg (\lg N)}{\lg N}$ bitとなる。

table

最後にpopcountで参照するtable。blockのbitパターンは計 $2^b = \sqrt{N}$ 通り。各bitパターンについて、1の数は $b$ を上回らないため $\lg b = \lg (\lg N) - 1$ bit整数で保持できる。

よって必要な空間計算量は $\sqrt{N} (\lg (\lg N) - 1)$ bitとなる。

以上を合計すると、簡潔ビットベクトルが必要とする空間計算量は

$\displaystyle N + \dfrac{N}{\lg N} + \dfrac{4 N \lg (\lg N)}{\lg N} + \sqrt{N} (\lg (\lg N) - 1)$ bit

となる。

$N$ 以外の項を $N$ で割った値、 $\dfrac{1}{\lg N}, \dfrac{4 \lg (\lg N)}{\lg N}, \dfrac{\lg (\lg N) - 1}{\sqrt{N}}$ はそれぞれ $N \rightarrow \infty$ にて $0$ に収束するため、これはランダウ記号で $N + o(N)$ bitと表すことができる。以上より簡潔性が示された。 $\square$

実装

上の仕組みが理解できれば後は実装だけなのだが、これが意外と難しい。以下では個人的に詰まりやすかった箇所について軽く解説する。

元データの持ち方

まず問題となるのが「元データの保持方法」である。

boolの配列として持てば一見無駄がないように思われる。しかしC++におけるboolのサイズは環境依存で、大抵は1byteらしい。これでは8倍損してしまう。

そこで今回は、unsigned char(8bit符号なし整数)の配列として元データを表現している。data[i] >> k & 1で元データの下から $8i + k$ bit目が求まる、といった次第である。ここの幅はblockに合わせると実装がしやすくなる。

初期化方法

bit列の初期化方法といっても、

上のように配列として与える
0と1からなる文字列として与える

などの様々な方法がある。今回はこれらのどれが来ても対応できるように、「初期化時点では長さのみ指定して、後からデータを書き換える」という方式にした。全部データを書き換え終わったら、build()を呼んで前計算をしてもらう。

ただこの方法には大きな難点がある。それが「build()の呼び忘れ」である。本当に呼び忘れるので使うときは注意してほしい。

index

実装していると否が応でも0-indexedか1-indexedか困ることになる。1-indexedにすると $access$ が実装しづらくなり、0-indexedにすると $rank$ が「先頭から $n + 1$ bitの1の数」と微妙な感じになる。

これを解決する方法として、以下では $rank$ を「半開区間 $[0, n)$ に含まれる1の数」としている。これなら全体を0-indexedに統一できて、かつ $rank$ の意味も自然になる。さらに $rank$ の実装も地味に楽になる。

実装例

その他諸々難しい箇所はあるが、そこは自分で実装して悩んでいただきたい。そうすることで、簡潔ビットベクトルへの理解がより深まるはずである。

以下は $N \lt 2^{16}$ を対象とした実装例である。参考までに。

#include <vector>

using namespace std;

using u8 = unsigned char;
using u16 = unsigned short;

class bitVector {
public:
    u16 length, cnum, bnum;
    // bit列の長さ、chunk数、chunk毎のblock数

    static int cw, bw;
    // chunkの幅、blockの幅 今回は256, 8

    vector<u8> bit;  // 元データ
    vector<u16> chunk;
    vector<vector<u8>> blocks;  // この辺りの型に注意
    static vector<u8> table;  // 理想は3bit整数だが、そんなものはないので妥協


    explicit bitVector(int N) : length(N) {
        cnum = (length + cw - 1) / cw;
        bnum = cw / bw;

        bit.assign(cnum * bnum, 0);
        chunk.assign(cnum + 1, 0);
        // +1は先頭の0が入っているやつの分
        blocks.assign(cnum, vector<u8>(bnum, 0));
    }

    /* --- pos桁目をbに変える --- */
    void set(int pos, int b) {
        int bpos = pos / bw;
        int offset = pos % bw;

        if (b == 0) {
            bit[bpos] &= ~(1 << offset);
        } else if (b == 1) {
            bit[bpos] |= (1 << offset);
        }
    }

    /* --- pos桁目のbitを返す --- */
    int access(int pos) const {
        int bpos = pos / bw;
        int offset = pos % bw;
        return bit[bpos] >> offset & 1;
    }

    /* --- numの立っているbit数を返す --- */
    u8 popCount(u8 num) const {
        return table[num];
    }

    /* --- 前計算 chunkとblocksを適切に埋める --- */
    void build() {
        chunk[0] = 0;
        for (int i = 0; i < cnum; ++i) {
            blocks[i][0] = 0;
            for (int j = 0; j < bnum - 1; ++j) {
                blocks[i][j + 1] = blocks[i][j] + popCount(bit[i * bnum + j]);
            }
            chunk[i + 1] = chunk[i] + blocks[i][bnum - 1] + popCount(bit[(i + 1) * bnum - 1]);
        }
    }

    /* --- [0, pos)にある1の数を返す --- */
    int rank(int pos) const {
        int cpos = pos / cw;
        int bpos = pos % cw / bw;
        int offset = pos % bw;

        // 余った部分から、必要な部分だけをbitmaskで抽出する
        u8 masked = (bit[cpos * bnum + bpos]) & ((1 << offset) - 1);
        return chunk[cpos] + blocks[cpos][bpos] + popCount(masked);
    }

    /* rank(idx)=numなる最小のidxを返す
       そもそも1がnum個なかったら-1を返す   */
    int select(int num) const {
        if (num == 0) return 0;
        if (rank(length) < num) return -1;

        int ok = length, ng = 0;
        while (ok - ng > 1) {
            int mid = (ok + ng) / 2;
            if (rank(mid) >= num) {
                ok = mid;
            } else {
                ng = mid;
            }
        }

        return ok;
    }
};

int bitVector::cw = 256;
int bitVector::bw = 8;
vector<u8> bitVector::table = {
    0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 
    4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8, 
};

参考リンク

以下、私が簡潔ビットベクトルを実装する上で参考にしたサイトを載せておく。

完備辞書（簡潔ビットベクトル）の解説
- 実装の上で一番参考にしたサイト。 $rank$ の解説がむちゃくちゃ丁寧。
簡潔データ構造超入門～つくって学ぶ簡潔ビットベクトル～
- popcountの項で触れた「シフト演算子でpopcountする方法」が載っている。__builtin_popcountがなければ参考にどうぞ。
完備辞書
- 記事の最後にある実装例が非常に参考になった。上の実装例も大半はこれに沿っている。