Mister雑記

競プロやります。

高速ゼータ・メビウス変換

ググれば大量の記事が出てくるが、知見の整理ということで自分でも書いてみる。

ゼータ・メビウス変換とは

「高速」ゼータ・メビウス変換とは名前の通りゼータ・メビウス変換を高速化したものなので、まずはそちらを理解する必要がある。

定義: ゼータ変換、メビウス変換

 Xを有限集合、 Gを加法が定義された集合1とする。

このとき写像 f: {\cal P}(X) \to Gから、以下を満たす写像 g: {\cal P}(X) \to Gを求めることをゼータ変換という2

 \displaystyle g(S) = \sum_{T \supseteq S} f(T)

なお、この式中の T \supseteq S T \subseteq Sにしたものも同様にゼータ変換と呼ぶことにする3

そしてこの逆変換、すなわち g: {\cal P}(X) \to Gから以下を満たす f: {\cal P}(X) \to Gを求めることをメビウス変換という4

 \displaystyle f(S) = \sum_{T \supseteq S} (-1)^{|T| - |S|}g(T)

この式中の T \supseteq S T \subseteq Sにすると、上で包含関係の向きを逆にした場合のゼータ変換の逆変換になる。

この変換がゼータ変換の逆変換であることは全く自明ではないが、この記事ではその証明は省略する5


以降、集合 Xとして \{ 1, 2, \cdots, n \}のみを考え、関数 fを長さ 2^nの配列として表現する。これと対応して、以降 f(\{ 2, 4, 5 \}) f(11010)のように表現する。

ゼータ変換の高速化

ここからゼータ変換、特に実装しやすい g(S) = \sum_{T \subseteq S} f(T)の方を高速に行うことを考える。メビウス変換も符号を少し変えるだけでほぼ同様に実現できる。

計算量の評価は、加算が行われる回数 C (Complexity)によって行う。

愚直解

まず愚直に求めようとすると、以下のようなコードになるだろう。

for (int s = 0; s < (1 << n); ++s) {
    for (int t = 0; t < (1 << n); ++t) {
        if ((s & t) == t) g[s] += f[t];
    }
}

 S T Xの部分集合全体を回り、 S \cap T = T、すなわち T \subseteq Sのときだけ g(S) f(T)を加算する。定義通りである。

 Xの部分集合は計 2^n個あることから、このコードの計算量は C = 2^n \cdot 2^n = 4^nとなる。

部分集合の列挙

しかしこのコードでは Tの探索に無駄が多すぎる。できることなら Xではなく Sの部分集合全体だけを回したいが、実はこれは以下のコードで実現できる。

for (int s = 0; s < (1 << n); ++s) {
    for (int t = s; t >= 0; t = ((t - 1) & s)) {
        g[s] += f[t];
    }
}

この手法はビット演算 (bit 演算) の使い方を総特集! 〜 マスクビットから bit DP まで 〜に載っている。

 Sの部分集合は 2^{|S|}個あり、 Xのサイズ kの部分集合は {}_n C_k個あることから、このコードの計算量は C = \sum_{k = 0}^{n} {}_n C_k \cdot 2^k = 3^nとなる6

高速ゼータ変換

さらに高速化して、 C = n \cdot 2^{n - 1}でゼータ変換を行うのが本記事の主題である高速ゼータ変換と呼ばれる手法である。

方針は「下位要素から順に拡張していく」というもの。 まず dp_{S, k}を「 1 \sim k Sの部分集合、それより上は Sと一致しているもの全部の総和」と定義する。言葉では分かりにくいので例を挙げると、

 \displaystyle dp_{\color{red}{110}\color{blue}{101}, 3} =
f(\color{red}{110}\color{blue}{000}) +
f(\color{red}{110}\color{blue}{001}) +
f(\color{red}{110}\color{blue}{100}) +
f(\color{red}{110}\color{blue}{101})

といった具合である。定義より、 dp_{S, 0} = f(S) dp_{S, n} = g(S)

これを k = 1から k = nまで順に更新していく。このとき、 k \in Sか否かで場合分けをする。具体例で考えてみると、

 \displaystyle dp_{\color{red}{10}\color{green}{0}\color{blue}{10}, 3} =
f(\color{red}{10}\color{green}{0}\color{blue}{00}) +
f(\color{red}{10}\color{green}{0}\color{blue}{10}) =
dp_{\color{red}{10}\color{green}{0}\color{blue}{10}, 2} \\
dp_{\color{red}{10}\color{green}{1}\color{blue}{10}, 3} =
f(\color{red}{10}\color{green}{0}\color{blue}{00}) +
f(\color{red}{10}\color{green}{0}\color{blue}{10}) +
f(\color{red}{10}\color{green}{1}\color{blue}{00}) +
f(\color{red}{10}\color{green}{1}\color{blue}{10}) =
dp_{\color{red}{10}\color{green}{0}\color{blue}{10}, 2} +
dp_{\color{red}{10}\color{green}{1}\color{blue}{10}, 2}

のように、 k \in Sの場合は dp_{S \backslash \{ k \}, k - 1}も足す必要がある。

これを実装に落とし込むと以下の通り。

for (int s = 0; s < (1 << n); ++s) {
    dp[s][0] = f[s];
}

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

for (int s = 0; s < (1 << n); ++s) {
    g[s] = dp[s][n];
}

集合及びDPの添字は1-indexedなのに対してbit演算は0-indexedなので、シフト数を1減らしてやる必要がある。

配列の使い回し

さらに更新について考えてみると、DPテーブルは1次元のものを使い回せることに気づく。すなわち、以下のような実装ができる。

for (int s = 0; s < (1 << n); ++s) {
    dp[s] = f[s];
}

for (int k = 1; k <= n; ++k) {
    // この時点でdp[s]は上の実装のdp[s][k-1]と一致
    for (int s = 0; s < (1 << n); ++s) {
        if ((s >> (k - 1)) & 1) {
            dp[s] += dp[s ^ (1 << (k - 1))];
        }
    }
}

for (int s = 0; s < (1 << n); ++s) {
    g[s] = dp[s];
}

ここで「 dp_{S \backslash \{ k \}} dp_{S}より先に更新されているが大丈夫なのか」となるが、 k \not \in S \backslash \{ k \}なので、この週ではそもそも dp_{S \backslash \{ k \}}は更新されていない。よって問題ない。

さらにkを1ずらしたりDP配列として直にgを使ったりして整理することで、他の記事でもよく見られるような実装になる。

for (int s = 0; s < (1 << n); ++s) {
    g[s] = f[s];
}

for (int k = 0; k < n; ++k) {
    for (int s = 0; s < (1 << n); ++s) {
        if ((s >> k) & 1) {
            g[s] += g[s ^ (1 << k)];
        }
    }
}

長い道のりだったが、ようやく高速ゼータ変換のコードに辿り着いた。

このコードの計算量は、各 S \subseteq Xについて |S|回加算更新が行われるため、 C = \sum_{k = 0}^{n} {}_n C_k \cdot k = n \cdot 2^{n - 1}となる7。加算が O(1)なら、 n = 18くらいなら安定して間に合うだろう( C \simeq 2.4 \times 10^6)。

高速ゼータ・メビウス変換の実装例まとめ

最後にまとめとして、各種変換の実装例を載せておく。 なお配列 f, gはコード外で宣言済みとする。

(注※) verifyはしていないので、コピペは自己責任で。実装ミスがあれば報告していただけると助かります。

高速ゼータ変換 上位集合版

 g(S) = \sum_{T \supseteq S} f(T)

for (int s = 0; s < (1 << n); ++s) {
    g[s] = f[s];
}

for (int k = 0; k < n; ++k) {
    for (int s = 0; s < (1 << n); ++s) {
        if (!((s >> k) & 1)) {
            g[s] += g[s ^ (1 << k)];
        }
    }
}

先程とは逆に、 k \not \in Sのときに dp_{S \cup \{ k \}}を加える。この遷移は Sが大きい方から小さい方に加えられるので、更新順序的に問題ない。

ゼータ変換 下位集合版

 g(S) = \sum_{T \subseteq S} f(T)

for (int s = 0; s < (1 << n); ++s) {
    g[s] = f[s];
}

for (int k = 0; k < n; ++k) {
    for (int s = 0; s < (1 << n); ++s) {
        if ((s >> k) & 1) {
            g[s] += g[s ^ (1 << k)];
        }
    }
}

先程長々と解説した通り。

メビウス変換 上位集合版

 f(S) = \sum_{T \supseteq S} (-1)^{|T| - |S|} g(T)

for (int s = 0; s < (1 << n); ++s) {
    f[s] = g[s];
}

for (int k = 0; k < n; ++k) {
    for (int s = 0; s < (1 << n); ++s) {
        if (!((s >> k) & 1)) {
            f[s] -= f[s ^ (1 << k)];
        }
    }
}

非自明ではあるが、実はゼータ変換の更新の符号を反転させるだけでいい。  kが抜けると dp_{S \cup \{k\}}の各項について |S|が1減るので符号が反転する、みたいなお気持ち。

メビウス変換 下位集合版

 f(S) = \sum_{T \subseteq S} (-1)^{|S| - |T|} g(T)

for (int s = 0; s < (1 << n); ++s) {
    f[s] = g[s];
}

for (int k = 0; k < n; ++k) {
    for (int s = 0; s < (1 << n); ++s) {
        if ((s >> k) & 1) {
            f[s] -= f[s ^ (1 << k)];
        }
    }
}

先と同じく符号を反転させるだけでいい。


  1. ゼータ変換では、 \sumを使うため Gは可換半群であるべきだと思われる。しかし代数学にはあまり自信がないので注釈に留めておく。

  2. 別な表現として、 g = \zeta fを満たすような写像 \zeta: G^{{\cal P}(X)} \to G^{{\cal P}(X)}をゼータ変換、(存在すれば)その逆写像 \zeta^{-1}メビウス変換と呼ぶこともできる。が、却って分かりにくくなりそうなので割愛。

  3. この辺は記事によってはメビウス変換になっていたりして曖昧な部分だが、ここでは高速ゼータ変換/高速メビウス変換に従うことにする。

  4. メビウス変換では逆元を必要とするため、 Gはより強い性質をもつアーベル群である必要があると思われる。

  5. 指数時間アルゴリズム入門の53ページに式変形が載っている。 f gが逆になっているので注意。

  6. 2つ目の等号は、二項定理で (1 + 2)^nを展開すると成り立つことが分かる。

  7. この2つ目の等号を示す方法として、例えば (1 + x)^nを二項定理を使って展開し、辺々を微分してから x = 1を代入するなどがある。