其他分享
首页 > 其他分享> > fft && ntt 模板

fft && ntt 模板

作者:互联网

fft参考blog
ntt参考blog

FFT模板
点击查看代码
#include<bits/stdc++.h>
using namespace std;
template<typename T>
inline void read(T &x){
	x=0;T fl=1;char tmp=getchar();
	while(tmp<'0'||tmp>'9')fl=tmp=='-'?-fl:fl,tmp=getchar();
	while(tmp>='0'&&tmp<='9')x=(x<<1)+(x<<3)+tmp-'0',tmp=getchar();
	x=x*fl;
}
const double Pi=acos(-1);
const int maxn=4e6;
struct FastFourierTransform{
	complex<double>omega[maxn],iomega[maxn];
	
	void init(const int n){
		for(int i=0;i<n;i++){
			omega[i]=complex<double>(cos(2*Pi/n*i),sin(2*Pi/n*i));
			iomega[i]=conj(omega[i]);
		}
	}
	
	void transform(complex<double>*a,int n,complex<double> *omega){
		int k=0;
		while((1<<k)<n)k++;
		for(int i=0;i<n;i++){
			int t=0;
			for(int j=0;j<k;j++)
				if(i&(1<<j))t|=1<<k-j-1;
			if(t<i)swap(a[i],a[t]);
		}
		for(int l=2;l<=n;l<<=1){
			int m=l/2;
			for(complex<double> *p=a;p!=a+n;p+=l){
				for(int i=0;i<m;i++){
					complex<double> t=omega[n/l*i]*p[i+m];
					p[i+m]=p[i]-t;
					p[i]=p[i]+t;
				}
			}
		}
	}
	
	void dft(complex<double>*a,const int n){
		transform(a,n,omega);
	}
	
	void idft(complex<double>*a,const int n){
		transform(a,n,iomega);
		for(int i=0;i<n;i++)
			a[i]/=n;
	}
}fft;
inline int multiply(int *a1,int n1,int *a2,const int n2){
	int n=1;
	while(n<n1+n2)n<<=1;
	complex<double>c1[maxn],c2[maxn];
	for(int i=0;i<n1;i++)
		c1[i]=a1[i];
	for(int i=0;i<n2;i++)
		c2[i]=a2[i];
	fft.init(n);
	fft.dft(c1,n),fft.dft(c2,n);
	for(int i=0;i<n;i++) c1[i]=c1[i]*c2[i];
	fft.idft(c1,n);
	n1=n1+n2-1;
	for(int i=0;i<n1+n2-1;i++) a1[i]=floor(c1[i].real()+0.5);
	return n1;
}
int n,m;
int a[maxn],b[maxn];
signed main(){
	cin>>n>>m;n++,m++;
	for(int i=0;i<n;i++)
		read(a[i]);
	for(int i=0;i<m;i++)
		read(b[i]);
	n=multiply(a,n,b,m);
	for(int i=0;i<n;i++)
		printf("%d ",a[i]);
	return 0;
}
NTT模板

求原根,可以先求出p-1的质因子。由于原根往往不大,可以从2到p-1依次枚举i,

\(若i^{(p-1)/p_i}≠1(mod p),{\forall}p_i为(p-1)的质因数,则i为p的一个原根。\)

如此将复数单位根替换成模p意义下的

点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=998244353;
const int maxn=2.5e6;
template<typename T>
inline void read(T &x){
	x=0;T fl=1;char tmp=getchar();
	while(tmp<'0'||tmp>'9')fl=tmp=='-'?-fl:fl,tmp=getchar();
	while(tmp>='0'&&tmp<='9')x=(x<<1)+(x<<3)+tmp-'0',tmp=getchar();
	x=x*fl;
}
inline ll pw(ll x,ll n,ll p){
	ll ans=1;
	while(n){
		if(n&1)ans=ans*x%p;
		x=x*x%p,n>>=1;
	}
	return ans;
}
inline ll root(const ll p){
	ll pri[60],cnt=0;
	ll x=p-1;
	for(int k=2;k*k<=p-1;k++){
		if(x%k==0){
			pri[++cnt]=k;
			while(x%k==0)x/=k;
		}
	}
	if(x>1)pri[++cnt]=x;
	int fl;
	for(int i=2;i<=p;i++){
		fl=0;
		for(int j=1;j<=cnt;j++){
			if(pw(i,(p-1)/pri[j],p)==1){
				fl=1;
				break;
			}
		}
		if(!fl)return i;
	}
	throw;
}
inline void exgcd(const ll a,const ll b,ll &x, ll &y){
	if(!b)x=1,y=0;
	else exgcd(b,a%b,y,x),y-=a/b*x;
}
inline ll inv(const ll a,const ll p){//must exist
	ll x,y;
	exgcd(a,p,x,y);
	return (x%p+p)%p;
	
}
struct NumberTheoreticTransform{
	ll omega[maxn],iomega[maxn];
	
	void init(const int n){
		ll g=root(mod),x=pw(g,(mod-1)/n,mod);
		omega[0]=iomega[0]=1;
		for(int i=1;i<n;i++){
			omega[i]=omega[i-1]*x;
			iomega[i]=inv(omega[i],mod);
		}
	}
	
	void transform(ll *a,const int n,ll *omega){
		int k=0;
		while((1<<k)<n)k++;
		for(int i=0;i<n;i++){
			int t=0;
			for(int j=0;j<k;j++)if(i&(1<<j))t|=1<<k-j-1;
			if(i>t)swap(a[i],a[t]);
		}
		for(int l=2;l<=n;l++){
			int m=l/2;
			for(ll *p=a;p!=a+n;p+=l)
				for(int i=0;i<m;i++){
					ll t=omega[n/l*i]*p[i+m]%mod;
					p[i+m]=(p[i]-t+mod)%mod;
					p[i]=(p[i]+t)%mod;
				}
		}
	}
	
	void dft(ll *a,const int n){
		transform(a,n,omega);
	}
	
	void idft(ll *a,const int n){
		transform(a,n,iomega);
		ll x=inv(n,mod);
		for(int i=0;i<n;i++)
			a[i]=a[i]*x%mod;
	}
}ntt;

标签:tmp,const,int,ll,fft,ntt,&&,fl,void
来源: https://www.cnblogs.com/xyc1719/p/16533688.html