https://en.wikipedia.org/wiki/Gaussian_elimination
https://en.wikipedia.org/wiki/Invertible_matrix
https://en.wikipedia.org/wiki/Determinant
高斯消元可以做3件事
系数矩阵行列式非0（满秩）
那么方程有解且唯一
系数矩阵行列式为0（不满秩）
那么方程要么无解，要么无穷多解
高斯消元
swap(a, b)
a和b可以是数组，即使数组不可以相互赋值
解方程 3 种情况
高斯消元不仅可以用于解方程，还可以用于行列式求值。
但是行列式满足的性质，和方程组略有不同
性质：
一行所有数字同时乘以k
第i行所有数字都减去第j行对应的数字
需要考虑精度问题，所以在把非零行换到第行时，一般换绝对值最大的。
参考代码 Luogu P3389
#include <bits/stdc++.h> using namespace std; int n; double a[120][120]; int main() { cin >> n; for (int i = 0; i < n; i++) { for (int j = 0; j <= n; j++) { cin >> a[i][j]; } // a[i][0] * x[0] + a[i][1] * x[1] + ... + a[i][n-1] * x[n-1] = a[i][n] } for (int i = 0; i < n; i++) { int t = i; for (int j = i; j < n; j++) { if (abs(a[t][i]) < abs(a[j][i])) { // 换绝对值最大的一行到第i行。 t = j; } } swap(a[t], a[i]); // 不同于赋值，swap可以交换两个数组。 double u = a[i][i]; if (fabs(u) < 1e-6) { printf("No Solution\n"); // 这里可能无解，或者无穷多解。 return 0; } for (int j = 0; j <= n; j++) { a[i][j] /= u; // 第i行 /= a[i][i]，使得a[i][i] = 1。 } for (int j = 0; j < n; j++) { if (i == j) { continue; } // 试图消去除自己之外的所有行。 u = a[j][i]; for (int k = 0; k <= n; k++) { a[j][k] -= u * a[i][k]; } } } for (int i = 0; i < n; i++) { printf("%.2f\n", a[i][n]); // 每一行只有a[i][i] = 1，所以 x[i] = a[i][n]。 } return 0; }
最容易实现的版本，只需要实现逆元即可。
参考代码 Luogu P4783
// 需要O2 #include <bits/stdc++.h> using namespace std; int n, a[401][801]; const int p = 1000000007; // 一定要const int，或者define，否则会超时！！ int pw(int x, int y) { int re = 1; for (; y > 0; y >>= 1) { if (y & 1) { re = (long long)re * x % p; } x = (long long)x * x % p; } return re; } int main() { cin >> n; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { cin >> a[i][j]; } a[i][n + i] = 1; // 右侧加一个n x n单位矩阵。 } for (int i = 0; i < n; i++) { int t = i; for (int j = i; j < n; j++) { if (a[j][i] > 0) { // 换任意非0行即可。 t = j; } } swap(a[t], a[i]); int u = pw(a[i][i], p - 2); if (u == 0) { printf("No Solution\n"); return 0; } for (int j = 0; j < 2 * n; j++) { // 变换时，范围是2n a[i][j] = (long long)a[i][j] * u % p; // u是a[i][i]的逆元。 } for (int j = 0; j < n; j++) { if (i == j) { continue; } u = a[j][i]; for (int k = 0; k < 2 * n; k++) { // 变换时，范围是2n a[j][k] = (a[j][k] - (long long)u * a[i][k]) % p; if (a[j][k] < 0) { a[j][k] += p; } } } } for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { printf("%d%c", a[i][j + n], j == n - 1 ? '\n' : ' '); } } // 输出右侧n x n的矩阵。 return 0; }
可以用bitset或者压位来优化常数。
参考代码
Luogu P3389
实数 高斯消元，模板题。
bzoj 1013
实数 高斯消元。
Luogu P4783
整数取模 高斯消元。
spoj JZPLIT
spoj JZPLIT2
bitset优化高斯消元解异或方程
https://github.com/wwwwodddd/Zukunft/blob/master/luogu/P4783.cpp
P4783 【模板】矩阵求逆
https://github.com/wwwwodddd/Zukunft/blob/master/luogu/P3389.cpp
【模板】高斯消元法
高斯消元
枚举第i行：
a[i][i] 可能是 0
找第j行，使得a[j][i] != 0 （对于实数的版本，把绝对值最大的换到第i行可以减少误差）
交换第i行和第j行（必须让a[i][i]不为0）
第i行同时乘以 a[i][i] 的逆元 （除以a[i][i]）（让a[i][i]变成1，这里这样写可以让代码更简单
将所有不是第i行的 第j行，减去若干倍的 第i行，使得a[j][i]变成0
（如果 只消去第i行之后的，那么在最后计算答案是，还需要代回求值）
x, y (x-x[1])^2+(y-y[1])^2 == (x-x[2])^2+(y-y[2])^2 - 2 x x[1] + x[1]^2 - 2 y y[1] + y[1]^2 = - 2 x x[2] + x[2]^2 - 2 y y[2] + y[2]^2 相邻2个方程相减，消去二次项
P3389 【模板】高斯消元法
P4783 【模板】矩阵求逆