ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

fft && ntt 模板

2022-07-30 00:06:00  阅读:175  来源: 互联网

标签:tmp const int ll fft ntt && fl void


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

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有