終于學(xué)會(huì)了$FFT$并無(wú)法自拔

啥是FFT

$FFT$指快速傅里葉變換

在$OI$中的應(yīng)用是$O(nlogn)$計(jì)算函數(shù)卷積

人話：多項(xiàng)式乘法

多項(xiàng)式毒瘤模板題的萬(wàn)惡之源

系數(shù)表達(dá)式與點(diǎn)值表達(dá)式

系數(shù)表達(dá)式就是平常的表示方法

$$f(x)=\sum _{i=0}^n{a}_i{x}^i$$

點(diǎn)值表達(dá)式就是取函數(shù)上$n$個(gè)不同的點(diǎn)來(lái)表示

如：多項(xiàng)式$f(x)=x^3?2x^2?x+1$

點(diǎn)值表達(dá)式可以為：$(?1,?1)(0,1)(1,?1)$

由代數(shù)基本定理，兩者等價(jià)

初步推導(dǎo)

為什么要引入點(diǎn)值表達(dá)式？

如果相乘的多項(xiàng)式選取了相同的$x$，就可以$O(n)$計(jì)算乘積

現(xiàn)在的問(wèn)題是：由系數(shù)轉(zhuǎn)點(diǎn)值，乘完后再轉(zhuǎn)系數(shù)

然而系數(shù)轉(zhuǎn)點(diǎn)值需要取$n$個(gè)$x$,然后求值，復(fù)雜度$O(n^2)$

更尷尬的是點(diǎn)值轉(zhuǎn)系數(shù)高斯消元$O(n^3)$

看來(lái)要在取數(shù)上做手腳了

復(fù)平面與單位根

定義$i^2=?1$，則$a+bi$表示一個(gè)復(fù)數(shù)

不需要糾結(jié)$i$是個(gè)啥，只需要保留，遇到$i^2$就化成$?1$就行了

定義復(fù)平面上$(a,b)$表示$a+bi$,也可看成向量

有個(gè)奇妙的性質(zhì)：復(fù)數(shù)相乘，相當(dāng)于轉(zhuǎn)角（$x$正半軸逆時(shí)針旋轉(zhuǎn)）相加，長(zhǎng)度相乘

定義$n$的單位根$w_n$滿足$w_n^n=1$

這樣的單位根有$n$個(gè)，為$w_n^0,w_n^1,w_n^2,...,w_n^{n?1}$

他們將平分整個(gè)復(fù)平面

折半引理

$w_n^k=w_{2n}^{2k}$

之所以叫折半而不叫倍增，主要是分?jǐn)?shù)太難看所以轉(zhuǎn)了一下

畫個(gè)圖就能證明

快速傅里葉變換

為表述方便，我們把不足的$n$增加次數(shù)湊成$2$的整數(shù)次冪。注意$n$為長(zhǎng)度，最高次為$n?1$

腦子一熱，不如我們把$n$個(gè)單位根作$x$吧！

開始推式子了

$$f(x)=\sum _{i=0}^{n?1}{a}_i{x}^i$$

閑著沒(méi)事按奇偶分類

$$f(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_1(x)=\sum _{i=0}^{\frac{n}{2}}{a}_{2i}{x}^i$$

$$f_2(x)=\sum _{i=0}^{\frac{n}{2}?1}{a}_{2i+1}{x}^i$$

那么有

$$f(x)=f_1(x^2)+x{f}_2(x^2)$$

對(duì)于$k\le \frac{n}{2}$,代入$x=w_n^k$

$$f(w_n^k)=f_1(w_n^{2k})+w_n^k{f}_2(w_n^{2k})$$

同時(shí)

$$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})$$

顯然$w_n^n=1,w_{2n}^n=?1$

所以

$$f(w_n^{k+\frac{n}{2}})=f_1(w_n^{2k})?w_n^k{f}_2(w_n^{2k})$$

我們發(fā)現(xiàn)兩個(gè)式子結(jié)構(gòu)是相同的

根據(jù)折半引理

$$f_1(w_n^{2k})=f_1(w_{\frac{n}{2}}^k)$$

$$f_2(w_n^{2k})=f_2(w_{\frac{n}{2}}^k)$$

故

$$f(w_n^k)=f_1(w_{\frac{n}{2}}^k)+w_n^k{f}_2(w_{\frac{n}{2}}^k)$$

$$f(w_n^{k+\frac{n}{2}})=f_1(w_{\frac{n}{2}}^k)?w_n^k{f}_2(w_{\frac{n}{2}}^k)$$

我們發(fā)現(xiàn)規(guī)模不斷減半，然后可以愉快的遞歸了。

復(fù)雜度$O(NlogN)$

逆變換

現(xiàn)在我們轉(zhuǎn)成了點(diǎn)值，進(jìn)行了乘法，該轉(zhuǎn)回系數(shù)了。

假設(shè)原多項(xiàng)式系數(shù)是$a$，將$N$個(gè)單位根代入得到了$N$個(gè)復(fù)數(shù)$b$

即：$b_k={\sum }_{i=0}^{n?1}{a}_i{w}_n^{ki}$

現(xiàn)在我們要求$a$

把$b$再插回去試試

我們倒著插單位根

即$c_k$為多項(xiàng)式$$\sum _{i=0}^{n?1}{b}_i{w}_n^{?ki}$$

代入$b$

$$c_k=\sum _{i=0}^{n?1}\sum _{j=0}^{n?1}{a}_j{w}_n^{ij}{w}_n^{?ki}$$

把$w_n$合并

$$c_k=\sum _{i=0}^{n?1}\sum _{j=0}^{n?1}{a}_j{w}_n^{ij?ki}$$

$$=\sum _{i=0}^{n?1}\sum _{j=0}^{n?1}{a}_j{w}_n^{i(j?k)}$$

$a_j$可以放前面來(lái)

$$c_k=\sum _{j=0}^{n?1}{a}_j\sum _{i=0}^{n?1}{w}_n^{i(j?k)}$$

考慮$$S=\sum _{i=0}^{n?1}{w}_n^{i(j?k)}$$

當(dāng)$j=k$,顯然$S=n$

當(dāng)$j?=k$

$$w_n^{j?k}S=\sum _{i=1}^n{w}_n^{i(j?k)}$$

$$(w_n^{j?k}?1)S=w_n^{n(j?k)}?w_n^{j?k}=w_n^{j?k}?w_n^{j?k}=0$$

所以$S=0$

這說(shuō)明只有$j=k$時(shí)對(duì)答案有貢獻(xiàn)

此時(shí)

$$c_k=a_{k}n$$

$$a_k=\frac{c_k}{n}$$

$FFT$至此就結(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)

什么？你沒(méi)過(guò)？

觀察一下原來(lái)的代碼 我們發(fā)現(xiàn)數(shù)組開小了

在遞歸函數(shù)中，每一層都會(huì)開一個(gè)數(shù)組，這造成效率極低。

我們觀察遞歸的全過(guò)程

（二進(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)

那我們可以模擬遞歸的過(guò)程，算出最后的下標(biāo)，再一層一層推回去

如何計(jì)算翻轉(zhuǎn)后的數(shù)？

①將第二位及以上的翻轉(zhuǎn)并右移1位

②將第一位移到最高位

可以$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ì)算最后的序列。<是為了避免換回來(lái)，寫>也一樣for (int len=2;len<=lim;len<<=1)//合并出的長(zhǎng)度{int mid=(len>>1);//被合并的長(zhǎng)度complex Wn(cos(Pi/mid),type*sin(Pi/mid));//單位根for (int s=0;s<=lim;s+=len)//枚舉每一段的長(zhǎng)度{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)算的，都可以考慮$FFT$

甚至可以用來(lái)匹配字符串

總結(jié)

以上是生活随笔為你收集整理的FFT:从入门到沉迷的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。