学习数学真是一件赛艇的事
FFT是我到目前OI的数学相关学过最难的了
其实理解以后发现并不是很难,只是需要的基础知识比较多
前置技能
主要包括复数相关,线性代数相关,分治基础
复数相关,复平面向量
可汗学院讲的真不错,强烈推荐看一看,内容比较全面,不过时间有点长
这里对复数做一点简单的总结
定义:
$$i^2=-1$$
这里的(i)是虚数单位(imaginary number)
对于任何数(z)都可以写成形如:
$$z=a+b*i$$
对于我们之前说的实数,放在这个形式里就是b=0
其中a称为实部,(bi)称为虚部
复数的四则运算
在这里我们只需要掌握加减乘三则运算就好
设(A=a+bi), (B=c+di)
加法:
(A+B=(a+c)+(b+d)i)
减法
(A+B=(a-c)+(b-d)i)
乘法
(A*B=(a+bi)(c+di)=ac+adi+bci+bdi^2 =ac-bd+adi+bci=(ac-bd)+(ad+bc)i)
共轭复数
若(z=a+bi,z'=a-bi)则称(z'是z的共轭复数)
复平面
首先定义一个复平面,x轴表示实部,y轴表示虚部
那么任何一个复数都可以在这个平面上以一个向量的形式表示出来
然后考虑四则运算的几何表示
加法和减法依然满足平行四边形法则
乘法根据棣莫弗定理要记住模长相乘,幅角相加
单位根
在复平面上以原点为圆心,1为半径作圆得到(单位圆)
设((omega_n)^n=1),称(omega_n)为n次单位根
在复平面上我们可以想象成用n条射线,从实轴开始,把圆均分成n部分
第一条射线与单位圆的交点所形成的向量
假设n=16,如图
((omega_n)^n)就是n个(omega_n)相乘(n-1)得到的向量,
因为要满足模长相乘,幅角相加
模长都为1不变,幅角=(frac{2pi*(n-1)}{n})
就相当于
整个平面被均分成了16份,所以旋转15次以后又回到(1,0),所以(omega_n^n=1).
对于单位根我们要知道它的几个性质
1. (omega_{2n}^{2k}=omega_n^k)
把平面分成2n格,旋转2k格=把平面分成n格,旋转k格.
2. (omega_n^{n+k}=omega_n^k)
指数函数的性质.
3. (omega_n^{frac{n}{2}+k}=-omega_n^k)
因为(omega_n^{frac{n}{2}}=-1).
4. (omega_n=cos)(frac{2pi}{n}+isin)({frac{2pi}{n}})
几何性质
线性代数相关
对于这一方面要求理解的并不是很多,只要了解矩阵乘法和矩阵的逆就好了
不了解也没有关系,接下来会详细说明
多项式乘法
例题:已知两个n次多项式(A=a_0+a_1x+a_2x^2+...+a_nx^n)和(B=b_0+b_1x+b_2x^2+...+b_nx^n),求A*B的各项系数.
首先看到这道题的做法就是(O(n^2))的暴力,将A的每个系数与B的每个系数相乘
想一想有没有可以优化的地方?
然而并没有......
接下来就需要FFT的操作了
系数表示法与点值表示法
对于一个多项式,我们最常用的把他表示出来的方法是系数表示法
就是形如(A=a_0+a_1x+a_2x^2+...+a_nx^n)的式子
其实还有另外一种表示方法点值表示法
((x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))...(x_n,f(x_n)).)
就像两点确定一条直线,三点确定一条抛物线一样,n+1个点能确定一个n次多项式
我们发现对于多项式A和B
(A=(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))...(x_n,f(x_n))).
(B=(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2))...(x_n,g(x_n))).
(A*B=(x_0,f(x_0))(*)(g(x_0)),(x_1,f(x_1))(*)(g(x_1))...(x_n,f(x_n))(*)(g(x_n)).)
这个操作是(O(n))的
这个(O(n))给我们提供了一个很好的思路,我们可以通过某种方法将系数表示变成点值表示,
再O(n)计算A*B,最后再通过某种方法将点值表示变回系数表示
一张图
考虑第一个奇怪的方法,如果采用暴力赋值计算,复杂度还是(O(n^2))的
快速幂?naive 更慢!(O(n^2logn))
所以我们不得不采用一种特殊的方法
给你一个多项式
(A(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+...+a_{n-1}x^{n-1})
我们设
(A_0(x)=a_0+a_2x+a_4x^2+a_6x^3+...+a_{n-2}x^{frac{n}{2}})
(A_1(x)=a_1+a_3x+a_5x^2+a_7x^3+...+a_{n-1}x^{frac{n}{2}})
可以发现
(A(x)=A_0(x^2)+xA_1(x^2))
然后就是一步骚操作,令(x=omega_{n}^k)
left[
egin{matrix}
a_0
a_1
a_2
vdots
a_{n-1}
end{matrix}
ight]
left[
egin{matrix}
y_0
y_1
y_2
vdots
y_{n-1}
end{matrix}
ight]
left[
egin{matrix}
(x_0)^0 & (x_0)^1 &(x_0)^2 &cdots & (x_0)^{n-1}
(x_1)^0 & (x_1)^1 &(x_1)^2 &cdots & (x_1)^{n-1}
(x_2)^0 & (x_2)^1 &(x_2)^2 &cdots & (x_2)^{n-1}
vdots & vdots &vdots & ddots & vdots
(x_{n-1})^0 & (x_{n-1})^1 &(x_{n-1})^2 &cdots & (x_{n-1})^{n-1}
end{matrix}
ight]*left[
egin{matrix}
a_0
a_1
a_2
vdots
a_{n-1}
end{matrix}
ight]=left[
egin{matrix}
y_0
y_1
y_2
vdots
y_{n-1}
end{matrix}
ight]
left[
egin{matrix}
(omega_n0)0 & (omega_n0)1 &(omega_n0)2 &cdots & (omega_n0){n-1}
(omega_n1)0 & (omega_n1)1 &(omega_n1)2 &cdots & (omega_n1){n-1}
(omega_n2)0 & (omega_n2)1 &(omega_n2)2 &cdots & (omega_n2){n-1}
vdots & vdots &vdots & ddots & vdots
(omega_n{n-1})0 & (omega_n{n-1})1 &(omega_n{n-1})2 &cdots & (omega_n{n-1}){n-1}
end{matrix}
ight]*left[
egin{matrix}
a_0
a_1
a_2
vdots
a_{n-1}
end{matrix}
ight]=left[
egin{matrix}
y_0
y_1
y_2
vdots
y_{n-1}
end{matrix}
ight]
在刚才的变换中,我们已知(B)向量,然后通过分治算法求得(C)向量.
现在我们通过(O(n))时间的乘法得到新的(C)向量,也就是新的点值表示,我们要重新求回原来的(B)向量,怎么办
考虑逆矩阵.
$$A{-1}*A*B=A{-1}*C$$
$$B=A^{-1}*C$$
我们只要求出(A)的逆矩阵就可以求出(B)向量了
(A)的逆矩阵并不好求,我们只需要知道它是什么就好了
(A)矩阵是一个特殊的范德蒙德矩阵,范德蒙德矩阵就是指每一行的元素为一个等比数列.
对于这个矩阵,我们有
$$A^{-1}=frac{1}{n}*overline{A}$$
其中(overline{A})是(A的共轭矩阵),共轭矩阵就是指矩阵内的所有元素都取共轭复数得到的矩阵.
具体证明最后再讲.
有了这个性质我们重新看看最开始的式子:
$$B=A^{-1}CB=frac{1}{n}overline{A}*C$$
是不是很神奇?原来的只要将我们分治的时候代入的(omega_n^x)换成(omega_n^{-x})就可以求出(B)向量,也就是系数表示了.
现在证明
$$A^{-1}=frac{1}{n}*overline{A}$$
设
$$frac{1}{n}overline{A}A=C$$
我们就是要证明(C=U),(U)是单位矩阵
(C_{i,j}=frac{1}{n}sum_{k=0}^{n-1}A_{i,k}*overline{A}_{k,j})
当(i=j)时
(A_{i,k}*overline{A}_{k,i}=1),因为模长相乘,幅角相加,模长为一的两个共轭复数相乘以后等于1.
(C_{i,j}=frac{1}{n}sum_{k=0}^{n-1}1=1)
当(i eq j)时
(C_{i,j}=frac{1}{n}sum_{k=0}^{n-1}A_{i,k}*overline{A}_{k,j})
(=frac{1}{n}sum_{k=0}^{n-1}(omega_n^i)^k*(omega_n^{-k})^j)
(=frac{1}{n}sum_{k=0}^{n-1}(omega_n^{i-j})^k)
(omega_n^{i-j})可以视为常数项,就变成了等比数列求和公式.
(=frac{1}{n}frac{omega_n^{i-j}[1-(omega_n^{i-j})^n]}{1-omega^{i-j}})
因为((omega_n^{i-j})^n=1,i!=j).
(C_{i,j}=0).
所以(C=U)就是我们要证的单位矩阵.
对我们刚才的那些操作取个名字,第一步系数->点值的操作叫DFT,
第三步点值->系数的操作叫IDFT
所有操作合起来叫做FFT.
最后看一看代码
洛谷上有模板题.
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define maxn 5000010
const double Pi=acos(-1.0);
inline int read(){
int ret=0,ff=1;
char ch=getchar();
while(ch<'0'||ch>'9'){
if(ch=='-') ff=-ff;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
ret=ret*10+ch-'0';
ch=getchar();
}
return ret*ff;
}
int lim=1,l=0;
int r[maxn];
struct Complex{
double x,y;
Complex(double _x=0,double _y=0){x=_x,y=_y;}
Complex operator+(Complex t){ return Complex(x+t.x,y+t.y);}
Complex operator-(Complex t){ return Complex(x-t.x,y-t.y);}
Complex operator*(Complex t){ return Complex(x*t.x-y*t.y,x*t.y+y*t.x);}
}a[maxn],b[maxn];
void fft(Complex *A,int op){
for(int i=0;i<lim;++i) if(i<r[i]) swap(A[i],A[r[i]]);
//用之前的r[i]存的顺序对A重新排列
for(int Mid=1;Mid<lim;Mid<<=1){
//披着倍增外套的分治
Complex Wn(cos(Pi/Mid),op*sin(Pi/Mid));//因为Mid是我们枚举的中点,本来就是要求区间的2倍,所以把2约掉
for(int R=Mid<<1,j=0;j<lim;j+=R){
Complex w(1,0);
for(int k=j;k<Mid+j;++k,w=w*Wn){
Complex t1=A[k],t2=w*A[k+Mid];
A[k]=t1+t2,A[k+Mid]=t1-t2;
}
}
}
}
int main(){
// freopen("mod.in","r",stdin);
int n=read(),m=read();
for(int i=0;i<=n;++i) a[i].x=read();
for(int j=0;j<=m;++j) b[j].x=read();
while(lim<=n+m) lim<<=1,++l;
for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
//这一步就是二进制的转置部分,r[i]表示i的二进制转置,比如说r[6(110)]=3(011)
//r[i]由r[i/2]递推得来,对比i和i/2的二进制规律,我们发现i=(i>>2)<<1+(i&1)
//因为r[i]是i的倒序,所以也应该是倒序递推
fft(a,1),fft(b,1);
for(int i=0;i<lim;++i) a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<=n+m;++i) printf("%d ",int(a[i].x/lim+0.5));
return 0;
}