有标号荒漠计数

有标号荒漠计数

考虑随意选择一个点为根，则仙人掌的\(\text{EGF}\)考虑用以下方式递归生成

令树边为二元环，则一个点周围的点都是都是与它直接相连的环

断开这个点，对于周围断开的环，环上每个点下面认为是一个仙人掌，设某个环断开之后的大小为\(c\)

当\(c=1\)时，不需要考虑排列重复，即为\(F(x)\)

当\(c>1\)时，考虑环正反排列，即为\(\cfrac{F^c(x)}{2}\)

那么就容易得到\(\displaystyle F(x)=x \cdot \text{exp}(F(x)+\sum _{i\ge 2}\frac{F^i(x)}{2})\)

变一下就是\(\displaystyle F=x\cdot \text{exp}(\frac{F^2}{2-2F}+F)=x\cdot \text{exp}(\frac{2F-F^2}{2-2F})\)

是的，我们要解这个方方方方方方程。。。牛顿迭代代代代代代代

\(\displaystyle f(F(x))=x\cdot \text{exp}(\frac{2F-F^2}{2-2F})-F=0\)

\(\displaystyle f(z)=x\cdot \text{exp}(\frac{2z-z^2}{2-2z})-z\)

\(\displaystyle f'(z)=x\cdot \text{exp}(\frac{2z-z^2}{2-2z})(1+\frac{2z-z^2}{2z^2-4z+2})-1\)

\(\displaystyle =x\cdot \text{exp}(\frac{2z-z^2}{2-2z})(\frac{1}{2}+\frac{1}{2z^2-4z+2})-1\)

设上一层的迭代结果为\(G(x)\)，带入牛顿迭代结论\(\displaystyle F(x)=G(x)-\frac{f(G)}{f'(G)}\)

设\(\displaystyle H=x\cdot \text{exp}(\frac{2G-G^2}{2-2G})\)，那么得到Luogu题解里\(\text{N}\color{red}{\text{aCl_Fish}}\)一样的式子(还要没有推错)

\(\displaystyle F=G-\frac{2H-2G}{H(1+\frac{1}{(1-G)^2})-2}\)

最后还要变成无根，除掉\(n\)即可

仙人掌转荒漠您只需要一个\(\text{exp}\)就好了

const int L=18,N=1<<L|10,P=998244353;

ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

int I[N],J[N];
int rev[N],w[N];
void Init(){
	w[1<<(L-1)]=1;
	int t=qpow(3,(P-1)>>L);
	rep(i,(1<<(L-1))+1,1<<L) w[i]=1ll*w[i-1]*t%P;
	drep(i,(1<<(L-1))-1,1) w[i]=w[i<<1];
	rep(i,J[0]=1,N-1) J[i]=1ll*J[i-1]*i%P;
	I[N-1]=qpow(J[N-1]);
	drep(i,N-1,1) I[i-1]=1ll*I[i]*i%P;
}
int Init(int n){
	int R=1,c=-1;
	while(R<=n) R<<=1,c++;
	rep(i,0,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
	return R;
}
void NTT(int n,V &A,int f) {
	static ull a[N];
	if((int)A.size()<n) A.resize(n);
	rep(i,0,n-1) a[rev[i]]=A[i];
	for(int i=1;i<n;i<<=1) {
		int *e=w+i;
		for(int l=0;l<n;l+=i*2) {
			for(int j=l;j<l+i;++j) {
				int t=a[j+i]*e[j-l]%P;
				a[j+i]=a[j]+P-t;
				a[j]+=t;
			}
		}
	}
	rep(i,0,n-1) A[i]=a[i]%P,Mod2(A[i]);
	if(f==-1) {
		reverse(A.begin()+1,A.end());
		ll base=1ll*I[n]*J[n-1]%P;
		rep(i,0,n-1) A[i]=A[i]*base%P;
	}
}

V operator + (V a,const V &b) {
	if(a.size()<b.size()) a.resize(b.size());
	rep(i,0,b.size()-1) a[i]+=b[i],Mod1(a[i]);
	return a;
}
V operator - (V a,const V &b) {
	if(a.size()<b.size()) a.resize(b.size());
	rep(i,0,b.size()-1) a[i]-=b[i],Mod2(a[i]);
	return a;
}
V operator * (V a,V b) {
	int n=a.size()-1,m=b.size()-1;
	int R=Init(n+m);
	NTT(R,a,1),NTT(R,b,1);
	rep(i,0,R-1) a[i]=1ll*a[i]*b[i]%P;
	NTT(R,a,-1),a.resize(n+m+1);
	return a;
}
V operator * (V a,const int &x) {
	for(int &i:a) i=1ll*i*x%P;
	return a;
}
V operator * (const int &x,V a) { return a*x; }

void println(const V &a){
	for(int i:a) printf("%d ",i);
	puts("");
}
V read(int n){
	V A(n);
	rep(i,0,n-1) A[i]=rd();
	return A;
}
V operator ~ (V a) {
	int n=a.size(),m=(n+1)>>1;
	if(n==1) return {(int)qpow(a[0])};
	V b=a; b.resize(m),b=~b;
	int R=Init(n*2);
	NTT(R,a,1),NTT(R,b,1);
	rep(i,0,R-1) a[i]=(P+2-1ll*a[i]*b[i]%P)*b[i]%P;
	NTT(R,a,-1),a.resize(n);
	return a;
}
V Deriv(V a) {
	rep(i,0,a.size()-2) a[i]=1ll*(i+1)*a[i+1]%P;
	a.pop_back();
	return a;
}
V Integ(V a){
	a.pb(0);
	drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*J[i-1]%P*I[i]%P;
	a[0]=0;
	return a;
}
V Ln(V a){
	int n=a.size();
	a=Deriv(a)*~a;
	return a.resize(n-1),Integ(a);
}
V Exp(V a) {
	if(a.size()==1) return assert(a[0]==0),V{1};
	int n=a.size();
	V b=a; b.resize((n+1)/2),b=Exp(b),b.resize(n);
	a=a-Ln(b),a[0]++;
	a=a*b,a.resize(n);
	return a;
}

V operator << (V a,int x) {
	a.resize(a.size()+x);
	drep(i,a.size()-1,x) a[i]=a[i-x];
	rep(i,0,x-1) a[i]=0;
	return a;
}
V operator >> (V a,int x) {
	if((int)a.size()<=x) return V{};
	rep(i,x,a.size()-1) a[i-x]=a[i];
	a.resize(a.size()-x);
	return a;
}

V Newton(int n){
	if(n==1) return V{0};
	if(n==2) return V{0,1};
	V G=Newton((n+1)/2); G.resize(n);
	V IG=~(V{1}-G);
	V H=(2*G-G*G); H.resize(n),H=H*IG*((P+1)/2),H.resize(n),H=Exp(H)<<1;
	V F=IG*IG; F.resize(n),F[0]++;
	F=H*F,F.resize(n),F[0]-=2,Mod2(F[0]);
	F=G-2*(H-G)*~F;
	return F.resize(n),F;
}

int main(){
	int n=rd()+1; Init();
	V F=Newton(n); 
	rep(i,1,F.size()-1) F[i]=1ll*F[i]*I[i]%P*J[i-1]%P;
	F=Exp(F);
	printf("%d\n",int(1ll*F.back()*J[n-1]%P));
}