prime's diary

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

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

ビット演算 競技プログラミング

はじめに

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

www.adventar.org

の1日目の記事です。

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

2021/01/03: 高速ゼータ変換/メビウス変換/畳み込みに関する部分を全面的に書き換えました。

表現方法

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) {
  // ここに処理を書く
}

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

部分集合と和に関するテクニック

2021/01/03 高速ゼータ変換/メビウス変換/畳み込みに関する部分を全面的に書き換えました。

  • 部分集合について和を取るもの

  • {\displaystyle{ g(U) = \sum_{T \supset U} f(T)}} （Uを包含するTについて和を取る: sum over supersets）

  • {\displaystyle{ g(U) = \sum_{T \subset U} f(T)}} （Uに包含されるTについて和を取る: sum over subsets）

  • 逆変換、要素数ごとに符号が反転する

  • {\displaystyle{ g(U) = \sum_{T \supset U} (-1)^{|U-T|}f(T)}} （Uを包含するTについて和を取る）

  • {\displaystyle{ g(U) = \sum_{T \subset U} (-1)^{|U-T|}f(T)}} （Uに包含されるTについて和を取る）

これらの計算は、上述の「ある部分集合の部分集合全体を渡る」テクニックにより、O(3N)で計算することができます（計算量の証明は二項定理による）。 しかし、さらに効率よく計算するアルゴリズムが存在します。これらはしばしば高速ゼータ変換/高速メビウス変換と呼ばれています。

あらかじめdp[T]にf(T)の値を入れておくことにします。 1は

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)];
    }
  }
}

2は

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

3は

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)];
    }
  }
}

4は

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

というコードでO(N2N)でできます。((T >> i) & 1) が0か1か、Tにビットを付け足すのか取り除くのか、演算が+なのか-なのかが違うだけです。 また、ここでは加算演算を利用しましたが、可換で結合的で、単位元がある（つまり、可換モノイド）であれば1と2が、さらに逆元が存在（つまり、可換群（アーベル群））すれば3,4を適用できます。

畳み込み

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

fm = ranked_sum_over_subsets(f);
gm = ranked_sum_over_subsets(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 = inverse_ranked_sum_over_subsets(conv);

ランクkのsum over subsetsは{\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) == 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) == 1) {
        dp[k][T] -= dp[k][T ^ (1 << i)];
      }
    }
  }
}

でdp[popcount(T)][T]に結果が代入されるように計算できる。 また、ここでは加算と乗算を用いましたが、他の可換環でも同様の計算を行うことができます。 (max, +)を用いるトロピカル半環の場合は、元の集合が十分小さい場合にはexp(x)(底は適切に取る)に対して(+,×)の畳み込みを行えばできます。

参考文献

さいごに

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

www.adventar.org

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