[模板] 快速傅里叶变换(FFT)

多项式

假设有$n$次多项式:
$$\sum _{i=0}^n{a}_i{x}^i$$
如何表示这个多项式.
首先想到的是用所有系数表示，即:
{ a n , a n − 1 , a n − 2 , . . . , a 2 , a 1 , a 0 } \{a_n,a_{n-1},a_{n-2},...,a_2,a_1,a_0\} {an​,an−1​,an−2​,...,a2​,a1​,a0​}
我们将这样的表示方法为系数表示法.
但是，系数表示法在表示多项式相加和相乘很不方便.
譬如又有$n$次表达式:
$$\sum _{i=0}^n{b}_i{x}^i$$
那么将两式子相乘的时间复杂度是$O(n^2)$的.
如果$|a|$=$10^5$，这种算法就不是那么优越.
所幸，我们还有另一种表示方法——点值表示法.
设想，如果我们用足够多的值对(x,y) 表示这个函数:
$$y=\sum _{i=0}^n{a}_i{x}^i$$
表达式的系数是可以确定的.
更美妙的是，在点值表达式下，多项式加法和乘法是可以在$O(n)$实现的!
比如，有点对
$$(a_1,X_1)(a_2,X_2)...(a_n,X_n)$$
$$(a_1,Y_1)(a_2,Y_2)...(a_n,Y_n)$$
那么相乘得
$$(a_1,X_1{Y}_1)(a_2,X_2{Y}_2)...(a_n,X_n{Y}_n)$$
妙啊

转换

但是，我们在实际生活中并不经常使用点值表示法，而是经常使用系数表示法.
于是，便有了两者之间的转换.
系数表示法和点值表示法可以用如下矩阵乘法转换:

$$\left[\begin{array}{c}1x_1{x}_1^2...x_1^n\\ \\ 1x_2{x}_2^2...x_2^n\\ \\ .....\\ \\ 1x_n{x}_n^2...x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ \\ a_1\\ \\ ...\\ \\ a_n\end{array}\right]=\left[\begin{array}{c}{X}_1\\ \\ X_2\\ \\ ...\\ \\ X_n\end{array}\right]$$
这种变换，叫作离散傅里叶变换(DFT).
反之，由点值表示法退回系数表示法，叫作离散傅里叶逆变换(IDFT).
花了好久点个赞呗
但是发现，这样转换时间复杂度还是$O(n^2)$的!
所以，就有了快速傅里叶变换(FFT).

快速傅里叶变换

又称(Fast Fourier Transform)(FFT)，详情见Mr.Du.

快速傅里叶变换是离散傅里叶变换的快速版本.这还用说吗
它可以以优秀的$O(nlogn)$解决这个问题.
举个例子:
我们代入$a=1$这个特殊值，那么可以得到$X={\sum }_{i=0}^n{a}_i$.
我们代入$a=-1$、$a=0$也可以得到特殊值.
但是，这样的特殊值总是有限的.
所以，我们将眼光转向别处——复数域.
比如说有一个复数坐标系:

这个圆叫作单位圆没错我们见过，上面的点都可以乘方到1.
如果想要乘方$n$次，那么这些值就是这样分布的:

这是$n=8$的情况.最好画了不是吗
//下面所有$n$都是$2$的幂，注意！
我们称这样的数列为 { ω 1 , ω 2 , . . . , ω n } \{\omega_1,\omega_2,...,\omega_n\} {ω1​,ω2​,...,ωn​}.(当然是乘方n次后为1的 )
并称这些为 $1$的$n$次单位根 .

铺垫

由神奇而美妙的欧拉公式:

$e^{ix}=\cos x+i\sin x$

还有复数的极角表达式:

$(a,{\theta }_1)(b,{\theta }_2)=(ab,{\theta }_1+{\theta }_2)$

以及复数运算满足交换律和结合律.
我们就可以以这些为铺垫，进行接下来的研究了 ?

定理

下面的定理不再证明，可以结合奆老博客和Mr.Du证明:
//提示:欧拉定理和复平面(上图)是个好东西.

${\omega }_n^1=\cos (2\pi /n)+i\sin (2\pi /n)$
对于${\omega }_n$序列，${\omega }_k={\omega }_1^k$(很实用)
${\omega }_{2n}^{2k}={\omega }_n^k$
${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$
${\omega }_n^{k+n}={\omega }_n^k$

算法构建

经过这么一波操作之后，好像什么都没做…
因为直到现在，我们还是没有在这个数据上构造优化算法.
构造那么神奇的数据，不实现又有何用
要不我们先将一个${\omega }_n$代入试试.
比如我们设定
$$A(n)=\sum _{i=0}^n{a}_n{\omega }_n^i$$
看上去没什么变化…
但是我把它的奇数指数和偶数指数分开，分别写成:
$$A_1(n)=a_0+a_2{\omega }_n+...+a_n{\omega }_n^{\frac{n}{2}}$$
和
$$A_2(n)=a_1{\omega }_n+a_3{\omega }_n^2+...+a_{n-1}{\omega }_n^{\frac{n}{2}}$$
//$n$是$2$的幂.
所以我们震惊的发现，竟然有这样的表达式！
$$A(x)=A_1(x^2)+{\omega }_n{A}_2(x^2)$$
哇我真棒
但是，这还是没有用…
在仔细研究为什么需要使用复数.
因为，${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$!
所以，我们取$x={\omega }_n^{k+\frac{n}{2}}$得:
$$A(x)=A_1(x^2)-{\omega }_n{A}_2(x^2)$$
整理一下得($k<\frac{n}{2}$):
$$A({\omega }_n^k)=A_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k)$$
$$A({\omega }_n^{k+\frac{n}{2}})=A_1({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k)$$
出现了!FFT!

IFFT

看似好像不那么容易.
但是我们已经会FFT了！
在分支向下的时候，保存每一项的系数,
具体怎么保存呢?
在转换的时候，我们乘的${\omega }_n$.
转换回来的时候，就乘上它的共轭复数 What?你不会?
于是，我们就很快活

递归版FFT&IFFT

题目传送门
//为什么叫BF很快揭晓
先随便打个真·暴力上去…

#include<bits/stdc++.h>
using namespace std;
const int N=1e6+5;
long long a[N],b[N],c[N];
int n,m;
int main()
{
	//True BF
	//O(mn)...
    scanf("%d%d",&m,&n);
    for (int i=0;i<=m;i++) scanf("%d",&a[i]);
    for (int i=0;i<=n;i++) scanf("%d",&b[i]);
    for (int i=0;i<=m;i++)
        for (int j=0;j<=n;j++)
            c[i+j]+=a[i]*b[j];
    for (int i=0;i<=n+m;i++)
        printf("%d ",c[i]);
    return 0;
}

Result:

//哇竟然能有55分！
而上面的诸多公式告诉我们，FFT可以完美解决这好像就是一道FFT模板题…
根据我们的FFT算法，很容易想到递归实现FFT:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double pi=acos(-1.0); 
const int N=5e6+5;
struct complex
{
	double x,y;
	complex(double _x=0,double _y=0) {x=_x,y=_y;}
}a[N],b[N];
complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void FFT(int n,complex *a,int inv)
{
	if (n==1) return;
	int k=n>>1;
	complex a1[k],a2[k];
	for (int i=0;i<n;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];//A1和A2
	FFT(k,a1,inv);//
	FFT(k,a2,inv);//分治
	complex Wn=complex(cos(2.*pi/n),inv*sin(2.*pi/n)),w=complex(1,0);//omegaN和单位元
	for (int i=0;i<k;i++,w=w*Wn)
	{
		a[i]=a1[i]+w*a2[i];
		a[i+k]=a1[i]-w*a2[i];//一次推两个
	}
}
int main()
{
	int n,m;
	scanf("%d%d",&n,&m);
	for (int i=0;i<=n;i++) scanf("%lf",&a[i].x);
	for (int i=0;i<=m;i++) scanf("%lf",&b[i].x);
	int p=1;
	while (p<=n+m) p<<=1;
	FFT(p,a,1),FFT(p,b,1);//系数转点值
	for (int i=0;i<=p;i++) a[i]=a[i]*b[i];
	FFT(p,a,-1);//点值转系数
	for (int i=0;i<=n+m;i++) printf("%d ",(int)(abs)(a[i].x/p+0.5));
	return 0;
}

然而我竟然AC了!这是怎么回事?
Result:

这还怎么提出迭代版FFT!别提了

迭代版FFT&IFFT

接下来，我们将它优化一下:
这里是经过FFT分治后的序列

初始序列(id)	0	1	2	3	4	5	6	7
1次操作	0	2	4	6	1	3	5	7
2次操作	0	4	2	6	1	5	3	7
2进制编码	000	100	010	110	001	101	011	111
2进制逆序编码	000	001	010	011	100	101	110	111

我真聪明
可以发现，最后的编码是按逆序升序排列的.
换句话说:每个下标位置处理后的位置是它二进制逆序后的位置!

用这个短小精悍的代码就可以解决一切问题:

//名叫Rader算法
int p=1,L=0;
while (p<=n+m) p<<=1,L++;
for (int i=0;i<p;i++)
	order[i]=(order[i>>1]>>1)|((i&1)<<(L-1));
	//(i&1)<<(L-1) 是为了考虑后面的1<<(L-1)个数

真有趣 ?
接下来，我们就更快活开心了 ?

蝴蝶效应

属于小优化的范畴.
主要是因为复数乘法过于缓慢，所以可以预先处理好.
而已qwq.不要想到龙卷风那里去了

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=5e6+5;
const double pi=acos(-1.);
struct complex 
{
	double x,y;
	complex(double _x=0,double _y=0) {x=_x,y=_y;}
}a[N],b[N];
complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
int order[N];
int p=1,L;
void FFT(complex *a,int inv)
{
	for (int i=0;i<p;i++) 
		if (i<order[i]) swap(a[i],a[order[i]]);
	for (int l=1;l<p;l<<=1)
	{
		complex Wn(cos(pi/l),inv*sin(pi/l));
		for (int R=l<<1,j=0;j<p;j+=R)
		{
			complex w(1,0);
			for (int k=0;k<l;k++,w=w*Wn)
			{
				complex x=a[j+k],y=w*a[j+l+k];//B
				a[j+k]=x+y;
				a[j+l+k]=x-y;
			}
		}
	}	
} 
int main()
{
	int n,m;
	scanf("%d%d",&n,&m);
	for (int i=0;i<=n;i++) scanf("%lf",&a[i].x);
	for (int i=0;i<=m;i++) scanf("%lf",&b[i].x);
	while (p<=n+m) p<<=1,L++;
	for (int i=0;i<p;i++)
		order[i]=(order[i>>1]>>1)|((i&1)<<(L-1));
	FFT(a,1),FFT(b,1);
	for (int i=0;i<=p;i++) a[i]=a[i]*b[i];
	FFT(a,-1);
	for (int i=0;i<=n+m;i++)
		printf("%d ",(int)(a[i].x/p+0.5));
	return 0;
}

Result:

感谢奆老关注 qwq ?

后记

FFT，又称Fast-Fast-TLE…
这里为什么不直接从$<complex>$这个库里调用$complex$?
因为…太慢了…