前言：由于作者未系统过学习线性代数，故下文肯定有不严谨的成分，不过应付算法竞赛是绰绰有余了QAQ。
\[\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} \]高斯消元通常用来求解线性方程组，实现方式大致分为两种，普通高斯消元以及约旦-高斯消元。据说后者精度误差更小，故本文先介绍后者。
假设该方程组存在唯一解。由于有 \(n\) 个方程，\(n\) 个未知数，大体思路是通过加减消元的方式，将原方程组的系数化为 \(a_{i,j}=0 \ (i \ne j)\) 的形式，这样就能直接计算方程组的解了。
因此先枚举 \(i\)，然后依次执行以下步骤：
- 首先需要选择一个方程放在第 \(i\) 行。这个方程的意义是将除这个方程以外的其他方程的 \(x_{i}\) 的系数全部化为 \(0\)。这个方程只能从第 \(i\) 行到第 \(n\) 行中选择，因为前面的 \(i-1\) 行的方程的位置已经固定了。我们尽量选择 \(x_{i}\) 的系数的绝对值最大的。据说这样做精度误差更小（
又来），代码实现如下：
lb Mx = 0;
int id = 0;
for(int j = i; j <= n; j++) { //第j个方程
if(Abs(a[j][i]) > Mx) {
Mx = a[j][i];
id = j;
}
}
for(int j = i; j <= n + 1; j++) swap(a[i][j], a[id][j]); //将选好的方程换到第 i 行
- 接着就可以利用第 \(i\) 个方程将其他方程消元了。假设正在消第 \(j\) 个方程，只需要执行这个操作即可： \(\textcircled{j} \leftarrow \textcircled{j}-\textcircled{i} \times \frac{a_{j,i}}{a_{i,i}}\)。
for(int j = 1; j <= n; j++) { //第j个方程
if(j == i) continue;
lb Num = 1.0 * a[j][i] / a[i][i];
for(int k = 1; k <= n + 1; k++) a[j][k] -= a[i][k] * Num;
}
我们就完成了对存在唯一解的方程组的求解。时间复杂度 \(O(n^3)\)。
但有时候一个方程组会出现无解或者无数解的情况，如何判断呢？可以发现，在进行加减消元前，如果找不到主元（关于 \(x_{i}\) 的系数全为 \(0\)），那么一定不是唯一解的情况。那到底是无解还是无数解呢？
这里需要对上面的求解方法稍加修改，记一个 \(r\) 表示当前有效的方程的数量，如果找到能找到主元，就把这个方程放到有效方程的后面。求解完后，如果存在一个无效方程的 \(b\) 值不为 \(0\)，就一定无解（因为此时这些方程的系数肯定全为 \(0\)），否则就是无数解。
以下就是 P2455 [SDOI2006] 线性方程组 的代码：
#include<bits/stdc++.h>
#define int long long
#define N 110
#define lb long double
using namespace std;
int n;
lb a[N][N], b[N]; //a[i][j]表示第i个方程,第j个系数
lb x[N];
lb Abs(lb x) {
return x > 0 ? x : -x;
}
void My_work() { //约旦消元QAQ
int r = 1; //r表示有效的方程
for(int i = 1; i <= n; i++) { //每个未知数
lb Mx = 0;
int id = 0;
for(int j = r; j <= n; j++) { //第j个方程
if(Abs(a[j][i]) > Mx) {
Mx = a[j][i];
id = j;
}
}
if(id == 0) continue;
for(int j = i; j <= n + 1; j++) swap(a[r][j], a[id][j]);
for(int j = 1; j <= n; j++) { //第j个方程
if(a[j][i] == 0 || j == r) continue;
lb Num = 1.0 * a[j][i] / a[r][i];
for(int k = 1; k <= n + 1; k++) a[j][k] -= a[r][k] * Num;
}
/*for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n + 1; j++) printf("%.2Lf ", a[i][j]);
printf("\n");
}
printf("\n");*/
r++;
}r--;
//printf("%lld\n",r);
if(r < n) {
for(int i = r + 1; i <= n; i++) {
if(Abs(a[i][n + 1]) > 1e-9) {
cout << -1 << endl;
return;
}
}
cout << 0 << endl;
return;
}
for(int i = 1; i <= n; i++)
printf("x%lld=%.2Lf \n", i, a[i][n + 1] / a[i][i]);
}
signed main() {
//ios::sync_with_stdio(0);
//cin.tie(0); cout.tie(0);
scanf("%lld", &n);
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n + 1; j++) scanf("%Lf", &a[i][j]);
}
My_work();
return 0;
}