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

HINT

Source

[ Submit][ Status][ Discuss] 


分析：
其实我们之前遇到过类似式子： bzoj3527
题目给出的式子实在太丑了，美化一下：

$C_k={\sum_{i=k}^{n-1}a_ib_{i-k}}$

化式子吧

$c_k={\sum_{i=k}^{n-1}a_ib_{i-k}}$

$={\sum_{i-k=0}^{n-k-1}a_ib_{i-k}}$

$={\sum_{j=1}^{n-k-1}a_{j+k}b_{j}}$ $(j=i-k)$

设$N=n-k-1,c'_N=c_k$

$c'_N={\sum_{j=1}^{N}a_{n-1-(N-j)}b_{j}}$

设$a'_{N-j}=a_{n-1-(N-j)},a'_i=a_{n-1-i}$

$c'_N={\sum_{j=0}^{N}a'_{N-j}b_{j}}$

标准卷积，上FFT

tip

好久不写FFT
注意下标从0开始
传入的长度是fn（2的整数幂）

#include<bits/stdc++.h>

using namespace std;

const double pi=acos(-1.0);
const int N=400005;
struct node{
    double x,y;
    node (double xx=0,double yy=0) {
        x=xx; y=yy;
    }
};
node a[N],b[N],omega[N],a_omega[N];
int n,m,fn;

node operator +(const node &a,const node &b) {return node(a.x+b.x,a.y+b.y);}
node operator -(const node &a,const node &b) {return node (a.x-b.x,a.y-b.y);}
node operator *(const node &a,const node &b) {return node(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

void init(int n) {
    for (int i=0;i<n;i++) {
        omega[i]=node(cos(2.0*i*pi/n),sin(2.0*i*pi/n));
        a_omega[i]=node(cos(2.0*i*pi/n),-sin(2.0*i*pi/n));
    }
}

void FFT(int n,node *a,node *w) {
    int i,j=0,k;
    for (i=0;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (i=2;i<=n;i<<=1) {
        int m=i>>1;
        for (j=0;j<n;j+=i) 
            for (k=0;k<m;k++) {
                node z=a[j+k+m]*w[n/i*k];
                a[j+k+m]=a[j+k]-z;
                a[j+k]=a[j+k]+z;
            }
    }
}

double p[N],q[N],ans[N];

int main()
{
    scanf("%d",&n);
    init(fn);
    for (int i=0;i<n;i++) scanf("%lf%lf",&p[i],&q[i]);

    for (int i=0;i<n;i++) {  //从0开始 
        a[i].x=p[n-i-1];     //i--->n-1-i
        b[i].x=q[i];
    }

    fn=1;
    while (fn<=2*n) fn<<=1;
    init(fn);
    FFT(fn,a,omega);
    FFT(fn,b,omega);
    for (int i=0;i<=fn;i++) a[i]=a[i]*b[i];
    FFT(fn,a,a_omega);

    for (int i=0;i<n;i++) 
        ans[n-i-1]=a[i].x/fn;
    for (int i=0;i<n;i++) 
        printf("%d\n",(int)(ans[i]+0.5));

    return 0;
}