FFT算法学习笔记
寫在前邊
1.辣雞RRRR_wys之前csdn的博客，千年不更。。。還很水。。。于是開了這個Blog。。。妄圖拯救一下自己
2.最近接觸了一些多項式理論。于是翹掉了愉快的高頻自控，通過《算導(dǎo)》稍稍學(xué)習(xí)了一下
3.算法競賽中，FFT主要解決多項式的乘法等問題
FFT基本概念
1.FFT即快速傅里葉變換, 是離散傅里葉變換的加速算法。可以在o(nlogn)的時間內(nèi)，完成DFT和DFT-1
2.DFT即離散傅里葉變換, 這里主要就是多項式的系數(shù)向量轉(zhuǎn)換成點值表示的過程
多項式的表示
1.系數(shù)表達：對于一個次數(shù)上界為n的多項式?A(x)=∑(j=0 - n-1) aj*xj?而言，系數(shù)表達式便是由系數(shù)組成的向量a = (a0, a1, a2, a3, ... , an-1)
? ? ? ? ?(1)計算A(x0)的值即：A(x0) = a0+x0(a1+x0(a2+x0(a3+x0(a4+...+x0(an-2+x0*an-1))...))) ?時間為o(n)
?(2)通過系數(shù)向量計算兩個多項式的乘法：設(shè)輸出向量為c, 而a,b長度n=max(na,nb)那么顯然兩個c的長度為2n-1，計算方法為：cj =?∑(k=0 ~ j)?ak*bj-k，
? ? ?即向量a,b的卷積
2.點值表達式：對于一個次數(shù)上界為n的多項式?A(x)=∑(j=0 ~ n-1)?aj*xj?而言，點值表達式便是由n個點值所構(gòu)成的集合{ (x0,y0),?(x1,y1),?(x2,y2), ... , (xn-1,yn-1)?}
? ?(1)運用之前提到的方法，一個個計算點值時間為o(n2)，之后可以看到如果我們巧妙地選取點xk可以使時間變?yōu)閛(nlogn)
? ?(2)將點值表達式轉(zhuǎn)換為系數(shù)表達式：將已知的點值帶入，多項式可已得到一個有n個未知數(shù)的n個線性方程
?a0 + a1*x0 + a2*x02 + a3* x03 + ... + an-1*x0n-1 = y0?(1)
a0 + a1*x1 + a2*x12 + a3* x13 + ... + an-1*x1n-1 = y1?(2)
.......
a0?+ a1*xn-1?+ a2*xn-12?+ a3* xn-13?+ ... + an-1*xn-1n-1?= yn-1?(n)
? 可以寫成矩陣形式，既為一個范德蒙矩陣，通過線代知識解這個方程就可以求出系數(shù)向量，復(fù)雜度o(n3)。另一種方法是采用拉格朗日插值法求解，復(fù)雜度o(n2)。因此，n個點的求值運算和插值運算，是定義完備的互逆運算。
(3)通過點值計算兩個多項式的乘積，c =?{ (x0,yA0*yB0),?(x1,yA1*yB1),?(x2,yA2*yB2), ... , (xn-1,yA(n-1)*yB(n-1))?}, 可以看到復(fù)雜度為o(n)
??(4)計算多項式乘法的過程：1）o(nlogn)分別求出兩個多項式的點值表達式 2）點值乘法o(n) 3）插值o(nlogn) 4）求出答案多項式的系數(shù)表達式
DFT與FFT
下面開始討論，如何在o(nlogn)的時間內(nèi)，完成求值和插值運算
1.單位復(fù)數(shù)根
??? n次單位復(fù)數(shù)根滿足ωn=1的復(fù)數(shù)ω，n此單位復(fù)數(shù)根的數(shù)目恰好有n個，這些根是：ωnk?= e2πik/n 這n個根均勻的分布在復(fù)平面上，其他n次單位根都是ωn的冪次。n個n次單位根在乘法意義下形成一個群，該群與模n意義下的加法群，有相同的結(jié)構(gòu)，即ωni ωnj = ωn(i + j) mod n
(1)?消去引理：ωndkd =ωnk ?(n>=0, k>=0, d>0)
(2)?ωnn/2 =ω2 =-1
(3)?折半引理：如果n>0為偶數(shù)，那么n個n次單位復(fù)數(shù)根的平方的集合就是n/2個n/2次單位復(fù)數(shù)根的集合
(4)?求和引理：∑(j=0 ~ n-1) (ωnk)j = 0 (n>=1, 不能被n整除的非負整數(shù)k)
2.DFT
我們希望計算ωn0, ωn1, ωn2, .... , ωnn-1處的值（在多項式乘法中其實是2n個點值）。假設(shè)給定系數(shù)a=(a0, a1, a2, ... , an-1), 即對于k=0, 1, 2, 3, ... , n-1求出yk?= A(ωnk)，記向量y=(y1, y2, y3, ... ,yn-1)就是系數(shù)向量a的離散傅里葉變換DFT，記為 y?= DFTn(a)
3. FFT
通過快速傅里葉變換FFT，利用單位復(fù)數(shù)根的性質(zhì)，我們可以在o(nlogn)的時間內(nèi)計算DFTn(a)。首先假設(shè)次數(shù)n恰好是2的整數(shù)冪，不足在高次位置項添0。
定義兩個新的次數(shù)界為n/2的多項式?A[0]?= a0?+ a2*x?+ a4*x2?+ a6* x3?+ ... + an-2*xn/2-1 , A[1]?= a1?+ a3*x?+ a5*x2?+ a7* x3?+ ... + an-1*xn/2-1
A[0]包含所有偶數(shù)下標(biāo)的系數(shù)，A[1]包含所有奇數(shù)下標(biāo)的系數(shù)，于是有A(x) = A[0](x2) + x*A[1](x2)?。
所以，求解A(x)的DFN，轉(zhuǎn)換為了求解多項式A[0], A[1]在(ωn0)2, (ωn1)2, (ωn2)2?, ... , (ωnn-1)2?處的取值。
根據(jù)折半引理，(ωn0)2, (ωn1)2, (ωn2)2?, ... , (ωnn-1)2?并不是由n個不同的值組成，而是由n/2個n/2次單位根復(fù)數(shù)根所組成，每個根正好出現(xiàn)兩次。因此，我們遞歸的對次數(shù)界為n/2的多項式A[0], A[1]在n/2個n/2次單位復(fù)數(shù)根處進行求值。這樣我們求出了次數(shù)界為n的次多項式在n次單位復(fù)數(shù)根處的值。
注意到除了遞歸調(diào)用的時間外，每次所需時間為o(n)，其中n是輸入向量的長度。因此有如下遞歸式：T(n) = 2T(n/2) + o(n) = o(nlogn)，時間為o(nlogn)。
計算在單位復(fù)數(shù)根處插值，求解系數(shù)向量。根據(jù)點值表達式，我們可以列出線性方程組：
a0?+ a1?+ a2?+ a3+ ... + an-1 = y0 (1)
a0?+ a1* (ωn)1+ a2*(ωn)2?+ a3* (ωn)3?+ ... + an-1*(ωn)n-1 = y1?(2)
a0?+ a1*?(ωn)2+ a2*(ωn)4?+ a3*?(ωn)6?+ ... + an-1*(ωn)2(n-1)? = ?y2 ?(3)
......
a0?+ a1*?(ωn)n-1+ a2*(ωn)2(n-1)?+ a3*?(ωn)3(n-1)?+ ... + an-1*(ωn)(n-1)(n-1)? = ?yn-1?(4)
根據(jù)這個線性方程組，我們可以令單位復(fù)數(shù)根組成的系數(shù)矩陣為Vn，運用yVn-1求解系數(shù)。可以證明：Vn-1的(j, k)處的值為(ωn)-kj/n ?(j,k = 0, 1, 2, ... , n-1)
有了上面的式子我們能推導(dǎo)出DFT-1：aj =(1/n)*?∑(k=0 ~ n-1) yk(ωn)-kj ?(j = 0, 1, 2, ... , n-1)
觀察可得，相對于DFT，我們把a與y替換，用ωn-1代替ωn，并將每個元素除以n。這樣，我們也可以在o(nlogn)的時間內(nèi)計算出DFT-1
卷積定理：a(卷積)b = DFT2n-1( DFT2n(a)*DFT2n(b) ) , 其中向量a, b用0填充使其長度達到2n。
[UR#34]多項式乘法
遞歸實現(xiàn)：參考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; }
然而，實際應(yīng)用中，需要算法盡量地快速，因此，我們引出一種小常數(shù)FFT的迭代實現(xiàn)。
高效FFT實現(xiàn)
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.考慮遞歸調(diào)用的過程，將每次輸入向量構(gòu)成一顆樹，類似于zkw線段樹，如果可以實現(xiàn)自底向上更新，相比遞歸常數(shù)就很優(yōu)秀了。
?(1)首先，我們成對的取出輸入向量樹最底層的元素，利用蝴蝶操作計算出每對的DFT，這樣我們就有了n/2個二元素DFT，同理，下一步取出這n/2個DFT，計算四元素DFT，最終形成n元素DFT
?(2)那么如何確定最底層的元素排列呢？嘗試一下就能發(fā)現(xiàn)，元素出現(xiàn)的順序就是是一個位逆序置換，即：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]多項式乘法
迭代實現(xiàn)：參考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)
可是，復(fù)數(shù)運算的常數(shù)依然是很巨大的，而且精度不高。當(dāng)題目給的模數(shù)很特殊的時候，我們就可以利用數(shù)論變換來加速算法。
?
轉(zhuǎn)載于:https://www.cnblogs.com/RRRR-wys/p/8868386.html
總結(jié)
- 上一篇： 学徒表哥的电脑坏了寄过来维修学徒表哥的电
- 下一篇： Full_of_Boys训练1总结