前言
懒得写前言。希望学到后面的数学知识不要放弃，写在这里督促自己。（作为一个没有接触过复数的萌新）
感觉再不学点东西就真的可以退役了。
什么是 FFT？
FFT 是一种将多项式转换为其 点值表示 的算法。区别于暴力代入的 \(\mathcal{O}(n^2)\) 的不优秀复杂度，其时间复杂度为 \(\mathcal{O}(n\log n)\)。
P.S. 点值表示 是一个用 \(n+1\) 个点来表示一个 \(n\) 次多项式的表示方法。
FFT 有什么用？
作用：加速多项式乘法。
即有两个 \(n\) 次多项式 \(A\) 和 \(B\)，你要将它们相乘得到多项式 \(C\)。
朴素的做法是 \(\mathcal{O}(n^2)\) 的，即枚举 \(A\) 和 \(B\) 的每一项并相乘，但这样的复杂度显然是不优秀的。
那么我们可以通过什么样的方法来优化这个呢？这时候我们想到了点值表示。点值表示的好处是什么？它的好处在于 \(A\) 和 \(B\) 各取 \(2n\) 个点出来，相应地一乘，就可以得到 \(C\) 的点值表示，并且这个过程的时间复杂度是 \(\mathcal{O}(n)\) 的！
那么我们接下来需要解决的问题就是优化将多项式转化成点值表示及其逆操作，这就是 FFT 的作用。
DFT（离散傅里叶变换）
DFT 是 FFT 的朴素版算法。
复数相关数学知识
P.S. 因为笔者之前没有接触过复数相关知识，所以如果写得有纰漏，欢迎在评论区指出，不甚感激。
如果你之前学过复数，下面的就都不用看了。
你可以将复数理解成一个 向量。复数有 实部 和 虚部，例如复数 \(5+3i\) 中，\(5\) 就是实部，\(3\) 就是虚部，其中 \(i=\sqrt{-1}\)。复数 \(5+3i\) 就像向量 \((5,4)\) 一样。
但是复数也是一个数，可以进行运算，下面给出 复数相乘 的规则：模长相乘，幅角相加。什么意思？模长就是向量的模长，幅角就是一根射线从 \(x\) 轴正方向开始，逆时针旋转直至与向量重合所转过的角度。
C++ STL 中的食用方法：用 complex<double> x; 以定义，可以对其直接进行加减乘除。（P.S. 笔者不建议使用 STL 中的复数，建议手写复数，毕竟又不长，万一有毒瘤出题人卡 STL 呢）
算法介绍
那么复数在 DFT 中有什么用呢？
当然是傅里叶灵机一动，把这 \(n\) 个点值表示的复数取成 \(n\) 个模长为 \(1\) 的复数。
但是这 \(n\) 个复数并不是随便取的，而是在一个单位圆上 \(n\) 等分取出来的。一个点 \((x,y)\) 表示 \(x\) 为实部，\(y\) 为虚部构成的复数。如下图所示：（虽然是盗的图，嘻嘻。）
设点 \((1,0)\) 开始，逆时针将这 \(n\) 个点从 \(0\) 开始编号，第 \(k\) 个记为 \(\omega_n^k\)。显然 \(w_n^k=(w_n^1)^k\)，由 模长相乘，幅角相加 可以得到。称 \(w_n^1\) 为 \(n\) 次单位根。
根据第 \(\omega_n^k\) 的幅角为 \(\frac{k}{n}2\pi\) 可得它的向量为 \(\cos(\frac{k}{n}2\pi),\sin(\frac{k}{n}2\pi)\)，所以其对应的复数为 \(\cos(\frac{k}{n}2\pi)+i\cdot\sin(\frac{k}{n}2\pi)\)。
话说回来，DFT 是什么呢？实际上 DFT 就是把 \(\omega_n^0,\omega_n^1,...,\omega_n^{n-1}\) 代入得到的一种特殊的点值表示。
单位根的性质
P.S. 这两个性质都相当显然。
性质 1：\(\omega_{2n}^{2k}=\omega_n^k\)
证明：它们对应的向量相同。
性质 2：\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)
证明：它们对应的向量等大反向。
DFT 的性质及结论
相信你读到这里自然会产生疑问：为什么要选取单位根代入呢？当然是因为 DFT 有一些特殊的性质啦！
设 \((y_0,y_1,...,y_{n-1})\) 为 \(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\) 的 DFT 结果。我们在设一个多项式 \(B(x)=y_0+y_1x+y_2x^2+...+y_{n-1}x^{n-1}\)，将单位根的 倒数 代入得 \((z_0,z_1,...,z_{n-1})\)。
下面开始推柿子（其实你不推也没事的）：
根据 DFT 可知 \(y_k=\sum_{i=0}^{n-1}a_i(\omega_n^k)^i\)。
\[\begin{aligned} z_k&=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^i)^j)(\omega_n^{-k})^i\\ &=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i) \end{aligned} \]下面是最神奇的一步，考虑对 \(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\) 求值。若 \(j=k\)，即 \(j-k=0\) 时，它等于 \(n\)。其他情况下，可以用等比数列求和，即 \(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i=\frac{(\omega_n^{j-k})^n-1}{\omega_{n}^{j-k}-1}=\frac{(\omega_n^n)^{j-k}-1}{\omega_{n}^{j-k}-1}=\frac{1^{j-k}-1}{\omega_{n}^{j-k}-1}=0\)。
这说明了什么？这说明了除了 \(j=k\) 的情况，其它情况都是 \(0\)，都没有用！所以 \(z_k\) 就等于 \(na_k\)！，也就是 \(\displaystyle a_k=\frac{z_k}{n}\)。
那么我们就得到了一个结论：
把多项式 \(A(x)\) 的 DFT 结果作为另一个多项式的系数 \(B(x)\)，取单位根的倒数即 \((w_n^0,w_n^{-1},...,w_n^{-(n-1)})\) 为 \(x\) 代入 \(B(x)\) 后，得到的结果除以 \(n\) 就又得到了多项式 \(A(x)\) 的系数。
我们再来看一看这个结论。求 \(A(x)\) 的 DFT 是 \(\mathcal{O}(n)\) 的，然后将单位根倒数代入也是 \(\mathcal{O}(n)\) 的，后面的操作也是 \(\mathcal{O}(n)\) 的，总共就是 \(\mathcal{O}(n)\) 的了！那么这就实现了傅里叶逆变换。
P.S. 这里笔者对着前面讲的东西好好想了一会儿，理清楚了 DFT 的定义及其奇妙之处，希望读者也能理清楚这些东西，争取不要重头再来。
快速傅里叶变换（FFT）
虽然我们已经处理了 将点值表示转化为多项式 这一操作，但是我们还未处理 将多项式转化为点值表示，所以总共的时间复杂度仍然是 \(\mathcal{O}({n^2})\) 的。
此时，一种名为 快速傅里叶变换 的算法横空出世，这是一种 分治 算法。
P.S. 后面的数学证明看起来有一点点复杂，但是只要每一步都搞懂地学下去也不算难，希望公式恐惧症患者们不要放弃。还有，为了方便，我们下面都设 \(n\) 为 2 的整数幂，如果不是的话补 0 即可。
首先，我们还是设 \(A(x)=a_0+a_1x+a_2x^2+...a_{n-1}x^{n-1}\)。然后考虑按照下标地奇偶性分成奇偶两类：
\[A(x)=(a_0+a_2x^2+...a_{n-2}x^{n-2})+(a_1x+a_3x^3+...a_{n-1}x^{n-1}) \]设两个多项式 \(A_1(x)=(a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1})\) 和 \(A_2(x)=(a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}\)。
则可以用这两个多项式来表示 \(A(x)\)：
\[A(x)=A_1(x^2)+xA_2(x^2) \]那么假设要把 \(x=\omega_n^k\) 代入，其中 \(k<\frac{n}{2}\)，得到（P.S. 最后一行由 性质 1 得到）：
\[\begin{aligned} A(x)&=A(\omega_n^k)\\ &=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_{\frac{n}{2}}^k) \end{aligned} \]如果要分治，则必然要考虑求 \(x=\omega_n^{k+\frac{n}{2}}\) 的值，那么（P.S. 倒数第二行由 性质 1 得到，最后一行由 性质 2 得到）：
\[\begin{aligned} A(x)&=A(\omega_n^{k+\frac{n}{2}})\\ &=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})\\ &=A_1(\omega_n^{2k}\cdot\omega_n^n)+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k}\cdot\omega_n^n)\\ &=A_1(\omega_{n}^{2k})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^k)+\omega_{n}^{k+\frac{n}{2}}A_2(\omega_{\frac{n}{2}}^k)\\ &=A_1(\omega_{\frac{n}{2}}^k)-\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^k) \end{aligned} \]此时我们发现，如果我们知道 \(A_1(x)\) 和 \(A_2(x)\) 在 \((\omega_{\frac{n}{2}}^0,\omega_{\frac{n}{2}})^1,...,\omega_{\frac{n}{2}}^{\frac{n}{2}-1})\) 处的点值表示，就可以在 \(\mathcal{O}(n)\) 的时间内算出 \(A(x)\) 在 \((\omega_n^0,\omega_n^1,...,\omega_n^{n-1})\) 的点值表示了。这是一个 分治 的算法，总时间复杂度为 \(\mathcal{O}(n\log n)\)！
至此，对 FFT 基础理论部分的讲解已经结束。（P.S. 应该不算太难吧）
递归 FFT 的实现
利用分治的思想，我们可以写出一个递归 FFT。边界为 \(n=1\) 时直接 return。
P3803 【模板】多项式乘法（FFT） 代码如下：
#include<bits/stdc++.h>
#define int long long
#define repe(i,l,r) for(int (i)=l;(i)<=r;(i)++)
#define rept(i,n) for(int (i)=1;(i)<=n;(i)++)
#define rep(i,n) for(int (i)=0;(i)<n;(i)++)
#define FOR(i,r,l) for(int (i)=r;(i)>=l;(i)--)
#define pb push_back
#define mpr make_pair
#define lwb lower_bound
#define upb upper_bound
#define fi first
#define se second
#define sz(x) ((int)(x).size())
#define lowbit(x) x&-x
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
template<typename T>void die(T s){cout<<s<<endl;exit(0);}
template<typename T>void chmax(T& a,T b){if(a<b)a=b;return;}
template<typename T>void chmin(T& a,T b){if(a>b)a=b;return;}
int read(){int sum=0,f=1;char c;c=getchar();while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}while(c>='0'&&c<='9'){sum=(sum<<1)+(sum<<3)+(c-'0');c=getchar();}return sum*f;}
void out(int x){if(x<0){x=-x;putchar('-');}if(x>=10)out(x/10);putchar(x%10+'0');}
int fast(int a,int b,int P){int res=1;if(P<=0){while(b){if(b&1)res=res*a;a=a*a;b>>=1;}}else{while(b){if(b&1)res=res*a%P;a=a*a%P;b>>=1;}}return res;}
const int N=4e6+10;
const double pi=acos(-1.0);
int n,m;
struct cp{ //手写复数
double x,y;
cp(double _x=0,double _y=0){x=_x,y=_y;}
}a[N],b[N];
cp operator+(cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
cp operator-(cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
cp operator*(cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void fft(int lim,cp *arr,int typ){ //FFT，其中 typ 表示是 FFT 还是其逆变换
if(lim==1)return;
cp a1[lim>>1],a2[lim>>1];
for(int i=0;i<lim;i+=2){ //奇偶性分组
a1[i>>1]=arr[i];
a2[i>>1]=arr[i+1];
}
fft(lim>>1,a1,typ); // 分治
fft(lim>>1,a2,typ); // 分治
cp basen=cp(cos(2.0*pi/lim),typ*sin(2.0*pi/lim)),omg=cp(1.0,0.0); // basen 为单位根，omg 为当前的 omega 值
rep(i,lim>>1){
arr[i]=a1[i]+omg*a2[i]; // 根据之前推的柿子算出作左半边
arr[i+(lim>>1)]=a1[i]-omg*a2[i]; // 根据之前推的柿子算出作右半边
omg=omg*basen; // 更新当前 omega 值
}
}
signed main(){
n=read(),m=read();
rep(i,n+1)a[i].x=(double)read();
rep(i,m+1)b[i].x=(double)read();
int lim=1;
while(lim<=n+m)lim<<=1; // 补成 2 的整数次幂
fft(lim,a,1); // 正变换
fft(lim,b,1); // 正变换
rep(i,lim+1)a[i]=a[i]*b[i]; // 点值表示成绩
fft(lim,a,-1); // 逆变换
rep(i,n+m+1){
out((int)(a[i].x/lim+0.5)); // 按照刚才推的柿子，要除以 n
putchar(' ');
}
return 0;
}
//https://www.luogu.com.cn/problem/P3803
将这个代码提交到洛谷上，就可以得到 \(\color{green}{\mathtt{AC}}\) 的好成绩！
P.S. 笔者看到有些人说递归版 FFT 过不了，实践可过，并且不用开 O2。但是即使如此，递归版 FFT 还是很慢的，这个代码最慢跑到了 1.01s，仍需要优化。
优化：非递归版 FFT
首先，FFT 是每一次按奇偶性分组，例如：
| 状态 | 序列 |
|---|---|
| 初始 | \((0,1,2,3,4,5,6,7)\) |
| 第一轮操作后 | \((0,2,4,6\mid1,3,5,7)\) |
| 第二轮操作后 | \((0,4\mid2,6\mid1,5\mid 3,7)\) |
| 第三轮操作后 | \((0,4,2,6,1,5,3,7)\) |
经过一顿超出人类智商范围的找规律后得到：每个数最后的位置为其二进制反转后的位置。什么意思？例如 \(6(110_2)\) 最后的位置在 \(3(011_2)\)，\(4(100_2)\) 最后的位置在 \(1(001_2)\)。
这样就可以先按照最终的顺序摆放，然后每次迭代合并即可。但是还有一个问题，就是之前的 \(A_1(\omega_{\frac{n}{2}}^k)-\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^k)\) 和 \(A_1(\omega_{\frac{n}{2}}^k)-\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^k)\) 如何处理。容易发现，这个分别变成了分成两组的对应位置。这是什么意思呢？比如说序列 \((0,1,2,3,4,5,6,7)\) 中，\(0\) 和 \(1\) 分别在奇偶组的第一个位置，那么对应到迭代的序列中也是一样的，即代码中的 i+j 和 i+j+mid。这个操作被称为 蝴蝶操作。
下面呈上模板题的代码：
#include<bits/stdc++.h>
#define int long long
#define repe(i,l,r) for(int (i)=l;(i)<=r;(i)++)
#define rept(i,n) for(int (i)=1;(i)<=n;(i)++)
#define rep(i,n) for(int (i)=0;(i)<n;(i)++)
#define FOR(i,r,l) for(int (i)=r;(i)>=l;(i)--)
#define pb push_back
#define mpr make_pair
#define lwb lower_bound
#define upb upper_bound
#define fi first
#define se second
#define sz(x) ((int)(x).size())
#define lowbit(x) x&-x
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
template<typename T>void die(T s){cout<<s<<endl;exit(0);}
template<typename T>void chmax(T& a,T b){if(a<b)a=b;return;}
template<typename T>void chmin(T& a,T b){if(a>b)a=b;return;}
int read(){int sum=0,f=1;char c;c=getchar();while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}while(c>='0'&&c<='9'){sum=(sum<<1)+(sum<<3)+(c-'0');c=getchar();}return sum*f;}
void out(int x){if(x<0){x=-x;putchar('-');}if(x>=10)out(x/10);putchar(x%10+'0');}
int fast(int a,int b,int P){int res=1;if(P<=0){while(b){if(b&1)res=res*a;a=a*a;b>>=1;}}else{while(b){if(b&1)res=res*a%P;a=a*a%P;b>>=1;}}return res;}
const int N=4e6+10;
const double pi=acos(-1.0);
int n,m,r[N];
struct cp{
double x,y;
cp(double _x=0,double _y=0){x=_x,y=_y;}
}a[N],b[N];
cp operator+(cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
cp operator-(cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
cp operator*(cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void fft(int lim,cp *arr,int typ){
rep(i,lim)if(i<r[i])swap(arr[i],arr[r[i]]); //交换得到最终的序列
for(int mid=1;mid<lim;mid<<=1){ //枚举迭代长度
cp basen(cos(pi/mid),typ*sin(pi/mid)); //单位根
for(int i=0,R=mid<<1;i<lim;i+=R){ //枚举分组
cp omg(1.0,0.0); //当前的 omega 值
rep(j,mid){
cp x=arr[i+j],y=omg*arr[i+mid+j]; //上文所提到的蝴蝶操作
arr[i+j]=x+y; //按照公式更新
arr[i+mid+j]=x-y; //按照公式更新
omg=omg*basen; //更新 omega 值
}
}
}
}
signed main(){
n=read(),m=read();
rep(i,n+1)a[i].x=(double)read();
rep(i,m+1)b[i].x=(double)read();
int lim=1,cnt=0;
while(lim<=n+m){
lim<<=1;
cnt++;
}
rep(i,lim)r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1)); //预处理出每个二进制倒转得到的二进制
fft(lim,a,1);
fft(lim,b,1);
rep(i,lim+1)a[i]=a[i]*b[i];
fft(lim,a,-1);
rep(i,n+m+1){
out((int)(a[i].x/lim+0.5));
putchar(' ');
}
return 0;
}
//https://www.luogu.com.cn/problem/P3803 优化版
该代码最慢已经被优化至 0.764s。
总结
相信读到这里的读者们已经对 FFT 有了清晰的了解，体会到了 FFT 的奇妙之处。
在这里，再来梳理一下 FFT 的思考步骤。
首先，解决的问题是多项式相乘（将 \(x=10\) 代入可以优化高精度乘法）。发现点值表示后需要解决两个问题：多项式 \(\Rightarrow\) 点至表示 以及 点值表示 \(\Rightarrow\) 多项式。
分别解决这两个问题都是无与伦比的，上面也都介绍了这两个问题的解决方案，确实十分巧妙。
至此，就是 FFT 的全部
参考资料
Rabbit_Hu 的博客（主要参考，图也是这里的）
attack 的博客（参考了下实现方法）