首页 > 代码库 > [BZOJ2194]快速傅立叶之二

[BZOJ2194]快速傅立叶之二

2024-07-12 20:55:03   222人阅读

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


联系
我们