高斯消元法

· · 题解

题目~题目~~：P3389 【模板】高斯消元法

经过我与lhf小朋友将近一个小时的研究，终于成功解读实现了高斯消元法

写一篇博客来纪念一下下

//别粘，粘了肯定错，但是可以看一下思路（其实本篇分析不太多，主要是注释比较可爱）好吧本质是因为我太菜了然后和机房小伙伴琢磨了好久不会实现还请大佬们移步不要重度吐槽

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <bits/stdc++.h>
#include <cmath>
using namespace std;
const int N = 110;
const double eps = 1e-6;
int n;
double a[N][N];
int gauss()
{
    int nowy, nowx;
    for (nowy = 0, nowx = 0; nowy < n; nowy ++ ){
        int pivot = nowx;
        for (int i = nowx; i < n; i ++ )
            if (fabs(a[i][nowy]) > fabs(a[pivot][nowy]))
                pivot = i;
        if (fabs(a[pivot][nowy]) < eps) continue;
        for (int i = nowy; i < n + 1; i ++ ) 
            swap(a[pivot][i], a[nowx][i]);
        for (int i = n; i >= nowy; i -- ) 
            a[nowx][i] /= a[nowx][nowy];
        for (int i = nowx + 1; i < n; i ++ )
            if (fabs(a[i][nowy]) > eps)
                for (int j = n; j >= nowy; j -- )
                    a[i][j] -= a[nowx][j] * a[i][nowy];
        nowx ++ ;
    }
    if (nowx < n){
        for (int i = nowx; i < n; i ++ )
            if (fabs(a[i][n]) > eps)
                return 2;
        return 1;
    }
    for (int i = n - 1; i >= 0; i -- )
        for (int j = i + 1; j < n; j ++ )
            a[i][n] -= a[j][n] * a[i][j];

    return 0;
}

int main()
{
    cin >> n;
    for (int i = 0; i < n; i ++ )
        for (int j = 0; j < n + 1; j ++ )
            cin >> a[i][j];
    int t = gauss();
    if (t == 0)
        for (int i = 0; i < n; i ++ ) printf("%.2lf\n", a[i][n]);   
    else if (t == 1) puts("Infinite group solutions");
    else puts("No Solution\n");
    return 0;
}

正文部分开始...

上三角形矩阵：

  1. 每行第一个元编号严格大于前一行
  2. 可以有若干行左边不含未知数，右边为常数的“方程”

想了一下，感觉这个题目还不太完善，不容易表达，因此又拿来了另一个题目

高斯消元2

（不保证未知数和方程的数量相当）

然后，为保证代码可读性，我打算标一堆注释...大业开始...

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>

#define lld long long
#define db double

using namespace std;
const int mod = 1e9 + 7;
const int N = 509;

lld qpow(lld a, lld b) {//快速幂无需解释
    lld Ans = 1;
    while (b) {
        if (b & 1) {
            Ans = (Ans * a) % mod;
        }
        a = (a * a) % mod;
        b >>= 1;
    }
    return Ans;
}//快速幂无需解释结束

int n, m;//n个方程，m-1个未知数
lld a[N][N], x[N];//a[][]来存放系数， x来存放未知数解
bool Unique = 1, No_Ans = 0;//立两个flag来处理不同情况（看英文就懂了）

void Swap(int i, int j) {
//实现矩阵的交换行
    for (int k = 1; k <= m; k++) swap(a[i][k], a[j][k]);
}

void gauss() {//高斯定理正文内容
//gauss预处理部分
    for (int j = 1; j < m; j++) {//从第一个未知数开始走，一直到最后一个未知数
        if (a[j][j] == 0) {//如果说，三角形出现了缺位（当然是在保证三角形消元的前提下）
            bool flag = 1;//先立一个事关唯一解的flag（现在假设解不唯一）
            for (int i = j + 1; i <= n; i++) {//从这一行的后面开始走，一直到最后一行，寻找系数绝对值更大的一行，并且进行交换
                if (a[i][j]) {//如果找到了下方有未知数的位置不在它后面
                    Swap(i, j);//就进行这两行的交换（以保证三角形成立）
                    flag = 0;//并且设定唯一解
                    break;//跳出寻找更大系数的循环
                }
            }
            if (flag)//如果说flag没有被推翻，那么解不唯一
                Unique = 0;
        }
//gauss预处理部分结束                      
        //从下一行开始（下一个方程），一直到最后一个方程
        for (int i = j + 1; i <= n; i++) {
            //不是因为矩阵运算不能够做除法，而是因为这个题需要取模，因此我们寻找逆元（顺便搞了负号）
            lld k = -(a[i][j] * qpow(a[j][j], mod - 2)) % mod;  //找逆元
            //找完逆元之后，我们就可以进行减法运算了

            for (int r = 1; r <= m; r++) //以后每行的每一个未知数前的系数都要减去当前（第i行）的对应未知数前的系数（就是消元的具体操作）
                a[i][r] = (a[i][r] + a[j][r] * k) % mod;
        }
    }

    bool zero = 1;//先立一个flag代表左边未知数系数为0（代表“空”，没有待求解未知数）
    //这样的方程组是否有解只与最后左边不含未知数，右边为常数的“方程”有关(方程有解当且仅当这些“方程”等式的右边为0，前面的方程只要满足第一个条件，无论如何都不会导致方程无解
    for (int i = n; i >= 1; i--) {//从最后一行开始走到第一行
        for (int j = 1; j <= m - 1; j++)//从第一个未知数开始走到最后一个未知数
            if (a[i][j]) {//如果说这一行的某个未知数的系数不是零（就是说我们找到了这一行里待解的未知数）
                zero = 0;//更改一下flag
                break;//跳出循环
            }
        if (!zero)//如果说，查到了左边带系数的未知数，就可以跳出循环了。（仔细看大括号）
            break;
        if (a[i][m] != 0) {//如果右边常数不是0，那就没有解了，更改一下flag，跳出就可以了（毕竟能走到这一步，那必须得是前方的系数有不为零的）
            No_Ans = 1;
            break;
        }
    }
}

void solve() {//求具体解
    for (int i = m - 1; i >= 1; i--) {//从最后一个未知数开始走，到第一个未知数结束

        x[i] = a[i][m];//第i个未知数的解x[i]先用这一行等号右边的常数赋值
        for (int j = i + 1; j <= m - 1; j++)//从下一个未知数（的位置）开始，到最后一个未知数（的位置）结束
            x[i] = (x[i] - a[i][j] * x[j]) % mod;//对本行未知数求解

 //把以上过程总结一下的话就是为了求解本行的未知数，我们用等号右边的常数减去本方程其余所有的量,获得本行未知数前的最后系数（此时的顺序尤其重要）
        //不要忘了我们是在解未知数，所以还要把x[i](也就是未知数前的系数)除走。选择求逆元，并在最后记得取模。
        x[i] = (x[i] * qpow(a[i][i], mod - 2)) % mod;
        x[i] = (x[i] + mod) % mod;
        }
}

void print() {//输出x1, x2, x3, x4,...,x_(m-1)
    for (int i = 1; i < m; i++) printf("%lld\n", x[i]);
}

int main() {
    scanf("%d%d", &n, &m);//有n个方程，m-1个未知数
    for (int i = 1; i <= n; i++)//输入第i行
        for (int j = 1; j <= m; j++)//输入每行的系数及最后的常数项
            scanf("%lld", &a[i][j]);
    gauss();//进行高斯定理（被搞死
    if (No_Ans)//noans初值设置是0，变为1说明没有答案
        puts("No solution.");
    else if (!Unique)//unique初值设定为1，变为0说明有多解
        puts("Many Solutions.");
    else {//剩下的就是正解部分了
        puts("Unique solution.");
        solve();//解方程（没错只有到了正解部分我们才是需要求解方程具体解的）
        print();//输出
    }
    return 0;
}