【LOJ6569】仙人掌计数（生成函数）（多项式全家桶）（牛顿迭代）

传送门

• 题意：$n$ 个点带标号的仙人掌个数
• 求一个有根的个数，令其 $EGF$ 为 $f(x)$，那么一个根可以在多个环，或是只连出去一条边，对于环上的点，是以这个点为根的有序拼接（这样就不用管环上的顺序了），但是直接无序拼接的话正反会被算一遍，然后我们将每个环的贡献 $\exp$ 起来可以得到
  $$f(x)=x\exp (f(x)+\sum _{i\ge 2}\frac{f^i(x)}{2})=x\exp (\frac{2f(x)-f^2(x)}{2-2f(x)})$$
  解方程，牛顿迭代，假设求得 $f_0$ 使得 $H(f_0)\equiv 0(mod{x}^n)$
  则 $H(f)\equiv H(f_0)+H^{\prime }(f_0)(f-f_0)(mod{x}^{2n})$，即 $f=f_0-\frac{H(f_0)}{H^{\prime }(f_0)}$
  $$\begin{gathered} H^{\prime }(f)=(\frac{2f-f^2}{2-2f}{)}^{\prime }x\exp (\frac{2f-f^2}{2-2f})-1 \\ (\frac{2f-f^2}{2-2f}{)}^{\prime }=\frac{(2-2f)(2-2f)+2(2f-f^2)}{(2-2f{)}^2}=1+\frac{4f-2f^2}{(2-2f{)}^2} \end{gathered}$$
  然后就可以迭代了，实现得不是很精细所以应该有点慢

#include<bits/stdc++.h>
#define cs const
#define pb push_back
#define poly vector<int>
using namespace std;
cs int Mod = 998244353;
int add(int a, int b){ return a + b >= Mod ? a + b - Mod : a + b; }
int dec(int a, int b){ return a - b < 0 ? a - b + Mod : a - b; }
int mul(int a, int b){ return 1ll * a * b % Mod; }
int ksm(int a, int b){ int as=1; for(;b;b>>=1,a=mul(a,a)) if(b&1) as=mul(as,a); return as; }
void Add(int &a, int b){ a = add(a,b); }
void Mul(int &a, int b){ a = mul(a,b); }
void Dec(int &a, int b){ a = dec(a,b); }
cs int M = 1 << 20 | 5, C = 20;
cs int N = 131200;
int fac[M], ifac[M], iv[M];
void output(poly a){
	for(int i=0; i<a.size(); i++) cout<<a[i]<<" ";
	puts("");
}
void pre_work(int n){
	fac[0]=fac[1]=ifac[0]=ifac[1]=iv[0]=iv[1]=1;
	for(int i=2; i<=n; i++) iv[i]=mul(Mod-Mod/i,iv[Mod%i]);
	for(int i=2; i<=n; i++) fac[i]=mul(fac[i-1],i);
	for(int i=2; i<=n; i++) ifac[i]=mul(ifac[i-1],iv[i]);
}
poly w[C+1];
void NTT_init(){
	for(int i=1; i<=C; i++) w[i].resize(1<<(i-1));
	int wn=ksm(3,(Mod-1)/(1<<C)); w[C][0]=1;
	for(int i=1; i<(1<<(C-1)); i++) w[C][i]=mul(w[C][i-1],wn);
	for(int i=C-1;i;i--) for(int j=0; j<(1<<(i-1)); j++) w[i][j]=w[i+1][j<<1];
}
int up, bit; poly rev;
void init(int deg){
	up=1; bit=0; while(up<deg) up<<=1,++bit; rev.resize(up);
	for(int i=0; i<up; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void NTT(poly &a, int typ=1){
	for(int i=0; i<up; i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=1,l=1;i<up;i<<=1,++l)
	for(int j=0; j<up; j+=(i<<1))
	for(int k=0; k<i; k++){
		int x=a[k+j], y=mul(w[l][k],a[k+j+i]);
		a[k+j]=add(x,y); a[k+j+i]=dec(x,y);
	}
	if(typ==-1){
		reverse(a.begin()+1,a.end());
		for(int i=0; i<up; i++) Mul(a[i],iv[up]);
	}
}
poly operator * (poly a, poly b){
	int deg=a.size()+b.size()-1; init(deg); 
	a.resize(up); b.resize(up); NTT(a); NTT(b);
	for(int i=0; i<up; i++) Mul(a[i],b[i]); 
	NTT(a,-1); a.resize(deg); return a;
}
poly inv(poly a, int deg){
	poly b(1,ksm(a[0],Mod-2)),c;
	for(int lim=4; (lim>>2)<deg; lim<<=1){
		c.resize(lim>>1); init(lim);
		for(int i=0; i<(lim>>1); i++) c[i]=i<(int)a.size()?a[i]:0;
		c.resize(up); b.resize(up); NTT(b); NTT(c);
		for(int i=0; i<up; i++) Mul(b[i],dec(2,mul(b[i],c[i])));
		NTT(b,-1); b.resize(lim>>1); 
	} b.resize(deg); return b;
}
poly deriv(poly a){
	for(int i=0; i+1<(int)a.size(); i++) a[i]=mul(i+1,a[i+1]);
	a.pop_back(); return a;
}
poly integ(poly a){
	a.pb(0);
	for(int i=a.size()-1;i;i--) a[i]=mul(iv[i],a[i-1]);
	a[0]=0; return a;
}
poly ln(poly a, int deg){
	a = integ(inv(a,deg)*deriv(a)); 
	a.resize(deg); return a;
}
poly Exp(poly a, int deg){
	poly b(1,1), c;
	for(int lim=2; (lim>>1)<deg; lim<<=1){
		c=ln(b,lim); Dec(c[0],1);
		for(int i=0; i<lim; i++) c[i]=dec(i<(int)a.size()?a[i]:0,c[i]);
		b=b*c; b.resize(lim); 
	} b.resize(deg); return b;
}
poly operator * (int coe, poly a){
	for(int i=0; i<(int)a.size(); i++)
	Mul(a[i],coe); return a;
}
poly operator + (poly a, poly b){
	int deg = max(a.size(),b.size()); a.resize(deg); b.resize(deg);
	for(int i=0; i<deg; i++) Add(a[i],b[i]); return a;
}
poly operator - (poly a, poly b){
	int deg = max(a.size(),b.size()); a.resize(deg); b.resize(deg);
	for(int i=0; i<deg; i++) Dec(a[i],b[i]); return a;
}
poly Mulx(poly a){ a.insert(a.begin(),0); return a; }
poly Newton(int n){
	if(n==1) return poly(1,1);
	poly f0 = Newton((n+1)>>1);
	poly t = Mulx(Exp((2*f0-f0*f0)*inv(poly(1,2)-2*f0,n),n));t.resize(n);
	poly f = f0 - ((2*t-2*f0)*inv((poly(1,1)+inv((f0-poly(1,1))*(f0-poly(1,1)),n))*t-poly(1,2),n)); 	
	f.resize(n); return f;
}
int main(){
	pre_work(1<<C); NTT_init();
	poly f = Newton(N); 
	int T, n; scanf("%d",&T);
	while(T--){
		scanf("%d",&n); 
		cout << mul(f[n],fac[n-1]) << '\n';
	} return 0;
}