P5900 无标号无根树计数 题解

· · 题解

无标号无根树计数

这道题题解区有很多神仙题解，讲解得很全面，本篇题解仅做一些细节补充和代码实现方面的讲解。

题意

求 n 个点的无标号无根树数量，答案对 998,244,353 取模。(1\le n\le 2\times10^5)

题解

注意到无根树这个限制比较恶心，因为对于一棵树我们总要找个根才好处理各种东西。所以考虑先计算无标号有根树的数量，然后减去根不是重心的方案。

而对于有根树，钦定完一个根之后剩下 n-1 个点会分成若干个非空集合（且非空集合之间顺序无所谓，这也是为什么后文用的 OGF 而不是 EGF），每个集合对应一棵子树，而每棵子树又是类似的结构，即具有递归的形式。这样如果我们设答案的 OGF 为：

F(x)=\sum_{n\ge 0}f_nx^n

则它应该满足等式：

[x^n]F(x)=[x^{n-1}]\prod_{k\ge 1}\left(\sum_{i\ge 0}x^{ik}\right)^{f_k}

其实就是把选择集合的过程当做背包，每种集合的大小当成一种物品，会贡献 f_i 的方案数，结合上面说的无标号有根树的性质很容易得出这个等式。根据：

\left(\dfrac{1}{1-x^k}\right)^{f_k}=\left(\sum_{i\ge 0}x^{ik}\right)^{f_k}

能把原式的 \sum 干掉：

[x^n]F(x)=[x^{n-1}]\prod_{k\ge 1}\left(\dfrac{1}{1-x^k}\right)^{f_k}

这里再简单乘上一个 x 就能把 [x^n],[x^{n-1}] 干掉了：

F(x)=x\prod_{k\ge 1}\left(\dfrac{1}{1-x^k}\right)^{f_k}

看来我们问题的关键就落在了 \prod_{k\ge 1}\left(\frac{1}{1-x^k}\right)^{f_k} 这个东西上面。注意 \prod 是很难处理的，所以考虑取对数转化为 \sum：

\begin{aligned}\prod_{k\ge 1}\left(\dfrac{1}{1-x^k}\right)^{f_k}&=\exp\left(\sum_{k\ge 1}\ln\left(\dfrac{1}{1-x^k}\right)^{f_k}\right)\\&=\exp\left(\sum_{k\ge1}-f_k\ln(1-x^k)\right)\end{aligned}

现在压力来到了 \ln(1-x^k) 这边。因为 \ln 很难处理，但它的导数很好处理，所以考虑求导：

\ln(1-x^k)=\int\dfrac{\mathrm{d}}{\mathrm{d}x}\ln(1-x^k)=\int\dfrac{-kx^{k-1}}{1-x^k}

我们能看到现在积分号里面的就是一个我们知道对应幂级数的封闭形式了：

\int\dfrac{-kx^{k-1}}{1-x^k}=-\int\sum_{n\ge 0}kx^{(n+1)k-1}=-\int\sum_{n\ge 1}kx^{nk-1}

现在这个式子就很好看了，只需要把积分号算进去：

-\int\sum_{n\ge 1}kx^{nk-1}=-\sum_{n\ge 1}\dfrac{kx^{nk-1+1}}{nk-1+1}=-\sum_{n\ge 1}\dfrac{x^{nk}}{n}

好，这东西的另一种形式显然更便于式子的化简，我们代回去：

\begin{aligned}\exp\left(\sum_{k\ge1}-f_k\ln(1-x^k)\right)&=\left(\exp\sum_{k\ge1}f_k\sum_{n\ge 1}\dfrac{x^{nk}}{n}\right)\\&=\exp\left(\sum_{n\ge 1}\dfrac{1}{n}\sum_{k\ge 1}f_kx^{nk}\right)\\&=\exp\left(\sum_{n\ge 1}\dfrac{F(x^n)}{n}\right)\end{aligned}

啊这就好多了，注意到我们有等式：

F(x)=x\exp\left(\sum_{n\ge 1}\dfrac{F(x^n)}{n}\right)

还是类似 \ln，遇到 \exp 考虑求导：

\begin{aligned}F'(x)&=\left(x\exp\left(\sum_{n\ge 1}\dfrac{F(x^n)}{n}\right)\right)'\\&=x'\exp\left(\sum_{n\ge 1}\dfrac{F(x^n)}{n}\right)+x\left(\exp\left(\sum_{n\ge 1}\dfrac{F(x^n)}{n}\right)\right)'\\&=\dfrac{F(x)}{x}+x\exp\left(\sum_{n\ge 1}\dfrac{F(x^n)}{n}\right)\left(\sum_{n\ge 1}\dfrac{F(x^n)}{n}\right)'\\&=\dfrac{F(x)}{x}+F(x)\sum_{n\ge 1}\dfrac{F'(x^n)nx^{n-1}}{n}\\&=\dfrac{F(x)}{x}+F(x)\sum_{n\ge 1}F'(x^n)x^{n-1}\end{aligned}

两边同时乘上 x，会让式子更好看：

xF'(x)=F(x)+F(x)\sum_{n\ge 1}F'(x^n)x^n

注意到这个等式哪哪都好，就是最后一个求和式很恶心。考虑定义 G(x)=\sum_{n\ge 1}F'(x^n)x^n，并研究一下 [x^i]G(x) 的取值。让我们把 G(x) 展开：

\begin{aligned}G(x)&=\sum_{n\ge 1}x^n\sum_{j\ge 1}jf_jx^{n(j-1)}\\&=\sum_{n\ge 1}\sum_{j\ge 1}jf_jx^{nj}\end{aligned}

这样的话，对于 x^i 做贡献的所有 jf_j 应该满足 j|i，即：

[x^i]G(x)=\sum_{j|i}f_jj

把 G(x) 带回原式：

xF'(x)-F(x)=F(x)G(x)

而我们发现 xF'(x) 其实就是：

\sum_{n\ge 1}nf_nx^n

所以：

[x^n]F(x)=\dfrac{[x^n]F(x)G(x)}{n-1}

即：

f_n=\dfrac{\sum_{k=1}^{n-1}(f_kg_{n-k})}{n-1},f_1=1

半在线卷积求解即可。不过这里的半在线卷积与谷谷 分治 FFT 板子 的半在线卷积并不是一种。注意那道题 g 是给出的，而这里 g 是依赖于前面求出的 f 的。在那道题中，我们用到的半在线卷积依赖于这个式子：

f_{l\sim mid}\times g_{0\sim r-l}\rightarrow f_{mid\sim r}

来算左边对右边的贡献。但这里因为只求出了 l\sim mid 的 f，g_{mid,\sim r-l} 的值还没有计算，这题这样做显然是错误的。而正确的策略是，当 l=0 时：

f_{l\sim mid}\times g_{0\sim mid}\rightarrow f_{mid\sim r}

这时候显然有贡献漏掉了。考虑当 l>0 时补回来：

\begin{cases}g_{l\sim mid}\times f_{0,r-l}\\f_{l\sim mid}\times g_{0,r-l}\end{cases}\rightarrow f_{mid\sim r}

这样就能补回来漏掉的贡献了。

现在求出来 F(x) 了，不要忘了我们要计数的是无根树不是有根树。所以要减掉根不是重心的方案。当 n 是奇数时，显然原树只能有一个重心。而重心的定义是所有子树大小均小于等于 \lfloor\frac{n}{2}\rfloor，所以我们只需要钦定一棵子树大小大于 \lfloor\frac{n}{2}\rfloor 即可，要减去的方案数为：

\sum_{k=\lfloor\frac{n}{2}\rfloor+1}^{n-1}f_kf_{n-k}

即看成两棵有根树分别处理，然后把大小为 k 的那个接到另一个的根上。而当 n 为偶数时，麻烦了，因为可能出现两个重心的情况，且其中一个还是根，即存在一棵子树大小恰等于 \frac{n}{2}。而现在还是有两种情况，第一种是这棵子树和去掉这棵子树后的原树（类比奇数情况下 f_k 和 f_{n-k}）恰好同构，那这个就不会重复计数了，因为这俩根本质上是一样的。而其他情况就会重复计数一次，减掉的方案数即在所有可能的形态中选俩接起来：

\dbinom{f_{\frac{n}{2}}}{2}

这样这道题就解决了，时间复杂度 \mathcal{O}(n\log^2n)。在刚刚的推导中我把 f_0 这一项都忽略了，但在实现里面我是算上了的，但显然影响不大。

#include <cstdio>
#include <algorithm>
const int N = 1e6 + 10, mod = 998244353; typedef long long ll;
int f[N], g[N], A[N], B[N], rev[N], lim, n, m;
inline int ksm(int a, int b)
{
    int ret = 1;
    while (b)
    {
        if (b & 1) ret = (ll)ret * a % mod;
        a = (ll)a * a % mod; b >>= 1;
    }
    return ret;
}
inline void init(int n)
{
    lim = 1; m = 0; while (lim <= n) lim <<= 1, ++m;
    for (int i = 0; i < lim; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (m - 1));
}
inline void NTT(int* f, int len, int on)
{
    for (int i = 0; i < len; ++i) if (i < rev[i]) std::swap(f[i], f[rev[i]]);
    for (int h = 2; h <= len; h <<= 1)
    {
        int gn = ksm(3, (ll)on * ((mod - 1) / h) % (mod - 1));
        for (int j = 0; j < len; j += h)
            for (int k = j, g = 1; k < j + h / 2; ++k, g = (ll)g * gn % mod)
            {
                int u = f[k], t = (ll)g * f[k + h / 2] % mod;
                f[k] = (u + t) % mod; f[k + h / 2] = ((u - t) % mod + mod) % mod;
            }
    }
    if (on == mod - 2) for (int i = 0, inv = ksm(len, on); i < len; ++i) f[i] = (ll)f[i] * inv % mod;
}
inline void mul(int* f, int lF, int* g, int lG)
{   
    for (int i = 0; i < lF; ++i) A[i] = f[i];
    for (int i = lF; i < lim; ++i) A[i] = 0;
    for (int i = 0; i < lG; ++i) B[i] = g[i];
    for (int i = lG; i < lim; ++i) B[i] = 0;
    NTT(A, lim, 1); NTT(B, lim, 1); 
    for (int i = 0; i < lim; ++i) A[i] = (ll)A[i] * B[i] % mod;
    NTT(A, lim, mod - 2);
}
inline void cdq(int l, int r)
{
    if (r == 1) return ;
    if (l + 1 == r) 
    {
        if (l != 1) f[l] = (ll)f[l] * ksm(l - 1, mod - 2) % mod;
        for (int i = l; i < n; i += l) (g[i] += (ll)f[l] * l % mod) %= mod;
        return ;
    }
    int mid = (l + r) >> 1; cdq(l, mid); init(r - l);
    if (l > 0) 
    {
        mul(f, r - l, g + l, mid - l);
        for (int i = mid; i < r; ++i) (f[i] += A[i - l]) %= mod;
    }
    mul(f + l, mid - l, g, !l ? mid : r - l);
    for (int i = mid; i < r; ++i) (f[i] += A[i - l]) %= mod;
    cdq(mid, r);
}
signed main()
{
    int m; scanf("%d", &m); n = 1; while (n <= m) n <<= 1;
    f[1] = 1; cdq(0, n); int ans = f[m];
    for (int i = m / 2 + 1; i < m; ++i) (ans += (ll)(mod - f[i]) * f[m - i] % mod) %= mod;
    if (!(m & 1)) (ans += (mod - (ll)f[m / 2] * (f[m / 2] - 1) / 2 % mod)) %= mod;
    printf("%d\n", ans); return 0;
}

题做完了，我们来看看从这道题我们能得到什么。注意到刚刚我们处理的式子：

\prod_{k\ge 1}\left(\frac{1}{1-x^k}\right)^{f_k}

的组合意义是把原集合分成若干个非空子集，且顺序无所谓，大小为 i 的组方案数为 f_i，对应的总方案数。这个东西跟 EGF 中的 \exp 似乎很相似，所以我们考虑给它单独起一个名字：Euler 变换，F(x) 的 Euler 变换通常记为 \mathcal{E}F(x)。刚刚已经得出了它的两种形式：

\mathcal{E}F(x)=\prod_{k}\left(\frac{1}{1-x^k}\right)^{f_k}=\exp\left(\sum_{n}\dfrac{F(x^n)}{n}\right)

个人认为算是 OGF 版的 “\exp”，不考虑分成子集的前后顺序。

ps. 半在线卷积的部分参考 \tt command\_ block 的 博客。

pps. 本题还可以用牛顿迭代在 \mathcal{O}(n\log n) 的时间内求解，但那个方法涉及到 \exp，所以常数会比较爆炸，具体可以看其他神仙的题解。

ppps. 关于最后扔出的 Euler 变换第二个形式的证明，题解区还给出了用群论方法证明的思路。

pppps. 我的代码常数大得要命，看看理解咋实现可以，借鉴具体实现方式就算了。