bzoj2194 快速傅立叶之二

传送门

上面骗人的

请计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n ，并且有 n < = 10 ^ 5。 a,b中的元素均为小于等于100的非负整数。

fft详解(这次是真的)

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

这题是上面形式反过来

我们可以再返回去

把a,c数组翻转完发现上面的式子就是c[n-k] = sigma(a[n-i] * b[i-k])

所以fft板子粘上就行

注意lim开到n<<1

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);
13 ll read() {
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 }
26 //head
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() {
53     n = read();
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 
57     rep(i,0,n-1) a[i] = read(),b[i] = read();
58     reverse(a,a+n);
59 
60     fft(a,1),fft(b,1);
61     rep(i,0,lim-1) a[i] = a[i] * b[i];
62     fft(a,-1);
63     per(i,n-1,0) printf("%d
", int(a[i].x / lim + 0.5));
64     return 0;
65 }