【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加速。

 

递归：

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
#define Maxn 100010*4
const double pi=3.141592653;

struct P
{
    double x,y;
    P() {x=y=0;}
    P(double x,double y):x(x),y(y){}
    friend P operator + (P x,P y) {return P(x.x+y.x,x.y+y.y);}
    friend P operator - (P x,P y) {return P(x.x-y.x,x.y-y.y);}
    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);}
}a[Maxn],b[Maxn];

void fft(P *s,int n,int f)
{
    if(n==1) return;
    P a0[n>>1],a1[n>>1];
    for(int i=0;i<=n;i+=2) a0[i>>1]=s[i],a1[i>>1]=s[i+1];
    fft(a0,n>>1,f);fft(a1,n>>1,f);
    P wn(cos(2*pi/n),f*sin(2*pi/n)),w(1,0);
    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];
}

int main()
{
    int n;
    scanf("%d",&n);n--;
    for(int i=0;i<=n;i++) scanf("%lf%lf",&a[i].x,&b[i].x);
    for(int i=0;i<=n/2;i++) swap(b[i].x,b[n-i].x);
    int nn=1;
    while(nn<=2*n) nn<<=1;
    fft(a,nn,1);fft(b,nn,1);
    for(int i=0;i<=nn;i++) a[i]=a[i]*b[i];
    fft(a,nn,-1);
    for(int i=n;i<=n+n;i++) printf("%d\n",(int)(a[i].x/nn+0.5));
    return 0;
}

 

2017-04-13 10:39:50