FFT 模板

Posted 2020-07-07 Pat

FFT（Fast Fourier Transformation/快速傅立叶变换），确切地说应该称之为FDFT（Fast Discrete Fourier Transformation/快速离散傅立叶变换），因为FFT是为解DFT问题而设计的一种快速算法。在深入讨论之前，有必要特别指出这一点。

DFT问题：

给定一个复数域上的n-1次多项式A(x)的系数表示（cofficient representation）(a[0], a[1], ..., a[n-1])，求A(x)的某个点-值表示（point-value representation）。

递归FFT的实现（《算法导论》）

 1 #include <bits/stdc++.h>
 2 #define rep(i, l, r) for(int i=l; i<r; i++)
 3 using namespace std;
 4 const double PI(acos(-1));
 5 
 6 struct P{
 7     double r, i;
 8     P(double r, double i):r(r), i(i){}
 9     P(int n):r(cos(2*PI/n)), i(sin(2*PI/n)){}    //!!error-prone
10     P():r(0), i(0){}    //default constructor
11     void operator *=(const P &a){
12         double R=r*a.r-i*a.i, I=r*a.i+a.r*i;
13         r=R, i=I;
14     }
15     P operator+(const P a){
16         return P(r+a.r, i+a.i);
17     }
18     P operator-(const P a){
19         return P(r-a.r, i-a.i);
20     }
21     P operator*(const P a){
22         return P(r*a.r-i*a.i, r*a.i+a.r*i);
23     }
24     void out(){
25         // cout<<r<<‘ ‘<<i<<endl;
26     }
27 };
28 
29 
30 typedef vector<P> vec;
31 
32 vec FFT(vec a, int t){
33     int n=a.size();    //n is a power of 2
34     if(n==1) return a;
35     P wn(t*n), w(1, 0);
36     vec b[2], y[2];
37     for(int i=0; i<n; i++)
38         b[i%2].push_back(a[i]);
39 
40     y[0]=FFT(b[0], t), y[1]=FFT(b[1], t);
41     vec res(n);
42     for(int i=0, h=n>>1; i<h; i++){
43         res[i]=y[0][i]+w*y[1][i];    
44         res[i+h]=y[0][i]-w*y[1][i];
45         w*=wn;
46     }
47     return res;
48 }
49 
50 const int N(5e4+5);
51 char s[N], t[N];
52 
53 int trans(int x){
54     int i=0;
55     for(; x>1<<i; i++);
56     return 1<<i;
57 }
58 
59 vec a, b, Y1, y2, res;
60 
61 int ans[N];
62 
63 int main(){
64     for(; ~scanf("%s%s", s, t); ){
65         int n=strlen(s), m=strlen(t), l=trans(n+m-1);
66         a.clear(), b.clear();    //error-prone
67         for(int i=0; i<n; i++) a.push_back(P(s[n-1-i]-‘0‘, 0)); a.resize(l);    //error-prone
68         for(int i=0; i<m; i++) b.push_back(P(t[m-1-i]-‘0‘, 0)); b.resize(l);
69         Y1=FFT(a, 1), y2=FFT(b, 1);
70         rep(i, 0, l) Y1[i]*=y2[i];
71         res=FFT(Y1, -1);
72         rep(i, 0, l) res[i].r/=l;    //error-prone
73         rep(i, 0, l) ans[i]=(int)(res[i].r+0.5); ans[l]=0;    //error-prone
74         rep(i, 0, l) ans[i+1]+=ans[i]/10, ans[i]%=10;
75         int p=l;
76         for(;p && !ans[p]; --p);
77         for(; ~p; putchar(ans[p--]+‘0‘));    //error-prone
78         puts("");
79     }
80     return 0;
81 }

Comment:

1. 此代码是为HDU1402写的。代码中，凡注释error-prone处，都应特别小心。我犯的最傻逼的错误是第9行，应当是2*PI，我写成PI了。

2. 这个FFT的递归实现完全是参照《算法导论》(3rd. Ed. Chapter 30, Polynomials and the FFT) 的，但这个实现常数大，空间复杂度高，TLE了

3. 迭代实现后面会补上

4. FFT的数值稳定性（精度）问题，还有待考虑。