高斯消元

高斯消元是用于求解 nn 元线性方程组的算法。这里不讲别的,只讲算法流程。

首先,可以用一个 nn(n+1)(n + 1) 列的矩阵来表示一个 nn 元线性方程组:

{a1,1x1+a1,2x2++a1,nxn=a1,n+1a2,1x1+a2,2x2++a2,nxn=a2,n+1an,1x1+an,2x2++an,nxn=an,n+1\begin{cases} a_{1, 1}x_1 + a_{1, 2}x_2 + \cdots + a_{1, n} x_n = a_{1, n + 1}\\ a_{2, 1}x_1 + a_{2, 2}x_2 + \cdots + a_{2, n} x_n = a_{2, n + 1}\\ \cdots\\ a_{n, 1}x_1 + a_{n, 2}x_2 + \cdots + a_{n, n} x_n = a_{n, n + 1} \end{cases}

可以用下面这个矩阵来表示:

{a1,1a1,2a1,na1,n+1a2,1a2,2a2,na2,n+1an,1an,2an,nan,n+1}\begin{Bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} & | & a_{1, n + 1} \\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} & | & a_{2, n + 1} \\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} & | & a_{n, n + 1} \end{Bmatrix}

高斯消元会将它通过三个矩阵的初等变换变成一个阶梯型矩阵,下面是三种初等变换:

  1. 交换矩阵的两行。
  2. 将矩阵的一行元素全部乘上 λ\lambda ( λ0\lambda \neq 0 )。
  3. 将矩阵的一行减去 / 加上另一行。

阶梯型矩阵是矩阵的一种类型。它的基本特征是:所给矩阵为行阶梯型矩阵,则矩阵中每一行的第一个不为零的元素的左边及其所在列以下全为零。

高斯消元最后形成的阶梯型矩阵形如:

{a1,1a1,2a1,3a1,na1,n+10a2,2a2,3a2,na2,n+100a3,3a3,na3,n+1000an,nan,n+1}\begin{Bmatrix} a_{1, 1} & a_{1, 2} & a_{1, 3} & \cdots & a_{1, n} & | & a_{1, n + 1} \\ 0 & a_{2, 2} & a_{2, 3} & \cdots & a_{2, n} & | & a_{2, n + 1} \\ 0 & 0 & a_{3, 3} & \cdots & a_{3, n} & | & a_{3, n + 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots & | & \vdots \\ 0 & 0 & 0 & \cdots & a_{n, n} & | & a_{n, n + 1} \end{Bmatrix}

所以得到了消元之后的矩阵,只需要从最后一行不断往上代,就能分别求出 x1,x2,,xnx_1, x_{2}, \cdots, x_{n}

同时还有约旦·高斯消元,这种算法最后消出来的矩阵为:

{a1,1000a1,n+10a2,200a2,n+100a3,30a3,n+1000an,nan,n+1}\begin{Bmatrix} a_{1, 1} & 0 & 0 & \cdots & 0 & | & a_{1, n + 1} \\ 0 & a_{2, 2} & 0 & \cdots & 0 & | & a_{2, n + 1} \\ 0 & 0 & a_{3, 3} & \cdots & 0 & | & a_{3, n + 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots & | & \vdots \\ 0 & 0 & 0 & \cdots & a_{n, n} & | & a_{n, n + 1} \\ \end{Bmatrix}

这种消元不用往回代,只需要用最后一列的值除以对应行的系数就能得到 xix_i

矩阵求逆

定义:

如果矩阵 AA 存在 A1A^{-1},满足 A1×A=A×A1=EA^{-1} \times A = A \times A^{-1} = E,其中 EE 是单位矩阵。那么 A1A^{-1} 就是 AA 的逆矩阵,矩阵 AA 是可逆的。

求法:

可以通过下面的公式求得:

A1×[AE]=[EA1]A^{-1} \times [AE] = [EA^{-1}]

其中 [AB][AB] 表示 AABB 两个矩阵拼起来得到的矩阵。

可以用高斯·约旦消元求出结果,这里举个例子:

[210121012]\begin{bmatrix}2 & -1 & 0 \\\\-1 & 2 & -1 \\\\0 & -1 & 2\end{bmatrix}

的逆矩阵,首先写出 [AE][AE]

[210100121010012001]\begin{bmatrix}2 & -1 & 0 & 1 & 0 & 0 \\\\-1 & 2 & -1 & 0 & 1 & 0 \\\\0 & -1 & 2 & 0 & 0 & 1 \end{bmatrix}

对左边进行消元可得:

[21010003211210004313231]\begin{bmatrix}2 & -1 & 0 & 1 & 0 & 0 \\ 0 & \frac{3}{2} & -1 & \frac{1}{2} & 1 & 0 \\0 & 0 & \frac{4}{3} & \frac{1}{3} & \frac{2}{3} & 1 \end{bmatrix}

此时已消成上三角矩阵,高斯消元开始回代,用约旦会消成对角矩阵:

[200321120320343234004313231]\begin{bmatrix}2 & 0 & 0 & \frac{3}{2} & 1 & \frac{1}{2} \\0 & \frac{3}{2} & 0 & \frac{3}{4} & \frac{3}{2} & \frac{3}{4} \\0 & 0 & \frac{4}{3} & \frac{1}{3} & \frac{2}{3} & 1\end{bmatrix}

最后每行除以系数:

[10034121401012112001141234]\begin{bmatrix}1 & 0 & 0 & \frac{3}{4} & \frac{1}{2} & \frac{1}{4} \\0 & 1 & 0 & \frac{1}{2} & 1 & \frac{1}{2} \\ 0 & 0 & 1 & \frac{1}{4} & \frac{1}{2} & \frac{3}{4} \end{bmatrix}

此时右半边即为所求。

这里要判 AA 是否是一个可逆矩阵,其实就是在高斯·约旦算法消元时判断是否无解。


例题

[JSOI2008] 球形空间产生器

显然有 dist1=dist2==distn+1dist_1 = dist_2 = \cdots = dist_{n + 1},所以有 :

{(a1,1x1)2+(a1,2x2)2++(a1,nxn)2=(a2,1x1)2+(a2,2x2)2++(a2,nxn)2=(an+1,1x1)2+(an+1,2x2)2++(an+1,nxn)2{2(a2,1a1,1)2(a2,2a1,2)2(a2,na1,n)i=1n(a2,i2a1,i2)2(a3,1a1,1)2(a3,2a1,2)2(a3,na1,n)i=1n(a3,i2a1,i2)2(an+1,1a1,1)2(an+1,2a1,2)2(an+1,na1,n)i=1n(an+1,i2a1,i2)}\begin{aligned} &\begin{cases} (a_{1, 1} - x_1)^2 + (a_{1, 2} - x_2)^2 + \cdots + (a_{1, n} - x_n)^2 = \\ (a_{2, 1} - x_1)^2 + (a_{2, 2} - x_2)^2 + \cdots + (a_{2, n} - x_n)^2 = \\ \cdots\\ (a_{n+1, 1} - x_1)^2 + (a_{n+1, 2} - x_2)^2 + \cdots + (a_{n+1, n} - x_n)^2 \end{cases}\\ \Rightarrow& \begin{Bmatrix} 2(a_{2, 1} - a_{1, 1}) & 2(a_{2, 2} - a_{1, 2}) & \cdots & 2(a_{2, n} - a_{1, n}) & | & \sum\limits_{i=1}^n(a_{2, i}^2 - a_{1, i}^2)\\ 2(a_{3, 1} - a_{1, 1}) & 2(a_{3, 2} - a_{1, 2}) & \cdots & 2(a_{3, n} - a_{1, n}) & | & \sum\limits_{i=1}^n(a_{3, i}^2 - a_{1, i}^2)\\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ 2(a_{n+1, 1} - a_{1, 1}) & 2(a_{n+1, 2} - a_{1, 2}) & \cdots & 2(a_{n+1, n} - a_{1, n}) & | & \sum\limits_{i=1}^n(a_{n+1, i}^2 - a_{1, i}^2)\\ \end{Bmatrix} \end{aligned}

接着就可以用高斯消元解决此题。

[USACO09NOV] Lights G

根据异或的性质,显然一个点操作两次及以上是没有意义的。

所以可以写成一个异或方程组的形式:

{(V(1,1)×x1)(V(1,2)×x2)(V(1,n)×xn)=1(V(2,1)×x1)(V(2,2)×x2)(V(2,n)×xn)=1(V(n,1)×x1)(V(n,2)×x2)(V(n,n)×xn)=1\begin{cases} (V(1, 1) \times x_{1}) \oplus (V(1, 2) \times x_{2}) \oplus \cdots \oplus (V(1, n) \times x_{n}) = 1 \\ (V(2, 1) \times x_{1}) \oplus (V(2, 2) \times x_{2}) \oplus \cdots \oplus (V(2, n) \times x_{n}) = 1 \\ \cdots \\ (V(n, 1) \times x_{1}) \oplus (V(n, 2) \times x_{2}) \oplus \cdots \oplus (V(n, n) \times x_{n}) = 1 \\ \end{cases}

这上面 xix_i 表示点 ii 是否被操作过,Vi,jV_{i, j} 表示 ii 被操作是否会影响到 jj 的状态,也就是 i,ji, j 之间是否有连边(在这种定义下,需要钦定 Vi,i=1V_{i, i} = 1)。

这种方程组也可以用高斯消元消掉,如果有自由元就要深搜一下自由元状态,然后这题就没了。

[HNOI2011] XOR和路径

首先考虑 dp,设 fuf_u 表示以 uu 为起点的答案,有:

fu=1degu(u,v,w)E(fvw)f_{u} = \frac{1}{deg_u} \cdot \sum\limits_{(u, v, w) \in E} \left( f_{v} \oplus w\right)

这里 degudeg_u 表示 uu 点的度数(注意自环只算一次)。

但是这道题会有环,所以正常 dp 会有后效性,于是考虑用高斯消元来解决后效性问题。上面的式子还不能直接用高斯消元,因为有异或,考虑把异或给转化掉。原式等于:

degu×fu=(u,v,w)E(fvw)deg_u \times f_u = \sum\limits_{(u, v, w) \in E} (f_{v} \oplus w)

接着由于期望的线性性,所以可以拆位计算异或,拆位之后,式子变成:

degu×fu=(u,v,w)Ew=0fv+(u,v,w)Ew=1(1fv)deg_u \times f_u = \sum\limits_{(u, v, w') \in E \land w' = 0} f_v + \sum\limits_{(u, v, w') \in E \land w' = 1} (1 - f_v)

其中 ww'ww 拆位之后的值(值只能是 0011),再移项:

(u,v,w)Ew=0fv(u,v,w)Ew=1fvdegu×fu=(u,v,w)Ew=11\sum\limits_{(u, v, w') \in E \land w' = 0} f_v - \sum\limits_{(u, v, w') \in E \land w' = 1} f_v - deg_u \times f_u = - \sum\limits_{(u, v, w') \in E \land w' = 1} 1

对于 1(n1)1 \sim (n - 1) 都有一个这样的约束,对于 nn,有 fn=0f_n = 0。所以用高斯消元就可以求出 f1f_1,答案就是 k=0logV2kf1\sum\limits_{k = 0}^{\log V} 2^k \cdot f_1'。时间复杂度 O(n3logV)\mathcal{O}(n^3\log V)

[HNOI2013] 游走

首先可以发现因为期望的可加性,所以答案是每条边的编号乘上每条边的期望经过次数。现在编号是由我们自己定的,每条边的期望经过次数是不变的,由于要最小化答案所以我们一定要让大的乘小的、小的乘大的。现在就考虑求边的期望经过次数,有点难以求得,不如先求点的期望经过次数,再通过点的经过次数得到边的:假设一条边 (u,v)(u, v) ( u,vnu, v \ne n ),设 fuf_u 是点 uu 的期望经过次数,那么 (u,v)(u, v) 这条边的期望经过次数就等于 fu×1degu+fv×1degvf_u \times \frac{1}{deg_u} +f_v \times \frac{1}{deg_v}

接着就考虑求 ffff 有朴素转移式子:

fu=(u,v)Evnfv×1degv+[u=1]f_{u} = \sum\limits_{(u, v) \in E \land v \ne n} f_{v} \times \frac{1}{deg_v} + [u = 1]

并且有 fn=1f_n = 1,这是一个显然有后效性的东西(因为原图不是 DAG),所以考虑用高斯消元求:

fu(u,v)Evnfv×1degv=[u=1]f_u - \sum\limits_{(u, v) \in E \land v \ne n} f_{v} \times \frac{1}{deg_v} = [u = 1]

最后算出边的期望就行了。

First Knight

这道题题意就是将 (n,m)(n, m) 大小的矩阵,每个点有概率往上下左右走,且不可能走出矩阵。要求求出 (1,1)(1, 1)(n,m)(n, m) 的期望步数。

很显然,有一个 dp,设 fi,jf_{i, j} 表示从 (i,j)(i, j)(n,m)(n, m) 的期望步数,那么有:

fi,j={0(i=nj=m)U(i,j)fi1,j+D(i,j)fi+1,j+L(i,j)fi,j1+R(i,j)fi,j+1otherwisef_{i, j} = \begin{cases} 0 & (i = n \land j = m) \\ U(i, j) \cdot f_{i - 1, j} + D(i, j) \cdot f_{i + 1, j} + L(i, j) \cdot f_{i, j - 1} + R(i, j) \cdot f_{i, j + 1} & \text{otherwise} \end{cases}

然后就可以移项,变成一个高斯消元,看起来是 O(n6)\mathcal{O}(n^6) 的。但是实际上,在消其它项的系数的时候,可以发现真正被消的只有 O(1)\mathcal{O}(1) 个系数,所以实际上是 O(n4n5)\mathcal{O}(n^4 \sim n^5) 的(太蒻了不会算)。

[USACO10HOL] Driving Out the Piggies G

这题和游走非常相似啊,设 fif_i 表示炸弹在 ii 爆炸的概率,有:

fu=pq×([u=1]+(u,v)Efv×1degv×(1pq))f_u = \frac{p}{q} \times \left( [u = 1] + \sum\limits_{(u, v) \in E} f_{v} \times \frac{1}{deg_v} \times (1 - \frac{p}{q} ) \right)

接着移项然后高斯消元就行了。


矩阵树定理

行列式的定义

行列式是一个函数定义,取值是一个标量。

对一个 n×nn \times n 的矩阵 AA,其 nn 阶行列式写作 det(A)\det(A) 或者 A| A |,定义为:

det(A)=A=P(1)τ(P)i=1nAi,Pi\det(A) = |A| = \sum\limits_{P}(-1)^{τ(P)}\prod\limits_{i = 1}^n A_{i, P_i}

PP 是一个排列,τ(P)\tau(P) 表示一个排列的逆序对数。

数学上的定义是:行列式可以理解为所有列向量所夹的几何体的有向体积,这样可以结合从几何直观出发理解为线性变换的伸缩因子。

但是 OI 中,不用管那么多数学意义。只需求出行列式的值即可。

行列式的性质

行列式的性质决定了高斯消元求行列式的可行性,性质如下:

  • **交换行列式的两行,行列式变号。**这个很好证:交换了两行,相当于对于每个排列,都交换了两个元素。排列中交换两个元素会使逆序对数奇偶性改变,首先交换排列中两个相邻元素会使奇偶性改变,而交换两个元素相当于交换 2x12x - 1 次相邻元素。(1)(1)
  • 交换排列的一行和一列,行列式不变。 (2)(2)
  • 如果 CC 中有一行是 AA 对应的行与 BB 对应的行之和,其余元素都是相同的。有 det(C)=det(A)+det(B)\det(C) = \det(A) + \det(B)(3)(3)

[a1,1a1,2a1,nai,1+bi,1ai,2+bi,2ai,n+bi,nan,1an,2an,n]=[a1,1a1,2a1,nai,1ai,2ai,nan,1an,2an,n]+[a1,1a1,2a1,nbi,1bi,2bi,nan,1an,2an,n]\begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\a_{i, 1} + b_{i, 1} & a_{i, 2} + b_{i, 2} & \cdots & a_{i, n} + b_{i, n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\ \end{bmatrix} = \begin{bmatrix} 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{bmatrix} + \begin{bmatrix}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{bmatrix}

  • 带入式子时,用乘法分配律不难证明。
  • (3)(3) 可以推出,当矩阵的某一行集体乘上 λ\lambda 时,行列式也乘上 λ\lambda(4)(4)
  • 如果矩阵中,有一行与另一行成比例,那么行列式的值为 00(5)(5)

[a1,1a1,2a1,nλa1,1λa1,2λa1,nan,1an,2an,n]=0\begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ \lambda \cdot a_{1, 1} & \lambda \cdot a_{1, 2} & \cdots & \lambda \cdot a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\ \end{bmatrix} = 0

  • 证明:首先根据 (4)(4)λ\lambda 提出来。再根据 (1)(1) 将相同的两行交换,此时矩阵应该不变,并且由 (1)(1) 行列式的值应该变号。由于 det(A)=det(A)\det(A) = -\det(A),解得 det(A)=0\det(A) = 0
  • 将矩阵的 一行集体乘上 λ\lambda 的值 加到矩阵的另一行,行列式不变(6)(6)

[a1,1a1,2a1,nai,1+λa1,1ai,2+λa1,2ai,n+λa1,nan,1an,2an,n]=[a1,1a1,2a1,nai,1ai,2ai,nan,1an,2an,n]\begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{i, 1} + \lambda \cdot a_{1, 1} & a_{i, 2} +\lambda \cdot a_{1, 2} & \cdots & a_{i, n} + \lambda \cdot a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\ \end{bmatrix} = \begin{bmatrix} 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{bmatrix}

  • 证明:由 (3)(3),可以将式子变成:

[a1,1a1,2a1,nai,1+λa1,1ai,2+λa1,2ai,n+λa1,nan,1an,2an,n]=[a1,1a1,2a1,nai,1ai,2ai,nan,1an,2an,n]+[a1,1a1,2a1,nλa1,1λa1,2λa1,nan,1an,2an,n]\begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{i, 1} + \lambda \cdot a_{1, 1} & a_{i, 2} +\lambda \cdot a_{1, 2} & \cdots & a_{i, n} + \lambda \cdot a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\ \end{bmatrix} = \begin{bmatrix} 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{bmatrix} + \begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ \lambda \cdot a_{1, 1} & \lambda \cdot a_{1, 2} & \cdots & \lambda \cdot a_{1, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\ \end{bmatrix}

  • 不难发现后面那坨东西的行列式是 00

行列式求值

暴力:O(n2×n!)\mathcal{O}(n^2 \times n!)

高斯消元:

由于行列式的性质 (1),(4),(6)(1), (4), (6),可以用高斯消元将其消成上三角矩阵:

[a1,1a1,2a1,3a1,n0a2,2a2,3a2,n00a3,3a3,n000an,n]\begin{bmatrix} a_{1, 1} & a_{1, 2} & a_{1, 3} & \cdots & a_{1, n} \\ 0 & a_{2, 2} & a_{2, 3} & \cdots & a_{2, n} \\ 0 & 0 & a_{3, 3} & \cdots & a_{3, n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & a_{n, n} \end{bmatrix}

这种特殊行列式的值是 i=1nai,i\prod\limits_{i = 1}^n a_{i, i}

因为考虑 i=1nai,Pi\prod\limits_{i = 1}^n a_{i, P_i} 中一旦某个 ai,Pia_{i, P_i} 选到 00 整个就为 )) 了。所以:Pn=1(n1)P_{n} = 1 \sim (n - 1) 时乘积都为 00,所以 Pn=nP_{n} = n,而 Pn1P_{n - 1} 由于 PnP_n 取了 nn 了只能取 n1n - 1,以此类推。对答案有贡献的排列 PP 只有 Pi=iP_i = i 这一个。

所以用高斯消元将其化为上三角式子后,再判断一下符号就能求出行列式的值了。

但是高斯消元会用到除法,有的时候会发现模数不是个质数,这样有很多数会找不到逆元。此时就要用辗转相除法。也就是像 gcd\gcd 那样,不断地交换并相减,根据 (1),(6)(1), (6),这样只会使符号改变,解决了模数(或者精度)问题。


代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
template <class T>
T power(T a, long long b) {
T res = 1;
for (; b; b >>= 1, a *= a) {
if (b & 1)
res *= a;
} return res;
}
template <long long mod>
class ModLL {
public:
long long n;
static long long Mod;
constexpr ModLL() : n{} {}
constexpr ModLL(long long x) : n(norm(x % getmod())) {}
constexpr long long norm(long long x) {
if (x >= getmod()) x %= getmod();
if (x <= -getmod()) x %= getmod();
if (x < 0) x += getmod();
return x;
}
constexpr long long getmod() {return (mod > 0 ? mod : Mod);}
explicit constexpr operator long long() const {return n;}
constexpr ModLL operator -() const {ModLL a; a.n = norm(getmod() - n); return a;}
constexpr ModLL inv() {assert(n != 0); return power((*this), getmod() - 2);}
constexpr ModLL &operator += (ModLL w) & {n = norm( n + w.n); return (*this);}
constexpr ModLL &operator -= (ModLL w) & {n = norm( n - w.n); return (*this);}
constexpr ModLL &operator *= (ModLL w) & {n = norm( 1LL * n * w.n % getmod()); return (*this);}
constexpr ModLL &operator /= (ModLL w) & {return (*this) *= w.inv();}
friend constexpr ModLL operator + (ModLL a, ModLL b) {ModLL res = a; res += b; return res;}
friend constexpr ModLL operator - (ModLL a, ModLL b) {ModLL res = a; res -= b; return res;}
friend constexpr ModLL operator * (ModLL a, ModLL b) {ModLL res = a; res *= b; return res;}
friend constexpr ModLL operator / (ModLL a, ModLL b) {ModLL res = a; res /= b; return res;}
friend constexpr bool operator == (ModLL a, ModLL b) {return (a.n == b.n);}
friend constexpr bool operator != (ModLL a, ModLL b) {return (a.n != b.n);}
friend constexpr istream &operator >> (istream &is, ModLL &a) {
long long x = 0; is >> x;
a = ModLL(x); return is;
}
friend constexpr ostream &operator << (ostream &os, const ModLL &a) {return (os << (a.n));}
} ;
template <int mod>
class ModInt {
public:
int n;
static int Mod;
constexpr ModInt() : n{} {}
constexpr ModInt(int x) : n(norm(x % getmod())) {}
constexpr int norm(int x) {
if (x >= getmod()) x %= getmod();
if (x <= -getmod()) x %= getmod();
if (x < 0) x += getmod();
return x;
}
constexpr static int getmod() {return (mod > 0 ? mod : Mod);}
explicit constexpr operator int() const {return n;}
constexpr ModInt operator -() const {ModInt a; a.n = norm(getmod() - n); return a;}
constexpr ModInt inv() const {assert(n != 0); return power((*this), getmod() - 2);}
constexpr ModInt &operator += (ModInt w) & {n = norm( n + w.n); return (*this);}
constexpr ModInt &operator -= (ModInt w) & {n = norm( n - w.n); return (*this);}
constexpr ModInt &operator *= (ModInt w) & {n = norm( 1LL * n * w.n % getmod()); return (*this);}
constexpr ModInt &operator /= (ModInt w) & {return (*this) *= w.inv();}
friend constexpr ModInt operator + (ModInt a, ModInt b) {ModInt res = a; res += b; return res;}
friend constexpr ModInt operator - (ModInt a, ModInt b) {ModInt res = a; res -= b; return res;}
friend constexpr ModInt operator * (ModInt a, ModInt b) {ModInt res = a; res *= b; return res;}
friend constexpr ModInt operator / (ModInt a, ModInt b) {ModInt res = a; res /= b; return res;}
friend constexpr bool operator == (ModInt a, ModInt b) {return (a.n == b.n);}
friend constexpr bool operator != (ModInt a, ModInt b) {return (a.n != b.n);}
friend constexpr istream &operator >> (istream &is, ModInt &a) {
int x = 0; is >> x;
a = ModInt(x); return is;
}
friend constexpr ostream &operator << (ostream &os, const ModInt &a) {return (os << (a.n));}
} ;
template <>
long long ModLL <0> :: Mod = (long long)(1E18) + 9;
template <>
int ModInt <0> :: Mod = 998244353;
using P = ModInt <0>;
const int N = 605;
int n;
P *a[N], _a[N][N];
P det(P **a) {
bool fh = 0; P ans = 1;
for (int i = 1; i <= n; ++i) {
for (int j = i + 1; j <= n; ++j) {
while (a[i][i].n > 0) {
int div = a[j][i].n / a[i][i].n;
for (int k = 1; k <= n; ++k)
a[j][k] -= P(div) * a[i][k];
swap(a[i], a[j]); fh ^= 1;
}
swap(a[i], a[j]); fh ^= 1;
}
}
for (int i = 1; i <= n; ++i)
ans *= a[i][i];
if (fh) ans *= P(-1);
return ans;
}
signed main(void) {
ios :: sync_with_stdio(false);
cin.tie(0); cout.tie(0);
cin >> n; int p; cin >> p; P :: Mod = p;
for (int i = 1; i <= n; ++i) a[i] = _a[i];
for (int i = 1; i <= n; ++i) for (int j = 1; j <= n; ++j)
cin >> a[i][j];
cout << det(a) << '\n';
return 0;
}

下面引用了 Achtoria 的这篇文章

矩阵树定理定义

  1. 对应无向图,定义 D(G)D(G) 为图 GG 的度数矩阵,其中:

D(G)(i,j)={degi(i=j)0(ij)D(G)(i, j) = \begin{cases}deg_i & (i = j) \\0 & (i \ne j)\end{cases}

  1. 对于有向图,定义 Din(G)D^{in}(G)GG 的入度矩阵,Dout(G)D^{out}(G) 为图 GG 的出度矩阵,其中:

Din(G)(i,j)={ini(i=j)0(ij),Dout(G)(i,j)={outi(i=j)0(ij)D^{in}(G)(i, j) = \begin{cases}in_i & (i = j) \\ 0 & (i \ne j)\end{cases}, D^{out}(G)(i, j) = \begin{cases}out_i & (i = j) \\ 0 * (i \ne j)\end{cases}

  1. 定义 A(G)A(G) 为图 GG 的邻接矩阵,其中:

A(G)(i,j)=the numebr of the edges from i to jA(G)(i, j) = \text{the numebr of the edges from } i \text{ to } j

  1. 定义图 GG 的基尔霍夫矩阵 K(G)=D(G)A(G)K(G) = D(G) - A(G)
  2. 定义无向图 GG 的生成树数量为 t(G)t(G),有向图内向树生成树(根为 uu)数量为 tout(G,u)t^{out}(G, u),有向图外向树数量为 tin(G,u)t^{in}(G, u)

矩阵树定理与扩展

  • 经典:
    给出一个无向无权图 GG,令 K(G)K'(G)K(G)K(G) 去掉第 kk 行和第 kk 列(kk 任意)的 n1n - 1 阶主子式。
    det(K)\det(K')GG 的生成树数量。
    不会证明,有兴趣可以看看上面引用的博客。
  • 加权扩展:
    根据乘法原理,对于某种生成树的形态,其对答案的贡献是每条边重的次数的乘积
    所以带权边 (u,v,w)(u, v, w) 相当于是建了 ww(u,v)(u, v) 这条边。注意这里度数矩阵也要边。
  • 有向扩展:
    有向图上也是可以对外向树或内向树计数的,这里有定理:
    定义 Kin(G)=Din(G)A(G)K^{in}(G) = D^{in}(G) - A(G)Kout(G)K^{out}(G) 同理;
    则:

i[1,n]tout(G,i)=det(Kout(G)(1,2,,i1,i+1,,n1,2,,i1,i+1,,n))i[1,n]tin(G,i)=det(Kin(G)(1,2,,i1,i+1,,n1,2,,i1,i+1,,n))\begin{aligned} \forall i \in [1, n] t^{out}(G, i) &= \det\left(K^{out}(G) \binom{1, 2, \cdots, i - 1, i + 1, \cdots, n}{1, 2, \cdots, i - 1, i + 1, \cdots, n}\right) \\ \forall i \in [1, n] t^{in}(G, i) &= \det\left(K^{in}(G) \binom{1, 2, \cdots, i - 1, i + 1, \cdots, n}{1, 2, \cdots, i - 1, i + 1, \cdots, n}\right) \\ \end{aligned}


例题