记录一些多项式代码
FFT
const double PI=acos(-1.0);
struct Complex
{
double R,I;
Complex(){}
Complex(double _R,double _I)
{
R=_R;
I=_I;
}
};
Complex operator + (Complex x,Complex y)
{
return Complex(x.R+y.R,x.I+y.I);
}
Complex operator - (Complex x,Complex y)
{
return Complex(x.R-y.R,x.I-y.I);
}
Complex operator * (Complex x,Complex y)
{
return Complex(x.R*y.R-x.I*y.I,x.R*y.I+x.I*y.R);
}
void FFT(Complex *a,int n,double type)
{
for(int i=0,j=0;i<n;i++)
{
if(i<j)
{
swap(a[i],a[j]);
}
int k=n>>1;
while(true)
{
j^=k;
if(j>=k)
{
break;
}
else
{
k>>=1;
}
}
}
for(int m=1;m<n;m<<=1)
{
Complex w=Complex(cos(PI/m),sin(PI/m)*type);
for(int k=0;k<n;k+=(m<<1))
{
Complex wi=Complex(1.0,0.0);
for(int i=k;i<k+m;i++,wi=wi*w)
{
Complex t0=a[i],t1=a[i+m]*wi;
a[i]=t0+t1;
a[i+m]=t0-t1;
}
}
}
if(type==1)
return;
for(int i=0;i<n;i++)
{
a[i].R/=n;
a[i].I/=n;
}
}
Complex AA[MAXN],BB[MAXN];
int N;
void FFT(int n,int *a,int m,int *b,int *ans)
{
N=1;
for(int i=0;i<=n;i++)
{
AA[i].R=a[i];
}
for(int i=0;i<=m;i++)
{
BB[i].R=b[i];
}
while(N<=n+m)
{
N<<=1;
}
FFT(AA,N,1);
FFT(BB,N,1);
for(int i=0;i<=N;i++)
{
AA[i]=AA[i]*BB[i];
}
FFT(AA,N,-1);
for(int i=0;i<=N;i++)
{
ans[i]=int(AA[i].R+0.5);
AA[i].R=AA[i].I=BB[i].R=BB[i].I=0;
}
}
\(n\) 是 \(a\) 多项式最高次大小，\(m\) 是 \(b\) 多项式最高次大小，\(ans\) 为答案多项式。
NTT
const int P=998244353,_3=332748118;
long long Pow(long long a,long long b)
{
int ans=1;
while(b)
{
if(b&1)
ans=ans*a%P;
a=a*a%P;
b>>=1;
}
return ans;
}
void NTT(long long *a,int n,int type)
{
for(int i=0,j=0;i<n;i++)
{
if(i<j)
{
swap(a[i],a[j]);
}
int k=n>>1;
while(true)
{
j^=k;
if(j>=k)
{
break;
}
else
{
k>>=1;
}
}
}
for(int m=1;m<n;m<<=1)
{
long long w=Pow(type==1?3:_3,(P-1)/(m<<1));
for(int k=0;k<n;k+=(m<<1))
{
long long wi=1;
for(int i=k;i<k+m;i++,wi=wi*w%P)
{
long long t0=a[i],t1=a[i+m]*wi%P;
a[i]=(t0+t1)%P;
a[i+m]=(t0-t1+P)%P;
}
}
}
if(type==1)
return;
int _n=Pow(n,P-2);
for(int i=0;i<n;i++)
{
a[i]=a[i]*_n%P;
}
}
long long AA[MAXN],BB[MAXN];
int N;
void NTT(int n,int *a,int m,int *b,int *ans)
{
memset(AA,0,sizeof(AA));
memset(BB,0,sizeof(BB));
N=1;
for(int i=0;i<=n;i++)
{
AA[i]=a[i];
}
for(int i=0;i<=m;i++)
{
BB[i]=b[i];
}
while(N<=n+m)
{
N<<=1;
}
NTT(AA,N,1);
NTT(BB,N,1);
for(int i=0;i<=N;i++)
{
AA[i]=AA[i]*BB[i]%P;
}
NTT(AA,N,-1);
for(int i=0;i<=N;i++)
{
ans[i]=AA[i];
}
}
_3 是 \(3\) 关于 \(P\) 的逆元，\(P\) 可以更改。用法同 FFT。
比较拉的全家桶。
#include <bits/stdc++.h>
using namespace std;
const int Gen=3,P=998244353,MAXN=1e6+50;
void Add1(int &a,int b)
{
a+=b;
if(a>=P)
a-=P;
}
int Add2(int a,int b)
{
static int Sum;
Sum=a+b;
if(Sum>=P)
return Sum-P;
else
return Sum;
}
void Sub1(int &a,int b)
{
a-=b;
if(a<0)
a+=P;
}
int Sub2(int a,int b)
{
static int Sum;
Sum=a-b;
if(Sum<0)
return Sum+P;
else
return Sum;
}
void Mul1(int &a,int b)
{
a=1ll*a*b%P;
}
int Mul2(int a,int b)
{
return 1ll*a*b%P;
}
void CopyArray(int *a,int *b,int l,int r)
{
for(int i=l;i<=r;i++)
{
a[i]=b[i];
}
}
void ClearArray(int *a,int l,int r)
{
for(int i=l;i<=r;i++)
{
a[i]=0;
}
}
void AddArray(int *a,int *b,int l,int r)
{
for(int i=l;i<=r;i++)
{
Add1(a[i],b[i]);
}
}
void MulArray(int *a,int *b,int l,int r)
{
for(int i=l;i<=r;i++)
{
Mul1(a[i],b[i]);
}
}
void RevArray(int *a,int l,int r)
{
static int b[MAXN];
for(int i=l;i<=r;i++)
{
b[i]=a[i];
}
for(int i=l,j=r;i<=r;i++,j--)
{
a[i]=b[j];
b[j]=0;
}
}
int Pow(int a,int b=P-2)
{
int ans=1;
while(b)
{
if(b&1)
{
Mul1(ans,a);
}
Mul1(a,a);
b>>=1;
}
return ans;
}
const int InvGen=Pow(Gen);
void NTT(int *a,int n,int type)
{
for(int i=0,j=0;i<n;i++)
{
if(i<j)
{
swap(a[i],a[j]);
}
int k=n>>1;
while(true)
{
j^=k;
if(j>=k)
{
break;
}
else
{
k>>=1;
}
}
}
for(int m=1;m<n;m<<=1)
{
int w=Pow(type==1?3:InvGen,(P-1)/(m<<1));
for(int k=0;k<n;k+=(m<<1))
{
int wi=1;
for(int i=k;i<k+m;i++,Mul1(wi,w))
{
int t0=a[i],t1=Mul2(a[i+m],wi);
a[i]=Add2(t0,t1);
a[i+m]=Sub2(t0,t1);
}
}
}
if(type==1)
return;
int _n=Pow(n);
for(int i=0;i<n;i++)
{
Mul1(a[i],_n);
}
}
void MulPoly(int lim,int *a,int *b,int Mod=1e9)
{
static int BB[MAXN];
NTT(a,lim,1);
CopyArray(BB,b,0,lim);
NTT(BB,lim,1);
MulArray(a,BB,0,lim);
NTT(a,lim,-1);
ClearArray(a,Mod+1,lim);
ClearArray(BB,0,lim);
}
void MulPoly(int n,int *a,int m,int *b,int Mod=1e9)
{
static int BB[MAXN];
int lim=1;
while(lim<=n+m)
{
lim<<=1;
}
NTT(a,lim,1);
CopyArray(BB,b,0,m);
NTT(BB,lim,1);
MulArray(a,BB,0,lim);
NTT(a,lim,-1);
ClearArray(a,Mod+1,lim);
ClearArray(BB,0,lim);
}
void InvPoly(int n,int *a)
{
static int AA[MAXN],BB[MAXN],ans[MAXN];
ans[0]=Pow(a[0],P-2);
int len,N;
for(len=1;len<=(n<<1);len<<=1)
{
N=len<<1;
CopyArray(AA,a,0,len-1);
CopyArray(BB,ans,0,len-1);
NTT(AA,N,1);
NTT(BB,N,1);
for(int i=0;i<=N;i++)
{
ans[i]=Mul2(Sub2(2,Mul2(AA[i],BB[i])),BB[i]);
}
NTT(ans,N,-1);
ClearArray(ans,len,N);
}
ClearArray(AA,0,len);
ClearArray(BB,0,len);
CopyArray(a,ans,0,n);
ClearArray(ans,0,len);
}
void DivPoly(int n,int *a,int m,int *b)
{
static int AA[MAXN],BB[MAXN];
int Len=n-m;
RevArray(b,0,m);
CopyArray(AA,b,0,Len);
RevArray(b,0,m);
RevArray(a,0,n);
CopyArray(BB,a,0,Len);
RevArray(a,0,n);
InvPoly(Len,AA);
MulPoly(Len,AA,Len,BB);
RevArray(AA,0,Len);
MulPoly(m,b,Len,AA);
for(int i=0;i<m;i++)
{
b[i]=(a[i]-b[i]+P)%P;
}
ClearArray(b,m,m+Len);
CopyArray(a,AA,0,Len);
ClearArray(a,Len+1,n);
}
void DaoPoly(int n,int *a)
{
for(int i=1;i<n;i++)
{
a[i-1]=Mul2(a[i],i);
}
a[n]=0;
}
int Inv[MAXN];
void Init(int n)
{
Inv[1]=1;
for(int i=2;i<=n;i++)
{
Inv[i]=Mul2(Inv[P%i],P-P/i);
}
}
void JiPoly(int n,int *a)
{
for(int i=n;i>=1;i--)
{
a[i]=Mul2(a[i-1],Inv[i]);
}
a[0]=0;
}
void LnPoly(int n,int *a)
{
static int AA[MAXN];
CopyArray(AA,a,0,n);
InvPoly(n,AA);
DaoPoly(n,a);
MulPoly(n,a,n,AA,n);
JiPoly(n-1,a);
ClearArray(AA,0,n);
}
void ExpPoly(int n,int *a)
{
static int AA[MAXN],ans[MAXN];
ans[0]=1;
int len=2,N=1;
while(N<=n)
{
N<<=1;
}
for(;len<=N;len<<=1)
{
CopyArray(AA,ans,0,len>>1);
LnPoly(len,AA);
for(int i=0;i<len;i++)
{
AA[i]=Sub2(a[i],AA[i]);
}
Add1(AA[0],1);
MulPoly(len,ans,len,AA,len<<1);
ClearArray(ans,len,len<<1);
ClearArray(AA,len,len<<1);
}
ClearArray(AA,0,len);
CopyArray(a,ans,0,n);
ClearArray(ans,0,len);
}
/*
void ModPoly(int n,int *a,int m,int *b)
{
int Invb=Pow(b[m],P-2);
for(int i=n;i>=m;i--)
{
int Shang=Mul2(a[i],Invb);
for(int j=i,k=m;j>=i-m;j--,k--)
{
Sub1(a[j],Mul2(b[k],Shang));
}
}
}
*/
void ModPoly(int n,int *a,int m,int *b)
{
static int AA[MAXN],BB[MAXN],CC[MAXN];
int Len=n-m;
RevArray(b,0,m);
CopyArray(AA,b,0,Len);
RevArray(b,0,m);
RevArray(a,0,n);
CopyArray(BB,a,0,Len);
RevArray(a,0,n);
InvPoly(Len,AA);
MulPoly(Len,AA,Len,BB);
RevArray(AA,0,Len);
CopyArray(CC,b,0,m);
MulPoly(m,CC,Len,AA);
for(int i=0;i<m;i++)
{
Sub1(a[i],CC[i]);
}
ClearArray(a,m,n);
ClearArray(CC,0,n);
ClearArray(AA,0,Len<<1);
ClearArray(BB,0,Len<<1);
}
int F[20][MAXN<<2];
vector<int>Vec[MAXN<<1];
int tot,Root;
int ls[MAXN<<1],rs[MAXN<<1];
void Solve1(int &u,int l,int r,int *X)
{
static int A[MAXN<<2],B[MAXN<<2],C[MAXN<<2];
u=++tot;
int Len=1;
while(Len<=(r-l+1))
{
Len<<=1;
}
Vec[u].resize(Len);
if(l==r)
{
Vec[u][0]=X[l]?Sub2(P,X[l]):0;
Vec[u][1]=1;
return;
}
int Mid=l+r>>1;
Solve1(ls[u],l,Mid,X);
Solve1(rs[u],Mid+1,r,X);
int Size=Vec[ls[u]].size();
for(int i=0;i<Size;i++)
{
A[i]=Vec[ls[u]][i];
}
ClearArray(A,Size,Len);
Size=Vec[rs[u]].size();
for(int i=0;i<Size;i++)
{
B[i]=Vec[rs[u]][i];
}
ClearArray(B,Size,Len);
MulPoly(Len,A,B);
for(int i=0;i<Len;i++)
{
Vec[u][i]=A[i];
}
ClearArray(A,0,Len);
ClearArray(B,0,Len);
}
int G[MAXN<<2];
void Solve2(int depth,int u,int l,int r,int n,int *ans,int *X)
{
int Len=r-l+1;
for(int i=0;i<=Len;i++)
{
G[i]=Vec[u][i];
}
CopyArray(F[depth],F[depth-1],0,max(n,Len));
if(n>=Len)
{
ModPoly(n,F[depth],Len,G);
}
ClearArray(G,0,Len);
if(Len<=400)
{
for (int i=l;i<=r;i++)
{
int Now=1;
for (int j=0;j<=min(n,Len);j++)
{
Add1(ans[i],Mul2(Now,F[depth][j]));
Mul1(Now,X[i]);
}
}
return;
}
int Mid=l+r>>1;
depth++;
Solve2(depth,ls[u],l,Mid,Len-1,ans,X);
Solve2(depth,rs[u],Mid+1,r,Len-1,ans,X);
}
void DuoDianQiuZhi(int n,int *a,int m,int *X,int *ans)
{
for(int i=0;i<=n;i++)
{
F[0][i]=a[i];
}
for(int i=0;i<=m;i++)
{
G[i]=Vec[Root][i];
}
Solve2(1,Root,1,m,n,ans,X);
}
int DeltaDao[MAXN],ans[MAXN];
int HH[20][MAXN<<2];
void FenZhiNTT2(int depth,int u,int l,int r,int *Y)
{
static int AA[MAXN],BB[MAXN];
if(l==r)
{
HH[depth][0]=Y[l];
return;
}
int Mid=l+r>>1;
int Size;
FenZhiNTT2(depth+1,ls[u],l,Mid,Y);
Size=Vec[rs[u]].size();
for(int i=0;i<Size;i++)
{
AA[i]=Vec[rs[u]][i];
}
CopyArray(HH[depth],HH[depth+1],0,Mid-l+1);
ClearArray(HH[depth+1],0,Mid-l+1);
MulPoly(Mid-l+1,HH[depth],Size,AA);
ClearArray(AA,0,Size);
FenZhiNTT2(depth+1,rs[u],Mid+1,r,Y);
Size=Vec[ls[u]].size();
for(int i=0;i<Size;i++)
{
AA[i]=Vec[ls[u]][i];
}
CopyArray(BB,HH[depth+1],0,r-Mid);
ClearArray(HH[depth+1],0,r-Mid);
MulPoly(r-Mid,BB,Size,AA);
AddArray(HH[depth],BB,0,r-l+1);
ClearArray(AA,0,r-l+1);
ClearArray(BB,0,r-l+1);
}
void FastLag(int n,int *x,int *y,int *ans)
{
Solve1(Root,1,n,x);
int Size=Vec[Root].size();
for(int i=0;i<Size;i++)
{
DeltaDao[i]=Vec[Root][i];
}
DaoPoly(Size,DeltaDao);
DuoDianQiuZhi(n,DeltaDao,n,x,ans);
for(int i=1;i<=n;i++)
{
Mul1(y[i],Pow(ans[i]));
}
FenZhiNTT2(1,Root,1,n,y);
for(int i=0;i<=n-1;i++)
{
ans[i]=(HH[1][i]%P+P)%P;
}
}
int main()
{
}
Wonder Egg Priority
本文来自博客园，作者：0htoAi，转载请注明原文链接：https://www.cnblogs.com/0htoAi/p/16421006.html
