闲话 22.12.9

闲话

我说我今天主要刷的是网络流题你信嘛（
主要是在上网课 然后老师讲课：“一说荒漠就会想到仙人掌”
就在群里发这题
浅看了一下发现不会 然后写了

是……辣个男人！他又来了（

有标号荒漠计数

P5434

我们定义一张图是荒漠，当且仅当其每个极大联通分量都是仙人掌。
给定 \(n\)，求有多少 \(n\) 个点的不同的荒漠。结果对 \(998244353\) 取模。

\(n \le 10^5\)。

首先要计数有标号荒漠，可以先计数有标号仙人掌。最后作 \(\exp\) 就能得到答案。
注意答案是无根的，但是仙人掌显然是有根好统计一些。因此我们记 \(G\) 为有标号有根仙人掌计数的生成函数。考虑将一系列根连成环接在指定的一个根上，并由环没有方向，有

\[\begin{aligned} G & = x\exp\left(G + \frac 12 \sum_{i \ge 2} G^i\right) \\ & = x\exp\left(G + \frac 12 \left(\frac{1}{1 - G} - (1 + G) \right) \right) \\ & = x\exp\left(G + \frac{G^2}{2(1 - G)} \right) \\ & = x\exp\left( \frac{2G - G^2}{2(1 - G)} \right) \end{aligned}\]

无根考虑答案即为 \([x^n]\exp\left(xG'(x)\right)\)。
由于 \(\exp\) 里面东西太多，我们考虑记

\[h(x) = \frac{2x - x^2}{2(1-x)} \]

设 \(F\) 为答案计数的生成函数，有 \(F = x\left(\exp(h(x))\right)'\)。
到这一步的话自然想到 Lagrange 反演。设 \(F = g(f(x))\)，有 \(G = x\left(g(f(x))\right)' = x\ f'(x)\ g'(f(x))\)。现在只需要把另一部分化成同一形式就行了。

\[\begin{aligned} G' & = \exp(h(G)) + x\left(\exp(h(G))\right)' \\ & = \exp(h(G)) + x\ G'\ h'(G)\ \exp(h(G)))' \\ & = \exp(h(G))(1 + x\ G'\ h'(G)) \\ xG'& = x\exp(h(G))(1 + x\ G'\ h'(G)) \\ & = G(1 + x\ G'\ h'(G)) \\ & = G + x\ G\ G'\ h'(G) \\ G & = xG' - x\ G\ G'\ h'(G) \\ & = xG'(1 - G \ h'(G) ) \end{aligned}\]

自然导出了 \(f = G, g'(x) = 1 - xh'(x)\)。记 \(H(x) = \int h(x) dx\)，有

\[\begin{aligned} g(x) & = \int1 - xh'(x) \text dx \\ & = x - \int xh'(x) \text dx \\ & = x - \left(x h(x) - \int h(x) \text dx\right) \\ & = x - x h(x) + H(x) \end{aligned}\]

这样我们需要的所有拼图都齐全了。拉格朗日反演得到

\[[x^n] e^{F(x)} = [x^n] e^{g(f(x))} = \frac 1n [x^{n-1}] g'(x) e^{g(x)} \left(\frac {x}{f^{⟨-1⟩}(x)}\right) \]

于是可以分别讨论三个（两个？）部分如何做。

\(g'(x)\)

显然

下面需要 \(h'(x)\)。根据 Wolfram|Alpha，

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

于是

\[\begin{aligned} g'(x) = 1 - xh'(x) = 1 - \frac{x^3 - 2x^2 - 2x}{2(1 - x)^2} = \frac{-\frac 12 x^3 + 2x^2 - 3x + 1}{(1 - x)^2} \end{aligned}\]

\(x / f^{⟨-1⟩}(x)\)

观察 \(f^{⟨-1⟩}\) 是个啥。其一是 \(f^{⟨-1⟩} (f(x)) = x\)，其二是 \(f = x\exp\left( \frac{2f- f^2}{2(1 - f)} \right)\)。我们发现，假设 \(f^{⟨-1⟩}\) 中是连乘形式，且其中有 \(x\) 项，将 \(f(x)\) 带入后，只需要将后面的 \(\exp\) 消除即可。
不妨考虑 \(f^{⟨-1⟩}(x) = x\exp\left(\frac{x^2 - 2x}{2(1 - x)}\right)\)，容易发现此函数满足需要。

因此我们能得到 \(x / f^{⟨-1⟩}(x) = \exp\left(\frac{2x - x^2}{2(1 - x)}\right)\)。
因为上面还有一个 \(n\) 次方，该部分的贡献是 \(\exp\left(\frac{n(2x - x^2)}{2(1 - x)}\right)\) 。

\(e^{g(x)}\)

下面需要 \(H(x)\)。根据 Wolfram|Alpha，

\[H(x) = \frac 14\left((x-1)^2 - 2\ln(1- x)\right) \]

于是

\[\begin{aligned} g(x) & = x - x h(x) + H(x) \\ & = x + \left(\frac{x - 1}2\right)^2 - \frac 12 \ln (1 - x) - \frac{x(2x - x^2)}{2(1 - x)} \\ & = \frac{2x + x^2} 4 - \frac 12 \ln (1 - x) - \frac{x(2x - x^2)}{2(1 - x)} \\ e^{g(x)} & = \exp\left( \frac{2x + x^2} 4 - \frac 12 \ln (1 - x) - \frac{x(2x - x^2)}{2(1 - x)} \right) \\ [后两项] & = \exp\left( \frac{2x + x^2} 4 - \frac 12 \ln (1 - x) - \frac{x(2x - x^2)}{2(1 - x)} + \frac{n(2x - x^2)}{2(1 - x)} \right) \\ [后两项] & = (1 - x)^{- 1/2} \exp\left( \frac{2x + x^2} 4 - \frac{x(2x - x^2)}{2(1 - x)} + \frac{n(2x - x^2)}{2(1 - x)} \right) \\ [后两项] & = (1 - x)^{- 1/2} \exp\left( \frac{2x + x^2} 4 + (n - x)\frac{(2x - x^2)}{2(1 - x)} \right) \\ [整个的] & = \frac{-\frac 12 x^3 + 2x^2 - 3x + 1}{(1 - x)^{2/5}} \exp\left( \frac{2x + x^2} 4 + (n - x)\frac{(2x - x^2)}{2(1 - x)} \right) \end{aligned}\]

根据化简过程中不很简的部分，这玩意微分有限，存在线性递推，因此可以 \(O(n) / O(\sqrt n \log n)\) 地算。

具体地，我们首先考察 \(a(x) = \exp\left( \frac{2x + x^2} 4 + (n - x)\frac{(2x - x^2)}{2(1 - x)} \right)\)。求导后化简可以得到

\[2(1 - x)^2 a'(x) = \left((2n + 1) - (2n + 5) x + (n + 4)x^2 - x^3\right) a(x) \]

过程 dirty，不细说了（

随后考察 \(b(x) = (1 - x)^{-5/2} a(x)\)。同理有

\[2(1 - x)^2 b'(x) = \left((2n + 6) - (2n + 10) x + (n + 4)x^2 - x^3\right) b(x) \]

这自然导出了递推方案。答案即为

\[[x^n / n!]\ b(x) \left(-\frac 12 x^3 + 2x^2 - 3x + 1\right) \]

初值抄鰰的（
代码全屏好看（

code

#include <bits/stdc++.h>
using namespace std; using pii = pair<int,int>; using vi = vector<int>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
template <typename T> T rand(T l, T r) { return uniform_int_distribution<T>(l, r)(rnd); }
template <typename T1, typename T2> T1 min(T1 a, T2 b) { return a < b ? a : b; }
template <typename T1, typename T2> T1 max(T1 a, T2 b) { return a > b ? a : b; }
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define ufile(x) 
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define Aster(i, s) for (int i = head[s], v; i; i = e[i].next)
#define all(s) s.begin(), s.end()
#define eb emplace_back
#define pb pop_back
#define em emplace
const int N = 1e5 + 10, mod = 998244353, ans[] = {1, 1, 2, 8, 57};
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
int n, inv[N], fac;

signed main() {
    ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cin >> n;
    if (n <= 4) cout << ans[n], exit(0);
    inv[0] = inv[1] = fac = 1;
    rep(i,2,n) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod, fac = 1ll * fac * i % mod;
    vector<int> coef{ 0, 2 * n + 6, mod - (2 * n + 10), n + 4, mod - 1 },
    g{ 1, n + 3, ((1ll * (n + 7) * n >> 1) + 5) % mod, ((1ll * (n + 12) * n % mod + 42) * n % mod + 38) * inv[6] % mod };
    g.resize(n + 1);
    rep(i,4,n) {
        g[i] = (1ll * coef[1] * g[i - 1] % mod + 1ll * coef[2] * g[i - 2] % mod +
                1ll * coef[3] * g[i - 3] % mod + 1ll * coef[4] * g[i - 4] % mod) * inv[2] % mod;
        g[i] = (g[i] + 1ll * ((i - 1) << 1) * g[i - 1] % mod - 1ll * (i - 2) * g[i - 2] % mod + mod) % mod * inv[i] % mod;
    }
    fac = 1ll * fac * (g[n] - 3ll * g[n - 1] % mod + 2ll * g[n - 2] - 1ll * inv[2] * g[n - 3] % mod + 2 * mod) % mod;
    cout << (fac + mod) % mod << endl;
}

upd：使用 ODE+\(O(\sqrt n \log n)\) 自动机，可以让计算机帮我们在 \(O(\sqrt n \log n)\) 复杂度内完成。Submission

code

signed main() {
	int n; cin >> n;
	ODE B(2);
	B[0].redegree(3); B[1].redegree(2);
	B[0][3] = 1, B[0][2] = mod - n - 4, B[0][1] = 2 * n + 10, B[0][0] = mod - 2 * n - 6;
	B[1][2] = 2, B[1][1] = mod - 4, B[1][0] = 2;
	poly f(4); f[3] = mod - ginv(2), f[2] = 2, f[1] = mod - 3, f[0] = 1;
	ODE A = ode_power(1).composite(pfrac(topoly(f)));
	A = A * B;
	A.getPRec();
	auto l = A.trans_poly();
	int d = 0;
	for (int i = 0; i < l.size(); ++ i) d = max(d, l[i].degree());
	cout << 1ll * gfac(n) * p_recur(n, l.size() - 1, d, { 1, n, 
		int(1ll * "x^2 + x - 4"_p.eval(n) * ginv(2) % mod), 
		int(1ll * "x^3 + 3x^2 - 9x - 19"_p.eval(n) * ginv(6) % mod), 
		int(1ll * "x^4 + 6x^3 - 9x^2 - 88x - 87"_p.eval(n) * ginv(24) % mod), 
		int(1ll * "x^5 + 10x^4 + 5x^3 - 220x^2 - 685x - 452"_p.eval(n) * ginv(120) % mod), 
		int(1ll * "x^6 + 15x^5 + 45x^4 - 365x^3 - 2715x^2 - 5517x - 2600"_p.eval(n) * ginv(720) % mod)
		}, l.data()) % mod << endl;
}