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

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

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

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

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

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

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

g(S)=∑T⊇Sf(T)${\displaystyle g(S)=\sum _{T\supseteq S}f(T)}$

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

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

f(S)=∑T⊇S(−1)|T|−|S|g(T)${\displaystyle f(S)=\sum _{T\supseteq S}(-1{)}^{|T|-|S|}g(T)}$

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

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

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

ゼータ変換の高速化

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

計算量の評価は、加算が行われる回数C$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$S$とT$T$はX$X$の部分集合全体を回り、S∩T=T$S\cap T=T$、すなわちT⊆S$T\subseteq S$のときだけg(S)$g(S)$にf(T)$f(T)$を加算する。定義通りである。

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

部分集合の列挙

しかしこのコードではT$T$の探索に無駄が多すぎる。できることならX$X$ではなくS$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$S$の部分集合は2|S|$2^{|S|}$個あり、X$X$のサイズk$k$の部分集合はnCk${}_n{C}_k$個あることから、このコードの計算量はC=∑nk=0nCk⋅2k=3n$C=\sum _{k=0}^n{}_n{C}_k\cdot 2^k=3^n$となる⁶。

高速ゼータ変換

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

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

dp110101,3=f(110000)+f(110001)+f(110100)+f(110101)${\displaystyle d{p}_{110101,3}=f(110000)+f(110001)+f(110100)+f(110101)}$

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

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

dp10010,3=f(10000)+f(10010)=dp10010,2dp10110,3=f(10000)+f(10010)+f(10100)+f(10110)=dp10010,2+dp10110,2$$\begin{gathered} {\displaystyle d{p}_{10010,3}=f(10000)+f(10010)=d{p}_{10010,2}} \\ d{p}_{10110,3}=f(10000)+f(10010)+f(10100)+f(10110)=d{p}_{10010,2}+d{p}_{10110,2} \end{gathered}$$

のように、k∈S$k\in S$の場合はdpS∖{k},k−1$d{p}_{S\mathrm{\setminus }\{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];
}

ここで「dpS∖{k}$d{p}_{S\mathrm{\setminus }\{k\}}$はdpS$d{p}_S$より先に更新されているが大丈夫なのか」となるが、k∉S∖{k}$k\notin S\mathrm{\setminus }\{k\}$なので、この週ではそもそもdpS∖{k}$d{p}_{S\mathrm{\setminus }\{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⊆X$S\subseteq X$について|S|$|S|$回加算更新が行われるため、C=∑nk=0nCk⋅k=n⋅2n−1$C=\sum _{k=0}^n{}_n{C}_k\cdot k=n\cdot 2^{n-1}$となる⁷。加算がO(1)$O(1)$なら、n=18$n=18$くらいなら安定して間に合うだろう(C≃2.4×106$C\simeq 2.4\times {10}^6$)。

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

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

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

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

g(S)=∑T⊇Sf(T)$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∉S$k\notin S$のときにdpS∪{k}$d{p}_{S\cup \{k\}}$を加える。この遷移はS$S$が大きい方から小さい方に加えられるので、更新順序的に問題ない。

ゼータ変換 下位集合版

g(S)=∑T⊆Sf(T)$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)=∑T⊇S(−1)|T|−|S|g(T)$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$k$が抜けるとdpS∪{k}$d{p}_{S\cup \{k\}}$の各項について|S|$|S|$が1減るので符号が反転する、みたいなお気持ち。

メビウス変換 下位集合版

f(S)=∑T⊆S(−1)|S|−|T|g(T)$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$G$は可換半群であるべきだと思われる。しかし代数学にはあまり自信がないので注釈に留めておく。↩

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

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

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

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

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

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