首页 > 其他 > 详细

BZOJ 2194 [快速傅里叶变换 卷积]

时间:2017-02-11 12:16:17      阅读:318      评论:0      收藏:0      [点我收藏+]

2194: 快速傅立叶之二

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 1168  Solved: 673
[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]。 


 

卷积 (f x g)(n)=∑{f(i)*g(n-i):i=0...n-1}

多项式乘法就是一个系数向量的卷积

可以用FFT快速计算卷积

遇到和不是定值的情况可以反转一个向量

如本题反转a向量后

  c[k]=∑(a[n-i-1]*b[i-k]) k<=i<=n-1

更换求和指标 i=i-k

  c[k]=∑(a[n-i-k-1]*b[i]) 0<=i<=n-k-1

把-k-1消去,令t=n-k-1

  c[n-t-1]=∑(a[t-i]*b[i]) 0<=i<=t

这样就是标准的卷积形式啦

 

#include <iostream>
#include <cstdio>
#include <string>
#include <algorithm>
#include <cmath>
using namespace std;
const int N=3e5+5;
const double PI=acos(-1);
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<0||c>9){if(c==-)f=-1;c=getchar();}
    while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
    return x*f;
}
struct Vector{
    double x,y;
    Vector(double a=0,double b=0):x(a),y(b){}
};
typedef Vector CD;
Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator *(Vector a,Vector b){return Vector(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
Vector operator /(Vector a,int b){return Vector(a.x/b,a.y/b);}
Vector conj(Vector a){return Vector(a.x,-a.y);}

struct FastFourierTransform{
    int n;
    CD omega[N],omegaInv[N];
    void ini(){
        for(int k=0;k<n;k++) 
            omega[k]=CD(cos(2*PI/n*k),sin(2*PI/n*k)),
            omegaInv[k]=conj(omega[k]);
    }
    void transform(CD *a,CD *omega){
        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            if(t<i) swap(a[t],a[i]);
        }

        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            for(CD *p=a;p!=a+n;p+=l)
                for(int k=0;k<m;k++){
                    CD t=omega[n/l*k]*p[k+m];
                    p[k+m]=p[k]-t;
                    p[k]=p[k]+t;
                }
        }
    }
    void DFT(CD *a,int flag){
        if(flag==1) transform(a,omega);
        else{
            transform(a,omegaInv);
            for(int i=0;i<n;i++) a[i]=a[i]/n;
        }
    }
    void MUL(CD *a,CD *b,int m){
        n=1;
        while(n<m) n<<=1;
        ini();
        DFT(a,1);DFT(b,1);
        for(int i=0;i<n;i++) a[i]=a[i]*b[i];
        DFT(a,-1);
    }
}fft;
int n;
CD A[N],B[N];
int main(){
    freopen("in","r",stdin);
    n=read();
    for(int i=0;i<n;i++) scanf("%lf%lf",&A[n-i-1].x,&B[i].x);
    fft.MUL(A,B,n+n-1);
    for(int i=0;i<n;i++) printf("%d\n",int(A[n-i-1].x+0.5));
}

 

 

 

BZOJ 2194 [快速傅里叶变换 卷积]

原文:http://www.cnblogs.com/candy99/p/6388946.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!