【CTSC2010】【洛谷P4191】性能优化(混合基FFT求循环卷积)

传送门


Solution:

听说Bluestein过不了这道题就没写,去学习了一下混合基FFT。

一般的FFT是把每一层分为两组,然后利用单位根的性质来合并。

由于题目说了 n n n能够拆分成非常小的质因数,于是就能够用混合基FFT了,其实就是把分治策略改一下就行了。

推荐唐jz神仙的题解:https://blog.csdn.net/skywalkert/article/details/51737272


代码:

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const
#define gc get_char

namespace IO{
	inline char get_char(){
		static cs int Rlen=1<<22|1;
		static char buf[Rlen],*p1,*p2;
		return (p1==p2)&&(p2=(p1=buf)+fread(buf,1,Rlen,stdin),p1==p2)?EOF:*p1++;
	}
	
	template<typename T>
	inline T get(){
		re char c;
		while(!isdigit(c=gc()));re T num=c^48;
		while(isdigit(c=gc()))num=(num+(num<<2)<<1)+(c^48);
		return num;
	}
	inline int getint(){return get<int>();}
}
using namespace IO;

using std::cerr;
using std::cout;

int mod;
inline int add(int a,int b){return (a+b>=mod)?a+b-mod:a+b;}
inline int dec(int a,int b){return (a<b)?a-b+mod:a-b;}
inline void Inc(int &a,int b){(a+=b)>=mod?(a-=mod):a;}
inline void Dec(int &a,int b){(a-=b)<0?a+=mod:a;}
inline int mul(int a,int b){return (ll)a*b>=mod?(ll)a*b%mod:(ll)(a*b);}
inline int power(int a,int b,int res=1){
	while(b){
		if(b&1)res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}

int prime[100],pcnt;
inline void factor(int n){
	for(int re i=2;i*i<=n;++i)
	while(n%i==0)n/=i,prime[++pcnt]=i;
	if(n>1)prime[++pcnt]=n;
}

inline int get_g(int n){
	factor(n);
	for(int re i=2;;++i){
		bool flag=true;
		for(int re j=1;j<=pcnt;++j)if(power(i,n/prime[j])==1){
			flag=false;break;
		}
		if(flag)return i;
	}
}

cs int N=5e5+5;

int n,m;
int tp[N]; 

inline void reverse(int *A){
	for(int re i=pcnt,block=n;i;--i){
		for(int re idx=0,k=0;k<n;k+=block)
		for(int re j=0;j<prime[i];++j)
		for(int re l=0;l<block;l+=prime[i])
		tp[idx++]=A[k+j+l];
		for(int re k=0;k<n;++k)A[k]=tp[k];
		block/=prime[i];
	}
}

int pos[N];
int wn[N];
inline void DFT(int *A,int typ){
	for(int re i=0;i<n;++i)tp[i]=A[pos[i]];
	for(int re i=0;i<n;++i)A[i]=tp[i];
	for(int re i=1,block=1;i<=pcnt;++i){
		cs int mid=block,wi=wn[n/(block*=prime[i])];
		memset(tp,0,sizeof(int)*n);
		
		for(int re j=0;j<n;j+=block){
			int wk=1;
			for(int re k=0;k<block;++k){
				for(int re l=k%mid,w=1;l<block;l+=mid,w=mul(w,wk))
				Inc(tp[j+k],mul(w,A[j+l]));
				wk=mul(wk,wi);
			}
		}
		for(int re j=0;j<n;++j)A[j]=tp[j];
	}
	if(typ==-1){
		std::reverse(A+1,A+n);
		for(int re i=0,inv=power(n,mod-2);i<n;++i)A[i]=mul(A[i],inv);
	}
}

int a[N],b[N];
signed main(){
//	freopen("program.in","r",stdin);
	n=getint(),m=getint();mod=n+1;
	wn[0]=1;wn[1]=get_g(n);
	for(int re i=0;i<n;++i)a[i]=getint();
	for(int re i=0;i<n;++i)b[i]=getint();
	for(int re i=2;i<n;++i)wn[i]=mul(wn[i-1],wn[1]);
	for(int re i=1;i<n;++i)pos[i]=i;reverse(pos);
	DFT(a,1),DFT(b,1);
	for(int re i=0;i<n;++i)a[i]=mul(a[i],power(b[i],m));
	DFT(a,-1);
	for(int re i=0;i<n;++i)cout<<a[i]<<"\n";
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值