代码

#include <bits/stdc++.h>
typedef long long unsigned ll;
constexpr uint LOG(21),N(1<<LOG),MOD(998244353),PR(3),LIM(50),im(86583718),A(2309898375);
uint norm(uint x){return x>=MOD?x-MOD:x;}
uint reduce(ll x){return x%MOD;}
//uint reduce(ll x){x-=MOD*(((x>>28)*A)>>33);return x>=MOD?x-MOD:x;}
uint mul(uint x,uint y){return reduce((ll)x*y);}
uint ksm(uint x,uint y){uint ans(1);for(;y;y>>=1,x=mul(x,x))(y&1)&&(ans=mul(ans,x));return ans;}
uint inv(uint x){return ksm(x,MOD-2);}
typedef std::vector<uint> poly;
typedef const poly& ref;
inline uint sz(ref f){return f.size();}
uint pr[LOG+1][N];
struct _g
{
	_g(void)
	{
		uint w0(ksm(PR,(MOD-1)>>(LOG)));
		for(uint p(LOG);~p;--p,w0=mul(w0,w0))
		{
			pr[p][0]=1;
			for(uint i(1);i<(1u<<p);++i)
				pr[p][i]=mul(pr[p][i-1],w0);
		}
	}
}__g;
void DIF(uint *f,uint n,uint b)
{
	for(uint i(n>>1),*j,*k,*buf,p(b),nx,ny;i;i>>=1,--p)
	for(j=f;j<f+n;j+=i<<1)
	for(k=j,buf=pr[p];k<j+i;++k,++buf)
	{
		nx=norm(*k+k[i]),ny=norm(*k-k[i]+MOD);
		*k=nx,k[i]=mul(ny,*buf);
	}
}
void DIT(uint *f,uint n,uint)
{
	for(uint i(1),*j,*k,*buf,p(1),nx,ny;i<n;i<<=1,++p)
	for(j=f;j<f+n;j+=i<<1)
	for(k=j,buf=pr[p];k<j+i;++k,++buf)
	{
		nx=*k,ny=mul(k[i],*buf);
		*k=norm(nx+ny),k[i]=norm(nx-ny+MOD);
	}
	std::reverse(f+1,f+n);
	uint iv(MOD-(MOD-1)/n);
	for(uint *i(f);i<f+n;++i)
		*i=mul(iv,*i);
}
poly mul(ref f,uint k)
{
	poly g(f);
	for(auto &i:g)
		i=mul(i,k);
	return g;
}
poly DFT(ref f,uint b,bool flag)
{
	uint n(1<<b);
	static uint a[N];
	memcpy(a,f.data(),sz(f)*sizeof(uint));
	memset(a+sz(f),0,(n-sz(f))*sizeof(uint));
	(flag?DIT:DIF)(a,n,b);
	return poly(a,a+n);
}
void dft(poly &f,uint b,bool flag)
{
	uint n(1<<b);
	f.resize(n);
	(flag?DIT:DIF)(f.data(),n,b);
}
poly resize(ref f,uint n)
{
	if(sz(f)>=n)
		return poly(f.begin(),f.begin()+n);
	poly g(f);
	g.resize(n);
	return g;
}
poly mul(ref f,ref g)
{
	uint p(sz(f)+sz(g)-1);
	if(sz(f)<=10||sz(g)<=10||p<=LIM)
	{
		poly h(p);
		for(uint i(0),n(sz(f)),m(sz(g));i<n;++i)
		for(uint j(0);j<m;++j)
			h[i+j]=reduce(h[i+j]+(ll)f[i]*g[j]);
		return h;
	}
	uint b(0),n(1);
	while(n<p)
		++b,n<<=1;
	poly f0(DFT(f,b,0)),g0(DFT(g,b,0));
	for(uint i(0);i<n;++i)
		f0[i]=mul(f0[i],g0[i]);
	dft(f0,b,1);
	f0.resize(p);
	return f0;
}
void poly_multiply(uint *a,int n,uint *b,int m,uint *c)
{
	std::vector<uint> A(a,a+n+1),B(b,b+m+1),C(mul(A,B));
	memcpy(c,C.data(),(n+m+1)*sizeof(uint));
}

评测结果