寫在前邊

　　1.辣雞RRRR_wys之前csdn的博客，千年不更。。。還很水。。。于是開了這個Blog。。。妄圖拯救一下自己

　　2.最近接觸了一些多項式理論。于是翹掉了愉快的高頻自控，通過《算導》稍稍學習了一下

　　3.算法競賽中，FFT主要解決多項式的乘法等問題

FFT基本概念

　　1.FFT即快速傅里葉變換, 是離散傅里葉變換的加速算法。可以在o(nlogn)的時間內，完成DFT和DFT^-1

　　2.DFT即離散傅里葉變換, 這里主要就是多項式的系數向量轉換成點值表示的過程

多項式的表示

　　1.系數表達：對于一個次數上界為n的多項式?A(x)=∑_{(j=0 - n-1)} a_j*x^j?而言，系數表達式便是由系數組成的向量a = (a₀, a₁, a₂, a₃, ... , a_n-1)

? ? ? ? ?(1)計算A(x₀)的值即：A(x₀) = a₀+x₀(a₁+x₀(a₂+x₀(a₃+x₀(a₄+...+x₀(a_n-2+x₀*a_n-1))...))) ?時間為o(n)

　　 ?(2)通過系數向量計算兩個多項式的乘法：設輸出向量為c, 而a,b長度n=max(na,nb)那么顯然兩個c的長度為2n-1，計算方法為：c_j =?∑_{(k=0 ~ j)}?a_k*b_j-k，

　　 ? ? ?即向量a,b的卷積

　　2.點值表達式：對于一個次數上界為n的多項式?A(x)=∑_{(j=0 ~ n-1)}?a_j*x^j?而言，點值表達式便是由n個點值所構成的集合{ (x₀,y₀),?(x₁,y₁),?(x₂,y₂), ... , (x_n-1,y_n-1)?}

　　 ? ?(1)運用之前提到的方法，一個個計算點值時間為o(n²)，之后可以看到如果我們巧妙地選取點x_{k可以使時間變為o(nlogn)}

　　 ? ?(2)將點值表達式轉換為系數表達式：將已知的點值帶入，多項式可已得到一個有n個未知數的n個線性方程

　　　　?a₀ + a₁*x₀ + a₂*x₀² + a₃* x₀³ + ... + a_n-1*x₀^n-1 = y_0?(1)

　　　　 a₀ + a₁*x₁ + a₂*x₁² + a₃* x₁³ + ... + a_n-1*x₁^n-1 = y_1?(2)

　　　　 .......

　　　　 a₀?+ a₁*x_n-1?+ a₂*x_n-1²?+ a₃* x_n-1³?+ ... + a_n-1*x_n-1^n-1?= y_n-1?(n)

?　　　　 可以寫成矩陣形式，既為一個范德蒙矩陣，通過線代知識解這個方程就可以求出系數向量，復雜度o(n³)。另一種方法是采用拉格朗日插值法求解，復雜度o(n²)。因此，n個點的求值運算和插值運算，是定義完備的互逆運算。

　　　 (3)通過點值計算兩個多項式的乘積，c =?{ (x₀,y_A0*y_B0),?(x₁,y_A1*y_B1),?(x₂,y_A2*y_B2), ... , (x_n-1,y_A(n-1)*y_B(n-1))?}, 可以看到復雜度為o(n)

　　 ??(4)計算多項式乘法的過程：1）o(nlogn)分別求出兩個多項式的點值表達式 2）點值乘法o(n) 3）插值o(nlogn) 4）求出答案多項式的系數表達式

DFT與FFT

　　下面開始討論，如何在o(nlogn)的時間內，完成求值和插值運算

　　1.單位復數根

??? 　　n次單位復數根滿足ωⁿ=1的復數ω，n此單位復數根的數目恰好有n個，這些根是：ω_n^k?= e^2πik/n 這n個根均勻的分布在復平面上，其他n次單位根都是ω_n的冪次。n個n次單位根在乘法意義下形成一個群，該群與模n意義下的加法群，有相同的結構，即ω_nⁱ ω_n^j = ω_n^{(i + j) mod n}

　　　(1)?消去引理：ω_nd^kd =ω_n^{k ?}(n>=0, k>=0, d>0)

　　　(2)?ω_n^n/2 =ω₂ =-1

　　　(3)?折半引理：如果n>0為偶數，那么n個n次單位復數根的平方的集合就是n/2個n/2次單位復數根的集合

　　　(4)?求和引理：∑_{(j=0 ~ n-1)} (ω_n^k)^j = 0 (n>=1, 不能被n整除的非負整數k)

　　2.DFT

　　我們希望計算ω_n⁰, ω_n¹, ω_n², .... , ω_n^n-1處的值（在多項式乘法中其實是2n個點值）。假設給定系數a=(a₀, a₁, a₂, ... , a_n-1), 即對于k=0, 1, 2, 3, ... , n-1求出y_k?= A(ω_n^k)，記向量y=(y₁, y₂, y₃, ... ,y_n-1)就是系數向量a的離散傅里葉變換DFT，記為 y?= DFT_n(a)

　　3. FFT

　　通過快速傅里葉變換FFT，利用單位復數根的性質，我們可以在o(nlogn)的時間內計算DFT_n(a)。首先假設次數n恰好是2的整數冪，不足在高次位置項添0。

　　定義兩個新的次數界為n/2的多項式?A_[0]?= a₀?+ a₂*x?+ a₄*x²?+ a₆* x³?+ ... + a_n-2*x^n/2-1 , A_[1]?= a₁?+ a₃*x?+ a₅*x²?+ a₇* x³?+ ... + a_n-1*x^n/2-1

　　A_[0]包含所有偶數下標的系數，A_[1]包含所有奇數下標的系數，于是有A(x) = A_[0](x2) + x*A_[1](x²)?。

　　所以，求解A(x)的DFN，轉換為了求解多項式A_[0], A_[1]在(ω_n⁰)², (ω_n¹)², (ω_n²)²?, ... , (ω_n^n-1)²?處的取值。

　　根據折半引理，(ω_n⁰)², (ω_n¹)², (ω_n²)²?, ... , (ω_n^n-1)²?并不是由n個不同的值組成，而是由n/2個n/2次單位根復數根所組成，每個根正好出現兩次。因此，我們遞歸的對次數界為n/2的多項式A_[0], A_[1]在n/2個n/2次單位復數根處進行求值。這樣我們求出了次數界為n的次多項式在n次單位復數根處的值。

　　注意到除了遞歸調用的時間外，每次所需時間為o(n)，其中n是輸入向量的長度。因此有如下遞歸式：T(n) = 2T(n/2) + o(n) = o(nlogn)，時間為o(nlogn)。

　　計算在單位復數根處插值，求解系數向量。根據點值表達式，我們可以列出線性方程組：

　　a₀?+ a₁?+ a₂?+ a₃+ ... + a_n-1 = y0 (1)

　　a₀?+ a₁* (ω_n)¹+ a₂*(ω_n)²?+ a₃* (ω_n)³?+ ... + a_n-1*(ω_n)^n-1 = y_1?(2)

　　a₀?+ a₁*?(ω_n)²+ a₂*(ω_n)⁴?+ a₃*?(ω_n)⁶?+ ... + a_n-1*(ω_n)^2(n-1)? = ?y_{2 ?}(3)

　　......

　　a₀?+ a₁*?(ω_n)^n-1+ a₂*(ω_n)^2(n-1)?+ a₃*?(ω_n)^3(n-1)?+ ... + a_n-1*(ω_n)^(n-1)(n-1)? = ?y_n-1?(4)　

　　根據這個線性方程組，我們可以令單位復數根組成的系數矩陣為V_n，運用yV_n^-1求解系數。可以證明：V_n^-1的(j, k)處的值為(ω_n)^-kj/n ?(j,k = 0, 1, 2, ... , n-1)

　　有了上面的式子我們能推導出DFT^-1：a_j =(1/n)*?∑_{(k=0 ~ n-1)} y_k(ω_n)^{-kj ?}(j = 0, 1, 2, ... , n-1)

　　觀察可得，相對于DFT，我們把a與y替換，用ω_n^-1代替ω_n，并將每個元素除以n。這樣，我們也可以在o(nlogn)的時間內計算出DFT^-1

　　卷積定理：a(卷積)b = DFT_2n^-1( DFT_2n(a)*DFT_2n(b) ) , 其中向量a, b用0填充使其長度達到2n。

[UR#34]多項式乘法

遞歸實現：參考hzwer的代碼

#include <bits/stdc++.h> #define pi acos(-1.0) typedef long long ll; using namespace std; int n,m; complex<double> a[300000],b[300000]; void fft(complex<double> *x, int n, int ty){if(n==1)return;complex<double> l[n>>1],r[n>>1];for(int i=0;i<=n;i+=2)l[i>>1]=x[i],r[i>>1]=x[i+1];fft(l,n>>1,ty);fft(r,n>>1,ty);complex<double> wn(cos(2*pi/n),sin(ty*2*pi/n)),w(1,0),t;for(int i=0;i<n>>1;i++,w*=wn)t=w*r[i],x[i]=l[i]+t,x[i+(n>>1)]=l[i]-t; } int main() {scanf("%d%d",&n,&m);for(int i=0,x;i<=n;++i)scanf("%d",&x),a[i]=x;for(int i=0,x;i<=m;++i)scanf("%d",&x),b[i]=x;m=n+m;for(n=1;n<=m;n<<=1);fft(a,n,1);fft(b,n,1);for(int i=0;i<=n;++i)a[i]=a[i]*b[i];fft(a,n,-1);for(int i=0;i<=m;++i)printf("%d ",(int)(a[i].real()/n+0.5));return 0; }

　　

　　然而，實際應用中，需要算法盡量地快速，因此，我們引出一種小常數FFT的迭代實現。

高效FFT實現

　　1.蝴蝶操作：

for(int i=0;i<n>>1;i++,w*=wn)t=w*r[i],x[i]=l[i]+t,x[i+(n>>1)]=l[i]-t;

　　2.考慮遞歸調用的過程，將每次輸入向量構成一顆樹，類似于zkw線段樹，如果可以實現自底向上更新，相比遞歸常數就很優秀了。

　　 ?(1)首先，我們成對的取出輸入向量樹最底層的元素，利用蝴蝶操作計算出每對的DFT，這樣我們就有了n/2個二元素DFT，同理，下一步取出這n/2個DFT，計算四元素DFT，最終形成n元素DFT

　　 ?(2)那么如何確定最底層的元素排列呢？嘗試一下就能發現，元素出現的順序就是是一個位逆序置換，即：n = 8時，{0, 4, 2, 6, 1, 5, 3, 7}, 寫成二進制形式，即{000, 100, 010, 110, 001, 101, 011, 111}，將所有元素二進制位逆序可得{000, 001, 010, 011, 100, 101, 110, 111} = {0, 1, 2 ,3, 4, 5, 6, 7}

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

?

[UR#34]多項式乘法

迭代實現：參考hzwer的代碼

#include <bits/stdc++.h> #define pi acos(-1.0) typedef complex<double> E; using namespace std; int n,m,L,R[300000]; E a[300000],b[300000]; void fft(E *a,int f){for(int i=0;i<n;++i)if(i<R[i])swap(a[i],a[R[i]]);for(int i=1;i<n;i<<=1){E wn(cos(pi/i),f*sin(pi/i));for(int p=(i<<1),j=0;j<n;j+=p){E w(1,0);for(int k=0;k<i;++k,w*=wn){E x=a[j+k],y=w*a[j+k+i];a[j+k]=x+y;a[j+k+i]=x-y;}}} } int main() {scanf("%d%d",&n,&m);for(int i=0,x;i<=n;++i)scanf("%d",&x),a[i]=x;for(int i=0,x;i<=m;++i)scanf("%d",&x),b[i]=x;m=n+m;for(n=1;n<=m;n<<=1)L++;for(int i=0;i<n;++i)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));fft(a,1);fft(b,1);for(int i=0;i<=n;++i)a[i]=a[i]*b[i];fft(a,-1);for(int i=0;i<=m;++i)printf("%d ",(int)(a[i].real()/n+0.5));return 0; }

?

NTT(Number Theory Transform)

　　可是，復數運算的常數依然是很巨大的，而且精度不高。當題目給的模數很特殊的時候，我們就可以利用數論變換來加速算法。

　　

?

轉載于:https://www.cnblogs.com/RRRR-wys/p/8868386.html

總結

以上是生活随笔為你收集整理的FFT算法学习笔记的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。