前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >拆系数FFT

拆系数FFT

作者头像
饶文津
发布2020-06-02 15:55:09
8950
发布2020-06-02 15:55:09
举报
文章被收录于专栏:饶文津的专栏饶文津的专栏

1.基本介绍

代码语言:javascript
复制
//495ms
#include <bits/stdc++.h>
#define rep(i,l,r) for(int i=l,ed=r;i<ed;++i)
typedef long long ll;
const double PI = acos(-1);
const int N = 1<<20;
const int BUF_SIZE=33554431;
using namespace std;

struct buf{
    char a[BUF_SIZE],b[BUF_SIZE],*s,*t;
    buf():s(a),t(b){a[fread(a,1,sizeof a,stdin)]=0;}
    ~buf(){fwrite(b,1,t-b,stdout);}
    operator int(){
        int x=0;
        while(*s<48)++s;
        while(*s>32)
            x=x*10+*s++-48;
        return x;
    }
    void out(int x){
        static char c[12];
        char*i=c;
        if(!x)*t++=48;
        else{
            while(x){
                int y=x/10;
                *i++=x-y*10+48,x=y;
            }
            while(i!=c)*t++=*--i;
        }
        *t++=10;
    }
}it;
struct cp{
    double x,y;
    cp(double _x=0,double _y=0):x(_x),y(_y){}
    cp operator +(const cp&amp; b)const{return cp(x+b.x,y+b.y);}
    cp operator -(const cp&amp; b)const{return cp(x-b.x,y-b.y);}
    cp operator *(const cp&amp; b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
    cp operator !()const{return cp(x,-y);}
}w[N];
void fft(cp p[],int n){
    for(int i=0,j=0;i<n;++i){
        if(i>j)swap(p[i],p[j]);
        for(int l=n>>1;(j^=l)<l;l>>=1);
    }
    for(int i=2;i<=n;i<<=1)
    for(int j=0,m=i>>1;j<n;j+=i)
        rep(k,0,m){
            cp b=w[n/i*k]*p[j+m+k];
            p[j+m+k]=p[j+k]-b;
            p[j+k]=p[j+k]+b;
        }
}
void conv(int n,ll *x,ll *y,ll *z){
    static cp p[N],q[N],h(0,-0.25);
    rep(i,0,n){
        w[i]=cp(cos(2*PI*i/n),sin(2*PI*i/n));
        p[i]=cp(x[i],y[i]);
    }
    fft(p,n);
    rep(i,0,n){
        int j=i?(n-i):0;
        q[j]=(p[i]*p[i]-!p[j]*!p[j])*h;
    }
    fft(q,n);
    rep(i,0,n)z[i]=q[i].x/n+0.5;
}
int n,m,p;
ll a[N],b[N],c[N];
int main(){
    n=it+1;m=it+1;
    rep(i,0,n) a[i]=it;
    rep(i,0,m) b[i]=it;
    for(n+=m-1,p=1;p<n;p<<=1);
    conv(p,a,b,c);
    rep(i,0,n)it.out(c[i]);
    return 0;
}

2.更快的卷积

代码语言:javascript
复制
//325ms
#include <bits/stdc++.h>
#define rep(i,l,r) for(int i=l,ed=r;i<ed;++i)
typedef long long ll;
const double PI = acos(-1);
const int N = 1<<20;
const int BUF_SIZE=33554431;
using namespace std;

struct buf{
    char a[BUF_SIZE],b[BUF_SIZE],*s,*t;
    buf():s(a),t(b){a[fread(a,1,sizeof a,stdin)]=0;}
    ~buf(){fwrite(b,1,t-b,stdout);}
    operator int(){
        int x=0;
        while(*s<48)++s;
        while(*s>32)
            x=x*10+*s++-48;
        return x;
    }
    void out(int x){
        static char c[12];
        char*i=c;
        if(!x)*t++=48;
        else{
            while(x){
                int y=x/10;
                *i++=x-y*10+48,x=y;
            }
            while(i!=c)*t++=*--i;
        }
        *t++=10;
    }
}it;
struct cp{
    double x,y;
    cp(double _x=0,double _y=0):x(_x),y(_y){}
    cp operator +(const cp&amp; b)const{return cp(x+b.x,y+b.y);}
    cp operator -(const cp&amp; b)const{return cp(x-b.x,y-b.y);}
    cp operator *(const cp&amp; b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
    cp operator *(double b)const{return cp(b*x,b*y);}
    cp operator !()const{return cp(x,-y);}
}w[N];
void fft(cp *p,int n){
    for(int i=0,j=0;i<n;++i){
        if(i>j)swap(p[i],p[j]);
        for(int l=n>>1;(j^=l)<l;l>>=1);
    }
    for(int i=2;i<=n;i<<=1)
    for(int j=0,m=i>>1;j<n;j+=i)
        rep(k,0,m){
            cp b=w[n/i*k]*p[j+m+k];
            p[j+m+k]=p[j+k]-b;
            p[j+k]=p[j+k]+b;
        }
}
void conv(int n,ll *x,ll *y,ll *z){
    static cp p[N],q[N],a[N];
    rep(i,0,n){
        (i&amp;1?p[i>>1].y:p[i>>1].x)=x[i];
        (i&amp;1?q[i>>1].y:q[i>>1].x)=y[i];
    }
    rep(i,0,n>>=1)w[i]=cp(cos(2*PI*i/n),sin(2*PI*i/n));
    fft(p,n);fft(q,n);
    rep(i,0,n){
        int j=i?n-i:0;
        a[j]=p[i]*q[i]-((cp(1,0)+w[i])*(p[i]-!p[j])*(q[i]-!q[j]))*0.25;
    }
    fft(a,n);
    rep(i,0,n)z[i<<1]=a[i].x/n+0.5,z[i<<1|1]=a[i].y/n+0.5;
}
int n,m,p;
ll a[N],b[N],c[N];
int main(){
    n=it+1;m=it+1;
    rep(i,0,n) a[i]=it;
    rep(i,0,m) b[i]=it;
    for(n+=m-1,p=2;p<n;p<<=1);
    conv(p,a,b,c);
    rep(i,0,n)it.out(c[i]);
    return 0;
}

3.拆系数FFT

代码语言:javascript
复制
//933ms
#include <bits/stdc++.h>
#define rep(i,l,r) for(int i=l,ed=r;i<ed;++i)
typedef long long ll;
const double PI = acos(-1);
const int N = 1<<20;
const ll mod = 1e9+7;
const int BUF_SIZE=33554431;
using namespace std;

struct buf{
    char a[BUF_SIZE],b[BUF_SIZE],*s,*t;
    buf():s(a),t(b){a[fread(a,1,sizeof a,stdin)]=0;}
    ~buf(){fwrite(b,1,t-b,stdout);}
    operator int(){
        int x=0;
        while(*s<48)++s;
        while(*s>32)
            x=x*10+*s++-48;
        return x;
    }
    void out(int x){
        static char c[12];
        char*i=c;
        if(!x)*t++=48;
        else{
            while(x){
                int y=x/10;
                *i++=x-y*10+48,x=y;
            }
            while(i!=c)*t++=*--i;
        }
        *t++=10;
    }
}it;
struct cp{
    double x,y;
    cp(double _x=0,double _y=0):x(_x),y(_y){}
    cp operator +(const cp&amp; b)const{return cp(x+b.x,y+b.y);}
    cp operator -(const cp&amp; b)const{return cp(x-b.x,y-b.y);}
    cp operator *(const cp&amp; b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
    cp operator !()const{return cp(x,-y);}
}w[N];
void fft(cp p[],int n){
    for(int i=0,j=0;i<n;++i){
        if(i>j)swap(p[i],p[j]);
        for(int l=n>>1;(j^=l)<l;l>>=1);
    }
    for(int i=2;i<=n;i<<=1)
    for(int j=0,m=i>>1;j<n;j+=i)
        rep(k,0,m){
            cp b=w[n/i*k]*p[j+m+k];
            p[j+m+k]=p[j+k]-b;
            p[j+k]=p[j+k]+b;
        }
}
void conv(int n,ll *x,ll *y,ll *z){
    static cp p[N],q[N],a[N],b[N],c[N],d[N];
    static cp r(0.5,0),h(0,-0.5),o(0,1);
    rep(i,0,n){
        w[i]=cp(cos(2*PI*i/n),sin(2*PI*i/n));
        x[i]=(x[i]+mod)%mod,y[i]=(y[i]+mod)%mod;
        p[i]=cp(x[i]>>15,x[i]&amp;32767),q[i]=cp(y[i]>>15,y[i]&amp;32767);
    }
    fft(p,n);fft(q,n);
    rep(i,0,n){
        int j=i?(n-i):0;
        static cp ka,ba,kb,bb;
        ka=(p[i]+!p[j])*r;
        ba=(p[i]-!p[j])*h;
        kb=(q[i]+!q[j])*r;
        bb=(q[i]-!q[j])*h;
        a[j]=ka*kb;b[j]=ka*bb;
        c[j]=kb*ba;d[j]=ba*bb;
    }
    rep(i,0,n){
        p[i]=a[i]+b[i]*o;
        q[i]=c[i]+d[i]*o;
    }
    fft(p,n);fft(q,n);
    rep(i,0,n){
        ll a,b,c,d;
        a=(ll)(p[i].x/n+0.5)%mod;
        b=(ll)(p[i].y/n+0.5)%mod;
        c=(ll)(q[i].x/n+0.5)%mod;
        d=(ll)(q[i].y/n+0.5)%mod;
        z[i]=((a<<30)+((b+c)<<15)+d)%mod;
    }
}
int n,m,p;
ll a[N],b[N],c[N];
int main(){
    n=it+1;m=it+1;
    rep(i,0,n) a[i]=it;
    rep(i,0,m) b[i]=it;
    for(n+=m-1,p=1;p<n;p<<=1);
    conv(p,a,b,c);
    rep(i,0,n)it.out((c[i]+mod)%mod);
    return 0;
}
本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2018-09-18 ,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1.基本介绍
  • 2.更快的卷积
  • 3.拆系数FFT
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档