代码

#pragma GCC optimize("O3,unroll-loops")
#include<bits/stdc++.h>
#define LL long long
#define LLL __int128
#define uint unsigned
#define ldb long double
#define uLL unsigned long long
using namespace std;
typedef vector<int> poly;
typedef vector<int> Vec;
typedef vector<Vec> Mat;
typedef tuple<poly,poly> Vec2;
typedef tuple<poly,poly,poly,poly> Mat2;
mt19937 rng(chrono::system_clock::now().time_since_epoch().count());
const int Mod=998244353,G=3,img=86583718,winv=ULLONG_MAX%Mod+1;
const int Lim=64,rinv=2-Mod;
const LL Mod2=(LL)Mod*Mod;
poly Grt({1}),Grt0,Grt1,frc({1,1}),inv({0,1}),ivf({1,1});vector<poly>Tr,Ts;
inline int fadd(const int&x,const int&y){int z=x+y;return z<0?z+Mod:z-Mod;}
inline int fsub(const int&x,const int&y){int z=x-y;return z<0?z+Mod:z-Mod;}
inline int reduce(const LL&x){return (x>>32)-(((LL)((int)(x*rinv))*Mod)>>32);}
inline int qpow(int x,int y,int z=1){
	for(;y;(y>>=1)&&(x=(LL)x*x%Mod))if(y&1)z=(LL)z*x%Mod;return z;
}
inline void Init(const int&n){
	for(int i=frc.size();i<=n;++i)
		frc.emplace_back((LL)frc.back()*i%Mod),
		inv.emplace_back(Mod-Mod/i*(LL)inv[Mod%i]%Mod),
		ivf.emplace_back((LL)ivf.back()*inv.back()%Mod);
}
inline int Binom(const int&n,const int&m){
	if(n<m||m<0)return 0;
	return Init(n),(LL)frc[n]*ivf[m]%Mod*ivf[n-m]%Mod;
}
inline int BSGS(int x){
	unordered_map<int,int>mp;
	int S=sqrtl(Mod)+1,j=1,q=1;
	for(int i=0;i<S;++i,j=(LL)j*G%Mod)mp[(LL)j*x%Mod]=i;
	for(int i=0;i<=S;++i,q=(LL)q*j%Mod)if(mp.count(q)&&i*S>=mp[q])return i*S-mp[q];
	return -1;
}
inline poly Add(poly P,poly Q){
    if(P.size()<Q.size())swap(P,Q);
    for(int i=Q.size();i--;)(P[i]+=Q[i])>=Mod&&(P[i]-=Mod);
    return P;
}
inline poly Sub(poly P,poly Q){
    if(P.size()<Q.size())P.resize(Q.size());
    for(int i=Q.size();i--;)(P[i]-=Q[i])<0&&(P[i]+=Mod);
    return P;
}
inline poly Neg(poly P){
	for(int i=P.size();i--;)P[i]&&(P[i]=Mod-P[i]);
	return P;
}
inline bool Empty(poly&P){
	for(;!P.empty()&&!P.back();P.pop_back());
	return P.empty();
}
inline int Eval(poly&P,int x){
	int z=0;
	for(int i=P.size();i--;)z=((LL)z*x+P[i])%Mod;
	return z;
}
inline Mat operator*(const Mat&x,const Mat&y){
	const int n=x.size(),m=y.size(),q=y[0].size();
	Mat z(n,Vec(q));
	for(int i=0;i<n;++i)for(int j=0;j<m;++j)if(x[i][j])
		for(int k=0;k<q;++k)z[i][k]=(z[i][k]+(LL)x[i][j]*y[j][k])%Mod;
	return z;
}
template<int T>inline void dft(poly&P){
	const int m=P.size(),n=m/T,n2=n/2;
    int L=Grt.size();
	if(L*2<n)for(Grt.resize(n2);L*2<n;L<<=1){
		const int cr=qpow(G,Mod/(L<<2));
		for(int i=0;i<L;++i)Grt[i+L]=(LL)cr*Grt[i]%Mod;
	}
	if(Grt0.size()<Grt.size()){
		const int pr=Grt0.size(),ed=Grt.size();
		Grt0.resize(ed),Grt1.resize(ed);
		for(int i=pr;i<ed;++i)
			Grt0[i]=reduce((LL)winv*Grt[i]),Grt1[i]=(uint)rinv*Grt0[i];
	}
	const auto dfs=[&](auto&&self,int i,int k){
		if(k==1)return;
		const int R=k*T,R4=R/4;
		int*P0=&P[i*R],*P1=P0+R4,*P2=P1+R4,*P3=P2+R4;
		const int g20=Grt1[i+i],gg0=Grt0[i+i];
		const int g21=Grt1[i+i+1],gg1=Grt0[i+i+1];
		const int g2=Grt1[i],gg=Grt0[i];
		for(int j=0,z,m,z0,z1,z2,z3,s0,d0,s1,d1;j<R4;++j){
			z=P2[j],m=(uint)g2*z;
			z0=((LL)z*gg-(LL)m*Mod)>>32;
			z=P3[j],m=(uint)g2*z;
			z1=((LL)z*gg-(LL)m*Mod)>>32;
			s0=fadd(P0[j],z0),d0=fsub(P0[j],z0);
			s1=P1[j]+z1,d1=P1[j]-z1;
			z=s1,m=(uint)g20*z;
			z2=((LL)z*gg0-(LL)m*Mod)>>32;
			z=d1,m=(uint)g21*z;
			z3=((LL)z*gg1-(LL)m*Mod)>>32;
			P0[j]=fadd(s0,z2);
			P1[j]=fsub(s0,z2);
			P2[j]=fadd(d0,z3);
			P3[j]=fsub(d0,z3);
		}
		self(self,i*4+0,k/4);
		self(self,i*4+1,k/4);
		self(self,i*4+2,k/4);
		self(self,i*4+3,k/4);
	};
	if(__lg(n)&1){
		const int R=n*T,R2=n2*T;
		for(int i=0;i<R2;++i){
			const int j=i+R2,v=fsub(P[i],P[j]);
			P[i]=fadd(P[i],P[j]),P[j]=v;
		}
		dfs(dfs,0,n2),dfs(dfs,1,n2);
	}
	else dfs(dfs,0,n);
    for(int&i:P)(i<0&&(i+=Mod));
}
template<int T>inline void idft(poly&P){
	const int m=P.size(),n=m/T,ni=reduce((LL)qpow(n,Mod-2)*winv),n2=n/2;
	poly rev(n);
	for(int i=1;i<n;++i)rev[i]=(rev[i>>1]>>1)+(i&1)*n2;
	for(int&i:P)i=reduce((LL)i*ni);
	for(int i=1;i<n;++i)if(i<rev[i])swap_ranges(P.begin()+T*i,P.begin()+T*(i+1),P.begin()+T*rev[i]);
	for(int i=1;i<n2;++i)swap_ranges(P.begin()+T*i,P.begin()+T*(i+1),P.begin()+T*(n-i));
	dft<T>(P);
	for(int i=1;i<n;++i)if(i<rev[i])swap_ranges(P.begin()+T*i,P.begin()+T*(i+1),P.begin()+T*rev[i]);
}
template<int T>inline void mul_mod(poly&P,poly&Q){
	const int m=P.size(),n=m/T;
	dft<T>(P),dft<T>(Q);
	for(int i=n;i--;){
		vector<uLL>ans(T*2);
		for(int j=0;j<T;++j){
			if(j>9&&j%9==8)for(int k=j;k<j+T-10;++k)
				ans[k]-=Mod2*9*(ans[k]>>63);
			for(int k=0;k<T;++k)
				ans[j+k]+=(LL)P[i*T+j]*Q[i*T+k];
		}
		const int c=((i&1)?Mod-Grt[i>>1]:Grt[i>>1]);
		for(int j=T;j--;)P[i*T+j]=(ans[j]+ans[j+T]%Mod*c)%Mod;
	}
	idft<T>(P);
}
template<size_t... m>inline void Mul_solve(index_sequence<m...>,int n,poly&P,poly&Q){
	static void (*ptrs[])(poly&,poly&)={&mul_mod<m+1>...};
	ptrs[n-1](P,Q);
}
inline poly Mul(poly P,poly Q){
    if(P.empty()||Q.empty())return poly();
	const int pn=P.size(),qn=Q.size(),rn=pn+qn-1;
	int b=1;while((Lim<<b)<rn)++b;
	const int T=((rn-1)>>b)+1,m=T<<b;
	P.resize(m),Q.resize(m);
	Mul_solve(make_index_sequence<Lim>{},T,P,Q);
	return P.resize(rn),P;
}
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
	poly P=Mul(poly(a,a+n+1),poly(b,b+m+1));
    copy(P.begin(),P.end(),c);
}

评测结果