前言
这篇博客很长,因为考虑到读者很多都是对FFT的了解除了名字再无其它的大白,所以(并不)详细地解释了许多,只要你耐心地把它看完,你一定会对FFT有更深刻的理解。
准备工作
首先我们要从泰勒公式说起。
泰勒公式是指若一个函数在
x=x0
x
=
x
0
处有n阶导数,我们可以用含这n阶导数和
(x−x0)
(
x
−
x
0
)
的多项式来模拟函数在
(x−Δx,x+Δx)(其中Δx无限接近于0)
(
x
−
Δ
x
,
x
+
Δ
x
)
(
其
中
Δ
x
无
限
接
近
于
0
)
的取值。用公式表示就是:
(此图来源于百度百科)
然后,我们可以根据泰勒公式推导出欧拉公式,即
(证明请自行查阅百度百科)
现在我们来考虑方程
xn=1
x
n
=
1
的解,由此引出n次单位根的概念,即
e2πi/n
e
2
π
i
/
n
,单位根的0——n-1次方都是这个方程的解
多项式
多项式有两种常用的表示方法,系数表示法和点值表示法,系数表示法就是我们常用的
y=kx+b
y
=
k
x
+
b
,点值表示法是说我们可以用至少n个点来表示n-1次多项式(求出各个系数)。
而这两种表示法是可以互相转化的,从系数到点值我们称为离散傅里叶变换(DFT),从点值到系数我们称为逆离散傅里叶变换(IDFT)。
快速傅里叶变换(FFT)
结合上面两点,我们取单位根的0——n-1次方为取样点,即可在
nlogn
n
l
o
g
n
的时间内快速完成离散傅里叶变换。(在操作之前我们需要先将阶数补为2的t次方)
对于一个多项式:
A(x)=a0+a1∗x+a2∗x2+...+an−1∗xn
A
(
x
)
=
a
0
+
a
1
∗
x
+
a
2
∗
x
2
+
.
.
.
+
a
n
−
1
∗
x
n
我们取
e2πi/n
e
2
π
i
/
n
的0到n-1次方为取样点,用A[k]表示第k个取样点的结果
然后我们采取分支的思路,将A分为O和E
O(x)=a1+a3∗x+a5∗x2+...+an−1∗xn/2−1
O
(
x
)
=
a
1
+
a
3
∗
x
+
a
5
∗
x
2
+
.
.
.
+
a
n
−
1
∗
x
n
/
2
−
1
E(x)=a0+a2∗x+a4∗x2+...+an−2∗xn/2−1
E
(
x
)
=
a
0
+
a
2
∗
x
+
a
4
∗
x
2
+
.
.
.
+
a
n
−
2
∗
x
n
/
2
−
1
我们会发现
A(x)=x∗O(x2)+E(x2)
A
(
x
)
=
x
∗
O
(
x
2
)
+
E
(
x
2
)
即
A[k]=x∗O[k]+E[k]
A
[
k
]
=
x
∗
O
[
k
]
+
E
[
k
]
(因为A[k]代表
x
x
取,O[k],E[k]代表
x2
x
2
取
e2π(k−1)/(n/2)=e4π(k−1)/n
e
2
π
(
k
−
1
)
/
(
n
/
2
)
=
e
4
π
(
k
−
1
)
/
n
,两者刚好相同)
这样我们就成功地把问题的规模缩小到了1/2,即在
nlogn
n
l
o
g
n
的时间内完成了DFT。
在DFT的过程中,我们相当于给初始系数矩阵A乘上了一个n*n的矩阵,如下图
那么我们做IDFT的话,只需要求这个矩阵的逆矩阵即可。
可以证明,对于该矩阵每一项取倒数再除以n就是该矩阵的逆矩阵。
这样时间复杂度就和DFT完全一样,是
nlogn
n
l
o
g
n
至此我们完成了FFT
FFT的写法
从上面的叙述中,FFT显然可以通过递归来实现分治,但这样会花费大量空间用于创建和维护数组,而我们可以使用迭代法(循环)避免这一部分的空间使用。(NOI实际评测时理论上会开无限栈,即栈的大小算在总内存里面,所以递归写法应该也能过)
我们以0——7为例来研究这样分治的特点
0,1,2,3,4,5,6,7
(0,2,4,6)(1,3,5,7)
(0,4)(2,6)(1,5)(3,7)
0,4,2,6,1,5,3,7
写成二进制:000,100,010,110,001,101,011,111
如果每个数都从右往左读的话,我们就会(惊讶地)发现,这是按从小往大的顺序排序的,现在我们只要能将数按照这种方式排序,就能够用循环来写FFT了。
于是我们递推预处理出rev数组
void get_rev(int bit)
{
for(int i=0;i<(1<<bit);i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
然后根据rev对a重新排序,就可以进行迭代FFT了,一个元素的DFT为其本身,两两合并得到两个元素的DFT,两两合并得到四个元素的DFT,进行下去直到得到整串的DFT。
例题
洛谷板子题
两多项式相乘,先各做一遍DFT,每个点值相乘,再做一遍IDFT即可。
代码:
#include <bits/stdc++.h>
using namespace std;
typedef complex<double> cpx;
const double pi=acos(-1);
cpx a[3000000],b[3000000];
int rev[3000000];
int n,m;
void get_rev(int bit)
{
for(int i=0;i<(1<<bit);i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void fft(cpx *a,int len,int dft)
{
for(int i=0;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int step=1;step<len;step<<=1)
{
cpx wn=exp(cpx(0,pi*dft/step));
for(int j=0;j<len;j+=(step<<1))
{
cpx wnk(1,0);
for(int k=j;k<j+step;k++)
{
cpx x=a[k];
cpx y=wnk*a[k+step];
a[k]=x+y;
a[k+step]=x-y;
wnk*=wn;
}
}
}
if(dft==-1) for(int i=0;i<len;i++) a[i]/=len;
}
int main()
{
scanf("%d%d",&n,&m);
n++;m++;
for(int i=0;i<n;i++)
{
double kk;
scanf("%lf",&kk);
a[i]=cpx(kk,0);
}
for(int i=0;i<m;i++)
{
double kk;
scanf("%lf",&kk);
b[i]=cpx(kk,0);
}
int bit=1,s=2;
for(;(1<<bit)<n+m;bit++) s<<=1;
get_rev(bit);
fft(a,s,1);
fft(b,s,1);
for(int i=0;i<s;i++) a[i]*=b[i];
fft(a,s,-1);
for(int i=0;i<n+m-1;i++) printf("%d ",(int)(a[i].real()+0.5)); printf("\n");
return 0;
}