题解 P5900 【无标号无根树计数】
Weng_Weijie
2020-01-05 20:29:13
## 前言
虽然 $\Theta(n\log^2n)$ 也能过，但是没什么意义，只是**跑的就是比 $\Theta(n\log n)$ 快**。
在我的博客里看可能公式会更好看。
## Euler 变换简介
为了方便，这里介绍一下幂级数的 $\text{Euler}$ 变换。
$\text{Euler}$ 变换的其中一个定义式如下：
$$\mathcal E(F(x))=\prod_i{\dfrac{1}{(1-x^i)^{f_i}}}$$
这个变换类似于 $\text{EGF}$ 中的 $\exp$。
$$\exp F(x)=\sum_i\dfrac{F^i(x)}{i!}$$
$\text{EGF}$ 的 $\exp$ 具有的组合意义是：将 $1\sim n$ 分成若干**非空**组，大小为 $i$ 的组内部具有 $f_i$ 中方案，问最后的总方案数。
而 $\text{Euler}$ 变换就是**去除了这个标号**的方案数，去掉标号会导致「如果两个组**大小和内部方案均相同**，则它们不可区分」。
因此 $\text{Euler}$ 的第一个定义式就易懂了。（即大小为 $i$ 的每种内部方案都可以选若干个，每种内部方案的生成函数都是 $(1-x^i)^{-1}$。）
现在有两种方法可以得到 $\text{Euler}$ 变换的第二个定义式。
### 方法一：使用指数、对数方法
$$\mathcal E(F(x))=\exp\left(\sum_i\ln\dfrac{1}{(1-x^i)^{f_i}}\right)
$$
$$=\exp\left(\sum_i-f_i\ln(1-x^i)\right)$$
$$=\exp\left(\sum_if_i\sum_{j}\dfrac{x^{ij}}j\right)$$
$$=\exp\left(\sum_j\dfrac 1j\sum_if_ix^{ij}\right)$$
$$=\exp\left(\sum_j\dfrac{F(x^j)}j\right)$$
因此我们得到了 $\text{Euler}$ 变换的第二个定义式。
### 方法二：使用 $\text{Burnside}$ 引理 / $\text{Pólya}$ 定理
首先枚举分成几个组，设有 $n$ 个。
这样的话，容易发现定理中的置换群即为 **$n$ 元对称群 $S_n$**。
**定理内容即：对所有的 $n$ 元置换，求在该置换作用下不动点个数的平均值。**
这时候考虑 $k$ 个位置构成 $k$ 阶循环的**指数生成函数**：
此时这 $k$ 个位置必须选择同样的方案，因此可以看做一个某个方案大小扩大了 $k$ 倍。同时 $k$ 个位置构成循环方案数即圆排列为 $(k-1)!$。
因此这个 $\text{EGF}$ 为
$$C_k(x)=\dfrac{1}{k!}(k-1)!\cdot F(x^k)=\dfrac{F(x^k)}{k}$$
根据 $\exp$ 的组合意义，可以想到将所有 $C_k(x)$ 相加，再做 $\exp$ 得到所有置换的贡献。
事实上这是正确的，原因是定理中的除以置换群大小（即 $n!$）与 $\text{EGF}$ 得到方案数时乘的 $n!$ 抵消了。（难理解的话可以**多设一个变元表示置换的大小**以更好地理解）
同样得到了那个定义式：
$$\mathcal E(F(x))=\exp\left(\sum_k\dfrac{F(x^k)}k\right)$$
## 回到原问题
### 解决无标号有根树计数问题
首先解决**无标号有根树**的问题。
容易发现，给定根以后，所有的子树可以理解成构成大小为 $n-1$ 的若干个组。
因此我们需要解的方程即：
$$F(x)=x\cdot\mathcal E(F(x))$$
#### 方法一：先求导再使用分治 $\text{FFT}$ 解决
将上式两边求导再同时乘上 $x$（即使用 $\vartheta$ 算子，$x^n$ 项的系数乘以 $n$）。
这里略去了一些中间步骤。
$$\vartheta F(x)=F(x)+F(x)\left(\sum_{k\geq 1}F'(x^k)x^k\right)$$
设 $G(x)=\sum\limits_kF'(x^k)x^k$。
可以观察到 $g_n=\sum\limits_{d|n}f_d\cdot d$。
因此 $f_n=\dfrac{1}{n-1}\left(\displaystyle\sum_{k}f_kg_{n-k}\right)$，$f_1=1$。
此时便可以使用分治 $\text{FFT}$ 解决这个问题。
#### 方法二：$\text{Newton}$ 迭代法
我们需要解的方程是：
$$G(F(x))=F(x)-x\cdot\mathcal E(F(x))$$
假设当前已经求出了 $F(x)\bmod{x^n}$，需要求 $F(x)\bmod{x^{2n}}$。
容易发现的是对 $k\geq 2$，我们已经知道了 $F(x^k)\mod {x^{2n}}$ 的值。
这时我们可以用一个常数来代替它。（强调一下**这不代表它原来是一个常数**，也就是说 $F(x^2)$ 对 $F(x)$ 求导并不是 $0$）。
设
$$P(x)=\left(x\exp\left(\sum_{k\geq 2}\dfrac{F(x^k)}k\right)\right)\bmod{x^{2n}}$$
由于 $P(x)$ 是我们规定的一个常数，则
$$F(x)\equiv P(x)\exp F(x)\pmod{x^{2n}}$$
这个形式就易于 $\text{Newton}$ 迭代了，具体过程不详细介绍了。
事实上这个方法需要做 $\exp$，导致常数比较大，实现不优秀很有可能比前一个方法慢。
### 解决无标号无根树计数问题
这个问题思路大概是：用有根树的方案减去根不是重心的方案。
**当 $n$ 是奇数时**：
如果根不是重心，必然存在恰好一个子树，它的大小超过 $\left\lfloor\dfrac n2\right\rfloor$（设它的大小为 $k$），那么这个子树和这棵子树以外的其他点的方案数恰好为 $f_k\cdot f_{n-k}$（可以看成两棵有根树）。
因此答案为
$$f_n-\sum_{k=\lfloor\frac n2\rfloor+1}^{n-1}f_k\cdot f_{n-k}$$
**当 $n$ 是偶数时**：
出现的问题是，有可能存在两个重心，且其中一个是根（即存在一棵子树大小恰为 $\dfrac n2$）。
那么如果这个子树和剥离该子树后的树完全相同，这样的方案在 $f_n$ 只会计算一次，即不需要减去。
如果这个子树和另一棵树不相同，会被统计两次，需要减掉一次。
因此偶数时额外减掉的方案数为
$$\binom{f_{\frac n2}}2$$
**至此，这个问题在 $O(n\log n)$ 或 $O(n\log^2n)$ 的时间复杂度内解决**。
### 参考代码：
```cpp
#include <bits/stdc++.h>
typedef long long LL;
typedef unsigned long long ULL;
const int mod = 998244353, N = 262144;
int wn[N], lim, s, rev[N], w[N];
void reduce(int &x) { x += x >> 31 & mod; }
int pow(int x, int y, int ans = 1) {
for (; y; y >>= 1, x = (LL) x * x % mod)
if (y & 1) ans = (LL) ans * x % mod;
return ans;
}
void fftinit(int len) {
wn[0] = lim = 1, s = -1; while (lim < len) lim <<= 1, ++s;
for (int i = 0; i < lim; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << s;
const int g = pow(3, (mod - 1) / lim);
for (int i = 1; i < lim; ++i) wn[i] = (LL) wn[i - 1] * g % mod;
}
void fft(int *A, int typ) {
static ULL tmp[N];
for (int i = 0; i < lim; ++i) tmp[rev[i]] = A[i];
for (int i = 1; i < lim; i <<= 1) {
for (int j = 0, t = lim / i / 2; j < i; ++j) w[j] = wn[j * t];
for (int j = 0; j < lim; j += i << 1)
for (int k = 0; k < i; ++k) {
const ULL x = tmp[k + j + i] * w[k] % mod;
tmp[k + j + i] = tmp[k + j] + mod - x, tmp[k + j] += x;;
}
}
for (int i = 0; i < lim; ++i) A[i] = tmp[i] % mod;
if (!typ) {
const int il = pow(lim, mod - 2); std::reverse(A + 1, A + lim);
for (int i = 0; i < lim; ++i) A[i] = (LL) A[i] * il % mod;
}
}
int inv[N];
int f[N], g[N], n;
int a[N], b[N];
void init(int n) {
inv[0] = 1, inv[1] = 1;
for (int i = 2; i < n; ++i)
inv[i] = (LL) (mod - mod / i) * inv[mod % i] % mod;
}
void solve(int l, int r) {
if (r - l <= 32) {
for (int i = l; i < r; ++i) {
for (int j = l; j < i; ++j) {
f[i] = (f[i] + (LL) f[j] * g[i - j]) % mod;
if (l > 1) f[i] = (f[i] + (LL) g[j] * f[i - j]) % mod;
}
f[i] = (LL) f[i] * inv[i - 1] % mod;
int v = (LL) f[i] * i % mod;
for (int p = i; p <= n; p += i) reduce(g[p] += v - mod);
}
return;
}
int mid = l + r + 1 >> 1;
solve(l, mid), fftinit(r - l);
auto update = [&] (int *f, int la, int *g, int lb) {
std::memcpy(a, f, la << 2), std::memset(a + la, 0, lim - la << 2);
std::memcpy(b, g, lb << 2), std::memset(b + lb, 0, lim - lb << 2);
fft(a, 1), fft(b, 1);
for (int i = 0; i < lim; ++i) a[i] = (LL) a[i] * b[i] % mod;
fft(a, 0);
for (int i = mid; i < r; ++i) reduce(::f[i] += a[i - l] - mod);
};
if (l > 1) update(f + l, mid - l, g, r - l);
update(g + l, mid - l, f, l == 1 ? mid : r - l);
solve(mid, r);
}
int main() {
std::ios::sync_with_stdio(0), std::cin.tie(0);
std::cin >> n, f[1] = 1, init(n), solve(1, n + 1);
int ans = f[n];
for (int i = n / 2 + 1; i < n; ++i)
ans = (ans + (LL) (mod - f[i]) * f[n - i]) % mod;
if (n % 2 == 0) ans = (ans + mod - (LL) f[n / 2] * (f[n / 2] - 1) / 2 % mod) % mod;
std::cout << ans << '\n';
return 0;
}
```