P5434 有标号荒漠计数

P5434 有标号荒漠计数

数学、组合计数、生成函数

---

其实森林也是沙漠

挺难的一道题，并非独立切。

设 \(f_n\) 为 \(n\) 个点仙人掌数量，\(F(x)\) 为 \(\{f_n\}\) 的 EGF，答案即为 \([\frac{x^n}{n!}]\exp(F(x))\)。

图论计数，一般考虑设函数方程或者将一个组合对象拆分为多个小组合对象，这里考虑设函数方程。

发现 \(F(x)\) 其实不能转移，考虑设

\[G(x)=\sum\limits_{i=0}^{+\infty}\frac{f_i}{(i-1)!} \]

其组合意义为有根仙人掌的生成函数。

对一个仙人掌，考虑任取一个点 \(u\) 作为根。与 \(u\) 相连的桥，其贡献为 \(G(x)\)；与 \(u\) 相连的环，其贡献为 \(\frac{1}{2}\sum\limits_{i=2}^{+\infty}G^i(x)=\frac{G^2(x)}{2(1-G(x))}\) 这里除以 \(2\)、不除阶乘是因为与 \(u\) 相连的环是可翻转的排列。
将桥和环任意组合（无序无数量限制，即 \(\exp\)）可以得到关于 \(G(x)\) 的函数方程。

\[G(x)=x\exp\left(G(x)+\frac{G^2(x)}{2\left(1-G(x)\right)}\right) \]

牛顿迭代即可，时间复杂度 \(\Theta(n\log n)\)。

---

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
typedef pair<ll,ll> pll;
const int MOD=998244353,G=3,INV2=(MOD+1)>>1;
void add(int &x,int y){
	x+=y;
	if(x>=MOD) x-=MOD;
}
int qpow(int a,int b){
	int mul=1;
	while(b){
		if(b&1) mul=(ll)mul*a%MOD;
		a=(ll)a*a%MOD;
		b>>=1;
	}
	return mul;
}
const int N=800005;
int fact[N],invfact[N],inv[N],stp[N],invstp[N];
void init(){
	fact[0]=1;
	for(int i=1;i<N;i++) fact[i]=(ll)i*fact[i-1]%MOD;
	invfact[N-1]=qpow(fact[N-1],MOD-2);
	for(int i=N-1;i>0;i--) invfact[i-1]=(ll)i*invfact[i]%MOD;
	for(int i=1;i<N;i++) inv[i]=(ll)invfact[i]*fact[i-1]%MOD;
	for(int w=1;w<N;w<<=1) stp[w]=qpow(G,(MOD-1)/w),invstp[w]=qpow(stp[w],MOD-2);
}
namespace Poly{
	void ntt(vector<int> &f,bool flag=false){
		int n=f.size();
		vector<int> rev(n);
		for(int i=1;i<n;i++){
			rev[i]=rev[i>>1]>>1;
			if(i&1) rev[i]|=(n>>1);
		}
		for(int i=1;i<n;i++) if(rev[i]<i) swap(f[rev[i]],f[i]);
		for(int w=1;w<n;w<<=1){
			int step=stp[w<<1];
			if(flag) step=invstp[w<<1];
			for(int i=0;i<n;i+=(w<<1)){
				int cur=1;
				for(int j=i;j<i+w;j++,cur=(ll)cur*step%MOD){
					int a=f[j],b=(ll)cur*f[j+w]%MOD;
					f[j]=(a+b)%MOD,f[j+w]=(a-b+MOD)%MOD;
				}
			}
		}
		if(flag){
			int inv=::inv[n];
			for(int i=0;i<n;i++) f[i]=(ll)inv*f[i]%MOD;
		}
	}
	vector<int> inv(vector<int> f,int n=-1){
		if(n==1) return {qpow(f[0],MOD-2)};
		if(n==-1) n=f.size();
		else f.resize(n);
		vector<int> h=inv(f,n>>1);
		f.resize(n<<1),h.resize(n<<1);
		ntt(f),ntt(h);
		for(int i=0;i<(n<<1);i++) f[i]=(2*h[i]%MOD-(ll)h[i]*h[i]%MOD*f[i]%MOD+MOD)%MOD;
		ntt(f,1);
		f.resize(n);
		return f;
	}
	vector<int> ln(vector<int> f){
		int n=f.size();
		vector<int> h=inv(f);
		for(int i=0;i+1<n;i++) f[i]=(ll)(i+1)*f[i+1]%MOD;
		f[n-1]=0;
		f.resize(n<<1),h.resize(n<<1);
		ntt(f),ntt(h);
		for(int i=0;i<(n<<1);i++) f[i]=(ll)f[i]*h[i]%MOD;
		ntt(f,1);
		for(int i=n-1;i>0;i--) f[i]=(ll)::inv[i]*f[i-1]%MOD;
		f[0]=0;
		f.resize(n);
		return f;
	}
	vector<int> exp(vector<int> f,int n=-1){
		if(n==1) return {1};
		if(n==-1) n=f.size();
		else f.resize(n);
		vector<int> h=exp(f,n>>1);
		h.resize(n);
		vector<int> g=ln(h);
		f.resize(n<<1),g.resize(n<<1),h.resize(n<<1);
		ntt(f),ntt(g),ntt(h);
		for(int i=0;i<(n<<1);i++) f[i]=(ll)h[i]*(1+MOD-g[i]+f[i])%MOD;
		ntt(f,1);
		f.resize(n);
		return f;
	}
	vector<int> mul(vector<int> f,vector<int> g,int n){
		f.resize(n<<1),g.resize(n<<1);
		ntt(f),ntt(g);
		for(int i=0;i<(n<<1);i++) f[i]=(ll)f[i]*g[i]%MOD;
		ntt(f,1);
		f.resize(n);
		return f;
	}
	vector<int> calc(int n){
		if(n==1) return {0};
		vector<int> f=calc(n>>1);
		f.resize(n);
		vector<int> g=f;
		for(int i=0;i<n;i++) g[i]=(ll)(MOD-2)*g[i]%MOD;
		add(g[0],2);
		g=mul(mul(f,f,n),inv(g),n);
		for(int i=0;i<n;i++) add(g[i],f[i]);
		vector<int> A=exp(g),tmp=f;
		for(int i=n-1;i>0;i--) A[i]=A[i-1];
		A[0]=0;
		A.resize(n<<2),f.resize(n<<2);
		ntt(A),ntt(f);
		vector<int> B(n<<2),C(n<<2);
		for(int i=0;i<(n<<2);i++){
			int tmp=(ll)(1+MOD-f[i])*(1+MOD-f[i])%MOD;
			B[i]=(ll)(A[i]-f[i]+MOD)*tmp%MOD;
			C[i]=((ll)A[i]*(tmp+(ll)f[i]*(1+MOD-f[i])%MOD+(ll)f[i]*f[i]%MOD*INV2%MOD)%MOD-tmp+MOD)%MOD;
		}
		ntt(B,1),ntt(C,1);
		B.resize(n),C.resize(n);
		B=mul(B,inv(C),n);
		B.resize(n);
		for(int i=0;i<n;i++) add(tmp[i],MOD-B[i]);
		return tmp;
	}
}
int n;
int main(){
	init();
	scanf("%d",&n);
	int len;
	for(len=1;len<=n;len<<=1);
	vector<int> f=Poly::calc(len);
	for(int i=1;i<len;i++) f[i]=(ll)inv[i]*f[i]%MOD;
	f=Poly::exp(f);
	printf("%lld\n",(ll)fact[n]*f[n]%MOD);
	return 0;
}