前言：由於作者未系統過學習線性代數，故下文肯定有不嚴謹的成分，不過應付演算法競賽是綽綽有餘了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;
}