多项式全家桶

已更新乘法、求逆、开平方、ln&exp、快速幂、带余除法。常数平均在洛谷前两页。

#include<bits/stdc++.h>
#define RG register
#define R RG int
#define ll long long
#define I inline
using namespace std;
const int bufsz=1<<20; 
char buf[bufsz],*p1,*p2;
#define GC (p1==p2&&(p2=(p1=buf)+fread(buf,1,bufsz,stdin),p1==p2)?0:*p1++)
I int read()
{
	char ch=GC;int x=0;
	while(ch<'0'||ch>'9')ch=GC;
	while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=GC;}
	return x;
}
const ll P=998244353,G=3,iG=332748118;
const int _N=1<<17,N=_N<<2;
I int pls(int x,int y){return x+y>=P?x+y-P:x+y;}
I int mns(int x,int y){return x-y<0?x-y+P:x-y;}
I int div2(int x)
{
	if(x&1)return x+P>>1;
	return x>>1;
}
int update(int x)
{
	int l;
	for(l=1;l<x;l<<=1);
	return l;
}
ll qpow(ll a,ll b=P-2)
{
	ll res=1;
	for(;b;b>>=1,a=a*a%P)
		if(b&1)res=res*a%P;
	return res;
}
namespace Cipolla
{
	ll U;
	struct Complex
	{
		ll a,b;
		Complex(ll r,ll i):a(r),b(i){}
		I void operator*=(Complex &x)
		{
			ll k=(a*x.a+U*b%P*x.b)%P;
			b=(a*x.b+b*x.a)%P;a=k;
		}
		friend Complex cpow(Complex a,ll b)
		{
			Complex res(1,0);
			for(;b;b>>=1,a*=a)
				if(b&1)res*=a;
			return res;
		}
	};
	ll sqrt(ll n)
	{	
		if(n<=1)return n;
		ll a;
		do a=rand()%P;while(qpow(U=(a*a-n)%P,P-1>>1)==1);
		ll ans=cpow(Complex(a,1),P+1>>1).a;
		if(ans<0)ans+=P;
		return min(ans,P-ans);
	}
}
namespace Poly
{
	int iv[N]={0,1},mx=1,w[N]={1},r[N],nr=0,tx[N],ty[N],tz[N],to[N];
	void inv_init(int n)
	{
		if(mx>=n)return;
		for(R i=mx+1;i<=n;i++)iv[i]=(P-P/i)*iv[P%i]%P;
		mx=n;
	}
	void init(int n)
	{
		if(nr==n)return;nr=n;
		for(R i=1;i<n;i++)
		{
			r[i]=r[i>>1]>>1;
			if(i&1)r[i]|=n>>1;
		}
	}
	void NTT(int *s,int n,bool f)
	{
		for(R i=0;i<n;i++)if(i<r[i])swap(s[i],s[r[i]]);
		for(R i=1;i<n;i<<=1)
		{
			ll ng=qpow(f?iG:G,(P-1)/(i<<1));
			for(R j=i-2>>1;~j;j--)w[j<<1|1]=(w[j<<1]=w[j])*ng%P;
			for(R j=0;j<n;j+=i<<1)
				for(R k=j;k<j+i;k++)
				{
					R x=s[k],y=1ll*w[k-j]*s[k+i]%P;
					s[k]=pls(x,y);
					s[k+i]=mns(x,y);
				}
		}
		if(f)
		{
			ll v=qpow(n);
			for(R i=0;i<n;i++)s[i]=s[i]*v%P;
		}
	}
	void mul(int *a,int *b,int *c,int p,int n=N,int m=N)
	{
		n=min(n,p);m=min(m,p);
		int o=n+m-1,l=update(o);
		for(R i=0;i<l;i++)
		{
			tx[i]=i<n?a[i]:0;
			ty[i]=i<m?b[i]:0;
		}
		init(l);
		NTT(tx,l,0);NTT(ty,l,0);
		for(R i=0;i<l;i++)c[i]=1ll*tx[i]*ty[i]%P;
		NTT(c,l,1);
		for(R i=l;i<p;i++)c[i]=0;
	}
	void inv(int *a,int *b,int p,int n)
	{
		R l;
		b[0]=qpow(a[0]);b[1]=0;
		for(l=1;l<p;l<<=1)
		{
			for(R i=0;i<l<<1;i++)tx[i]=i<n?a[i]:0;
			for(R i=l<<1;i<l<<2;i++)b[i]=tx[i]=0;
			init(l<<2);
			NTT(tx,l<<2,0);NTT(b,l<<2,0);
			for(R i=0;i<l<<2;i++)b[i]=(2+P-1ll*b[i]*tx[i]%P)*b[i]%P;
			NTT(b,l<<2,1);
			for(R i=l<<1;i<l<<2;i++)b[i]=0;
		}
	}
	int div(int *a,int *b,int *c,int *d,int n,int m)
	{
		if(n<m)
		{
			for(R i=0;i<n;i++)d[i]=a[i];
			return n;
		}
		reverse_copy(a,a+n,tz);
		reverse_copy(b,b+m,to);
		int p=n-m+1;
		inv(to,ty,p,m);
		mul(tz,ty,c,p,n,p);
		reverse(c,c+p);
		mul(b,c,d,m-1,m,p);
		for(R i=0;i<m-1;i++)d[i]=mns(a[i],d[i]);
		return m-1;
	}
	void sqrt(int *a,int *b,int p,int n)
	{
		R l;
		b[0]=Cipolla::sqrt(a[0]);
		for(l=1;l<p;l<<=1)
		{
			inv(b,ty,l<<1,l);
			for(R i=0;i<l<<1;i++)tx[i]=i<n?a[i]:0;
			for(R i=l<<1;i<l<<2;i++)ty[i]=tx[i]=0;
			init(l<<2);
			NTT(tx,l<<2,0);NTT(ty,l<<2,0);
			for(R i=0;i<l<<2;i++)tx[i]=1ll*tx[i]*ty[i]%P;
			NTT(tx,l<<2,1);
			for(R i=0;i<l;i++)b[i]=div2(pls(b[i],tx[i]));
			for(R i=l;i<l<<1;i++)b[i]=div2(tx[i]);
		}
	}
	void derv(int *a,int *b,int n)
	{
		for(R i=1;i<n;i++)b[i-1]=1ll*a[i]*i%P;
		b[n-1]=0;
	}
	void intg(int *a,int *b,int n)
	{
		inv_init(n-1);
		for(R i=n-1;i;i--)b[i]=1ll*a[i-1]*iv[i]%P;
		b[0]=0;
	}
	void ln(int *a,int *b,int p,int n=N)
	{
		n=min(n,p);
		inv(a,b,p,n);derv(a,tx,n);
		R l=update(p);
		for(R i=n;i<l<<1;i++)tx[i]=0;
		for(R i=p;i<l<<1;i++)b[i]=0;
		init(l<<1);
		NTT(tx,l<<1,0);NTT(b,l<<1,0);
		for(R i=0;i<l<<1;i++)b[i]=1ll*tx[i]*b[i]%P;
		NTT(b,l<<1,1);
		intg(b,b,p);
	}
	void exp(int *a,int *b,int p,int n=N)
	{
		n=min(n,p);
		R l;
		b[0]=1;b[1]=0;
		for(l=1;l<p;l<<=1)
		{
			ln(b,ty,l<<1,l);
			for(R i=0;i<l<<1;i++)tx[i]=mns(i<n?a[i]:0,ty[i]);
			for(R i=l<<1;i<l<<2;i++)b[i]=tx[i]=0;
			tx[0]++;
			init(l<<2);
			NTT(b,l<<2,0);NTT(tx,l<<2,0);
			for(R i=0;i<l<<2;i++)b[i]=1ll*b[i]*tx[i]%P;
			NTT(b,l<<2,1);
			for(R i=l<<1;i<l<<2;i++)b[i]=0;
		}
	}
	void power(int *a,int *b,int c1,int c2,int p,int n=N)
	{
		n=min(n,p);R t,ut;
		for(t=0;!a[t]&&t<n;t++);
		if(t==n||1ll*t*c1>=p)return;
		ut=t*c1;
		ll u=qpow(a[t],c2),v=qpow(a[t]);
		for(R i=0;i+t<n;i++)a[i]=1ll*a[i+t]*v%P;
		ln(a,tz,p-ut,n-t);
		for(R i=0;i<p-ut;i++)tz[i]=1ll*tz[i]*c1%P;
		exp(tz,b,p-ut,p-ut);
		for(R i=p-1;i-ut>=0;i--)b[i]=1ll*b[i-ut]*u%P;
		for(R i=ut-1;~i;i--)b[i]=0;
	}
}
int a[N],b[N],c[N],d[N];
int main()
{
	int n=read()+1,m=read()+1;
	for(R i=0;i<n;i++)a[i]=read();
	for(R i=0;i<m;i++)b[i]=read();
	int len=Poly::div(a,b,c,d,n,m);
	for(R i=0;i<n-m+1;i++)printf("%d ",c[i]);puts("");
	for(R i=0;i<len;i++)printf("%d ",d[i]);
	return 0;
}
  • FFT+MTT(4FFT版)+任意模数求逆(跑了洛谷第一!!)
#include<bits/stdc++.h>
#define RG register
#define R RG int
#define ll long long
#define D double
#define I inline
using namespace std;
const int bufsz=1<<20; 
char buf[bufsz],*p1,*p2;
#define GC (p1==p2&&(p2=(p1=buf)+fread(buf,1,bufsz,stdin),p1==p2)?0:*p1++)
I int read()
{
	char ch=GC;int x=0;
	while(ch<'0'||ch>'9')ch=GC;
	while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=GC;}
	return x;
}
const int _N=1e5+3,N=1<<18,T=15,TT=T<<1;
const ll S=1<<T,P=1e9+7;
const D pi=acos(-1);
ll qpow(ll a,ll b=P-2)
{
	int res=1;
	for(;b;b>>=1,a=a*a%P)
		if(b&1)res=res*a%P;
	return res;
}
struct Cp
{
	D re,im;
	Cp(D r=0,D i=0):re(r),im(i){};
	I Cp operator+(const Cp &a){return Cp(re+a.re,im+a.im);}
	I Cp operator-(const Cp &a){return Cp(re-a.re,im-a.im);}
	I Cp operator*(const Cp &a){return Cp(re*a.re-im*a.im,re*a.im+im*a.re);}
};
Cp w[N];
int r[N];
void init_r(int n)
{
	for(R i=1;i<n;i++)
	{
		r[i]=r[i>>1]>>1;
		if(i&1)r[i]|=n>>1;
	}
}
void FFT(Cp *s,int n,bool f)
{
	for(R i=0;i<n;i++)if(i<r[i])swap(s[i],s[r[i]]);
	for(R i=1;i<n;i<<=1)
	{
		const Cp wn(cos(pi/i),f?-sin(pi/i):sin(pi/i));
		for(R j=i-2>>1;~j;j--)w[j<<1|1]=(w[j<<1]=w[j])*wn;
		for(R j=0;j<n;j+=i<<1)
			for(R k=j;k<j+i;k++)
			{
				Cp x=s[k],y=s[k+i]*w[k-j];
				s[k]=x+y;s[k+i]=x-y;
			}
	}
	if(f)
	{
		const D in=1.0/n;
		for(R i=0;i<n;i++)s[i].re*=in,s[i].im*=in;
	}
}
void MTT(int *a,int *b,int *res,int n,int m,int l,int x,int y)
{
	static Cp A[N],B[N],c0[N],c1[N];
	for(R i=0;i<n;i++)A[i]=Cp(a[i]&S-1,a[i]>>T);for(R i=n;i<l;i++)A[i]=Cp();
	for(R i=0;i<m;i++)B[i]=Cp(b[i]&S-1,b[i]>>T);for(R i=m;i<l;i++)B[i]=Cp();
	FFT(A,l,0);FFT(B,l,0);
	c0[0]=Cp(A[0].re)*B[0];
	c1[0]=Cp(A[0].im)*B[0];
	for(R i=1,j=l-1;i<l;i++,j--)
	{
		c0[i]=Cp((A[i].re+A[j].re)*0.5,(A[i].im-A[j].im)*0.5)*B[i];
		c1[i]=Cp((A[i].im+A[j].im)*0.5,(-A[i].re+A[j].re)*0.5)*B[i];
	}
	FFT(c0,l,1);FFT(c1,l,1);
	for(R i=x;i<y;i++)
	{
		ll c00=c0[i].re+0.5,c01=c0[i].im+0.5,c10=c1[i].re+0.5,c11=c1[i].im+0.5;
		res[i-x]=(c00+((c01+c10)%P<<T)+(c11%P<<TT))%P;
	}
}
void inv(int *a,int *res,int n)
{
	static int tp[N];
	res[0]=qpow(a[0]);
	for(R l=1;l<n;l<<=1)
	{
		init_r(l<<1);
		MTT(a,res,tp,l<<1,l,l<<1,l,l<<1);
		for(R i=0;i<l;i++)tp[i]=P-tp[i];
		MTT(res,tp,res+l,l,l,l<<1,0,l);
	}
}
int a[_N],b[_N];
int main()
{
	w[0].re=1;
	int n;
	scanf("%d",&n);
	for(R i=0;i<n;i++)a[i]=read();
	inv(a,b,n);
	for(R i=0;i<n;i++)printf("%d ",b[i]); 
	return 0;
}
posted @ 2021-01-21 01:01  nkxjlym  阅读(54)  评论(0)    收藏  举报