[BZOJ2194]快速傅立叶之二

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

 

Solution

  FFT表示从上学期开始到昨晚才完全搞明白。。。

  分两步:蝶形变换,DFT。具体过程见http://picks.logdown.com/posts/177631-fast-fourier-transform

  自己手推一遍终于弄懂了= =

 

 1 /**************************************************************
 2     Problem: 2194
 3     User: wjy1998
 4     Language: C++
 5     Result: Accepted
 6     Time:2660 ms
 7     Memory:19788 kb
 8 ****************************************************************/
 9  
10 #include<cstdio>
11 #include<cmath>
12 const double PI=3.14159265359;
13 struct P{double x,y;};
14 P operator+(const P&a,const P&b){return (P){a.x+b.x,a.y+b.y};}
15 P operator-(const P&a,const P&b){return (P){a.x-b.x,a.y-b.y};}
16 P operator*(const P&a,double p){return (P){a.x*p,a.y*p};}
17 P operator*(const P&a,const P&b){double d=a.x*b.x,e=a.y*b.y,f=(a.x+a.y)*(b.x+b.y);return (P){d-e,f-e-d};}
18 P operator/(const P&a,double p){return (P){a.x/p,a.y/p};}
19 P operator/(const P&a,const P&b){double x=b.x*b.x+b.y*b.y;return (P){(a.x*b.x+a.y*b.y)/x,(a.y*b.x-a.x*b.y)/x};}
20  
21 int a[270000],b[270000],n,m;
22 P w[2][270000],x[270000],y[270000];
23 void FFT(P*x,int k,int v)
24 {
25     int i,j,l;P tmp;
26     for(i=j=0;i<k;i++)
27     {
28         if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
29         for(l=k>>1;(j^=l)<l;l>>=1);
30     }
31     for(i=2;i<=k;i<<=1)
32     for(j=0;j<k;j+=i)
33     for(l=0;l<i>>1;l++)
34     {
35         tmp=x[j+l+(i>>1)]*w[v][k/i*l];
36         x[j+l+(i>>1)]=x[j+l]-tmp;
37         x[j+l]=x[j+l]+tmp;
38     }
39 }
40 int main(){
41     scanf("%d",&m);int i;
42     for(i=0;i<m;i++)scanf("%d%d",&a[i],&b[m-i]);
43     for(n=1;n<m<<1;n<<=1);
44     for(i=0;i<=n;i++)w[0][i]=(P){cos(PI*2*i/n),sin(PI*2*i/n)};
45     for(i=0;i<=n;i++)w[1][i]=w[0][n-i];
46     for(i=0;i<n;i++)x[i]=(P){a[i],0};FFT(x,n,0);
47     for(i=0;i<n;i++)y[i]=(P){b[i],0};FFT(y,n,0);
48     for(i=0;i<n;i++)x[i]=x[i]*y[i];FFT(x,n,1);
49     for(i=0;i<m;i++)printf("%d\n",(int)(x[m+i].x/n+0.5));
50 }
View Code
posted @ 2014-07-19 08:28  n+e  阅读(422)  评论(0编辑  收藏  举报