[拉格朗日反演][FFT][NTT][多项式大全]详解

1、目录

• • • 1、目录
    • 2、多项式的两种表示法
    • 3、复数
      • 1.复数、复平面的定义
      • 2.欧拉公式
      • 3.复数的运算
      • 4.单位复数根
    • 4、DFT&FFT
      • 1.DFT
      • 2.FFT
    • 5、IDFT&IFFT
    • 6、NTT
      • 1.定义
      • 2.原根
      • 3.公式
    • 7、其他扩展
    • 8、代码实现
    • 9、拉格朗日反演

2、多项式的两种表示法

1.系数表示法

我们最常用的多项式表示法就是系数表示法，一个次数界为n$n$的多项式S(x)$S(x)$可以用一个向量s=(s0,s1,s2,⋯,sn−1)$s=(s_0,s_1,s_2,\cdots ,s_n-1)$系数表示如下：

S(x)=∑k=0n−1skxk$$S(x)=\sum _{k=0}^{n-1}{s}_k{x}^k$$

系数表示法很适合做加法，可以在O(n)$O(n)$的时间复杂度内完成，表达式为：

S(x)=A(x)+B(x)=∑k=0n−1(ak+bk)xk$$S(x)=A(x)+B(x)=\sum _{k=0}^{n-1}(a_k+b_k)x^k$$

当中

sk=ak+bk$$s_k=a_k+b_k$$

但是，系数表示法不适合做乘法，时间复杂度为O(n2)$O(n^2)$，表达式为：

S(x)=A(x)B(x)=∑k=0n−1(∑j=0n−1ajbk−j)xk$$S(x)=A(x)B(x)=\sum _{k=0}^{n-1}\left(\sum _{j=0}^{n-1}{a}_j{b}_{k-j}\right)x^k$$

当中

sk=∑j=0kajbk−j$$s_k=\sum _{j=0}^k{a}_j{b}_{k-j}$$

这就是卷积的一般形式，记s=a⊗b$s=a\otimes b$，我们要想办法加速这个过程。

2.点值表示法

顾名思义，点值就是多项式在一个点处的值。多项式A(x)$A(x)$的点值表达是一个集合：

{(x0,y0),(x1,y1),(x2,y2),⋯,(xn−1,yn−1)}$$\{(x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots ,(x_{n-1},y_{n-1})\}$$

使得对于k=0,1,2,⋯,n−1$k=0,1,2,\cdots ,n-1$有xk$x_k$两两不相同且yk=A(xk)$y_k=A(x_k)$。

n$n$个点可以确定唯一一个n$n$次多项式。

点值表达有很多优良的性质，加法和乘法都可以在O(n)$O(n)$的时间复杂度内完成。

现有A(x)$A(x)$的点值表达

{(x0,y0),(x1,y1),(x2,y2),⋯,(xn−1,yn−1)}$$\{(x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots ,(x_{n-1},y_{n-1})\}$$

和B(x)$B(x)$的点值表达

{(x0,y′0),(x1,y′1),(x2,y′2),⋯,(xn−1,y′n−1)}$$\{(x_0,y_0^{\prime }),(x_1,y_1^{\prime }),(x_2,y_2^{\prime }),\cdots ,(x_{n-1},y_{n-1}^{\prime })\}$$

则C(x)=A(x)+B(x)$C(x)=A(x)+B(x)$的点值表达为：

{(x0,y0+y′0),(x1,y1+y′1),(x2,y2+y′2),⋯,(xn−1,yn−1+y′n−1)}$$\{(x_0,y_0+y_0^{\prime }),(x_1,y_1+y_1^{\prime }),(x_2,y_2+y_2^{\prime }),\cdots ,(x_{n-1},y_{n-1}+y_{n-1}^{\prime })\}$$

C(x)=A(x)B(x)$C(x)=A(x)B(x)$的点值表达为：

{(x0,y0y′0),(x1,y1y′1),(x2,y2y′2),⋯,(xn−1,yn−1y′n−1)}$$\{(x_0,y_0{y}_0^{\prime }),(x_1,y_1{y}_1^{\prime }),(x_2,y_2{y}_2^{\prime }),\cdots ,(x_{n-1},y_{n-1}{y}_{n-1}^{\prime })\}$$

可见，点值表示可以帮助我们更快地进行卷积，可是如何在系数表示法和点值表示法之间相互转化呢？

3、复数

当x$x$为实数时，无法很好地对转换方法进行优化。为了优化计算xn$x^n$所浪费的时间，我们需要x$x$有循环的性质。但点值表示法需要n$n$个两两不同的值，而在实数域中只有1$1$和−1$-1$，因此，我们需要复数的帮助。

1.复数、复平面的定义

我们把形如a+bi$a+bi$的数称为复数z$z$，其中a$a$为实部(Real)，记为Rz$\mathrm{\Re }z$；b$b$为虚部(Imaginary)，记为Iz$\mathrm{\Im }z$。

每一点都对应唯一复数的平面叫复平面，相当于一个把Rz$\mathrm{\Re }z$作为横坐标，把Iz$\mathrm{\Im }z$作为纵坐标的笛卡尔坐标系。如图：

模长：复平面上原点到复数z$z$的距离，记为|z|$|z|$。根据勾股定理，|z|=|a+bi|=a2+b2−−−−−−√$|z|=|a+bi|=\sqrt{a^2+b^2}$

辐角：复平面上x$x$轴与复数z$z$所对应向量之间的夹角，在(−π2,π2)$(-\frac{\pi }{2},\frac{\pi }{2})$之间的记为辐角主值argz$\arg z$。

2.欧拉公式

大名鼎鼎的欧拉公式：

eit=cost+isint$$e^{it}=\cos t+i\sin t$$

根据三角函数在单位圆上的几何意义，公式是容易理解的。

几何意义：

当中θ$\theta$为角度，t$t$为弧长。

则根据欧拉公式，可将一个复数表示为一个二元组(a,θ)$(a,\theta )$，即模长和辐角（相当于复平面上极坐标系的表示方法）。值为:a(cosθ+isinθ)$a(\cos \theta +i\sin \theta )$

特殊情况：欧拉恒等式eiπ+1=0$e^{i\pi }+1=0$

3.复数的运算

（1）复数加法

运算规则：实部、虚部分别相加

(a+bi)+(c+di)=a+c+bi+di=(a+c)+(b+d)i$$(a+bi)+(c+di)=a+c+bi+di=(a+c)+(b+d)i$$

几何意义：如图

结果相当于两个向量所构成的平行四边形的对角线。如果把一个复数所对应的向量视为一个移动的变换，那么向量加法就是连续运用这两个变换相当于的新变换。

（2）复数乘法

运算规则：展开

(a+bi)(c+di)=ac+adi+bci+bdi2=(ac−bd)+(ad+bc)i$$(a+bi)(c+di)=ac+adi+bci+bd{i}^2=(ac-bd)+(ad+bc)i$$

几何意义：如图

如图，arga+argb=arga×b$\arg a+\arg b=\arg a\times b$,|a|×|b|=|a×b|$|a|\times |b|=|a\times b|$

总结就是：模长相乘，辐角相加。

因此，如果模长为1$1$，那么它的n$n$次方一定还在单位圆上。

证明：

根据欧拉公式，已知x=(a1,θ1)=a1(cosθ1+isinθ1),y=(a2,θ2)=a2(cosθ2+isinθ2)$x=(a_1,{\theta }_1)=a_1(\cos {\theta }_1+i\sin {\theta }_1),y=(a_2,{\theta }_2)=a_2(\cos {\theta }_2+i\sin {\theta }_2)$

则

x×y=a1a2(cosθ1+isinθ1)(cosθ2+isinθ2)=a1a2[(cosθ1cosθ2−sinθ1sinθ2)+i(cosθ1sinθ2+sinθ1cosθ2)]=a1a2[(cos(θ1+θ2)+cos(θ1−θ2)2+cos(θ1+θ2)−cos(θ1−θ2)2)+i(sin(θ1+θ2)−sin(θ1−θ2)2+sin(θ1+θ2)+sin(θ1−θ2)2)]=a1a2[cos(θ1+θ2)+isin(θ1+θ2)](积化和差公式)$$\begin{array}{rl}x\times y& =a_1{a}_2(\cos {\theta }_1+i\sin {\theta }_1)(\cos {\theta }_2+i\sin {\theta }_2)\\ & =a_1{a}_2\left[(\cos {\theta }_1\cos {\theta }_2-\sin {\theta }_1\sin {\theta }_2)+i(\cos {\theta }_1\sin {\theta }_2+\sin {\theta }_1\cos {\theta }_2)\right]\\ & =a_1{a}_2[\left(\frac{\cos ({\theta }_1+{\theta }_2)+\cos ({\theta }_1-{\theta }_2)}{2}+\frac{\cos ({\theta }_1+{\theta }_2)-\cos ({\theta }_1-{\theta }_2)}{2}\right)\\ \text{(积化和差公式)}& & +i\left(\frac{\sin ({\theta }_1+{\theta }_2)-\sin ({\theta }_1-{\theta }_2)}{2}+\frac{\sin ({\theta }_1+{\theta }_2)+\sin ({\theta }_1-{\theta }_2)}{2}\right)]& =a_1{a}_2\left[\cos ({\theta }_1+{\theta }_2)+i\sin ({\theta }_1+{\theta }_2)\right]\end{array}$$

∴|x×y|=|x|×|y|,arg(x×y)=argx+argy$\therefore |x\times y|=|x|\times |y|,\arg (x\times y)=\arg x+\arg y$

证毕。

4.单位复数根

（1）基本性质

单位复数根是方程ωn=1${\omega }^n=1$的解，第k$k$个解记为ωkn${\omega }_n^k$(这里的k$k$事实上是乘方的含义)

n=16$n=16$的解在复平面上的位置如下：

可以看到，n$n$个解把单位圆分成了n$n$等弧，交点即为根。而且，ωkn${\omega }_n^k$实际上是ωn${\omega }_n$的n$n$次方，模长仍为1$1$,辐角翻倍！

为什么呢？

∵|xn|=|x|n,argxn=nargx$\because |x^n|=|x{|}^n,\arg x^n=n\arg x$

∴|ω|n=|ωn|,argωn=nargω$\therefore |\omega {|}^n=|{\omega }^n|,\arg {\omega }^n=n\arg \omega$

∴|ω|n=1(|ω|∈R+),argω=360∘n$\therefore |\omega {|}^n=1(|\omega |\in {\mathbb{R}}^{+}),\arg \omega =\frac{{360}^{\circ }}{n}$

∴|ω|=1,argω=360∘n$\therefore |\omega |=1,\arg \omega =\frac{{360}^{\circ }}{n}$

这就很明显了。

所以，ωkn${\omega }_n^k$事实上表示的是ωn${\omega }_n$的k$k$次幂。为什么选择单位复数根呢？因为它有循环的优良性质，即ωnn=1${\omega }_n^n=1$。由于其他的都可以由ω1n${\omega }_n^1$得到，因此称为主n$n$次单位根，又记为ωn${\omega }_n$。

根据单位复数根的平分圆的意义和欧拉公式，ωkn=e2πikn=cos2πkn+isin2πkn${\omega }_n^k=e^{\frac{2\pi ik}{n}}=\cos \frac{2\pi k}{n}+i\sin \frac{2\pi k}{n}$。

（2）计算引理

显然，由于单位复数根循环(ωznn=e2πiz=[(eπi)2]z=1z=1${\omega }_n^{zn}=e^{2\pi iz}={\left[{\left(e^{\pi i}\right)}^2\right]}^z=1^z=1$)，有变换恒等式：

ωkn=ωk+wnn(w∈Z)$${\omega }_n^k={\omega }_n^{k+wn}(w\in \mathbb{Z})$$

每一份再分成k$k$份，编号也变成k$k$倍，位置自然不变(ωdkdn=e2πidkdn=e2πikn=ωkn${\omega }_{dn}^{dk}=e^{\frac{2\pi idk}{dn}}=e^{\frac{2\pi ik}{n}}={\omega }_n^k$)，所以有消去引理：

ωdkdn=ωkn$${\omega }_{dn}^{dk}={\omega }_n^k$$

由于过了n/2$n/2$就会绕过半圈(ωn/2n=eπinn=eπi=−1${\omega }_n^{n/2}=e^{\frac{\pi in}{n}}=e^{\pi i}=-1$)，所以有折半引理：

ωkn=−ωk±n/2n$${\omega }_n^k=-{\omega }_n^{k\pm n/2}$$

对单位复数根求和，根据几何级数(等差数列求和公式)，可得(∑k=0n−1(ωkn)j=(ωnn)j−1ω1n−1=0$\sum _{k=0}^{n-1}({\omega }_n^k{)}^j=\frac{({\omega }_n^n{)}^j-1}{{\omega }_n^1-1}=0$)，即有求和引理（要注意公式的使用条件）：

∑k=0n−1(ωkn)j=0,ωkn≠1$$\sum _{k=0}^{n-1}({\omega }_n^k{)}^j=0,{\omega }_n^k\ne 1$$

4、DFT&FFT

1.DFT

DFT就是求多项式A(x)$A(x)$在点(ω0n,ω1n,ω2n,⋯,ωn−1n)$({\omega }_n^0,{\omega }_n^1,{\omega }_n^2,\cdots ,{\omega }_n^{n-1})$处取值的过程。即：

yk=A(ωkn)=∑j=0n−1ajωkjn$$y_k=A({\omega }_n^k)=\sum _{j=0}^{n-1}{a}_j{\omega }_n^{kj}$$

结果y=(y0,y1,y2,⋯,yn−1)$y=(y_0,y_1,y_2,\cdots ,y_{n-1})$就是a$a$的离散傅里叶变换(DFT),记为y=DFTn(a)$y=DF{T}_n(a)$

2.FFT

（1）递归

DFT的O(n2)$O(n^2)$算法太慢了，FFT使用分治策略优化速度到O(nlogn)$O(n\log n)$。

考虑奇偶分治。

现在，我们假设n=2t$n=2^t$，设原系数a=(a0,a1,a2,⋯,an−1)$a=(a_0,a_1,a_2,\cdots ,a_{n-1})$，分治为偶数部分a1=(a0,a2,a4,⋯,an−2)$a_1=(a_0,a_2,a_4,\cdots ,a_{n-2})$，奇数部分a2=(a1,a3,a5,⋯,an−1)$a_2=(a_1,a_3,a_5,\cdots ,a_{n-1})$，已经递归求出y1=DFTn/2(a1)$y_1=DF{T}_{n/2}(a_1)$，y2=DFTn/2(a2)$y_2=DF{T}_{n/2}(a_2)$，现在我们要合并y1,y2$y_1,y_2$，得到y=DFTn(a)$y=DF{T}_n(a)$(蝴蝶操作)。

对于n=1$n=1$的边界情况，结果是显然的：因为k=0$k=0$，故ω01=1${\omega }_1^0=1$，即结果等于原系数。

对于n>1$n>1$，现在我们枚举k∈[1,n]$k\in [1,n]$要合并出yk$y_k$：

yk=A(ωkn)=∑j=0n−1ajωkjn=∑j=0n/2−1a2jω2kjn+∑j=0n/2−1a2j+1ω2kj+kn=∑j=0n/2−1a2jω2kjn+ωkn∑j=0n/2−1a2j+1ω2kjn=∑j=0n/2−1a1jωkjn/2+ωkn∑j=0n/2−1a2jωkjn/2(消去引理)$$\begin{array}{rl}{y}_k& =A({\omega }_n^k)\\ & =\sum _{j=0}^{n-1}{a}_j{\omega }_n^{kj}\\ & =\sum _{j=0}^{n/2-1}{a}_{2j}{\omega }_n^{2kj}+\sum _{j=0}^{n/2-1}{a}_{2j+1}{\omega }_n^{2kj+k}\\ & =\sum _{j=0}^{n/2-1}{a}_{2j}{\omega }_n^{2kj}+{\omega }_n^k\sum _{j=0}^{n/2-1}{a}_{2j+1}{\omega }_n^{2kj}\\ \text{(消去引理)}& & =\sum _{j=0}^{n/2-1}{a}_{1_j}{\omega }_{n/2}^{kj}+{\omega }_n^k\sum _{j=0}^{n/2-1}{a}_{2_j}{\omega }_{n/2}^{kj}\end{array}$$

• 对于k<n/2$k<n/2$：
  

  yk=y1k+ωkny2k$$\begin{array}{rl}{y}_k& =y_{1_k}+{\omega }_n^k{y}_{2_k}\end{array}$$

• 对于k≥n/2$k\ge n/2$：
  

  yk=∑j=0n/2−1a1jω(k−n/2)jn/2+ωkn∑j=0n/2−1a2jω(k−n/2)jn/2=y1k−n/2+ωkny2k−n/2=y1k−n/2−ωk−n/2ny2k−n/2(变换恒等式)(折半引理)$$\begin{array}{rl}\text{(变换恒等式)}& y_k& =\sum _{j=0}^{n/2-1}{a}_{1_j}{\omega }_{n/2}^{(k-n/2)j}+{\omega }_n^k\sum _{j=0}^{n/2-1}{a}_{2_j}{\omega }_{n/2}^{(k-n/2)j}& =y_{1_{k-n/2}}+{\omega }_n^k{y}_{2_{k-n/2}}\\ \text{(折半引理)}& & =y_{1_{k-n/2}}-{\omega }_n^{k-n/2}{y}_{2_{k-n/2}}\end{array}$$

我们用k+n/2$k+n/2$替代k$k$，就得到yk+n/2=y1k−ωkny2k$y_{k+n/2}=y_{1_k}-{\omega }_n^k{y}_{2_k}$

结合在一起就得到{ykyk+n/2=y1k+ωkny2k=y1k−ωkny2k$\{\begin{array}{ll}{y}_k& =y_{1_k}+{\omega }_n^k{y}_{2_k}\\ y_{k+n/2}& =y_{1_k}-{\omega }_n^k{y}_{2_k}\end{array}$
这样我们就可以把两个n/2$n/2$长的序列合并为一个n$n$长的序列了。

以下图的递归序，就可以在O(nlogn)$O(n\log n)$的时间复杂度内完成求解了。

（2）迭代

递归方法消耗内存、时间过大，无法承受。我们每次把下标分为奇数部分和偶数部分，是否有办法直接求出最后的递归运算顺序，以避免递归呢？

这样想：
第一次奇偶划分，我们按照二进制的倒数第一位排序；
第二次奇偶划分，我们按照二进制的倒数第二位排序；
第n$n$次奇偶划分，我们按照二进制的倒数第n$n$位排序；
依此类推。

因此，结果顺序就是原序列按照二进制位翻转的大小排序的结果。只要依次交换ak,arev(k)$a_k,a_{rev(k)}$，求出序列，就可以用迭代方法相邻归并实现快速傅里叶变换。

或者，我们也可以用更加代数的方法来发现这个结论。
已知现在位置为k=(b1b2b3⋯bn)2$k=(b_1{b}_2{b}_3\cdots b_n{)}_2$，按照奇偶重排：

• 若bn=0$b_n=0$，则位置变为k2=(0b1b2⋯bn−1)2=(bnb1b2⋯bn−1)2$\frac{k}{2}=(0b_1{b}_2\cdots b_{n-1}{)}_2=(b_n{b}_1{b}_2\cdots b_{n-1}{)}_2$，即为把最后一位提到第一位。
• 若bn=1$b_n=1$，则位置变为2n−1−1+k+12=2n+k−12=(1b1b2⋯bn−10)22=(bnb1b2⋯bn−1)2$2^{n-1}-1+\frac{k+1}{2}=\frac{2^n+k-1}{2}=\frac{(1b_1{b}_2\cdots b_{n-1}0{)}_2}{2}=(b_n{b}_1{b}_2\cdots b_{n-1}{)}_2$，同样是把最后一位提到第一位。

则反复n$n$次之后，就相当于二进制反转了。

如n=8$n=8$时，求出二进制：

0	1	2	3	4	5	6	7
0002${000}_2$	0012${001}_2$	0102${010}_2$	0112${011}_2$	1002${100}_2$	1012${101}_2$	1102${110}_2$	1112${111}_2$

翻转：

0	1	2	3	4	5	6	7
0002${000}_2$	1002${100}_2$	0102${010}_2$	1102${110}_2$	0012${001}_2$	1012${101}_2$	0112${011}_2$	1112${111}_2$

按翻转后的值排序：

0	4	2	6	1	5	3	7
0002${000}_2$	0012${001}_2$	0102${010}_2$	0112${011}_2$	1002${100}_2$	1012${101}_2$	1102${110}_2$	1112${111}_2$

这样就可以把奇偶合并转化为左右归并，迭代实现了。

5、IDFT&IFFT

何把点值表达变回系数表达呢？如果把求值写成矩阵形式，就是：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢ω0nω0nω0nω0n⋮ω0nω0nω1nω2nω3n⋮ωn−1nω0nω2nω4nω6n⋮ω2(n−1)nω0nω3nω6nω9n⋮ω3(n−1)n⋯⋯⋯⋯⋱⋯ω0nωn−1nω2(n−1)nω3(n−1)n⋮ω(n−1)2n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢a0a1a2a3⋮an−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢y0y1y2y3⋮yn−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{cccccc}{\omega }_n^0& {\omega }_n^0& {\omega }_n^0& {\omega }_n^0& \cdots & {\omega }_n^0\\ {\omega }_n^0& {\omega }_n^1& {\omega }_n^2& {\omega }_n^3& \cdots & {\omega }_n^{n-1}\\ {\omega }_n^0& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \cdots & {\omega }_n^{2(n-1)}\\ {\omega }_n^0& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \cdots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ {\omega }_n^0& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ ⋮\\ y_{n-1}\end{array}\right]$$

如果我们要把y$y$变成a$a$，就需要求出第一个矩阵的逆，即：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢a0a1a2a3⋮an−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢ω0nω0nω0nω0n⋮ω0nω0nω1nω2nω3n⋮ωn−1nω0nω2nω4nω6n⋮ω2(n−1)nω0nω3nω6nω9n⋮ω3(n−1)n⋯⋯⋯⋯⋱⋯ω0nωn−1nω2(n−1)nω3(n−1)n⋮ω(n−1)2n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥−1⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢y0y1y2y3⋮yn−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]={\left[\begin{array}{cccccc}{\omega }_n^0& {\omega }_n^0& {\omega }_n^0& {\omega }_n^0& \cdots & {\omega }_n^0\\ {\omega }_n^0& {\omega }_n^1& {\omega }_n^2& {\omega }_n^3& \cdots & {\omega }_n^{n-1}\\ {\omega }_n^0& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \cdots & {\omega }_n^{2(n-1)}\\ {\omega }_n^0& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \cdots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ {\omega }_n^0& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right]}^{-1}\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ ⋮\\ y_{n-1}\end{array}\right]$$

这个范德蒙德矩阵极为特殊，它的逆矩阵是：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢a0a1a2a3⋮an−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥=1n⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢ω0nω0nω0nω0n⋮ω0nω0nω−1nω−2nω−3n⋮ω−(n−1)nω0nω−2nω−4nω−6n⋮ω−2(n−1)nω0nω−3nω−6nω−9n⋮ω−3(n−1)n⋯⋯⋯⋯⋱⋯ω0nω−(n−1)nω−2(n−1)nω−3(n−1)n⋮ω−(n−1)2n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢y0y1y2y3⋮yn−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]=\frac{1}{n}\left[\begin{array}{cccccc}{\omega }_n^0& {\omega }_n^0& {\omega }_n^0& {\omega }_n^0& \cdots & {\omega }_n^0\\ {\omega }_n^0& {\omega }_n^{-1}& {\omega }_n^{-2}& {\omega }_n^{-3}& \cdots & {\omega }_n^{-(n-1)}\\ {\omega }_n^0& {\omega }_n^{-2}& {\omega }_n^{-4}& {\omega }_n^{-6}& \cdots & {\omega }_n^{-2(n-1)}\\ {\omega }_n^0& {\omega }_n^{-3}& {\omega }_n^{-6}& {\omega }_n^{-9}& \cdots & {\omega }_n^{-3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ {\omega }_n^0& {\omega }_n^{-(n-1)}& {\omega }_n^{-2(n-1)}& {\omega }_n^{-3(n-1)}& \cdots & {\omega }_n^{-(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ ⋮\\ y_{n-1}\end{array}\right]$$

只是把每项取倒数，并将结果除以n$n$即可。

证明：

记原矩阵为Ann$A_{nn}$,我们给出的矩阵为Bnn$B_{nn}$。

∴Aij=ωijn,Bij=1nω−ijn(0≤i,j≤n−1)$\therefore A_{ij}={\omega }_n^{ij},B_{ij}=\frac{1}{n}{\omega }_n^{-ij}(0\le i,j\le n-1)$

ABij=∑k=0n−1AikBkj=1n∑k=0n−1ωiknω−kjn=1n∑k=0n−1ωik−kjn=1n∑k=0n−1(ωi−jn)k$$\begin{array}{rl}A{B}_{ij}& =\sum _{k=0}^{n-1}{A}_{ik}{B}_{kj}\\ & =\frac{1}{n}\sum _{k=0}^{n-1}{\omega }_n^{ik}{\omega }_n^{-kj}\\ & =\frac{1}{n}\sum _{k=0}^{n-1}{\omega }_n^{ik-kj}\\ & =\frac{1}{n}\sum _{k=0}^{n-1}({\omega }_n^{i-j}{)}^k\end{array}$$

• 当i=j$i=j$时：

ABij=1n∑k=0n−11k=1n×n=1$$\begin{array}{rl}A{B}_{ij}& =\frac{1}{n}\sum _{k=0}^{n-1}{1}^k\\ & =\frac{1}{n}\times n\\ & =1\end{array}$$

• 当i≠j$i\ne j$时：

ABij=1n∑k=0n−1(ωi−jn)k=0(求和引理)$$\begin{array}{rl}A{B}_{ij}& =\frac{1}{n}\sum _{k=0}^{n-1}({\omega }_n^{i-j}{)}^k\\ \text{(求和引理)}& & =0\end{array}$$

综上所述，AB=In$AB=I_n$，即B=A−1$B=A^{-1}$

以上，离散傅里叶逆变换(IDFT)的表达式为：

ak=1n∑j=0n−1ykω−kjn$$a_k=\frac{1}{n}\sum _{j=0}^{n-1}{y}_k{\omega }_n^{-kj}$$

记a=IDFTn(y)$a=IDF{T}_n(y)$。

同理，可以用相同的方法把IDFT加速到O(nlogn)$O(n\log n)$,称为IFFT。

6、NTT

1.定义

有时候我们会想要模素数p$p$意义下的多项式乘法。此时，由次数界为n$n$的多项式A(x),B(x)$A(x),B(x)$的系数表达a,b$a,b$求S(x)=A(x)B(x)$S(x)=A(x)B(x)$的系数表达s$s$的公式为：

sk=∑j=0kajbk−jmodp$$s_k=\sum _{j=0}^k{a}_j{b}_{k-j}\mathrm{mod}p$$

FFT无能为力，我们需要一种新的DFT，以数论的办法进行，这就是快速数论变换(NTT)。

2.原根

（1）定义

我们需要一种有数论循环性质的新数，原根恰好满足我们的要求。

设m$m$为正整数，a$a$为整数，若a$a$模m$m$的阶等于φ(m)$\phi (m)$，则称a$a$为模m$m$的一个原根。

假设g$g$是素数p$p$的原根，有1<g<p$1<g<p$，且对于k=0,1,2,⋯,p−1$k=0,1,2,\cdots ,p-1$，有gkmodp$g^k\mathrm{mod}p$的结果两两不同，且gp−1≡1(modp)$g^{p-1}\equiv 1(\mathrm{mod}p)$。

可以发现，原根同样有循环的性质。因此，我们类比ωkn${\omega }_n^k$的定义，把原来的ωkn=e2πikn${\omega }_n^k=e^{\frac{2\pi ik}{n}}$替换为gk(p−1)n$g^{\frac{k(p-1)}{n}}$。

（2）性质

我们来证明一些类似单位复数根的性质。

变换恒等式：
因为：

gp−1≡1$$g^{p-1}\equiv 1$$

所以：

g(k+n)(p−1)n≡gk(p−1)ngp−1≡gk(p−1)n$$g^{\frac{(k+n)(p-1)}{n}}\equiv g^{\frac{k(p-1)}{n}}{g}^{p-1}\equiv g^{\frac{k(p-1)}{n}}$$

---

消去引理：
显然有：

gdk(p−1)dn≡gk(p−1)n$$g^{\frac{dk(p-1)}{dn}}\equiv g^{\frac{k(p-1)}{n}}$$

---

折半引理：
因为：

g(k+n2)(p−1)n≡gk(p−1)ngp−12$$g^{\frac{(k+\frac{n}{2})(p-1)}{n}}\equiv g^{\frac{k(p-1)}{n}}{g}^{\frac{p-1}{2}}$$

所以若要使：

gk(p−1)n+g(k+n2)(p−1)n≡0(modp)$$g^{\frac{k(p-1)}{n}}+g^{\frac{(k+\frac{n}{2})(p-1)}{n}}\equiv 0(\mathrm{mod}p)$$

成立，即：

gk(p−1)2(1+gp−12)≡0$$g^{\frac{k(p-1)}{2}}\left(1+g^{\frac{p-1}{2}}\right)\equiv 0$$

只需证：

gp−12≡p−1$$g^{\frac{p-1}{2}}\equiv p-1$$

由于：

gk≡0,1,2,⋯,p−1$$g^k\equiv 0,1,2,\cdots ,p-1$$

那么我们可以设：

gp−12≡x,x=0,1,2,⋯,p−1$$g^{\frac{p-1}{2}}\equiv x,x=0,1,2,\cdots ,p-1$$

因为：

(gp−12)2≡gp−1≡1$${\left(g^{\frac{p-1}{2}}\right)}^2\equiv g^{p-1}\equiv 1$$

所以：

x2−1≡0$$x^2-1\equiv 0$$

即：

(x+1)(x−1)≡0$$(x+1)(x-1)\equiv 0$$

又因为p$p$为素数，所以有：

x+1≡0或x−1≡0$$x+1\equiv 0或x-1\equiv 0$$

所以：

x=1,p−1$$x=1,p-1$$

又因为：

gp−1≡1,gkmodp两两不同$$g^{p-1}\equiv 1,g^k\mathrm{mod}p两两不同$$

所以：

x=p−1$$x=p-1$$

即：

gp−12≡p−1$$g^{\frac{p-1}{2}}\equiv p-1$$

得证：

gk(p−1)n+g(k+n2)(p−1)n≡0$$g^{\frac{k(p-1)}{n}}+g^{\frac{(k+\frac{n}{2})(p-1)}{n}}\equiv 0$$

---

求和引理：

∑k=0n−1gkj(p−1)n≡∑k=0n−1(gj(p−1)n)k≡gnj(p−1)n−1gj(p−1)n−1≡(gp−1)j−1gj(p−1)n−1≡0$$\sum _{k=0}^{n-1}{g}^{\frac{kj(p-1)}{n}}\equiv \sum _{k=0}^{n-1}{\left(g^{\frac{j(p-1)}{n}}\right)}^k\equiv \frac{g^{\frac{nj(p-1)}{n}}-1}{g^{\frac{j(p-1)}{n}}-1}\equiv \frac{{\left(g^{p-1}\right)}^j-1}{g^{\frac{j(p-1)}{n}}-1}\equiv 0$$

这样，我们就证明了原根由于复数单位根相同的性质。

3.公式

我们用原根替换复数单位根，得到：

yk=A(gk(p−1)n)=∑j=0n−1akgkj(p−1)n$$y_k=A(g^{\frac{k(p-1)}{n}})=\sum _{j=0}^{n-1}{a}_k{g}^{\frac{kj(p-1)}{n}}$$

即y=NTTn(a)$y=NT{T}_n(a)$。逆变换：

ak=1n∑j=0n−1ykg−kj(p−1)n$$a_k=\frac{1}{n}\sum _{j=0}^{n-1}{y}_k{g}^{\frac{-kj(p-1)}{n}}$$

即a=INTTn(y)$a=INT{T}_n(y)$。

7、其他扩展

1.任意模数FFT

有的时候我们需要卷积后模上一个数，这个数不是NTT模数，甚至可能不是个质数。那我们该怎么做呢？

这里使用拆系数FFT，本质是以时间换精度。

现在给定次数界为m$m$的两个多项式A(x),B(x)$A(x),B(x)$，要求A(x)B(x)modP$A(x)B(x)\mathrm{mod}P$。

首先，令M=⌊P−−√⌋$M=\lfloor \sqrt{P}\rfloor$，再对于每个ai$a_i$或bi$b_i$，把其变为kiM+ri(ri<M)$k_{i}M+r_i(r_i<M)$的形式。这样，ki$k_i$和ri$r_i$就都小于等于M$M$了。

那么卷积就可以写成：

ci=∑k=0iakbi−k=∑k=0i(kaiM+rai)(kbi−kM+rbi−k)=M2∑k=0ikaikbi−k+M∑k=0i(kairbi−k+raikbi−k)+∑k=0irairbi−k$$c_i=\sum _{k=0}^i{a}_k{b}_{i-k}=\sum _{k=0}^i(k_{a_i}M+r_{a_i})(k_{b_{i-k}}M+r_{b_{i-k}})=M^2\sum _{k=0}^i{k}_{a_i}{k}_{b_{i-k}}+M\sum _{k=0}^i(k_{a_i}{r}_{b_{i-k}}+r_{a_i}{k}_{b_{i-k}})+\sum _{k=0}^i{r}_{a_i}{r}_{b_{i-k}}$$

那么我们看到，ci$c_i$可以由k$k$和r$r$合并得到。那么我们对ka,kb,ra,rb$k_a,k_b,r_a,r_b$分别做FFT，再对应算出x1=kakb,x2=karb+rakb,x3=rarb$x_1=k_a{k}_b,x_2=k_a{r}_b+r_a{k}_b,x_3=r_a{r}_b$，对它们再分别做IFFT，就可以由c=M2x1+Mx2+x3$c=M^2{x}_1+M{x}_2+x_3$得到答案了。

这么做的好处究竟在哪里呢？事实上，计算后的最大值由mP$mP$变为了m⌊P−−√⌋$m\lfloor \sqrt{P}\rfloor$，避免了溢出。

时间复杂度：O(nlogn)$O(n\log n)$

2.多项式求逆

现在我们有一个次数界为n$n$的多项式A(x)$A(x)$，要求B(x)$B(x)$满足A(x)B(x)≡1(modxn)$A(x)B(x)\equiv 1(\mathrm{mod}{x}^n)$。

我们考虑倍增实现。

• 当n=1$n=1$时，直接求逆元求得B(x)≡A(x)p−2$B(x)\equiv A(x{)}^{p-2}$。
• 当n>1$n>1$时，已有A(x)G(x)≡1(modxn2)$A(x)G(x)\equiv 1(\mathrm{mod}{x}^{\frac{n}{2}})$：

因为：

A(x)B(x)≡1(modxn)$$A(x)B(x)\equiv 1(\mathrm{mod}{x}^n)$$

又因为除了0$0$次项之外到n−1$n-1$次都为0$0$，因此到n2−1$\frac{n}{2}-1$次项也为零：

A(x)B(x)≡1(modxn2)$$A(x)B(x)\equiv 1(\mathrm{mod}{x}^{\frac{n}{2}})$$

又

A(x)G(x)≡1(modxn2)$$A(x)G(x)\equiv 1(\mathrm{mod}{x}^{\frac{n}{2}})$$

两式相减：

A(x)[B(x)−G(x)]≡0(modxn2)$$A(x)[B(x)-G(x)]\equiv 0(\mathrm{mod}{x}^{\frac{n}{2}})$$

因为：

A(x)≠0$$A(x)\ne 0$$

所以：

B(x)−G(x)≡0(modxn2)$$B(x)-G(x)\equiv 0(\mathrm{mod}{x}^{\frac{n}{2}})$$

既然0$0$至n2−1$\frac{n}{2}-1$次项都为零，那么平方之后0$0$至n−1$n-1$次项也为零：

B(x)2−2B(x)G(x)+G(x)2≡0(modxn)$$B(x{)}^2-2B(x)G(x)+G(x{)}^2\equiv 0(\mathrm{mod}{x}^n)$$

又

A(x)B(x)≡1(modxn)$$A(x)B(x)\equiv 1(\mathrm{mod}{x}^n)$$

两边同时乘上A(x)$A(x)$，得：

B(x)−2G(x)+A(x)G(x)2≡0(modxn)$$B(x)-2G(x)+A(x)G(x{)}^2\equiv 0(\mathrm{mod}{x}^n)$$

移项，得：

B(x)≡2G(x)−A(x)G(x)2(modxn)$$B(x)\equiv 2G(x)-A(x)G(x{)}^2(\mathrm{mod}{x}^n)$$

这样就可以了。

时间复杂度：O(nlogn)$O(n\log n)$

3.多项式开根

前置：多项式求逆。

现在我们有一个次数界为n$n$的多项式A(x)$A(x)$，要求B(x)$B(x)$满足B(x)2≡A(x)(modxn)$B(x{)}^2\equiv A(x)(\mathrm{mod}{x}^n)$。

还是倍增。

• 当n=1$n=1$时，B(x)$B(x)$等于A(x)$A(x)$在模意义下的开根。
• 当n>1$n>1$时，已有G(x)2≡A(x)(modxn2)$G(x{)}^2\equiv A(x)(\mathrm{mod}{x}^{\frac{n}{2}})$：

因为：

G(x)2≡A(x)(modxn2)$$G(x{)}^2\equiv A(x)(\mathrm{mod}{x}^{\frac{n}{2}})$$

移项，得：

G(x)2−A(x)≡0(modxn2)$$G(x{)}^2-A(x)\equiv 0(\mathrm{mod}{x}^{\frac{n}{2}})$$

两边平方，同理可得：

G(x)4−2G(x)2A(x)+A(x)2≡0(modxn)$$G(x{)}^4-2G(x{)}^{2}A(x)+A(x{)}^2\equiv 0(\mathrm{mod}{x}^n)$$

所以：

[G(x)2+A(x)]2−4G(x)2A(x)≡0(modxn)$$[G(x{)}^2+A(x){]}^2-4G(x{)}^{2}A(x)\equiv 0(\mathrm{mod}{x}^n)$$

即：

[G(x)2+A(x)]2≡4G(x)2A(x)(modxn)$$[G(x{)}^2+A(x){]}^2\equiv 4G(x{)}^{2}A(x)(\mathrm{mod}{x}^n)$$

除过去：

[G(x)2+A(x)]24G(x)2≡A(x)(modxn)$$\frac{[G(x{)}^2+A(x){]}^2}{4G(x{)}^2}\equiv A(x)(\mathrm{mod}{x}^n)$$

得到：

A(x)≡(G(x)2+A(x)2G(x))2≡B(x)2(modxn)$$A(x)\equiv {\left(\frac{G(x{)}^2+A(x)}{2G(x)}\right)}^2\equiv B(x{)}^2(\mathrm{mod}{x}^n)$$

即：

B(x)≡G(x)2+A(x)2G(x)≡A(x)2G(x)+G(x)2(modxn)$$B(x)\equiv \frac{G(x{)}^2+A(x)}{2G(x)}\equiv \frac{A(x)}{2G(x)}+\frac{G(x)}{2}(\mathrm{mod}{x}^n)$$

这就可以了。

时间复杂度：O(nlogn)$O(n\log n)$

4.多项式求导

根据导数的可加性和幂函数求导法则d(cxk)dx=ckxk−1$\frac{\mathbb{d}(c{x}^k)}{\mathbb{d}x}=ck{x}^{k-1}$，有多项式的导数为：

d(∑k=0n−1akxk)dx=∑k=0n−1d(akxk)dx=∑k=1n−1kakxk−1=∑k=0n−2(k+1)ak+1xk$$\frac{\mathbb{d}(\sum _{k=0}^{n-1}{a}_k{x}^k)}{\mathbb{d}x}=\sum _{k=0}^{n-1}\frac{\mathbb{d}(a_k{x}^k)}{\mathbb{d}x}=\sum _{k=1}^{n-1}k{a}_k{x}^{k-1}=\sum _{k=0}^{n-2}(k+1)a_{k+1}{x}^k$$

时间复杂度：O(n)$O(n)$

5.多项式积分

根据积分的可加性和幂函数的不定积分∫cxkdx=ckxk+1${\displaystyle \int c{x}^k\mathbb{d}x=\frac{c}{k}{x}^{k+1}}$，有多项式的不定积分为：

∫∑k=0n−1akxkdx=∑k=0n−1∫akxkdx=∑k=0n−1akkxk+1=∑k=1nak−1k−1xk$$\int \sum _{k=0}^{n-1}{a}_k{x}^k\mathbb{d}x=\sum _{k=0}^{n-1}\int a_k{x}^k\mathbb{d}x=\sum _{k=0}^{n-1}\frac{a_k}{k}{x}^{k+1}=\sum _{k=1}^n\frac{a_{k-1}}{k-1}{x}^k$$

时间复杂度：O(n)$O(n)$

6.多项式求ln

前置：多项式求导&积分&求逆

现在已知多项式A(x)$A(x)$，要求B(x)=lnA(x)$B(x)=\ln A(x)$。我们两边微分，得到：

B′(x)=d(lnA(x))dA(x)dA(x)dx$$B^{\prime }(x)=\frac{\mathbb{d}(\ln A(x))}{\mathbb{d}A(x)}\frac{\mathbb{d}A(x)}{\mathbb{d}x}$$

B′(x)=A′(x)A(x)$$B^{\prime }(x)=\frac{A^{\prime }(x)}{A(x)}$$

再两边积分，就得到：

B(x)=∫A′(x)A(x)dx$$B(x)=\int \frac{A^{\prime }(x)}{A(x)}\mathbb{d}x$$

因此，我们直接多项式求导+求逆+积分解决。

时间复杂度：O(nlogn)$O(n\log n)$

7.多项式求exp

前置：多项式求ln

现在，我们已知多项式A=A(x)$A=A(x)$（这样写是因为在这里把A$A$视为是与x$x$无关的），求F(x)=exp(A)=eA$F(x)=\exp (A)=e^A$。只要我们设G(x)=lnx−A$G(x)=\ln x-A$，就得到：

G(F(x))=lnF(x)−A=0$$G(F(x))=\ln F(x)-A=0$$

我们考虑用牛顿迭代法来倍增解这个方程。

对于牛顿迭代法的初始解，即结果的常数项，我们并不知道具体值。但是如果不对的话，也只是缺少了一个常数罢了，那我们不妨设F(x)=1$F(x)=1$。

倍增：现在设我们已经求出了F(x)$F(x)$的前n$n$项F0(x)$F_0(x)$，即：

F0(x)≡F(x)(modxn)$$F_0(x)\equiv F(x)(\mathrm{mod}{x}^n)$$

根据牛顿迭代法，我们求出下一个近似解为：

F(x)≡F0(x)−G(F0(x))G′(F0(x))≡F0(x)−lnF0(x)−AF0(x)−1≡F0(x)(1−lnF0(x)+A(x))$$F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G^{\prime }(F_0(x))}\equiv F_0(x)-\frac{\ln F_0(x)-A}{F_0(x{)}^{-1}}\equiv F_0(x)(1-\ln F_0(x)+A(x))$$

如此，就可以倍增实现了。

时间复杂度：O(nlogn)$O(n\log n)$

8.多项式求幂

已知多项式A(x)$A(x)$和指数k$k$，求A(x)k$A(x{)}^k$。

在幂数很大的时候，为了确定出系数，时间复杂度为O(knlog(nk))$O(kn\log (nk))$，我们必须进行优化。

我们发现：

A(x)k=(elnA(x))k=eklnA(x)=exp(klnA(x))$$A(x{)}^k={\left(e^{\ln A(x)}\right)}^k=e^{k\ln A(x)}=\exp (k\ln A(x))$$

于是我们可以用多项式求ln+多项式求exp在O(nlogn)$O(n\log n)$的时间内求出。

8、代码实现

1.二进制反转

可以用一种类似dp的思想计算。

边界：0$0$的二进制翻转为0$0$

递归式：对于a$a$，已经算出了rev⌊a2⌋$re{v}_{\lfloor \frac{a}{2}\rfloor }$，a$a$就是除去最后一位的二进制翻转(rev⌊a2⌋$re{v}_{\lfloor \frac{a}{2}\rfloor }$)向后移动一位再补上第一位。即：

rev[a]=(rev[a>>1]>>1)|((a&1)<<(l-1))

l$l$为要翻转的位数。

2.复数类

套公式即可。

struct complex{
    double re,im;
    complex(double re=0,double im=0):re(re),im(im){}
    complex operator+(complex &b)const{return complex(re+b.re,im+b.im);}
    complex operator-(complex &b)const{return complex(re-b.re,im-b.im);}
    complex operator*(complex &b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
};

3.FFT

• 1：根据二进制翻转交换ar$a_r$和arev(r)$a_{rev(r)}$
• 2：枚举归并步长i∈[1,n)$i\in [1,n)$，i$i$为二的幂；
  
  • 2.1： 根据欧拉公式求出ω1i=cosπi+isinπi${\omega }_i^1=\cos \frac{\pi }{i}+i\sin \frac{\pi }{i}$
  • 2.2：枚举归并位置j$j$，归并[j,j+i)$[j,j+i)$和[j+i,j+2i)$[j+i,j+2i)$，步长为2i$2i$
    
    • 2.2.1：枚举x$x$的幂数k∈[0,i)$k\in [0,i)$进行蝴蝶操作计算y$y$，根据单位根计算ωki${\omega }_i^k$
      
      • 2.2.1.1：yj+k=yj+k+ωkiyj+k+i$y_{j+k}=y_{j+k}+{\omega }_i^k{y}_{j+k+i}$
      • 2.2.1.2：yj+k+i=yj+k−ωkiyj+k+i$y_{j+k+i}=y_{j+k}-{\omega }_i^k{y}_{j+k+i}$

注意：由于在C++中值会被更改，因此需要引入临时变量。

void FFT(complex c[]){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}

4.IFFT

由于公式中只差一个负号而已，因此引入一个参数type$type$，在欧拉公式的地方乘上去，再除以n$n$就可以了。

void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=m;i++)a[i].re/=n;
}

5.多项式乘法

模板：洛谷3083

注意：

1、由于FFT要求n$n$为2$2$的幂且结果的次数界较大，所以要把两个因式的系数补到l$l$位，l$l$满足l=2t$l=2^t$且l/2$l/2$大于等于因式的次数界。

2、FFT$FFT$虽然在数学上精准，但在C++中误差巨大，因此虚部不会为0$0$，忽略即可。实部也不为正数，可以加上0.1$0.1$再向下取整。

代码：

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const double pi=acos(-1);
struct complex{
    double re,im;
    complex(double re=0,double im=0):re(re),im(im){}
    complex operator+(complex &b)const{return complex(re+b.re,im+b.im);}
    complex operator-(complex &b)const{return complex(re-b.re,im-b.im);}
    complex operator*(complex &b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
}a[2097153],b[2097153];
int n,m,l,r[2097153];
void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=m;i++)a[i].re/=n;
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%lf",&a[i].re);
    for(int i=0;i<=m;i++)scanf("%lf",&b[i].re);
    m+=n;
    for(n=1;n<=m;n<<=1)l++;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a);
    FFT(b);
    for(int i=0;i<=n;i++)a[i]=a[i]*b[i];
    IFFT(a);
    for(int i=0;i<=m;i++)printf("%d ",int(a[i].re+0.5));
}

6.(I)FFT+(I)NTT

给出一个n$n$次多项式和一个m$m$次多项式的系数表达(1≤n,m≤4000000$1\le n,m\le 4000000$)，求乘积。type=0$type=0$时，直接计算；type=1$type=1$时，结果取模998244353$998244353$(原根为3$3$)。

注：为了方便阅读，代码效率不高。若要提速，可以把单位根打表。而且，由于gp−1n$g^{\frac{p-1}{n}}$中p−1n$\frac{p-1}{n}$必须为整数，故仅在第一个比n+m$n+m$大的2$2$的整数次幂是p−1$p-1$的约数时才可行，此处998244353−1=998244352=119×223$998244353-1=998244352=119\times 2^{23}$，223=8388608>n+m$2^{23}=8388608>n+m$。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
const double pi=acos(-1);
const int p=998244353;
typedef long long ll;
struct complex{
    double re,im;
    complex(double re=0,double im=0):re(re),im(im){}
    complex operator+(complex b)const{return complex(re+b.re,im+b.im);}
    complex operator-(complex b)const{return complex(re-b.re,im-b.im);}
    complex operator*(complex b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
}af[131073],bf[131073];
struct modp{
    int a;
    modp(int a=0):a(a){}
    modp operator+(modp b)const{return (a+b.a)%p;}
    modp operator-(modp b)const{return (a-b.a+p)%p;}
    modp operator*(modp b)const{return ll(a)*b.a%p;}
    modp operator/(modp b)const{return (b^(p-2))*a;}
    modp operator^(int b)const{
        modp ans=1,bs=a;
        while(b){
            if(b&1)ans=ans*bs;
            bs=bs*bs;
            b>>=1;
        }
        return ans;
    }
}an[131073],bn[131073];
const modp g=3;
int n,m,l,r[131073],type;
modp gn[18];
void init(){
    m+=n;
    for(n=1;n<=m;n<<=1)l++;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<=17;i++)gn[i]=g^((p-1)/(1<<i));
}
void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=n;i++)c[i].re/=n;
}
void NTT(modp c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1,id=1;i<n;i<<=1,id++){
        modp ggn=gn[id]^(tp==1?1:p-2);
        for(int j=0;j<n;j+=(i<<1)){
            modp gg=1;
            for(int k=0;k<i;k++,gg=gg*ggn){
                modp x=c[j+k],y=gg*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void INTT(modp c[]){
    NTT(c,-1);
    for(int i=0;i<=n;i++)c[i]=c[i]/n;
}
int main(){
    scanf("%d%d%d",&n,&m,&type);
    if(type==0){
        for(int i=0;i<=n;i++)scanf("%lf",&af[i].re);
        for(int i=0;i<=m;i++)scanf("%lf",&bf[i].re);
    }else{
        for(int i=0;i<=n;i++)scanf("%d",&an[i].a);
        for(int i=0;i<=m;i++)scanf("%d",&bn[i].a);
    }
    init();
    if(type==0){
        FFT(af);
        FFT(bf);
        for(int i=0;i<=n;i++)af[i]=af[i]*bf[i];
        IFFT(af);
        for(int i=0;i<=m;i++)printf("%d ",int(af[i].re+0.5));
    }else{
        NTT(an);
        NTT(bn);
        for(int i=0;i<=n;i++)an[i]=an[i]*bn[i];
        INTT(an);
        for(int i=0;i<=m;i++)printf("%d ",an[i].a);
    }
    printf("\n");
}

常数只有上面那个的三分之一的NTT（正式考试请务必采用这种写法）：

PS:有一道题上面那个3700ms,这个1000ms

inline ll pow(ll a,int b){
    ll ans=1;
    while(b){
        if(b&1)ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
inline ll add(ll a,ll b){return a+b>p?a+b-p:a+b;}
inline ll cut(ll a,ll b){return a-b<0?a-b+p:a-b;}
void init(){
    for(n=1;n<=s;n<<=1);
    nn=n;
    gn[0][0]=gn[1][0]=1;
    gn[0][1]=pow(g,(p-1)/(n<<1));
    gn[1][1]=pow(gn[0][1],p-2);
    for(int i=2;i<(n<<1);i++){gn[0][i]=gn[0][i-1]*gn[0][1]%p;gn[1][i]=gn[1][i-1]*gn[1][1]%p;}
    inv[1]=1;
    for(int i=2;i<=(n<<1);i++)inv[i]=inv[p%i]*(p-p/i)%p;
}
void NTT(ll c[],int n,int tp=1){
    for(int i=0;i<n;i++){
        r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
        if(i<r[i])swap(c[i],c[r[i]]);
    }
    for(int i=1;i<n;i<<=1){
        for(int j=0;j<n;j+=(i<<1)){
            for(int k=0;k<i;k++){
                ll x=c[j+k],y=gn[tp!=1][nn/i*k]*c[j+k+i]%p;
                c[j+k]=add(x,y);
                c[j+k+i]=cut(x,y);
            }
        }
    }
}
void INTT(ll c[],int n){
    NTT(c,n,-1);
    for(int i=0;i<n;i++)c[i]=c[i]*inv[n]%p;
}

7.任意模数FFT

给定n,m,P$n,m,P$，再给定次数界为n$n$的第一个多项式和次数界为m$m$的第二个多项式，求两个多项式的卷积模P$P$。

注意：拆系数FFT精度损失极大，需要使用long double保证正确性。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
typedef long long ll;
int x,n,m,l,r[262145],pm,p;
const long double pi=acos(-1);
struct complex{
    long double re,im;
    complex(long double re=0,long double im=0):re(re),im(im){}
    complex operator+(const complex &b)const{return complex(re+b.re,im+b.im);}
    complex operator-(const complex &b)const{return complex(re-b.re,im-b.im);}
    complex operator*(const complex &b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
}k1[262145],r1[262145],k2[262145],r2[262145],c1[262145],c2[262145],c3[262145];
void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=n;i++)c[i].re/=n;
}
void init(){
    m+=n;
    for(n=1;n<=m;n<<=1)l++;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
int main(){
    scanf("%d%d%d",&n,&m,&p);
    pm=sqrt(p);
    for(int i=0;i<=n;i++){
        scanf("%d",&x);
        k1[i]=x/pm;
        r1[i]=x%pm;
    }
    for(int i=0;i<=m;i++){
        scanf("%d",&x);
        k2[i]=x/pm;
        r2[i]=x%pm;
    }
    init();
    FFT(k1);
    FFT(r1);
    FFT(k2);
    FFT(r2);
    for(int i=0;i<=n;i++){
        c1[i]=k1[i]*k2[i];
        c2[i]=k1[i]*r2[i]+r1[i]*k2[i];
        c3[i]=r1[i]*r2[i];
    }
    IFFT(c1);
    IFFT(c2);
    IFFT(c3);
    for(int i=0;i<=m;i++){
        ll s1=ll(c1[i].re+0.5)%p*pm%p*pm%p,s2=ll(c2[i].re+0.5)%p*pm%p,s3=ll(c3[i].re+0.5)%p;
        printf("%lld ",((s1+s2)%p+s3)%p);
    }
}

7.多项式的运算

依赖关系：

直接按照公式打就好了。

我们先修改NTT：

void NTT(modp c[],int n,int tp=1){
    for(int i=0;i<n;i++){
        r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
        if(i<r[i])swap(c[i],c[r[i]]);
    }
    for(int i=1,id=1;i<n;i<<=1,id++){
        modp ggn=gn[id]^(tp==1?1:p-2);
        for(int j=0;j<n;j+=(i<<1)){
            modp gg=1;
            for(int k=0;k<i;k++,gg=gg*ggn){
                modp x=c[j+k],y=gg*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void INTT(modp c[],int n){
    NTT(c,n,-1);
    for(int i=0;i<n;i++)c[i]=c[i]/n;
}

求逆：

void inverse(modp c[],int n=n){
    static modp t[262145],tma[262145];
    t[0]=c[0]^(p-2);
    for(int k=2;k<=n;k<<=1){
        for(int i=0;i<(k<<1);i++)tma[i]=(i<k?c[i]:0);
        for(int i=(k>>1);i<(k<<1);i++)t[i]=0;
        NTT(tma,k<<1);
        NTT(t,k<<1);
        for(int i=0;i<(k<<1);i++)t[i]=t[i]*2-t[i]*t[i]*tma[i];
        INTT(t,k<<1);
    }
    memcpy(c,t,sizeof(int)*n);
}

开根(c[0]=1$c[0]=1$)：

void sqrt(modp c[],int n=n){
    static modp t[262145],tma[262145],tmb[262145];
    t[0]=1;
    for(int k=2;k<=n;k<<=1){
        for(int i=0;i<k;i++)tma[i]=t[i]*2;
        inverse(tma,k);
        for(int i=0;i<(k<<1);i++)tmb[i]=(i<k?c[i]:0);
        NTT(tma,k<<1);
        NTT(tmb,k<<1);
        for(int i=0;i<(k<<1);i++){
            modp tmp=tma[i];
            tma[i]=t[i];
            t[i]=tmp*tmb[i];
        }
        INTT(t,k<<1);
        for(int i=0;i<(k<<1);i++)t[i]=(i<k?t[i]+tma[i]/2:0);
    }
    memcpy(c,t,sizeof(int)*n);
}

求导：

void derivative(modp c[],int n=n){for(int i=0;i<n;i++)c[i]=c[i+1]*(i+1);}

积分：

void integrate(modp c[],int n=n){for(int i=n-1;i>=1;i--)c[i]=c[i-1]*inv[i];c[0]=0;}

求ln：

void ln(modp c[],int n=n){
    static modp t[262145];
    for(int i=0;i<(n<<1);i++)t[i]=(i<n?c[i]:0);
    derivative(t,n);
    inverse(c,n);
    NTT(t,n<<1);
    NTT(c,n<<1);
    for(int i=0;i<(n<<1);i++)c[i]=c[i]*t[i];
    INTT(c,n<<1);
    for(int i=n;i<(n<<1);i++)c[i]=0;
    integrate(c,n);
}

求exp：

void exp(modp c[]){
    static modp t[262145],ta[262145];
    t[0]=1;
    for(int k=2;k<=n;k<<=1){
        for(int i=0;i<(k<<1);i++)ta[i]=t[i];
        ln(ta,k);
        for(int i=0;i<k;i++)ta[i]=c[i]-ta[i];
        ta[0].a++;
        NTT(t,k<<1);
        NTT(ta,k<<1);
        for(int i=0;i<(k<<1);i++)t[i]=t[i]*ta[i];
        INTT(t,k<<1);
        for(int i=k;i<(k<<1);i++)t[i]=0;
    }
    memcpy(c,t,sizeof(modp)*n);
}

求幂：
我们在多项式求exp中假定了c[0]=1$c[0]=1$，那么如果常数项不是1$1$的话我们就把常数项变为1$1$在运算后再用快速幂变回来即可。

void pow(modp c[],int k){
    ln(c);
    for(int i=0;i<n;i++)c[i]=c[i]*k;
    exp(c);
}

9、拉格朗日反演

1.形式幂级数

对于任意一个域F$F$我们定义在其上的形式幂级数为：

f(x)=∑k=0∞akxk,ak∈F$$f(x)=\sum _{k=0}^{\mathrm{\infty }}{a}_k{x}^k,a_k\in F$$

记所有的形式幂级数为F[[x]]$F[[x]]$。

2.反演公式

拉格朗日反演是求关于函数方程的幂级数展开系数非常重要的工具，可以用于组合计数函数的系数提取。

公式内容：

这里[xn]f(x)$[x^n]f(x)$指取f(x)$f(x)$中xn$x^n$的系数。

若f(x),g(x)∈F[[x]]$f(x),g(x)\in F[[x]]$且f(g(x))=x$f(g(x))=x$，则称f(x)$f(x)$为g(x)$g(x)$的复合逆。满足：

[xn]g(x)=1n[x−1]1f(x)n(1)$$\begin{array}{}\text{(1)}& [x^n]g(x)=\frac{1}{n}[x^{-1}]\frac{1}{f(x{)}^n}\end{array}$$

特别的，如果f(x)=xϕ(x)$f(x)={\displaystyle \frac{x}{\varphi (x)}}$，那么有：

[xn]g(x)=1n[xn−1]ϕ(x)n(2)$$\begin{array}{}\text{(2)}& [x^n]g(x)=\frac{1}{n}[x^{n-1}]\varphi (x{)}^n\end{array}$$

公式的推导：
由式f(x)=xϕ(x)$f(x)={\displaystyle \frac{x}{\varphi (x)}}$，得ϕ(x)=xf(x)$\varphi (x)={\displaystyle \frac{x}{f(x)}}$，代入(2)$(2)$式可得：

[xn]g(x)=1n[xn−1](xf(x))n$$[x^n]g(x)=\frac{1}{n}[x^{n-1}]{\left(\frac{x}{f(x)}\right)}^n$$

公式的推广：

[xn]h(g(x))=1n[xn−1]h′(x)(xf(x))n$$[x^n]h(g(x))=\frac{1}{n}[x^{n-1}]h^{\prime }(x){\left(\frac{x}{f(x)}\right)}^n$$