【学习笔记】快速傅里叶变换(FFT)
【学习笔记】快速傅里叶变换
学习之前先看懂这个
浅谈范德蒙德(Vandermonde)方阵的逆矩阵的求法以及快速傅里叶变换(FFT)中IDFT的原理——gzy
hhh开个玩笑。
讲一下\(FFT\) 的流程，我也不准备长篇大论地分析\(FFT...\)
- 将系数表示法转换为点值表示法 \(O(n \log n)\)
- 对于点值表示法直接进行操作 \(O(n)\)
- 将点值表示法转换为系数表示法 \(O(n \log n)\)
这样的流程，最终复杂度是\(O(n \log n)\) 的，现在我们从最简单的点值表示法讲起。
对于点值表示法直接进行操作 :
- 一个\(n\)次多项式可以被\(n+1\)个点唯一确定。
证明：显然。设\(f(x)=\Sigma_i^{n+1} a_ix^{i-1}\)，直接通过待定系数法解即可。
- 将两个多项式相加，(从函数的角度上理解)，就是\(x\)不变，\(y\)相加，减法同理
- 将两个多项式相乘，(从函数的角度上理解)，就是\(x\)不变，\(y\)相乘。除法求逆
那么假如我们已经有了多项式的点值表示法了，怎么办乘起来呢?
for(register int t=0;t< len;++t) A[t]*=B[t];
就是这样简单。
将系数表示法转换为点值表示法 \((DFT)\)
朴素的想法是带入\(O(n^2)\)求值，还是用了秦六韶的情况下...
为什么我们不能考虑优化，是因为不同的\(x\)的值是不一样的(不然呢)
但是考虑一种特殊情况：
是不是就一样了这为我们优化复杂度创造了条件，引入复数的概念。利用复数中的"单位负根"的性质我们可以直接分治求值。
详细证明我最后放连接。
将点值表示法转换为系数表示法\((IDFT)\)
还记得开头提到的范德蒙德(Vandermonde)方阵吗...
如果不知道请转到《工程数学——线性代数(第六版)》高等教育出版社第\(18\)面，或者浅谈范德蒙德(Vandermonde)方阵的逆矩阵的求法以及快速傅里叶变换(FFT)中IDFT的原理——gzy。
考虑拿个系数右乘这个矩阵
实际上就是把\(x_0->x_{n-1}\)这么多变量带入式子。
范德蒙德这个矩阵是可逆的，于是我们有
假如我们可以快速求出范德蒙德矩阵的逆矩阵，然后顺便对于这个特殊的矩阵乘法快速求值，这样我们就可以得到原来的系数是吧。
同样利用上面\(IDT\)的方法，可以做到\(n\log n\)，甚至只要在某个地方把\(1->-1\)，数学是不是很神奇？所以最后的代码很短，大概是这样的：
FFT(A,1);
FFT(B,1);
for(register int t=0;t<=len;++t) A[t]*=B[t];
FFT(A,-1);
是不是很神奇？
我在这里讲得很省略，主要是为了给自己梳理思路(方便记忆)，\(FFT\)有许多有趣的数学推理，大家可以自己去查阅资料，限于篇幅和笔者水平这里就略去了咕咕咕。
注意！如果你就看了我的解析你绝对看不懂下面的代码，除非你是神仙(那你还是看吧,这位神仙)
#include<bits/stdc++.h>
using namespace std; typedef long long ll;
template < class ccf > inline ccf qr(ccf ret){ ret=0;
register char c=getchar();
while(c<48||c>57) c=getchar();
while(c>=48&&c<=57) ret=ret*10+c-48,c=getchar();
return ret;
}
inline int qr(){return qr(1);}
const int maxn=1<<21|1;
const double pi=acos(-1.0);
const double eps=1e-5;
int n,m,len(1);
typedef complex<double> is;
is A[maxn];
is B[maxn];
int r[maxn];
inline void FFT(is*a,double ty){
is w,temp,wn,*a1,*a0;
for(register int t=0;t< len;++t)
if(t<r[t]) swap(a[t],a[r[t]]);
for(register int t=1;t<len;t<<=1){
wn=is(cos(pi/t),sin(pi/t)*ty);
for(register int j=0;j<len;j+=t<<1){
a1=t+(a0=a+j); w=is(1,0);
for(register int k=0;k<t;++k,++a0,++a1,w*=wn)
temp=*a1*w, *a1=*a0-temp, *a0=*a0+temp;
}
}
}
int main(){
#ifndef ONLINE_JUDGE
freopen("in.in","r",stdin);
freopen("out.out","w",stdout);
#endif
n=qr();m=qr();
int cnt=0;
while(len<=n+m)len<<=1,++cnt;
for(register int t=0;t<=n;++t) A[t]=qr();
for(register int t=0;t<=m;++t) B[t]=qr();
for(register int t=0;t< len;++t) r[t]=(r[t>>1]>>1)|((t&1)<<(cnt-1));
FFT(A,1); FFT(B,1);
for(register int t=0;t< len;++t) A[t]*=B[t];
FFT(A,-1);
for(register int t=0;t<=m+n;++t) printf("%.0lf%c",fabs((A[t].real())/len+eps),t==m+n?'\n':' ');
return 0;
}
参考资料
假如你不想看那些具体而不够抽象也不够严谨的讲解并且想一次完全明白FFT(数学基础要好)
