高斯消元

转载自https://www.luogu.com.cn/blog/pks-LOVING/p3389-mu-ban-gao-si-xiao-yuan-fa

emmmm 这个消元方式其实严格来说是可行性算法……而不是优化性算法……不过话说由于我太蒟蒻了，所以并不知道什么更优的算法（#滑稽）

嗯，其实这个算法是 O(n^3)O(n3) 的算法，需要一些矩阵及行列式的知识……那么由本蒟蒻来记录一下这个算法吧！

那么假设有一个线性方程组是长这样的：

egin{Bmatrix} 3x & + & 2y &+& z &=&10 \5x & + & y &+& 6z &=&25 \2x & + & 3y &+& 4z &=&20end{Bmatrix}⎩⎨⎧​3x5x2x​+++​2yy3y​+++​z6z4z​===​102520​⎭⎬⎫​

emmmemmm 这就是一个很简单的三元一次方程，让我们想想常规方法该怎么做（先不谈 codecode ）

初中老师说过：我们可以加减消元或者代入消元,但是我们需要在程序里实现的时候，需要一种有规律可循的算法。所以我们选择加减消元，但用代入消元回带。

整体思路就是我们可以先在某一个式子里，用这个式子的 xx 消去其他式子里的 xx ,然后在剩下的两个式子里再选择一个式子里的 yy ，用这个 yy 消去最后剩下的式子里的 yy 。那么现在最后一个方程里就只有一个未知数 zz 了。倘若 zz 的系数是 11 ，那么我们就可以直接得出答案来了（别觉得这句话是废话）。

比如刚才这个方程，我们用第二个式子去消1、3式里的 xx ：

egin{Bmatrix} 0 imes x & + & frac{7}{5}y &+& (-frac{13}{5}z) &=&-5 \5x & + & y &+& 6z &=&25 \0 imes x & + & frac{13}{5}y &+& frac{8}{5}z &=&10end{Bmatrix}⎩⎨⎧​0×x5x0×x​+++​57​yy513​y​+++​(−513​z)6z58​z​===​−52510​⎭⎬⎫​

整理之后再用第三个式子里的 yy 消去第一个式子里的 yy (注意，由于第二个式子作为消元用式，所以接下来的运算不再考虑二式):

egin{Bmatrix}0 imes y &+& (-frac{225}{65}z) &=&-frac{135}{13} \ frac{13}{5}y &+& frac{8}{5}z &=&10end{Bmatrix}{0×y513​y​++​(−65225​z)58​z​==​−13135​10​}

那么我们发现在 11 式中只剩下一个未知数了，那么就可解得：z=3z=3

带回三式里解出y=2y=2

再将 xx 、 yy 带回最早被消掉的二式里，解得x=1x=1

好像这个方法再数学逻辑上讲是特别麻烦的，但是却是一个通用性强的方法 qwqqwq

那么放在程序实现上来讲，我们可以先用一个 n imes (n+1)n×(n+1) 的矩阵来记录每一个式子的系数以及结果。譬如上式就可以用矩阵表示成这样：

egin{Bmatrix} 3& 2 & 1 &|& 10 \5 & 1 & 6 &|& 25 \2 & 3 & 4 &|&20end{Bmatrix}⎩⎨⎧​352​213​164​∣∣∣​102520​⎭⎬⎫​

左边记录系数，右边记录每个式子的结果。

那么首先我们需要用第一列中（所有的 xx 中）系数最大的来消其他两个式子。而这里为了方便起见，我们将这个选中的系数置为 11 ，方便上例中地不断带回原式的操作（这样在回带的时候就可以不考虑原本的系数了）。

由于最多也只能用 doubledouble 型存储，所以必然会有精度误差。但如果我们每次都选用最大系数的来消掉其他系数，就可以最大程度地来减小误差。以下是一种不严谨地、适合意会的证明（选读）：

---

---

假设我们现在在处理第 nn 个未知数，此时在众多的未知数 nn 中，他们的系数分别是 k_1 k_2 k_3 k_ 4k1​k2​k3​k4​ …… kmkm ,那么考虑，在选完 k_iki​ 之后，下面我们要进行的是把 k_iki​ 消成 11 。那么此时对于第 ii 行的其他的系数以及结果我们都要除以 k_iki​ 。

之后呢？之后我们要进行的操作是用这个式子来消掉其他式子里的该未知数啊 qwqqwq 。如果要这么操作肯定会让其他式子别的未知数的系数，减去当前式子的别的未知数的系数乘上某个值（事实上假设选择含 k_iki​ 的式子,则对于每个式子 jj 而言，每个系数减去当前系数的倍数，这个倍数应该为 k_jkj​ ）

那么这样看来，对于当次用来消元的式子的每个系数 q_{i1}q_{i2}q_{i3}q_{i4}qi1​qi2​qi3​qi4​ …… q_{iw}qiw​ (假设当前元的系数是 q_{i1}qi1​ )而言，对于每一个其他式子的该项系数 q_{jw}qjw​ ,都需要让 q_{jw}qjw​ 变成

q_{jw}-frac{q_{j1}}{q_{i1}} imes q_{iw}qjw​−qi1​qj1​​×qiw​

那么我们观察这个式子， q_{i1}qi1​ 越大， frac{q_{j1}}{q_{i1}}qi1​qj1​​ 期望越小，那么我们考虑，这个值越小，我们就约可以把它看作一个“基本单位”。从而我们就使得减出来的值失精程度越低，最后即可保证数据是从期望上来讲最精确。

嗯，讲的很麻烦，大家挑重点看吧（或者只看最后一个自然段）

(逃

---

---

在置为 11 之后，我们需要来用这个式子去消其他的式子（别忘了每个式子的结果也要消）。那么在最后，我们只需要将这个矩阵的最右下角（也就是最后一个元的实际值）不断回带即可。

代码长这个样子：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
double map[111][111];
double ans[111];
double eps=1e-7;
int main(){
    int n;
    cin>>n;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&map[i][j]);
    for(int i=1;i<=n;i++){
        int r=i;
        for(int j=i+1;j<=n;j++)
            if(fabs(map[r][i])<fabs(map[j][i]))
                r=j;//find_the_biggest_number_of_the_first_column（at present) 
        if(fabs(map[r][i])<eps){
            printf("No Solution");
            return 0;
        }
        if(i!=r)swap(map[i],map[r]);//对换一行或一列,属于找最大当前系数的其中一步。（这样就可以只处理当前行的系数啦！） 
        double div=map[i][i];
        for(int j=i;j<=n+1;j++)
            map[i][j]/=div;
        for(int j=i+1;j<=n;j++){
            div=map[j][i];
            for(int k=i;k<=n+1;k++)
                map[j][k]-=map[i][k]*div;
        }
    }
    ans[n]=map[n][n+1];
    for(int i=n-1;i>=1;i--){
        ans[i]=map[i][n+1];
        for(int j=i+1;j<=n;j++)
            ans[i]-=(map[i][j]*ans[j]);
    }//回带操作
    for(int i=1;i<=n;i++)
        printf("%.2lf
",ans[i]);
}