无标号有度数限制基环树森林计数

无标号有度数限制基环树森林计数

写这个的原因是前几天模拟赛考了这个，这个计数还是很有意思的，过程比较有启发性。

思路和代码均借鉴该博客，本文更多的是思路方向和细节上的一些补充。

我们要求的东西是给定 \(n,k\)，求 \(n\) 个点，每个点入度不超过 \(k\) 的无标号内向基环树森林个数。

大体的思路是先求出 \(dp_n\) 表示 \(n\) 个点，度数限制下无标号有根树的方案数，再挂到环上算出 \(n\) 个点的基环树的方案数，最后合并成森林的方案数。

求无标号有根树方案数

先不考虑度数限制，那问题的关键是去重。暴力容斥什么的还是太吃操作了，考虑直接上生成函数来刻画这个情况。直接枚举子树大小，对于每种子树大小为了防止重复，在一起统计方案数。那么对于大小为 \(i\) 的子树枚举选它的个数，实际上方案数就是 \(x^0+x^i+x^{2i}+\cdots=\frac1{1-x^i}\)，因此我们有：

\[dp_n=[x^{n-1}]\prod_{i=1}^\infty(\frac1{1-x^i})^{dp_i} \]

考虑加上度数限制。这个东西直接加一维硬算。那么有：

\[dp_n=\sum_{k\le K}[x^{n-1}y^k]\prod_{i=1}^\infty(\frac1{1-x^iy})^{dp_i} \]

加和的形式还是太吃操作了，考虑给这个东西乘上 \(\frac1{1-y}\) 实际上就能把所有 \(\le K\) 的 \(k\) 系数为 \(1\) 地贡献到 \(K\)，因此有：

\[dp_n=[x^{n-1}y^k]\prod_{i=1}^\infty(\frac1{1-x^iy})^{dp_i}\frac1{1-y} \]

现在的问题是转移的系数不好算。考虑如何算 \([x^t](\frac1{1-x^iy})^{dp_i}\) 。记 \(w=x^iy\)，展开之后就是 \((w^0+w^1+\cdots)^{dp_i}\) 状物，那么求第 \(t\) 项系数就是等价于求 \(dp_i\) 个非负整数的和等于 \(t\)，这个东西是 \(\binom{t+dp_i-1}{t}\)。

那 \(dp_i\) 这个东西很大，取模带进去算对吗？啊是对的，考虑卢卡斯定理有 \(\binom{n}{m}=\binom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\binom{n\bmod p}{m\bmod p}\)，而下指标 \(t\) 很小，有 \(\lfloor\frac tp\rfloor=0\)，因此实际上直接带进去取模就是对的。

然后背包维护就是容易的了，对于复杂度，实际上单次转移 \(x\) 的信息量只有 \(n\log n\)，这个东西是 \(\sum_{i=1}^n\frac ni\) 的经典复杂度，因此总的复杂度是 \(O(n^3\log n)\)。

需要注意的是最后的根要多留一条边给基环树。

求无标号内向基环树方案数

这个东西要在环上做，那直接上 Burnside 引理。

\[|X/G|=\frac1{|G|}\sum_{g\in G}|X^g| \]

考虑记一个 \(F(n,m)\) 表示环上 \(n\) 个位置，每个位置的大小是 \(a_i\)，且 \(\sum a_i=m\) 的方案数。那么枚举位置数量 \(h\) 和旋转置换环旋转的个数 \(t\)，那么实际要统计的位置数是 \(\gcd(h,t)\)，放置的点数是 \(\frac{n}{\frac{h}{\gcd(h,t)}}=\frac{n\gcd(h,t)}{h}\)，于是上引理有：

\[tr_n=\sum_{h=1}^n\frac1h\sum_tF(\gcd(h,t),\frac{n\gcd(h,t)}{h}) \]

那 \(F\) 的方案数也好求了：

\[F(n,m)=\sum_{i=1}^mF(n-1,m-i) \]

Final Part

然后就是把若干基环树合成森林。这个合成的过程和一开始求树的方案数的方法是一致的：

\[ans_n=[x^n]\prod_{i=1}^n(\frac1{1-x^i})^{tr_i} \]

直接暴力做一下就好了。

下面给出代码：

#include <bits/stdc++.h>
using namespace std;
const int N = 255, mod = 1e9 + 7;
inline void add(int &x, int y) {x = x + y >= mod ? x + y - mod : x + y;
}
int qpow(int x, int y) {int ans = 1;while (y) {if (y & 1) ans = 1ll * ans * x % mod;x = 1ll * x * x % mod;y >>= 1;}return ans;
}
int fac[N], inv[N];
int C(int n, int m) {if (n < m || n < 0 || m < 0) return 0;int ans = 1;for (int i = m + 1; i <= n; i++) ans = 1ll * ans * i % mod;return 1ll * ans * inv[n - m] % mod;
}int wn, wk;
int dp[N], f[N][N], num[N];int main() {fac[0] = 1;for (int i = 1; i < N; i++) fac[i] = 1ll * fac[i - 1] * i % mod;inv[N - 1] = qpow(fac[N - 1], mod - 2);for (int i = N - 2; ~i; --i) inv[i] = 1ll * inv[i + 1] * (i + 1) % mod;cin >> wn >> wk;for (int i = 0; i <= wk; i++) f[0][i] = 1;dp[1] = 1;for (int n = 1; n <= wn; n++) {vector<pair<int, int>>v;int r = n, t = 0;while (r <= wn) ++t, v.emplace_back(t, C(t + dp[n] - 1, dp[n] - 1)), r += n;for (int x = wn; ~x; --x)for (int y = wk; ~y; --y)for (auto p : v) {int i = p.first, w = p.second;if (i * n <= x && i <= y) add(f[x][y], 1ll * f[x - i * n][y - i] * w % mod);}dp[n + 1] = f[n][wk];}for (int i = 1; i <= wn; i++) dp[i] = f[i - 1][wk - 1];memset(f, 0, sizeof f);f[0][0] = 1;for (int i = 1; i <= wn; i++)for (int j = 1; j <= wn; j++)for (int k = 1; k <= j; k++)add(f[i][j], 1ll * f[i - 1][j - k] * dp[k] % mod);for (int n = 1; n <= wn; n++)for (int h = 1; h <= n; h++) {int iv = qpow(h, mod - 2);for (int t = 1; t <= h; t++) {int g = __gcd(h, t);if (n * g % h == 0) add(num[n], 1ll * f[g][n * g / h] * iv % mod);}}memset(dp, 0, sizeof dp);dp[0] = 1;for (int n = 1; n <= wn; n++) {vector<pair<int, int>>v;int r = n, t = 0;while (r <= wn) ++t, v.emplace_back(t, C(t + num[n] - 1, num[n] - 1)), r += n;for (int x = wn; ~x; --x)for (auto p : v) {int i = p.first, w = p.second;if (i * n <= x) add(dp[x], 1ll * dp[x - i * n] * w % mod);}}cout << dp[wn] << '\n';return 0;
}