zoukankan      html  css  js  c++  java
  • (各种)FFT模板

    先来一个标准的归并版FFT 2881ms

    #include<math.h>
    #include<stdio.h>
    #include<string.h>
    #include<algorithm>
    #define D double
    using namespace std;
    struct Z{ D x,y; } a[280010],b[280010],z[280010];
    inline Z operator+ (Z x,Z y){ return (Z){x.x+y.x,x.y+y.y}; }
    inline Z operator- (Z x,Z y){ return (Z){x.x-y.x,x.y-y.y}; }
    inline Z operator* (Z x,Z y){ return (Z){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x}; }
    void FFT(Z* a,Z* b,int l,int r,int g){
    	if(l==r) return;
    	int m=l+r>>1,n=r-l+1,i,j=l,k=m+1;
    	Z w=(Z){cos(2*M_PI/n),g*sin(2*M_PI/n)},z=(Z){1,0};
    	for(i=l;i<=r;++i) b[~i&1?j++:k++]=a[i];
    	FFT(b,a,l,m,g); FFT(b,a,m+1,r,g);
    	for(i=l,j=m-l+1;i<=m;++i,z=z*w){
    		a[i]=b[i]+z*b[i+j];
    		a[i+j]=b[i]-z*b[i+j];
    	}
    }
    int n,m,M=1,s[280010];
    int main(){
    	scanf("%d%d",&n,&m);
    	for(int i=0;i<=n;++i) scanf("%lf",&a[i].x);
    	for(int i=0;i<=m;++i) scanf("%lf",&b[i].x);
    	while(M<=max(n,m)) M<<=1; M<<=1;
    	FFT(a,z,0,M-1,1);
    	FFT(b,z,0,M-1,1);
    	for(int i=0;i<M;++i) a[i]=a[i]*b[i];
    	FFT(a,z,0,M-1,-1);
    	for(int i=0;i<M;++i) s[i]+=int(a[i].x/M+0.5);
    	for(int i=0;i<=n+m;++i) printf("%d ",s[i]);
    }
    让后就可以写一个蝴蝶变换版本的 1097ms

    #include<math.h>
    #include<stdio.h>
    #include<string.h>
    #include<algorithm>
    #define D double
    using namespace std;
    struct Z{ D x,y; } a[280010],b[280010],z[280010];
    inline Z operator+ (Z x,Z y){ return (Z){x.x+y.x,x.y+y.y}; }
    inline Z operator- (Z x,Z y){ return (Z){x.x-y.x,x.y-y.y}; }
    inline Z operator* (Z x,Z y){ return (Z){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x}; }
    inline Z rt(double X){ return (Z){cos(X),sin(X)}; }
    void FFT(Z* a,Z* b,int n,int g){
    	for(int i=0,j,k,t;i<n;++i){
    		for(j=0,k=i,t=n-1;t;t>>=1,k>>=1) j=(j<<1)|(k&1);
    		b[j]=a[i];
    	}
    	for(int m=2;m<=n;m<<=1){
    		Z w=rt(g*M_PI/m),z=(Z){1,0},u,v;
    		for(int i=0,k=m>>1;i<k;++i){
    			for(int j=i;j<n;j+=m){
    				u=b[j]; v=b[j+k]*z;
    				b[j]=u+v; b[j+k]=u-v;
    			}
    			z=z*w;
    		}
    	}
    	for(int i=0;i<n;++i) a[i]=b[i];
    }
    int n,m,M=1,s[280010];
    int main(){
    	scanf("%d%d",&n,&m);
    	for(int i=0;i<=n;++i) scanf("%lf",&a[i].x);
    	for(int i=0;i<=m;++i) scanf("%lf",&b[i].x);
    	while(M<=max(n,m)) M<<=1; M<<=1;
    	FFT(a,z,M,2);
    	FFT(b,z,M,2);
    	for(int i=0;i<M;++i) a[i]=a[i]*b[i];
    	FFT(a,z,M,-2);
    	for(int i=0;i<M;++i) s[i]+=int(a[i].x/M+0.5);
    	for(int i=0;i<=n+m;++i) printf("%d ",s[i]);
    }
    加速蝴蝶变换版本的 620ms

    #include<math.h>
    #include<stdio.h>
    #include<string.h>
    #include<algorithm>
    #define D double
    using namespace std;
    int n,m,M=1,L,s[280010],r[280010];
    struct Z{ D x,y; } a[280010],b[280010],z[280010];
    inline Z rt(double X){ return (Z){cos(X),sin(X)}; }
    inline Z operator+ (Z x,Z y){ return (Z){x.x+y.x,x.y+y.y}; }
    inline Z operator- (Z x,Z y){ return (Z){x.x-y.x,x.y-y.y}; }
    inline Z operator* (Z x,Z y){ return (Z){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x}; }
    void FFT(Z* a,Z* b,int n,int g){
    	for(int i=0;i<n;++i) b[r[i]]=a[i];
    	for(int m=2;m<=n;m<<=1){
    		Z w=rt(g*M_PI/m),u,v,z;
    		for(int i,k=m>>1,j=0;j<n;j+=m){
    			for(z=(Z){1,i=0};i<k;++i,z=z*w){
    				u=b[i+j]; v=z*b[i+j+k];
    				b[i+j]=u+v; b[i+j+k]=u-v;
    			}
    		}
    	}
    	memcpy(a,b,n*sizeof(Z));
    }
    int main(){
    	scanf("%d%d",&n,&m);
    	for(int i=0;i<=n;++i) scanf("%lf",&a[i].x);
    	for(int i=0;i<=m;++i) scanf("%lf",&b[i].x);
    	for(;M<=max(n,m);++L) M<<=1; M<<=1;
    	for(int i=0;i<M;++i) r[i]=(r[i>>1]>>1)|((i&1)<<L);
    	FFT(a,z,M,2);
    	FFT(b,z,M,2);
    	for(int i=0;i<M;++i) a[i]=a[i]*b[i];
    	FFT(a,z,M,-2);
    	for(int i=0;i<M;++i) s[i]+=int(a[i].x/M+0.5);
    	for(int i=0;i<=n+m;++i) printf("%d ",s[i]);
    }

    最终版本:490ms

    #include<cctype>
    #include<math.h>
    #include<stdio.h>
    #include<string.h>
    #include<algorithm>
    #define D double
    #define Rg register
    #define N 280010
    using namespace std;
    int n,m,M=1,L,s[N],r[N];
    struct Z{ D x,y; } a[N],b[N],z[N];
    inline Z rt(D X){ return (Z){cos(X),sin(X)}; }
    inline Z operator+ (Z x,Z y){ return (Z){x.x+y.x,x.y+y.y}; }
    inline Z operator- (Z x,Z y){ return (Z){x.x-y.x,x.y-y.y}; }
    inline Z operator* (Z x,Z y){ return (Z){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x}; }
    inline int rd(){  
    	int X=0,w=0; char ch=0;
        while(!isdigit(ch)) {w|=ch=='-';ch=getchar();}
        while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
        return w?-X:X; 
    } 
    void FFT(Z* a,Z* b,int n,int g){
    	for(int i=0;i<n;++i) b[r[i]]=a[i];
    	for(int m=2;m<=n;m<<=1){
    		Rg Z w=rt(g*M_PI/m),u,v,z;
    		for(Rg int i,k=m>>1,j=0;j<n;j+=m){
    			for(z=(Z){1,i=0};i<k;++i,z=z*w){
    				u=b[i+j]; v=z*b[i+j+k];
    				b[i+j]=u+v; b[i+j+k]=u-v;
    			}
    		}
    	}
    	memcpy(a,b,n*sizeof(Z));
    }
    int main(){
    	n=rd(); m=rd();
    	for(int i=0;i<=n;++i) a[i].x=rd();
    	for(int i=0;i<=m;++i) b[i].x=rd();
    	for(;M<=max(n,m);++L) M<<=1; M<<=1;
    	for(int i=0;i<M;++i) r[i]=(r[i>>1]>>1)|((i&1)<<L);
    	FFT(a,z,M,2);
    	FFT(b,z,M,2);
    	for(int i=0;i<M;++i) a[i]=a[i]*b[i];
    	FFT(a,z,M,-2);
    	for(int i=0;i<M;++i) s[i]+=int(a[i].x/M+0.5);
    	for(int i=0;i<=n+m;++i) printf("%d ",s[i]);
    }

    数论版(NTT) 649ms

    #include<cctype>
    #include<math.h>
    #include<stdio.h>
    #include<string.h>
    #include<algorithm>
    #define N 280010
    #define Rg register
    #define LL long long
    #define M 1004535809 
    using namespace std;
    int n,m,T,L,r[N],a[N],b[N],z[N]; LL W[N],inv;
    inline LL pow(LL x,LL k,LL& s){
    	for(s=1;k;x=x*x%M,k>>=1) k&1?s=s*x%M:0;
    }
    inline int rd(){  
    	int X=0; char ch=0;
        while(!isdigit(ch)) ch=getchar();
        while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
        return X; 
    }
    void NTT(int* a,int* b,int n,int g){
    	for(int i=0;i<n;++i) b[r[i]]=a[i];
    	for(int m=2;m<=n;m<<=1){
    		LL w=W[(g?(n/m):(n-n/m))],u,v,z;
    		for(int i,k=m>>1,j=0;j<n;j+=m)
    			for(z=1,i=0;i<k;++i,z=z*w%M){
    				u=b[i+j]; v=z*b[i+j+k]%M;
    				b[i+j]=(u+v)%M; b[i+j+k]=(u-v+M)%M;
    			}
    	}
    	memcpy(a,b,n<<2);
    }
    int main(){
    	n=rd(); m=rd();
    	for(int i=0;i<=n;++i) a[i]=rd();
    	for(int i=0;i<=m;++i) b[i]=rd();
    	T=n+m; m=max(n,m); n=1;
    	for(;n<=m;++L) n<<=1; n<<=1;
    	for(int i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<L);
    	pow(3,(M-1)/n,W[*W=1]);
    	for(int i=2;i<=n;++i) W[i]=W[i-1]*W[1]%M;
    	NTT(a,z,n,1);
    	NTT(b,z,n,1);
    	for(int i=0;i<n;++i) a[i]=(LL)a[i]*b[i]%M;
    	NTT(a,z,n,0);
    	pow(n,M-2,inv);
    	for(int i=0;i<=T;++i) printf("%lld ",a[i]*inv%M);
    }


  • 相关阅读:
    sfs2x 连接 mongodb
    java websocket
    webstorm 4.0 注册码
    解决 sfs2 admin tool 找不到扩展
    window 注册表五大类
    opengl 学习第二日
    java google Protobuf
    扩展 java sencha touch PhonegapPlugin
    sencha touch2 kryonet socket phonegap 通信 作者:围城
    sencha touch2 layout 笔记
  • 原文地址:https://www.cnblogs.com/Extended-Ash/p/8449237.html
Copyright © 2011-2022 走看看