FFT入门——学习笔记

FFT入门

给一个非常好的入门视频：

快速傅里叶变换

复数与单位根

定义：$i^2=-1$为虚数单位，我们称形如$a+bi(a,b\in R)$的数为复数。

我们可以用复数在复平面上表示点$(0,0)->(a,b)$的向量，我们称$x$的正轴与该向量的夹角为幅角，$\sqrt{a^2+b^2}$为模长。
下文中，我们默认$n,m$为2的幂。
我们以原点为起点，单位圆的$n$等分点为终点，设幅角为正且最小的向量所对应的复数为$w_n^0$，表示$n$次单位根，设其余$n-1$个单位根为$w_n^1,w_n^2\dots$。
这里显然有：$w_n^n=w_n^0=1$。

引出欧拉公式：

$$w_n^k=\cos \frac{2k\pi }{n}+i\sin \frac{2k\pi }{n}$$

那么由向量的运算法则（模长相乘，幅角相加），显然有$w_n^k=(w_n^1{)}^k$。

那么几个小性质：

1. $w_n^0=w_n^n=1$
2. $w_{dn}^{dk}=w_n^k$
3. $w_n^k=(w_n^1{)}^k$
4. $w_n^{k+\frac{n}{2}}=-w_n^k$

关于性质4的证明，可以将$k+\frac{n}{2}$看作旋转180度的结果，自然为复。

给张图

点值表示法

对于多项式$A(x)=a_0+a_{1}x+a_2{x}^2+\dots +a_{n-1}{x}^{n-1}$,$(a_0,a_1,a_2\dots )$为其系数表示法。

而带入$n$个不同的$x$，得到的$n$个二元组$(x_0,y_0),(x_1,y_1),(x_2,y_2)\dots \dots$即为多项式的点值表示法。

对于两个多项式$A(x),B(x)$，计算他们的乘积$C(x)=A(x)B(x)$的点值表示就只需要将$A(x),B(x)$对应的点值表示的$y$值相乘即可，也即$(x_0,y_{0,A}\times y_{0,B})\dots \dots$，但我们需要$n+m$个点。

可见点值表示法求多项式的积是非常方便的，这引出了一个极为伟大的思想：求出$A,B$的点值表示->求出$A(x)\times B(x)$的点值表示->将点值表示化为系数表示

其中第二步可以$O(n+m)$做到。其中表示$A,B$分别为$n,m$阶多项式。

因为$C$是$n+m$阶多项式，所以可以将$A,B$不足的位补0。下面我们来尝试优化第一个和第三个步骤，这就是FFT所做的事情。

FFT流程

我们现在讨论如何求出$A$的点值表达。将单位根$w_n^k$带入，得到：

$A(w_n^k)=\sum _{k^{\prime }=0}^{n-1}{a}_k{w}_n^{k{k}^{\prime }}$。
因为$n$是$2$的幂。考虑奇偶分组：

$A1(x)=a_0+a_2{x}^2+a_4{x}^4\dots ,A2(x)=a_{1}x+a_3{x}^3+a_5{x}^5\dots$。

得到：$A(x)=A1(x^2)+xA2(x^2)$。

将$w_n^k(k<\frac{n}{2})$带入得到：

$$A(w_n^k)=A1(w_n^{2k})+w_n^{k}A2(w_n^{2k})$$

将$w_n^{k+\frac{n}{2}}$带入得到：

$$A(w_n^{k+\frac{n}{2}})=A1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}A2(w_n^{2k+n})=A1(w_n^n{w}_n^{2k})-w_n^{k}A2(w_n^{2k})$$

$k,k+\frac{n}{2}$正好取遍$[0,n-1]$。
这个式子告诉我们，只需要我们处理左半区间的$A1,A2$的点值表达式，就可以快速得到$A$的整个区间的点值表达式。故可以分治处理，复杂度$O(n\log n)$。

可以写出代码：

//node 是复数类
db pi=acos(-1.0);
void solve(int limit,node a[],int tag){
    if(limit==1)return ;//单个常数没有用 
    node a1[(limit>>1)+5],a2[(limit>>1)+5];
    for(int i=0;i<=limit;i+=2){
    	a1[i>>1]=a[i];
    	a2[i>>1]=a[i+1];
	}//处理出系数 
	solve(limit>>1,a1,tag);
	solve(limit>>1,a2,tag);
	node wn={cos(2.0*pi/limit),tag*sin(2.0*pi/limit)},w={1.0,0};//pi是圆周率
	for(int i=0;i<(limit>>1);i++,w=w*wn){
		a[i]=a1[i]+w*a2[i];
		a[i+(limit>>1)]=a1[i]-w*a2[i];
	}
}
//tag=1

a数组就给出了$x=w_n^0,w_n^1,w_n^2\dots$的$A$的点值表示。
至于我们为什么要搞一个tag，会告诉你答案的。

点值化系数

直觉告诉我们，因为上文的solve在对其进行系数化点值的时候，我们搞了个变量tag=1，不难想到它的逆过程就是tag=-1，就可以让点值化系数了！

但是，这个系数不是我们想要的。

设上述FFT把$A(x)$的点值表示求出为$(y_0,y_1,y_2\dots )$，逆过程求出来的系数是$(c_0,c_1,c_2\dots )$，则$c_k=\sum _{j=0}^{n-1}{y}_j{w}_n^{-jk}$。

展开：

$$c_k=\sum _{j=0}^{n-1}\left(\sum _{i=0}^{n-1}{a}_i{w}_n^{ij}\right)w_n^{-jk}=\sum _{i=0}^{n-1}{a}_i\sum _{j=0}^{n-1}{w}_n^{j(i-k)}$$

将$w_n^{i-k}$看作整体，运用等比数列求和公式，得到：$c_k=\sum _{i=0}^{n-1}{a}_i\frac{w_n^{(i-k{)}^n}-1}{w_n^{i-k}-1}$

因为$w_n^{(i-k{)}^n}=(w_n^n{)}^{i-k}=1$，故分子为0，分母显然不为0，但当$i=k$时，得到其为：$n{a}_i$显然可以给出：$c_k=n{a}_k$，故可以给出：

$$a_k=\frac{c_k}{n}$$

这就是FFT求逆。点值化系数。

这有什么用呢？我们在求出来$A,B$的点值表示并且相乘后，再求逆得出$(c_0,c_1\dots )$，再对每一个$c$除以$n+m$就可以得出$C$的系数了！

所以整个过程是这样的。盗用一张图：

代码实现：

code

迭代换递归

不难看出，递归的解法虽然简洁明了，但却有极大的空间和常数开销，我们考虑将其从递归转化为迭代：

观察这一张图。

我们将底层的数的下标写作二进制数：

原序列	000	001	010	011	100	101	110	111
重排后	000	100	010	110	001	101	011	111

乍一看，我们好像发现：这个二进制数好像上下对应是反过来了。考虑证明这个性质，实际也不难，每一次划分就对一个二进制位进行了交换，自然会反过来。

将一个l位的二进制数$x$倒置，可以这样做：设r[x]为$x$倒置后的结果，则有r[x]=(r[x>>1]>>1)|((x&1)<<(l-1))，其中$limit=2^l$。

所以我们可以预处理出这个倒置的数组r，即可处理出合并的顺序，然后将长度为二的幂的区间自大到小合并即可，就省去了自顶向下带来的巨大常数和空间开销。

code

NTT

首先，类比FFT，我们利用了单位根的以下性质：

1. $w_n^0=w_n^n=1$
2. $w_{dn}^{dk}=w_n^k$
3. $w_n^k=(w_n^1{)}^k$
4. $w_n^{k+\frac{n}{2}}=-w_n^k$

而NTT是解决在模意义下的多项式乘法。

为什么我们需要NTT?因为FFT它炸精了！

而NTT的重要思想就是在整数域，模意义下的同样具有以上性质的整数，这让我们发现了——原根！

前置知识

阶

定义：
若$a,p$互素，那么满足$a^n\equiv 1(modp)$的最小正整数$n$即为$a$模$p$的阶，记作${\delta }_p(a)$

例如${\delta }_7(2)=3$

原根

设$p\in {\mathbb{N}}^{\mathbb{+}},a\in \mathbb{Z}$，若满足${\delta }_p(a)=\phi (p)$，则称$a$为模$p$的一个原根。

注意，对于模数$p$，如果它有原根，那么它的原根数量是$\phi (\phi (p))$。

对于$m$来说，存在模$m$的原根的重要条件是：$m=2,4,p^a,2p^a(p\in Prime,a\in {\mathbb{N}}^{\mathbb{+}})$

性质

对于原根，存在一个非常重要的定理：

设$p$是素数，$g$是模$p$的一个原根，那么$g^{i}modp(1<g<p,0\le i<p-1)$互不相同。

用原根代替单位根

这里因为$n$是2的幂，所以我们对$p$有一定要求，$p=a{2}^x+1$，常见的有：

$998244353=119\times 2^{23}+1,1004535809=479\times 2^{21}+1$。$3$是他们的原根之一。

设$g_n^n\equiv 1(modp)$且$g_n^1,g_n^1\dots \dots g_n^{n-1}$在模$p$下互不相同，则$g_n\equiv g^{\frac{p-1}{n}}(modp)$等价于$w_n^1$。

证明：

1. $g_n^n\equiv 1(modp)$

根据定义显然

2. $g_n^1,g_n^1\dots \dots g_n^{n-1}$在模$p$下互不相同

根据定义显然

3. $w_n^{k+\frac{n}{2}}=-w_n^k,w_n^2=w_{\frac{n}{2}}^1$

由于$g_n^1=g^{\frac{p-1}{n}}$，设$m=\frac{p-1}{n}$，则$nm=p-1$，当$n^{\prime }=\frac{n}{2}$时，$m^{\prime }=2m$。所以$g_n^2=g^{2m}=g_{\frac{n}{2}}^1$，这样我们就证明了后面一条定理
然后对于$g_n^{k+\frac{n}{2}}=g_n^k·g_n^{\frac{n}{2}}$，$g_n^{\frac{n}{2}}=g^{\frac{p-1}{2}}$，根据费马小定理，可以得到$g^{\frac{p-1}{2}}=1\text{或}-1$。然后因为$g^0\equiv 1(modp)$，根据性质2可以得到$g^{\frac{p-1}{2}}=-1$。所以带入即可得证。

4. $\sum _{j=0}^{n-1}{g}_n^{j(i-k)}$当且仅当在$i-k=0$时为$n$，否则为$0$

同理当$i-k=0$时显然为$n$
当$i\ne k$时，根据等比数列求和公式可以得到$\frac{g_n^{(i-k{)}^n}-1}{g_n^{i-k}-1}$，根据原根的定义和费马小定理：$g_n^n=g^{p-1}\equiv 1(modp)$，所以也可以得到分子为0。

综上，$g_n^1=g^{\frac{p-1}{n}}$为一个可替代$w$的选择。

在上面的FFT代码中，我们仅仅需要写个快速幂，再更改几行：

code

不过为什么我的NTT跑不过FFT啊，还慢得一批。
附上一张常用表：
//（$g$ 是$mod(r{2}^k+1)$的原根）

素数	r	k	g
3	1	1	2
5	1	2	2
17	1	4	3
97	3	5	5
193	3	6	5
257	1	8	3
7681	15	9	17
12289	3	12	11
40961	5	13	3
65537	1	16	3
786433	3	18	10
5767169	11	19	3
7340033	7	20	3
23068673	11	21	3
104857601	25	22	3
167772161	5	25	3
469762049	7	26	3
998244353	119	23	3
1004535809	479	21	3
2013265921	15	27	31
2281701377	17	27	3
3221225473	3	30	5
75161927681	35	31	3
77309411329	9	33	7
206158430209	3	36	22
2061584302081	15	37	7
2748779069441	5	39	3
6597069766657	3	41	5
39582418599937	9	42	5
79164837199873	9	43	5
263882790666241	15	44	7
1231453023109121	35	45	3
1337006139375617	19	46	3
3799912185593857	27	47	5
4222124650659841	15	48	19
7881299347898369	7	50	6
31525197391593473	7	52	3
180143985094819841	5	55	6
1945555039024054273	27	56	5
4179340454199820289	29	57	3
上面这一坨都不用管，只需要记下以下三个质数即可
$998244353,1004535809,469762049$，他们都存在原根$3$且都在int范围内。

任意模数NTT

它大概是这样的，由于算法竞赛中常见的模数是${10}^9$次方级别的，最常用的是998244353,1e9+7这两个，于是对于多项式乘法$A(x)B(x)$，设两个多项式分别是$n,m(n<m)$阶的，那么乘法所产生的最大值是${10}^9\times {10}^{9}m$，由于$m$不会大于${10}^6$级别，所以答案的值不会超过${10}^{24}$。可以这样考虑，我们选择3个1e9级别的模数，比如上文所说的$998244353,1004535809,469762049$，这三者的乘积显然是大于1e24的，我们分别对于这些数跑NTT，这很方便，但需要9次NTT。然后我们便得到了如下这个方程组：

$$\{\begin{array}{rl}x& \equiv a_1(mod{p}_1)\\ x& \equiv a_2(mod{p}_2)\\ x& \equiv a_3(mod{p}_3)\end{array}$$

因为直接合并三个质数会爆long long，于是我们需要科技。

根据CRT的方法，可以先合并前两个式子，这样$p_1{p}_2$不会爆long long，得到：

$$\{\begin{array}{rl}x& \equiv A(modP)\\ x& \equiv a_3(mod{p}_3)\end{array}$$

其中$P=p_1{p}_2,A=a_1{p}_2{p}_2^{-1}+a_2{p}_1{p}_1^{-1}$,$p_1^{-1},p_2^{-1}$分别是在模$p_2,p_1$时的逆元。

然后就有：$x=kP+A\equiv a_3(mod{p}_3)$，立即可推得：$k\equiv (a_3-A)\times P^{-1}(mod{p}_3)$。

然后就给出了答案：$x=kP+A(modp)$。

这里在做$P^{-1}$的时候long long会挂掉。而计算逆元用费马小定理，故我们需要龟速乘！奇淫技巧！
当然偷懒可以直接__int128。奇淫技巧类似于：

ll mul(ll a,ll b,ll P){
	a=(a%P+P)%P,b=(b%P+P)%P;
	return ((a*b-(ll)((long double)a/P*b+1e-6)*P)%P+P)%P;
}

建议背下。
于是稍加更改便可以写出三模数NTT代码：
细节蛮多的，我会尽量在代码中标注。

code

事实上：

简直慢的要死

MTT

同理解决任意模数的问题。

不妨将 $A(x),B(x)$ 分别的系数拆掉，也即设 $A1(x)\times 2^{15}+A2(x)=A(x),B(x)$ 同理。

然后我们来考虑加速，朴素的想法是：

$A(x)\times B(x)=A1(x)B1(x)\times 2^{30}+2^{15}\times (A2(x)\times B1(x)+A1(x)\times B2(x))+A2(x)B2(x)$

那么不妨利用 $FFT$ 的特性，令 $f_x=a{1}_x+i·a{2}_x,g_x=b{1}_x+i·b{2}_x$。

因此 $f\times g=A1(x)\times B1(x)-A2(x)\times B2(x)+i·(A1(x)\times B2(x)+A2(x)\times B1(x))$。

然后再有取 $f^{\prime }=a{1}_x+i·(-a{2}_x)$，做 $f^{\prime }\times g$。

得到 $A1(x)\times B1(x)+A2(x)\times B2(x)+i·(A1(x)\times B2(x)-A2(x)\times B1(x))$。

分列方程，即可解出所需四个乘法。

也即 $f+f^{\prime }$ 得到了 $a1b1,a1b2$，$f-f^{\prime }$ 得到了 $a2b2,-a2b1$。

一共五次FFT。非常优秀。

注意复数除以实数： $\frac{a+i·b}{c}=\frac{a}{c}+i·\frac{b}{c}$

code

分治FFT

我们一般的卷积式子是这样的：

$$a_i=\sum _{j=0}^i{b}_j{c}_{i-j}$$

这时候就是裸的FFT/NTT

但有时候会遇到如下毒瘤的式子：

$$a_i=\sum _{j=1}^i{a}_{i-j}{b}_i$$

我们发现，TA和正常的卷积并无什么差别，但问题是$a$出现在了式子的两边？？？

这时候，正常的FFT/NTT不能满足我们的要求，于是考虑改变这个算法，食用分治法：

假设我们已经知道了$[l,mid]$中的$a$值，那么左半区间$[l,mid]$对右半区间$[mid+1,r]$的贡献$F$可以计算为：

$$F_i=\sum _{j=l}^{mid}{a}_j{b}_{i-j}$$

在每轮分治求这个卷积的值的时候，将辅助数组$A$填充$a$的$[l,mid]$部分，$B$填充$b$的$[1,r-l]$部分。
这时候就又是一个多项式乘法的模版了，求出$F$之后累加到$a$上即可。复杂度$O(n{\log }^{2}n)$

code

多项式求逆

注意所有的多项式初等函数基本满足高位对低位无影响，所以可以直接先把 $n$ 倍增到 $2^k-1$，方便计算。

问题描述，给定一个多项式$F(x)$，求一个多项式$G(x)$，满足 $F(x)\ast G(x)\equiv 1(mod{x}^n)$（对$998244353$取模），也即 $F(x)\ast G(x)$只有常数项为1，其余项均为0。

在这里，我们设$F(x)=f_0+f_{1}x+f_2{x}^2\dots f_{n-1}{x}^{n-1},G(x)=g_0+g_{1}x+g_2{x}^2+g_3{x}^3\dots g_{n-1}{x}^{n-1}$，其中$f$是已知的。

那么由于仅有常数项为1，可以得到：$g_0=\frac{1}{f_0}$，于是我们可以求出$g_0$。

方法1—分治FFT

展开设$A(x)=F(x)\ast G(x)$，系数为$a_0,a_1\dots$，得到：$a_0=1,a_i=0(i\in [1,n-1])$。

于是：

$$a_i=\sum _{k=0}^i{f}_k{g}_{i-k}=0$$

移项变式得到：

$$f_0{g}_i=-\sum _{k=1}^i{f}_k{g}_{i-k}\text{可以给出}{g}_i=-\sum _{k=1}^i{g}_{i-k}\frac{f_k}{f_0}$$

令$f_k^{\prime }=-\frac{f_k}{f_0}(k>0)$，此题就变成了分治FFT，时间复杂度$O(n{\log }^{2}n)$。

代码的话仅仅需要上文分治FFT的主函数改为：

int main(){
//	freopen("data.in","r",stdin);
//	freopen("data.out","w",stdout);
	ios::sync_with_stdio(false);
	cin>>n;
	inv_g=power(g,p-2);
	for(int i=0;i<n;i++)cin>>b[i];
	b[0]%=p;
	a[0]=power(b[0],p-2);
	int INV=power(b[0],p-2);
	for(int i=1;i<n;i++)b[i]=-1ll*b[i]%p*INV%p,b[i]=(b[i]%p+p)%p; 
	cdq(0,n-1);
	for(int i=0;i<n;i++)cout<<a[i]<<" "; 
}

然后你去luogu上交了一发，发现：啊啊啊啊啊我怎么只有75啊啊啊啊啊

方法2—FFT+倍增

考虑倍增，设我们已经求解出来$F(x)\ast H(x)\equiv 1(mod{x}^{\lceil \frac{n}{2}\rceil })$，并由此推出$G(x)$。

$$F(x)\ast G(x)\equiv 1(mod{x}^n)⟹F(x)\ast G(x)\equiv 1(mod{x}^{\lceil \frac{n}{2}\rceil })$$

故二者做差可得：

$F(x)\ast (G(x)-H(x))\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil })$

所以可以得到：$G(x)-H(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil })⟹G(x{)}^2+H(x{)}^2-2G(x)H(x)\equiv 0(mod{x}^n)$。

同时乘上$F(x)$得到$F(x)G(x{)}^2+F(x)H(x{)}^2-2F(x)G(x)H(x)\equiv 0(mod{x}^n)$。

根据定义，得到：

$$G(x)+F(x)H(x{)}^2-2H(x)\equiv 0(mod{x}^n)⟹G(x)\equiv 2H(x)-F(x)H(x{)}^2(mod{x}^n)⟹G(x)\equiv H(x)(2-F(x)H(x))(mod{x}^n)$$

用FFT处理一下两个乘法即可。然后根据摊还分析可得复杂度为：$T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)$。

这里我们的所有推导都是基于常数项对于$x^n$有逆元，故多项式可逆的充要条件是其常数项是否可逆。也即，其常数项不为零。若在模意义下，则需要常数项存在逆元

code

这里其实我们可以发现，分治FFT和多项式求逆是可以相互转化的，比如上文的分治FFT的例题可以这样解决：

题目：求$a_i=\sum _{k=0}^{i-1}{a}_k{b}_{i-k},a_0=1$。

设数列$a,b$的生成函数分别是$A(x)=a_0+a_{1}x+a_2{x}^2\dots +a_n{x}^n+\dots ,B(x)=b_0+b_{1}x+b_2{x}^2+\dots b_n{x}^n+\dots$，设$b_0=0$。
设：$A(x)B(x)=H(x)$。

则有：

$$h_i=\sum _{k=0}^i{a}_k{b}_{i-k}⟹h_i=\sum _{k=0}^{i-1}{a}_k{b}_{i-k}+a_i{b}_0⟹h_i=a_i(i\ge 1)$$

而$h_0=a_0{b}_0=0$。

故可以得到：

$$A(x)B(x)=A(x)-1$$

移项化简可得：

$$A(x)=\frac{1}{1-B(x)}$$

而$1-B(x)$可以看作生成函数$C(x)=x_0·1+\sum _{i=1}^{i}0·x^i$减去$B(x)$，故$b_i^{\prime }=-b_i(i\ge 1),b_0^{\prime }=1$。

多项式求逆即可解决。

多项式 $ln$

多项式导数积分及 $ln$ 公式

必须有 $f_0=1$

根据函数求导法则可得：多项式 $F(x)=\sum _i{x}^i{f}_i$ 的求导法则为：

$F^{\prime }(x)=\sum i{f}_i{x}^{i-1}$，当然这必须要求 $f_0=0$

取 $ln$ 可以看作先求导再乱搞然后积分。

积分公式是 $F^{\prime }(x)=\sum \frac{x^i{f}_{i-1}}{i}$。

我们对 $\ln F(x)$ 先求导再积分。

也即 $(\ln F(x){)}^{\prime }=\frac{1}{F(x)}\ln F(x)$

	vector<int> derivation(vector<int>&f){
		vector<int>g(f.size()+1);
		for(int k=1;k<f.size();k++)g[k-1]=f[k]*k%p;
		return g;
	}
	vector<int> quadrature(vector<int>&f){
		vector<int>g(f.size()+1);
		for(int k=1;k<f.size();k++)g[k]=f[k-1]*inv[k]%p;g[0]=0;
		return g;
	}
	vector<int> getln(vector<int>f,int n){
		vector<int>f0=derivation(f);//求导
		f=getinv(f);//求逆
		for(int i=n;i<f.size();i++)f[i]=0;
		for(int i=n;i<f0.size();i++)f0[i]=0;
		int limit=1,d=0;
		while(limit<n+n)limit<<=1,++d;
		init(d);//不要忘了 
		ntt(f0,limit,1);ntt(f,limit,1);
		for(int i=0;i<limit;i++)f[i]=f[i]*f0[i]%p;
		ntt(f,limit,0);int inv=power(limit,p-2);
		for(int i=0;i<limit;i++)f[i]=f[i]*inv%p;
		return quadrature(f);//积分
	}

多项式平移

给定 $F(x)=\sum f_i{x}^i$，求 $F(x+n)$

$F(x+n)=\sum f_i(x+n{)}^i=\sum f_i\sum (\genfrac{}{}{0ex}{}{i}{j})x^j{n}^{i-j}=\sum _i{x}^i\sum (\genfrac{}{}{0ex}{}{k}{i})n^{k-i}{f}_k$

所以仅需计算 $\sum x^i\sum k!f_k\times \frac{n^{k-i}}{(k-i)!}$

后者是一个差卷积，直接计算即可。

利用该方法可以倍增计算第一类斯特林数行：$x(x+1)\dots \dots (x+n-1)=\sum x_{i}S1(n,i)$

多项式 $exp$

计算 $\exp F(x)$，必须有 $f_0=0$

先求导，得到：

$(\exp F(x){)}^{\prime }=F^{\prime }(x)\exp F(x)(mod{x}^n)$

而我们又知道，$f_i^{\prime }=(i+1)f_{i+1}$

因此有：

$$[x^i](\exp F(x){)}^{\prime }=\sum ([x^j]F^{\prime }(x))\times ([x^{i-j}]\exp F(x))$$

然后根据导数性质得到：

$$[x^{i+1}](i+1)(\exp F(x))=\sum (j+1)([x^{j+1}]F(x))\times ([x^{i-j}\exp F(x)])$$

不妨设 $\exp F(x)=\sum h_i{x}^i$

就可以得到：

$$h_{i+1}(i+1)=\sum (j+1)f_{j+1}{h}_{i-j}⟹h_i\times i=\sum j\times f_j{h}_{i-j}$$

整理一下也就是 $n·h_n=\sum j·f_j·h_{n-j}$

令 $w_j=j·f_j$，这就是个卷积了。

分治 FFT 即可。

改进的多项式 $exp$

假定 $n=2^k$

利用牛顿迭代法可以得到以下结论：

记 $F_0(x)=\exp G(x)(mod{x}^{2^k})$

则 $F(x)=\exp G(x)(mod{x}^{2^{k+1}})\equiv F_0(x)·(1+G(x)-\ln F_0(x))$

玄学错误：需要多倍增一次。

#include<bits/stdc++.h>
using namespace std;
#define N 1050500
#define int long long 
const int p=998244353,g=3,invg=(p+2)/3;
int n,m,r[N],invs[N];//数逆元
void init(int d){
    for(int i=0;i<(1<<d);i++)r[i]=(r[i>>1]>>1|((i&1)<<d-1));
}
int power(int a,int b){
    int res=1;
    while(b){
        if(b&1)res=res*a%p;
        a=a*a%p;b>>=1;
    }
    return res;
}
namespace NTT{
    void ntt(int f[],int d,int tag){
        for(int i=0;i<(1<<d);i++)if(i<r[i])swap(f[i],f[r[i]]);
        for(int len=1;len<(1<<d);len<<=1){
            int gn=power((tag?g:invg),(p-1)/(len<<1));
            for(int j=0;j<(1<<d);j+=(len<<1)){
                int g0=1;
                for(int k=0;k<len;k++,g0=g0*gn%p){
                    int x=f[j+k],y=f[j+k+len]*g0%p;
                    f[j+k]=(x+y)%p;f[j+k+len]=(x-y+p)%p;
                }
            }
        }
        if(!tag){
            int inv=power(1<<d,p-2);
            for(int i=0;i<(1<<d);i++)f[i]=f[i]*inv%p;
        }
    }
}
void Wf(int a[],int b[],int d){
    for(int i=0;i<(1<<d);i++)a[i]=b[i+1]*(i+1)%p;
}
void Jf(int a[],int b[],int d){
    for(int i=1;i<(1<<d);i++)a[i]=b[i-1]*invs[i]%p;a[0]=0;
}
namespace Inv{
    int A[N],B[N],C[N],S[N];
    void Inv(int a[],int b[],int d){
        S[0]=power(b[0],p-2);
        S[1]=0;
        for(int len=2,c=1;len<(1<<d);len<<=1,++c){
            int limit=len<<1;
            for(int i=0;i<len;i++)A[i]=b[i],B[i]=S[i];
            for(int i=len;i<limit;i++)A[i]=B[i]=0;
            init(c+1);
            NTT::ntt(B,c+1,1);NTT::ntt(A,c+1,1);
            for(int i=0;i<limit;i++)S[i]=2ll*B[i]%p+p-B[i]*B[i]%p*A[i]%p,S[i]%=p;
            NTT::ntt(S,c+1,0);
            for(int i=len;i<limit;i++)S[i]=0;
        }
        for(int i=0;i<(1<<d);i++)a[i]=S[i];
    }   
}
namespace Ln{
    int A[N];
    void Ln(int a[],int b[],int d){
        //LN(F(X))=F'(X)/F(X)
        for(int i=0;i<(1<<d);i++)A[i]=b[i];
        Wf(a,b,d);
        Inv::Inv(A,b,d);
        init(d);
        NTT::ntt(A,d,1);
        NTT::ntt(a,d,1);
        for(int i=0;i<(1<<d);i++)A[i]=a[i]*A[i]%p;
        NTT::ntt(A,d,0);
        // for(int i=(1<<d);i<(1<<d+1);i++)a[i]=0;
        Jf(a,A,d);
    }
}
int G[N],F[N];
namespace Exp{
    int A[N],B[N],C[N],S[N];
    void Exp(int g[],int f[],int d){//G(x)=exp(F(x)) mod x^d
    //G(x)=G0(X)*(1+f(x)-ln(G0(X)))
        for(int i=0;i<(1<<d);i++)S[i]=0;
        S[0]=1,S[1]=0;
        for(int len=2,c=1;len<(1<<d);len<<=1,++c){
            int limit=len<<1;
            for(int i=0;i<len;i++)A[i]=f[i],B[i]=S[i];
            for(int i=len;i<limit;i++)A[i]=B[i]=C[i]=0;
            Ln::Ln(C,B,c);
            for(int i=0;i<len;i++)C[i]=(p-C[i]+A[i])%p;
            C[0]+=1;
            init(c+1);
            NTT::ntt(B,c+1,1);
            NTT::ntt(C,c+1,1);
            for(int i=0;i<limit;i++)S[i]=B[i]*C[i]%p;
            NTT::ntt(S,c+1,0);
            for(int i=len;i<limit;i++)S[i]=0;
            // cout<<len<<"\n";
            // for(int i=0;i<len;i++)cout<<S[i]<<" ";cout<<"\n";
        }
        for(int i=0;i<(1<<d);i++)g[i]=S[i];
    }
}
signed main(){
    ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
    cin>>n;
    for(int i=0;i<n;i++)cin>>F[i];
    for(int i=0;i<(n<<2);i++)invs[i]=power(i,p-2);
    int l=0;
    while((1<<l)<=(n<<1))++l;
    Exp::Exp(G,F,l);
    for(int i=0;i<n;i++)cout<<G[i]<<" ";
}