快速傅里叶变换(FFT) 学习笔记

合集 - 学习笔记(8)

1.Manacher 学习笔记2024-05-042.李超线段树 学习笔记2024-07-053.闵可夫斯基和 学习笔记2024-08-024.2-SAT 学习笔记2024-08-105.KMP&exKMP2024-08-106.高斯消元 学习笔记2024-08-127.拉格朗日插值 学习笔记2024-08-12

8.快速傅里叶变换(FFT) 学习笔记2024-08-19

收起

快速傅里叶变换(FFT) 学习笔记

简介：

卷积是形如 $C_k=\sum _{i\oplus j==k}{A}_i\ast B_j$ 的式子，其中 $\oplus$ 为表示某种运算。

而快速傅里叶变换（FFT）用于在 $O(n\log n)$ 的时间复杂度内求加法卷积。

对于一个 $k$ 次多项式，如果我们知道了它在 $k+1$ 个点上的值，那么我们可以求出这个多项式。而两个多项式在同一个点上的值相乘的话，就可以得到 这两个多项式相乘后的多项式 在这个点上的值。FFT 的主要思路就是先将两个多项式在若干个点上的取值，然后相乘，再通过点值求出两个多项式相乘的结果。

单位根

先考虑求点值。我们肯定不能暴力选择点去代入求值，这样的时间复杂度是 $O(n^2)$ 的，要找一些有特殊性质的点去代入。

于是我们在复数域找到了单位根，它是在复平面上的单位圆上的点，表示为 ${\omega }_n^k$ ，它的意义是将单位圆分割成 $n$ 份，按顺指针方向的第 $k$ 个点所表示的值。两个复数相乘时，它们的辐角相加（从实轴的正半轴开始，按逆时针旋转），模长相乘（该点到原点的距离）。于是我们可以得到单位根的以下性质：

$$\begin{gathered} {\omega }_n^i\ast {\omega }_n^j={\omega }_n^{i+j} \\ {\omega }_n^k={\omega }_n^{k\mathrm{\%}n} \\ {\omega }_{n\ast i}^{k\ast i}={\omega }_n^k \\ {\omega }_n^{k+n/2}=-{\omega }_n^k（\mathrm{当}n\mathrm{为}\mathrm{偶}\mathrm{数}\mathrm{时}） \end{gathered}$$

对于最后一条，我们可以将 $n/2$ 看做旋转半个圆周，所以就是取相反数。

对于一个 $n$ 次多项式，我们看看将单位根代入会有什么性质：

先将多项式 $F[x]$ 拆成两个部分： $FL[x]=x^0+x^2+x^4...,FR[x]=x^1+x^3+x^5...$ 。然后就可以将多项式 $F[x]$ 这样表示：

$$F[x]=FL[x^2]+x\ast FR[x^2]$$

再分别将 ${\omega }_n^k$ 和 ${\omega }_n^{k+n/2}$ 代入：

$$\begin{array}{rl}F[{\omega }_n^k]& =FL[({\omega }_n^k{)}^2]+{\omega }_n^k\ast FR[({\omega }_n^k{)}^2]\\ & =FL[{\omega }_{n/2}^k]+{\omega }_n^k\ast FR[{\omega }_{n/2}^k]\\ F[{\omega }_n^{k+n/2}]& =FL[({\omega }_n^{k+n/2}{)}^2]+{\omega }_n^{k+n/2}\ast FR[({\omega }_n^{k+n/2}{)}^2]\\ & =FL[{\omega }_{n/2}^{k+n/2}]-{\omega }_n^k\ast FR[{\omega }_{n/2}^{k+n/2}]\\ & =FL[{\omega }_{n/2}^k]-{\omega }_n^k\ast FR[{\omega }_{n/2}^k]\end{array}$$

发现上下两式只有后面一项的符号不同，所以如果我们求出了 $FL[x]$ 和 $FR[x]$ 在 ${\omega }_{n/2}^{0..n/2-1}$ 上的值，我们就可以求出 $F[x]$ 在 ${\omega }_n^{0...n-1}$ 上的值了，这样做一次的时间复杂度是 $O(n)$ 的。同时我们可以发现，对于 $FL[x]$ 和 $FR[x]$ ，其实可以看做一个相同的子任务，再次递归下去求解。一共需要递归 $logn$ 次。

但是此时我们又有一个问题：如果某个 $n$ 不能被整除了怎么办？我们可以给它在后面补位，系数设为 0 就行。再观察一下我们递归到最下面一层时，系数的编号与第一层相比有什么特点。可以发现它是二进制位的反转，所以我们可以用数位DP的方法求出反转后的编号，这样就不用一遍一遍递归下去求了。

tr[i]=(tr[i>>1]>>1)|(i&1)?limit>>1:0

求多项式

通过上面的步骤，我们已经将 $n$ 个点的值求出来了，接下来就要通过这 $n$ 个点的值求多项式。

设点值构成的序列为 $G[x]$ 。然后我们可以列出以下式子：

$$G[x]=\sum _{i=0}^{n-1}{\omega }_n^i\ast F[i]$$

然后可以得到以下结论：

$$F[x]\ast n=\sum _{i=0}^{n-1}{\omega }_n^{-i}\ast G[i]$$

这个式子的证明就是将上面的式子代入。此时我们发现，这个式子和我们求点值的式子十分相似，只是在单位根上有点区别。我们可以发现，在单位圆上， ${\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;
}