学习数学真是一件赛艇的事
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{2\pi*(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{2\pi}{n}+i\sin\)\({\frac{2\pi}{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\)
\[A(\omega_{n}^k)=A_0(\omega_{n}^{2k})+\omega_{n}^kA_1(\omega_{n}^{2k})\]\[=A_0(\omega_{\frac{n}{2}}^{k})+\omega_{n}^kA_1(\omega_{\frac{n}{2}}^{k})\]又令\(x=\omega_n^{\frac{n}{2}+k}\)\[A(\omega_n^{\frac{n}{2}+k})=A_0(\omega_{n}^{n+2k})+\omega_n^{\frac{n}{2}+k}A_1(\omega_{n}^{n+2k})\]\[=A_0(\omega_{n}^{2k})-\omega_{n}^kA_1(\omega_{n}^{2k})\]\[=A_0(\omega_{\frac{n}{2}}^{k})-\omega_{n}^kA_1(\omega_{\frac{n}{2}}^{k})\]比较上面两个式子,我们发现
假设我们已经知道了\(A_0(\omega_{\frac{n}{2}}^{k})\)和\(A_1(\omega_{\frac{n}{2}}^{k})\),我们就能同时知道\(A(\omega_{n}^k)\)和\(\omega_n^{\frac{n}{2}+k}\).
这样就把问题缩小了一半,对于每个问题都缩小一半,复杂度就从\(O(n^2)\)变成了\(O(nlogn)\).
就可以利用一种类似于线段树的操作来计算,如果没有理解就看看这张盗来的图
每一层向上转移是\(O(n)的\),因为树高只有\(logn\)层,总复杂度就变得很小
此时我们的第一步由系数到点值已经结束了,不过还要补充一点
观察树的最底层
序列为\(0,4,2,6,1,5,3,7\)
原序列\(0,1,2,3,4,5,6,7\)
转化成二进制来发现规律
新序列\(000,100,010,110,101,011,111\)
原序列\(000,001,010,011,101,110,111\)
我们发现新序列的每个数是原序列的二进制反转
现在我们已经知道了每个序列的下标,我们就可以利用下标实现
这种实现方法更像是倍增,而不是分治,虽然算法的本质还是分治
如果要看代码在最后面...
IDFT
对于把点值表示转化成系数表示,我们已经知道了用分治的方法快速求.
抛开分治以及log算法不谈,我们可以从另一方面理解刚才的操作
线性代数
把所有系数放在一起组成一个系数列向量:
\[
\left[
\begin{matrix}
a_0 \\
a_1 \\
a_2 \\
\vdots \\
a_{n-1} \\
\end{matrix}
\right]
\]
最终的n个点值也可以放在一起组成一个点值列向量
\[
\left[
\begin{matrix}
y_0 \\
y_1 \\
y_2 \\
\vdots \\
y_{n-1} \\
\end{matrix}
\right]
\]
系数列向量可以通过左乘一个矩阵得到点值列向量
\[
\left[
\begin{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}
\right]*\left[
\begin{matrix}
a_0 \\
a_1 \\
a_2 \\
\vdots \\
a_{n-1} \\
\end{matrix}
\right]=\left[
\begin{matrix}
y_0 \\
y_1 \\
y_2 \\
\vdots \\
y_{n-1} \\
\end{matrix}
\right]
\]
因为我们代入的\(x_0,x_1,x_2...x_{n-1}\)分别是\(\omega_n^0,\omega_n^1\omega_n^2...\omega_n^{n-1}\).
矩阵就变成了\[
\left[
\begin{matrix}
(\omega_n^0)^0 & (\omega_n^0)^1 &(\omega_n^0)^2 &\cdots & (\omega_n^0)^{n-1} \\
(\omega_n^1)^0 & (\omega_n^1)^1 &(\omega_n^1)^2 &\cdots & (\omega_n^1)^{n-1} \\
(\omega_n^2)^0 & (\omega_n^2)^1 &(\omega_n^2)^2 &\cdots & (\omega_n^2)^{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}
\right]*\left[
\begin{matrix}
a_0 \\
a_1 \\
a_2 \\
\vdots \\
a_{n-1} \\
\end{matrix}
\right]=\left[
\begin{matrix}
y_0 \\
y_1 \\
y_2 \\
\vdots \\
y_{n-1} \\
\end{matrix}
\right]
\]
设这个矩阵为
\[A*B=C\]
在刚才的变换中,我们已知\(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}*C\\B=\frac{1}{n}*\overline{A}*C\]
\[
\left[
\begin{matrix}
a_0 \\
a_1 \\
a_2 \\
\vdots \\
a_{n-1} \\
\end{matrix}
\right]= \frac{1}{n} \left[
\begin{matrix}
(\omega_n^0)^0 & (\omega_n^0)^1 &(\omega_n^0)^2 &\cdots & (\omega_n^0)^{n-1} \\
(\omega_n^{-1})^0 & (\omega_n^{-1})^1 &(\omega_n^{-1})^2 &\cdots & (\omega_n^{-1})^{n-1} \\
(\omega_n^{-2})^0 & (\omega_n^{-2})^1 &(\omega_n^{-2})^2 &\cdots & (\omega_n^{-2})^{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}
\right]*\left[
\begin{matrix}
y_0 \\
y_1 \\
y_2 \\
\vdots \\
y_{n-1} \\
\end{matrix}
\right]
\]
是不是很神奇?原来的只要将我们分治的时候代入的\(\omega_n^x\)换成\(\omega_n^{-x}\)就可以求出\(B\)向量,也就是系数表示了.
现在证明
\[A^{-1}=\frac{1}{n}*\overline{A}\]
设
\[\frac{1}{n}*\overline{A}*A=C\]
我们就是要证明\(C=U\),\(U\)是单位矩阵
\[\frac{1}{n}\left[
\begin{matrix}
(\omega_n^0)^0 & (\omega_n^0)^1 &(\omega_n^0)^2 &\cdots & (\omega_n^0)^{n-1} \\
(\omega_n^1)^0 & (\omega_n^1)^1 &(\omega_n^1)^2 &\cdots & (\omega_n^1)^{n-1} \\
(\omega_n^2)^0 & (\omega_n^2)^1 &(\omega_n^2)^2 &\cdots & (\omega_n^2)^{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}
\right]\left[
\begin{matrix}
(\omega_n^0)^0 & (\omega_n^0)^1 &(\omega_n^0)^2 &\cdots & (\omega_n^0)^{n-1} \\
(\omega_n^{-1})^0 & (\omega_n^{-1})^1 &(\omega_n^{-1})^2 &\cdots & (\omega_n^{-1})^{n-1} \\
(\omega_n^{-2})^0 & (\omega_n^{-2})^1 &(\omega_n^{-2})^2 &\cdots & (\omega_n^{-2})^{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}
\right]=\left[
\begin{matrix}
1 & 0 &0 &\cdots & 0 \\
0 & 1 &0 &\cdots & 0 \\
0 & 0 &1 &\cdots & 0 \\
\vdots & \vdots &\vdots & \ddots & \vdots \\
0 & 0 &0 &\cdots & 1 \\
\end{matrix}
\right]\]
\(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\neq 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;
}