有标号仙人掌计数(图的计数)
有标号仙人掌计数(图的计数)
解题思路
设 \(F\) 表示仙人掌的生成函数，有
\[F = x \exp(F+\frac 12\sum_{n\ge2}F^i)
\]
具体来说就是枚举是环上边和非环上边，环上边要枚举环的大小，环上挂着仙人掌是一个排列但正反只能算一遍
\[F = x \exp(F+\frac 12\sum_{n\ge2}F^i)\\
= x \exp(F+\frac {F^2}{2-2F})\\
= x \exp(\frac {2F-F^2}{2-2F})
\]
然后直接牛顿迭代即可
\[\large G(F) = x\exp(\frac {2F-F^2}{2-2F})-F\\
\large F = F_0-\frac {G(F_o)}{G'(F_o)}\\
\large F = F_0-\frac {x\exp(\frac {2F_0-F_0^2}{2-2F_o})-F_0}{x\exp(\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 {x\exp(\frac {2F_0-F_0^2}{2-2F_o})-F_0}{x\exp(\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 {2x\exp(\frac {2F_0-F_0^2}{2-2F_o})-2F_0}{x\exp(\frac {2F_0-F_0^2}{2-2F_o})(1+\frac{1}{(1-F_0)^2})-2}\\
\]
为了简化式子，设 \(G_0 = x\exp(\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 = '\n')
{
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;
}
