数论——高斯消元

---

Surprise!

介绍

何为高斯消元？就是拿来解多元一次方程（组）的一个算法。
什么是多元一次方程？就是有多个未知数，每一个未知数的次数皆为1的方程。多个多元一次方程构成多元一次方程组。如下是一个三元一次方程组：
$$\{\begin{array}{l}3x_1+2x_2+x_3=3\\ 2x_1+x_2+3x_3=6\\ x_1+3x_2+2x_3=9\end{array}$$
容易通过代入消元法或加减消元法获得唯一解：$$\{\begin{array}{l}{x}_1=-1\\ x_2=2\\ x_3=2\end{array}$$

具体题目

题目（只有5个n,m=3的数据，凑合着看看吧）

讲解

还是那个方程：$$\{\begin{array}{l}3x_1+2x_2+x_3=3①\\ 2x_1+x_2+3x_3=6②\\ x_1+3x_2+2x_3=9③\end{array}$$
我们可以直接解：
②$\ast \frac{3}{2}-$①，③$\ast 3-$①得
$$\{\begin{array}{l}3x_1+2x_2+x_3=3①\\ -\frac{1}{2}{x}_2+\frac{7}{2}{x}_3=6④\\ 7x_2+5x_3=24⑤\end{array}$$
所以由③$\ast (-\frac{1}{14})-$④可以得到
$$\{\begin{array}{l}3x_1+2x_2+x_3=3①\\ -\frac{1}{2}{x}_2+\frac{7}{2}{x}_3=6④\\ \frac{27}{7}{x}_3=\frac{54}{7}⑥\end{array}$$
于是就可以得到$x_3=2$。
代入④中，得$-\frac{1}{2}{x}_2+7=6\Rightarrow x_2=2$
代入①中，得$3x_1+4+2=3\Rightarrow x_1=-1$

于是得解。实际上，我们刚才已经模拟了一个高斯消元的过程。虽然这个过程在数学上看来是很麻烦的，但是计算机并不会觉得麻烦，而且这也是一个计算机中的可行算法。模拟手算中的代入消元法，这就是高斯消元的真谛。

矩阵

先不管矩阵是什么，看着就明白了。
如下，原方程$$\{\begin{array}{l}3x_1+2x_2+x_3=3①\\ 2x_1+x_2+3x_3=6②\\ x_1+3x_2+2x_3=9③\end{array}$$可以化为这样一个增广矩阵
$$\left[\begin{array}{ccccc}3& 2& 1& |& 3\\ 2& 1& 3& |& 6\\ 1& 3& 1& |& 9\end{array}\right]$$
按照这个道理，可以将$$\{\begin{array}{l}3x_1+2x_2+x_3=3①\\ -\frac{1}{2}{x}_2+\frac{7}{2}{x}_3=6④\\ \frac{27}{7}{x}_3=\frac{54}{7}⑥\end{array}$$也转化成一个矩阵，就是$$\left[\begin{array}{ccccc}3& 2& 1& |& 3\\ 0& -\frac{1}{2}& \frac{7}{2}& |& 6\\ 0& 0& \frac{27}{7}& |& \frac{54}{7}\end{array}\right]$$
可以看到，除去最后一行常数，系数构成了一个类似三角形的矩阵，叫做阶梯矩阵（虽然我个人还是习惯叫三角矩阵）
所以我们的目标就是将普通的矩阵化为三角矩阵，这里就是纯粹的模拟了。具体不懂的看注释吧……

#define Abs(x) ( x > 0 ? x : -x )
const double EPS = 1e-6;//这个东西是拿来预防double误差的,一般来说,一个数的绝对值<EPS则这个数为0
int n, m;//n个方程,m个未知数
double a[401][402];//矩阵
double ans[401];//储存答案,即所有x的解
bool fr[401];//自由元，是自由元为0，否则为1
register int now, k;
for(k = 1, now = 1; k <= n && now <= m; k++, now++) {//要遍历每一个方程,每一个未知数
    int j = k;
    //因为前面的都已经化成了三角矩阵,所以要在还未进行消元的部分,
    //即还没有化为三角矩阵的部分进行消元
    for(register int i = k + 1; i <= n; i++)
        if(Abs(a[j][now]) < Abs(a[i][now])) j = i;
        //找出当前这一个未知数中最大的系数，因为这样可以尽可能地避免误差
    if(Abs(a[j][now]) < EPS) { k--; continue; }
    //如果当前这一个未知数所有的系数都为0,因为最大值为0,则遍历下一个未知数
    if(j != k) swap(a[j], a[k]);//不用担心，是不会CE的
    //直接把消元的当前这一行放到前面的一行，这样就可以构造一个三角形的矩阵
    //处理当前这一行的系数,以构造三角形矩阵
    for(register int i = 1; i <= n; i++) {//对每一个方程消元
        if(Abs(a[i][now]) < EPS || i == k) continue;//若当前未知数的系数为0,当然就跳过啦
        double D = a[i][now] / a[k][now];
        //例题是浮点型的,所以用除法消元,如果是整形,就应该用gcd(a[i][now], a[k][now])
        for(register int x = now; x <= m + 1; x++) a[i][x] -= D * a[k][x];//消这个未知数
    }
}

可以明显看出，高斯消元的复杂度是$O(n^3$+巨大常数$)$的。（啊这真是一个快速的算法）

特殊情况处理

无解

例如有如下方程：
$$\{\begin{array}{l}x+y=1\\ x+y=2\end{array}$$
显然冲突，无解。这种情况该如何处理呢？
很显然，如果有这样的一行$\left[\begin{array}{ccccccc}0& 0& 0& \cdots & 0& |& b\end{array}\right](b̸=0)$，那这个方程显然无解，于是可以用这样的一行来判断一下是否有解：

for(reg int i = k; i <= n; i++)//为什么从k开始?因为前k开始系数都是被消掉的全是0了,这样才满足条件
    if(Abs(a[i][m + 1]) > EPS) return !printf("No solution");
//因为有消元了过后,前面全是0,而后面常数项非0的,这种情况当然无解

自由元

自由元是什么？
如下：$$\{\begin{array}{l}4x_1+2x_2+x_3=3①\\ 2x_1+x_2+x_3=6②\\ x_3=9③\end{array}$$
可以解出$x_3=9$，但是$x_1,x_2$都没有确定的解，这个就叫做自由元
因为消了$k$个元，所以如果未知数的个数$m>k$，那就可能存在自由元。

if(k <= m) {//存在自由元
    for(reg int i = k - 1; i >= 1; i--) {
        int Unknown = 0, idx = 0;//统计该方程中未出现的未知数的个数和那个数的下标
        for(reg int j = 1; j <= m; j++)//找每一个系数
            if(!fr[j] && Abs(a[i][j]) > EPS) Unknown++, idx = j;
        //非0,且不是已经有解了的未知数
        if(Unknown > 1) continue;//此方程暂时无解
        double tmp = a[i][m + 1];
        for(reg int j = 1; j <= m; j++)
            if(j != idx && Abs(a[i][j]) > EPS) tmp -= a[i][j] * ans[j];//模拟回代不解释
        //倒着推回来，进行代入消元
        ans[idx] = tmp / a[i][idx];//得到这个未知数的解,还要记住可能存在的系数
        fr[idx] = 1;//这个未知数肯定不是自由元啦~
    }
}

普通情况

for(reg int i = m; i >= 1; i--) {//m个元，m个含系数的未知式
    double tmp = a[i][m + 1];
    for(reg int j = i + 1; j <= m; j++)
        tmp -= a[i][j] * ans[j];//代入……
    ans[i] = tmp / a[i][i];//除以系数,得解
    fr[i] = 1;//不是自由元
}//普通的回带不解释

代码总结

#define reg register
#define Abs(a) ( a > 0 ? a : -a )
const double EPS = 1e-6;
int n, m;
double a[401][402];
double ans[401];
bool fr[401];//自由元
inline bool Gauss() {//返回1表示有解，返回0表示无解
    register int now, k;
    for(k = 1, now = 1; k <= n && now <= m; k++, now++) {//要遍历每一个方程,每一个未知数
        int j = k;
        //因为前面的都已经化成了三角矩阵,所以要在还未进行消元的部分,
        //即还没有化为三角矩阵的部分进行消元
        for(register int i = k + 1; i <= n; i++)
            if(Abs(a[j][now]) < Abs(a[i][now])) j = i;
            //找出当前这一个未知数中最大的系数，因为这样可以尽可能地避免误差
        if(Abs(a[j][now]) < EPS) { k--; continue; }
        //如果当前这一个未知数所有的系数都为0,因为最大值为0,则遍历下一个未知数
        if(j != k) swap(a[j], a[k]);//不用担心，是不会CE的
        //直接把消元的当前这一行放到前面的一行，这样就可以构造一个三角形的矩阵
        //处理当前这一行的系数,以构造三角形矩阵
        for(register int i = 1; i <= n; i++) {//对每一个方程消元
            if(Abs(a[i][now]) < EPS || i == k) continue;//若当前未知数的系数为0,当然就跳过啦
            double D = a[i][now] / a[k][now];
            //例题是浮点型的,所以用除法消元,如果是整形,就应该用gcd(a[i][now], a[k][now])
            for(register int x = now; x <= m + 1; x++) a[i][x] -= D * a[k][x];//消这个未知数
        }
    }
    for(reg int i = k; i <= n; i++)//为什么从k开始?因为前k开始系数都是被消掉的全是0了,这样才满足条件
        if(Abs(a[i][m + 1]) > EPS) return 0;
    //因为有消元了过后,前面全是0,而后面常数项非0的,这种情况当然无解
    if(k <= m) {//存在自由元
        for(reg int i = k - 1; i >= 1; i--) {
            int Unknown = 0, idx = 0;//统计该方程中未出现的未知数的个数和那个数的下标
            for(reg int j = 1; j <= m; j++)//找每一个系数
                if(!fr[j] && Abs(a[i][j]) > EPS) Unknown++, idx = j;
            //非0,且不是已经有解了的未知数
            if(Unknown > 1) continue;//此方程暂时无解
            double tmp = a[i][m + 1];
            for(reg int j = 1; j <= m; j++)
                if(j != idx && Abs(a[i][j]) > EPS) tmp -= a[i][j] * ans[j];//模拟回代不解释
            //倒着推回来，进行代入消元
            ans[idx] = tmp / a[i][idx];//得到这个未知数的解,还要记住可能存在的系数
            fr[idx] = 1;//这个未知数肯定不是自由元啦~
        }
    return 1;
    }
    for(reg int i = m; i >= 1; i--) {//m个元，m个含系数的未知式
        double tmp = a[i][m + 1];
        for(reg int j = i + 1; j <= m; j++)
            tmp -= a[i][j] * ans[j];//代入……
        ans[i] = tmp / a[i][i];//除以系数,得解
        fr[i] = 1;//不是自由元
    }//普通的回带不解释
    return 1;
}
int main() {
    scanf("%d%d", &n, &m);
    for(reg int i = 1; i <= n; i++)
        for(reg int j = 1; j <= m + 1; j++)
            scanf("%lf", &a[i][j]);
    bool sol = Gauss();
    if(!sol) {
        printf("No solution");
        return 0;
    }
    for(reg int i = 1; i <= m; i++) {
        if(!fr[i]) printf("X[%d] not determined\n", i);
        else {
            printf("X[%d] = ", i);
            if(Abs(ans[i]) < EPS)
                ans[i] = 0;
            printf("%.4lf\n", ans[i]);
        }
    }
}

附上：上方代码没有头文件的原因

#include<map>
#include<set>
#include<list>
#include<queue>
#include<deque>
#include<stack>
#include<ctime>
#include<cmath>
#include<vector>
#include<bitset>
#include<cctype>
#include<string>
#include<sstream>
#include<fstream>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<iomanip>
#include<iostream>
#include<algorithm>
#include<tr1/unordered_map>
using namespace std;
using namespace tr1;