On implementing O(nlog(n)^2) algorithm of FPS composition

#	User	Rating
1	tourist	4009
2	jiangly	3823
3	Benq	3738
4	Radewoosh	3633
5	jqdai0815	3620
6	orzdevinwang	3529
7	ecnerwala	3446
8	Um_nik	3396
9	ksun48	3390
10	gamegame	3386

#	User	Contrib.
1	cry	166
2	maomao90	163
2	Um_nik	163
4	atcoder_official	161
5	adamant	159
6	-is-this-fft-	158
7	awoo	157
8	TheScrasse	154
9	nor	153
9	Dominater069	153

Hi codeforces! I want to show that how to implement the algorithm from noshi91 who showed that we could modify the Bostan--Mori's algorithm to make use in FPS composition and compositional inverse.

noshi91's blog: FPS の合成と逆関数、冪乗の係数列挙 Θ(n (log(n))^2) in Japanese.

UPD: my submission. relative code from line 874 to 923.

Compositional Inverse

The key idea is consider the bivariate formal power series

$$$ \dfrac{1}{1-yf(x)}\in\mathbb{C}\left\lbrack\left\lbrack x,y\right\rbrack\right\rbrack $$$

Let $$$f(0)=0$$$, we have

$$$ \begin{aligned} \left\lbrack x^k\right\rbrack\frac{1}{1-yf(x)}=\left\lbrack x^k\right\rbrack\sum_{i\geq 0}y^if(x)^i \end{aligned} $$$

Apply Bostan--Mori's algorithm here:

$$$ \begin{aligned} \frac{P(x,y)}{Q(x,y)}&=\frac{P(x,y)Q(-x,y)}{Q(x,y)Q(-x,y)} \\ &=\frac{U_{\text{even}}(x^2,y)+xU_{\text{odd}}(x^2,y)}{V(x^2,y)} \end{aligned} $$$

$$$ \left\lbrack x^k\right\rbrack\frac{P(x,y)}{Q(x,y)}= \begin{cases} \left\lbrack x^{k/2}\right\rbrack &\dfrac{U_{\text{even}}(x,y)}{V(x,y)},\space &\text{if }k\equiv 0\pmod{2} \\ \left\lbrack x^{(k-1)/2}\right\rbrack &\dfrac{U_{\text{odd}}(x,y)}{V(x,y)},\space &\text{if }k\equiv 1\pmod{2} \end{cases} $$$

The degree of $$$y$$$ doubled after each time of iteration, but the degree of $$$x$$$ halved.

$$$ \begin{array}{ll} &\textbf{Algorithm }\operatorname{Enum-}k\operatorname{th-Term-of-Power}(f,k,n)\text{:} \\ &\textbf{Input}\text{: } f\in x\mathbb{C}\left\lbrack x\right\rbrack,k,n\in\mathbb{Z}_{\geq 0}\text{.} \\ &\textbf{Output}\text{: }\left(\left\lbrack x^k\right\rbrack f^0,\left\lbrack x^k\right\rbrack f,\left\lbrack x^k\right\rbrack f^2,\dots ,\left\lbrack x^k\right\rbrack f^{n-1}\right)\text{.} \\ 1&P(x,y)\gets 1 \\ 2&Q(x,y)\gets 1-yf \\ 3&\textbf{while }\deg_y(P)<n-1\textbf{ or }k\neq 0\textbf{ do}\text{:} \\ 4&\qquad U_0(x^2,y)+xU_1(x^2,y)\gets P(x,y)Q(-x,y) \\ 5&\qquad V(x^2,y)\gets Q(x,y)Q(-x,y) \\ 6&\qquad P\gets U_{k\bmod{2}}(x,y) \bmod{x^{\left\lfloor k/2\right\rfloor +1}} \\ 7&\qquad Q\gets V(x,y) \bmod{x^{\left\lfloor k/2\right\rfloor +1}} \\ 8&\qquad k\gets \left\lfloor k/2\right\rfloor \\ 9&\textbf{end while} \\ 10&A(y)=\sum_{j=0}^{\deg_y(P)}a_jy^j\gets \left\lbrack x^0\right\rbrack P(x,y) \\ 11&\textbf{return }(a_0,\dots ,a_{n-1}) \end{array} $$$

Combined with Lagrange Inversion Formula, we have the following algorithm:

$$$ \begin{array}{ll} &\textbf{Algorithm }\operatorname{Compositional-Inverse}(f,n)\text{:} \\ &\textbf{Input}\text{: }f\in x\mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack, f'(0)\neq 0,n\in\mathbb{Z}_{\geq 2}\text{.} \\ &\textbf{Output}\text{: }g\pmod{x^n} \text{ such that }f(g)\equiv g(f)\equiv x\pmod{x^n}\text{.} \\ 1&t\gets f'(0) \\ 2&F(x)\gets f\left(t^{-1}x\right) \\ 3&(a_0,\dots ,a_{n-1})\gets \operatorname{Enum-}k\operatorname{th-Term-of-Power}(F,n-1,n) \\ &\because \left\lbrack x^{n-1}\right\rbrack F^k=\frac{k}{n-1}\left\lbrack x^{n-1-k}\right\rbrack \left(F^{\langle -1\rangle}/x\right)^{-(n-1)} \\ &\therefore \frac{n-1}{k}a_{k}=\left\lbrack x^{n-1-k}\right\rbrack\left(F^{\langle -1\rangle}/x\right)^{-(n-1)},\space (k>0) \\ 4&G(x)\gets \sum_{k=1}^{n-1}\frac{n-1}{k}a_{k}x^{n-1-k} \\ 5&F^{\langle -1\rangle}/x\gets \left(G(x)^{1/(n-1)}\right)^{-1}\mod{x^{n-1}} \\ &\because (F^{\langle -1\rangle}/x)(0)=1 \\ &\therefore F^{\langle -1\rangle}/x=\exp\left(\frac{1}{1-n}\log G\right) \\ &\because F^{\langle -1\rangle}\circ f\circ \left(t^{-1}x\right)=x \\ &\therefore \left((t^{-1}x) \circ F^{\langle -1\rangle}\right)\circ f=\left((t^{-1}x) \circ F^{\langle -1\rangle}\right)\circ f \circ \left((t^{-1}x)\circ (tx)\right)=x \\ 6&\textbf{return }\left((t^{-1}x) \circ F^{\langle -1\rangle}\right) \end{array} $$$

Composition

Given $$$f=\sum_{k\geq 0}f_kx^k\in\mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack,g\in x\mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack$$$ then

$$$ f(g)=\sum_{k\geq 0}f_kg^k $$$

Consider again

$$$ \dfrac{f(y^{-1})}{1-yg(x)}=\sum_{k\geq 0}(\cdots +f_ky^{-k}+\cdots)y^kg(x)^k $$$

Our goal is computing

$$$ \left\lbrack y^0\right\rbrack \dfrac{f(y^{-1})}{1-yg(x)} $$$

Apply Bostan--Mori's algorithm here:

$$$ \begin{aligned} \frac{P(y)}{Q(x,y)}&=\frac{P(y)}{Q(x,y)Q(-x,y)}Q(-x,y) \\ &=\frac{P(y)}{V(x^2,y)}Q(-x,y) \end{aligned} $$$

We want to compute $$$\dfrac{P(y)}{Q(x,y)Q(-x,y)}$$$ first, and we only need to track finitly many coefficients of $$$y^j$$$ for $$$j\leq 0$$$.

$$$ \begin{array}{ll} &\textbf{Algorithm }\operatorname{Composition-Subprocedure}(P,Q,n)\text{:} \\ &\textbf{Input}\text{: }P=\sum_{0\leq j\leq n}p_jy^{-j}\in\mathbb{C}((y)),Q\in\mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack\left\lbrack y\right\rbrack ,Q(0)=1,[x^0y^k]Q=0,(\forall k>0)\text{.} \\ &\textbf{Output}\text{: }\left\lbrack y^{\left\lbrack -\deg_y(Q)+1,0\right\rbrack}\right\rbrack\dfrac{P}{Q}\bmod{x^{n+1}}\text{.} \\ 1&d\gets \deg_y(Q)\\ 2&\textbf{if }n=0\textbf{ then return }\left(\left\lbrack y^{-d+1}\right\rbrack P,\dots ,\left\lbrack y^0\right\rbrack P\right) \\ &n=0\implies \deg_x(Q)=0 \\ &\because Q(0)=1,[x^0y^k]Q=0,(\forall k>0) \\ &\therefore Q\equiv 1\pmod{x^1},\space \left(\left\lbrack y^{-d+1}\right\rbrack P,\dots ,\left\lbrack y^0\right\rbrack P\right)=\left(\left\lbrack y^{-d+1}\right\rbrack P/Q,\dots ,\left\lbrack y^0\right\rbrack P/Q\right) \\ 3&V(x^2,y)\gets Q(x,y)Q(-x,y)\bmod{x^{n+1}} \\ 4&(t_{-2d+1},\dots ,t_0)\gets \operatorname{Composition-Subprocedure}\left(P,V(x,y),\left\lfloor n/2\right\rfloor\right) \\ 5&T(x,y)\gets \sum_{j=-2d+1}^0t_jy^j \\ 6&U(x,y)=\sum_{j=-2d+1}^d u_jy^j\gets T(x^2,y)Q(-x,y)\bmod{x^{n+1}} \\ 7&\textbf{return }\left(u_{-d+1},\dots ,u_0\right) \end{array} $$$

and the composition algorithm is just

$$$ \begin{array}{ll} &\textbf{Algorithm }\operatorname{Composition}(f,g,n)\text{:} \\ &\textbf{Input}\text{: }f,g\in\mathbb{C}\left\lbrack x\right\rbrack ,n\in\mathbb{Z}_{>0}\text{.} \\ &\textbf{Output}\text{: }f(g)\bmod{x^n}\text{.} \\ 1&c\gets g(0) \\ 2&F(x)\gets f(x+c)\bmod{x^n} \\ 3&G(x)\gets (g-c)\bmod{x^n} \\ 4&(u)\gets\operatorname{Composition-Subprocedure}\left(F\left(y^{-1}\right),1-yG,n-1\right) \\ &f(g)=f\circ (x+c)\circ G=F(G) \\ 5&\textbf{return }(u) \end{array} $$$

One can try with composition_of_formal_power_series_large at Library Checker.

Reference

Alin Bostan, Ryuhei Mori. A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence.
noshi91. FPS の合成と逆関数、冪乗の係数列挙 Θ(n (log(n))^2).

hly1204's blog

Compositional Inverse

Composition

Reference