时间戳

首页 > 技术文章 > 【模板】多项式全家桶

【模板】多项式全家桶

BAJimH 2019-03-21 22:07 原文

FFT

struct Z
{
	double x,y;
	Z(double _x=0,double _y=0):x(_x),y(_y){};
};
Z operator +(Z x,Z y) {return Z(x.x+y.x,x.y+y.y);}
Z operator -(Z x,Z y) {return Z(x.x-y.x,x.y-y.y);}
Z operator *(Z x,Z y) {return Z(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
int bit[M],n,m;
Z wi[M+1],a[M],b[M];
void prp()
{
	fo(i,0,M-1) bit[i]=(bit[i>>1]>>1)|((i&1)<<(L-1));
	fo(i,0,M) wi[i]=Z(cos(2*pi*i/M),sin(2*pi*i/M));
}
void DFT(Z *a,int pd)
{
	fo(i,0,M-1) if(i<bit[i]) swap(a[i],a[bit[i]]);
	int half=1,lim=(M>>1);
	Z v;
	for(int m=2;m<=M;half=m,m<<=1,lim>>=1)
	{
		fo(i,0,half-1)
		{
			Z w=(pd==1)?wi[i*lim]:wi[M-i*lim];
			for(int j=i;j<M;j+=m)
			{
				v=w*a[j+half];
				a[j+half]=a[j]-v;
				a[j]=a[j]+v;
			}
		}
	}
	if(pd<0) fo(i,0,M-1) a[i].x/=M;
}
int main()
{
	cin>>n>>m;
	fo(i,0,n) scanf("%lf",&a[i].x);
	fo(i,0,m) scanf("%lf",&b[i].x);
	prp();
	DFT(a,1),DFT(b,1);
	fo(i,0,M-1) a[i]=a[i]*b[i];
	DFT(a,-1);
	fo(i,0,n+m) 
	{
		int v=(int)(a[i].x+0.5);
		printf("%d ",v);
	}
}

NTT

inline LL md(LL x)
{
	return(x<0)?(x+mo):((x>=mo)?x-mo:x);
}
void NTT(LL *a,int pd,int num)
{
	prp(num);
	fo(i,0,num-1) if(i<bit[i]) swap(a[i],a[bit[i]]);
	int lim=num>>1,half=1;
	LL v;
	for(int m=2;m<=num;half=m,m<<=1,lim>>=1)
	{
		fo(i,0,half-1)
		{
			LL w=(pd==1)?wi[i*lim]:wi[num-i*lim];
			for(int j=i;j<num;j+=m)
			{
				v=(w*a[j+half])%mo;
				a[j+half]=md(a[j]-v);
				a[j]=md(a[j]+v);
			}
		}
	}
	if(pd<0) fo(i,0,num-1) a[i]=a[i]*ny[num]%mo;
}

多项式求逆

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstring>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define M 262144
#define mo 998244353
#define LL long long
#define N 100005
using namespace std;
LL a[M+1],b[M+1],c[M+1],wi[M+1],wg[M+1],ny2,d[M+1];
int bit[M+1],n;
LL ksm(LL k,LL n)
{
	LL s=1;
	for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;
	return s;
}
void prp(int num,int L)
{
	ny2=ksm(num,mo-2);
	fo(i,0,num) wi[i]=wg[i*(M/num)];
	fo(i,0,num-1) bit[i]=(bit[i>>1]>>1)|((i&1)<<(L-1));
}
void NTT(LL *a,bool pd,int num)
{
	fo(i,0,num-1) if(bit[i]<i) swap(a[bit[i]],a[i]);
	int half=1,lim=num>>1;
	LL v;
	for(int m=2;m<=num;half=m,m<<=1,lim>>=1)
	{
		for(register int i=0;i<half;i++)
		{
			LL w=(!pd)?wi[i*lim]:wi[num-i*lim];
			for(register int j=i;j<num;j+=m)
			{
				v=a[j+half]*w%mo;
				a[j+half]=(a[j]-v+mo)%mo;
				a[j]=(a[j]+v)%mo;
			}
		}
	}	
	if(pd) fo(i,0,num-1) a[i]=a[i]*ny2%mo;
}
void make(int l,LL *a,LL *b)
{
	b[0]=ksm(a[0],mo-2);
	for(int m=1,t=2,num=4,cnt=2;m<l;m=t,t=num,num<<=1,cnt++)
	{
		prp(num,cnt);
		fo(i,0,m-1) c[i]=a[i],d[i]=b[i];
		fo(i,m,t-1) c[i]=a[i];
		fo(i,t,num-1) c[i]=0;
		NTT(c,0,num),NTT(b,0,num);
		fo(i,0,num-1) b[i]=b[i]*b[i]%mo*c[i]%mo;
		NTT(b,1,num);
		fo(i,0,t-1) b[i]=((LL)2*d[i]-b[i]+mo)%mo;
		fo(i,t,num-1) b[i]=0;
	}		
}
int main()
{
	cin>>n;
	fo(i,0,n-1) scanf("%lld",&a[i]);
	wg[0]=1;
	LL c=ksm(3,(mo-1)/M);
	fo(i,1,M) wg[i]=wg[i-1]*c%mo;
	make(n,a,b);
	fo(i,0,n-1) printf("%lld ",b[i]);
}

多项式除法（多项式取模）

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstring>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define M 524288
#define mo 998244353
#define LL long long
#define N 100005
using namespace std;
LL a[M+1],b[M+1],c[M+1],wi[M+1],wg[M+1],ny2,d[M+1],r[M+1],q[M+1],l2[M+1],cf[20],f[M+1];
int bit[M+1],n,m;
LL ksm(LL k,LL n)
{
	LL s=1;
	for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;
	return s;
}
void prp(int num)
{
	ny2=ksm(num,mo-2);
	fo(i,0,num) wi[i]=wg[i*(M/num)];
	fo(i,0,num-1) bit[i]=(bit[i>>1]>>1)|((i&1)<<(l2[num]-1));
}
void NTT(LL *a,bool pd,int num)
{
	fo(i,0,num-1) if(bit[i]<i) swap(a[bit[i]],a[i]);
	int half=1,lim=num>>1;
	LL v;
	for(int m=2;m<=num;half=m,m<<=1,lim>>=1)
	{
		for(register int i=0;i<half;i++)
		{
			LL w=(!pd)?wi[i*lim]:wi[num-i*lim];
			for(register int j=i;j<num;j+=m)
			{
				v=a[j+half]*w%mo;
				a[j+half]=(a[j]-v+mo)%mo;
				a[j]=(a[j]+v)%mo;
			}
		}
	}	
	if(pd) fo(i,0,num-1) a[i]=a[i]*ny2%mo;
}
void make(int l,LL *a,LL *b)
{
	b[0]=ksm(a[0],mo-2);
	for(int m=1,t=2,num=4;m<l;m=t,t=num,num<<=1)
	{
		prp(num);
		fo(i,0,m-1) c[i]=a[i],d[i]=b[i];
		fo(i,m,t-1) c[i]=a[i];
		fo(i,t,num-1) c[i]=0;
		NTT(c,0,num),NTT(b,0,num);
		fo(i,0,num-1) b[i]=b[i]*b[i]%mo*c[i]%mo;
		NTT(b,1,num);
		fo(i,0,t-1) b[i]=((LL)2*d[i]-b[i]+mo)%mo;
		fo(i,t,num-1) b[i]=0;
	}		
}
void rev(int num,LL *a,LL *b)
{
	fo(i,0,num-1) b[i]=a[num-i-1];
}
void div(LL *a,LL *b,LL *d,LL *r)
{
	rev(m+1,b,r);
	make(n-m+1,r,d);
	int num=cf[l2[n+m+1]];
	prp(num);
	fo(i,0,n) f[i]=a[n-i];
	fo(i,n-m+1,num-1) f[i]=0,d[i]=0;
	NTT(f,0,num),NTT(d,0,num);
	fo(i,0,num-1) d[i]=d[i]*f[i]%mo;
	NTT(d,1,num);
	fo(i,n-m+1,num-1) d[i]=0;
	fo(i,0,(n-m)>>1) swap(d[i],d[n-m-i]);
	fo(i,0,num-1) r[i]=d[i],f[i]=b[i];
	num=cf[l2[n+1]];
	prp(num);
	NTT(f,0,num),NTT(r,0,num);
	fo(i,0,num-1) r[i]=r[i]*f[i]%mo;
	NTT(r,1,num);
	fo(i,0,num-1) r[i]=(a[i]-r[i]+mo)%mo;
}
int main()
{
	cf[0]=1;
    fo(i,1,19) cf[i]=cf[i-1]*2,l2[cf[i]]=i;
	cin>>n>>m;
	fo(i,0,n) scanf("%lld",&a[i]);
	fo(i,0,m) scanf("%lld",&b[i]);
	wg[0]=1;
	LL c=ksm(3,(mo-1)/M);
	fod(i,M-1,1) if(!l2[i]) l2[i]=l2[i+1];
	fo(i,1,M) wg[i]=wg[i-1]*c%mo;
	div(a,b,q,r);
	fo(i,0,n-m) printf("%lld ",q[i]);
	printf("\n");
	fo(i,0,m-1) printf("%lld ",r[i]);
}

多点求值（洛谷 P5050）

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define M 262144
#define mo 998244353
#define N 100005
#define LL long long
using namespace std;
LL a[M+1],b[M+1],c[M+1],u1[M+1],wg[M+1],u2[M+1],wi[M+1],r[M+1],q[M+1],cf[19],l2[M+1],ny[M+1],a1[20][M+1],ans[N],s1[20][M+1];
int bit[M+1],n,m,fi[20][N+1],le[20][N+1],lf;
LL ksm(LL k,LL n)
{
	LL s=1;
	for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;
	return s;
}
void prp(int num)
{
	fo(i,0,num) wi[i]=wg[i*(M/num)];
	fo(i,0,num-1) bit[i]=(bit[i>>1]>>1)|((i&1)<<(l2[num]-1));
}
void NTT(LL *a,bool pd,int num)
{
	fo(i,0,num-1) if(i<bit[i]) swap(a[i],a[bit[i]]);
	LL w,v;
	for(int lim=num>>1,m=2,half=1;m<=num;half=m,lim>>=1,m<<=1)
	{
		fo(i,0,half-1)
		{
			w=(!pd)?wi[i*lim]:wi[num-i*lim];
			for(int j=i;j<num;j+=m)
			{
				v=w*a[j+half]%mo;
				a[j+half]=(a[j]-v+mo)%mo;
				a[j]=(a[j]+v)%mo;
			}
		}
	}
	if(pd) fo(i,0,num-1) a[i]=a[i]*ny[num]%mo;
}
void make(int l,LL *a,LL *b)
{
	b[0]=ksm(a[0],mo-2);
	for(int m=1,t=2,num=4;m<=l;m=t,t=num,num<<=1)
	{
		prp(num);
		fo(i,0,m-1) c[i]=a[i],u2[i]=b[i];
		fo(i,m,t-1) c[i]=a[i],u2[i]=b[i]=0;
		fo(i,t,num-1) c[i]=u2[i]=b[i]=0;
		NTT(c,0,num),NTT(b,0,num);
		fo(i,0,num-1) b[i]=b[i]*b[i]%mo*c[i]%mo;
		NTT(b,1,num);
		fo(i,0,t-1) b[i]=((LL)2*u2[i]-b[i]+mo)%mo;
		fo(i,t,num-1) b[i]=0;
	}
}
void rev(int num,LL *a,LL *b)
{
	fo(i,0,num-1) b[i]=a[num-i-1];
}
void div(int n,int m,LL *a,LL *b,LL *d,LL *r)
{
	rev(m,b,r);
	fo(i,m,n-m+1) r[i]=0;
	make(n-m+1,r,d);
	int num=cf[l2[n]]*2;
	prp(num),rev(n,a,u2);
	fo(i,n-m+1,num-1) d[i]=u2[i]=0;
	
	NTT(d,0,num),NTT(u2,0,num);
	fo(i,0,num-1) d[i]=d[i]*u2[i]%mo;
	NTT(d,1,num);
	fo(i,n-m+1,num-1) d[i]=0;
	fo(i,0,(n-m)>>1) swap(d[i],d[n-m-i]);

	num=cf[l2[n]],prp(num);
	fo(i,0,m-1) r[i]=d[i],u2[i]=b[i];
	fo(i,m,num-1) r[i]=d[i],u2[i]=0;
	NTT(r,0,num),NTT(u2,0,num);
	fo(i,0,num-1) r[i]=r[i]*u2[i]%mo;
	NTT(r,1,num);
	fo(i,0,m-1) r[i]=(a[i]-r[i]+mo)%mo;
	fo(i,m,num-1) r[i]=0;
}
void prd(int t,int l,int r)
{
	if(l==r)
	{
		fi[t][l]=lf+1;
		le[t][l]=2;
		a1[t][fi[t][l]]=-b[l],a1[t][fi[t][l]+1]=1;
		lf+=2;
		return;
	}
	int mid=(l+r)>>1;
	prd(t+1,l,mid),prd(t+1,mid+1,r);
}
void doit(int t,int l,int r)
{
	if(l==r) return;
	int mid=(l+r)>>1;
	doit(t+1,l,mid),doit(t+1,mid+1,r);
	fi[t][l]=fi[t+1][l];
	int num=cf[l2[le[t+1][l]+le[t+1][mid+1]-1]];
	prp(num);
	fo(i,0,num-1)
	{
		u1[i]=(i<le[t+1][l])?a1[t+1][fi[t+1][l]+i]:0;
		u2[i]=(i<le[t+1][mid+1])?a1[t+1][fi[t+1][mid+1]+i]:0;
	}
	NTT(u1,0,num),NTT(u2,0,num);
	fo(i,0,num-1) u1[i]=u1[i]*u2[i]%mo;
	NTT(u1,1,num);
	le[t][l]=le[t+1][l]+le[t+1][mid+1]-1;
	fo(i,0,le[t][l]) a1[t][fi[t][l]+i]=u1[i];
}
void query(int t,int l,int r,int n)
{
	if(l==r) ans[l]=s1[t][0];
	else
	{
		int mid=(l+r)>>1;
		div(n,le[t+1][l],s1[t],a1[t+1]+fi[t+1][l],q,s1[t+1]);
		query(t+1,l,mid,le[t+1][l]-1);
		div(n,le[t+1][mid+1],s1[t],a1[t+1]+fi[t+1][mid+1],q,s1[t+1]);
		query(t+1,mid+1,r,le[t+1][mid+1]-1);
	}
}
int main()
{
	cf[0]=1;
	fo(i,1,18) cf[i]=cf[i-1]<<1,l2[cf[i]]=i;
	fod(i,M-1,1) if(!l2[i]) l2[i]=l2[i+1];
	ny[1]=1;
	fo(i,2,M) ny[i]=(-ny[mo%i]*(LL)(mo/i)%mo+mo)%mo;
	wg[0]=1;
	LL c=ksm(3,(mo-1)/M);
	fo(i,1,M) wg[i]=wg[i-1]*c%mo;
	cin>>n>>m;
	fo(i,0,n) scanf("%lld",&a[i]);
	fo(i,0,m-1) scanf("%lld",&b[i]);
	prd(0,0,m-1);
	doit(0,0,m-1);
	fo(i,0,n) s1[0][i]=a[i];
	query(0,0,m-1,n+1);
	fo(i,0,m-1) printf("%lld\n",ans[i]);
}

多项式开根（洛谷P5205）

#include <bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define M 524288
#define LL long long
#define mo 998244353
using namespace std;
int l2[M+1],cf[22],n,t,bit[M+1];
LL wi[M+1],wg[M+1],b[M+1],u1[M+1],u2[M+1],u3[M+1],ny[M+1],c[M+1];
void prp(int num)
{
	fo(i,0,num) bit[i]=(bit[i>>1]>>1)|((i&1)<<(l2[num]-1)),wi[i]=wg[i*(M/num)];
}
LL ksm(LL k,LL n)
{
	LL s=1;
	for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;
	return s;
}
void NTT(LL *a,bool pd,int num)
{
	fo(i,0,num-1) if(i<bit[i]) swap(a[i],a[bit[i]]);
	LL v;	
	for(register int lim=num>>1,h=1,m=2;m<=num;h=m,lim>>=1,m<<=1)
	{
		for(register int j=0;j<num;j+=m)
		{
			fo(i,0,h-1)
			{
				v=((!pd)?wi[i*lim]:wi[num-i*lim])*a[i+j+h]%mo;	
				a[i+j+h]=(a[i+j]-v+mo)%mo;
				a[i+j]=(a[i+j]+v)%mo;
			}
		}
	}
	if(pd) fo(i,0,num-1) a[i]=a[i]*ny[num]%mo;
}
void inv(int l,LL *a,LL *b)
{
	b[0]=ksm(a[0],mo-2);
	for(int m=1,t=2,num=4;m<l;m=t,t=num,num<<=1)
	{
		prp(num);
		fo(i,0,m-1) u1[i]=a[i],u2[i]=b[i];
		fo(i,m,t-1) u1[i]=a[i],u2[i]=b[i]=0;
		fo(i,t,num-1) u1[i]=0,u2[i]=b[i]=0;
		NTT(u1,0,num),NTT(b,0,num);
		fo(i,0,num-1) b[i]=b[i]*b[i]%mo*u1[i]%mo,u1[i]=0;
		NTT(b,1,num);
		fo(i,0,t-1) b[i]=((LL)2*u2[i]%mo-b[i]+mo)%mo,u2[i]=0;
		fo(i,t,num-1) b[i]=0;
	}		
}
void make(int l,LL *a,LL *b)
{
	b[0]=1;
	for(int m=1,t=2,num=4;m<l;m=t,t=num,num<<=1)
	{
		inv(t,b,u3);
		fo(i,0,t-1) u1[i]=a[i];
		fo(i,t,num-1) u1[i]=u3[i]=0;
		prp(num);
		NTT(u1,0,num),NTT(u3,0,num);
		fo(i,0,num-1) u1[i]=u1[i]*u3[i]%mo,u3[i]=0;
		NTT(u1,1,num);
		fo(i,0,t-1) b[i]=(u1[i]+b[i]+mo)%mo*ny[2]%mo,u1[i]=0;
		fo(i,t,num-1) b[i]=0,u1[i]=0;
	}
}
int main()
{
	cf[0]=1;
	fo(i,1,19) l2[cf[i]=cf[i-1]<<1]=i;
	fod(i,M-1,2) if(!l2[i]) l2[i]=l2[i+1];
	wg[0]=1,wg[1]=ksm(3,(mo-1)/M);
	fo(i,2,M) wg[i]=wg[i-1]*wg[1]%mo;
	ny[1]=1;
	fo(i,2,M) ny[i]=(-ny[mo%i]*(LL)(mo/i)%mo+mo)%mo;
	
	int n;
	cin>>n;
	fo(i,0,n-1) scanf("%lld",&b[i]);
	make(cf[l2[n]],b,c);
	fo(i,0,n-1) printf("%lld ",c[i]);
}

多项式取ln（洛谷P4725）

#include <bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define mo 998244353
#define LL long long
#define M 262144
#define N 100005
#define L 18
using namespace std;
LL wi[M+1],wg[M+1],a[M+1],b[M+1],ny[M+1];
int bit[M+1],n,cf[L+1],l2[M+1];
LL ksm(LL k,LL n)
{
	LL s=1;
	for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;
	return s;
}
void prp(int num)
{
	fo(i,0,num)
	{
		bit[i]=(bit[i>>1]>>1)|((i&1)<<(l2[num]-1));
		wi[i]=wg[i*(M/num)];
	}
	ny[num]=ksm(num,mo-2);
}
void NTT(LL *a,bool pd,int num)
{
	fo(i,0,num-1) if(bit[i]<i) swap(a[i],a[bit[i]]);
	LL v;
	for(int h=1,m=2,lim=num>>1;m<=num;h=m,m<<=1,lim>>=1)
	{
		for(int j=0;j<num;j+=m)
		{
			fo(i,0,h-1)
			{
				v=((!pd)?wi[i*lim]:wi[num-i*lim])*a[i+j+h]%mo;
				a[i+j+h]=(a[i+j]-v+mo)%mo;
				a[i+j]=(a[i+j]+v)%mo;
			}
		}
	}
	if(pd) fo(i,0,num-1) a[i]=a[i]*ny[num]%mo;
}
void inv(int n,LL *a,LL *b)
{
	static LL u1[M+1],u2[M+1];
	fo(i,0,n-1) b[i]=0;
	b[0]=ksm(a[0],mo-2);
	for(int m=1,t=2,num=4;m<n;m=t,t=num,num<<=1)
	{
		prp(num);
		fo(i,0,num-1) u1[i]=u2[i]=0;
		fo(i,0,m-1) u1[i]=b[i];
		fo(i,0,t-1) u2[i]=a[i];
		NTT(u1,0,num),NTT(u2,0,num);
		fo(i,0,num-1) u1[i]=u1[i]*u1[i]%mo*u2[i]%mo;
		NTT(u1,1,num);
		fo(i,0,t-1) b[i]=((LL)2*b[i]-u1[i]+mo)%mo;
	}
}
void getln(int n,LL *a,LL *b)
{
	static LL u3[M+1],u4[M+1];
	int num=cf[l2[2*n+1]];
	inv(n+1,a,u3);
	fo(i,n+1,num) u3[i]=0;
	prp(num);
	fo(i,0,num-1) u4[i]=0;
	fo(i,0,n-1) u4[i]=a[i+1]*(LL)(i+1)%mo;
	NTT(u4,0,num),NTT(u3,0,num);
	fo(i,0,num-1) u4[i]=u3[i]*u4[i]%mo;
	NTT(u4,1,num);
	b[0]=0;
	fo(i,1,n-1) b[i]=u4[i-1]*ny[i]%mo;
}
int main()
{
	cin>>n;
	cf[0]=1;
	fo(i,1,L) l2[cf[i]=cf[i-1]<<1]=i;
	fod(i,M-1,2) if(!l2[i]) l2[i]=l2[i+1];
	ny[1]=1;
	fo(i,2,M) ny[i]=(-ny[mo%i]*(LL)(mo/i)%mo+mo)%mo; 
	wg[0]=1,wg[1]=ksm(3,(mo-1)/M);
	fo(i,2,M) wg[i]=wg[i-1]*wg[1]%mo;
	fo(i,0,n-1) scanf("%lld",&a[i]);
	getln(n,a,b);
	fo(i,0,n-1) printf("%lld ",b[i]);
}

推荐阅读