FFT感性瞎扯
为了自己以后再用FFT时不再一脸懵X,本蒟蒻决定感性理解一下FFT。
FFT可以干啥?
把两个多项式乘在一起。具体地说,对于两个多项式$f:f(x)=\Sigma_{i=0}^n f_i x^i$和$g:g(x)=\Sigma_{i=0}^m g_i x^i$,得到一个多项式$h:h(x)=\Sigma_{x=0}^{n+m}\Sigma_{i=0}^x f(i)g(x-i)$。
显然,如果暴力的话,是$O(n^2)$的。
但是,FFT可以做到$O(nlog_2n)$!
我们一点一点把它扯个不明不白。
Part1.复数
只需要记住一点结论:
复数相乘时,模长相乘,辐角相加。
比如说,对于两个极角表示法的复数$u(u_l,u_{\theta})$与$v(v_l,v_{\theta})$,有$w=uv=(u_l v_l,u_{\theta}+v_{\theta})$
Part2.单位根
(以下,默认$n$为$2$的整数幂)
我们把单位圆上的某个复数$(cos(2k\pi/n),sin(2k\pi/n))$称作单位根$\omega_n^k$。换句话说,就是极角表示法下的复数$(1,2k\pi/n)$。当然,因为是在单位圆上,我们可以只关注它的辐角$(2k\pi/n)$。
如图。
![8次单位根](https://cdn.luogu.com.cn/upload/image_hosting/otcevc5c.png)
结论1:$\omega_n^k=\omega_n^{k+n}$。
类比正负角。可以感性理解一下,就等于绕了一整圈,度数$>2\pi$的角。有$\omega_n^k=2k\pi/n=2k\pi/n+2\pi=2k\pi/n+2n\pi/n=2(k+n)\pi/n=\omega_n^{k+n}$
结论2:$\omega_n^k=-\omega_n^{k+n/2}$
等于绕了半圈的角。
结论3:当$k$为偶数时,$\omega_n^k=\omega_{n/2}^{k/2}$。这是最重要的结论。
有$\omega_n^k=2k\pi/n=2(k/2)\pi/(n/2)=\omega_{n/2}^{k/2}$
Part3.DFT
显然,$n+1$个点可以唯一确定一个$n$次多项式。
而我们如果知道$f$和$g$上各$n+m+1$个点,就可以在$O(n)$时间內把它们乘在一起。对于每个点$(x,f(x))$与$(x,g(x))$,新点即为$(x,f(x)g(x))$。
DFT解决的就是将多项式$f$由系数表达转为点值表达。
开始推式子:
$f(x)$
$=a_0+a_1x+a_2x^2+...+a_nx^n$
$=(a_0+a_2x^2+...+a_nx^n)+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})$
$=(a_0+a_2x^2+...+a_nx^n)+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})$
令$f1(x)=a_0+a_2x+...+a_nx^{n/2}$,$f2(x)=a_1+a_3x+...+a_{n-1}x^{n/2}$
则$f(x)=f1(x^2)+xf2(x^2)$。
令$x=\omega_n^k(k<n/2)$。
则$f(\omega_n^k)$
$=f1((\omega_n^k)^2)+\omega_n^kf2((\omega_n^k)^2)$
$=f1(\omega_n^{2k})+\omega_n^kf2(\omega_n^{2k})$
$=f1(\omega_{n/2}^k)+\omega_n^kf2(\omega_{n/2}^k)$
再令$x=\omega_n^{k+n/2}(k<n/2)$。
则$f(\omega_n^{k+n/2})$
$=f1((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}f2((\omega_n^{k+n/2})^2)$
$=f1(\omega_n^{2k+n})+\omega_n^{k+n/2}f2(\omega_n^{2k+n})$
$=f1(\omega_n^{2k})+\omega_n^{k+n/2}f2(\omega_n^{2k})$
$=f1(\omega_{n/2}^k)-\omega_n^kf2(\omega_{n/2}^k)$
注意到什么了吗?
**$f(\omega_n^k)$与$f(\omega_n^{k+n/2})$,结果只有最后一项的符号不同!!!!!!**
然后就可以分治了。只要知道$f1(\omega_{n/2}^k)$与$f2(\omega_{n/2}^k)$的值,就可以直接得出$f(\omega_n^k)$与$f(\omega_n^{k+n/2})$!!!
因此我们可以令$f1(\omega_{n/2}^k)$成为新的$f(x)$,进一步递推。
Part4.IDFT
只需要将DFT中的$\omega_n^k$改为$\omega_n^{-k}$即可。别忘了最后记得除上数组长度!!!
Part5.代码:
#include<bits/stdc++.h> using namespace std; const int MAXN=4000000; int n,m,lim=1,t,res[MAXN]; const double pi=acos(-1); struct cp{ double x,y; cp(double u=0.0,double v=0.0){ x=u,y=v; } friend cp operator +(const cp &lv,const cp &rv){ return cp(lv.x+rv.x,lv.y+rv.y); } friend cp operator -(const cp &lv,const cp &rv){ return cp(lv.x-rv.x,lv.y-rv.y); } friend cp operator *(const cp &lv,const cp &rv){ return cp(lv.x*rv.x-lv.y*rv.y,lv.x*rv.y+lv.y*rv.x); } }f[MAXN],g[MAXN]; void FFT(cp *a,int sz,int inv){ if(sz==1)return; int md=sz>>1; static cp b[MAXN]; for(int i=0;i<md;i++)b[i]=a[i<<1],b[i+md]=a[(i<<1)+1]; for(int i=0;i<sz;i++)a[i]=b[i]; FFT(a,md,inv),FFT(a+md,md,inv); for(int i=0;i<md;i++){ cp x=cp(cos(2*pi*i/sz),inv*sin(2*pi*i/sz)); b[i]=a[i]+x*a[i+md],b[i+md]=a[i]-x*a[i+md]; } for(int i=0;i<sz;i++)a[i]=b[i]; } int read(){ char c=getchar(); int x=0; while(c>'9'||c<'0')c=getchar(); while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar(); return x; } int main(){ n=read(),m=read(); for(int i=0,t;i<=n;i++)f[i].x=read(); for(int i=0,t;i<=m;i++)g[i].x=read(); while(lim<=n+m)lim<<=1; FFT(f,lim,1),FFT(g,lim,1); for(int i=0;i<lim;i++)f[i]=f[i]*g[i]; FFT(f,lim,-1); for(int i=0;i<lim;i++)res[i]=(int)(f[i].x/lim+0.5); for(int i=0;i<=n+m;i++)printf("%d ",res[i]); return 0; }
Part6.非递归
在FFT时,一个位置上的数的最终位置,是把其位置的二进制表达翻转后的新位置。
即有:
$(7)_{10}=(111)_2 \rightarrow(111)_2=(7)_{10}$
$(13)_{10}=(1101)_2\rightarrow(1011)_2=(11)_{10}$
$(22)_{10}=(10110)_2\rightarrow(01101)_2=(13)_{10}$
显然,当长度不同时,一个数可能翻转到不同的位置。
而它的终点位置,可以如此递推:
for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
其中,$lim$是长度,$lg=log_2(lim)$。
Part7.最终非递归代码:
#include<bits/stdc++.h> using namespace std; const int MAXN=4000000; int n,m,lim=1,t,res[MAXN],lg,rev[MAXN]; const double pi=acos(-1); struct cp{ double x,y; cp(double u=0.0,double v=0.0){ x=u,y=v; } friend cp operator +(const cp &lv,const cp &rv){ return cp(lv.x+rv.x,lv.y+rv.y); } friend cp operator -(const cp &lv,const cp &rv){ return cp(lv.x-rv.x,lv.y-rv.y); } friend cp operator *(const cp &lv,const cp &rv){ return cp(lv.x*rv.x-lv.y*rv.y,lv.x*rv.y+lv.y*rv.x); } }f[MAXN],g[MAXN]; int read(){ char c=getchar(); int x=0; while(c>'9'||c<'0')c=getchar(); while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar(); return x; } void FFT(cp *a,int tp){ for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]); for(int md=1;md<lim;md<<=1){ cp rt=cp(cos(pi/md),tp*sin(pi/md)); for(int stp=md<<1,pos=0;pos<lim;pos+=stp){ cp w=cp(1,0); for(int i=0;i<md;i++,w=w*rt){ cp x=a[pos+i],y=w*a[pos+md+i]; a[pos+i]=x+y; a[pos+md+i]=x-y; } } } } int main(){ n=read(),m=read(); while(lim<=(n+m))lim<<=1,lg++; for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1)); for(int i=0;i<=n;i++)f[i].x=read(); for(int i=0;i<=m;i++)g[i].x=read(); FFT(f,1),FFT(g,1); for(int i=0;i<lim;i++)f[i]=f[i]*g[i]; FFT(f,-1); for(int i=0;i<lim;i++)res[i]=(int)(f[i].x/lim+0.5); for(int i=0;i<=n+m;i++)printf("%d ",res[i]); return 0; }
完结撒傅里叶~~~
![](https://cdn.luogu.com.cn/upload/image_hosting/4l5cvdy8.png)