「KDOI-07」能量场 题解
Xseventh
·
·
题解
题目要求的是生成树计数,但是是基环树。
考虑使用矩阵树定理。
定义矩阵 D 和 A,其中 D_{i,i} = \sum_{j=1}^n a_i+a_j = na_i + \sum_{j=1}^na_j,A_{i,j}=a_i+a_j。显然如果我们求的是正常的生成树,则 det(D-A) 即为答案。
注意到基环树的环缩起来也是一棵树。即对于基环树上环中点集 S,在矩阵 D - A 上删除与 S 有关的行列并求出行列式,本质与将环缩为一个点,只考虑环外贡献的答案相同。
因此我们有了基环树计数的初步方法,即枚举环上集合 S,将环上贡献与环外贡献分开计算,环外贡献使用矩阵树定理,环内可以直接状压打牌解决,时间复杂度 O(2^nn^3)。
矩阵 D - A 本身的生成具有一定性质,考虑是否可以利用。对于矩阵 A,将其第 2 到 n 行都减去第 1 行,即其原本的值是这样:
\begin{matrix}
a_1+a_1 & a_1+a_2 & a_1+a_3 & \dots & a_1+a_n \\
a_2+a_1 & a_2+a_2 & a_2+a_3 & \dots & a_2+a_n \\
a_3+a_1 & a_3+a_2 & a_3+a_3 & \dots & a_3+a_n \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
a_n+a_1 & a_n+a_2 & a_n+a_3 & \dots & a_n+a_n \\
\end{matrix}
操作后变为:
\begin{matrix}
a_1+a_1 & a_1+a_2 & a_1+a_3 & \dots & a_1+a_n \\
a_2-a_1 & a_2-a_1 & a_2-a_1 & \dots & a_2-a_1 \\
a_3-a_1 & a_3-a_1 & a_3-a_1 & \dots & a_3-a_1 \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
a_n-a_1 & a_n-a_1 & a_n-a_1 & \dots & a_n-a_1 \\
\end{matrix}
发现这几行每行值都相等,即矩阵 A 的秩至多为 2。
把行列式按照定义展开,则 det(D-A) = \sum_p(-1)^{inv(p)} \prod_{i=1}^n(D-A)_{i,p_i} = \sum_p(-1)^{inv(p)} \prod_{i=1}^n(D_{i,p_i} - A_{i,p_i})。显然 D_{i,p_i} 和 A_{i,p_i} 中必然有一个为 0,又因 i=p_i 时 D_{i,p_i} 才有贡献,则矩阵 A 的贡献就类似于其子矩阵的行列式。之前已知 A 的秩至多为 2,那么 A 有贡献时子矩阵大小必然不超过 2。
考虑直接枚举这个子矩阵大小 k,情况如下:
$k=1$时,$det(D - A) = \sum_p(-1)^{inv(p)}(-A_{u,p_u})\prod_{i=1,i \ne u}^n D_{i,p_i} = (-2a_u)\prod_{i=1,i \ne u}^n D_{i,i}$。
$k=2$时,$det(D - A) = (-a_u^2 + 2a_ua_v - a_v^2)\prod_{i=1,i \ne u, i \ne v}^n D_{i,i}$。
于是我们学会了 $O(n^2)$ 计算行列式,非常高级。
考虑优化指数级复杂度的部分,即对于环的集合的枚举。令当前环上点集为 $S$,环内贡献为 $\sum_p \prod_{i=1}^{|S|} (a_{p_i} + a_{p_{i+1}})$,此处 $p$ 的含义是环上点的一种排列,$p_{i+1}$ 若超出范围则首尾相接。考虑到 $a_{p_i}$ 只会和与自己相邻的点同时出现,则最多出现两次。对式子直接展开,则 $a_{p_i}$ 最多出现二次幂,即式子中只会出现 $a_{p_i}^0$,$a_{p_i}^1$,$a_{p_i}^2$。
针对此性质进行优化。我们直接用展开的式子进行计算,即将 $\sum_p \prod_{i=1}^{|S|} (a_{p_i} + a_{p_{i+1}})$ 展开后若干个形式为 $\prod_{i=1}^{|S|}a_{p_i}^{t_i}$ 的结构拿来计算,其中 $t_i$ 的值为 $0$,$1$ 或者 $2$。换言之,我们直接把一个环拆做多份贡献。那么对于一个 $\prod_{i=1}^{|S|}a_{p_i}^{t_i}$ 形式的贡献,可能会有多种环的形态满足可以跑出这样的贡献,假设 $t_i=0$ 的 $i$ 有 $a$ 个,假设 $t_i=1$ 的 $i$ 有 $b$ 个,假设 $t_i=2$ 的 $i$ 有 $c$ 个,我们要求出能够满足以上三个条件的环的数量。
此时编号没有意义,我们只需要使用 $a$,$b$,$c$ 进行计算。回到原本环上 $\sum_p \prod_{i=1}^{|S|} (a_{p_i} + a_{p_{i+1}})$ 的模型,可以视作对环上每条边进行定向,最后每个点 $i$ 的入度即对应转化后点 $i$ 的指数 $t_i$。
问题转化为已知 $a$ 个点入度为 $0$,$b$ 个点入度为 $1$,$c$ 个点入度为 $2$,求可以构造出的环的数量。
显然环上的形态是度数为 $0$ 的点和度数为 $2$ 的点交替出现,中间间隔若干个度数为 $1$ 的点,则必须满足 $a=c$,否则方案数为 $0$。
先将度数为 $0$ 的点和度数为 $2$ 的点交替放,方案数为 $\frac{(a-1)!a!}{2}$。再将度数为 $1$ 的点随意插入,总方案数为 $\frac{(a-1)!a!(2a+b-1)!}{2(2a-1)!}$,但是特判 $a = 0$ 的情况,计算是简单的。
考虑打牌,定义状态 $dp_{i,j,a,b,c}$ 表示当前遍历到第 $i$ 个点,度数为 $0$,$1$,$2$ 的分别有 $a$,$b$,$c$ 个,基环树上环内 $A$ 矩阵没有贡献,贡献了一个 $j$ 或者贡献完了的方案数,还要考虑优化。
对于环内部分,状态中的 $j$ 实际的形态是 $(-a_u^2 + 2a_ua_v - a_v^2)$ 的。当枚举到 $u$ 时,实际可以将 $u$ 和 $v$ 独立开来直接计入答案而不是再开一维。即先加入 $a_u^2$ 的贡献,对于 $2a_u$ 的部分另开一倍空间记录下来同步更新,枚举到 $v$ 时再拿来计算并累加。时间复杂度变为 $O(n^4)$。
继续优化。回到 $D_{i,i} = \sum_{j=1}^n a_i+a_j = na_i + \sum_{j=1}^na_j$ 的定义上。我们在计算中多次使用了 $D_{i,i}$,环上的这一部分限制死了我们 $a$,$b$,$c$ 的状态是去不掉一个的。考虑将 $\sum_{j=1}^na_j$ 视作常数 $r$,则 $D_{i,i}=r+na_i$,即 $\prod_{i=1}^nD_{i,i} =\prod_{i=1}^n (r+na_i)$。可以发现我们能够用前面将环的贡献拆开的方法来拆这个式子的贡献,其中 $a_i$ 的幂次最大为 $1$。干脆将 $(-a_u^2 + 2a_ua_v - a_v^2)$ 这一部分也拆出来,无非就是 $a_i$ 的最大幂次变为了 $2$,似乎也可以用同种方法算。
尝试将环上这部分贡献直接并入 $a$,$b$,$c$ 中,这样我们又可以省下一维状态。
重新定义状态 $dp_{i,a,b,c}$ 表示已遍历前 $i$ 个点,其中部分可能在环内,部分可能在环外,贡献形式为 $a_i^0$ 的有 $a$ 个,贡献形式为 $a_i^1$ 的有 $b$ 个,贡献形式为 $a_i^2$ 的有 $c$ 个,此时权值相乘的值之和。由于此时 $a+b+c=i$,去除一维 $c$。显然不考虑附带的常数,将环内部分拆开乘上环外部分是没有问题的,不过这要求我们对于环内和环外的多种附带权值方案进行额外的计算。即计算一个数组可以表示已知 $a$,$b$,$c$ 时,将这三种贡献的形式分配至环内或是环外,所有这样的方案附带的权值之和。
定义 $f_{i,j,k}$ 表示 $i$ 个数贡献形式为 $a_i^0$,$j$ 个数贡献形式为 $a_i^1$,$k$ 个数贡献形式为 $a_i^2$,环内环外的所有系数。同理,$i + j + k = n$,去除一维 $k$。
具体的计算步骤上,先花 $O(n^2)$ 的复杂度枚举**有向环上**度为 $0$ 和 $1$ 的点数 $x$ 与 $y$,显然度为 $2$ 的点数 $z=x$;再枚举矩阵 $A$ 贡献了 $0$ 个,$1$ 个,还是 $2$ 个数;于是我们可以算出剩下在环外集合 $T$ 的大小以及具体的幂次分布(显然只有 $0$ 次和 $1$ 次),最后枚举集合 $T$ 中 $1$ 次幂的数量,即 $na_i$,就可以算出具体附带的权值。
$dp_{i,a,b}$ 可以滚动,最后的空间复杂度为 $O(n^2)$,时间复杂度为 $O(n^3)$,但是究极小常。
```cpp
#include <bits/stdc++.h>
using namespace std;
const int N = 1010;
const int Mod = 998244353;
int n, r;
int a[N], fact[N], infact[N];
int C[N][N], PowN[N], Powr[N];
int dp[2][N][N], f[N][N];
int Ad(int a, int b) {
return (a + b >= Mod ? a + b - Mod : a + b);
}
int qmi(int a, int b, int p) {
int res = 1;
while (b) {
if (b & 1) res = 1ll * res * a % p;
a = 1ll * a * a % p;
b >>= 1;
}
return res;
}
void init() {
C[0][0] = 1;
for (int i = 1; i < N; i ++) {
C[i][0] = 1;
for (int j = 1; j <= i; j ++) C[i][j] = Ad(C[i - 1][j], C[i - 1][j - 1]);
}
PowN[0] = 1;
for (int i = 1; i <= n; i ++) PowN[i] = 1ll * PowN[i - 1] * n % Mod;
Powr[0] = 1;
for (int i = 1; i <= n; i ++) Powr[i] = 1ll * Powr[i - 1] * r % Mod;
fact[0] = 1;
for (int i = 1; i <= n; i ++) fact[i] = 1ll * fact[i - 1] * i % Mod;
infact[n] = qmi(fact[n], Mod - 2, Mod);
for (int i = n; i >= 1; i --) infact[i - 1] = 1ll * infact[i] * i % Mod;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; i ++) scanf("%d", &a[i]), r = Ad(r, a[i]);
init();
dp[0][0][0] = 1;
for (int i = 1; i <= n; i ++) {
for (int A = 0; A < i; A ++) {
for (int B = 0; A + B < i; B ++) {
int F = dp[(i - 1) & 1][A][B], w0 = 1, w1 = a[i], w2 = 1ll * a[i] * a[i] % Mod;
dp[i & 1][A + 1][B] = Ad(dp[i & 1][A + 1][B], F * w0);
dp[i & 1][A][B + 1] = Ad(dp[i & 1][A][B + 1], 1ll * F * w1 % Mod);
dp[i & 1][A][B] = Ad(dp[i & 1][A][B], 1ll * F * w2 % Mod);
dp[(i - 1) & 1][A][B] = 0;
}
}
}
for (int i = 0; i <= n; i ++) {
for (int j = 0; i + j + i <= n; j ++) {
int Cir = i + j + i;
if (Cir < 3) continue;
for (int A = 0; A <= 3; A ++) {
int I = i, J = j, K = i, Pre = 1;
if (I == 0) Pre = fact[J - 1];
else Pre = 1ll * fact[I - 1] * fact[I] % Mod * fact[I * 2 + J - 1] % Mod * infact[2] % Mod * infact[2 * I - 1] % Mod;
if (A == 1) J ++, Pre = 1ll * Pre * (Mod - 2) % Mod;
if (A == 2) I ++, K ++, Pre = 1ll * Pre * (Mod - 1) % Mod;
if (A == 3) J += 2, Pre = 1ll * Pre * 2 % Mod;
if (I + J + K > n) continue;
for (int T = 0; T <= n - (I + J + K); T ++) {
int t0 = T, t1 = n - (I + J + K) - T;
int G = 1ll * Pre * Powr[t0] % Mod * PowN[t1] % Mod;
int C_ = 1ll * C[i + t0][i] * C[j + t1][j] % Mod;
if (A == 0) f[I + t0][J + t1] = Ad(f[I + t0][J + t1], 1ll * G * C_ % Mod);
if (A == 1) f[I + t0][J + t1] = Ad(f[I + t0][J + t1], 1ll * G * C_ % Mod * (J + t1) % Mod);
if (A == 2) f[I + t0][J + t1] = Ad(f[I + t0][J + t1], 1ll * G * C_ % Mod * (I + t0) % Mod * K % Mod);
if (A == 3) f[I + t0][J + t1] = Ad(f[I + t0][J + t1], 1ll * G * C_ % Mod * C[J + t1][2] % Mod);
}
}
}
}
int ans = 0;
for (int A = 0; A <= n; A ++)
for (int B = 0; A + B <= n; B ++)
ans = Ad(ans, 1ll * dp[n & 1][A][B] * f[A][B] % Mod);
printf("%d", ans);
return 0;
}
```