bzoj2194 快速傅立叶之二 FFT

2194: 快速傅立叶之二

Time Limit: 10 Sec Memory Limit: 259 MB
Submit: 1697 Solved: 1004
[Submit][Status][Discuss]

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

卷积的形式

这样的形式化出来才可以fft

我们把b数组翻转，i-->-i
c[i+j]+=a[i]*b[j]

 1 #pragma GCC optimize(2)
 2 #pragma G++ optimize(2)
 3 #include<cstdio>
 4 #include<cmath>
 5 #include<iostream>
 6 #include<algorithm>
 7 #include<cstring>
 8 
 9 #define N 300007
10 #define pi acos(-1)
11 using namespace std;
12 inline int read()
13 {
14     int x=0,f=1;char ch=getchar();
15     while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
16     while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
17     return x*f;
18 }
19 
20 int n,m,L;
21 int rev[N];
22 struct comp
23 {
24     double r,v;
25     inline comp operator+(const comp &a){return (comp){r+a.r,v+a.v};}
26     inline comp operator-(const comp &a){return (comp){r-a.r,v-a.v};}
27     inline comp operator*(const comp &a){return (comp){r*a.r-v*a.v,r*a.v+v*a.r};}
28 }a[N],b[N];
29 
30 void FFT(comp *a,int f)
31 {
32     for (int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
33     for (int i=1;i<n;i<<=1)
34     {
35         comp wn=(comp){cos(pi/i),f*sin(pi/i)};
36         for (int j=0;j<n;j+=(i<<1))
37         {
38             comp w=(comp){1,0};
39             for (int k=0;k<i;k++,w=w*wn)
40             {
41                 comp x=a[j+k],y=w*a[j+k+i];
42                 a[j+k]=x+y,a[j+k+i]=x-y;
43             }
44         }
45     }
46     if(f==-1)for (int i=0;i<n;i++)a[i].r/=n;
47 }
48 int main()
49 {
50     n=read()-1,m=n*2;
51     for(int i=0;i<=n;i++)
52         a[i].r=(double)read(),b[n-i].r=(double)read();
53     for(n=1;n<=m;n<<=1,L++);
54     if(L)L--;
55     for(int i=0;i<n;i++)
56         rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
57     FFT(a,1),FFT(b,1);
58     for (int i=0;i<n;i++)a[i]=a[i]*b[i];
59     FFT(a,-1);
60     for (int i=m/2;i<=m;i++)
61         printf("%d
",(int)(a[i].r+0.5));
62 }