Image Description
hotman78
  • 2020-08-14T18:08

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

\newcommand{\d}{\mathrm{d}}

はじめに

形式的冪級数は、無限に続く収束性を考えない関数 f(x)f(x) を考えて、数え上げなどの問題を解くテクニックです。

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

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

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

わからないところや間違っていることがあればtwitterまでお願いします

前提知識

フーリエ変換
形式的冪級数の基礎

前編で解説する内容

1/f(x)1/f(x) verify
ef(x)e^{f(x)} verify
f(x)\sqrt{f(x)} verify
微分積分
log⁡(f(x))\log(f(x)) verify
f(x)kf(x)^k verify
商と余り
f(x+c)f(x+c) verify
f(g(x))f(g(x)) verify
部分和問題 verify

定義

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

によって定義されます

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

記号的ニュートン法

概要

形式的冪級数におけるニュートン法は
先頭11項を求める
先頭22項を求める
先頭44項を求める
⋮\vdots
先頭2n2^n項を求める

を繰り返すことで

O(1log⁡1+2log⁡2+4log⁡4+⋯+(N/2)log⁡(N/2)+Nlog⁡N)=O(Nlog⁡N)\mathcal{O}(1\log 1+2\log 2+4\log 4+\dots+(N/2)\log(N/2)+N\log N)=\mathcal{O}(N\log N)

でNN項目まで正確な値を求めるテクニックです

逆数の場合

gn(x)g_n(x)を0∼2n−10\sim 2^n-1項目までの1/f(x)1/f(x)の値とします

gn+1(x)=gn(x)(2−gn(x)×f(x) mod x2n+1) mod x2n+1g_{n+1}(x)=g_{n}(x)(2-g_{n}(x)\times f(x) \bmod x^{2^{n+1}}) \bmod x^{2^{n+1}}
が成り立ちます

  • 証明

定義より gn(x)×f(x) mod x2n=1g_{n}(x)\times f(x) \bmod x^{2^n}=1 なので、
1+hn(x)x2n=gn(x)×f(x) mod x2n1+h_{n}(x)x^{2^n}=g_{n}(x)\times f(x) \bmod x^{2^n} と置くことができ

gn+1(x)×f(x) mod x2n+1=gn(x)×f(x)(2−gn(x)×f(x) mod x2n+1) mod x2n+1=(1+h(x)x2n)(2−(1+h(x)x2n)) mod x2n+1=1−h(x)2x2n+1 mod x2n+1=1\begin{align*} g_{n+1}(x)\times f(x)\bmod x^{2^{n+1}}&=g_{n}(x)\times f(x)(2-g_{n}(x)\times f(x) \bmod x^{2^{n+1}})\bmod x^{2^{n+1}}\\ &=(1+h(x)x^{2^n})(2-(1+h(x)x^{2^n})) \bmod x^{2^{n+1}}\\ &=1-h(x)^2x^{2^{n+1}} \bmod x^{2^{n+1}}\\ &=1 \end{align*}

よって数学的帰納法により示せました

色々な場合

同様の理屈によって

求めたいもの 初項 漸化式
1/f(x)1/f(x) g0(x)=1/([x0]f(x))g_0(x)=1/([x^0]f(x)) gn+1(x)=gn(x)(2−gn(x)∗f(x)) mod x2n+1g_{n+1}(x)=g_n(x)(2-g_n(x)*f(x)) \bmod x^{2^{n+1}}
ef(x)e^{f(x)} g0(x)=1g_0(x)=1 gn+1(x)=gn(x)(f(x)+1−log⁡(gn(x))) mod x2n+1g_{n+1}(x)=g_n(x)(f(x)+1-\log(g_n(x)))\bmod x^{2^{n+1}}
f(x)\sqrt{f(x)} g0(x)=[x0]f(x)g_0(x)=\sqrt{[x^0]f(x)} gn+1(x)=(gn(x)+f(x)/gn(x))/2 mod x2n+1g_{n+1}(x)=(g_n(x)+f(x)/g_n(x))/2\bmod x^{2^{n+1}}

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

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

微分積分

難しそうに見える両者ですが、実は簡単です。
微分は

ddx(∑k=0∞xk)=∑k=0∞k×xk−1\frac{\mathrm{d}}{\mathrm{d} x}(\sum_{k=0}^{\infty}x^k)=\sum_{k=0}^{\infty}k\times x^{k-1}

積分は

∫(∑k=0∞xk)dx=∑k=0∞xk+1k+1\int (\sum_{k=0}^{\infty}x^k)dx=\sum_{k=0}^{\infty}\frac{x^{k+1}}{k+1}

によって定義されるため、そのとおりに配列を計算するだけで、 O(N)O(N) です

log

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

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

(多項式としての)商/余り

(1/f(x))×g(x)(1/f(x))\times g(x) との違いは切り捨てる方向です
(1/f(x))×g(x)(1/f(x))\times g(x) は xNx^N 項目より大きい項を切り捨てる割り算と考えられますが、
g(x)g(x) を f(x)f(x) で割った商は x0x^0 項目より小さい項を切り捨てる割り算と考えられます
よって x−1x^{-1} を代入して (1/f(x−1))xdeg(f)−deg(g)×g(x−1)(1/f(x^{-1})) x^{deg(f)-deg(g)}\times g(x^{-1}) を計算することによって商は手に入ります。
実装上はx−kx^{-k}を扱えないため、適当に下駄を履かせます
余りは商を p(x)p(x) とすれば、 g(x)−p(x)×f(x)g(x)- p(x) \times f(x) です
f(x)÷g(x)=h(x)+(A1x−1+A2x−2+A3x−3+⋯ )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)(A1x−1+A2x−2+A3x−3+⋯ )f(x)=h(x)g(x)+g(x)(A_1x^{-1}+A_2x^{-2}+A_3x^{-3}+\cdots)
とした時に、
g(x)(A1x−1+A2x−2+A3x−3+⋯ )g(x)(A_1x^{-1}+A_2x^{-2}+A_3x^{-3}+\cdots) の deg\mathrm{deg}はdeg(g)\mathrm{deg}(g)未満なので余りの要件を満たすと考えると
f(1/x)/g(1/x)=h(1/x)+(A1x+A2x2+A3x3+⋯ )f(1/x)/g(1/x)=h(1/x)+(A_1 x+A_2 x^2+A_3 x^3+\cdots)
となるため正当性がわかりやすいと思います

累乗

繰り返し二乗法により O(Nlog⁡2N)\mathcal{O}(N\log^2N) ですが、 O(Nlog⁡N)\mathcal{O}(N\log N) でもできます

f(x)k=exp⁡(k×log⁡(f(x)))f(x)^k = \exp(k\times \log(f(x)))

を使います
exp⁡\exp と log⁡\log の適用条件に気をつけてください

f(x+c)

ff を NN 次多項式とします
f(x)=∑k=0NAkxkf(x)=\sum_{k=0}^{N}A_k x^k とすると、

f(x+c)=∑k=NNAk(x+c)k=∑k=0N∑i=0kAk(ki)ck−ixi=∑k=0N∑i=0kAkk!ck−i(k−i)!xii!=∑j=0N∑i=0NAkk!cjj!xii!\begin{align*} 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!} \end{align*}

(但し、 j=k−ij=k-i とおいた)
ここで

(∑k=0NAkk!xk)×(∑j=0Ncjj!x−j)=∑k=0N∑j=0N(Akk!cjj!)xk−j(\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}

より、これを計算し、最後に 1i!\frac{1}{i!} をかけ合わせることで求める事ができます。
x−jx^{-j} の項に関しては実装上は xNx^N などを掛けて下駄を履かせてあげるといいです

f(g(x))

テイラー展開

テイラーの定理より多項式 f,定数 a に対して、
f(x+a)=f(a)+f′(a)x+f′′(a)2!x2+⋯f(x+a)=f(a)+f'(a)x+\frac{f''(a)}{2!}x^2+\cdots
であるが、実は xx と aa に形式的冪級数を代入することで、形式的冪級数に拡張することが出来、
g(x)=p(x)+q(x)xkg(x)=p(x)+q(x)x^k として、 aa に p(x)p(x) 、 xx に q(x)xkq(x)x^k を代入すると、 f(p(x)+q(x)xk)=f(p(x))+f′(p(x))q(x)xk+f′′(p(x))2!q2(x)x2k+⋯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}+\cdots
式の形より ⌈Nk⌉\lceil \frac{N}{k}\rceil 項目まで計算すれば十分である事がわかります
ここでは k=N/log⁡Nk=\sqrt{N/\log N} とします

分割統治法(f(p(x))の計算)

f(x)=s(x)+t(x)x∣f(x)∣2f(x)=s(x)+t(x)x^{\frac{|f(x)|}{2}} とすると
f(p(x))=s(p(x))+p(x)∣f(x)∣2t(p(x))f(p(x))=s(p(x))+p(x)^{\frac{|f(x)|}{2}}t(p(x))
となり、これを再帰的に計算することで O((Nlog⁡N)32)\mathcal{O}((N\log N)^\frac{3}{2}) で計算できます

f(p(x))から f'(p(x))を計算する

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

部分和問題

集合 S に含まれる要素を 1 つずつ選んで M を作る方法の場合の数を O(∣S∣+Mlog⁡M)\mathcal{O}(|S|+M\log M) で求められます。(四則演算が O(1)\mathcal{O}(1) で出来るという仮定の元)
[xM]∏e∈S(1+xe)[x^M]\prod_{e \in S}(1+x^e) が求めたい数ですが、
これは [xM]e∑e∈Slog⁡(1+xe)[x^M]e^{\sum_{e \in S}\log(1+x^e)} です
ここで、 log⁡\log の和の計算が大変そうに見えますが、log⁡(1+xe)\log(1+x^e) は xkex^{ke} の項以外は 00 なので重複要素を同一視すれば調和級数で抑えられます
よって求められました。

おわりに

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

参考資料

形式的冪級数(Formal-Power-Series) | Luzhiled's memo
僕の形式的冪級数はこれをパク参考にしています
Operations on Formal Power Series - Codeforces

謝辞

記事の公開後、maspy さん/えびちゃんさんの指摘をもとに記事の怪しいところ等を加筆修正しました
ありがとうございます

追記

中編・後編の進捗が停滞している間にmaspy さんの神まとめが出来ていました、さらに勉強したい人にはおすすめです

;