一次线性方程组 高斯消元
Link_Cut_Y
·
·
算法·理论
高斯消元原理
高斯消元用来解如下形式的方程组：
\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}
首先将所有系数写成系数矩阵：
\begin{bmatrix}
a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} & b_{1}\\
a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} & b_{2}\\
\vdots & \vdots & \ddots & \vdots & \vdots\\
a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} & b_{n}
\end{bmatrix}
接下来可以进行初等行列变换，将矩阵消成下三角矩阵。类似这样：
\begin{bmatrix}
a'_{1, 1} & a'_{1, 2} & \cdots & a'_{1, n} & b'_{1}\\
0 & a'_{2, 2} & \cdots & a'_{2, n} & b'_{2}\\
0 & 0 & \cdots & a'_{3, n} & b'_{3}\\
\vdots & \vdots & \ddots & \vdots & \vdots\\
0 & 0 & \cdots & a'_{n, n} & b'_{n}
\end{bmatrix}
然后最后一行就表示 a_{n, n} x_n = a_{n, n + 1}，可以解出 x_n。
然后往回带，解出 x_{n - 1}，以此类推即可。时间复杂度 O(n ^ 3)。
void gauss() {
rep(i, 1, n) {
rep(j, i, n) if (fabs(a[j][i]) > eps) {
swap(a[j], a[i]); break;
} if (fabs(a[i][i]) <= eps) continue;
rep(j, i + 1, n) {
if (fabs(a[j][i]) <= eps) continue;
db inv = (db)a[j][i] / a[i][i];
rep(k, i, n + 1)
a[j][k] -= inv * a[i][k];
}
}
for (int i = n; i; i -- ) {
f[i] = a[i][n + 1] / a[i][i];
for (int j = i - 1; j; j -- )
a[j][n + 1] -= a[j][i] * f[i];
}
}
高斯消元应用
题意：给定 n 点 m 边无向简单图。要求给每条边重新编号 1 \sim m，使得每条边编号不同且从 1 到 n 经过路径的权值和期望尽可能小。
发现边数很多，所以从点数入手。设 $f_i$ 表示第 $i$ 个点被经过的期望次数，那么有：
$$
\begin{cases}
f_1 = \sum \limits_{(1, v) \in E, v \ne n}^{} \dfrac{f_v}{deg_v} + 1 \\
f_u = \sum \limits_{(u, v) \in E, v \ne n}^{} \dfrac{f_v}{deg_v} + 1 (u \ne 1)\\
\end{cases}
$$
当然，$n$ 号点被我们排除在外，因为到了 $n$ 就停止了。
然后发现，上面的柿子形成了一个 $n - 1$ 元方程。用高斯消元解方程就可以得到 $f_u$。
然后可以得到每条边被经过的期望次数：$g_{(u, v)} = \dfrac{f_u}{deg_u} + \dfrac{f_v}{deg_v}$。根据边被经过的期望次数贪心赋权即可。
其他用高斯消元求解 $dp$ 的例题：[P4321 随机漫游
](https://www.luogu.com.cn/problem/P4321)
- 矩阵求逆
矩阵 $A$ 的逆矩阵定义为 $A^{-1}$，满足 $A^{-1} \times A = I$ 且 $A ^ {-1} \times A = I$。其中 $I$ 表示单位矩阵。
求法：将原矩阵 $A$ 右面拼上一个单位矩阵 $I$。然后对左边一半也就是 $A$ 做初等行列变换消成单位矩阵，同时对右边做左边同样的操作（左边行加右边也加，左边同乘右边也同乘）。最后得到的右半边就是 $A ^ {-1}$。
证明真的不会/kk
```cpp
int gauss() {
rep(i, 1, n) {
rep(j, i, n)
if (a[j][i] != 0) {
swap(a[i], a[j]); break;
}
if (a[i][i] == 0) return 0;
rep(j, i + 1, n) {
if (a[j][i] == 0) continue;
int inv = a[j][i] * qpow(a[i][i]) % mod;
rep(k, i, 2 * n) a[j][k] -= inv * a[i][k] % mod,
a[j][k] = (a[j][k] % mod + mod) % mod;
}
}
rep(i, 1, n) {
int inv = qpow(a[i][i]);
rep(j, 1, 2 * n) a[i][j] = a[i][j] * inv % mod;
}
dep(i, n, 1) {
rep(j, 1, i - 1) {
int inv = a[j][i];
rep(k, 1, 2 * n) a[j][k] -= a[i][k] * inv % mod,
a[j][k] = (a[j][k] % mod + mod) % mod;
}
} return 1;
}
signed main() {
read(n); int B = (10 * n + 5) * sizeof(LL);
rep(i, 0, n + 1) a[i] = (LL*)malloc(B);
rep(i, 1, n) rep(j, 1, n) read(a[i][j]);
rep(i, 1, n) a[i][i + n] = 1; int s = gauss();
if (!s) return puts("No Solution"), 0;
for (int i = 1; i <= n; i ++ , pc('\n'))
rep(j, 1, n) write(' ', a[i][j + n]);
return 0;
}
```