蒟蒻 Nanjo_Qi 前天考了一次試……第一題就華麗麗地爆零了。
解一次方程組我會啊，但是解一千個有百來八十個未知數的……棄了棄了orz。
考完了才知道有高斯消元這個神奇的東西，於是就去簡單了解了一下。
高斯消元法是線性代數規划中的一個算法，可用來為線性方程組求解，還可以求出矩陣的秩，以及求出可逆方陣的逆矩陣。消元法就是將方程組中的一方程的未知數用含有另一未知數的代數式表示，並將其代人到另一方程中，這就消去了一未知數，得到一解；或將方程組中的一方程倍乘某個常數加到另外一方程中去，也可達到消去一未知數的目的（其實就是中小學最常用的解二元，三元一次方程組的思想嘛）。我們把構成n個未知數的，含n個方程的方程組，亂搞一下，一個一個未知數求解，求出一個就再代入回去求其它的，這樣就能解出所有的來啦（也有可能會無解或有多組解）！
那該怎么搞呢？
我們舉個栗子，有這樣一個三元一次方程組：
5x + 4y + 7z = 8
3x + 9y – 1z = -2
8x + 3y + 2z = 9
我們默認每一列上的未知數是相同的，如果沒有就補個0，這樣拋除未知數，我們就得到了一個n*(n+1)的矩陣，然后以后我們規定未知數用x1(x)，x2(y)，x3(z)等表示：
5 4 7 8
3 9 -1 -2
8 3 2 9
下面到了算法的核心：
如果我們正在求解xi，那么首先我們選取一個方程j>=i，使a[j][i]不為0，然后把這個方程和第i行方程互換一下。隨后對於每一行j>i，我們把這個方程通過加減方程i的倍數，將其中的xi的系數化為0，即進行消元操作。
這個地方建議手動拿例子模擬一下，代碼長這樣：
1 int p = i; 2 for(int j=i; j<=n; ++j) if( fabs(a[p][i]) ) { p = j; break; } //選取一個方程使a[j][i]不為0 3 for(int j=i; j<=n+m; ++j) swap(a[i][j], a[p][j]); 4 for(int j=i+1; j<=n; ++i) { //向下消元 5 if( !a[j][i] ) continue; //如果系數已經是0了，不用消元 6 for(int k=i+1; k<=n+m; ++k) 7 a[j][k] -= a[i][k] / a[i][i] * a[j][i]; 8 a[j][i] = 0; 9 }
然后外面套一層循環i = 1 ~ n，完成每一個消元的操作，這時對矩陣的操作就完成啦。
這時這個矩陣被搞成了這樣，去掉常數那一邊的話會是一個嚴格上三角矩陣（保留兩位小數）：
在高斯消元算法中，當且僅當它是嚴格上三角矩陣時，才是有唯一解的。
5.00 4.00 7.00 8.00
0.00 6.60 -5.20 -6.80
0.00 0.00 -11.88 -7.30
也就是：
5x + 4y + 7z = 8
6.6y – 5.2z = -6.8
-11.88z = -7.3
此時從后往前看，我們容易解出z！然后就容易解出y！然后代入解出x，就做完了。
其實還是很簡單的，這一步可以這樣解決：
1 x[n] = a[n][n+i] / a[n][n]; //未知數x[n] 2 for(int j=n-1; j>=1; --j) { //往前推 3 x[j] = a[j][n+i]; 4 for(int k=j+1; k<=n; ++k) 5 x[j] -= a[j][k] * x[k]; 6 x[j] /= a[j][j]; 7 }
這樣，高斯消元法求解線性方程組的時間是O(n^3)的。
那無解的時候怎么辦呢？
只需要判斷是否出現某時對xi進行消元時，能否找到找到一個xi的系數不為零就可以了，就是：
1 if( fabs(a[p][i])<eps ) return 0;
下面是考試題目：
Pear 上數學課（equation.pas/c/cpp）
題目描述：
Pear 被一群長得很像的一次方程組圍住了，快去拯救他！
Pear 最近的作業之一，就是練習解一次方程組，老師給了他 m 個方程組，每個方程組未知數個數一樣（n 個），未知數的系數一樣，唯一不同的是每個方程的常數項。
比如 n=2，m=3，pear 手上有這樣的 3 個方程組：
𝑥 + 𝑦 = 3 𝑥 + 𝑦 = 20 𝑥 + 𝑦 = 4
𝑥 − 𝑦 = 5 𝑥 − 𝑦 = 10 𝑥 − 𝑦 = 0
那他會分別解得：
𝑥 = 4 𝑥 = 15 𝑥 = 2
𝑦 = −1 𝑦 = 5 𝑦 = 2
現在這個任務就交給你了！
保證方程組有解，允許誤差范圍在 1e-3 以內。
輸入格式：
第一行兩個正整數 n，m；
之后 n 行每行 n 個整數，第 i 行的第 j 個整數 ai,j 表示第 i 個方程中 xi 的系數；
之后 m 行每行 n 個整數，第 i 行的第 j 個整數 bi,j 表示第 m 個方程組的第 j 個方程的常數項（約定方程的左邊寫未知數，右邊寫常數項）。
輸出格式：
輸出 m 行 n 個實數，第 i 行第 j 個實數表示方程組 i 中的 xj 的值。
輸入樣例：
2 3
1 1
1 -1
3 5
20 10
4 0
輸出樣例：
4 -1
15 5
2 2
數據規模：
對於 10%的數據 n=2；
對於 30%的數據 n,m<=20；
對於 50%的數據 n<=50,m<=100；
對於 100%的數據 n<=100,m<=1000。
HINT
本題包含 special judge。
有所不同的是題目中含有m個方程組，一個個做的時間是O(n^3*m)，會T。
很方便的是每一組數據所有方程組系數是一樣的，只有常數不一樣，而可以發現在各種對矩陣的操作中實際上常數項並不會影響什么，只是隨之加加減減而已，所以不妨將這些常數全放到矩陣對應位置，我們一起操作，這樣的時間是O(n^2*(n+m))的，大概是可過的范圍。
最后貼一個代碼：
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 using namespace std; 5 6 const double eps = 1e-10; 7 const int maxn = 100 + 10; 8 const int maxm = 1000 + 10; 9 int n, m; double a[maxn][maxn+maxm], x[maxn]; 10 11 inline double Abs(double x) { return x>=0 ? x : -x; } 12 13 inline int Gaussian() { 14 for(int i=1; i<=n; ++i) { 15 int p = i; 16 for(int j=i+1; j<=n; ++j) 17 if( Abs(a[j][i]) > Abs(a[p][i]) ) p = j; 18 //這樣選可以有所避免精度問題 19 if( Abs(a[p][i])<eps ) return 0; 20 for(int j=i; j<=n+m; ++j) swap(a[i][j], a[p][j]); 21 for(int j=i+1;j<=n;j++) { 22 if( a[j][i]==0 ) continue; 23 for(int k=i+1; k<=n+m; ++k) 24 a[j][k] -= a[i][k] / a[i][i] * a[j][i]; 25 a[j][i] = 0; 26 } 27 } 28 return 1; 29 } 30 31 int main() { 32 scanf("%d%d", &n, &m); 33 for(int i=1; i<=n; ++i) 34 for(int j=1; j<=n; ++j) 35 scanf("%lf", &a[i][j]); 36 for(int i=1; i<=m; ++i) 37 for(int j=1; j<=n; ++j) 38 scanf("%lf", &a[j][n+i]); 39 40 int flag = Gaussian(); 41 if( !flag ) { printf("-1\n"); return 0; } 42 43 for(int i=1; i<=m; ++i) { //對於m組方程組分別求解 44 x[n] = a[n][n+i] / a[n][n]; 45 for(int j=n-1; j>=1; --j) { 46 x[j] = a[j][n+i]; 47 for(int k=j+1; k<=n; ++k) 48 x[j] -= a[j][k] * x[k]; 49 x[j] /= a[j][j]; 50 } 51 for(int j=1; j<=n; ++j) 52 printf("%.9lf ", x[j]); 53 printf("\n"); 54 } 55 return 0; 56 }