首页 > 代码库 > FFT —— 快速傅里叶变换
FFT —— 快速傅里叶变换
问题:
已知A[], B[], 求C[],使:
定义C是A,B的卷积,例如多项式乘法等。
朴素做法是按照定义枚举i和j,但这样时间复杂度是O(n2).
能不能使时间复杂度降下来呢?
点值表示法:
我们把A,B,C看作表达式。
即:
A(x)=a0 + a1* x + a2 * x2 +...
将A={(x1,A(x1)), (x2,A(x2)), (x3,A(x3))...}叫做A的点值表示法。
那么使用点值表示法做多项式乘法就很简单了:对应项相乘。
那么,如何将A和B转换成点值表示法,再将C转化回系数表示法(即最初的表示方法)呢?
如果任取n个点,按照定义计算,那么还是O(n2)的。
这样就要用到快速傅里叶变换。
快速傅里叶变换:
既然任取n个点,按照定义计算太慢,就要找一些特殊点。
我们用n个n次单位复数根(1的n次方根,涉及到复数,1的方根不止1和-1)来计算:
1的n次方根是,其中i是虚数单位。
我们定义wn = e^(2∏i)是主n次单位根,那么所有n次单位复数根都是它的次方。
我们要求出A(wnk),就要采用分治思想。
我们将奇偶系数分离(先假设n为偶数),即定义
A1(x)=a0 + a2* x + a4 * x2 +...
A2(x)=a1 + a3* x + a5 * x2 +...
那么A(x)=A1(x2) + xA2(x2)。
要计算A(wnk)=A1((wnk)2) + wnkA2((wnk)2),
就要用到(wnk)2 = wn/2k mod (n / 2)(证略)。
所以A(wnk)=A1(wn/2k mod (n / 2)) + wnkA2(wn/2k mod (n / 2))
我们发现A1,A2都是n/2项的,且只需要算wn/2k的值,那么这就和开始的问题一样了,可以分治。
边界也很容易:n=1的时候A1本身就是值。
合并解。
A(wnk)=A1(wn/2k mod (n / 2)) + wnkA2(wn/2k mod (n / 2))
那么可以A(wnk)和A(wnk+n/2)一起算(0<=k<n/2):
设u = A1(wn/2k), t = wnkA2(wn/2k),
那么A(wnk) = u + t
A(wnk+n/2) = A1(wn/2k) + wnk + n/2 A2(wn/2k)
= A1(wn/2k) + wnk wnn/2 A2(wn/2k)
= A1(wn/2k) - wnk A2(wn/2k)
= u - t
所以这样就能算出A的点值表示法。
一个问题:分治要求n是2的幂,不是怎么办? 补0, 直到是2的幂。
剩下的问题:如何把C转化回系数表示法。
快速傅里叶逆变换:
我们把C做一遍快速傅立叶变换,只是求的是wnn, wnn-1, ..., wn1的值而不是wn0, wn1, ..., wnn-1的值,最后每一项除以n即可。
证略。
1 void Rader(complex y[],int len) 2 { 3 int i,j,k; 4 for(i = 1, j = len/2;i < len-1;i++) 5 { 6 if(i < j)swap(y[i],y[j]); 7 k = len/2; 8 while( j >= k) 9 { 10 j -= k; 11 k /= 2; 12 } 13 if(j < k)j += k; 14 } 15 } 16 void FFT(complex y[],int len,int on) //on = 1 快速傅里叶变换, on = 0 快速傅里叶逆变换 17 { 18 Rader(y,len); 19 for(int h = 2;h <= len;h <<= 1) 20 { 21 complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h)); //e^ki = cosk + isink 22 for(int j = 0;j < len;j += h) 23 { 24 complex w(1,0); 25 for(int k = j;k < j+h/2;k++) 26 { 27 complex u = y[k]; 28 complex t = w*y[k+h/2]; 29 y[k] = u+t; 30 y[k+h/2] = u-t; 31 w = w*wn; 32 } 33 } 34 } 35 if(on == -1) 36 for(int i = 0;i < len;i++) 37 y[i].r /= len; 38 }
39 //复数实现略
FFT —— 快速傅里叶变换