FFT —— 快速傅里叶变换

Posted 2020-08-22

问题：

　　已知A[], B[], 求C[],使：

　　　　

　　定义C是A,B的卷积，例如多项式乘法等。

 

　　朴素做法是按照定义枚举i和j，但这样时间复杂度是O(n²).

　　能不能使时间复杂度降下来呢？

 

点值表示法：

　　我们把A，B，C看作表达式。

　　即：

　　　　A(x)=a₀ + a₁* x + a₂ * x² +...

　　将A={(x₁,A(x₁)), (x₂,A(x₂)), (x₃,A(x₃))...}叫做A的点值表示法。

　　那么使用点值表示法做多项式乘法就很简单了：对应项相乘。

　　那么，如何将A和B转换成点值表示法，再将C转化回系数表示法（即最初的表示方法）呢？

　　如果任取n个点，按照定义计算，那么还是O(n²)的。

　　这样就要用到快速傅里叶变换。

 

快速傅里叶变换：

　　既然任取n个点，按照定义计算太慢，就要找一些特殊点。

　　我们用n个n次单位复数根（1的n次方根，涉及到复数，1的方根不止1和-1）来计算：

　　　　1的n次方根是,其中i是虚数单位。

　　我们定义w_n = e^(2∏i)是主n次单位根，那么所有n次单位复数根都是它的次方。

　　我们要求出A(w_n^k)，就要采用分治思想。

　　我们将奇偶系数分离（先假设n为偶数），即定义

　　　　　　A1(x)=a₀ + a₂* x + a₄ * x² +...

　　　　　　A2(x)=a₁ + a₃* x + a₅ * x² +...

　　那么A(x)=A1(x²) + xA2(x²)。

　　要计算A(w_n^k)=A1((w_n^k)²) + w_n^kA2((w_n^k)²),

　　就要用到(w_n^k)² = w_n/2^{k mod (n / 2)}(证略)。

　　所以A(w_n^k)=A1(w_n/2^{k mod (n / 2)}) + w_n^kA2(w_n/2^{k mod (n / 2)})

　　我们发现A1,A2都是n/2项的，且只需要算w_n/2^k的值，那么这就和开始的问题一样了，可以分治。

　　边界也很容易：n=1的时候A₁本身就是值。

　　合并解。

　　A(w_n^k)=A1(w_n/2^{k mod (n / 2)}) + w_n^kA2(w_n/2^{k mod (n / 2)})

　　那么可以A(w_n^k)和A(w_n^k+n/2)一起算(0<=k<n/2)：

　　　　设u = A1(w_n/2^k), t = w_n^kA2(w_n/2^k),

　　　　那么A(w_n^k) = u + t

　　　　　　A(w_n^k+n/2) = A1(w_n/2^k) + w_n^{k + n/2} A2(w_n/2^k)

　　　　　　　　　　　　= A1(w_n/2^k) + w_n^k  w_n^n/2 A2(w_n/2^k)

　　　　　　　　　　　　= A1(w_n/2^k) - w_n^k A2(w_n/2^k)

　　　　　　　　　　　　= u - t

　　所以这样就能算出A的点值表示法。

　　一个问题：分治要求n是2的幂，不是怎么办？ 补0， 直到是2的幂。

 

　　剩下的问题：如何把C转化回系数表示法。

快速傅里叶逆变换：

　　我们把C做一遍快速傅立叶变换，只是求的是w_nⁿ, w_n^n-1, ..., w_n¹的值而不是w_n⁰, w_n¹, ..., w_n^n-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 //复数实现略