はじめに

この記事は、「bit列を用いたbit全探索」を前提知識としています。それ以外は基本的に記事中で説明するように心がけています。
もし不安がある場合は、この記事を読むと理解が進む可能性があります。

また、記事中のサンプルコードは基本的にC#で書かれていますが、できる限り平易な書き方を心がけています。
本記事で説明するテクニックで落ちる計算量は、現実的には高々$20$倍程度の改善にしかなりません。このオーダーだと定数倍レベルが大事になってくると感じているため、そこも意識した記事構成としています。ご了承ください。

bit全探索

突然ですが、bit全探索を用いて以下の問題を解いてみましょう。

問題

$A_1,A_2 \cdots A_N$という数列があります。この数列から数をいくつか選び、その総和が$M$に最も近くなるような選び方をします。そのときのMとの絶対値を求めて下さい。
【制約】
$1\leqq N\leqq 20$
$1\leqq A_i,M\leqq {10}^8$

この問題はNが十分に小さいため、$2^N$個¹存在する選び方全てを列挙することで解くことができます。それぞれの集合において$N$bitを走査して和を計算するため、計算量は$O(N\cdot 2^N)$となります。

int res = INF;
//全ての状態を列挙
for (int i = 0; i < (1 << N); i++)
{
    int sum = 0;
    //現在の状態iにおいての和を計算
    for (int j = 0; j < N; j++)
        if (((i >> j) & 1) == 1) sum += A[j];
    //より良かったらresを更新
    res = Min(res, Abs(M - sum));
}

$0 \sim 2^N-1$までの全ての数字を全探索することで、$N$要素の集合の存在する選び方の列挙をしています。
これの計算量の$N$が落とせて$O(2^N)$にできるよ、というのが本記事の目的です。

計算量を落とす

上のコードでは、それぞれの状態の時に$N$ステップで和を計算しています。上手いことをするとbitが立っている部分だけを抽出して平均$\frac{N}{2}$ステップで列挙が可能になったりする²んですが、今回はそういう話ではありません。

例えば、$100$という状態から$101$という状態に変化するときを考えてみましょう。この時、変わったのは最下位の1bitのみです。これは$A_0$を示しています。
$x$という状態の和を計算する処理を$CalcSum(x)$とすると、$CalcSum(100)+A_0$と$CalcSum(101)$は同一となります。
また、$101$という状態の和からも$100$の状態の和を求めることができます。これは$CalcSum(101)-A_0$とすれば良く、減算が加算の逆演算となっているため成立しています。また、加算は交換法則が成立するため、操作のすぐ後でなくとも逆元を適用することが可能になっています。
このように前の結果を利用し、次の結果を導出することを差分更新と呼ぶこととすると、

bitの違う部分が定数($d$)個
bitの異なる位置が高速に求まる

という条件で状態を全列挙できれば、$O(d)$で差分更新が可能になることがわかります。
ここで、定数$d$はどんなに桁数が大きくても$1$とできます。これが最も効率的であることは明らかなため、今後はこの$d=1$のケースのみ考えることとします。
先程のコードの状態の列挙順は通常の2進数と同じでしたが、上手い列挙方法を使うことで差分更新を高速で可能とします。この列挙順は複数ありますが、そのうちの一つをグレイコードと呼びます。

グレイコード

とはなにか

上にも述べた通り、グレイコードとは2進数における数字の表し方の一つで、複数の嬉しい性質を持っています。[^4]具体的には、i番目のグレイコードはi^(i>>1)というbit演算による結果と対応します。ここで、^はxor、>>はbitを1つ右にシフトすることを示します。

XOR、ビットシフトについて

XORは「2つのbitが異なるならば1、そうでないならば0」という動作をします。具体的には、`1100^1010==0110`となります。また、ビットシフトはその名の通り「すべてのbitを右/左に指定個ずらす」という動作をします。余ったところに何を埋めるのか、正負を表す符号のbitはどうするか等で複数種類存在しますが、今回は対象が必ず正なのであまり意識しなくて良いです。具体的には、`101101>>1==010110`となります。以上の解説はざっくりとしたものなので、より正確なものは[Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%93%E3%83%83%E3%83%88%E6%BC%94%E7%AE%97)等を見てください。これらの演算はあまり馴染みがないかもしれません。しかし、コンピュータでは、これらの演算は加算と同程度、またはそれ以上に高速に行える演算になっています。

[^4]: 先程は「列挙順」と書きましたが、これは数一つ一つにコードを対応づけることと意味としては等しいです。

通常の2進数では、

…000
…001
…010
…011
…100

と数字が続きますが、グレイコードでは

…000
…001
…011
…010
…110

と続きます。通常の2進数と比較をし、上記の一般項i^(i>>1)がしっかりと成立していること、「隣り合うコードでは、立っているbitは1しか変わらない」という性質を確認してください。
更に、i番目とi+1番目のグレイコードの異なっているbitは、i+1という2進数の最下位bit(Least Significant Bit,LSB)³であるという性質を持ちます。

[](

以上の性質の証明

\ \

\ )

実装

ソースコードの通りです。できる限りすべての処理に注釈を入れるようにしましたが、分からないところがあれば聞いてください。

//初期状態はすべてが0の状態
int sum = 0;
//初期状態のスコアであるところの、Abs(M-0)を格納
int res = M;
//iを1からはじめている。
//なぜならば、i==0の時のグレイコードは0となり、差分更新するべき場所が存在しないから。
for (int i = 1; i < (1 << N); i++)
{
    //i番目のグレイコード
    int code = i ^ (i >> 1);

    int bitpos = 0;
    //下から見ていき、その桁が0であれば一桁見る場所を上げる、を繰り返す。
    //iがbitを含まないのはi==0の時のみなので、このループは停止する。
    while (((1 << bitpos) & i) == 0) bitpos++;

    //変化したことが確実なbitの現在の状況を観測することで、どのように変化したかは分かる。
    //グレイコード上でそのbitが立っていた場合は(0→1に変化した場合なため)演算を、
    if (((1 << bitpos) & code) != 0) sum += A[bitpos];
    //立っていなかった場合は(1→0に変化した場合なため)逆演算を適用する。
    else sum -= A[bitpos];

    //より良かったらresを更新
    res = Min(res, Abs(M - sum));
}

aのLSBの場所についてはループで愚直に求めていますが、この方法でも定数時間の平均計算量⁴で求めることができます⁵。

まとめ

bit全列挙が可能な場合、

(高速に行える)逆演算

さえ存在していれば、グレイコードを用いた全列挙で高速化をすることがが可能になります。

余談1 : LSBの場所の計算方法について

今回はループを用いて平均定数時間で計算を行いましたが、この他にも実装方法は複数存在します。

方針1 : CPU組み込み関数を使う

C++などの言語では、一部のCPUに備わった特殊な命令を使用することができる場合があります。
最も代表的なのはGCCのbuiltin関数です。これには、最下位bitから連なっている0の数(ctz:count trailing zeros)を求めるctzという命令が存在します。これは0-indexedにおけるLSBの位置と同じことが分かります。よって、以下のように書くことができます。

int bitpos = __builtin_ctz(i);

この関数はCPUの一命令で実行されるため、定数時間となります。

方針2 : 頭のいいビット演算の黒魔術を使う

LSBを求める

まず、LSB自体は簡単なビット演算で求めることができます。ある2進数iのLSBは、2の補数を用いて負数を表現している場合はi & -iとすることで求めることができます。

LSBの位置を求める

今、このLSBはbitが一つだけ立っている状態、つまり1 << Nという表記で表すことができます。
また、このNが求めたいLSBの位置であることも分かります。
ここで、a * (1 << n) == a << nという事実が実はあります。2進数の掛け算も通常の筆算のように行うことができます。1 << nは一桁しか1の桁のがないため繰り上がりが起きません。また、その桁とaを掛けたものは、aをnビットシフトしたものと同じであることが分かります。詳しくは以下の図を見てください。

さて、De Bruijn Sequenceというものがこの世の中には存在しています。追い追いこれについても記事を一本書きたいと思っていますが、どういったものかをとりあえず簡単に解説します。

突然ですが問題です。
0,1の文字を並べて文字列を作ります。長さNの連続した部分文字列が、有り得る$2^N$通りの文字列と過不足なく対応しているような文字列を構築してください。
サンプル1
N=3
0010111000
全ての連続した部分文字列は001,010,101,011,111,110,100,000となり、条件を満たします。

これの答えが(次数2,サイズN)のDe Bruijn Sequenceです。^8
これの何が嬉しいかというと、これを桁溢れギリギリまで詰めた場合、$0\sim 2^N-1$個左シフトした時の上位bit$N$桁が互いに相違なものになります。
とりあえず、上のサンプルで試してみます。これは次数2,サイズ3のDe Bruijn Sequenceです。ビットシフトは不足した桁を0埋めするため、末尾の0は不要となります。これにより、0010_1110と$2^3$bitに収めることができました。これの上位3bitを見てみます。

0010_1110 << 0 == (001)0_1110
0010_1110 << 1 == (010)1_1100
0010_1110 << 2 == (101)1_1000
0010_1110 << 3 == (011)1_0000
0010_1110 << 4 == (111)0_0000
0010_1110 << 5 == (110)0_0000
0010_1110 << 6 == (100)0_0000
0010_1110 << 7 == (000)0_0000

括弧でくくったところが注目すべき場所で、ここは先程もいった通り互いに異なっています。シフトで溢れた元々あったさらに上位のbitは切り詰められるため、この上位bitは連続する部分文字列の列挙と一致するのです。後はこれらを右にいくらかシフトしてあげれば、下位ビットも切り詰めることができます。(以下は切り詰めた一例です)

(0010_1110 << 3) >> 5 == 011

さて、これが一意に定まったということは、シフトされた数字からシフト幅を求めることができるということです。方法は簡単で、この結果を添え字とした配列を作って上げればよいだけです。この配列は、

//DeBruijnSequenceの次数
int Order = 3;
for(int i = 0; i < (1 << Order); i++) DeBruijnBitPosition[(DeBruijnSequence << i) >> ((1 << Order) - Order)] = i;

といった形で生成することができます。
よって、最初に示したa * (1 << n) == a << nという情報より、

//サイズ5のDe Bruijn Sequence
const uint DeBruijnSequence = 0b0000111011010111110011000101001U;
//DeBruijnBitPosition[i] := 上のDe Bruijn Sequenceにて、対応しているシフト幅
int[] DeBruijnBitPosition = { 0, 1, 23, 2, 29, 24, 19, 3, 30, 27, 25, 11, 20, 8, 4, 13, 31, 22, 28, 18, 26, 10, 7, 12, 21, 17, 9, 6, 16, 5, 15, 14 };

uint lsb = (uint)(i & -i);
int bitpos = DeBruijnBitPosition[(DeBruijnSequence * lsb) >> 27];

とすることでシフト量、つまりはLSBの場所を定数時間で計算することができます。

余談2 : 再帰の実装

今回再帰は出てきませんでしたが、再帰を用いると逆演算がなくとも探索が可能です。計算量は$O(2^N)$ですが、実はこの再帰関数は$2^{N+1}-1$回呼ばれます⁶。

res = DFS(0, 0);
int DFS(int sum, int depth)
{
    if (depth == N)
        return Abs(M - sum);
    return Min(DFS(sum, depth + 1), DFS(sum + A[depth], depth + 1));
}

再帰で逆演算が不要なのは、値をコピーして渡しているためです。これにより、逆演算による状態の復元が不要となっています。⁷
なので、値が高速にコピーできない場合(配列になっている等)は計算量が余計に掛かってしまうことがあり、グレイコードを用いた方が有利です。また、再帰を用いると関数呼び出しのオーバーヘッドなどの影響で少し遅くなりがちです。⁸

	愚直bit全探索	グレイコードbit全探索	再帰
オーダー	$O(N\cdot 2^N)$	$O(2^N)$	$O(copy\cdot 2^N)$
定数倍	軽い	ちょっと重い	重い
必須性質	なし	可逆性,交換法則	コピー可能

どれを使うかはケースバイケースだと思います。

それぞれの要素を選ぶ/選ばないで2通りあり、要素がN個存在しているので$2^N$となります。
これを実現するコードは以下の通りです。 ↩
少し後にも出てきますが、最下位bitは高速に取得することが可能です。最下位bitを取得し、そこを消す…を繰り返すと、bitの数だけのイテレーションで済むようになります。bitは半分の確率で立っているため、計算量は二分の一となります。 ↩
bit毎に見た際に、最も右側(小さい)場所にある1のbitのこと。0101であれば0001、1100であれば0100等。
以上の性質を用いて、差分更新を実装していきます。 ↩
一回の演算についての最悪ケースでは定数時間に収まらない可能性がありますが、クエリがランダムに飛んでくると仮定した場合に、十分な試行回数を積めば平均で一定の時間に収まる、という意味です。 ↩
それぞれのbitについて、下からi番目のbitを探索が通過する(=そこまでがすべて0)確率は$2^{-i}$で、終了に必要な探索回数は$1$。$1+{\sum_{k=1}^{inf} 2^{-k}}=2$より、定数に収束します。
よって、この実装の計算量は$O(2^N)$となります。 ↩
末端ノード($depth=N$となる場所)が$2^N$個、それ以外のノードが$2^N-1=\sum_{k=0}^{N-1} 2^k$個なため。セグ木の要素数が$2^{(N+1)}-1$となる理由と同一です。 ↩
再帰関数を抜ける際に逆元を適用して状態を元に戻すことにより、再帰による実装でもコピーコストを掛けずに$O(2^N)$での計算が可能になります。 ↩
コピーが高速であることを利用した$O(2^N)$化は、再帰以外の実装でも可能です。配列に全てを保存しておき、LSBを除いたところから更新をすれば良いです。
以上をまとめると、以下の表のようになります。 ↩