快速傅里叶变换(FFT) 学习笔记
快速傅里叶变换(FFT) 学习笔记
简介：
卷积是形如 \(C_k=\sum_{i\oplus j==k} A_i*B_j\) 的式子，其中 \(\oplus\) 为表示某种运算。
而快速傅里叶变换（FFT）用于在 \(O(n\log n)\) 的时间复杂度内求加法卷积。
对于一个 \(k\) 次多项式，如果我们知道了它在 \(k+1\) 个点上的值，那么我们可以求出这个多项式。而两个多项式在同一个点上的值相乘的话，就可以得到 这两个多项式相乘后的多项式 在这个点上的值。FFT 的主要思路就是先将两个多项式在若干个点上的取值，然后相乘，再通过点值求出两个多项式相乘的结果。
单位根
先考虑求点值。我们肯定不能暴力选择点去代入求值，这样的时间复杂度是 \(O(n^2)\) 的，要找一些有特殊性质的点去代入。
于是我们在复数域找到了单位根，它是在复平面上的单位圆上的点，表示为 \(\omega_n^k\) ，它的意义是将单位圆分割成 \(n\) 份，按顺指针方向的第 \(k\) 个点所表示的值。两个复数相乘时，它们的辐角相加（从实轴的正半轴开始，按逆时针旋转），模长相乘（该点到原点的距离）。于是我们可以得到单位根的以下性质：
对于最后一条，我们可以将 \(n/2\) 看做旋转半个圆周，所以就是取相反数。
对于一个 \(n\) 次多项式，我们看看将单位根代入会有什么性质：
先将多项式 \(F[x]\) 拆成两个部分： \(FL[x]=x^0+x^2+x^4...,FR[x]=x^1+x^3+x^5...\) 。然后就可以将多项式 \(F[x]\) 这样表示：
再分别将 \(\omega_n^k\) 和 \(\omega_n^{k+n/2}\) 代入：
发现上下两式只有后面一项的符号不同，所以如果我们求出了 \(FL[x]\) 和 \(FR[x]\) 在 \(\omega_{n/2}^{0..n/2-1}\) 上的值，我们就可以求出 \(F[x]\) 在 \(\omega_n^{0...n-1}\) 上的值了，这样做一次的时间复杂度是 \(O(n)\) 的。同时我们可以发现，对于 \(FL[x]\) 和 \(FR[x]\) ，其实可以看做一个相同的子任务，再次递归下去求解。一共需要递归 \(log\ n\) 次。
但是此时我们又有一个问题：如果某个 \(n\) 不能被整除了怎么办？我们可以给它在后面补位，系数设为 0 就行。再观察一下我们递归到最下面一层时，系数的编号与第一层相比有什么特点。可以发现它是二进制位的反转，所以我们可以用数位DP的方法求出反转后的编号，这样就不用一遍一遍递归下去求了。
tr[i]=(tr[i>>1]>>1)|(i&1)?limit>>1:0
求多项式
通过上面的步骤，我们已经将 \(n\) 个点的值求出来了，接下来就要通过这 \(n\) 个点的值求多项式。
设点值构成的序列为 \(G[x]\) 。然后我们可以列出以下式子：
然后可以得到以下结论：
这个式子的证明就是将上面的式子代入。此时我们发现，这个式子和我们求点值的式子十分相似，只是在单位根上有点区别。我们可以发现，在单位圆上， \(\omega_n^{-1}\) 所在的点和 \(\omega_n^1\) 所在的点，它们关于实轴对称。于是我们也可以用同样的方法来求上面的式子。最后记得除以 \(n\) 就是我们最终的多项式。
Code
#include<cstdio>
#include<cmath>
using namespace std;
const int N=3e6+5;const double PI=6.28318530717958646;
struct Date{
double x,y;
Date operator*(Date v) {
Date res;res.x=x*v.x-y*v.y;res.y=x*v.y+y*v.x;return res;
};
Date operator+(Date v) {
Date res;res.x=x+v.x;res.y=y+v.y;return res;
};
Date operator-(Date v) {
Date res;res.x=x-v.x;res.y=y-v.y;return res;
};
};Date f[N],g[N],w[N];
int tr[N],n,m;
void swap(double &x,double &y)
{
double t=x;x=y;y=t;
}
void FFT(int flag)
{
int i,j,len;Date t;
w[0].x=1;w[0].y=0;
w[1].x=cos(PI/n);w[1].y=sin(PI/n)*flag;
for(i=1;i<n;i++)
w[i]=w[i-1]*w[1];
for(i=0;i<n;i++)
{
if(i<tr[i])
{
swap(f[i].x,f[tr[i]].x);swap(f[i].y,f[tr[i]].y);
}
}
for(len=1;len<n;len=len<<1)
{
for(i=0;i<n;i+=2*len)
{
for(j=0;j<len;j++)
{
t=f[i+j+len]*w[n/len/2*j];
f[i+j+len]=f[i+j]-t;f[i+j]=f[i+j]+t;
}
}
}
}
int main()
{
int i,h,len;
scanf("%d%d",&n,&m);
for(i=0;i<=n;i++)
scanf("%lf",&f[i].x);
for(i=0;i<=m;i++)
scanf("%lf",&g[i].x);
h=0;len=n+m;
while((1<<h)<n+m+1)
h++;
n=1<<h;
for(i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
FFT(1);
for(i=0;i<n;i++)
{
swap(f[i].x,g[i].x);swap(f[i].y,g[i].y);
}
FFT(1);
for(i=0;i<n;i++)
f[i]=f[i]*g[i];
FFT(-1);
for(i=0;i<=len;i++)
printf("%d ",int(f[i].x/n+0.49));
return 0;
}
