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);
    }


  • 相关阅读:
    机器学习-初学者入门
    安装.cer证书并将证书从.cer格式转化为.pem格式
    字符串反转C#的实现
    Linux系统下远程文件拷贝scp命令
    【Django】ESRTful APi
    数据结构-栈跟队列基础部分
    数据结构-排序
    数据分析--Matplotlib的基本使用
    数据分析--pandas的基本使用
    数据分析--numpy的基本使用
  • 原文地址:https://www.cnblogs.com/Extended-Ash/p/8449237.html
Copyright © 2011-2022 走看看