在圆上作 n 个点将单位圆等分成 n 份（以单位圆与实轴正半轴的交点为一个等分点）。然后以原点为起点，圆上的这 n 个等分点为终点，作出 n 个向量。其中幅角最小的那个的向量被称为 n 次单位根，记做 \omega_n^1，其余的向量分别为 \omega_n^2，\omega_n^3，\dots，\omega_n^n。
单位根的代数意义是 x^n=1 的解，一定会有 n 个解，这 n 个解正是 \omega_n^1，\omega_n^2，\omega_n^3，\dots，\omega_n^n，其中实数解只有 1，有时有 -1，其余都是复数解。
所以除 1 以外的幅角最小的那个解一定是以 \frac{360^\circ}{n} 为幅角， 1 为模长的向量，这个解也就是单位根 \omega_n^1。同时，单位根幅角的倍数也一定是方程的解，因为当 \theta \times n=360^\circ，就有 k\times \theta \times n=k\times360^\circ，也是 360^\circ的倍数，也能在 n 次幂后会落在向量 1 上。 所有以单位根幅角的倍数为幅角，模长为 1 的向量一共有 n 个，那正是单位圆上的 n 等分点为终点的向量。
有了幅角和模长，怎么具体的把这个复数写出来呢？三角函数是个很有用的公式，高中知识告诉我们 sin\theta 表示当幅角为 \theta， 模长为 1 时的 y 坐标，cos\theta 表示 x 坐标。复平面中 x 和 y 分别表示 a 和 b 的值，所以我们就可以用三角函数将 n 次单位根表示出来，\omega_n^1=cos\frac{360^\circ}{n}+i\times sin\frac{360^\circ}{n}，同理 \omega_n^k=cos\frac{360^\circ\times k}{n}+i\times sin\frac{360^\circ\times k}{n}，因为之前说过 k 次方就是将幅角乘 k。我们还需要把公式换成弧度制，因为 C++ 的内置函数是弧度制的。弧度制很简单，只需要记住 \pi=180^\circ 即可。于是得到 \omega_n^k=cos\frac{2k\pi}{n}+i\times sin\frac{2k\pi}{n}，这个式子叫做欧拉公式，好像也可以用泰勒展开得到，让人感叹数学的魔幻。
上图中，OM 为 cos 值，ON 为 sin 值。
因为 FFT 中需要正负成对出现的值，所以我们还得看看单位根的正负关系。复数 a+bi 的相反数是 -a-bi，其几何意义是 xy 坐标同时取反，也就是以原点为中心，旋转 180^\circ。用向量表示就是模长相同，方向相反。这很重要，因为我们发现单位圆 n 等分，当 n 为偶数时，就能得到 \frac{n}{2} 对这样互为相反数的点，而它们都可以用单位根表示出来。
我们惊奇地发现这 \frac{n}{2} 个向量平方以后重新将圆等分为了 \frac{n}{2} 份，当然前提是 n 为偶数。如何理解？因为 n 为偶数，所以这前 \frac{n}{2} 个向量本身平分了 x 轴上方的半个单位圆。那么平方以后，幅角翻倍，就变成平分整个单位圆了。所以如果 \frac{n}{2} 还是偶数，就能继续分治下去，因为既然是平分单位圆，那就还是正负成对的。
上加一线表示共轭，也就是复数含 i 的项取反，前面的 \frac{1}{n} 表示矩阵内的每个项都除以 n。总体来说，其逆矩阵就是每项取共轭，再除以 n。
点值矩阵乘以这个逆矩阵得到的矩阵每项再除以 n 就是系数矩阵，但这个工作也是 O(n^2) 的，如何加速？细心的读者可能已经发现了，这个过程和傅里叶快速变换是一样的，完全可以用相同的分治来加速，因为共轭实质上就是按 x 轴翻折，单位圆上单位根的性质不变。代码实现中两者只有一行语句的差别，也就是单位根含 i 项的正负号。
代码实现
洛谷 P3803 FFT
#include <bits/stdc++.h>
#define N 5000000
#define err 1e-6
using namespace std;
int n, m;
const double pi = acos(-1.0);
complex<double> a[N], b[N];
void fft(complex<double> *a,int len,int inv){
if(len==1)return;
complex<double> e[len/2],o[len/2];
for(int i=0;i<len;i+=2)
e[i/2]=a[i],o[i/2]=a[i+1];
fft(e,len/2,inv); fft(o,len/2,inv);
complex<double> wn(cos(2*pi/len),sin(2*pi/len)*inv),wk(1,0);
for(int i=0;i<len/2;i++,wk*=wn)
a[i]=e[i]+wk*o[i],
a[i+len/2]=e[i]-wk*o[i];
}
int main(){
cin>>n>>m;
for(int i=0,buf;i<=n;i++)cin>>buf,a[i].real(buf);
for(int i=0,buf;i<=m;i++)cin>>buf,b[i].real(buf);
int len=1; while(len<n+m+1)len<<=1;
fft(a,len,1); fft(b,len,1);
for(int i=0;i<len;i++)a[i]*=b[i];
fft(a,len,-1);
for (int i=0;i<=n+m;++i)
printf("%.0f ",a[i].real()/len);
return 0;
}
参考资料
Fast Fourier Transform (2020), Youtube [online]. Available from:
https://www.youtube.com/watch?v=h7apO7q16V0&feature=emb_logo
[Accessed on June 1 2022]
一小时学会快速傅里叶变换 (2022), zhihu [online]. Available from:
https://zhuanlan.zhihu.com/p/31584464
[Accessed on May 31 2022]
快速傅里叶变换 (2022), zhihu [online]. Available from: https://zhuanlan.zhihu.com/p/347091298
[Accessed on May 31 2022]