【競技プログラミング】形式的冪級数の応用テクニック(前編)

# 0.はじめに
形式的冪級数は、収束性を考えない無限に続く多項式（級数） $f(x)$ であり、競技プログラミングではこれを利用して数え上げなどの問題を解くテクニックの総称としても扱われます。

最近はyukicoderなどでの出題が増えており、atcoderで出されることがある重要なテクニックになりつつあります。

今回はそんな形式的冪級数の応用的なテクニックをまとめていきます。

以下$[x^k]f(x)$を$f(x)$の$x^k$の係数の意味で使います

わからないところや間違っていることがあれば[twitter](https://twitter.com/hotmanww/status/1227520715659694080)までお願いします

### 0.1. 前提知識
フーリエ変換
形式的冪級数の基礎

### 0.2.　前編で解説する内容
$1/f(x)$[verify](https://judge.yosupo.jp/problem/inv_of_formal_power_series)
$e^{f(x)}$[verify](https://judge.yosupo.jp/problem/exp_of_formal_power_series)
$\sqrt{f(x)}$[verify](https://judge.yosupo.jp/problem/sqrt_of_formal_power_series)
微分積分
$\log(f(x))$[verify](https://judge.yosupo.jp/problem/log_of_formal_power_series)
$f(x)^k$[verify](https://judge.yosupo.jp/problem/pow_of_formal_power_series)
商と余り
$f(x+c)$[verify](https://judge.yosupo.jp/problem/polynomial_taylor_shift)
$f(g(x))$[verify](https://judge.yosupo.jp/problem/composition_of_formal_power_series)
部分和問題[verify](https://judge.yosupo.jp/problem/sharp_p_subset_sum)

### 0.3. 定義
**形式的冪級数においてはf(x)について、f(c)や、g(x)の定数項が0でない時のf(g(x))が定義されないです(この制約のおかげで無限項の計算を有限項に丸めても問題がなくなります)**
また、四則演算は$\mathcal{O}(1)$ で行われると仮定します
$[x^k]f(x)$:$f(x)$の$x^k$の項
$1/f(x)$ :$g(x)\times f(x)=1$となるような$g(x)$
$e^{f(x)}:=\sum_{k=0}^{\infty}\frac{f(x)^k}{k!}$($[x^0]f(x)=0$の時のみ定義されます)
$\sqrt{f(x)}$:$g(x)\times g(x)=f(x)$となる$g(x)$(一般には定数項が1のf(x)以外には定義されないらしいのですが、今回はこのように定義します)
$\log(1+f(x)):=-\sum_{k=1}^{\infty}\frac{(-1)^k}{k}f(x)^k$($[x^0]f(x)=0$の時のみ定義されます)
$f(x+c):=\sum_{k=0}^{\infty}\frac{[x^k]f(x)}{k!}(x+c)^k$(形式的冪級数においてではなく多項式において定義されます)
$f(g(x)):=\sum_{k=0}^{\infty}\frac{[x^k]f(x)}{k!}g(x)^k$(形式的冪級数においては$[x^0]g(x)=0$を要請します(多項式においては要請されないです))
$f'(x):=\sum_{k=0}^{\infty}k([x^k]f(x)) x^{k-1}$
$\int f(x)dx:=\sum_{k=0}^{\infty}\frac{x^{k+1}}{k+1}$

によって定義されます

形式的冪級数における定義は基本的に$f(x)$に対する操作について操作後の$x^0 \sim x^N$項目までが、操作前の$f(x)$の$x^0 \sim x^N$項目までによって一意に定まるように定義されています

# 1.記号的ニュートン法
### 1.0.概要
形式的冪級数におけるニュートン法は
先頭$1$項を求める
先頭$2$項を求める
先頭$4$項を求める
$\vdots$
先頭$2^n$項を求める
を繰り返すことで
$O(1\log1+2\log2+4\log4+\dots+(N/2)\log(N/2)+N\log N)=O(N\log N)$で
$N$項目まで正確な値を求めるテクニックです

### 1.1.逆数の場合
$g_n(x)$ を以下のように定義します。

$$
\begin{align}
g_0(x) &\equiv \frac{1}{[x^0]f(x)} \pmod{x}\\\\
g_{n+1}(x) &\equiv g_{n}(x)(2-g_{n}(x)\times f(x)) \pmod{x^{2^{n+1}}}~~(n \geq 0)
\end{align}
$$

すると、$g_n(x) \equiv \frac{1}{f(x)} \pmod{x^{2^{n}}}$ となることが分かります。

- 証明

定義より $g_{n}(x)\times f(x) \equiv 1 \pmod{x^{2^n}}$ なので、
$1+h_{n}(x)x^{2^n} \equiv g_{n}(x)\times f(x) \pmod{x^{2^{n+1}}}$ と置くことができ

$$
\begin{align}
g_{n+1}(x)\times f(x) &\equiv (g_{n}(x)(2-g_{n}(x)\times f(x)))\times f(x)\\\\
&\equiv (g_{n}(x)\times f(x)) \times (2-g_{n}(x)\times f(x))\\\\
&\equiv (1+h(x)x^{2^n}) \times (2-(1+h(x)x^{2^n}))\\\\
&\equiv 1-h^2(x)x^{2^{n+1}} \\\\
&\equiv 1 \pmod{x^{2^{n+1}}}
\end{align}
$$
よって数学的帰納法により示せました

- 追記
分かりやすく補足されている記事が出ていました
[形式的冪級数の逆元を求める式はなんであんな形なんですか？ - Qiita](https://qiita.com/Suu0313/items/ee9995e24101baab1026)

### 1.2　色々な場合
同様の理屈によって

|求めたいもの|初項|漸化式|
|:---:|:---:|:---:|
|$1/f(x)$|$g_0(x)=1/([x^0]f(x))$|$g_{n+1}(x)=g_n(x)(2-g_n(x)*f(x)) \bmod x^{2^{n+1}}$|
|$e^{f(x)}$|$g_0(x)=1$|$g_{n+1}(x)=g_n(x)(f(x)+1-\log(g_n(x)))\bmod x^{2^{n+1}}$|
|$\sqrt{f(x)}$|$g_0(x)=\sqrt{[x^0]f(x)}$|$g_{n+1}(x)=(g_n(x)+f(x)/g_n(x))/2\bmod x^{2^{n+1}}$|

が成り立ちます
これは所謂普通のニュートン法と同様の式になっています
注意点として
$e^{f(x)}=\sum_{k=0}^{\infty}f(x)^k/k!$
より、$[x^0]f(x)$が$0$でなければ任意の自然数$N$に対し$[x^0]e^{f(x)}$が$f(x)$の先頭$N$項から決まらないので$e^{f(x)}$は定義されないです
また、$[x^0]f(x)=0$の時、この方法では$\sqrt{f(x)}$が計算出来ないので、$\sqrt{x^2f(x)}=x\sqrt{f(x)}$を使って定数項を0以外にしてください

$e^{f(x)}$の漸化式に出てくる$\log$の計算方法は後述します

# 2.微分積分
難しそうに見える両者ですが、実は簡単です。
微分は$\frac{d(\sum_{k=0}^{\infty}x^k)}{dx}=\sum_{k=0}^{\infty}k\times x^{k-1}$
積分は$\int (\sum_{k=0}^{\infty}x^k)dx=\sum_{k=0}^{\infty}\frac{x^{k+1}}{k+1}$
によって定義されるため、そのとおりに配列を計算するだけで、$O(N)$です

# 3.log
今までの知識でできます
$\log(f(x))=\int \frac{f'(x)}{f(x)}dx$
は形式的冪級数においても成り立つので、これを計算すれば$O(N\log N)$です

注意点としては
定義より、$[x^0]f(x)=1$でないといけません

# 4.(多項式としての)商/余り
$(1/f(x))\times g(x)$との違いは切り捨てる方向です
$(1/f(x))\times g(x)$は$x^N$項目より大きい項を切り捨てる割り算と考えられますが、
$g(x)$を$f(x)$で割った商は$x^0$項目より小さい項を切り捨てる割り算と考えられます
よって$x^{-1}$を代入して$(1/f(x^{-1})) x^{deg(f)-deg(g)}\times g(x^{-1})$を計算することによって商は手に入ります。
実装上は$x^{-k}$を扱えないため、適当に下駄を履かせます
余りは $g(x)-(商)\times f(x)$です

$f(x)\div g(x)=h(x)+(A_1x^{-1}+A_2x^{-2}+A_3x^{-3}+\cdots)$
と表せると考えて
$f(x)=h(x)g(x)+g(x)(A_1x^{-1}+A_2x^{-2}+A_3x^{-3}+\cdots)$
とした時に、
$g(x)(A_1x^{-1}+A_2x^{-2}+A_3x^{-3}+\cdots)$の$deg$は$deg(g)$未満なので余りの要件を満たすと考えると
$f(1/x)/g(1/x)=h(1/x)+(A_1  x+A_2 x^2+A_3 x^3+\cdots)$
となるため正当性がわかりやすいと思います

# 5.累乗
繰り返し二乗法により$O(N\log^2N)$ですが、$O(N\log N)$でもできます
$e^{k\times \log(f(x))}$
を使います
$\exp$と$\log$の適用条件に気をつけてください

# 6. f(x+c)
$f$ を$N$ 次多項式とします
$f(x)=\sum_{k=0}^{N}A_k x^k$とすると、
$f(x+c)$
$=\sum_{k=N}^{N}A_k(x+c)^k$
$=\sum_{k=0}^{N}\sum_{i=0}^{k}A_k\binom{k}{i}c^{k-i}x^{i}$
$=\sum_{k=0}^{N}\sum_{i=0}^{k}A_k k! \frac{c^{k-i}}{(k-i)!}\frac{x^i}{i!}$
$=\sum_{j=0}^{N}\sum_{i=0}^{N}A_{k}k!\frac{c^{j}}{j!} \frac{x^i}{i!}$ (但し、 $j=k-i$ とおいた)
ここで $(\sum_{k=0}^{N}A_k k! x^k) \times (\sum_{j=0}^{N} \frac{c^j}{j!} x^{-j}) = \sum_{k=0}^{N}\sum_{j=0}^{N}(A_{k}k!\frac{c^{j}}{j!})x^{k-j}$
より、これを計算し、最後に $\frac{1}{i!}$ をかけ合わせることで求める事ができます。
$x^{-j}$の項に関しては実装上は$x^N$などを掛けて下駄を履かせてあげるといいです

# 7. f(g(x))
### 追記
現在この操作を行う実用上も高速かつシンプルな $\mathcal{O}(N\log^2 N)$ の[アルゴリズム](https://arxiv.org/abs/2404.05177)が Yasunori Kinoshita(noshi91) らによって提案されているため今から実装する場合はそちらをおすすめします！
### 7.1.テイラー展開
テイラーの定理より多項式f,定数aに対して、
$f(x+a)=f(a)+f'(a)x+\frac{f''(a)}{2!}x^2+...$
であるが、実は$x$と$a$に形式的冪級数に置き換えることで、形式的冪級数に拡張することが出来(証明略)、
$g(x)=p(x)+q(x)x^k$として、$a$に$p(x)$、$x$に$q(x)x^k$を代入すると、$f(p(x)+q(x)x^k)=f(p(x))+f'(p(x))q(x)x^k+\frac{f''(p(x))}{2!}q^2(x)x^{2k}+...$
式の形より$\lceil \frac{N}{k}\rceil$項目まで計算すれば十分である事がわかります
ここでは$k=\sqrt{N/\log N}$とします

### 7.2.分割統治法(f(p(x))の計算)
$f(x)=s(x)+t(x)x^{\frac{|f(x)|}{2}}$とすると
$f(p(x))=s(p(x))+p(x)^{\frac{|f(x)|}{2}}t(p(x))$
となり、これを再帰的に計算することで$O((N\log N)^\frac{3}{2})$で計算できます

### 7.3. f(p(x))からf'(p(x))を計算する
二項目以降も同様に計算していたらオーダーが大きくなってしまうのですが、
そこで、$(f(p(x)))'=f'(p(x))p'(x)$を利用して$f'(p(x))=(1/p'(x))*(f(p(x)))'$
とすることで高速化ができます
$O(\sqrt{N\log N})$回$O(N\log N)$の操作を行うので計算量は$O((N\log N)^\frac{3}{2})$です
よってf(g(x))は$O((N\log N)^{3/2})$で計算できました

# 8. 部分和問題
集合Sに含まれる要素を1つずつ選んでMを作る方法の場合の数を$O(|S|+M\log M)$
で求められます
$[x^M]\prod_{e \in S}(1+x^e)$が求めたい数ですが、
これは$[x^M]e^{\sum_{e \in S}\log(1+x^e)}$です
ここで、logの和の計算が大変そうに見えますが、$\log(1+x^e)$は$x^{ke}$の項以外は$0$なので重複要素を同一視すれば調和級数で抑えられます
よって求められました。

# 9. おわりに
今回は計算に焦点を当てて解説しました
形式的冪級数はとても奥が深くて楽しいので是非、色々研究してみてほしいです。
後編はスターリング数などの数え上げに関する数や、多項式補完やmultipoint evalutionなどの形式的冪級数から少し離れた重要な操作、Subset Convolutionや、包除原理などバラエティに富んだ話題を書きたいと考えています

# 10. 参考資料
[形式的冪級数(Formal-Power-Series) | Luzhiled's memo](https://ei1333.github.io/luzhiled/snippets/math/formal-power-series.html)
僕の形式的冪級数はこれを~~パク~~参考にしています
[Operations on Formal Power Series - Codeforces](https://codeforces.com/blog/entry/56422)

# 11. 謝辞
記事の公開後、maspyさん/えびちゃんさんの指摘をもとに記事の怪しいところ等を加筆修正しました
ありがとうございます
# 12 ps
中編・後編の進捗が停滞している間に[maspyさんの神まとめ](https://maspypy.com/%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%83%BB%E5%BD%A2%E5%BC%8F%E7%9A%84%E3%81%B9%E3%81%8D%E7%B4%9A%E6%95%B0-%E9%AB%98%E9%80%9F%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%A7%E3%81%8D%E3%82%8B%E3%82%82%E3%81%AE)が出来ていました、さらに勉強したい人にはおすすめです