有标号仙人掌计数(图的计数)

有标号仙人掌计数(图的计数)

解题思路

设 (F) 表示仙人掌的生成函数，有

[F = x exp(F+frac 12sum_{nge2}F^i) ]

具体来说就是枚举是环上边和非环上边，环上边要枚举环的大小，环上挂着仙人掌是一个排列但正反只能算一遍

[F = x exp(F+frac 12sum_{nge2}F^i)\ = x exp(F+frac {F^2}{2-2F})\ = x exp(frac {2F-F^2}{2-2F}) ]

然后直接牛顿迭代即可

[large G(F) = xexp(frac {2F-F^2}{2-2F})-F\ large F = F_0-frac {G(F_o)}{G'(F_o)}\ large F = F_0-frac {xexp(frac {2F_0-F_0^2}{2-2F_o})-F_0}{xexp(frac {2F_0-F_0^2}{2-2F_o})(frac {2F_0-F_0^2}{2-2F_o})'-1}\ ]

[(frac {2x-x^2}{2-2x})' = x + frac {x^2}{2-2x}=1+frac {2x(2-2x)+2x^2}{(2-2x)^2}\ =frac {4-4x+2x^2}{(2-2x)^2}=frac {x^2-2x+2}{2(1-x)^2} ]

[large F = F_0-frac {xexp(frac {2F_0-F_0^2}{2-2F_o})-F_0}{xexp(frac {2F_0-F_0^2}{2-2F_o})frac{F_0^2-2F_0+2}{2(1-F_0)^2}-1}\ large F = F_0-frac {2xexp(frac {2F_0-F_0^2}{2-2F_o})-2F_0}{xexp(frac {2F_0-F_0^2}{2-2F_o})(1+frac{1}{(1-F_0)^2})-2}\ ]

为了简化式子，设 (G_0 = xexp(frac {2F_0-F_0^2}{2-2F_o}))，有

[large F= F_0-frac{2G_0-2F_0}{(1+frac{1}{(1-F_0)^2})G_0-2} ]

可以暴力解了

#pragma GCC optimize(2)
#include <queue>
#include <vector>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define MP make_pair
#define ll long long
#define fi first
#define se second
using namespace std;

template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}

//template <typename T>
//void read(T &x) {
//    x = 0; bool f = 0;
//    char c = getchar();
//    for (;!isdigit(c);c=getchar()) if (c=='-') f=1;
//    for (;isdigit(c);c=getchar()) x=x*10+(c^48);
//    if (f) x=-x;
//}

template<typename F>
inline void write(F x, char ed = '
')
{
	static short st[30];short tp=0;
	if(x<0) putchar('-'),x=-x;
	do st[++tp]=x%10,x/=10; while(x);
	while(tp) putchar('0'|st[tp--]);
	putchar(ed);
}

template <typename T>
inline void Mx(T &x, T y) { x < y && (x = y); }

template <typename T>
inline void Mn(T &x, T y) { x > y && (x = y); }

const int N = 633400;
const int P = 998244353;
int r[N], I[N*20], lim, n;
ll A[N], B[N], C[N], D[N], g[20][N];
ll inv[N], fac[N], q;

ll fpw(ll x, ll mi) {
	ll res = 1;
	for (; mi; mi >>= 1, x = x * x % P)
		if (mi & 1) res = res * x % P;
	return res;
}

void init(int n) {
	inv[0] = inv[1] = fac[0] = fac[1] = 1;
	for (int i = 2;i <= n; i++)
		inv[i] = inv[P % i] * (P - P / i) % P, fac[i] = fac[i-1] * i % P;
	for (int i = 0;i <= 19; i++) I[1 << i] = i;
    for (int i = 0;i <= 19; ++i) {
        ll *G = g[i]; G[0] = 1;
        const int gi = G[1] = fpw(3, (P - 1) / (1 << (i+1)));
        for(int j = 2;j < 1 << i; ++j) G[j] = G[j-1] * gi % P;
    }
}


void NTT(ll*a,int f){
    for(int i=1;i<lim;++i)if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=0;i<I[lim];++i){
        const ll*G=g[i],c=1<<i;
        for(int j=0;j<lim;j+=c<<1)
        for(int k=0;k<c;++k){
            const int x=a[j+k],y=a[j+k+c]*G[k]%P;
            (a[j+k]+=y)%=P,a[j+k+c]=(x-y+P)%P;
        }
    }
    if(f==-1){
    	ll inv = fpw(lim, P - 2);
        for(int i=0;i<lim;++i)a[i]=a[i]*inv%P;
        std::reverse(a+1,a+lim);
    }
}


/*

F(G) - A = 0;
F(G0) - A = 0;
0 = T(G0)+T(G0)'(G-G0)
G = G0-T(G0)/T(G0)'

1/G0 - A = 0;

G = G0 + (1/G0-A)G0^2
G = G0 + G0 - G0^2A
G = G0(2 - A*G0)

*/

void Inv(ll *a, ll *b, int deg) {
	b[0] = fpw(a[0], P - 2); int len;
	for (len = 2;len < (deg << 1); len <<= 1) {
		lim = len << 1;
		for (int i = 1;i < lim; i++)
			r[i] = (r[i>>1]>>1) | ((i & 1) ? len : 0);
		for (int i = 0;i < len; i++) A[i] = a[i], B[i] = b[i];
		NTT(A, 1), NTT(B, 1);
		for (int i = 0;i < lim; i++)
			b[i] = (2 + (P - A[i]) * B[i]) % P * B[i] % P;
		NTT(b, -1);
		for (int i = len;i < lim; i++) b[i] = 0;
	}
    for (int i = 0;i < len; i++) A[i] = B[i] = 0;
}

void chick(ll *a, int deg) {
	for (int i = deg - 1;i >= 1; i--)
		a[i] = a[i-1] * inv[i] % P;
	a[0] = 0;
}

void door(ll *a, int deg) {
	for (int i = 1;i < deg; i++)
		a[i-1] = a[i] * i % P;
	a[deg - 1] = 0;
}

void Ln(ll *a, int deg) {
	Inv(a, C, deg), door(a, deg);
	NTT(a, 1), NTT(C, 1);
	for (int i = 0;i < lim; i++)
		a[i] = a[i] * C[i] % P;
	NTT(a, -1), chick(a, deg);
}

/*

F(G) - A = 0;
F(G0) - A = 0;
0 = T(G0)+T(G0)'(G-G0)
G = G0-T(G0)/T(G0)'

1/G0 - A = 0;

G = G0 + (LnG0-A)G0
G = -(-1+LnG0-A)G0

*/

void Exp(ll *a, ll *b, int deg) {
	b[0] = 1; int len;
	for (len = 2;len < (deg << 1); len <<= 1) {
		for (int i = 0;i < len; i++) D[i] = b[i];
		Ln(D, len); D[0] = (a[0] + 1 - D[0] + P) % P;
		for (int i = 1;i < len; i++)
			D[i] = (a[i] - D[i] + P) % P;
		NTT(D, 1), NTT(b, 1);
		for (int i = 0;i < lim; i++)
			b[i] = b[i] * D[i] % P;
		NTT(b, -1);
		for (int i = len;i < lim; i++) b[i] = 0;
	}
	for (int i = 0;i < len; i++) D[i] = 0;
}

ll T[N], F[N], F0[N], F1[N], F2[N], G[N], E[N], IF[N];
void newton(ll *f, int deg) {
//	f[1] = 1;
	int len; ll inv2 = (P + 1) >> 1;
	for (len = 2;len < (deg << 1); len <<= 1) {
		for (int i = 0;i < len; i++) F2[i] = f[i], T[i] = P - f[i], G[i] = 0, F0[i] = f[i];
		T[0]++, Inv(T, IF, len);
		NTT(IF, 1), NTT(F2, 1), NTT(F0, 1);
		for (int i = 0;i < lim; i++) {
			F2[i] = (2 * F0[i] - F2[i] * F2[i]) % P * IF[i] % P * inv2 % P;
			IF[i] = IF[i] * IF[i] % P;
		}
		NTT(F2, -1), Exp(F2, G, len), NTT(IF, -1), IF[0]++;
		for (int i = len - 1;i >= 1; i--) G[i] = G[i-1]; G[0] = 0;
		for (int i = len;i < lim; i++) IF[i] = 0;
		NTT(IF, 1), NTT(G, 1);
		for (int i = 0;i < lim; i++) IF[i] = IF[i] * G[i] % P;
		NTT(IF, -1), IF[0] = (IF[0] + P - 2) % P;
		for (int i = 0;i < lim; i++) G[i] = (2 * G[i] - 2 * F0[i] + 2 * P) % P;
		Inv(IF, F0, len), NTT(F0, 1);
		for (int i = 0;i < lim; i++) F0[i] = F0[i] * G[i] % P;
		NTT(F0, -1);
		for (int i = 0;i < len; i++) f[i] = (f[i] - F0[i] + P) % P;
//		for (int i = len;i < lim; i++) f[i] = 0;
//		cout << "OK" << endl;
	}
}

int main() {
	freopen ("hs.in","r",stdin);
	freopen ("hs.out","w",stdout); 
	read(n);
//	n = 131073;
	init(n + 1), newton(F, n + 1);
//	for (read(q); q; q--) 
//		read(n), write(F[n] * fac[n-1] % P);
	for (int i = 0;i <= n; i++) F[i] = F[i] * inv[i] % P;
	Exp(F, F1, n + 1);
	write(F1[n] * fac[n] % P);
	return 0;
}