読者です 読者をやめる 読者になる 読者になる

prime's diary

そすうの日々を垂れ流しちゃうやつだよ

ビット列による部分集合表現 【ビット演算テクニック Advent Calendar 2016 1日目】

はじめに

この記事はビット演算テクニック Advent Calendar 2016

www.adventar.org

の1日目の記事です。

この記事ではビット演算を使って有限集合の部分集合を表現、操作するテクニックを紹介します。 本文中で説明を省略した部分については、詳しい説明のある参考文献を紹介しておきます。

表現方法

S={A,B,C,D}という4要素の集合があったとします。 各要素を次のように2進数の数字と対応付けます。 A: 0001 B: 0010 C: 0100 D: 1000 こうすると、Sの部分集合は存在する要素のビットごとの論理和を取ることで4ビットのビット列として一意に表すことができ、例えば {A,B,C}: 0111 {B,D}: 1010 のようになります。

以後、部分集合とそのビット列による表現をしばしば同一視します。

基本操作

空集合、全体集合

それぞれ0000, 1111が対応します。

和集合

和集合はビットごとの論理和(|)で書くことができます。

T = {A,B,C}: 0111, U = {B,D}: 1010

T∪U = {A,B,C,D}: 1111 = 0111 | 1010

共通部分

共通部分はビットごとの論理積(&)で書くことができます。

T = {A,B,C}: 0111, U = {B,D}: 1010

T∩U = {B}: 0010 = 0111 & 1010

補集合

補集合はビット反転(と必要に応じて&でマスクを取る操作)で書くことができます

T = {A,B,C}: 0111

Tc = {D}: 1000 = ~0111

(上位のビットが無視できない場合は 1000 = ~0111 & 1111)

部分集合全体を渡る

numを全体集合の要素数とすると、

for (int i = 0; i < (1 << num); ++i) {
  //ここに処理を書く
}

(本当はシフト演算子のところの括弧は不要ですが、見やすさの為付けています)

包含関係を判定する

補集合と共通部分を持つかどうかで判定できるので、 TがUに含まれるかどうかは~U&T==0で書けます。

素数を数える

1になっているビットの数を数える操作(population count, popcount)が相当します。 最近のIntelのCPUならずばりPOPCNT命令がありますし、そうでない場合も高速に処理する方法が知られています。 32要素の集合Sの部分集合Tの要素数は、

T = (T & 0x55555555) + ((T >>  1) & 0x55555555);
T = (T & 0x33333333) + ((T >>  2) & 0x33333333);
T = (T & 0x0F0F0F0F) + ((T >>  4) & 0x0F0F0F0F);
T = (T & 0x00FF00FF) + ((T >>  8) & 0x00FF00FF);
T = (T & 0x0000FFFF) + ((T >> 16) & 0x0000FFFF);

で計算できます。 ここで、 0x5=0b0101 0x3=0b0011 となっていて、 n行目で2nビットごとのpopcountを計算しています。 さらに最適化して、

T -= (T >> 1) & 0x55555555;
T = (T & 0x33333333) + ((T >> 2) & 0x33333333);
T = (T + (T >> 4)) & 0x0F0F0F0F;
T = (T * 0x01010101) >> 24;

と書くこともできます。

応用操作

ある部分集合の部分集合全体を渡る

T⊃UなるUを全部列挙します。 全体集合の部分集合全体を渡って、Tに含まれているかを判定すればできますが、要素数が少ない場合は無駄が多いですね。 実は、次のようにすることでTの要素数をNとしたとき2N回でループを抜けることができます。

for (int U = (1 << num) - 1; U >= 0; --U) {
  U &= T;
  // ここに処理を書く
}

UをTでマスクしているのがポイントです。 1を引く(デクリメントする)という操作は1になっている一番下位のビットを0にして、それより下のビットをすべて1にする、という操作なので(0のときは除く)、結果的にもれなくかぶりなく列挙できます。

ある集合を部分集合に持つ集合全体を渡る

T⊂UなるUを全部列挙します。

for (int i = T; i < (1 << num); i=(i+1)|T) {
  // ここに処理を書く
}

仕組みは部分集合列挙とほとんど同じです。

素数kの部分集合全体を渡る

int x, y;
for (int T = (1 << k) - 1; T < (1 << N); x = T & -T, y = T + x, T = (((T & ~y) / x) >> 1) | y) {
  // ここに処理を書く
}

かなりトリッキーですね。 詳しくは参考文献の解説を読んでください。

高速ゼータ変換

Sの部分集合を引数に取る関数fがあったとして、 {\displaystyle{ g(U) = \sum_{T \supset U} f(T)}} をSの各部分集合Uに対して求めたいとします。 各Uについて個別に求めれば二項定理により全体の計算量はO(3N)となりますが、Tの取りうる数が2N通りしかないので無駄が多いです。 これを効率よく計算する高速ゼータ変換というアルゴリズムがあって、 あらかじめdp[T]にf(T)の値を入れておくと、

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if (((T >> i) & 1) == 0) {
      dp[T] += dp[T | (1 << i)];
    }
  }
}

というコードでO(N2N)でできます。

高速メビウス変換・逆変換

Sの部分集合を引数に取る関数fがあったとして、 [tex:{\displaystyle{ g(U) = \sum{T \subset U} f(T)}}] をメビウス変換(ゼータ変換と包含の向きが逆)、 [tex:{\displaystyle{ g(U) = \sum{T \subset U} (-1)^{|U-T|}f(T)}}] をメビウス逆変換といいます。(集合の引き算は差集合。) その名の通りメビウス変換して逆変換すると元に戻ります。

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if ((T >> i) & 1) {
      dp[T] += dp[T ^ (1 << i)];
    }
  }
}

逆変換は差集合の要素数に依存する符号があって複雑そうですが、実はこちらも似たようなコードでできます。(符号が逆になっただけ)

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if ((T >> i) & 1) {
      dp[T] -= dp[T ^ (1 << i)];
    }
  }
}

畳み込み

Sの部分集合を引数に取る関数f,gがあったとして、 {\displaystyle{ (f \circ g)(U) = \sum_{T \subset U} f(T)g(U-T)}} をSの各部分集合Uに対して求めたいとします。 (ランク付きの)高速メビウス変換・逆変換を用いると計算できます(多項式乗算におけるFFT/逆FFTと似た理屈)

fm = ranked_mobius(f);
gm = ranked_mobius(g);
for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < (1 << N); ++i) {
    for (int j = 0; j <= k; ++j) {
      conv[k][i] += fm[j][i]*gm[k-j][i];
    }
  }
}
result = ranked_mobius_inv(conv);

ランクkのメビウス変換は{\displaystyle{ \hat{f}(k,U) = \sum_{T \subset U, |T| = k} f(T)}}

実際のコードはこんな感じ(dp[popcount(T)][T]にはf(T)の値を入れておく)

for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < N; ++i) {
    for (int T = 0; T < (1 << N); ++T) {
      if ((T >> i) & 1) {
        dp[k][T] += dp[k][T ^ (1 << i)];
      }
    }
  }
}

メビウス逆変換は{\displaystyle{ f(U) = \sum_{T \subset U} (-1)^{|U-T|} \hat{f}(|U|,T)}}

dp[k][T]{\hat{f}(k,T)}の値を入れておくと、

for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < N; ++i) {
    for (int T = 0; T < (1 << N); ++T) {
      if ((T >> i) & 1) {
        dp[k][T] -= dp[k][T ^ (1 << i)];
      }
    }
  }
}

でdp[popcount(T)][T]に結果が代入されるように計算できる。(実は畳み込みしなくてももともと結果が入っている)

参考文献

さいごに

ビット演算テクニック Advent Calendar 2016

www.adventar.org

みんな書いて!おねがい!