# bzoj2194 快速傅立叶之二

fft详解(这次是真的)

fft用于快速计算多项式乘法 满足形如c[n] = sigma(a[i] * b[n-i])的式子都可以有效计算

Code:

 1 #include<cstdio>
2 #include<cstring>
3 #include<cmath>
4 #include<algorithm>
5 #define rep(i,a,n) for(int i = a;i <= n;++i)
6 #define per(i,n,a) for(int i = n;i >= a;--i)
7 #define inf 2147483647
8 #define ms(a,b) memset(a,b,sizeof a)
9 using namespace std;
10 typedef double D;
11 typedef long long ll;
12 const D PI = acos(-1);
14     ll as = 0,fu = 1;
15     char c = getchar();
16     while(c < '0' || c > '9') {
17         if(c == '-') fu = -1;
18         c = getchar();
19     }
20     while(c >= '0' && c <= '9') {
21         as = as * 10 + c- '0';
22         c = getchar();
23     }
24     return as * fu;
25 }
27 const int N = 1<<21;
28 struct cp {
29     D x,y;
30     cp(D X = 0,D Y = 0) {x = X,y = Y;}
31     cp operator + (const cp &o) const {return cp(x+o.x,y+o.y);}
32     cp operator - (const cp &o) const {return cp(x-o.x,y-o.y);}
33     cp operator * (const cp &o) const {return cp(x*o.x-y*o.y,x*o.y+y*o.x);}
34 } a[N],b[N];
35
36 int r[N],lim = 1,base;
37 void fft(cp *a,int d) {
38     rep(i,0,lim-1) if(i < r[i]) swap(a[i],a[r[i]]);
39     for(int i = 1;i < lim;i <<= 1) {
40         cp T(cos(PI/i),d * sin(PI/i));
41         for(int k = 0;k < lim;k += (i<<1)) {
42             cp t(1,0);
43             rep(l,0,i-1) {
44                 cp nx = a[k+l],ny = a[k+l+i] * t;
45                 a[k+l] = nx+ny,a[k+l+i] = nx-ny;
46                 t = t * T;
47             }
48         }
49     }
50 }
51 int n;
52 int main() {
54     while(lim <= n<<1) lim <<= 1,base++;
55     rep(i,0,lim-1) r[i] = (r[i>>1]>>1) | ((i&1) << base-1);
56
65 }