zoukankan      html  css  js  c++  java
  • 分治FFT

    分治FFT

    概述

    严格的分治FFT为“一般形式”,而我们把所有带有FFT的分治都称为分治FFT,所以这个名称并没有什么意义。
    分治FFT代码比较复杂,低阶的比赛里应该不会出现。

    一般形式

    严格的分治FFT负责解决类似(displaystyle b_k=sum^{k-1}_{i=1}{b_ia_{k-i}})的数列,其中整个式子的形状符合卷积,且所求数列的每一项都依赖于前面的项。

    已知数组(a),求数组(b),其中(displaystyle b_k=sum^{k-1}_{i=1}{b_ia_{k-i}})

    考虑分治,对于区间([l,r)),先计算([l,m))的答案。对于(kgeq m),将(b_k)拆成(displaystyle sum^{m-1}_{i=0}b_ia_{k-i})(displaystyle sum^{k-1}_{i=m}b_ia_{k-i})两部分,后者可以给数组下标添加偏移量之后递归计算,而前者可以通过普通的卷积计算。

    总时间复杂度为(Theta(n log^2 n))

    实现

    假设我们求出了(l o mid)的答案,要求这些点对(mid+1 o r)的影响,那么对右半边点(x)的贡献为:(displaystyle w_x=sum_{i=l}^{mid} f[i] * g[x-i])这部分可以利用卷积来快速计算。计算完以后,答案直接加到答案数组就可以了。需要注意的是,如果要求左边点对右边点的影响,首先整个区间以左对该区间的贡献应该先求出。所以分治过程为先分治左边,在求出中间,然后在递归右边。

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long LL;
    const int MAXN=1e6+10;
    const LL MOD=998244353,G=3,iG=332748118;
    inline LL fpm(LL base,LL p){
    	LL ret=1;
    	while(p){
    		if(p&1)
    			ret=ret*base%MOD;
    		base=base*base%MOD;
    		p>>=1;
    	}
    	return ret;
    }
    void exgcd(LL a,LL b,LL &x,LL &y){
    	if(!b){
    		x=1;
    		y=0;
    		return;
    	}
    	exgcd(b,a%b,y,x);
    	y-=x*(a/b);
    }
    inline LL inv(LL x){
    	LL ii,jj;
    	exgcd(x,MOD,ii,jj);
    	return (ii%MOD+MOD)%MOD;
    }
    int lim,rev[MAXN],N,L;
    inline void prec(int l,int r){
    	lim=1;
    	L=0;
    	while(lim<r-l){
    		L++;
    		lim<<=1;
    	}
    	for(int i=0;i<lim;i++)
    		rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
    }
    inline void NTT(LL *a,int type){
    	for(int i=0;i<lim;i++)
    		if(i<rev[i])
    			swap(a[i],a[rev[i]]);
    	for(int hf=1;hf<lim;hf<<=1){
    		int len=hf<<1;
    		LL Wn=fpm(type==1?G:iG,(MOD-1)/len);
    		for(int j=0;j<lim;j+=len){
    			LL w=1;
    			for(int k=0;k<hf;k++){
    				LL t1=a[j+k],t2=a[j+k+hf]*w%MOD;
    				a[j+k]=(t1+t2)%MOD;
    				a[j+k+hf]=(t1-t2)%MOD;
    				w=w*Wn%MOD;
    			}
    		}
    	}
    }
    LL f[MAXN],g[MAXN],A[MAXN],B[MAXN];
    inline void cdq(int l,int r){
    	if(l+1>=r||l>=N)
    		return;
    	int mid=(l+r)>>1;
    	cdq(l,mid);
    	prec(l,r);
    	memcpy(A,f+l,((r-l)/2)<<3);
    	memcpy(B,g,(r-l)<<3);
    	memset(A+(r-l)/2,0,((r-l)/2)<<3);
    	NTT(A,1);
    	NTT(B,1);
    	for(int i=0;i<r-l;i++)
    		A[i]=A[i]*B[i]%MOD;
    	NTT(A,-1);
    	int ii=inv(r-l);
    	for(int i=0;i<r-l;i++)
    		A[i]=A[i]*ii%MOD;
    	for(int i=(r-l)>>1;i<r-l;i++)
    		f[i+l]=(f[i+l]+A[i])%MOD;
    	cdq(mid,r);
    }
    int main(){
    	scanf("%d",&N);
    	for(int i=1;i<N;i++)
    		scanf("%lld",g+i);
    	prec(0,N);
    	f[0]=1;
    	cdq(0,lim);
    	for(int i=0;i<N;i++)
    		printf("%lld ",(f[i]+MOD)%MOD);
    	return 0;
    }
    

    例题

    P4721 【模板】分治 FFT

    见“实现”

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long LL;
    const int MAXN=1e6+10;
    const LL MOD=998244353,G=3,iG=332748118;
    inline LL fpm(LL base,LL p){
    	LL ret=1;
    	while(p){
    		if(p&1)
    			ret=ret*base%MOD;
    		base=base*base%MOD;
    		p>>=1;
    	}
    	return ret;
    }
    void exgcd(LL a,LL b,LL &x,LL &y){
    	if(!b){
    		x=1;
    		y=0;
    		return;
    	}
    	exgcd(b,a%b,y,x);
    	y-=x*(a/b);
    }
    inline LL inv(LL x){
    	LL ii,jj;
    	exgcd(x,MOD,ii,jj);
    	return (ii%MOD+MOD)%MOD;
    }
    int lim,rev[MAXN],N,L;
    inline void prec(int l,int r){
    	lim=1;
    	L=0;
    	while(lim<r-l){
    		L++;
    		lim<<=1;
    	}
    	for(int i=0;i<lim;i++)
    		rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
    }
    inline void NTT(LL *a,int type){
    	for(int i=0;i<lim;i++)
    		if(i<rev[i])
    			swap(a[i],a[rev[i]]);
    	for(int hf=1;hf<lim;hf<<=1){
    		int len=hf<<1;
    		LL Wn=fpm(type==1?G:iG,(MOD-1)/len);
    		for(int j=0;j<lim;j+=len){
    			LL w=1;
    			for(int k=0;k<hf;k++){
    				LL t1=a[j+k],t2=a[j+k+hf]*w%MOD;
    				a[j+k]=(t1+t2)%MOD;
    				a[j+k+hf]=(t1-t2)%MOD;
    				w=w*Wn%MOD;
    			}
    		}
    	}
    }
    LL f[MAXN],g[MAXN],A[MAXN],B[MAXN];
    inline void cdq(int l,int r){
    	if(l+1>=r||l>=N)
    		return;
    	int mid=(l+r)>>1;
    	cdq(l,mid);
    	prec(l,r);
    	memcpy(A,f+l,((r-l)/2)<<3);
    	memcpy(B,g,(r-l)<<3);
    	memset(A+(r-l)/2,0,((r-l)/2)<<3);
    	NTT(A,1);
    	NTT(B,1);
    	for(int i=0;i<r-l;i++)
    		A[i]=A[i]*B[i]%MOD;
    	NTT(A,-1);
    	int ii=inv(r-l);
    	for(int i=0;i<r-l;i++)
    		A[i]=A[i]*ii%MOD;
    	for(int i=(r-l)>>1;i<r-l;i++)
    		f[i+l]=(f[i+l]+A[i])%MOD;
    	cdq(mid,r);
    }
    int main(){
    	scanf("%d",&N);
    	for(int i=1;i<N;i++)
    		scanf("%lld",g+i);
    	prec(0,N);
    	f[0]=1;
    	cdq(0,lim);
    	for(int i=0;i<N;i++)
    		printf("%lld ",(f[i]+MOD)%MOD);
    	return 0;
    }
    
  • 相关阅读:
    BZOJ 1057 悬线法求最大01矩阵
    POJ 2248
    SPOJ
    51NOD
    2017-2018 ACM-ICPC, NEERC, Moscow Subregional Contest J. Judging the Trick
    POJ 1379 模拟退火
    POJ 2420 模拟退火
    Frontend 事后诸葛亮
    【Frontend】Alpha Review 展示博客
    ASE19 团队项目 alpha 阶段 Frontend 组 scrum5 记录
  • 原文地址:https://www.cnblogs.com/guoshaoyang/p/11282228.html
Copyright © 2011-2022 走看看