不想多说了，看网上的题解吧，我大概说下思路。
首先考察Stirling的意义，然后求出递推式，变成卷积的形式。
然后发现贡献是一定的，我们可以分治+NTT。
也可以直接求逆（我不会啊啊啊啊啊）
#include <map>
#include <cmath>
#include <queue>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)
#define ll long long
#define mp make_pair
#define md 998244353
#define maxn 400005
int n;
ll ans,a[maxn],b[maxn],f[maxn],fac[maxn],inv[maxn],rev[maxn];
ll ksm(ll a,ll b)
{
ll ret=1;
for (;b;b>>=1,a=a*a%md) if (b&1) ret=ret*a%md;
return ret;
}
void NTT(ll *a,int n,int flag)
{
F(i,0,n-1) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int m=2;m<=n;m<<=1)
{
int mid=m>>1;
ll wn=ksm(3,((md-1)/m*flag+md-1)%(md-1));
for (int i=0;i<n;i+=m)
{
ll w=1;
F(j,0,mid-1)
{
ll u=a[i+j],v=a[i+j+mid]*w%md;
a[i+j]=(u+v)%md; a[i+j+mid]=(u-v+md)%md;
w=w*wn%md;
}
}
}
if (flag==-1)
{
ll inv=ksm(n,md-2);
F(i,0,n-1) a[i]=a[i]*inv%md;
}
}
void solve(int l,int r)
{
if (l==r) return ;
int mid=(l+r)>>1;
solve(l,mid);
int n=mid-l+r-l,len=1;
while ((1<<len)<n) len++;
n=(1<<len);
F(i,0,n-1) a[i]=b[i]=0;
F(i,l,mid) a[i-l]=f[i];
F(i,1,r-l) b[i-1]=inv[i]*2%md;
F(i,1,n-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
NTT(a,n,1);NTT(b,n,1);
F(i,0,n-1) a[i]=a[i]*b[i]%md;
NTT(a,n,-1);
F(i,mid+1,r) f[i]=(f[i]+a[i-l-1])%md;
solve(mid+1,r);
}
int main()
{
scanf("%d",&n);
fac[0]=inv[0]=1;
F(i,1,n) fac[i]=fac[i-1]*i%md,inv[i]=ksm(fac[i],md-2);
f[0]=1;solve(0,n);
F(i,0,n) ans=(ans+f[i]*fac[i]%md)%md;
printf("%lld
",ans);
}
UPD：
搞出生成函数，然后把递推式化简一下，发现成了$F(x)=frac{1}{1-H(x)}$的样子了。
然后多项式求逆即可
#include <map>
#include <cmath>
#include <queue>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)
#define ll long long
#define mp make_pair
#define md 998244353
#define maxn 500005
#define g 3
int rev[maxn],a[maxn],b[maxn];
int ksm(int a,int b)
{
int ret=1;
for (;b;b>>=1,a=(ll)a*a%md) if (b&1) ret=(ll)ret*a%md;
return ret;
}
void NTT(int *x,int n,int flag)
{
F(i,0,n-1) if (rev[i]>i) swap(x[rev[i]],x[i]);
for (int m=2;m<=n;m<<=1)
{
int wn=ksm(g,((md-1)/m*flag+md-1)%(md-1));
for (int i=0;i<n;i+=m)
{
int w=1;
for (int j=0;j<(m>>1);++j)
{
int u=x[i+j],v=(ll)x[i+j+(m>>1)]*w%md;
x[i+j]=(u+v)%md; x[i+j+(m>>1)]=(u-v+md)%md;
w=(ll)w*wn%md;
}
}
}
if (flag==-1)
{
int inv=ksm(n,md-2);
F(i,0,n-1) x[i]=(ll)x[i]*inv%md;
}
}
void Get_Inv(int * a,int * b,int n)
{
static int tmp[maxn];
if (n==1){b[0]=ksm(a[0],md-2);return;}
Get_Inv(a,b,n>>1);
F(i,0,n-1) tmp[i]=a[i],tmp[n+i]=0;
int l=0;while(!(n>>l&1))l++;
F(i,0,(n<<1)-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
NTT(tmp,n<<1,1);NTT(b,n<<1,1);
F(i,0,(n<<1)-1)tmp[i]=(ll)b[i]*(2-(ll)tmp[i]*b[i]%md+md)%md;
NTT(tmp,n<<1,-1);F(i,0,n-1) b[i]=tmp[i],b[n+i]=0;
}
int fac[maxn],inv[maxn],n,m;
int main()
{
scanf("%d",&n); for (m=1;m<=n;m<<=1);
fac[0]=1;F(i,1,n)fac[i]=(ll)fac[i-1]*i%md;
inv[n]=ksm(fac[n],md-2);
D(i,n-1,0)inv[i]=(ll)inv[i+1]*(i+1)%md;
a[0]=1;F(i,1,n) a[i]=((md-inv[i])*2)%md;
Get_Inv(a,b,m);int ans=b[n];
for (int i=n;i;--i) ans=((ll)ans*i+b[i-1])%md;
printf("%d
",ans);
}