矩阵树定理 & BEST 定理

行列式

定义

对于 $n\times n$ 的方阵 $A$：

$$\begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \end{bmatrix}$$

定义其行列式 $\det(A) = \displaystyle \sum_{p}(-1)^{\pi(p)}\prod_{i=1}^n a_{i,p_i}$，其中 $p$ 是长度为 $n$ 的排列，$\pi(p)$ 表示 $p$ 的逆序对数。

记作 $\det(A)$，$|A|$ 或 $\begin{vmatrix}a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} \\\vdots & \vdots & \ddots & \vdots \\a_{n, 1} & a_{n, 2} & \cdots & a_{n, n}\end{vmatrix}$。

例如二阶行列式根据定义可得

$$\begin{vmatrix} a_{1,1} & a_{1,2}\\ a_{2,1} &a_{2,2} \end{vmatrix} = (-1)^0 a_{1,1}a_{2,2} + (-1)^1 a_{1,2}a_{2,1} = a_{1,1}a_{2,2}-a_{1,2}a_{2,1}$$

性质

1. 方阵转置后行列式不变。
   即 $\det(A) = \det\left(A^{\mathrm{T}}\right)$

   证明：

   $$\det(A) = \displaystyle \sum_{p}(-1)^{\pi(p)}\prod_{i = 1}^n a_{i,p_i}\\ \det\left(A^{\mathrm{T}}\right)=\displaystyle \sum_{p}(-1)^{\pi(p)}\prod_{i = 1}^n a_{p_i,i}$$

   令 $p'$ 为满足 $p'_{p_i}=i$ 的排列，只需证明 $\pi(p)=\pi(p')$ 即可。

   任意 $p$ 中的逆序对 $(i,j)$ 满足 $i<j$ 且 $p_i>p_j$。则在 $p'$ 中有 $p_i>p_j$ 且 $p'_{p_i}>p'_{p_j}$。

   因此 $p$ 中逆序对对应 $p'$ 的一个逆序对，同理 $p'$ 中逆序对对应 $p$ 的一个逆序对，故 $\pi(p)=\pi(p')$。

2. 行列式某一行的公因子可以提取到行列式外面。
   即：

   $$\begin{vmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ ka_{i, 1} & ka_{i, 2} & \cdots & ka_{i, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \end{vmatrix}= k \begin{vmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{i, 1} & a_{i, 2} & \cdots & a_{i, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \end{vmatrix}$$

   证明：

   $$\text{左边} = \displaystyle \sum_{p}(-1)^{\pi(p)}\prod_{i = 1}^n ka_{i,p_i} = \text{右边}\\$$

3. 将行列式中某一行可以拆成两个部分拆开求和，与原行列式相等。
   即：

   $$\begin{vmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ b_{i, 1}+c_{i, 1} & b_{i, 2}+c_{i, 2} & \cdots & b_{i, n}+c_{i, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \end{vmatrix}= \begin{vmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ b_{i, 1} & b_{i, 2} & \cdots & b_{i, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \end{vmatrix}+ \begin{vmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ c_{i, 1} & c_{i, 2} & \cdots & c_{i, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \end{vmatrix}$$

   证明：

   $$\begin{aligned} \text{左边}&= \displaystyle \sum_{p}(-1)^{\pi(p)}( b_{i,p_i}+c_{i,p_i})\prod_{1\le j \le n, j\ne i} a_{j,p_j}\\ &= \displaystyle \sum_{p}(-1)^{\pi(p)}b_{i,p_i}\prod_{1\le j \le n, j\ne i} a_{j,p_j} + \displaystyle \sum_{p}(-1)^{\pi(p)}c_{i,p_i}\prod_{1\le j \le n, j\ne i} a_{j,p_j}\\ &= \text{右边} \end{aligned}$$

4. 交换行列式两行，行列式值反号。

   证明：
   只需证明交换排列中两个数逆序对奇偶性改变即可。

   不难发现如果交换相邻元素排列逆序对变化 $\pm 1$，奇偶性一定改变。则若交换 $a_i$ 和 $a_j\ (i<j)$，相当于 $a_i$ 右移 $(j-i)$ 次，$a_j$ 左移 $(j-i-1)$ 次。共交换奇数次，奇偶性改变。

5. 若行列式中两行成比例，行列式为零。

   证明：
   先利用性质 2 将比例提出来，此时行列式中有两行相等。
   那么根据性质 4，交换这两行行列式的值反号，而行列式的值又不变，即 $\det(A)=-\det(A)$，故 $\det(A)=0$。

6. 将行列式某一行的 k 倍加到另一行上，行列式值不变。

   证明：
   先利用性质 3 把 $k$ 倍拆出来，根据性质 5 另外一部分为 $0$，剩下的就是原行列式。

7. 对行满足的性质对列同样满足。

   证明：
   利用性质 1，先转置再转置回来即可。

求法

根据定义求复杂度太高。事实上可以在 $O(n^3)$ 的时间内求出行列式。只需要将行列式化成上三角形式，行列式的值就等于主对角线元素乘积。

证明：
即证：

$$\begin{vmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ 0 & a_{2, 2} & \cdots & a_{2, n} \\ \vdots& \vdots & \ddots & \vdots \\ 0& 0 & \cdots& a_{n, n} \end{vmatrix} = \prod_{i=1}^na_{i,i}$$

根据定义，想要选出一个所有元素都非零的排列，则 $p_n=n$，$p_{n-1}$ 不能等于 $n$ 只能等于 $n-1$，依次类推得到只有 $p_i=i$ 的一个排列元素都非零，所以得证。

也可以对这个行列式进行代数余子式展开，得到的是一样的结果。

通过类似高斯消元的过程，运用性质 4 和性质 6，所有行列式都能化成上三角形式。
只需要不断消元，如果这一列全部都是零那么行列式就是零，可以直接返回。

（代码略，可参考下面任意模数行列式代码）

事实上，有些题目要取模。如果模数是质数有逆元，消元是简单的。
如果模数不是质数，我们也可以在 $O(n^3)$ 复杂度内求出行列式。

假设现在要消第 $i$ 行，有 $a_{i,i},a_{j,i}\ne 0$，想要让 $a_{i,i}\ne 0, a_{j,i}=0$。

对 $a_{i,i}$ 和 $a_{j,i}$ 进行辗转相除，即不断进行 $(a,b)\rightarrow(b\bmod a, a)$ 直至 $a=0$，最后一定能得到 $(0,k)$ 的形式，再交换回来就完成了对 $a_{j,i}$ 的消元。

$(a,b)\rightarrow(b\bmod a, a)$ 这个操作等价于对第 $j$ 行减去 $\lfloor b/a\rfloor$ 倍的第 $i$ 行，而这个减法是可以取模的，然后交换两行，第一个操作不改变行列式，第二个操作行列式变号。

记得要维护交换次数，来决定行列式的符号。

为了证明复杂度，我们需要知道，这个 $(a,b)\rightarrow(b\bmod a, a)$ 不会进行太多次。因为每次这么操作之后，相当于 $a_{i,i} = \gcd(a_{i,i},a_{j,i})$，而每次交换都至少会让 $a_{i,i}$ 减半，所以只会进行 $O(n^2 + n\log p)$ 次交换，

因此，复杂度是 $O(n^2 \log p)$ 的，不是瓶颈，复杂度就是 $O(n^3)$。

// 求 det(A) mod p
int det(vector<vector<int>> A,int p){
    int n = A.size();
    int res = 1; // 答案
    for(int i=0;i<n;i++){
        // 把第 i 列第 (i+1~n) 行全部消成 0
        for(int j=i+1;j<n;j++){
            // 把 a[i][i] 消掉
            while(A[i][i]){
                int div = A[j][i]/A[i][i]; // 辗转相除
                for(int k=i;k<n;k++) // 只会进行 O(n^2 + n log p) 次
                    A[j][k] = ((A[j][k] - (ll)div*A[i][k])%p+p)%p; 
                swap(A[i],A[j]); res=-res;
            }
            swap(A[i],A[j]); res=-res;
        }
    }
    for(int i=0;i<n;i++)res=(ll)res*A[i][i]%p;
    res=(res+p)%p;
    return res;
}

矩阵树定理

记号定义

参考自 OI-wiki。

以下均是不带权图，有重边、但无自环。

无向图

设 $G$ 是一个有 $n$ 个顶点的无向图。定义度数矩阵 $D(G)$ 为

$$D_{ii}(G) = \mathrm{deg}(i),\ D_{ij} = 0,\ i\neq j.$$

设 $\#e(i,j)$ 为点 $i$ 与点 $j$ 相连的边数，并定义邻接矩阵 $A$ 为

$$A_{ij}(G)=A_{ji}(G)=\#e(i,j),\ i\neq j.$$

定义 Laplace 矩阵（亦称 Kirchhoff 矩阵）$L$ 为

$$L(G) = D(G) - A(G).$$

记图 $G$ 的所有生成树个数为 $t(G)$。

有向图

设 $G$ 是一个有 $n$ 个顶点的有向图。定义出度矩阵 $D^{out}(G)$ 为

$$D^\mathrm{out}_{ii}(G) = \mathrm{deg}^\mathrm{out}(i),\ D^\mathrm{out}_{ij} = 0,\ i\neq j.$$

类似地定义入度矩阵 $D^\mathrm{in}(G)$。

设 $\#e(i,j)$ 为点 $i$ 指向点 $j$ 的有向边数，并定义邻接矩阵 $A$ 为

$$A_{ij}(G)=\#e(i,j),\ i\neq j.$$

定义出度 Laplace 矩阵 $L^\mathrm{out}$ 为

$$L^\mathrm{out}(G) = D^\mathrm{out}(G) - A(G).$$

定义入度 Laplace 矩阵 $L^\mathrm{in}$ 为

$$L^\mathrm{in}(G) = D^\mathrm{in}(G) - A(G).$$

记图 $G$ 的以 $k$ 为根的所有根向树形图个数为 $t^\mathrm{root}(G,k)$。所谓根向树形图，是说这张图的基图是一棵树，所有的边全部指向父亲。

记图 $G$ 的以 $k$ 为根的所有叶向树形图个数为 $t^\mathrm{leaf}(G,k)$。所谓叶向树形图，是说这张图的基图是一棵树，所有的边全部指向儿子。

带权有向图矩阵树定理证明

对于带权图的以上符号可以类推，改成边权求和即可。

定理叙述

对于任意的 $k$，都有

$$\sum_{T\in\mathcal T^\mathrm{root}(G,k)}w(T)=\det L^\mathrm{out}(G)_{[n]\setminus\{k\},[n]\setminus\{k\}}.$$

这里，$\mathcal T^\mathrm{root}(G,k)$ 是 $G$ 的以 $k$ 为根的根向树形图的集合，$w(T)$ 代表 $T$ 中所有边权的乘积。

Cauchy–Binet 公式

对于大小分别为 $n\times m$ 和 $m\times n$ 的矩阵 $A,B$，有：

$$\det(AB) = \sum_{1\le j_1 < j_2 < \ldots < j_n \le m} \det(A_{j_1,j_2,\ldots,j_n})\det(B_{j_1,j_2,\ldots,j_n}).$$

其中 $A_{j_1,j_2,\ldots,j_n}$ 表示包含 $A$ 第 $j_1,j_2,\ldots j_n$ 列的 $n\times n$ 矩阵，$B_{j_1,j_2,\ldots,j_n}$ 表示包含 $B$ 第 $j_1,j_2,\ldots j_n$ 行的 $n\times n$ 矩阵。

证明：

$$\begin{aligned} \det(AB) &= \sum_{p}(-1)^{\pi(p)}\prod_{i=1}^n (AB)_{i,p_i}\\ &=\sum_{p}(-1)^{\pi(p)}\prod_{i=1}^n \sum_{j=1}^m a_{i,j}b_{j,p_i}.\\ \end{aligned}$$

其中 $p$ 是长度为 $n$ 的排列，$\pi(p)$ 为 $p$ 逆序对个数。

考察上式最后的形式，展开后每一项都包含了 $\prod_{i=1}^n a_{i,k_i}$（$k_1,k_2,\ldots k_n$ 是每项均在 $[1,m]$ 的有序数组，可以重复），因此上式可以写作

$$\begin{aligned} &\sum_{1\le k_1,\ldots k_n \le m}\prod_{i=1}^n a_{i,k_i} \sum_{p}(-1)^{\pi(p)}\prod_{i=1}^n b_{k_i,p_i}\\ =&\sum_{1\le k_1,\ldots k_n \le m}\prod_{i=1}^n a_{i,k_i} \det(B_{k_1,k_2,\ldots k_n}).\\ \end{aligned}$$

令 $j_1,j_2,\ldots j_n$ 为 $k_1,k_2,\ldots k_n$ 排序后结果，不难得出 $\det(B_{k_1,k_2,\ldots k_n}) = (-1)^{\pi(k_1,k_2,\ldots,k_n)} \det(B_{j_1,j_2,\ldots,j_n}).$（考虑交换两行的次数）

那么原式

$$\begin{aligned} &=\sum_{1\le k_1,\ldots k_n \le m} (-1)^{\pi(k_1,k_2,\ldots,k_n)}\prod_{i=1}^n a_{i,k_i} \det(B_{j_1,j_2,\ldots j_n})\\ &=\sum_{1\le j_1\le j_2 \le \ldots \le j_n \le m}\det(B_{j_1,j_2,\ldots j_n})\sum_{k \text{是} j \text{重排的结果}} (-1)^{\pi(k_1,k_2,\ldots,k_n)}\prod_{i=1}^n a_{i,k_i} \\ &=\sum_{1\le j_1\le j_2 \le \ldots \le j_n \le m}\det(A_{j_1,j_2,\ldots j_n})\det(B_{j_1,j_2,\ldots j_n}). \end{aligned}$$

若有两个 $j$ 相同，那么

$$\det(A_{j_1,j_2,\ldots j_n}) = \det(B_{j_1,j_2,\ldots j_n}) = 0.$$

故有

$$\det(AB) = \sum_{1\le j_1< j_2 < \ldots < j_n \le m}\det(A_{j_1,j_2,\ldots j_n})\det(B_{j_1,j_2,\ldots j_n}).$$

图的关联矩阵与 Laplace 矩阵的关系

对于有向图 $G=(V,E)$，顶点数为 $n$（第 $i$ 个点记作 $v_i$），边数为 $m$（第 $i$ 条边记作 $e_i$）。由此，可以定义 $m\times n$ 阶出度关联矩阵

$$M^\mathrm{out}_{ij}=\begin{cases} \sqrt{w(e_i)},&e_i\text{是} v_j \text{的出边},\\ 0,&\textrm{otherwise}, \end{cases}$$

和 $m\times n$ 阶入度关联矩阵

$$M^\mathrm{in}_{ij}=\begin{cases} \sqrt{w(e_i)},&e_i\text{是} v_j \text{的入边},\\ 0,&\textrm{otherwise}. \end{cases}$$

则有 $A(G) = (M^\mathrm{out})^T M^\mathrm{in}$。

因为 $A(G)_{i,j} = \sum_{k=1}^m M^\mathrm{out}_{k,i}M^\mathrm{in}_{k,j} = \sum_{k=1}^m w(e_k)[e_k \text{是} i \text{到} j \text{的边}] .$

有 $D^\mathrm{in}(G) = (M^\mathrm{in})^T M^\mathrm{in}$。

因为 $D^\mathrm{in}(G) _{i,i} = \sum_{k=1}^m M^\mathrm{in}_{k,i}M^\mathrm{in}_{k,i} = \sum_{k=1}^m w(e_k)[e_k \text{是} i \text{的入边}].$
而 $D^\mathrm{in}(G) _{i,j} = \sum_{k=1}^m M^\mathrm{in}_{k,i}M^\mathrm{in}_{k,j}=0\ (i\ne j).$

故有 $L^\mathrm{in}(G) = (M^\mathrm{in}-M^\mathrm{out})^T M^\mathrm{in}$。

图的关联矩阵与图性质的关系

我们把关联矩阵和 Laplace 矩阵已经联系了起来，接下来我们将关联矩阵和图的结构联系在一起。

有如下引理：

对于 $G$ 的一个子图 $(W,S)$，若它满足 $|W|=|S|\le n$，则子图 $T=(V,S)$ 是一个以 $V\setminus W$ 为根的叶向森林，当且仅当

$$\det(M^\mathrm{in}_{S,W})\det(M^\mathrm{in}_{S,W}-M^\mathrm{ou}_{S,W}) = w(T).$$

其中 $w(T)$ 代表 $T$ 中所有边权的乘积。

证明：

把行列式里面每一行都提出一个 $\sqrt{w(e_i)}$，就能得到 $\det(M'^\mathrm{in}_{S,W})\det(M'^\mathrm{in}_{S,W}-M'^\mathrm{out}_{S,W}) = 1$，$M'^\mathrm{in}$ 和 $M'^\mathrm{out}$ 中都只有 $0$ 和 $1$。

易有 $\det(M'^\mathrm{in}_{S,W})\det(M'^\mathrm{in}_{S,W}-M'^\mathrm{out}_{S,W}) = 1$ 等价于两项均为 $1$，我们分别分析两项为 $1$ 能代表什么性质。

1. $\det(M'^\mathrm{in}_{S,W})=1.$
   首先不存在全 $0$ 的行，且不存在两行相同，所以 $W$ 中每个点都恰是 $S$ 中一条边的终点，且互不相同。

2. $\det(M'^\mathrm{in}_{S,W}-M^\mathrm{out}_{S,W}) = 1.$
   该式等价于图中不存在环。

   • （$\Rightarrow$）用反证法，假设存在环 $e_1 \to e_2 \to \ldots \to e_k$，将 $M'^\mathrm{in}_{S,W}-M'^\mathrm{out}_{S,W}$ 第 $e_1 ,e_2 , \ldots , e_k$ 行加起来，每一列都是 $0$，所以 $M'^\mathrm{in}_{S,W}-M'^\mathrm{out}_{S,W}$ 不满秩，行列式为 $0$，矛盾！
   • （$\Leftarrow$）如果没有环，总存在 $W$ 中一个点与 $V\setminus W$ 相连，这一列恰有一个 $1$，因此可以把其他的 $-1$ 都消去，递归处理总能消去所有的 $-1$，而这些 $+1$ 的位置就是各行对应边的终点，就是 $M'^\mathrm{in}_{S,W}$。

现在只需要证明，满足 $W$ 中每个点都恰是 $S$ 中一条边的终点，且互不相同的无环图是以 $V\setminus W$ 为根的叶向森林。

对于在 $W$ 中任意一点 $u$，沿着它唯一入边，总会跳到 $V\setminus W$ 中的节点（因为没有环），而 $V\setminus W$ 中的节点一定没有入边（因为 $|W|=|S|$），故得证。

最终证明

证明：
记 $W=[n]\setminus\{k\}$ 为除去 $k$ 点外的剩余顶点的集合。

根据 Cauchy–Binet 公式

$$\det L^\mathrm{in}(G)_{W,W} = \sum_{S\subset[m];~|S|=n-1}\det(M^\mathrm{in}_{S,W})\det(M^\mathrm{in}_{S,W}-M^\mathrm{out}_{S,W}).$$

遍历所有的 $S$，由引理 2，当且仅当 $T=(V,S)$ 构成一个以 $V\setminus W=\{k\}$ 为根的叶向森林时，亦即 $T$ 是一个以 $k$ 为根的叶向树形图时，右侧累加 $w(T)$。

故得证。

对于无权图只需要把边权看作 $1$ 即可。

矩阵树定理无向图形式

把无向图每个边都拆成两条，那么任意一个无向图的生成树都与以某个点为根的叶向树一一对应，因此得到了矩阵树定理的无向图形式。

$$\sum_{T\in\mathcal T(G)}w(T) = \det L(G)_{[n]\setminus\{k\},[n]\setminus\{k\}}.$$

这里，$\mathcal T(G)$ 是 $G$ 的生成树的集合。这也说明，$L(G)$ 的所有 $(n-1)$ 阶主子式都相等。

总结

Laplace 矩阵 $L=D-A$。

$D_{i,i}$ 代表的是这个点的度数：

• 无向图：与 $i$ 相连边的权值和；
• 有向图叶向树：入边权值和；
• 有向图根向树：出边权值和；

$A_{i,j}$ 表示 $i$ 和 $j$ 之间边数（$i\ne j$）：

• 无向图：$(i,j)$ 权值和；
• 有向图：$i\to j$ 权值和；

而 $L$ 的 $(n-1)$ 主子式代表了生成树的个数：

• 无向图：所有的 $(n-1)$ 主子式都相等，就表示生成树边权乘积和。
• 有向图：删去第 $k$ 行第 $k$ 列的主子式代表的就是以 $k$ 为根的叶向/根向树边权乘积和。

例题

P6178 【模板】Matrix-Tree 定理

模板题。

---

代码：

#include<bits/stdc++.h>
using namespace std;

using ll=long long;

const int MOD = 1e9+7;
// 求 det(A) mod p
int det(vector<vector<int>> A,int p){
    int n = A.size();
    int res = 1; // 答案
    for(int i=0;i<n;i++){
        for(int j=i+1;j<n;j++){
            while(A[i][i]){
                int div = A[j][i]/A[i][i]; 
                for(int k=i;k<n;k++)
                    A[j][k] = ((A[j][k] - (ll)div*A[i][k])%p+p)%p; 
                swap(A[i],A[j]); res=-res;
            }
            swap(A[i],A[j]); res=-res;
        }
    }
    for(int i=0;i<n;i++)res=(ll)res*A[i][i]%p;
    res=(res+p)%p;
    return res;
}
int main(){
    ios::sync_with_stdio(0);cin.tie(0);
    int n,m,t; cin>>n>>m>>t;
    vector<vector<int>> a(n-1, vector<int>(n-1));
    for(int i=0;i<m;i++){
        int u,v,w; cin>>u>>v>>w;
        u-=2, v-=2;
        if(v>=0)   a[v][v]=(a[v][v]+w)%MOD;
        if(u>=0)   if(!t)a[u][u]=(a[u][u]+w)%MOD;
        if(u>=0 && v>=0){
            a[u][v]=(a[u][v]-w+MOD)%MOD;
            if(!t)a[v][u]=(a[v][u]-w+MOD)%MOD;
        }
    }

    cout << det(a, MOD) << '\n';
    return 0;
}

P3317 [SDOI2014] 重建

求

$$\sum_{T\in\mathcal T(G)} \left(\prod_{e \in T}p_e\right)\left(\prod_{e \not\in T}(1-p_e)\right).$$

其中 $\mathcal T(G)$ 是 $G$ 的生成树的集合。

---

把上式改写成

$$\left(\prod_{e=1}^m (1-p_e)\right) \sum_{T\in\mathcal T(G)} \left(\prod_{e \in T}\dfrac{p_e}{1-{p_e}}\right)$$

即可。

这道题这里需要加一点 eps 扰动才能通过，但可能会被 Hack。

事实上，把 $p_e=1$ 的边提前加进图中缩点就能更好的解决本题。

也就是说，无向图钦定包含某些边的生成树计数也是一个可做问题。

代码：

#include<bits/stdc++.h>
using namespace std;

const double eps = 1e-8;
double det(vector<vector<double>> A){
    int n = A.size();
    double res = 1;
    for(int i=0;i<n;i++){
        int mx=i;
        for(int j=i+1;j<n;j++)if(abs(A[j][i])>abs(A[mx][i]))mx=j;
        if(abs(A[mx][i])<eps){
            res=0;
            break;
        }
        if(mx!=i)swap(A[i],A[mx]),res=-res;
        res *= A[i][i];
        for(int k=i+1;k<n;k++){
            double mul=A[k][i]/A[i][i];
            for(int j=i;j<n;j++)A[k][j]-=A[i][j]*mul;
        }
    }
    return res;
}

int main(){
    ios::sync_with_stdio(0);cin.tie(0);
    cout<<fixed;

    int n; cin>>n;
    vector<vector<double>> a(n-1, vector<double>(n-1));
    double ans = 1;
    for(int i=-1;i<n-1;i++){
        for(int j=-1;j<n-1;j++){
            double p; cin>>p;
            if(i>=j)continue;
            if(abs(p)<eps)p=eps;
            if(abs(p-1)<eps)p=1-eps;
            ans *= (1-p);
            double w = p/(1-p);
            if(i>=0)a[i][i]+=w;
            if(j>=0)a[j][j]+=w;
            if(i>=0 && j>=0)a[i][j]-=w,a[j][i]-=w;
        }
    }
    cout << det(a) * ans << '\n';
    return 0;
}

P5296 [北京省选集训2019] 生成树计数

求所有生成树权值和的 $k$ 次方之和。

---

我们可以发现，在矩阵树定理中，我们在证明时用到了边权的加法，乘法运算。
所以边权并不一定需要是实数，甚至可以是向量、矩阵，生成函数，只要能够满足这些运算就可以了。

众所周知等比数列 $[1,p,p^2,\ldots,p^n,\ldots]$ 的 egf 是 $e^{px}$，而 $e^{ax} \cdot e^{bx} = e^{(a+b)x}$，我们就把乘法变成了加法，最后只要求出所有边 egf 乘积第 $x^k$ 项系数乘上 $k!$ 即可。

需要多项式求逆，多项式乘法，都用暴力 $O(k^2)$ 实现。

复杂度 $O(n^3 k^2)$。

#include<bits/stdc++.h>
using namespace std;

template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <unsigned umod>
struct modint {
    static constexpr int mod = umod;
    unsigned v;
    modint() : v(0) {}
template <class T, must_int<T> = 0>
     modint(T _v) {int x = _v % (int)umod; v = x < 0 ? x + umod : x;}
     modint operator+() const { return *this; }
     modint operator-() const { return modint() - *this; }
     friend int raw(const modint &self) { return self.v; }
     friend ostream &operator<<(ostream &os, const modint &self) { return os << raw(self);}
     modint &operator+=(const modint &rhs) {v += rhs.v;if (v >= umod) v -= umod;return *this;}
     modint &operator-=(const modint &rhs) {v -= rhs.v;if (v >= umod) v += umod;return *this;}
     modint &operator*=(const modint &rhs) {v = 1ull * v * rhs.v % umod; return *this;}
     modint &operator/=(const modint &rhs) { return *this *= rhs.inv(); }
     modint inv() const {
         assert(v);
         static unsigned lim = 1 << 21;
         static vector<modint> inv{0, 1};
         if (v >= lim) return qpow(*this, mod - 2);
         inv.reserve(v + 1);
         while (v >= inv.size()) {
             int m = inv.size();
             inv.resize(m << 1);
             for (int i = m; i < m << 1; i++)inv[i] = (mod - mod / i) * inv[mod % i];
         }
         return inv[v];
     }
    template <class T, must_int<T> = 0>
    friend modint qpow(modint a, T b) {modint r = 1;for (; b; b >>= 1, a *= a)if (b & 1) r *= a;return r;}
    friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
    friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
    friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
    friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
    friend bool operator==(const modint &lhs, const modint &rhs) { return lhs.v == rhs.v; }
    friend bool operator!=(const modint &lhs, const modint &rhs) { return lhs.v != rhs.v; }
};
const int MOD = 998244353;
using mint = modint<MOD>;

const int K = 31;

int k;
using Poly = vector<mint>;
mint f[K],finv[K];

Poly exp_poly(mint w){ // e^{wx} 的 egf
    Poly res(k);
    mint now = 1;
    for(int i=0;i<k;i++)res[i]=finv[i]*now,now*=w;
    return res; 
}

Poly operator-(const Poly& a){
    Poly c(k);
    for(int i=0;i<k;i++)c[i]=-a[i];
    return c;
}
Poly operator+(const Poly& a,const Poly& b){
    Poly c(k);
    for(int i=0;i<k;i++)c[i]=a[i]+b[i];
    return c;
}
Poly operator-(const Poly& a,const Poly& b){
    Poly c(k);
    for(int i=0;i<k;i++)c[i]=a[i]-b[i];
    return c;
}
Poly operator*(const Poly& a,const Poly& b){
    Poly c(k);
    for(int i=0;i<k;i++){
        for(int j=0;j<=i;j++)c[i]+=a[j]*b[i-j];
    }
    return c;
}

Poly inv(const Poly& a){
    assert(a[0]!=0);
    Poly res(k);
    res[0] = 1/a[0];
    for(int i=1;i<k;i++){
        for(int j=0;j<i;j++)res[i]-=res[j]*a[i-j];
        res[i]*=res[0];
    }
    return res;
}

Poly det(vector<vector<Poly>> A){
    int n = A.size();
    Poly res(k);res[0]=1;
    for(int i=0;i<n;i++){
        int mx=-1;
        for(int j=i;j<n;j++)if(A[j][i][0]!=0){mx=j;break;}
        if(mx==-1)assert(0);
        if(mx!=i)swap(A[i],A[mx]),res=-res;
        
        res = res * A[i][i];
        auto invi = inv(A[i][i]);
        for(int q=i+1;q<n;q++){
            auto mul = A[q][i] * invi;
            for(int j=i;j<n;j++)A[q][j] = A[q][j] - A[i][j]*mul;
        }
    }
    return res;
}

int main(){
    ios::sync_with_stdio(0);cin.tie(0);
    int n;cin >> n >> k; ++k;
    f[0]=finv[0]=1;
    for(int i=1;i<k;i++)f[i]=f[i-1]*i, finv[i]=finv[i-1]/i;


    vector<vector<Poly>> a(n-1, vector<Poly>(n-1, Poly(k)));
    for(int i=-1;i<n-1;i++){
        for(int j=-1;j<n-1;j++){
            int w; cin>>w;
            if(i>=j)continue;
            Poly v = exp_poly(w);
            if(i>=0)a[i][i]=a[i][i]+v;
            if(j>=0)a[j][j]=a[j][j]+v;
            if(i>=0 && j>=0)a[i][j]=a[i][j]-v,a[j][i]=a[j][i]-v;
        }
    }

    cout << det(a)[k-1]*f[k-1] << '\n';
    return 0;
}

这告诉我们所有生成树边权和的总和也是可以在 $O(n^3)$ 内求出的。

练习

• P4111 [HEOI2015] 小 Z 的房间
• P2144 [FJOI2007] 轮状病毒
• P6624 [省选联考 2020 A 卷] 作业题

BEST 定理

定理内容

设 $G$ 是有向欧拉图，$k$ 为任意顶点，那么 $G$ 的不同欧拉回路总数 $\mathrm{ec}(G)$（循环同构不重复计算）是

由于欧拉图入度等于出度，下文的 $\deg$ 无需标明是入度还是出度。

$$\mathrm{ec}(G) = t^\mathrm{root}(G,k)\prod_{v\in V}(\deg (v) - 1)!.$$

推论：
由定理立刻得知：

1. $G$ 从 $s$ 出发的欧拉回路总数是

   $$\mathrm{ec}_s(G) = t^\mathrm{root}(G,k)(\deg(s))!\prod_{v\ne s}(\deg (v) - 1)!.$$

2. $G$ 中欧拉路径总数是

   $$\mathrm{ec}_s(G) = |E| t^\mathrm{root}(G,k)\prod_{v\in V}(\deg (v) - 1)!.$$

注意定理使用条件：

1. 是有向图。
2. 有欧拉回路。

并注意求的是根向树个数。

证明：
见 OI-wiki。

模板

P5807 【模板】BEST 定理 | Which Dreamed It

给定一个有向图，求从 $1$ 开始又回到 $1$ 的欧拉回路的个数。

---

先要保证原图是欧拉图。

用并查集，确保原图是一个大的连通分量加上一些孤立点（且没有自环），然后大的连通分量内每个点入度等于出度。

然后再用矩阵树定理，注意对于孤立点的那些行和列应该跳过，不参与行列式的计算。

#include<bits/stdc++.h>
using namespace std;


template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <unsigned umod>
struct modint {
    static constexpr int mod = umod;
    unsigned v;
    modint() : v(0) {}
template <class T, must_int<T> = 0>
     modint(T _v) {int x = _v % (int)umod; v = x < 0 ? x + umod : x;}
     modint operator+() const { return *this; }
     modint operator-() const { return modint() - *this; }
     friend int raw(const modint &self) { return self.v; }
     friend ostream &operator<<(ostream &os, const modint &self) { return os << raw(self);}
     modint &operator+=(const modint &rhs) {v += rhs.v;if (v >= umod) v -= umod;return *this;}
     modint &operator-=(const modint &rhs) {v -= rhs.v;if (v >= umod) v += umod;return *this;}
     modint &operator*=(const modint &rhs) {v = 1ull * v * rhs.v % umod; return *this;}
     modint &operator/=(const modint &rhs) { return *this *= rhs.inv(); }
     modint inv() const {
         assert(v);
         static unsigned lim = 1 << 21;
         static vector<modint> inv{0, 1};
         if (v >= lim) return qpow(*this, mod - 2);
         inv.reserve(v + 1);
         while (v >= inv.size()) {
             int m = inv.size();
             inv.resize(m << 1);
             for (int i = m; i < m << 1; i++)inv[i] = (mod - mod / i) * inv[mod % i];
         }
         return inv[v];
     }
    template <class T, must_int<T> = 0>
    friend modint qpow(modint a, T b) {modint r = 1;for (; b; b >>= 1, a *= a)if (b & 1) r *= a;return r;}
    friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
    friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
    friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
    friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
    friend bool operator==(const modint &lhs, const modint &rhs) { return lhs.v == rhs.v; }
    friend bool operator!=(const modint &lhs, const modint &rhs) { return lhs.v != rhs.v; }
};
const int MOD = 1e6+3;
using mint = modint<MOD>;


const int N = 105;


int fa[N],sz[N];
int find(int x){
    return fa[x]==x?x:fa[x]=find(fa[x]);
}
void unite(int i,int j){
    i=find(i),j=find(j);
    if(i==j)return;
    if(sz[i] > sz[j])swap(i,j);
    fa[i] = j;
    sz[j] += sz[i], sz[i]=0;
}

int in[N],out[N];
bool loop[N];

mint A[N][N];

int n;

mint det(){
    mint res = 1;
    for(int i=2;i<=n;i++){
        if(find(i) != find(1))continue;
        int mx=i;
        for(int j=i;j<=n;j++){
            if(find(j)!=find(1))continue;
            if(A[j][i]!=0){mx=j;break;}
        }
        if(A[mx][i] == 0){
            res=0;break;
        }
        if(mx!=i)swap(A[i],A[mx]),res=-res;
        res *= A[i][i];
        for(int k=i+1;k<=n;k++){
            if(find(k)!=find(1))continue;
            auto mul=A[k][i]/A[i][i];
            for(int j=i;j<=n;j++)A[k][j]-=A[i][j]*mul;
        }
    }
    return res;
}

mint f[MOD];
void solve(){
    cin>>n;
    for(int i=1;i<=n;i++)fa[i]=i,sz[i]=1,in[i]=out[i]=0,loop[i]=0;
    for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=0;


    for(int u=1;u<=n;u++){
        int cnt; cin >> cnt;
        while(cnt--){
            int v; cin>>v;
            out[u]++,in[v]++;
            A[u][v]-=1;A[u][u]+=1;
            unite(u,v);
            if(u==v)loop[u]=1;
        }
    }

    bool ok = 1;
    for(int i=1;i<=n;i++){
        if(in[i] != out[i]){
            ok = 0; break;
        }
        if(find(i) != find(1) && (sz[find(i)]>1 || loop[i])){
            ok = 0; break;
        }
    }

    if(!ok){
        cout << 0 << '\n';
        return;
    }

    mint ans = det() * f[out[1]];
    for(int i=2;i<=n;i++){
        if(out[i])ans*=f[out[i]-1];
    }

    cout << ans << '\n';
}


int main(){
    ios::sync_with_stdio(0);cin.tie(0);
    cout<<fixed;

    f[0]=1;
    for(int i=1;i<MOD;i++)f[i]=f[i-1]*i;

    int T; cin>>T;
    while(T--)solve();
    return 0;
}

求欧拉路径个数

在 $O(n^3)$ 时间内也能求出图的欧拉路径个数。

当图里有欧拉回路的时候，处理方法与上面相同，求出欧拉回路个数后乘以总边数就是欧拉路径的个数了。
当图里有欧拉通路的时候，因为起点和终点确定，在图中加入一条从终点到起点的辅助边，那么在新图中每个欧拉回路都经过额外加的边，从这里断开就对应一个欧拉通路。反之同理。因此原图的欧拉通路数与新图的欧拉回路数（循环同构不重复计算）相同。

例题：ABC336G - 16 Integers

给出 $16$ 个非负整数 $X_{i,j,k,l}\ (i,j,k,l\in\{0,1\})$，令 $N = \displaystyle\sum_{i,j,k,l\in\{0,1\}}X_{i,j,k,l}$。

求长度为 $N+3$ 的 $\{0,1\}$ 合法序列 $A_{1},A_{2},\ldots,A_{n+3}$ 的个数，使得对于任意 $(i,j,k,l)$，满足 $A_s=i,A_{s+1}=j,A_{s+2}=k,A_{s+3}=l$ 的 $s$ 恰有 $X_{i,j,k,l}$ 个。

$N\le 10^6$。

---

先考虑把原问题转化为欧拉路径计数。

建立 $8$ 个点 $(i,j,k)$，$A$ 中每个四元组 $(i,j,k,l)$ 都相当于在图上走了 $(i,j,k)\to (j,k,l)$ 的一条边，问题就转化成给定一张图，求经过每条边固定次数的路径个数。

把要经过的边建出来，就是欧拉路径计数。

上文已经叙述过如何求欧拉路径个数了。在求完欧拉路径个数时，还需注意本题的边是不区分的，所以答案还要除掉这些，也就是要除掉 $\prod (A_{i,j,k,l})!$。

复杂度 $O(B^3 + N)$，其中 $B=8$。

#include<bits/stdc++.h>
using namespace std;

template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <unsigned umod>
struct modint {
    static constexpr int mod = umod;
    unsigned v;
    modint() : v(0) {}
template <class T, must_int<T> = 0>
     modint(T _v) {int x = _v % (int)umod; v = x < 0 ? x + umod : x;}
     modint operator+() const { return *this; }
     modint operator-() const { return modint() - *this; }
     friend int raw(const modint &self) { return self.v; }
     friend ostream &operator<<(ostream &os, const modint &self) { return os << raw(self);}
     modint &operator+=(const modint &rhs) {v += rhs.v;if (v >= umod) v -= umod;return *this;}
     modint &operator-=(const modint &rhs) {v -= rhs.v;if (v >= umod) v += umod;return *this;}
     modint &operator*=(const modint &rhs) {v = 1ull * v * rhs.v % umod; return *this;}
     modint &operator/=(const modint &rhs) { return *this *= rhs.inv(); }
     modint inv() const {
         assert(v);
         static unsigned lim = 1 << 21;
         static vector<modint> inv{0, 1};
         if (v >= lim) return qpow(*this, mod - 2);
         inv.reserve(v + 1);
         while (v >= inv.size()) {
             int m = inv.size();
             inv.resize(m << 1);
             for (int i = m; i < m << 1; i++)inv[i] = (mod - mod / i) * inv[mod % i];
         }
         return inv[v];
     }
    template <class T, must_int<T> = 0>
    friend modint qpow(modint a, T b) {modint r = 1;for (; b; b >>= 1, a *= a)if (b & 1) r *= a;return r;}
    friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
    friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
    friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
    friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
    friend bool operator==(const modint &lhs, const modint &rhs) { return lhs.v == rhs.v; }
    friend bool operator!=(const modint &lhs, const modint &rhs) { return lhs.v != rhs.v; }
};
const int MOD = 998244353;
using mint = modint<MOD>;

int G[8][8];
int in[8],out[8];
int fa[8];
int find(int i){
    return fa[i]==i?i:fa[i]=find(fa[i]);
}
void unite(int i,int j){
    i=find(i),j=find(j);
    fa[i]=j;
}

void no_sol(){
    cout << 0 << '\n';
    exit(0);
}

bool vis[8];
mint A[8][8];
mint det(int skip){
    int n = 8;
    mint res = 1;
    for(int i=0;i<n;i++){
        if(!vis[i] || i==skip)continue;
        int mx=i;
        for(int j=i;j<n;j++){
            if(!vis[j] || j==skip)continue;
            if(A[j][i]!=0){mx=j;break;}
        }
        if(A[mx][i] == 0){res=0;break;}
        if(mx!=i)swap(A[i],A[mx]),res=-res;
        res *= A[i][i];
        for(int k=i+1;k<n;k++){
            if(!vis[k] || k==skip)continue;
            auto mul=A[k][i]/A[i][i];
            for(int j=i;j<n;j++)A[k][j]-=A[i][j]*mul;
        }
    }
    return res;
}

const int N = 1e6+5;
mint f[N],invf[N];
int main(){
    ios::sync_with_stdio(0);cin.tie(0);
    
    f[0]=invf[0]=1;
    for(int i=1;i<N;i++)f[i]=f[i-1]*i,invf[i]=invf[i-1]/i;
    for(int i=0;i<8;i++)fa[i]=i;
    int m = 0;
    for(int i=0;i<2;i++)for(int j=0;j<2;j++)for(int k=0;k<2;k++){
        for(int l=0;l<2;l++){
            int cnt; cin>>cnt;
            m += cnt;
            int u = (i<<2)|(j<<1)|k;
            int v = (j<<2)|(k<<1)|l;
            G[u][v] += cnt;
            out[u]+=cnt, in[v]+=cnt;
            unite(u,v);
        }
    }

    if(m == 0){
        cout << 8 << '\n';
        return 0;
    }

    int p0=0;
    while(!in[p0] && !out[p0])++p0;  
    int st = -1, ed = -1;

    for(int i=0;i<8;i++){
        if(!in[i] && !out[i])continue;
        if(find(i) != find(p0))no_sol();
        vis[i] = 1;
        if(in[i] == out[i])continue;
        else if(out[i] == in[i] + 1)
            if(st != -1)no_sol();
            else st = i;
        else if(in[i] == out[i] + 1)
            if(ed != -1)no_sol();
            else ed = i;
        else no_sol();
    }

    if((st==-1)^(ed==-1))no_sol();

    if(st!=-1){
        for(int i=0;i<8;i++)for(int j=0;j<8;j++){
            if(i==ed && j==st)++G[i][j],++out[i],++in[j]; // 加一条边
            A[i][j] -= G[i][j], A[i][i] += G[i][j];
        }

        auto ans = det(st);
        for(int i=0;i<8;i++){
            if(in[i])ans*=f[in[i]-1];
        }
        for(int i=0;i<8;i++){
            for(int j=0;j<8;j++){
                if(i==ed && j==st)--G[i][j]; // 加的那条边是要区分的
                ans *= invf[G[i][j]]; // 其他边不区分
            }
        }
        cout << ans << '\n';
    }
    else{
        for(int i=0;i<8;i++){
            for(int j=0;j<8;j++){
                A[i][j] -= G[i][j];
                A[i][i] += G[i][j];
            }
        }
        auto ans = det(p0);
        for(int i=0;i<8;i++){
            if(in[i])ans*=f[in[i]-1];
        }
        for(int i=0;i<8;i++){
            for(int j=0;j<8;j++){
                ans *= invf[G[i][j]];
            }
        }
        ans = ans * m;
        cout << ans << '\n';
    }
    return 0;
}

利用 BEST 定理转化后直接求解

用 BEST 定理时，不一定要搭配矩阵树定理使用。

例题：CF1919E Counting Prefixes

给定长度为 $n$ 的数列 $p_1,p_2,\ldots,p_n$，求有多少个长度为 $n$ 的数列 $a_1,a_2,\ldots,a_n$，满足：

• $\forall i \in [1,n], |a_i|=1$；
• $a_i$ 的前缀和数组排序后恰为 $p$。

对 $998244353$ 取模，$\sum n \le 5000$。

---

就是给出 $a_i$ 前缀和数组中每个数出现次数。

问题转化成，在形如下图的双向链上，从 $1$ 或 $-1$ 开始，满足每个点经过次数与给定数组相同的路径有多少种。

因为这个图是链状的，可以转化为欧拉路径。

我们先讨论在 $s\rightsquigarrow t$ 的欧拉路径中，每个点经过的次数：

1. 起点 $s$：
   $\textrm{deg}^+(s) = \textrm{deg}^-(s) + 1$，经过次数为 $\textrm{det}^+(s)$。

2. 终点 $t$：
   $\textrm{deg}^-(t) = \textrm{deg}^+(t) + 1$，经过次数为 $\textrm{det}^-(t)$。

3. 其他点 $u$：
   设点 $u$ 的度数为 $\textrm{deg}(u)$，那么 $u$ 经过的次数就是 $\textrm{deg}(u)$。

特别的，若 $s=t$，则起点（也是终点）经过的次数是 $\textrm{deg}(s)+1$。

那么，如果确定起点和终点，每个点的入度和出度都已经确定。又因为这个图是链状的，依次递推就能求出每条边经过次数。

现在要求出一条链（加一条边）中以某个点为根的根向树个数。这是容易在 $O(n)$ 时间复杂度内解决的，只需要选额外那条边的起点求根向树个数，额外的那条边就对答案无影响。

枚举起点和终点也是 $O(n)$ 的，故复杂度为 $O(n^2)$。

更详细的实现见代码注释。

#include<bits/stdc++.h>
using namespace std;


template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <unsigned umod>
struct modint {
    static constexpr int mod = umod;
    unsigned v;
    modint() : v(0) {}
template <class T, must_int<T> = 0>
     modint(T _v) {int x = _v % (int)umod; v = x < 0 ? x + umod : x;}
     modint operator+() const { return *this; }
     modint operator-() const { return modint() - *this; }
     friend int raw(const modint &self) { return self.v; }
     friend ostream &operator<<(ostream &os, const modint &self) { return os << raw(self);}
     modint &operator+=(const modint &rhs) {v += rhs.v;if (v >= umod) v -= umod;return *this;}
     modint &operator-=(const modint &rhs) {v -= rhs.v;if (v >= umod) v += umod;return *this;}
     modint &operator*=(const modint &rhs) {v = 1ull * v * rhs.v % umod; return *this;}
     modint &operator/=(const modint &rhs) { return *this *= rhs.inv(); }
     modint inv() const {
         assert(v);
         static unsigned lim = 1 << 21;
         static vector<modint> inv{0, 1};
         if (v >= lim) return qpow(*this, mod - 2);
         inv.reserve(v + 1);
         while (v >= inv.size()) {
             int m = inv.size();
             inv.resize(m << 1);
             for (int i = m; i < m << 1; i++)inv[i] = (mod - mod / i) * inv[mod % i];
         }
         return inv[v];
     }
    template <class T, must_int<T> = 0>
    friend modint qpow(modint a, T b) {modint r = 1;for (; b; b >>= 1, a *= a)if (b & 1) r *= a;return r;}
    friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
    friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
    friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
    friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
    friend bool operator==(const modint &lhs, const modint &rhs) { return lhs.v == rhs.v; }
    friend bool operator!=(const modint &lhs, const modint &rhs) { return lhs.v != rhs.v; }
};
const int MOD = 998244353;
using mint = modint<MOD>;


const int N = 5e4;
mint f[N],invf[N];

void solve(){
    int n;cin>>n;
    int B = n+3; // 因为有负数，所有下标加偏移量
    int pt  = 0; // 有效点数
    vector<int> cnt(B*2+1);

    // L,R 是有 cnt 的范围
    int L =  B-1, R = B+1;

    for(int i=1;i<=n;i++){
        int x; cin>>x;
        if(!cnt[x+B])++pt;
        cnt[x+B]++;
        L = min(L,x+B), R = max(R, x+B);
    }
    // 多往外拓展一个，用于判断无解的边（否则可能会有边往边界外指）
    --L,++R;

    if(n==1){
        if(cnt[B+1]||cnt[B-1])cout<<1<<'\n';
        else cout<<0<<'\n';
        return;
    }

    // 可以发现，在最终的图中，每个点的出度就是它的经过次数
    // 因此，我们提前计算出BEST定理中的度数阶乘项，减小常数
    mint res = 1;
    for(int i=L;i<=R;i++)if(cnt[i])res*=f[cnt[i]-1];

    vector<int> out(B*2+1);
    vector<int> cnt_L(B*2+1),cnt_R(B*2+1);

    mint ans = 0;
    for(int s=B-1;s<=B+1;s+=2){
        for(int t=L;t<=R;++t){
            // 求出每个点的出度（入度不需要）
            for(int i=L;i<=R;i++)out[i]=cnt[i]-(i==t);

            // 因为边不区分，要除以所有边的条数阶乘，提前计算
            mint r = 1;

            bool ok = 1;
            cnt_L[L] = 0;
            for(int i=L;i<=R;i++){
                cnt_R[i] = out[i] - cnt_L[i];
                // (i,i+1) 之间往右的边比往左的多一条
                if(s<=i && i<t)cnt_L[i+1]=cnt_R[i]-1;
                // (i,i+1) 之间往左的边等于往右的边
                else if(t<=i && i<s)cnt_L[i+1]=cnt_R[i]+1;
                // (i,i+1) 之间左右边数相等
                else cnt_L[i+1]=cnt_R[i];

                // 非法的 s 和 t
                if(cnt_L[i]<0 || cnt_R[i]<0){
                    ok = 0;
                    break;
                }
                r *= invf[cnt_L[i]]*invf[cnt_R[i]]; // 计算 r
            }

            if(!ok)continue;


            // 首先，这张图的点数是 pt，也就是 p 数组中独特元素个数
            // 容易证明，该图一定存在 s->t 的欧拉通路

            // 这个图里面还需要添加一条 t->s 的边才能用 BEST 定理
            // 但是因为欧拉图的基尔霍夫矩阵的主余子式全部相等
            // 我们求以 t 为根的根向树个数，就可以忽略这条边了
            
            // 所以根向树个数就是 Π(cnt_R[i]) Π(cnt_L[j])
            // 求出 t 往左/往右 (pt-1) 条边的乘积            
            vector<mint> prod1(pt),prod2(pt);
            prod1[0]=prod2[0]=1;
            for(int i=1;i<pt;i++){
                if(t-i<0 || !cnt_R[t-i])break;
                prod1[i]=prod1[i-1]*cnt_R[t-i];
            }
            for(int i=1;i<pt;i++){
                if(t+i>R || !cnt_L[t+i])break;
                prod2[i]=prod2[i-1]*cnt_L[t+i];
            }

            mint T = 0; // 根向树个数
            for(int i=0;i<pt;i++)T+=prod1[i]*prod2[pt-1-i];
            ans += T * res * r;
        }
    }
    cout << ans << '\n';
}


int main(){
    ios::sync_with_stdio(0);cin.tie(0);

    f[0]=invf[0]=1;
    for(int i=1;i<N;i++)f[i]=f[i-1]*i,invf[i]=invf[i-1]/i;
    int T; cin>>T;
    while(T--)solve();
    return 0;
}

参考资料

• 矩阵树定理及其无向图形式证明 - 洛谷专栏
• Cauchy-Binet Formula - ProofWiki
• 矩阵树定理 - OI Wiki