$\gdef\all#1{\lbrack #1\rbrack}$
$\gdef\powall#1{2 ^ {\all{#1}}}$
$\gdef\set#1{\lbrace#1\rbrace}$
$\gdef\F{\mathcal{F}}$
追記
本記事の内容に関してより統一的な見方をしていそうな記事を教えて頂きました。
https://codeforces.com/blog/entry/92183
追記 (2023/05/02): 上記事 (多項式 との合成) に関する内容を追加しました。
導入
非負整数 $n$ に対して $\all{n} \coloneqq \set{0,1,\ldots,n-1}$ と定義します。
集合関数 $f\colon \powall{n} \to \mathbb{K},\ g\colon \powall{n} \to \mathbb{K}$ の subset convolution $f\ast g$ を以下で定義します。
$$f\ast g\colon S\in \powall{n} \mapsto \sum _ {X\subseteq S} f(X) g(S\setminus X) \in\mathbb{K}.$$
与えられた $f,g$ に対して $f\ast g$ を $O(2 ^ n n ^ 2)$ 時間で計算するアルゴリズム が知られています*1 。
$R _ n$ を $\powall{n}$ から $\mathbb{K}$ への関数全体の集合として $(R _ n, +, \ast)$ は可換環 を成します。ここで、加法 $+$ は通常の加法で $f+g\colon S\in \powall{n} \mapsto f(S) + g(S)\in\mathbb{K}$ です。
加法単位元 $e _ 0$ および乗法単位元 $e _ 1$ はそれぞれ次で表されます。
$$\begin{aligned}
e _ 0 &\colon S\in \powall{n} \mapsto 0 \in\mathbb{K}, \\
e _ 1 &\colon S \in \powall{n} \mapsto \begin{cases} 1 & \text{if}\quad S =\emptyset \\ 0 & \text{otherwise} \end{cases}.
\end{aligned}$$
本記事では、この環上での様々な操作をダブリングによって計算する方法を与えます。
ranked zeta transform -> 各点に対する演算 -> ranked mobius transform とする計算方法 *2 もありますが、ダブリングによる実装は subset convolution さえ書いてしまえば fps exp / log / pow などを書かずに済むという利点があります。ただ、後述しますが一部のケースはダブリングで対応することが (おそらく) 難しいことに注意してください。
ダブリングの基本方針
計算結果 $h\in R _ n$ を得るアルゴリズム を以下の手順で得ます。
$h(\emptyset)$ を得る
$m\in\all{n}$ を固定し、$S\subseteq \all{m}$ なる全ての $S$ に対して $h(S)$ が求まっていると仮定して $m\in T\subseteq \all{m + 1}$ なる全ての $T$ に対する $h(T)$ を計算する方法を得る
これをダブリングと呼んでいるのは、$m=0,1,\ldots,n-1$ の順に 2 を行うことで既知の $h(\bullet)$ の個数が倍々に増えていくことが理由です。
記法
集合 $S, T$ に対して $S \sqcup T$ で $S\cap T = \emptyset$ であるような $S, T$ に対する $S\cup T$ を表します (disjoint であることを明示的に示すための記法)。
表記の簡単のため、$x\in\all{n}$ に対して $S + x$ で $S \sqcup \set{x}$ を表します。
$m\in\all{n}$ と $f\in R _ n$ に対して $\mathsf{lo} _ m (f)$ および $\mathsf{hi} _ m (f)$ を以下で定義します。
$$\begin{aligned}
\mathsf{lo} _ m (f) &\colon S\in\powall{m}\mapsto f(S) \in \mathbb{K}, \\
\mathsf{hi} _ m (f) &\colon S\in\powall{m}\mapsto f(S+m) \in \mathbb{K}.
\end{aligned}$$
inv
$f\in R _ n$ に対して $f ^ {-1}$ を $f\ast g = e _ 1$ を満たす $g\in R _ n$ と定義します。ただし $f(\emptyset) \neq 0$ を仮定します。
$f ^ {-1}(\emptyset)$
$f ^ {-1}(S+m)$ for all $S\subseteq \all{m}$
計算量を $T(n)$ とすると、$T(n) = T(n - 1) + O(2 ^ n n ^ 2)$ より $T(n) = O(2 ^ n n ^ 2)$ です。
exp
$\displaystyle \exp f \coloneqq \sum _ {k = 0} ^ \infty \dfrac{f ^ k}{k!}$ と定義されます。ただし $f(\emptyset) = 0$ を仮定します。($f(\emptyset) \neq 0$ であっても $\exp (f(\emptyset) )$ が $\mathbb{K}$ において定義されるならば正当そうなんですが、一旦正当性が簡単でかつ実用上でもよく用いられるであろう $f(\emptyset) = 0$ のケースのみを考えます)
$(\exp f)(\emptyset)$
$(\exp f)(S+m)$ for all $S\subseteq \all{m}$
計算量は inv のときと同様の解析により $O(2 ^ n n ^ 2)$ 時間です。
log
$\log$ は $\exp$ の逆関数 なので、次が成り立ちます。ただし $f(\emptyset) = 1$ を仮定します。
$$\begin{aligned}
(\log f)(\emptyset) &= \log (f(\emptyset)) = 0, \\
\mathsf{hi} _ m(\log f) &= \mathsf{hi} _ m(f) \ast \mathsf{lo} _ m(f) ^ {-1} = \mathsf{hi} _ m(f) \ast \mathsf{lo} _ m(f ^ {-1}).
\end{aligned}$$
計算量は $O(2 ^ n n ^ 2)$ 時間です。
pow
$k$ を非負整数として $f ^ k$ の計算を考えます。$k = 0$ のときは $f ^ k = e _ 1$ なので、以下では $k \gt 0$ を仮定します。
今までとは少し異なる方針で計算します。$S\subseteq\all{n - 1}$ に対して、次が成り立ちます。
$$\begin{aligned}
f ^ k(S + (n - 1)) &= \sum _ {X _ 1\sqcup X _ 2\sqcup \cdots \sqcup X _ k = S + (n - 1)} \prod _ {i = 1} ^ k f(X _ i) \\
&= k\sum _ {Y\subseteq S} f(Y + (n - 1)) \sum _ {X _ 1\sqcup X _ 2\sqcup \cdots \sqcup X _ {k - 1} = S \setminus Y} \prod _ {i = 1} ^ {k - 1} f(X _ i) \\
&= k\times (\mathsf{hi} _ {n - 1}(f) \ast \mathsf{lo} _ {n - 1}(f) ^ {k - 1}) (S).
\end{aligned}$$
つまり、$\mathsf{hi} _ {n - 1} (f ^ k) = k\times(\mathsf{hi} _ {n - 1}(f) \ast \mathsf{lo} _ {n - 1}(f) ^ {k - 1})$ です。一方で次が成り立ちます。
$$\mathsf{lo} _ {n - 1}(f ^ k) = \mathsf{lo} _ {n - 1}(f) ^ k = \mathsf{lo} _ {n - 1}(f) \ast \mathsf{lo} _ {n - 1}(f) ^ {k - 1}.$$
$\mathsf{lo} _ {n - 1}(f ^ k)$ および $\mathsf{hi} _ {n - 1}(f ^ k)$ に共通して $\mathsf{lo} _ {n - 1}(f) ^ {k - 1}$ が現れていることに注意します。再帰 的に $\mathsf{lo} _ {n - 1}(f) ^ {k - 1}$ を計算することにすれば計算量 $T(n)$ は $T(n) = T(n - 1) + O(2 ^ n n ^ 2)$ より $T(n) = O(2 ^ n n ^ 2)$ です。
sqrt
$f\in R _ n$ に対して $f ^ {1/2}$ を $g ^ 2 = f$ を満たす $g\in R _ n$ と定義します (一般に $g$ は複数存在し得ます)。ただし $f(\emptyset) \neq 0$ かつ $(f(\emptyset) ) ^ {1/2}$ の存在を仮定します。
$f ^ {1/2}(\emptyset)$
$f ^ {1/2}(S+m)$ for all $S\subseteq \all{m}$
以上より $O(2 ^ n n ^ 2)$ 時間で計算できます。
これを一般化して $1/k$ 乗を計算することもできます。詳細な計算式は省略しますが、次が成り立ちます。
$f ^ {1/k} (\emptyset) = (f(\emptyset) ) ^ {1/k}$
$\mathsf{hi} _ m(f ^ {1/k}) = \dfrac{1}{k}(\mathsf{hi} _ m(f) \ast (\mathsf{lo} _ m(f ^ {1/k}) ^ {-1}) ^ {k - 1} ) = \dfrac{1}{k}(\mathsf{hi} _ m(f) \ast \mathsf{lo} _ m(f ^ {1/k}) \ast \mathsf{lo} _ m (f ^ {-1}) )$
つまり、$1/k$ 乗も $O(2 ^ n n ^ 2)$ 時間で計算できます。有理数 乗も $f ^ {p / q} = (f ^ {1/q}) ^ p$ により $O(2 ^ n n ^ 2)$ 時間です。
$f(\emptyset) = 0$ の場合に関して、私は ranked zeta transform -> 各点 sqrt -> ranked mobius transform とする計算方法しか知りません。ダブリングでも解けたら教えていただけると嬉しいです。
追記(2023/05/02): 多項式 との合成
https://codeforces.com/blog/entry/92183 の内容です。
多項式 $\displaystyle a(x) = \sum _ {i = 0} ^ {m - 1} a _ i x ^ i$ と集合冪級数 $f$ の合成として得られる集合冪級数 $\displaystyle g \coloneqq a\circ f = \sum _ {i = 0} ^ {m - 1} a _ i f ^ i$ を計算します。
pow と同様の方針で再帰 式を立てます。$S\subseteq\all{n - 1}$ に対して、次が成り立ちます。
$$\begin{aligned}
g(S + (n - 1)) &= \sum _ {i = 0} ^ {m - 1} a _ i \sum _ {X _ 1\sqcup X _ 2\sqcup \cdots \sqcup X _ i = S + (n - 1)} \prod _ {j = 1} ^ i f(X _ j) \\
&= \sum _ {i = 1} ^ {m - 1} ia _ i \sum _ {Y\subseteq S} f(Y + (n - 1)) \sum _ {X _ 1\sqcup X _ 2\sqcup \cdots \sqcup X _ {i - 1} = S \setminus Y} \prod _ {j = 1} ^ {i - 1} f(X _ j) \\
&= \sum _ {Y\subseteq S} f(Y + (n - 1)) \sum _ {i = 1} ^ {m - 1} ia _ i \cdot f ^ {i - 1} (S \setminus Y) \\
&= \sum _ {Y\subseteq S} f(Y + (n - 1)) \cdot (a' \circ f)(S\setminus Y) \\
&= (\mathsf{hi} _ {n - 1}(f) \ast \mathsf{lo} _ {n - 1}(a' \circ f))(S) \\
&= (\mathsf{hi} _ {n - 1}(f) \ast (a' \circ\mathsf{lo} _ {n - 1}(f) ) )(S).
\end{aligned}$$
つまり、$\mathsf{hi} _ {n - 1}(a\circ f) = \mathsf{hi} _ {n - 1}(f) \ast (a' \circ \mathsf{lo} _ {n - 1}(f))$ が成り立ちます。また $\mathsf{lo} _ {n - 1}(a\circ f) = a\circ \mathsf{lo} _ {n - 1}(f)$ です。
再帰 式に微分 $a'$ が現れています。再帰 を展開していくと、様々な整数組 $(i, j)$ に対して合成 $\displaystyle \left(\dfrac{\mathrm{d} ^ j a}{\mathrm{d} x ^ j}\right) \circ \mathsf{lo} _ {i}(f)$ が計算できれば十分なことが分かります。
そこで、次のように dp テーブルを定義します。求める値は $\mathsf{dp}(n, 0)$ です。
$$\mathsf{dp}(i, j) \coloneqq \left(\dfrac{\mathrm{d} ^ j a}{\mathrm{d} x ^ j}\right) \circ \mathsf{lo} _ {i}(f).$$
先ほどの再帰 式より、このテーブルは次のように $i$ の昇順に埋めることができます。
$$\begin{aligned}
\mathsf{dp}(0, j) &= \left(S \in \lbrace \emptyset \rbrace \mapsto \dfrac{\mathrm{d} ^ j a}{\mathrm{d} x ^ j}(f(S)) \in \mathbb{K}\right), \\
\mathsf{lo} _ {i - 1}(\mathsf{dp}(i, j)) &= \mathsf{dp}(i - 1, j), & \quad (1\leq i\leq n)\\
\mathsf{hi} _ {i - 1}(\mathsf{dp}(i, j)) &= \mathsf{hi} _ {i - 1}(f) \ast \mathsf{dp}(i - 1, j + 1). & \quad (1\leq i\leq n)\\
\end{aligned}$$
また、見るべき $i, j$ は $0\leq i \leq n,\ 0\leq j \leq n - i$ を満たす部分で十分なことが分かります。
計算量を評価します。$\mathsf{dp}(0,\ast)$ の計算は $O(nm)$ 時間です。$i = 1, 2, \ldots, n$ に対して、$\mathsf{dp}(i,\ast)$ の計算は大きさ $2 ^ {i - 1}$ の subset convolution を $n - (i - 1)$ 回行う部分がボトルネック となり $O(2 ^ i i ^ 2 (n - i))$ 時間です。
以上より、テーブル全体を埋めるのに掛かる時間計算量は $O(nm + 2 ^ n n ^ 2)$ です ($\displaystyle \sum _ {i = 1} ^ n 2 ^ i i ^ 2 (n - i) = O(2 ^ n n ^ 2)$ に注意)。
なお、テーブルのデータとして ranked zeta transform を行った状態のものを持っておくことで、ボトルネック となる subset convolution の定数倍を改善できます ( 実装例 )。
おわりに
ここで示した式通りに計算するプログラムを書いて簡単な動作確認はしているので正しいと思っていますが、万が一間違いがあれば教えてください。
