2
4
2016
0

【bzoj2194】快速傅里叶之二

2194: 快速傅立叶之二

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 758  Solved: 433
[Submit][Status][Discuss]

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

题目给的式子和卷积的形式非常像,我们只要稍作改变就行了。

由于正常的卷积是sigma(a[i]*b[k-i]),而题目里是要求sigma(a[i]*b[i-k]),那么我们可以先把b数组翻转一下,变成求sigma(a[i]*b[n-1+k-i]),发现这等于c[n-1+k],c为fft后的数组,所以对于每个i,我们只要输出c[n-1+i]就可以了。

代码:

#include<cstdio>
#include<cmath>
#include<algorithm>
#define N 500010
using namespace std;
double pi=acos(-1);
int n,m;
struct cp{double x,y;};
cp a[N],b[N],cur[N];
cp operator *(cp x,cp y){return (cp){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x};}
cp operator +(cp x,cp y){return (cp){x.x+y.x,x.y+y.y};}
cp operator -(cp x,cp y){return (cp){x.x-y.x,x.y-y.y};}
void fft(cp *a,int n,int fl)
{
    for(int i=n>>1,j=1;j<n;j++)
    {
        if (i<j) swap(a[i],a[j]);
        int k=n>>1;
        for (;k&i;i^=k,k>>=1);i^=k;
    }
    for(int m=2;m<=n;m<<=1)
    {
        cp w=(cp){cos(2*pi*fl/m),sin(2*pi*fl/m)};
        cur[0]=(cp){1,0};
        for (int i=1;i<m;i++) cur[i]=cur[i-1]*w;
        for (int i=0;i<n;i+=m)
            for (int j=i;j<i+(m>>1);++j)
            {
                cp u=a[j],v=a[j+(m>>1)]*cur[j-i];
                a[j]=u+v;
                a[j+(m>>1)]=u-v;
            }
    }
}
int main()
{
    scanf("%d",&n);
    for (int i=0;i<n;i++) scanf("%lf%lf",&a[i].x,&b[n-i-1].x);
    m=1;
    while (m<=n*2) m*=2;
    fft(a,m,1);
    fft(b,m,1);
    for (int i=0;i<m;i++) a[i]=a[i]*b[i];
    fft(a,m,-1);
    for (int i=0;i<n;i++) printf("%d\n",(int)(a[n-1+i].x/m+0.5));
}
Category: FFT | Tags: | Read Count: 581

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com