因为是有标号计数，考虑用指数型生成函数来解决。
设有标号荒漠的生成函数为 (F(x))，有标号仙人掌的生成函数为 (f(x))，因为荒漠是仙人掌所组成的集合，得：
[F(x) = exp (f(x))
]
解决有标号仙人掌的计数问题后，就能解决有标号荒漠的计数问题。
先在仙人掌上确定一点来考虑，和该点相邻的边有在环上和不在环上两种情况，即与该点相邻的元素为独立边和环。独立边可看作单独的一个仙人掌，得独立边的生成函数为 (f(x))。环的大小至少为 (3)，环上每个点也都对应单独的一个仙人掌，正方向反方向都会计算一遍，所以要除 (2)，得环的生成函数为 (sumlimits_{i geqslant 2 } frac{f^i(x)}{2})。
该点相邻的所有元素可以看作一个集合，同时还要考虑确定的那个点，得：
[f(x) = x exp(f(x) + sum_{i geqslant 2 } frac{f^i(x)}{2}) = x exp(frac{f^2(x) - 2f(x)}{2f(x)-2})
]
通过牛顿迭代来求解。设已经求得 (f_0(x))，满足 (f_0(x) equiv f(x) pmod{x^n})，设 (g(f(x))= x exp(frac{f^2(x) - 2f(x)}{2f(x)-2}) - f(x),h(x) = x exp(frac{f_0^2(x) - 2f_0(x)}{2f_0(x)-2}))，得：
[egin{aligned} f(x) & equiv f_0(x) - frac{g(f_0(x))}{{g}'(f_0(x))} pmod{x^{2n}} \ & equiv f_0(x) - frac{h(x)-f_0(x)}{h(x){(frac{f_0^2(x) - 2f_0(x)}{2f_0(x)-2})}'-1} pmod{x^{2n}} \ & equiv f_0(x) - frac{2h(x)-2f_0(x)}{h(x)(frac{1}{(f_0(x)-1)^2}+1)-2} pmod{x^{2n}}end{aligned}
]
因为之前的推导是确定了仙人掌上的一个点，所以对于 (n) 个点的仙人掌，一个方案会计算 (n) 次，要除以 (n)。
(code:)
#include<bits/stdc++.h>
#define maxn 800010
#define P 998244353
#define G 3
using namespace std;
typedef long long ll;
template<typename T> inline void read(T &x)
{
x=0;char c=getchar();bool flag=false;
while(!isdigit(c)){if(c=='-')flag=true;c=getchar();}
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
if(flag)x=-x;
}
ll n;
ll rev[maxn],f[maxn],g[maxn],fac[maxn],iv[maxn],w[22][maxn];
ll qp(ll x,ll y)
{
ll v=1;
while(y)
{
if(y&1) v=v*x%P;
x=x*x%P,y>>=1;
}
return v;
}
int calc(int n)
{
int lim=1;
while(lim<=n) lim<<=1;
for(int i=0;i<lim;++i)
rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
return lim;
}
void NTT(ll *a,int lim,int type)
{
for(int i=0;i<lim;++i)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(int len=1,t=1;len<lim;len<<=1,++t)
{
for(int i=0;i<lim;i+=len<<1)
{
for(int j=0;j<len;++j)
{
ll x=a[i+j],y=w[t][j]*a[i+j+len]%P;
a[i+j]=(x+y)%P,a[i+j+len]=(x-y+P)%P;
}
}
}
if(type==1) return;
ll inv=qp(lim,P-2);
for(int i=0;i<lim;++i) a[i]=a[i]*inv%P;
reverse(a+1,a+lim);
}
void Inv(int deg,ll *a,ll *b)
{
static ll t[maxn];
if(deg==1)
{
b[0]=qp(a[0],P-2);
return;
}
Inv((deg+1)>>1,a,b);
int lim=calc(deg<<1);
for(int i=0;i<deg;++i) t[i]=a[i];
for(int i=deg;i<lim;++i) t[i]=b[i]=0;
NTT(t,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;++i) b[i]=b[i]*(2-t[i]*b[i]%P+P)%P;
NTT(b,lim,-1);
for(int i=deg;i<lim;++i) b[i]=0;
}
void Ln(int deg,ll *a,ll *b)
{
static ll inva[maxn],dera[maxn];
Inv(deg,a,inva);
for(int i=0;i<deg-1;++i) dera[i]=a[i+1]*(i+1)%P;
dera[deg-1]=0;
int lim=calc(deg<<1);
for(int i=deg;i<lim;++i) dera[i]=inva[i]=0;
NTT(dera,lim,1),NTT(inva,lim,1);
for(int i=0;i<lim;++i) b[i]=dera[i]*inva[i]%P;
NTT(b,lim,-1);
for(int i=deg-1;i>=1;--i) b[i]=b[i-1]*qp(i,P-2)%P;
b[0]=0;
for(int i=deg;i<lim;++i) b[i]=0;
}
void Exp(int deg,ll *a,ll *b)
{
static ll t[maxn],lnb[maxn];
if(deg==1)
{
b[0]=1;
return;
}
Exp((deg+1)>>1,a,b),Ln(deg,b,lnb);
int lim=calc(deg<<1);
for(int i=0;i<deg;++i) t[i]=(a[i]-lnb[i]+P)%P;
for(int i=deg;i<lim;++i) t[i]=b[i]=0;
NTT(t,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;++i) b[i]=b[i]*(1+t[i])%P;
NTT(b,lim,-1);
for(int i=deg;i<lim;++i) b[i]=0;
}
void work(int deg,ll *f)
{
static ll h[maxn],t1[maxn],t2[maxn],t3[maxn];
if(deg==1)
{
f[0]=0;
return;
}
work((deg+1)>>1,f);
int lim=calc(deg<<1);
for(int i=0;i<deg;++i) t1[i]=f[i],t2[i]=f[i]*2;
NTT(t1,lim,1);
for(int i=0;i<lim;++i) t1[i]=(t1[i]*t1[i]-2*t1[i]+P)%P;
t2[0]-=2,Inv(deg,t2,t3),NTT(t3,lim,1);
for(int i=0;i<lim;++i) t3[i]=t1[i]*t3[i]%P,t1[i]=t1[i]+1;
NTT(t3,lim,-1),Exp(deg,t3,h),NTT(t1,lim,-1),Inv(deg,t1,t2);
for(int i=deg-1;i>=1;--i) h[i]=h[i-1];
h[0]=0;
for(int i=0;i<deg;++i) t1[i]=(2*h[i]-2*f[i]+P)%P;
t2[0]++,NTT(t2,lim,1),NTT(h,lim,1);
for(int i=0;i<lim;++i) t2[i]=t2[i]*h[i]%P;
NTT(t2,lim,-1),t2[0]-=2,Inv(deg,t2,t3),NTT(t1,lim,1),NTT(t3,lim,1);
for(int i=0;i<lim;++i) t1[i]=t1[i]*t3[i]%P;
NTT(t1,lim,-1);
for(int i=0;i<deg;++i) f[i]=(f[i]-t1[i]+P)%P;
for(int i=deg;i<lim;++i) f[i]=0;
}
void init()
{
fac[0]=fac[1]=iv[1]=1;
for(int i=2;i<=n;++i)
fac[i]=fac[i-1]*i%P,iv[i]=(P-P/i)*iv[P%i]%P;
for(int i=1;i<=19;++i) w[i][0]=1,w[i][1]=qp(G,(P-1)/(1<<i));
for(int i=1;i<=19;++i)
for(int j=2;j<1<<i;++j)
w[i][j]=w[i][j-1]*w[i][1]%P;
}
int main()
{
read(n),init(),work(n+1,f);
for(int i=0;i<=n;++i) f[i]=f[i]*iv[i]%P;
Exp(n+1,f,g),printf("%lld",g[n]*fac[n]%P);
return 0;
}