快速傅里叶变换（Fast Fourier Transform,FFT）在算法竞赛中主要用于求卷积，或者说多项式乘法。如果我们枚举两数的各系数相乘，时间复杂度是\(O(n^2)\)，而FFT可以将这一过程优化到\(O(n\log n)\)。
流程
整个FFT算法分\(3\)个过程：
- 将\(2\)个多项式的系数表示法转换为点值表示法（FFT）。
- 将\(2\)个多项式的点值相乘，得到结果的点值表示。
- 将结果的点值表示转换为系数表示法（IFFT）。
其中系数表示法就是形如\(\sum\limits_{i=0}^{n-1}a_i x^i\)的表示。
将其转换为点值表示，就是选了\(n\)个具体的点\(x_0,x_1,\dots,x_{n-1}\)，对应在该点的值\(y_1,y_2,\dots,y_{n-1}\)。
可以证明，任意互不相同的\(n\)个点值可以确定一个唯一的系数表示。
其中，我们发现步骤\(2\)仅仅需要将每个\(x\)对应的的\(2\)个\(y\)相乘，就可以得到结果的点值表示了，这一步骤是\(O(n)\)的。
所以我们重点分析步骤\(1\)和\(3\)。
注：以下步骤假设\(n\)是\(2\)的正整数次幂，具体代码实现中，为了保证正常运行，也会将\(n\)补成\(2\)的正整数次幂。
系数转点值
现在我们要选取\(n\)个点，将\(A,B\)两个多项式转为点值表示，这里拿\(A\)说明。
如果我们直接代数进去，时间复杂度仍然是\(O(n^2)\)。
接下来就要上FFT的核心了。我们可以将\(A\)拆成\(1\)个偶函数和\(1\)个奇函数的和，也就是把\(A\)按奇偶次项划分成\(2\)部分。
\[\begin{aligned} A(x)&=a_0+a_1x^1+a_2x^2+a_3x^3+\dots+a_{n-1}x^{n-1}\\ &=(a_0+a_2x^2+a_4x^4+\dots+a_{n-2}x^{n-2})+(a_1+a_3x_3+a_5x^5+\dots+a_{n-1}x^{n-1}) \end{aligned}\]设
\[A_1(x)=a_0+a_2x+a_4x^2+\dots+a_{n-2}x^{\frac{n-2}{2}} \]\[A_2(x)=a_1+a_2x+a_5x^2+\dots+a_{n-1}x^{\frac{n-2}{2}} \]则
\[A(x)=A_1(x^2)+x\times A_2(x^2) \]\[A(-x)=A_1(x^2)-x\times A_2(x^2) \]为什么呢？因为\(A_1\)是偶函数，而\(x\times A_2\)是奇函数。
我们发现如果我们令\(x_0,x_1,x_2,\dots,x_{n-1}\)是正负配对的，其中\([x_0,\dots,x_{\frac{n}{2}-1}]\)和\([x_\frac{n}{2},\dots,x_{n-1}]\)对应相反，而它们的平方也对应相等……
我们就仅需算出\(A_1(x_0^2),A_1(x_1^2),\dots,A_1(x_{\frac{n}{2}-1}^2)\)和\(A_2(x_0^2),A_2(x_1^2),\dots,A_2(x_{\frac{n}{2}-1}^2)\)，然后用上面的式子就可以\(O(n)\)计算出每个\(A(x_0),A(x_1),\dots,A(x_{n-1})\)了。
这相当于我们将该问题分成了\(2\)个规模减半的子问题。对于每个子问题，再重复上面的操作。
如果这个方法可行，我们就能再\(O(n\log n)\)的时间复杂度内求出点值表达。
但我们发现对于初始状态，我们可以令\(x_0,x_1,x_2,\dots,x_{n-1}\)正负配对。但是在实数范围内它们的平方\([x_0^2,x_1^2,\dots,x_{\frac{n}{2}-1}^{2}]\)都是非负数，怎么能保证正负配对呢？
没错，就是复数。
如上图，上一层的平方是下一层，而下一层可用于推出上一层，和我们刚才的过程相同。
我们发现，最终要求的\(x_0,x_1,\dots,x_{n-1}\)，就是\(x^8=k\)在复数域上的解。
这里不妨取\(k=1\)，这样我们解出的\(8\)个根，就是“\(n\)次单位复数根”。
如上图，根据欧拉公式，有\(e^{i\theta}=\cos \theta+i\sin \theta\)。那么我们可以得出对于\(k\in[0,n)\)，\(n\)次单位根\(\omega_n^{k}=\omega^{2\pi i*\frac{k}{n}}\)（因为\(2\pi\)是\(360^{\circ}\)嘛，\(k\)表示的就是\(n\)份中的第几份）
单位根\(\omega\)有如下性质（可以结合下面的图来理解）：
\(\omega_{n}^{k}=-\omega_{n}^{k+\frac{n}{2}}\)
\(\omega_{an}^{ak}=\omega_{n}^{k}\)
发现了么？第\(1\)个性质和\(x\)的对应关系\(x_{i}=x_{i+\frac{n}{2}}\)正好相同！
而第\(2\)个性质说明平方之后，这\(\frac{n}{2}\)个单位根仍然平均分布在单位圆上。
也就是说，我们可以用\(\omega\)代替\(x\)：
至此，我们可以写出代码了。不过写之前我们发现有一个小优化：
如果我们每层额外用一个\(O(n)\)来进行奇偶分类，虽然不影响总时间复杂度\(O(n\log n)\)，但是常数较大。我们试着观察奇偶分类后的形态（图片来源自为风月马前卒的博客）。
奇偶分类后的序列，下标的二进制表示正是原序列下标二进制表示的翻转。所以我们根据这个规则，先把\(A\)的系数数组按这样的规则重排一下，然后按正常步骤进行就可以了。
点值转系数
其实我们的多项式也可以表示为矩阵乘的形式：
\[Y=X\times A \]即
\[\begin{bmatrix} y_0\\y_1\\y_2\\ \vdots\\y_{n-1} \end{bmatrix}=\begin{bmatrix} 1&x_0&x_0^1&\cdots&x_0^{n-1}\\ 1&x_1&x_1^1&\cdots&x_1^{n-1}\\ 1&x_2&x_2^1&\cdots&x_2^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&x_{n-1}&x_{n-1}^1&\cdots&x_{n-1}^{n-1} \end{bmatrix}\times \begin{bmatrix} a_0\\a_1\\a_2\\ \vdots\\a_{n-1} \end{bmatrix}\]带入\([x_0,x_1,\dots,x_{n-1}]=[1,\omega_n^1,\dots,\omega_n^{n-1}]\)：
\[\begin{bmatrix} y_0\\y_1\\y_2\\ \vdots\\y_{n-1} \end{bmatrix}=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega_n^1&\omega_n^2&\cdots&\omega_n^{n-1}\\ 1&\omega_n^2&\omega_n^4&\cdots&\omega_n^{2(n-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\cdots&\omega_n^{(n-1)(n-1)} \end{bmatrix}\times \begin{bmatrix} a_0\\a_1\\a_2\\ \vdots\\a_{n-1} \end{bmatrix}\]现在我们要倒推\(A\)，就两边同时左乘\(X^{-1}\)：
\[\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega_n^1&\omega_n^2&\cdots&\omega_n^{n-1}\\ 1&\omega_n^2&\omega_n^4&\cdots&\omega_n^{2(n-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\cdots&\omega_n^{(n-1)(n-1)} \end{bmatrix}^{-1}\times \begin{bmatrix} y_0\\y_1\\y_2\\ \vdots\\y_{n-1} \end{bmatrix}= \begin{bmatrix} a_0\\a_1\\a_2\\ \vdots\\a_{n-1} \end{bmatrix}\]中间的矩阵\(X\)是一个范德蒙德矩阵的转置矩阵，其行列式是\(\prod\limits_{1\le i < j \le n }(x_j-x_i)\)（因为是线性代数的范畴所以就不给证明了，逆阵也是）。
由于\(x_i\)互不相等，所以行列式不为\(0\)，存在逆阵：
\[X^{-1}=\frac{1}{n}\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega_n^{-1}&\omega_n^{-2}&\cdots&\omega_n^{-(n-1)}\\ 1&\omega_n^{-2}&\omega_n^{-4}&\cdots&\omega_n^{-2(n-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_n^{-(n-1)}&\omega_n^{-2(n-1)}&\cdots&\omega_n^{-(n-1)(n-1)} \end{bmatrix}^{-1} \]可以发现，与原矩阵的区别，就是每个\(\omega\)都取了倒数，最后再乘一个\(\frac{1}{n}\)。
所以代码不用重新写了，只需要额外加一个参数type。FFT提供1，IFFT提供-1。IFFT中，从\(\omega_n^0\)逆时针走\(1\)个单位，要变成顺时针。所以单位的复数部分需要\(\times\)type。具体见代码。
Code
点击查看代码
#include<bits/stdc++.h>
#define N 4000010
#define Pi acos(-1)
using namespace std;
struct complx{
double x,y;
complx(double xx=0,double yy=0){x=xx,y=yy;}
}a[N],b[N];
int l,r[N],limit=1;
inline complx operator+(complx a,complx b){return complx(a.x+b.x,a.y+b.y);}
inline complx operator-(complx a,complx b){return complx(a.x-b.x,a.y-b.y);}
inline complx operator*(complx a,complx b){return complx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void fft(complx *a,int type){
for(int i=0;i<limit;i++) if(i<r[i]) swap(a[i],a[r[i]]);
for(int mid=1;mid<limit;mid<<=1){
complx Wn(cos(Pi/mid),type*sin(Pi/mid));
for(int R=mid<<1,j=0;j<limit;j+=R){
complx w(1,0);
for(int k=0;k<mid;k++,w=w*Wn){
complx x=a[j+k],y=w*a[j+mid+k];
a[j+k]=x+y;
a[j+mid+k]=x-y;
}
}
}
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);
int n,m;
cin>>n>>m;
for(int i=0;i<=n;i++) cin>>a[i].x;
for(int i=0;i<=m;i++) cin>>b[i].x;
while(limit<=n+m) limit<<=1,l++;
for(int i=0;i<limit;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
fft(a,1);
fft(b,1);
for(int i=0;i<limit;i++) a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<=n+m;i++) cout<<(int)(a[i].x/limit+0.5)<<" ";
return 0;
}