300+行完成几乎所有多项式操作...
(不得不说指针真是个好东西)
#include<bits/stdc++.h>
#define Set(a,b) memset(a,b,sizeof(a))
#define Clear(a,_begin_,_end_) for(int i=_begin_;i<_end_;++i) a[i]=0
#define Input_Array(a,_begin_,_end_) for(int i=_begin_;i<_end_;++i) init(a[i])
#define __ NULL
using namespace std;
const int N=1e5+10,MAXN=N<<2,MAXM=1e7;
typedef vector<int> Poly;
template <typename T> inline void init(T&x){
x=0;char ch=getchar();bool t=0;
for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=1;
for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
if(t) x=-x;return;
}typedef double db;
typedef long long ll;
int Inv[MAXN],rader[MAXN];
const int mod=998244353,phi=998244352,SIZE=sizeof(rader),inv2=499122177;
template<typename T>inline void Inc(T&x,int y,const int mod=998244353){x+=y;if(x>=mod) x-=mod;return;}
template<typename T>inline void Dec(T&x,int y,const int mod=998244353){x-=y;if(x < 0) x+=mod;return;}
template<typename T>inline int fpow(int x,T k,const int mod=998244353){int ret=1;for(;k;k>>=1,x=(ll)x*x%mod) if(k&1) ret=(ll)ret*x%mod;return ret;}
inline int Sum(int x,int y,const int mod=998244353){x+=y;if(x>=mod) return x-mod;return x;}
inline int Dif(int x,int y,const int mod=998244353){x-=y;if(x < 0 ) return x+mod;return x;}
inline int Init(int n){int len=1,up=-1;while(len<=n) len<<=1,++up;for(int i=0;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);return len;}
inline void Calc_Inversion(){Inv[1]=1;for(int i=2;i<MAXN;++i) Inv[i]=(ll)(mod-mod/i)*Inv[mod%i]%mod;return;}
namespace FFT{
const db PI=acos(-1);
struct Complex{
db x,y;Complex(db _x=0.0,db _y=0.0) {x=_x,y=_y;}
inline Complex operator +(const Complex &B){return Complex(x+B.x,y+B.y);}
inline Complex operator -(const Complex &B){return Complex(x-B.x,y-B.y);}
inline Complex operator /(const db a){return Complex(x/a,y/a);}
inline Complex operator *(const Complex &B){return Complex(x*B.x-y*B.y,x*B.y+y*B.x);}
}w[MAXN];
inline void Calcw(int len){for(int i=1;i<len;i<<=1) for(int j=0;j<i;++j) w[len/i*j]=Complex(cos(PI/i*j),sin(PI/i*j));}
inline void FFT(Complex*A,int n,int f){
for(int i=0;i<n;++i) if(rader[i]>i) swap(A[i],A[rader[i]]);
for(int i=1;i<n;i<<=1)
for(int j=0,p=i<<1;j<n;j+=p)
for(int k=0;k<i;++k) {
Complex W= ~f ? w[n/i*k]:Complex(w[n/i*k].x,-w[n/i*k].y);
Complex X=A[j|k],Y=W*A[j|k|i];
A[j|k]=X+Y,A[j|k|i]=X-Y;
}if(!~f) for(int i=0;i<n;++i) A[i]=A[i]/(db)n;
}
inline void Mul(int*a,int*b,int*c,int n,int m) {
int L=n+m-1;int len=Init(L);static Complex A[MAXN],B[MAXN];
for(int i=0;i<len;++i) A[i]=Complex(),B[i]=Complex();
for(int i=0;i<n;++i) A[i]=Complex((db)a[i],0);for(int i=0;i<m;++i) B[i]=Complex((db)b[i],0);
Calcw(len);FFT(A,len,1),FFT(B,len,1);
for(int i=0;i<len;++i) A[i]=A[i]*B[i];
FFT(A,len,-1);for(int i=0;i<L;++i) c[i]=(int)(A[i].x+0.5);
return;
}
namespace MTT{
inline void Mul(int*A,int*B,int*C,int len,int mod){
static Complex A1[MAXN],A2[MAXN],B1[MAXN],B2[MAXN];int MO=sqrt(mod);
for(int i=0;i<len;++i) A1[i]=Complex(A[i]/MO,0.0),B1[i]=Complex(A[i]%MO,0.0),A2[i]=Complex(B[i]/MO,0.0),B2[i]=Complex(B[i]%MO,0.0);
FFT(A1,len,1),FFT(A2,len,1),FFT(B1,len,1),FFT(B2,len,1);
for(int i=0;i<len;++i) {Complex X;
X=A1[i]*A2[i],A2[i]=A2[i]*B1[i];
B1[i]=B1[i]*B2[i];B2[i]=B2[i]*A1[i];
A1[i]=X,A2[i]=A2[i]+B2[i];
}int MOD=MO*MO%mod;
FFT(A1,len,-1),FFT(B1,len,-1),FFT(A2,len,-1);
for(int i=0;i<len;++i) {
int X=(ll)(A1[i].x+0.5)%mod,Y=(ll)(B1[i].x+0.5)%mod,Z=(ll)(A2[i].x+0.5)%mod;
int ans=(ll)MOD*X%mod;Inc(ans,(ll)MO*Z%mod);Inc(ans,Y);C[i]=ans;
}return;
}
}
}
namespace Any_Length_DFT{
int n,mod,g,len;int wn[N],iwn[N];
inline void Getroot(){
int phi=mod-1;int x=phi;static int pri[50],cnt=0;
for(int i=2;i*i<=x;++i) if(x%i==0) {pri[++cnt]=i,x/=i;while(x%i==0) x/=i;}
for(g=2;;++g){
bool fl=1;
for(int i=1;i<=cnt;++i)if(fpow(g,phi/pri[i],mod)==1) {fl=0;break;}
if(fl)return;
}
}
inline void Prework(int _len,int _mod){
n=_len,mod=_mod;Getroot();
wn[0]=iwn[0]=1,wn[1]=fpow(g,(mod-1)/n,mod),iwn[1]=fpow(wn[1],mod-2,mod);
for(int i=2;i<n;++i) wn[i]=(ll)wn[i-1]*wn[1]%mod,iwn[i]=(ll)iwn[i-1]*iwn[1]%mod;
len=Init(3*n-3);FFT::Calcw(len);
return;
}
inline int Comb(int x){return (ll)x*(x-1)/2%n;}
inline void DFT(int*A,int f){
int m=2*n-1;static int F[MAXN],G[MAXN];
if(~f) {
for(int i=0;i<n;++i) F[i]=(ll)A[i]*iwn[Comb(i)]%mod;Clear(F,n,len);
for(int i=0;i<m;++i) G[i]=wn[Comb(i)];Clear(G,m,len);
}
else {
for(int i=0;i<n;++i) F[i]=(ll)A[i]*wn[Comb(i)]%mod;Clear(F,n,len);
for(int i=0;i<m;++i) G[i]=iwn[Comb(i)];Clear(G,m,len);
}reverse(F,F+n);FFT::MTT::Mul(F,G,F,len,mod);
for(int k=0,i=n-1;i<m;++i,++k) {
if(~f) A[k]=(ll)F[i]*iwn[Comb(k)]%mod;
else A[k]=(ll)F[i]* wn[Comb(k)]%mod;
}
if(!~f) for(int i=0,inv=fpow(n,mod-2,mod);i<n;++i) A[i]=(ll)A[i]*inv%mod;
return;
}
}
namespace NTT{
int wn[30],iwn[30];
inline void Calcw(){for(int i=0;i<30;++i) wn[i]=fpow(3,phi/(1<<i)),iwn[i]=fpow(wn[i],mod-2);}
inline void NTT(int*A,int n,int f){
for(int i=0;i<n;++i) if(rader[i]>i) swap(A[i],A[rader[i]]);
for(int i=1,h=1;i<n;++h,i<<=1){
int W=(~f)? wn[h]:iwn[h];
for(int j=0,p=i<<1;j<n;j+=p){
for(int w=1,k=0;k<i;++k,w=(ll)w*W%mod){
int X=A[j|k],Y=(ll)w*A[j|k|i]%mod;
A[j|k]=Sum(X,Y),A[j|k|i]=Dif(X,Y);
}
}
}if(!~f) for(int i=0;i<n;++i) A[i]=(ll)A[i]*Inv[n]%mod;
return;
}
inline void Mul(int*a,int*b,int*c,int n,int m) {
int L=n+m-1;int len=Init(L);static int A[MAXN],B[MAXN];
for(int i=0;i<n;++i) A[i]=a[i];for(int i=0;i<m;++i) B[i]=b[i];
Clear(A,n,len);Clear(B,m,len);NTT(A,len,1),NTT(B,len,1);
for(int i=0;i<len;++i) c[i]=(ll)A[i]*B[i]%mod;
NTT(c,len,-1);return;
}
inline void CDQ(int l,int r,int*F,int*G){//分治 FFT
if(l==r) return;int mid=(l+r)>>1;
CDQ(l,mid,F,G);static int R[MAXN];
Mul(F+l,G+1,R,mid-l+1,r-l+1);
for(int i=mid-l,j=mid+1;j<=r;++j,++i) Inc(F[j],R[i]);
CDQ(mid+1,r,F,G);return;
}
inline void Poly_Inv(int*F,int*I,int n){
if(n==1) {memset(I,0,SIZE);I[0]=fpow(F[0],mod-2);return;}
Poly_Inv(F,I,(n+1)>>1);int L=n<<1,len=Init(L);
static int A[MAXN];for(int i=0;i<n;++i) A[i]=F[i];Clear(A,n,len);
NTT(I,len,1);NTT(A,len,1);
for(int i=0;i<len;++i) I[i]=Dif(Sum(I[i],I[i]),(ll)I[i]*I[i]%mod*A[i]%mod);
NTT(I,len,-1);Clear(I,n,len);return;
}
inline void Direv(int*A,int n){for(int i=0;i<n;++i) A[i]=(ll)A[i+1]*(i+1)%mod;A[n-1]=0;return;}
inline void Integ(int*A,int n){for(int i=n-1;i;--i) A[i]=(ll)A[i-1]*Inv[i]%mod; A[0]=0;return;}
inline void Poly_Ln(int*F,int*L,int n) {
static int A[MAXN],B[MAXN];for(int i=0;i<n;++i) A[i]=F[i];
Direv(A,n);Poly_Inv(F,B,n);Mul(A,B,L,n,n);
int TL=(n<<1)-1;Clear(L,n,TL);Integ(L,n);return;
}
inline void Poly_Exp(int*F,int*E,int n){
if(n==1) {memset(E,0,SIZE);E[0]=1;return;} static int A[MAXN];
Poly_Exp(F,E,(n+1)>>1);Poly_Ln(E,A,n);
for(int i=0;i<n;++i) A[i]=Dif(F[i],A[i]);Inc(A[0],1);
Mul(E,A,E,n,n);int TL=(n<<1)-1;Clear(E,n,TL);return;
}
inline int Get_sqrt(int n){//二次剩余
if(n<=1) return n;if((int)sqrt(n)*(int)sqrt(n)==n) return (int)sqrt(n);
srand(time(NULL));int a,delta,D=(mod-1)>>1;
while(233) {a=rand()%mod;delta=Dif((ll)a*a%mod,n);if(fpow(delta,D)!=1) break;}
D=(mod+1)>>1;int b=1;int c=1,d=0;
while(D){
if(D&1) {int _c=c;c=((ll)a*c+(ll)b*d%mod*delta)%mod,d=((ll)a*d+(ll)b*_c)%mod;}
int _a=a;a=((ll)a*a+(ll)b*b%mod*delta)%mod;b=(2ll*b*_a)%mod;D>>=1;
}c=min(c,mod-c);
return c;
}
inline void Poly_Sqrt(int*F,int*S,int n){
if(n==1) {memset(S,0,SIZE);S[0]=Get_sqrt(F[0]);return;}
Poly_Sqrt(F,S,(n+1)>>1);static int A[MAXN],B[MAXN];Poly_Inv(S,A,n);
int len=Init(n<<1);
for(int i=0;i<n;++i) B[i]=F[i];Clear(B,n,len);
NTT(S,len,1);NTT(A,len,1);NTT(B,len,1);
for(int i=0;i<len;++i) S[i]=((ll)S[i]*S[i]%mod+B[i])*inv2%mod*A[i]%mod;
NTT(S,len,-1);Clear(S,n,len);
return;
}
inline void Poly_Fpow(int*F,int k,int n){//常数项为 0 就先给他处理一下
int base=1;static int R[MAXN];
if(F[0]!=1) {base=fpow(F[0],k);int inv=fpow(F[0],mod-2);for(int i=0;i<n;++i) F[i]=(ll)F[i]*inv%mod;}
Poly_Ln(F,F,n);for(int i=0;i<n;++i) F[i]=(ll)F[i]*k%mod;
Poly_Exp(F,R,n);for(int i=0;i<n;++i) F[i]=(ll)R[i]*base%mod;
return;
}
inline void Poly_Mod(int*A,int*B,int*Q,int*R,int n,int m){
if(n<m) {for(int i=0;i<=n;++i) R[i]=A[i];return;}
static int C[MAXN],D[MAXN],E[MAXN];
int r=(n-m)<<1;int len=1;while(len<=r)len<<=1;
for(int i=0;i<=n-m;++i) D[i]=A[n-i];Clear(D,n-m+1,len);
for(int i=0;i<=m;++i) E[i]=B[m-i];Clear(E,m+1,len);
Poly_Inv(E,C,n-m+1);len=Init(r);
NTT(C,len,1),NTT(D,len,1);
for(int i=0;i<len;++i) C[i]=(ll)C[i]*D[i]%mod;
NTT(C,len,-1);reverse(E,E+1+m),reverse(C,C+1+n-m);
if(Q) for(int i=0;i<=n-m;++i) Q[i]=C[i];
len=Init(n);Clear(C,n-m+1,len);Clear(E,m+1,len);
NTT(C,len,1),NTT(E,len,1);
for(int i=0;i<len;++i) C[i]=(ll)C[i]*E[i]%mod;
NTT(C,len,-1);for(int i=0;i<m;++i) R[i]=Dif(A[i],C[i]);
return;
}
const int img=86583718;
inline void Poly_Sin(int*F,int*Sin,int n){
static int L[MAXN],R[MAXN],XL[MAXN],XR[MAXN];
for(int i=0;i<n;++i) L[i]=(ll)F[i]*img%mod,R[i]=mod-L[i];
Poly_Exp(L,XL,n);Poly_Exp(R,XR,n);
int inv=fpow(2ll*img%mod,mod-2);
for(int i=0;i<n;++i) Sin[i]=(ll)Dif(XL[i],XR[i])*inv%mod;
return;
}
inline void Poly_Cos(int*F,int*Cos,int n){
static int L[MAXN],R[MAXN],XL[MAXN],XR[MAXN];
for(int i=0;i<n;++i) L[i]=(ll)F[i]*img%mod,R[i]=mod-L[i];
Poly_Exp(L,XL,n);Poly_Exp(R,XR,n);
for(int i=0;i<n;++i) Cos[i]=(ll)Sum(XL[i],XR[i])*inv2%mod;
return;
}
inline void Poly_Arcsin(int*F,int*Arcsin,int n){
static int A[MAXN],S[MAXN],I[MAXN];
for(int i=0;i<n;++i) A[i]=I[i]=F[i];Direv(A,n);
Mul(I,I,I,n,n);int TL=(n<<1)-1;Clear(I,n,TL);
for(int i=0;i<n;++i) I[i]=Dif(mod,I[i]);Inc(I[0],1);
Poly_Sqrt(I,S,n);Poly_Inv(S,I,n);
Mul(A,I,Arcsin,n,n);Integ(Arcsin,n);return;
}
inline void Poly_Arctan(int*F,int*Arctan,int n){
static int A[MAXN],B[MAXN],I[MAXN];
for(int i=0;i<n;++i) A[i]=B[i]=F[i];Direv(A,n);
Mul(B,B,B,n,n);int TL=(n<<1)-1;Clear(B,n,TL);
Inc(B[0],1);Poly_Inv(B,I,n);
Mul(A,I,Arctan,n,n);Integ(Arctan,n);return;
}
#define ls (u<<1)
#define rs (u<<1|1)
int POOL[MAXM],cnt=0;
int *P[MAXN],*F[MAXN],*G[MAXN],X[N],*val;
inline int* Neospace(int len){int*ret=&POOL[cnt];cnt+=len;return ret;}
inline int Calc(int*F,int n,const int x){int X=1,ret=0;for(int i=0;i<n;++i) Inc(ret,(ll)F[i]*X%mod),X=(ll)X*x%mod;return ret;}
inline void Divide(int u,int l,int r){P[u]=NULL;//分治+NTT , 保留了中间的结果
if(l==r) {P[u]=Neospace(2);P[u][0]=mod-X[l];P[u][1]=1;return;}
static int L[MAXN],R[MAXN];int mid=(l+r)>>1;int LS=ls,RS=rs;
Divide(LS,l,mid);Divide(RS,mid+1,r);
int n=r-l+2,nl=mid-l+2,nr=r-mid+1;
for(int i=0;i<nl;++i) L[i]=P[LS][i];
for(int i=0;i<nr;++i) R[i]=P[RS][i];
Mul(L,R,L,nl,nr);P[u]=Neospace(n);
for(int i=0;i<n;++i) P[u][i]=L[i];
return;
}
inline void Poly_Evaluate(int u,int l,int r){// 多点求值
if(r-l+1<=1000) {for(int i=l;i<=r;++i) val[i]=Calc(F[u],r-l+1,X[i]);return;}// 长度小的暴力计算
int mid=(l+r)>>1,LS=ls,RS=rs;
int n=r-l+1,nl=mid-l+1,nr=r-mid;
static int R[MAXN];
F[LS]=Neospace(nl),F[RS]=Neospace(nr);
Poly_Mod(F[u],P[LS],__,R,n-1,nl);
for(int i=0;i<nl;++i) F[LS][i]=R[i];
Poly_Mod(F[u],P[RS],__,R,n-1,nr);
for(int i=0;i<nr;++i) F[RS][i]=R[i];
Poly_Evaluate(LS,l,mid);
Poly_Evaluate(RS,mid+1,r);
return;
}
inline void Solve_Evaluation(int*A,int*_X,int*ans,int n,int m){
cnt=0;F[1]=Neospace(m);val=ans;
for(int i=1;i<=m;++i) X[i]=_X[i];Divide(1,1,m);
Poly_Mod(A,P[1],__,F[1],n-1,m);//注意最开始的多项式也要取一次模
Poly_Evaluate(1,1,m);return;
}
inline void Poly_Interpolate(int u,int l,int r){// 多项式快速插值
if(l==r) {G[u]=Neospace(1);G[u][0]=val[l];return;}
int mid=(l+r)>>1,LS=ls,RS=rs,n=r-l+1,nl=mid-l+1,nr=r-mid;
Poly_Interpolate(LS,l,mid);
Poly_Interpolate(RS,mid+1,r);
static int L[MAXN],R[MAXN];
G[u]=Neospace(n);
Mul(G[LS],P[RS],L,nl,nr+1);
Mul(G[RS],P[LS],R,nr,nl+1);
for(int i=0;i<n;++i) G[u][i]=Sum(L[i],R[i]);
return;
}
inline void Solve_Interpolation(int*_X,int*Y,int*ans,int n){
cnt=0;for(int i=1;i<=n;++i) X[i]=_X[i];
Divide(1,1,n);static int Val[N];val=Val;int hcnt=cnt;
F[1]=Neospace(n);for(int i=0;i<=n;++i) F[1][i]=P[1][i];
Direv(F[1],n+1);Poly_Evaluate(1,1,n);
Clear(POOL,hcnt,cnt);cnt=hcnt;
for(int i=1;i<=n;++i) Val[i]=(ll)Y[i]*fpow(Val[i],mod-2)%mod;
Poly_Interpolate(1,1,n);for(int i=0;i<n;++i) ans[i]=G[1][i];
return;
}
}
namespace FWT{
inline void FWT(int*A,int n,int f,int opt){
if(opt==1) {// and
for(int i=1;i<n;i<<=1)
for(int j=0,p=i<<1;j<n;j+=p)
for(int k=0;k<i;++k)
f? Inc(A[j|k|i],A[j|k]):Dec(A[j|k|i],A[j|k]);
}else if(opt==2){// or
for(int i=1;i<n;i<<=1)
for(int j=0,p=i<<1;j<n;j+=p)
for(int k=0;k<i;++k)
f? Inc(A[j|k],A[j|k|i]):Dec(A[j|k],A[j|k|i]);
}else {// xor
for(int i=1;i<n;i<<=1)
for(int j=0,p=i<<1;j<n;j+=p)
for(int k=0;k<i;++k){
int X=A[j|k],Y=A[j|k|i];
A[j|k]=Sum(X,Y),A[j|k|i]=Dif(X,Y);
}
if(!f) for(int i=0;i<n;++i) A[i]=(ll)A[i]*Inv[n]%mod;
}
}
inline void Mul(int*A,int*B,int*C,int n,int opt){
static int X[MAXN],Y[MAXN];
for(int i=0;i<n;++i) X[i]=A[i],Y[i]=B[i];
FWT(X,n,1,opt);FWT(Y,n,1,opt);
for(int i=0;i<n;++i) C[i]=(ll)X[i]*Y[i]%mod;
FWT(C,n,0,opt);return;
}
}
int main()
{
Calc_Inversion();NTT::Calcw();
return 0;
}