高斯消元原理
高斯消元用来解如下形式的方程组：
\[\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];
}
}
高斯消元应用
- 解线性方程组
这个不用我说了吧。。。
- 在 \(\text{dp}\) 中用于消元
这个用处非常大。举个栗子：P3232 [HNOI2013] 游走
题意：给定 \(n\) 点 \(m\) 边无向简单图。要求给每条边重新编号 \(1 \sim m\)，使得每条边编号不同且从 \(1\) 到 \(n\) 经过路径的权值和期望尽可能小。
\(n \le 500, m \le 125000\)。
发现边数很多，所以从点数入手。设 \(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 随机漫游
- 矩阵求逆
矩阵 \(A\) 的逆矩阵定义为 \(A^{-1}\)，满足 \(A^{-1} \times A = I\) 且 \(A ^ {-1} \times A = I\)。其中 \(I\) 表示单位矩阵。
求法：将原矩阵 \(A\) 右面拼上一个单位矩阵 \(I\)。然后对左边一半也就是 \(A\) 做初等行列变换消成单位矩阵，同时对右边做左边同样的操作（左边行加右边也加，左边同乘右边也同乘）。最后得到的右半边就是 \(A ^ {-1}\)。
证明真的不会/kk
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;
}