[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

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

 

/**************************************************************
    Problem: 2194
    User: wjy1998
    Language: C++
    Result: Accepted
    Time:2660 ms
    Memory:19788 kb
****************************************************************/
 
#include<cstdio>
#include<cmath>
const double PI=3.14159265359;
struct P{double x,y;};
P operator+(const P&a,const P&b){return (P){a.x+b.x,a.y+b.y};}
P operator-(const P&a,const P&b){return (P){a.x-b.x,a.y-b.y};}
P operator*(const P&a,double p){return (P){a.x*p,a.y*p};}
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};}
P operator/(const P&a,double p){return (P){a.x/p,a.y/p};}
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};}
 
int a[270000],b[270000],n,m;
P w[2][270000],x[270000],y[270000];
void FFT(P*x,int k,int v)
{
    int i,j,l;P tmp;
    for(i=j=0;i<k;i++)
    {
        if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
        for(l=k>>1;(j^=l)<l;l>>=1);
    }
    for(i=2;i<=k;i<<=1)
    for(j=0;j<k;j+=i)
    for(l=0;l<i>>1;l++)
    {
        tmp=x[j+l+(i>>1)]*w[v][k/i*l];
        x[j+l+(i>>1)]=x[j+l]-tmp;
        x[j+l]=x[j+l]+tmp;
    }
}
int main(){
    scanf("%d",&m);int i;
    for(i=0;i<m;i++)scanf("%d%d",&a[i],&b[m-i]);
    for(n=1;n<m<<1;n<<=1);
    for(i=0;i<=n;i++)w[0][i]=(P){cos(PI*2*i/n),sin(PI*2*i/n)};
    for(i=0;i<=n;i++)w[1][i]=w[0][n-i];
    for(i=0;i<n;i++)x[i]=(P){a[i],0};FFT(x,n,0);
    for(i=0;i<n;i++)y[i]=(P){b[i],0};FFT(y,n,0);
    for(i=0;i<n;i++)x[i]=x[i]*y[i];FFT(x,n,1);
    for(i=0;i<m;i++)printf("%d\n",(int)(x[m+i].x/n+0.5));
}