FFT:从入门到沉迷
終于學(xué)會(huì)了FFTFFTFFT并無法自拔
啥是FFT
FFTFFTFFT指快速傅里葉變換
在OIOIOI中的應(yīng)用是O(nlogn)O(nlogn)O(nlogn)計(jì)算函數(shù)卷積
人話:多項(xiàng)式乘法
多項(xiàng)式毒瘤模板題的萬惡之源
系數(shù)表達(dá)式與點(diǎn)值表達(dá)式
系數(shù)表達(dá)式就是平常的表示方法
f(x)=∑i=0naixif(x)=\sum_{i=0}^{n}a_{i}x^if(x)=i=0∑n?ai?xi
點(diǎn)值表達(dá)式就是取函數(shù)上nnn個(gè)不同的點(diǎn)來表示
如:多項(xiàng)式f(x)=x3?2x2?x+1f(x)=x^3-2x^2-x+1f(x)=x3?2x2?x+1
點(diǎn)值表達(dá)式可以為:(?1,?1)(0,1)(1,?1)(-1,-1)(0,1)(1,-1)(?1,?1)(0,1)(1,?1)
由代數(shù)基本定理,兩者等價(jià)
初步推導(dǎo)
為什么要引入點(diǎn)值表達(dá)式?
如果相乘的多項(xiàng)式選取了相同的xxx,就可以O(n)O(n)O(n)計(jì)算乘積
現(xiàn)在的問題是:由系數(shù)轉(zhuǎn)點(diǎn)值,乘完后再轉(zhuǎn)系數(shù)
然而系數(shù)轉(zhuǎn)點(diǎn)值需要取nnn個(gè)xxx,然后求值,復(fù)雜度O(n2)O(n^2)O(n2)
更尷尬的是點(diǎn)值轉(zhuǎn)系數(shù)高斯消元O(n3)O(n^3)O(n3)
看來要在取數(shù)上做手腳了
復(fù)平面與單位根
定義i2=?1i^2=-1i2=?1,則a+bia+bia+bi表示一個(gè)復(fù)數(shù)
不需要糾結(jié)iii是個(gè)啥,只需要保留,遇到i2i^2i2就化成?1-1?1就行了
定義復(fù)平面上(a,b)(a,b)(a,b)表示a+bia+bia+bi,也可看成向量
有個(gè)奇妙的性質(zhì):復(fù)數(shù)相乘,相當(dāng)于轉(zhuǎn)角(xxx正半軸逆時(shí)針旋轉(zhuǎn))相加,長度相乘
定義nnn的單位根wnw_nwn?滿足wnn=1w_n^n=1wnn?=1
這樣的單位根有nnn個(gè),為wn0,wn1,wn2,...,wnn?1w_n^0,w_n^1,w_n^2,...,w_n^{n-1}wn0?,wn1?,wn2?,...,wnn?1?
他們將平分整個(gè)復(fù)平面
折半引理
wnk=w2n2kw_n^k=w_{2n}^{2k}wnk?=w2n2k?
之所以叫折半而不叫倍增,主要是分?jǐn)?shù)太難看所以轉(zhuǎn)了一下
畫個(gè)圖就能證明
快速傅里葉變換
為表述方便,我們把不足的nnn增加次數(shù)湊成222的整數(shù)次冪。注意nnn為長度,最高次為n?1n-1n?1
腦子一熱,不如我們把nnn個(gè)單位根作xxx吧!
開始推式子了
f(x)=∑i=0n?1aixif(x)=\sum_{i=0}^{n-1}a_ix^if(x)=i=0∑n?1?ai?xi
閑著沒事按奇偶分類
f(x)=∑i=0n2a2ix2i+∑i=0n2?1a2i+1x2i+1f(x)=\sum_{i=0}^{\frac{n}{2}}a_{2i}x^{2i}+\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^{2i+1}f(x)=i=0∑2n??a2i?x2i+i=0∑2n??1?a2i+1?x2i+1
太丑了
定義
f1(x)=∑i=0n2a2ixif_1(x)=\sum_{i=0}^{\frac{n}{2}}a_{2i}x^{i}f1?(x)=i=0∑2n??a2i?xi
f2(x)=∑i=0n2?1a2i+1xif_2(x)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^{i}f2?(x)=i=0∑2n??1?a2i+1?xi
那么有
f(x)=f1(x2)+xf2(x2)f(x)=f_1(x^2)+xf_2(x^2)f(x)=f1?(x2)+xf2?(x2)
對(duì)于k≤n2k\leq\frac{n}{2}k≤2n?,代入x=wnkx=w_n^kx=wnk?
f(wnk)=f1(wn2k)+wnkf2(wn2k)f(w_n^k)=f_1(w_n^{2k})+w_n^kf_2(w_n^{2k})f(wnk?)=f1?(wn2k?)+wnk?f2?(wn2k?)
同時(shí)
f(wnk+n2)=f1(wn2k+n)+wnk+n2f2(wn2k+n)f(w_n^{k+\frac{n}{2}})=f_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}f_2(w_n^{2k+n})f(wnk+2n??)=f1?(wn2k+n?)+wnk+2n??f2?(wn2k+n?)
顯然wnn=1,w2nn=?1w_n^n=1,w_{2n}^{n}=-1wnn?=1,w2nn?=?1
所以
f(wnk+n2)=f1(wn2k)?wnkf2(wn2k)f(w_n^{k+\frac{n}{2}})=f_1(w_n^{2k})-w_n^kf_2(w_n^{2k})f(wnk+2n??)=f1?(wn2k?)?wnk?f2?(wn2k?)
我們發(fā)現(xiàn)兩個(gè)式子結(jié)構(gòu)是相同的
根據(jù)折半引理
f1(wn2k)=f1(wn2k)f_1(w_n^{2k})=f_1(w_{\frac{n}{2}}^{k})f1?(wn2k?)=f1?(w2n?k?)
f2(wn2k)=f2(wn2k)f_2(w_n^{2k})=f_2(w_{\frac{n}{2}}^{k})f2?(wn2k?)=f2?(w2n?k?)
故
f(wnk)=f1(wn2k)+wnkf2(wn2k)f(w_n^k)=f_1(w_{\frac{n}{2}}^{k})+w_n^kf_2(w_{\frac{n}{2}}^{k})f(wnk?)=f1?(w2n?k?)+wnk?f2?(w2n?k?)
f(wnk+n2)=f1(wn2k)?wnkf2(wn2k)f(w_n^{k+\frac{n}{2}})=f_1(w_{\frac{n}{2}}^{k})-w_n^kf_2(w_{\frac{n}{2}}^{k})f(wnk+2n??)=f1?(w2n?k?)?wnk?f2?(w2n?k?)
我們發(fā)現(xiàn)規(guī)模不斷減半,然后可以愉快的遞歸了。
復(fù)雜度O(NlogN)O(NlogN)O(NlogN)
逆變換
現(xiàn)在我們轉(zhuǎn)成了點(diǎn)值,進(jìn)行了乘法,該轉(zhuǎn)回系數(shù)了。
假設(shè)原多項(xiàng)式系數(shù)是aaa,將NNN個(gè)單位根代入得到了NNN個(gè)復(fù)數(shù)bbb
即:bk=∑i=0n?1aiwnkib_k=\sum_{i=0}^{n-1}a_iw_n^{ki}bk?=∑i=0n?1?ai?wnki?
現(xiàn)在我們要求aaa
把bbb再插回去試試
我們倒著插單位根
即ckc_kck?為多項(xiàng)式∑i=0n?1biwn?ki\sum_{i=0}^{n-1}b_iw_n^{-ki}i=0∑n?1?bi?wn?ki?
代入bbb
ck=∑i=0n?1∑j=0n?1ajwnijwn?kic_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_jw_n^{ij}w_n^{-ki}ck?=i=0∑n?1?j=0∑n?1?aj?wnij?wn?ki?
把wnw_nwn?合并
ck=∑i=0n?1∑j=0n?1ajwnij?kic_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_jw_n^{ij-ki}ck?=i=0∑n?1?j=0∑n?1?aj?wnij?ki?
=∑i=0n?1∑j=0n?1ajwni(j?k)\quad=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_jw_n^{i(j-k)}=i=0∑n?1?j=0∑n?1?aj?wni(j?k)?
aja_jaj?可以放前面來
ck=∑j=0n?1aj∑i=0n?1wni(j?k)c_k=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}w_n^{i(j-k)}ck?=j=0∑n?1?aj?i=0∑n?1?wni(j?k)?
考慮S=∑i=0n?1wni(j?k)S=\sum_{i=0}^{n-1}w_n^{i(j-k)}S=i=0∑n?1?wni(j?k)?
當(dāng)j=kj=kj=k,顯然S=nS=nS=n
當(dāng)j≠kj \neq kj??=k
wnj?kS=∑i=1nwni(j?k)w_n^{j-k}S=\sum_{i=1}^{n}w_n^{i(j-k)}wnj?k?S=i=1∑n?wni(j?k)?
(wnj?k?1)S=wnn(j?k)?wnj?k=wnj?k?wnj?k=0(w_n^{j-k}-1)S=w_n^{n(j-k)}-w_n^{j-k}=w_n^{j-k}-w_n^{j-k}=0(wnj?k??1)S=wnn(j?k)??wnj?k?=wnj?k??wnj?k?=0
所以S=0S=0S=0
這說明只有j=kj=kj=k時(shí)對(duì)答案有貢獻(xiàn)
此時(shí)
ck=aknc_k=a_knck?=ak?n
ak=ckna_k=\frac{c_k}{n}ak?=nck??
FFTFFTFFT至此就結(jié)束了。
代碼:
#include <iostream> #include <cstdio> #include <cmath> #include <cctype> #define MAXN 2000005 const double PI=acos(-1.0); using namespace std; inline int read() {int ans=0,f=1;char c=getchar();while (!isdigit(c)){if (c=='-')f=-1;c=getchar();}while (isdigit(c)){ans=(ans<<3)+(ans<<1)+(c^48);c=getchar();}return f*ans; } struct complex {double x,y;complex(double x=0,double y=0):x(x),y(y){} }a[MAXN],b[MAXN]; complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);} complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);} complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} void fft(int lim,complex *a,int type) {if (lim==1)return;complex a1[lim>>1],a2[lim>>1];for (int i=0;i<=lim;i+=2)a1[i>>1]=a[i],a2[i>>1]=a[i+1];fft(lim>>1,a1,type);fft(lim>>1,a2,type);complex wn(cos(2.0*PI/lim),type*sin(2.0*PI/lim)),w(1,0);for (int i=0;i<(lim>>1);i++,w=w*wn){a[i]=a1[i]+w*a2[i];a[i+(lim>>1)]=a1[i]-w*a2[i];} } int main() {int n,m;n=read(),m=read();for (int i=0;i<=n;i++)//注意下標(biāo) a[i].x=read();for (int i=0;i<=m;i++)b[i].x=read();int lim=1;while (lim<=n+m)lim<<=1;fft(lim,a,1);fft(lim,b,1);for (int i=0;i<=lim;i++) a[i]=a[i]*b[i];fft(lim,a,-1);for (int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/lim+0.5));return 0; }迭代實(shí)現(xiàn)
什么?你沒過?
觀察一下原來的代碼 我們發(fā)現(xiàn)數(shù)組開小了
在遞歸函數(shù)中,每一層都會(huì)開一個(gè)數(shù)組,這造成效率極低。
我們觀察遞歸的全過程
(二進(jìn)制)
000 001 010 011 100 101 110 111
000 010 100 110|001 011 101 111
000 100|010 110|001 101|011|111
我們發(fā)現(xiàn),遞歸完后實(shí)際上是將二進(jìn)制翻轉(zhuǎn)
那我們可以模擬遞歸的過程,算出最后的下標(biāo),再一層一層推回去
如何計(jì)算翻轉(zhuǎn)后的數(shù)?
①將第二位及以上的翻轉(zhuǎn)并右移1位
②將第一位移到最高位
可以O(N)O(N)O(N)求出
#include <iostream> #include <cstdio> #include <cstring> #include <cctype> #include <cmath> #define MAXN 5000005 const double Pi=acos(-1.0); using namespace std; inline int read() {int ans=0,f=1;char c=getchar();while (!isdigit(c)){if (c=='-') f=-1;c=getchar();}while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();return f*ans; } struct complex{double x,y;complex(const double& x=0,const double& y=0):x(x),y(y){}}a[MAXN],b[MAXN]; inline complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);} inline complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);} inline complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} int lim=1,l,r[MAXN]; void fft(complex* a,int type) {for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);// 計(jì)算最后的序列。<是為了避免換回來,寫>也一樣for (int len=2;len<=lim;len<<=1)//合并出的長度{int mid=(len>>1);//被合并的長度complex Wn(cos(Pi/mid),type*sin(Pi/mid));//單位根for (int s=0;s<=lim;s+=len)//枚舉每一段的長度{complex w(1,0);for (int k=0;k<mid;k++,w=w*Wn){complex x=a[s+k],y=w*a[s+mid+k];//一個(gè)小優(yōu)化,避免大量復(fù)數(shù)乘法運(yùn)算a[s+k]=x+y,a[s+mid+k]=x-y;}}} } int main() {int n,m;n=read(),m=read();for (int i=0;i<=n;i++) a[i].x=read();for (int i=0;i<=m;i++) b[i].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翻轉(zhuǎn)后的數(shù)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; }凡是兩個(gè)數(shù)組一一對(duì)應(yīng)然后進(jìn)行各種運(yùn)算的,都可以考慮FFTFFTFFT
甚至可以用來匹配字符串
總結(jié)
以上是生活随笔為你收集整理的FFT:从入门到沉迷的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 《机械战警:暴戾都市》:缺陷明显的“宠粉
- 下一篇: 手游王者荣耀音乐王子高渐离