矩阵的初等变换

以下的矩阵大小都是 $n \times n$ 的，即为方阵。

初等行变换

作用：交换第 $i,j$ 行的所有元素，即构造矩阵 $T$ 使得 $TA$ 交换了 $A$ 矩阵的第 $i,j$ 行。（$i<j$）

构造方式：

$$
T=\begin{bmatrix}
I_{i-1} & & & & \\
& 0 & &1 \\
& & I_{j-i-1} &\\
& 1 & &0 \\
& & & & I_{n-j} \\
\end{bmatrix}
$$

空的地方都是 $0$。

倍法变换

作用：让第 $i$ 行的元素都乘上系数 $k$。

构造方式：

$$
T=\begin{bmatrix}
I_{i-1} & & \\
& k & \\
& & I_{n-i} \\
\end{bmatrix}
$$

空的地方都是 $0$。

消法变换

作用：让第 $i$ 行的元素都乘上系数 $k$ 并且加到第 $j$ 行上。（$i<j$）

构造方式：

$$
T=\begin{bmatrix}
I_{i-1} & & & & \\
& 1 & &0 \\
& & I_{j-i-1} &\\
& k & &1 \\
& & & & I_{n-j} \\
\end{bmatrix}
$$

空的地方都是 $0$。

这样我们就证明了高斯消元中的所有初等变换都可以写成矩阵 $A$ 左乘某个矩阵 $T$ 的方式，即 $A \gets TA$。

高斯（约当）消元

基本流程

给定一组 $n$ 元一次方程，试求出这个方程的解。

一般情况下可以视作给定矩阵 $A=m \times n,B = 1 \times m$，求出一个 $X=n \times 1$，使得 $AX=B$。

很显然，根据初中数学学到的知识，对于第 $i$ 个方程，如果它第 $i$ 项的系数不为 $0$，则可以把其他方程第 $i$ 项的系数都消掉而且不影响其他项的非 $0$ 性。

如果第 $i$ 项的系数为 $0$，那就交换一行上来，如果任意方程第 $i$ 项的系数都为 $0$，就是下面要说的一种情况。

然后在理想情况下得到的方程最后是这个样子的：

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

那么这个时候直接用右边的值除以这个矩阵的系数就可以了。

但是可能会有下列特例：

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

即 $n-2,n-1,n$ 没有值的限制，那么这个时候我们把这些元素称作自由元，值得注意的是：自由元互相不影响，即如果每个数字有 $P$ 种取值的话，整个方程解得个数就是 $P^\text{自由元的数量}$。

如果我们要给出方程的一组解，那就需要先确定自由元的值，然后因为其它值可能不确定，就需要代入自由元的取值，从矩阵的第 $n$ 行倒推回去，这样子做肯定是对的，我们可以考虑消元的顺序。

还有一种情况，就是无解，无解表现和自由元非常相似，只是后面系数都为 $0$ 的行，常数却不为 $0$，这样子方程就产生了冲突，即无解。

时间复杂度：$O(N^3)$，空间复杂度：$O(N^2)$。

减少数量

如果我们可以用一组变量 $a_1 \times a_n$，表示所有其它的数，且满足其它的数的方程的话，我们就可以将方程和未知数的数量缩减到 $O(n)$，而并非原来的 $O(n^2)$ 及更多。

应用场合

通常只需要把方程列出来，然后减少数量之后求解。

注意到消元的过程中只有加减乘除运算，故对于某个数取模意义下的方程也是同样的操作。

特别是对于 $2$ 取模（bool 变量的方程），因为位运算中，单个位的加减法都是异或，乘除法都是与运算，我们可以使用 bitset 优化至 $O(\dfrac{N^3}{64})$ 的时间复杂度。

或者是解一些期望概率的题。

这里顺带提一下：

图上随机游走类问题，如果要求点的期望，那么方程可以很轻松的列出来，即：

$$
dp_u = \begin{cases}\sum_{v \to u} \dfrac{dp_v}{d_v}+1 & u=1 \\ \sum_{v \to u} \dfrac{dp_v}{d_v} & u \ne 1\end{cases}
$$

特别的，如果限制了终点，特判即可，因为到了终点就停止了，所以到终点的概率就是它的期望。

这是最重要的结论之一。

矩阵的逆

如果一个矩阵 $A = n \times n$ 存在另外一个 $B = n \times n$ 的矩阵使得 $A \times B = I$，则 $B$ 称作 $A$ 的逆。

那么有一个定理，将一个矩阵 $[AI]$ 通过初等行变换成一个阶梯形最简矩阵，那么如果这个矩阵的左边是单位矩阵 $I$，代表 $A$ 这个矩阵有逆，逆就是右边的矩阵，即化简后变成 $[IA^{-1}]$ 的形式。

证明

我们可以知道任意一种初等行变换都可以表示成矩阵 $A \gets TA$ 的形式，那么高斯消元实际上就是乘了很多个不同的 $T$，并且这些矩阵是有结合律的。

最后高斯（约当）消元是需要每行除以一个系数，变成：

$$
\begin{bmatrix}
1 & 0 & 0 & \cdots &0& | & x_1 \\
0 & 1 & 0 & \cdots &0& | & x_2 \\
0 & 0 & 1 & \cdots &0& | & x_3 \\
0 & 0 & 0 & \ddots &0& | & x_k \\
0 & 0 & 0 & \ddots &1& | & x_n \\
\end{bmatrix}
$$

右边的通过竖线分开的是增广矩阵，可以省略，然后左边就是单位矩阵了，这个时候我们相当于对原来的矩阵乘了一些 $T$，那么怎么找出这个 $T$ 呢？

于是我们在 $A$ 的右边添加一个矩阵 $I$，变换的时候顺便在 $I$ 上记录了乘上的矩阵 $T$，最后输出就可以了。

证毕。

如果高斯（约当）消元中碰到了无解的情况，即有一个数字被消完了，并且其每一个系数都为 $0$，则 $A$ 没有逆矩阵。

那么代码就很好写了，最后记得除以每一行的系数，总时间复杂度为 $O(n^3)$。

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

#include<bits/stdc++.h>
#define mod 1000000007
#define ll long long
#define N 805
using namespace std;
struct matrix{ll a[N][N],n1,n2;}maps,inv;
ll n,m1,m2,i,j,k,x,K,temp,ans=LLONG_MAX,pos;
ll qmi(ll a,ll b,ll p){
	ll res = 1%p,t=a;
	while(b){
		if(b&1) res=res*t%p;
		t=t*t%p;
		b>>=1;
	}
	return res;
}
int main(){
	ios::sync_with_stdio(false);
	cin>>n;
	for(i=1;i<=n;i++) for(j=1;j<=n;j++) cin>>inv.a[i][j];
	for(i=1;i<=n;i++) inv.a[i][i+n]=1;
	for(i=1;i<=n;i++){
		pos = i;
		for(j=i+1;j<=n;j++) if(inv.a[j][i]>inv.a[i][i]) pos=j;
		for(j=1;j<=2*n;j++) swap(inv.a[i][j],inv.a[pos][j]);
		if(inv.a[i][i]==0){
			cout<<"No Solution\n";
			return 0;
		}
		for(j=1;j<=n;j++){
			if(i==j) continue;
			temp = inv.a[j][i]*qmi(inv.a[i][i],mod-2,mod)%mod;
			for(k=1;k<=2*n;k++) inv.a[j][k]-=(inv.a[i][k]*temp)%mod,inv.a[j][k]=(inv.a[j][k]%mod+mod)%mod;
		}
	}
	for(i=1;i<=n;i++){
		ll temp = qmi(inv.a[i][i],mod-2,mod);
		for(j=1;j<=2*n;j++) inv.a[i][j] = inv.a[i][j]*temp%mod;
	}
	for(i=1;i<=n;i++){
		for(j=1;j<=n;j++){
			cout<<inv.a[i][j+n]<<" ";
		}
		cout<<endl;
	}
	return 0;
}
/*
Input:
4 2 2
1 3
2 4
1 1 1 1
0 0 0 0

Output:
4
*/

线性相关

如果设 $a$ 是一个长度为 $m$ 的向量，设其为 $a_1,a_2,\dots,a_m$，那么设 $ka={ka_1,ka_2,\dots,ka_m}$，如果若干个向量满足存在实数集 $S$ 满足：

$$
S_1a_1+S_2a_2+S_3a_3+\dots+s_na_n=0
$$

则称这 $n$ 个向量线性相关，特别的我们称一个向量空间（即一个满足运算封闭性的向量集合）的线性基，为选出最多的向量空间中的向量，满足这些向量线性无关。

可以证明这些向量可以表出这个向量空间内的所有数，不能表出任意一个向量空间外的数。

表出代表只使用加法和乘法对于向量和实数进行操作得到另外的向量。

矩阵的秩

矩阵有两种秩：行秩和列秩。行秩就是把矩阵的每一个行看成向量，然后它所属的向量空间的线性基的大小；列秩就是把矩阵的每一个列看成向量，然后它所属的向量空间的线性基的大小。

可以证明出来行秩等于列秩。

证明：

初等行变换不影响矩阵的行秩，最后行秩等于 $n$ 减去自由变元的数量。

初等行变换也不影响矩阵的列秩。

交换两行不影响矩阵的列秩，因为相当于交换了两个未知数。

一行同时乘上一个数不影响矩阵的列秩，因为相当于某个未知数扩大了 $k$ 倍。

一行乘上一个数加到另一行上也不影响矩阵的列秩，因为相当于一个方程 $ax+b$，并且通过前两者结合也可以证明。

证毕。

所以任意矩阵的秩可以记作 $\operatorname{rank}(A)$ 或者 $r(A)$，它既可以表示列秩，也可以表示行秩。

如果一个 $n \times n$ 的方阵 $A$，如果 $r(A)=n$，则称它满秩。

如果一个 $n \times m$ 的矩阵 $A$，如果 $r(A)=n$，则称它行满秩；如果 $r(A)=m$，则称它列满秩。

显然 $r(A) \le \min(n,m)$。

重要结论

• 矩阵满秩、矩阵可逆、线性方程组有唯一解可以相互转化。

• 系数矩阵的秩等于增广矩阵的秩的时候，线性方程组有解。

• 满秩矩阵可以拆分成若干个初等变换矩阵的乘积。（考虑高斯消元的操作）

• 任意矩阵乘以满秩矩阵，秩不改变。（用上一条结论证明）

• $r(A+B) \le r(A)+r(B)$

• $\max(r(A),r(B)) \le r([A|B]) \le r(A)+r(B)$

• $r(AB) \le \min(r(A),r(B))$

行列式

定义

行列式是一个 $n \times n$ 的方阵，每个方阵 $A$ 对应了一个值 $\det(A)$，这个值就叫做这个行列式的值。

如果我们设 $\pi(p)$ 表示长度为 $n$ 的排列 $p$ 的逆序对数量，那么就有：

$$
\det(A) = \sum_{p}(-1)^{\pi(p)}\prod_{i=1}^n A_{i,p_i}
$$

这是行列式的一个全排列方法定义，由此可得矩阵 $A$ 的行列式与矩阵 $A^{T}$ （$A$ 的转置）的行列式相等。

或者是一种递归的解法，任取一行 $i$，则：

$$
\det(A)=\sum_{j=1}^n (-1)^{i+j}A_{i,j}\det(M_{i,j})
$$

其中 $M_{i,j}$ 表示从 $A$ 中删去第 $i$ 行和第 $j$ 列的数的矩阵，显然它是一个 $n-1 \times n-1$ 的方阵。

或者我们可以任取一行 $j$，则：

$$
\det(A)=\sum_{i=1}^n (-1)^{i+j}A_{i,j}\det(M_{i,j})
$$

前者叫做按第 $i$ 行展开，后者叫做按第 $j$ 列展开，所以不论按哪一行展开或者哪一列展开，所得结果都是一样的。

于是有结论：

定理：行列式 $\det A$ 的某一行（或某一列）的元素与另外一行（或另外一列）对应元素的代数余子式的乘积之和等于 $0$。

换句话说，当 $i\neg j$ 时：

$$a_{1,i}M_{1,j}+a_{2,i}M_{2,j}+\cdots+a_{i,n}M_{j,n}=0$$

$$a_{1,i}M_{1,j}+a_{2,i}M_{2,j}+\cdots+a_{n,i}M_{n,j}=0$$

计算

暴力计算行列式的时间复杂度达到了指数级别，我们可以利用矩阵的初等变换来减少行列式的计算时间，从而让其降为 $O(n^3)$。

交换两行

如果我们交换两行，相当于对于每个排列 $p$，都交换了 $i,j$ 两个位置上的值，那么显然可得 $\pi(p)$ 会变成原来的相反数，绝对值不变。

所以交换两行用一个变量来记录取反了多少次即可。

一行乘上 $k$

显然，行列式的值也会乘上 $k$。

一行乘上 $k$ 加到另外一行

引理，如果：

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

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

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

那么 $\det(A)=\det(B)+\det(C)$，我们通过全排列的乘法分配律就可以得到。

所以我们把操作的行拆成原来的部分，和 $k$ 乘上另外一行的部分，原来的部分显然行列式的值不会变，那么我们需要考虑第 $j$ 行等于第 $i$ 行乘上 $k$ 的行列式的值。（$i \ne j$）

我们还是考虑全排列，因为交换两行的话 $\pi(p)$ 的奇偶性会改变，我们也可以把这一行的 $k$ 提出来，这样就有完全相同的两行，那么显然答案为 $0$。

所以这种变换答案不变。

注意：以上行的操作通过全排列的性质也可以对列通用。

过程

因此，我们可以使用高斯消元消成上三角矩阵，这个时候 $A_{i,j}=0(i>j)$，所以行列式的值显然等于主对角线的乘积。

代码如下，时间复杂度为 $O(n^2(n+\log p))$，并且是没有逆元的解法，可以使用辗转相减法避免除法。

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

#include<bits/stdc++.h>
#define ll long long
#define N 605
using namespace std;
ll n,a[N][N],i,j,k,mod,ans=1,has,res;
int main(){
	ios::sync_with_stdio(false);
	cin>>n>>mod;
	for(i=1;i<=n;i++) for(j=1;j<=n;j++) cin>>a[i][j];
	for(i=1;i<=n;i++){
		for(j=i+1;j<=n;j++){
			while(a[i][i]){
				res = a[j][i]/a[i][i];
				for(k=i;k<=n;k++) a[j][k]=(a[j][k]-a[i][k]*res%mod+mod)%mod;
				swap(a[i],a[j]),has^=1;
			}
			swap(a[i],a[j]),has^=1;
		}
	}
	if(i<=n){
		cout<<0<<endl;
		return 0;
	}
	for(i=1;i<=n;i++) ans=ans*a[i][i]%mod;
	if(has) ans=(mod-ans)%mod;
	cout<<ans<<endl;
	return 0;
} 
/*
Input:
2 998244353
1 4
1 5

Output:
1
*/

性质

凸包

设 $k$ 维空间中存在 $k$ 个点，则这 $k$ 个点和原点构成的凸包面积 $\times k!$ 就是：

$$
A=\begin{vmatrix}
x_{1,1} & x_{2,1} & \cdots & x_{k,1} \\
x_{1,2} & x_{2,2} & \cdots & x_{k,2} \\
\vdots & \vdots & \ddots & \vdots \\
x_{1,k} & x_{2,k} & \cdots & x_{k,k}
\end{vmatrix}
$$

行列式的值为 $\det(A)$，答案就是 $|\det(A)|$ 这是一个结论，二维的情况显然可以证明，三维及以上就需要想象力了。

如果是给定的 $k+1$ 个点构成的凸包，我们就需要任意指定一个点为原点建系即可。

柯西-比内 Cauchy–Binet 公式

设 $A \in R^{n \times m},B \in R^{m \times n},n \le m$，那么有：

$$
\det(AB) = \sum_{S} \det(A_S)\det(B_S)
$$

其中 $S$ 表示下标集合 ${1,2,\dots,m}$ 的 $n$ 元子集，$A_S$ 表示 $A$ 取 $S$ 的列构成的 $n$ 元子式，$B$ 恰好相反。

当 $A,B$ 为方阵的时候 $\det(AB)=\det(A)\det(B)$。

范德蒙德 Vandermonde 行列式

定义为：

$$
A = \begin{vmatrix}
1 & 1 & \dots & 1 \\
x_1 & x_2 & \dots & x_n \\
\vdots & \vdots & \cdots & \vdots \\
x_1^{n-1} & x_2^{n-1} & \dots & x_n^{n-1} \\
\end{vmatrix}
$$

有 $\det(A)=\prod_{1 \le i<j \le n}(x_j-x_i)$，证明可以用高斯消元直接代入求值。

LGV 引理

定义

Lindström–Gessel–Viennot lemma 适用范围：有向无环图。

首先设 $w(P)$ 表示 $P$ 这条路径上的边权的乘积，如果是计数类问题我们可以把所有边权设为 $1$。

我们设 $e(u,v)$ 表示 $u$ 到 $v$ 的每一条路径的 $w(P)$ 之和，$A$ 为长度为 $n$ 的起点集合，$B$ 为长度为 $n$ 的终点集合。

设一组 $A$ 到 $B$ 的不相交集合 $S$ 为 $S_i$ 与 $S_j$ 没有公共端点，并且每一条路径恰好连接了 $A$ 中的一个数和 $B$ 中的一个数，而每个顶点也恰好只有一条路径连接它，由此可见 $S$ 的长度一定为 $n$，也即第 $i$ 条路径 $S_i$ 的起点是 $A_i$，终点是 $B_{p_i}$，$p$ 是长度为 $n$ 的排列，设 $\pi(p)$ 表示 $p$ 中逆序对的个数。

则有：

$$ M = \begin{bmatrix}e(A_1,B_1)&e(A_1,B_2)&\cdots&e(A_1,B_n)\\
e(A_2,B_1)&e(A_2,B_2)&\cdots&e(A_2,B_n)\\
\vdots&\vdots&\ddots&\vdots\\
e(A_n,B_1)&e(A_n,B_2)&\cdots&e(A_n,B_n)\end{bmatrix}$$

$$\det(M)=\sum\limits_{S:A\rightarrow B}(-1)^{\pi(p)}\prod\limits_{i=1}^n w(S_i)$$
其中 $\sum\limits_{S:A\rightarrow B}$ 表示满足上文要求的 $A\rightarrow B$ 的每一组不相交路径 $S$。

意思是什么呢？就是上面这个矩阵的行列式的值，等于所有可能的 $w(S)$ 的带符号和，如果 $w$ 恒为 $1$，那么就是所有合法的 $S$ 集合的带符号和。

证明

首先，由定义得：

$$
\begin{aligned}
\det(M)&=\sum\limits_{P:A\rightarrow B}(-1)^{\pi(p)}\prod\limits_{i=1}^n w(S_i) \\
&=\sum\limits_{P:A\rightarrow B}(-1)^{\pi(p)}\prod\limits_{i=1}^n \sum_{P:a_i \to b_{p_i}}w(P) \\
\end{aligned}
$$

我们可以发现通过乘法分配律展开之后，因为 $w(P)$ 中的 $P$ 可以是多条路径组合，所以就是 $\prod\limits_{i=1}^n \sum_{P:a_i \to b_{p_i}}w(P)=\sum_{P:p}w(P)$，表示所有起点集合和终点集合与 $p$ 相同的 $P$ 的和，$P$ 是任意路径组，不要求不相交。

所以：

$$
\begin{aligned}
\det(M)
&=\sum\limits_{p:A\rightarrow B}(-1)^{\pi(p)}\prod\limits_{i=1}^n \sum_{P:a_i \to b_{p_i}}w(P) \\
&=\sum\limits_{p:A\rightarrow B}(-1)^{\pi(p)}\sum_{P:p}w(P) \\
&=\sum\limits_{P:A\rightarrow B}(-1)^{\pi(p)}\prod_{i=1}^n w(P_i)\\
\end{aligned}
$$

到此时 $P$ 也是任意路径组合。

我们设 $U$ 是要求的集合，而 $V$ 是有相交的路径集合，那么有：

$$
\begin{aligned}
&\sum\limits_{P:A\rightarrow B}(-1)^{\pi(p)}\prod_{i=1}^n w(P_i)\\
=&\sum\limits_{U:A\rightarrow B}(-1)^{\pi(u)}\prod_{i=1}^n w(U_i)+\sum\limits_{V:A\rightarrow B}(-1)^{\pi(v)}\prod_{i=1}^n w(V_i)\\
\end{aligned}
$$

我们又可以发现：若 $P$ 存在两条相交的路径 $a_i \to u \to a_j,b_i \to u \to b_j$，那么一定有 $a_i \to u \to b_j,b_i \to u \to a_j$，这两个项符号相反，绝对值相同。（按照之前行列式交换两行的方式证明即可）

所以 $\sum\limits_{V:A\rightarrow B}(-1)^{\pi(v)}\prod_{i=1}^n w(V_i)=0$，那么 $\sum\limits_{P:A\rightarrow B}(-1)^{\pi(p)}\prod_{i=1}^n w(P_i)
=\sum\limits_{U:A\rightarrow B}(-1)^{\pi(u)}\prod_{i=1}^n w(U_i)$，由此，结论得证。

运用

如果是一张普通的有向无环图，是不可以直接运用 LGV 引理的，因为带符号会让问题变得不好处理。

当然，有一些特殊的图满足以下性质的话我们就可以使用 LGV 引理直接求：

• 需要求出对于每个 $i$ 找到 $a_i \to b_i$ 的路径，并且路径两两没有公共点。

• 如果排列 $p$ 中有逆序对，那么这些路径一定有交点。

根据 $S$ 的定义，这些情况不会被统计到答案中，因此我们可以直接计算。

例题参见 LGV 引理，这个图就满足上述两种情况，可以直接利用行列式求解。

代码如下：

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

#include<bits/stdc++.h>
#define ll long long
#define mod 998244353
#define N 605
#define M 5000005
using namespace std;
ll T,n,m,a[N],b[N],i,j,jc[M],inv[M],maps[N][N];
inline ll qmi(ll a,ll b,ll p){
	ll res = 1%p,t=a;
	while(b){
		if(b&1) res=res*t%p;
		t=t*t%p;
		b>>=1;
	}
	return res;
}
inline ll C(ll n,ll m){
	if(n<m) return 0;
	return jc[n]*inv[m]%mod*inv[n-m]%mod;
}
inline ll solve(ll n){
	ll i,j,k,has=0,ans=1,res;
	for(i=1;i<=n;i++){
		for(j=i+1;j<=n;j++){
			if(maps[j][i]>maps[i][i]) swap(maps[i],maps[j]),has^=1;
			res = maps[j][i]*qmi(maps[i][i],mod-2,mod)%mod;
			for(k=i;k<=n;k++) maps[j][k]=maps[j][k]-maps[i][k]*res,maps[j][k]=(maps[j][k]%mod+mod)%mod;
		}
	}
	if(i<=n) return 0;
	for(i=1;i<=n;i++) ans=ans*maps[i][i]%mod;
	if(has) ans=(mod-ans)%mod;
	return ans;
}
int main(){
	ios::sync_with_stdio(false);
	jc[0] = 1;
	for(i=1;i<=5e6;i++) jc[i]=jc[i-1]*i%mod;
	inv[5000000] = qmi(jc[5000000],mod-2,mod);
	for(i=5e6;i>=1;i--) inv[i-1]=inv[i]*i%mod;
	cin>>T;
	while(T--){
		cin>>n>>m;
		for(i=1;i<=m;i++) cin>>a[i]>>b[i];
		for(i=1;i<=m;i++) for(j=1;j<=m;j++) maps[i][j]=C(b[j]-a[i]+n-1,n-1);
		cout<<solve(m)<<endl;
	}
	return 0;
}
/*
Input:
3
3 2
1 2
2 3
5 2
1 3
3 5
10 5
3 5
4 7
5 8
7 9
9 10

Output:
3
155
2047320
*/

BEST 定理

分析

对于一道包含欧拉回路的有向图，有 BEST 定理：

$$
cnt=T \prod_{i=1}^n(ot_i-1)
$$

其中 $ot_i$ 表示第 $i$ 个点的出度，$T$ 表示这幅图的内向生成树的个数（取任何一个节点为根都可以，答案都相等，这是包含欧拉回路的有向图的性质）。

注意：这是把循环同构的欧拉回路去除掉的答案！

而这道题要求从 $1$ 开始，从 $1$ 结束，所以我们需要乘上 $ot_1$，证明在下面。

证明

偏感性。

首先考虑对于任意一条欧拉回路保留除去起始点 $s$ 剩下点的最后一条出边，那么这些出边一定构成一个以 $s$ 为根的内向树。

如果成环，那就说明形成了一个 $\rho$ 形，交点处因为入度等于出度，一定还会有一条出边，于是不可能成环。

接下来考虑反向构造，如果固定了起点和终点为 $s$，那么从 $s$ 到外面走有 $ot_s!$ 种走法，对于每个点，因为固定了一条出边，所以都有 $(ot_i-1)!$ 种走法，乘起来即可。

如果不固定起点，在上面的方案中，每种欧拉回路都因为 $s$ 的最后一条出边是任意选择的，所以除掉 $ot_s$ 就变成了文章一开始的形式。（也可以理解为 $s$ 把任何路径划分成了 $ot_s$ 段，但是对于整幅图来说这 $ot_s$ 段循环同构，所以需要除一下 $ot_s$）

实现

首先需要判断这个有向图有没有欧拉回路，判定条件为去掉独立点之后从 $1$ 开始能不能遍历到所有点，并且不是独立点的点的入度必须等于出度。

如果含有欧拉回路，用矩阵树定理求 $T$ 的值即可。

注意，这里求 $T$ 如果选第 $i$ 行和第 $i$ 列删掉，一定要保证 $i$ 不是独立点，代码选的是 $i=1$，当然你也可以选择以 $3$ 为根的内向树，但是要判断 $3$ 是不是独立点。

代码

代码如下：

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

#include<bits/stdc++.h>
#define ll long long
#define N 105
#define M 200005
#define mod 1000003
using namespace std;
vector<ll> op[N];
ll T,n,m,i,j,x,y,a[N][N],ans,jc[M],in[N],ot[N],vis[N],vit;
ll qmi(ll a,ll b,ll p){
	ll res = 1%p,t=a;
	while(b){
		if(b&1) res=res*t%p;
		t=t*t%p;
		b>>=1;
	}
	return res;
}
void dfs(ll x){
	vis[x] = 1;
	for(ll i=0;i<op[x].size();i++){
		vit++;
		if(vis[op[x][i]]) continue;
		dfs(op[x][i]);
	}
}
bool check(){
	dfs(1);
	if(vit!=m) return 0;
	for(ll i=1;i<=n;i++) if(in[i]!=ot[i]) return 0;
	return 1;
}
ll solve(){
	ll ans=1,has=0;
	for(ll i=1;i<=n;i++){
		if(i==1||!vis[i]) continue;
		for(ll j=i+1;j<=n;j++){
			if(j==1||!vis[j]) continue;
			if(!a[i][i]) swap(a[i],a[j]),has^=1;
			ll res = a[j][i]*qmi(a[i][i],mod-2,mod)%mod;
			for(ll k=i;k<=n;k++){
				if(k==1||!vis[k]) continue;
				a[j][k] = (a[j][k]-a[i][k]*res%mod+mod)%mod;
			}
		}
	}
	for(ll i=1;i<=n;i++) if(i!=1&&vis[i]) ans=ans*a[i][i]%mod;
	if(has) ans=(mod-ans)%mod;
	return ans;
}
int main(){
	ios::sync_with_stdio(false);
	jc[0] = 1;
	for(i=1;i<=2e5;i++) jc[i]=jc[i-1]*i%mod;
	cin>>T;
	while(T--){
		vit=0,m=0;
		cin>>n;
		for(i=1;i<=n;i++){
			cin>>x;
			while(x--) cin>>y,op[i].push_back(y),a[y][i]--,a[i][i]++,in[y]++,ot[i]++,m++;
		}
		if(check()){
			ans = solve();
			for(i=1;i<=n;i++){
				if(!vis[i]) continue;
				if(i==1) ans=ans*jc[ot[i]]%mod;
				else ans=ans*jc[ot[i]-1]%mod;
			}
			cout<<ans<<endl;
		}
		else cout<<0<<endl;
		for(i=1;i<=n;i++) op[i].clear(),in[i]=0,ot[i]=0,vis[i]=0;
		for(i=1;i<=n;i++) for(j=1;j<=n;j++) a[i][j]=0;
	}
	return 0;
}

Matrix-Tree 定理

矩阵树定理。

定理内容

对于一张无向图，定义它的度数矩阵为：

$$
D=\begin{bmatrix}
d_1 & & & & \\
& d_2 & & & \\
& & \ddots & & \\
& & & d_{n-1} & \\
& & & & d_n \\
\end{bmatrix}
$$

$d_i$ 为 $i$ 节点的度数，同时定义它的邻接矩阵为 $E_{i,j}$ 为 $i \to j$ 的边数，那么它的 Kirchhoff 矩阵 $L$ 为 $E-D$。

那么我们在行列式中学到了代数余子式 $M_r$ 表示删去第 $r$ 行和第 $r$ 列的元素之后行列式的值，那么有性质 $L_1=L_2=\dots=L_n=ans$，其中 $ans$ 就是这张无向图中的生成树的个数。

于是我们只需要在 $O(n^3)$ 内求出一个行列式的值即可。

并且如果边有边权（不一定是整数，有可能是向量，参见 P5296 [北京省选集训2019] 生成树计数 题解），那么矩阵树定理可以求出所有生成树的边权乘积的和。

对于一张有向图，根节点为 $r$，我们要找到以 $r$ 为根的外向树（边从父亲指向子节点，如果是内向树，将边反向即可），我们定义度数矩阵：

$$
D=\begin{bmatrix}
d_1 & & & & \\
& d_2 & & & \\
& & \ddots & & \\
& & & d_{n-1} & \\
& & & & d_n \\
\end{bmatrix}
$$

$d_i$ 为 $i$ 点的入度，定义邻接矩阵为 $m_{i,j}$ 表示 $i \to j$ 的边数，那么它的 Kirchhoff 矩阵 $L$ 也为 $E-D$。

因为指定了根节点，所以我们只能输出代数余子式 $L_{root}$ 的值，其它代数余子式的值就是以 $i$ 为根的时候，外向树的个数。

证明

略，详见 OI-wiki。

代码

代码如下，当 $t=0$ 的时候代表是无向图，当 $t=1$ 的时候代表是有向图：

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

#include<bits/stdc++.h>
#define ll long long
#define N 605
using namespace std;
ll n,m,t,x,y,z,a[N][N],i,j,k,mod=1e9+7,ans=1,has,res;
int main(){
	ios::sync_with_stdio(false);
	cin>>n>>m>>t;
	for(i=1;i<=m;i++){
		cin>>x>>y>>z;
		x--,y--;
		if(x==y) continue;
		if(t==0){
			if(x>=1&&y>=1) a[x][y]-=z,a[y][x]-=z;
			if(x>=1) a[x][x]+=z;
			if(y>=1) a[y][y]+=z;
		}
		else{
			if(x>=1&&y>=1)a[x][y]-=z;
			if(y>=1) a[y][y]+=z;
		}
	}
	n--;
	for(i=1;i<=n;i++) for(j=1;j<=n;j++) a[i][j]=(a[i][j]%mod+mod)%mod;
	for(i=1;i<=n;i++){
		for(j=i+1;j<=n;j++){
			while(a[i][i]){
				res = a[j][i]/a[i][i];
				for(k=i;k<=n;k++) a[j][k]=(a[j][k]-a[i][k]*res%mod+mod)%mod;
				swap(a[i],a[j]),has^=1;
			}
			swap(a[i],a[j]),has^=1;
		}
	}
	if(i<=n){
		cout<<0<<endl;
		return 0;
	}
	for(i=1;i<=n;i++) ans=ans*a[i][i]%mod;
	if(has) ans=(mod-ans)%mod;
	cout<<ans<<endl;
	return 0;
} 
/*
Input:
2 998244353
1 4
1 5

Output:
1
*/

扩展

如果我们要求有向图以不同节点为根的外向生成树的数量总和，那么需要我们拿出基尔霍夫矩阵，然后把这个矩阵任意一行改成全 $1$。

注意：是任意一行，列上的数不用管它，证明尚不明确，先记住这个结论吧。

例题：D - xfh1我回来了 | NKOJ，代码如下：

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

#include<bits/stdc++.h>
#pragma GCC optimize("Ofast")
#define ll long long
#define mod 1000000007
#define N 17
#define M 505
using namespace std;
ll n,m,x[M],y[M],z[M],i,j,f[1<<N],a[N][N],ans;
inline ll solve(){
	ll has = 0,ans = 1,res,i,j,k;
	for(i=1;i<=n;i++){
		for(j=i+1;j<=n;j++){
			while(a[i][i]){
				res = a[j][i]/a[i][i];
				for(k=i;k<=n;k++) a[j][k]=(a[j][k]-a[i][k]*res%mod+mod)%mod;
				swap(a[i],a[j]),has^=1;
			}
			swap(a[i],a[j]),has^=1;
		}
	}
	if(i<=n) return 0;
	for(ll i=1;i<=n;i++) ans=ans*a[i][i]%mod;
	if(has) ans=(mod-ans)%mod;
	return ans;
}
int main(){
	ios::sync_with_stdio(false);
	cin>>n>>m;
	for(i=1;i<=m;i++) cin>>x[i]>>y[i]>>z[i];
	for(i=0;i<(1<<(n-1));i++){
		memset(a,0,sizeof(a));
		for(j=1;j<=m;j++) if((i>>z[j])&1) a[x[j]][y[j]]--,a[y[j]][y[j]]++;
		for(j=1;j<=n;j++) a[1][j]=1;
		f[i] = solve();
	}
	for(i=0;i<(1<<(n-1));i++){
		for(j=(i&(i-1));;j=(i&(j-1))){
			f[i] = (f[i]-f[j]+mod)%mod;
			if(j==0) break;
		}
		for(j=0;j<n-1;j++) if(!((i>>j)&1)) break;
		ans = (ans+f[i]*j)%mod;
	}
	cout<<ans<<endl;
	return 0;
}