高斯消元

时间:2019-11-08
本文章向大家介绍高斯消元，主要包括高斯消元使用实例、应用技巧、基本知识点总结和需要注意事项，具有一定的参考价值，需要的朋友可以参考一下。

ll逆元版

//by Saber Alter Official
//Erishikigal & Ishtar
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <queue>

typedef long long ll;
typedef unsigned long long ull;

const int N = 515, mod = 1000000007;

ll n, m;//n个方程， m个元 

//用来求逆元的快速幂，不多嗦了 
inline ll QuickPow(ll a, ll b) {
    ll Cur = 1;
    while(b) {
        if(b & 1) {
            Cur *= a;
            Cur %= mod;
        }
        a *= a;
        a %= mod;
        b >>= 1;
    }
    return Cur;
}
 
struct Matrix {
    ll A[N][N], x[N];//A存的是矩阵，也就是方程的系数  x是用来存解哒 
    ll Valid;//储存的是消元后，存储的有效方程的个数，也就是有主元的可解方程
    Matrix () {
        memset(A, 0, sizeof A);
        memset(x, 0, sizeof x);
    }
    
    inline void FormulaSwap(int i, int j) {//交换……这个不用多说了，1到m的系数全换一遍 
        for(int k = 1;k <= m;k++) {
            ll _temp = A[i][k];
            A[i][k] = A[j][k];
            A[j][k] = _temp;
        }
    }
    
    inline void Elimate(int j, int i) {//消元啊！！！ A[j][j]是我们找到的主元的系数位置，A[i][j]是我们要消掉那个方程的元的系数 
        //逆元来消元啊！不是很在意精度或者模数非素就用double直接除啊！ 
        ll Coe = (A[i][j] * QuickPow(A[j][j], mod - 2)) % mod;//我们要找到A[i][j] = k * A[j][j]！也就是k = A[i][j] * A[j][j]!
        Coe = (Coe + mod) % mod;//防止倍数变成负的 ！ 
        for(int k = j;k <= m;k++) {//从A[i][j]向后开始消元啊！ 
            A[i][k] = (A[i][k] - Coe * A[j][k]) % mod; //都减去对应的主元方程位置处的系数!也就是Coe(k) * A[j][j->m]！注意要全都减掉！我们是在两个方程之间加减！ 
        } 
    }
    
    inline void Gauss() {//高斯消元 
        for(int j = 1;j <= m - 1;j++) {//从第一行开始枚举主元！1到m-1！因为只有m-1个未知数！第m个是右侧常数哇 
            if(!A[j][j]) {//假如我们枚举的这个主元系数是0！那么就要开始向下面找到这个元的第一个系数不是0的方程！交换这两个！然后再进行下去！ 
                for(int i = j + 1;i <= n;i++) {//向下找啊！找到第一个该元系数不是零的位置啊！ 
                    if(A[i][j] != 0) {//找到了啊！ 
                        FormulaSwap(i, j);//交换啊 ！ 
                        break;//跳出来啊！ 交换完了！到下面去消元了！ 
                    }
                }
            }
            if(!A[j][j]) continue;//全走完了，还找不到啊！换吧！这个元不用消啊！ 
            for(int i = j + 1;i <= n;i++) {//从这一行向下消元啊！！！ 
                Elimate(j, i); //消元啊！注意顺序啊！j是我们枚举的主元的位置啊！我们要从主元向下消啊！ 
            }
        }
    }
    
    
    inline bool Nosolution() {//无解的判断啊！ 
        bool NoSol = true;//先假设无解啊！ 
        for(int i = n;i >= 1;i--) {//从最后一个方程开始看！ 
            for(int j = 1;j <= m - 1;j++) {//从左到右扫一遍方程左边的系数！ 
                if(A[i][j]) {//左边剩下系数！有解啊！ 
                    NoSol = false;//记下来！不是无解！ 
                    Valid = i;//记下来！这是消元之后第一次出现有效方程的地方！ 
                    break;//跳出来就好了！ 
                }
            }
            if(!NoSol) return false;//你看我们上面扫了一遍最后一样找到了可解方程啊！ 
            if(A[i][m]) return true;//我们扫完最后一行，没有进入到上面的判断！那就是左边莫得系数！但是！右边不为0!这是无解条件啊！ 
        }
        return false;//这句没啥大用……为了防止Wall报而已 
    }
    
    inline bool ManySolutions() {//求完上面的情况，发现有效方程的数目小于我们有的未知数个数！有自由元啊！ 
        return Valid < m - 1;//记住！我们未知数的个数只有m-1个！第m个是方程右边的常数！ 
    }
    
    inline void Solve() {//开始解方程 
        for(int i = m - 1;i >= 1;i--) {//从最后一个未知数开始 
            x[i] = A[i][m] % mod;//先把它设成右边的常数 
            for(int j = i + 1;j <= m - 1;j++) {//假如说右边有一个元，通过上三角阵形式可以知道，那么它一定已经得出来了解 
                x[i] = (x[i] - A[i][j] * x[j]) % mod;//移项求解就可以了 
                x[i] = (x[i] + mod) % mod;//换成正的 
            }
            x[i] = (x[i] * QuickPow(A[i][i], mod - 2)) % mod;//假如说没有右边的未知数，那就是一个Ax=B的方程，直接逆元解就好了 
        }
        
        for(int i = 1;i <= m - 1;i++) {
            printf("%lld\n", x[i]);
        }
    }
    
}Ans;

int main() {
    scanf("%lld %lld", &n, &m);
    for(int i = 1;i <= n;i++) {
        for(int j = 1;j <= m;j++) {
            scanf("%lld", &Ans.A[i][j]);
        }
    }
    Ans.Gauss();
    if(Ans.Nosolution()) {
        printf("No solution.");
    }
    else if(Ans.ManySolutions()) {
        printf("Many Solutions.");
    }
    else {
        printf("Unique solution.\n");
        Ans.Solve();
    }
    
    return 0;
}

double除法版

//by Saber Alter Official
//Erishikigal & Ishtar
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <queue>

//#define Rin

typedef long long ll;
typedef unsigned long long ull;

const int N = 115;

struct Matrix {
    int n, m;//n formula,(m - 1 = n)x, 1 constant
    double A[N][N];
    double x[N];
    int valid;
    
    inline void FormulaSwap(int j, int i) {
        for(int k = 1;k <= m;k++) {
            double _temp = A[j][k];
            A[j][k] = A[i][k];
            A[i][k] = _temp;
        }
    }
    
    inline void Elimate(int j, int i) {
        double coe = A[i][j] / A[j][j];
        for(int k = 1;k <= m;k++) {
            A[i][k] = A[i][k] - coe * A[j][k];
        }
    }
    
    inline void Gauss() {
        for(int j = 1;j <= n;j++) {//枚举方程编号 
            if(!A[j][j]) {//方程主元系数为0 
                for(int i = j + 1;i <= n;i++) {
                    if(A[i][j]) {
                        FormulaSwap(j, i);
                        break;
                    }
                }
            }
            if(!A[j][j]) continue;
            for(int i = j + 1;i <= n;i++) {
                Elimate(j, i);
            }
        }
    }
    
    inline bool NoSolution() {
        bool NoSolu = true;
        for(int i = n;i >= 1;i--) {
            for(int j = 1;j <= n;j++) {
                if(A[i][j]) {
                    NoSolu = false;
                    valid = i;
                    break;
                }
            }
            if(!NoSolu) return false;
            if(A[i][m]) return true;
        }
        return false;
    }
    
    inline bool ManySolutions() {
        return valid < n;
    }
    
    inline void Solve() {
        for(int i = n;i >= 1;i--) {//枚举方程 
            x[i] = A[i][m];//设为右边的常数 
            for(int j = i + 1;j <= n;j++) {//枚举右边的元 
                x[i] = x[i] - A[i][j] * x[j];
            }
            x[i] = x[i] / A[i][i];
        }
        
        for(int i = 1;i <= n;i++) {
            printf("%.2lf\n", x[i]);
        }
    }
    
}Ans;

int main() {
    scanf("%d", &Ans.n);
    Ans.m = Ans.n + 1;
    for(int i = 1;i <= Ans.n;i++) {
        for(int j = 1;j <= Ans.m;j++) {
            scanf("%lf", &Ans.A[i][j]);
        }
    }
    Ans.Gauss();
    if(Ans.NoSolution() || Ans.ManySolutions()) {
        #ifdef Rin
        printf("%d %d\n", Ans.NoSolution(), Ans.ManySolutions());
        printf("%d\n", Ans.valid);
        #endif
        printf("No Solution\n");
        return 0;
    }
    Ans.Solve();
    
    return 0;
}

原文地址：https://www.cnblogs.com/Fructose-Ryllis/p/11818274.html

随机文章
本站知识点必读
最近更新