【BZOJ 2194】2194: 快速傅立叶之二（FFT）

2194: 快速傅立叶之二
Time Limit: 10 Sec Memory Limit: 259 MB
Submit: 1273 Solved: 745
Description

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

Input

第一行一个整数N,接下来N行，第i+2..i+N-1行，每行两个数，依次表示a[i],b[i] (0 < = i < N)。

Output

输出N行，每行一个整数，第i行输出C[i-1]。

Sample Input

5
3 1
2 4
1 1
2 4
1 4

Sample Output

24
12
10
6
1

HINT

Source

【分析】

　　这个卷积也很容易看出来

　　C[k]=sigma(a[i]*b[i-k])

　　变形，先讲b数组倒置：

　　C[k]=sigma(a[i]*b[n-i+k])

　　令D[k]=sigma(a[i]*b[k-i])

　　则C[k]=D[n+k]

　　D数组就是标准卷积了，用FFT加速。

递归：

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cstring>
 4 #include<iostream>
 5 #include<algorithm>
 6 #include<cmath>
 7 using namespace std;
 8 #define Maxn 100010*4
 9 const double pi=3.141592653;
10 
11 struct P
12 {
13     double x,y;
14     P() {x=y=0;}
15     P(double x,double y):x(x),y(y){}
16     friend P operator + (P x,P y) {return P(x.x+y.x,x.y+y.y);}
17     friend P operator - (P x,P y) {return P(x.x-y.x,x.y-y.y);}
18     friend P operator * (P x,P y) {return P(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
19 }a[Maxn],b[Maxn];
20 
21 void fft(P *s,int n,int f)
22 {
23     if(n==1) return;
24     P a0[n>>1],a1[n>>1];
25     for(int i=0;i<=n;i+=2) a0[i>>1]=s[i],a1[i>>1]=s[i+1];
26     fft(a0,n>>1,f);fft(a1,n>>1,f);
27     P wn(cos(2*pi/n),f*sin(2*pi/n)),w(1,0);
28     for(int i=0;i<(n>>1);i++,w=w*wn) s[i]=a0[i]+w*a1[i],s[i+(n>>1)]=a0[i]-w*a1[i];
29 }
30 
31 int main()
32 {
33     int n;
34     scanf("%d",&n);n--;
35     for(int i=0;i<=n;i++) scanf("%lf%lf",&a[i].x,&b[i].x);
36     for(int i=0;i<=n/2;i++) swap(b[i].x,b[n-i].x);
37     int nn=1;
38     while(nn<=2*n) nn<<=1;
39     fft(a,nn,1);fft(b,nn,1);
40     for(int i=0;i<=nn;i++) a[i]=a[i]*b[i];
41     fft(a,nn,-1);
42     for(int i=n;i<=n+n;i++) printf("%d
",(int)(a[i].x/nn+0.5));
43     return 0;
44 }

View Code

2017-04-13 10:39:50