多项式全家桶
已更新乘法、求逆、开平方、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;
}
浙公网安备 33010602011771号