码迷,mamicode.com
首页 > 其他好文 > 详细

FFT 模板

时间:2016-05-17 22:33:04      阅读:238      评论:0      收藏:0      [点我收藏+]

标签:

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的数值稳定性(精度)问题,还有待考虑。

FFT 模板

标签:

原文地址:http://www.cnblogs.com/Patt/p/5503322.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!