对边不交树结构的解析组合分析
《本文介绍一种小常数 \(O(n \log n)\) 计算 A091320(n, k) 一行的方法。》
定义 \(n\) 个节点的平面有根树是满足如下条件的树:
- \(n\) 个节点按标号在平面上顺序排列为一个正 \(n\) 边形,根为 \(1\) 号节点。
- 边是连接两个节点的直线段,任意两条边不在除端点外的位置交叉。
- 任意两个节点联通,且不存在环。
定义一棵平面有根树的叶子为只与一条边连接且不为根的节点。
给定 \(n\),请对每个 \(1\le k < n\) 计数 \(n\) 个节点、\(k\) 个叶子的平面有根树,对 \(998244353\) 取模。
这是去年(2023)八月的点子,但是想今年(2024)六月后再出。
然后是大家都熟知的事情了:任意二次分式的复合方法被公布;xiaoziyao 老师出了一道集训队互测。于是这道题的 gf 处理部分就完全 well-known 了。
而推导 gf 的部分(或得到 gf 的方法)可以在神 flajolet 的《对不交结构的解析组合分析》中找到。数列就是 A091320。那这题的所有部分就都 well-known 了。/dk
公开了。/fn
首先考虑不使用生成函数的做法。
定义 \(f(a, b, c)\) 为顺次排列 \(a\) 个点,其中第 \(b\) 个为根,有 \(c\) 个叶子的方案数。转移考虑枚举最右子树的范围,枚举右边的根,枚举叶子个数,有 \(O(n^6)\) dp。前缀和优化得到 \(O(n^5)\),插值维护得到 \(O(n^4)\)。
选用 \(O(n^5)\) 的做法即可过前三个 subtask。
code
感谢 kaguya 的代码。
#include <bits/stdc++.h>
using namespace std;
#define sz 305
#define mod 998244353
int n;
int dp[sz][sz][sz], sum[sz][sz];
int main()
{
scanf ("%d", &n);
dp[1][1][1] = 1;
for (int i = 1; i <= n; ++i)
for (int k = 1; k <= i; ++k)
{
dp[i][1][k] = (dp[i][1][k] + sum[i - 1][k]) % mod;
for (int j = 1; j <= i; ++j)
{
if (j > (i + 1 >> 1))
{
dp[i][j][k] = dp[i][i - j + 1][k];
continue;
}
for (int h = i; h > j && h > 2; --h)
for (int s = 1; s < k; ++s)
dp[i][j][k] = (dp[i][j][k] + 1ll * dp[h - 1][j][k - s] * sum[i - h + 1][s]) % mod;
}
for (int j = 1; j <= i; ++j)
sum[i][k] = (sum[i][k] + dp[i][j][k]) % mod;
}
for (int i = 1; i < n; ++i)
printf ("%d ", dp[n][1][i]);
return 0;
}
随后是正解。
不妨先不考虑 \(k\) 的限制,我们计数 \(n\) 个点的平面树。
考虑根节点的其中一棵子树。一个直接的想法就是将这棵子树拆分成左右两个部分,两个部分在平面图上分布在这颗子树的根两侧,这样和保证边不交等价。假设一棵子树的形态是 \(R-\bullet-L\),其中 \(L, R\) 为两个由子树组成的序列。可以发现,\(R-\bullet\) 和 \(\bullet - L\) 的结构也是平面树。这自然引出了下述的一种递推法。
定义拆分对是由两棵平面树组成的有序二元组,其大小为两棵平面树的大小之和减一,组合意义就是将这两棵平面树的根节点合并为一个。令 \(B(z)\) 为拆分对的 OGF,\(T(z)\) 为平面树的 OGF,模拟有序选择两棵平面树并合并根节点可以得到 \(B = {T^2}/z\)。由前面的讨论知道,一个拆分对就对应着一棵子树,因此模拟选取一个拆分对序列并连在根上可以得到 \(T = z\times \mathrm{SEQ}(B)\)。经过一些化简可以知道 \(T^3 - zT + z^2 = 0\)。
一个直接的想法就是找到复合逆的形式,而这就需要我们把 \(z\) 的次数化成相同的。令 \(T = z + zy\),我们就有 \(z^3(1 + y)^3 - z^2(1 + y) + z^2 = 0\),即 \(z(1 + y)^3 = y\)。这样我们就知道了 \(y^{\langle -1 \rangle}(z) = \dfrac {z}{(1 + z)^3}\),由拉格朗日反演知
随后考虑 \(k\) 的限制,这时加一维 \(w\) 表示叶子个数。
根据以往经验可以知道,我们总可以将不考虑 \(k\) 时的叶子(\(z\))换为带权叶子(\(zw\))得到正确的 BGF。分析可知 \(T = z\times \mathrm{SEQ}(B)\) 是不变的,因为这个过程中不会产生新叶子。接着考察 \(B(z)\)。分析组合意义可以知道 \(B[0] = 0, B[1] = 1\),这时应该将 \(B[1]\) 项标记,由于这是唯一的叶子。因此对 BGF \(B\) 有 \(B = \dfrac{T^2}z - z + zw\)。类似地化简可以得到 \(T^3 + (z^2 w - z^2 - z)T + z^2 = 0\)。注意到这个式子已经出现在了 oeis 上。
考虑 \(T = z + zy\),化简得 \(y = z(y + 1)(y^2 + 2y + w)\)。接着类似地处理可以得到
提取系数后自然能导出一种 \(O(n)\) 对单个 \(k\) 计算的方法,答案即为
code
int _2[N];
signed main () {
cin >> n;
_2[0] = 1;
rep(i,1,n) {
_2[i] = (_2[i - 1] << 1);
if (_2[i] >= mod) _2[i] -= mod;
}
rep(k,1,n - 1) {
int ans = 0;
rep(i,0,k-1) ans = (ans + 1ll * gC(n - 1, i) * gC(n - 1 - k, k - 1 - i) % mod * _2[n - 2 * k + i]) % mod;
cout << 1ll * C(n - 1, k) * ans % mod * ginv(n - 1) % mod << ' ';
}
cout << '\n';
}
下面考虑如何迅速对一行 \(k\) 提取
一个直接的结论是,如果维护两个多项式 \(F, G\),则提取一项系数可以模拟卷积,复杂度是 \(O\left(\min\{\mathrm{deg} F, \mathrm{deg} G\}\right)\) 的。可以设置阈值 \(B\),并维护两个多项式 \(F, G\)。初始令 \(F = (z + 1)^n, G = 1\),顺次计算 \(k\),每当 \(\mathrm{deg} G > B\) 就将 \(G\) 乘入 \(F\)。这样有 \(n-1\) 次 \(O(B)\) 的求值,以及 \(O(n/B)\) 次 \(O(B)\) 次多项式乘 \(O(n)\) 次多项式,复杂度是 \(O\left(n\times B + n\log n\times\dfrac nB\right)\),取 \(B = \sqrt{n \log n}\) 得到复杂度 \(O(n\sqrt{n \log n})\)。
code
signed main () {
cin >> n;
int B = sqrt(n * log(n));
poly f = "x + 1"_p.redegree(n - 1).pow(n - 1), g = "1"_p;
rep(k,0,n - 2) {
int ans = 0;
rep(i,0,min(g.degree(), n - 2 - k)) ans = (ans + 1ll * g[i] * f[n - 2 - k - i]) % mod;
a[n - 1 - k] = ans;
if (k < n - 2) {
g *= "x + 2"_p;
if (g.degree() > B) f *= g, g = "1"_p;
}
}
rep(i,1,n - 1) a[i] = 1ll * a[i] * gC(n - 1, i) % mod * ginv(n - 1) % mod;
rep(i,1,n - 1) cout << a[i] << ' ';
cout << '\n';
}
回到计算 \(c_k = [z^{n - 2}] (z + 1)^{n - 1}(z^2 + 2z)^{n - 1 - k}\) 上来。
为方便,将 \(n\) 减一,并使 \(a_k = c_{n - k} = [z^{n - 1}] (z + 1)^n(z^2 + 2z)^{k}\)。
令 \(b_k = \dbinom{n}{k + 1}\),矩阵 \(A\) 满足 \(A[k, i] = [z^i](z^2 + 2z)^k\),则我们知道
考虑转置原理。设向量 \(\bm a',\bm b'\) 使得 \(\bm a' = A^{\mathsf T} \bm b'\),即
令 \(A'(z)\) 满足 \(A'[i] = a'_i\),\(B'(z)\) 满足 \(B'[i] = b'_i\),则可以知道
由于 \(B'(x^2 + 2x) = B' \circ (x - 1) \circ x^2 \circ (x + 1)\),可以知道这是 \(\bm a' =M_{+1}M_{^\wedge 2} M_{-1} \bm b'\),其中 \(M_{+1}\) 为对左乘的向量对应多项式复合 \(x + 1\),另两个同理。因此由转置原理有 \(\bm a = M_{-1}^{\mathsf T}M_{^\wedge 2}^{\mathsf T}M_{+1}^{\mathsf T} \bm b\)。下面需要分开讨论。
1. \(M_{+ 1}^{\mathsf T}\)
考虑
则 \(b_k = \sum_{i\ge k} \binom{i}{k} a_i\),对应 \(M_{-1}[k, i] = \binom{i}{k}\),转置得 \(M_{-1}^{\mathsf T} [k, i] = \binom{k}{i}\),于是转置变换即为
egf 卷 \(e^x\) 即可。
2. \(M_{^\wedge2}^{\mathsf T}\)
可以知道 \(M_{^\wedge2}[k, i] = [2i = k]\),转置得 \(M_{^\wedge 2}^{\mathsf T}[k, i] = [2k = i]\),即
3. \(M_{-1}^\mathsf T\)
与 \(1.\) 类似。
转置有
egf 卷 \(e^{-x}\) 即可。
总时间复杂度 \(O(n\log n)\)。
code
signed main () {
cin >> n; n --;
fac[0] = inv[0] = inv[1] = ifc[0] = 1;
rep(i,1,2 * n) fac[i] = 1ll * i * fac[i - 1] % mod;
rep(i,2,2 * n) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
rep(i,1,2 * n) ifc[i] = 1ll * inv[i] * ifc[i - 1] % mod;
poly f(2 * n + 1), ex(2 * n + 1);
rep(i,0,n - 1) f[i] = 1ll * C(n, i + 1) * ifc[i] % mod;
// stage 1: [x + 1]^T
ex.redegree(2 * n);
rep(i,0,2 * n) ex[i] = ifc[i];
f = f * ex;
f.redegree(2 * n);
// stage 2: [x^2]^T
rep(i,0,n) f[i] = 1ll * f[2 * i] * fac[2 * i] % mod * ifc[i] % mod;
f.redegree(n);
// stage 3: [x - 1]^T
rep(i,0,n) ex[i] = ((i & 1) ? mod - ifc[i] : ifc[i]);
f = f * ex;
rep(k,1,n) c[k] = 1ll * f[n - k] * fac[n - k] % mod;
n ++;
rep(k,1,n-1) cout << 1ll * inv[n - 1] * C(n - 1, k) % mod * c[k] % mod << ' ';
}
lyin 给了一种思路:
仍然考察一个非根结点能贡献 \((子树 +1)\) 的权,并且知道 \(\dfrac{\mathrm d G(1-G)^{-1}}{\mathrm d G} = \dfrac{1}{(1 - G)^2}\),可以得到 \(G\) 的式子。随后对根不带权地选择子树能得到上面的两个式子。
再往后推的话,可以知道
化简知道
于是有
玄妙。后面不想推了。
以下是博客签名,与正文无关。
请按如下方式引用此页:
本文作者 joke3579,原文链接:https://www.cnblogs.com/joke3579/p/non-crossing-tree.html。
遵循 CC BY-NC-SA 4.0 协议。
请读者尽量不要在评论区发布与博客内文完全无关的评论,视情况可能删除。