首页 > 其他分享 >高斯消元

高斯消元

时间：2023-01-29 16:33:29浏览次数：49  
标签：return int il cs ri 高斯消 define

简述

高斯消元法（Gauss-Jordan elimination）是求解线性方程组的经典算法，它在当代数学中有着重要的地位和价值，是线性代数课程教学的重要组成部分。
高斯消元法除了用于线性方程组求解外，还可以用于行列式计算、求矩阵的逆，以及其他计算机和工程方面。

主要步骤

高斯消元法在将增广矩阵化为最简形后对于自由未知量的赋值，需要掌握线性相关知识，且赋值存在人工经验的因素，使得在学习过程中有一定的困难，将高斯消元法划分为五步骤，从而提出五步骤法，内容如下：

  1. 增广矩阵行初等行变换为行最简形；
  2. 还原线性方程组；
  3. 求解第一个变量；
  4. 补充自由未知量；
  5. 列表示方程组通解。

解线性方程组

il int Gauss(int n,int m){
    ri int x=1,y=1;
    for(x=1;x<=n&&y<m;++x,++y){
        ri int mx=x;
        for(ri int i=x+1;i<=n;++i){
            if(abs(a[mx][y])<abs(a[i][y])) mx=i;
        }
        if(abs(a[mx][y])<eps){--x;continue;}
        if(mx!=x) swap(a[mx],a[x]);
        for(ri int i=y+1;i<=m;++i) a[x][i]/=a[x][y];
        a[x][y]=1;
        for(ri int i=x+1;i<=n;++i){
            for(ri int j=y+1;j<=m;++j){
                a[i][j]-=a[i][y]*a[x][j];
            }
            a[i][y]=0;
        } 
    }
    if(x<=n){
        for(ri int i=x;i<=n;++i){
            if(abs(a[i][m])>eps) return -1;
        }
        return n-x+1;		
    }
    for(ri int i=n;i;--i){
        xi[i]=a[i][m];
        for(ri int j=i-1;j;--j){
            a[j][m]-=xi[i]*a[j][i];
        }
    }
    return 0;
}
AC·code 
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
#define abs(a) (a<0?-a:a)
using namespace std;
il int rd(){
    ri int x=0,f=1;ri char c=getchar();
    while(c<'0'||c>'9') f=(c=='-')?-f:f,c=getchar();
    while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
    return x*f;
}
cs int N=105;cs double eps=1e-7;
int n,fl;double a[N][N],xi[N];
il int Gauss(int n,int m){
    ri int x=1,y=1;
    for(x=1;x<=n&&y<m;++x,++y){
        ri int mx=x;
        for(ri int i=x+1;i<=n;++i){
            if(abs(a[mx][y])<abs(a[i][y])) mx=i;
        }
        if(abs(a[mx][y])<eps){--x;continue;}
        if(mx!=x) swap(a[mx],a[x]);
        for(ri int i=y+1;i<=m;++i) a[x][i]/=a[x][y];
        a[x][y]=1;
        for(ri int i=x+1;i<=n;++i){
            for(ri int j=y+1;j<=m;++j){
                a[i][j]-=a[i][y]*a[x][j];
            }
            a[i][y]=0;
        } 
    }
    if(x<=n){
        for(ri int i=x;i<=n;++i){
            if(abs(a[i][m])>eps) return -1;
        }
        return n-x+1;		
    }
    for(ri int i=n;i;--i){
        xi[i]=a[i][m];
        for(ri int j=i-1;j;--j){
            a[j][m]-=xi[i]*a[j][i];
        }
    }
    return 0;
}
int main(){
    n=rd();
    for(ri int i=1;i<=n;++i){
        for(ri int j=1;j<=n+1;++j){
            a[i][j]=1.0*rd();
        }
    }
    fl=Gauss(n,n+1);
    if(fl) printf("No Solution");
    else{
        for(ri int i=1;i<=n;++i){
            printf("%.2lf\n",xi[i]);
        }
    }
    return 0;
}

解异或方程组

注意到异或方程组的增广矩阵是 \(01\) 矩阵（矩阵中仅含有 \(0\) 与 \(1\)），所以我们可以使用 \(C++\) 中的 \(std::bitset\) 进行优化，将时间复杂度降为 \(O(\dfrac{n^2m}{\omega})\)，其中 \(n\) 为元的个数，\(m\) 为方程条数，\(\omega\) 一般为 \(32\)（与机器有关）。

bitset<1010> a[2010];
il int Gauss(int n,int m){
	int res=0,mx=n;
	for(ri int i=1;i<=n;++i){
        res=i;
        for(ri int j=i;j<=m;++j){
            if(a[j][i]){
                res=j;
                mx=max(res,mx);
                break;
            }
        }
        if(!a[res][i]) return -1;//多解或无解 
        if(res!=i) swap(a[res],a[i]);
        for(ri int j=1;j<=m;++j){
            if(i!=j&&a[j][i]) a[j]^=a[i];
        }
    }
    return mx;//在第mx行可以确定唯一解
}
AC·code 
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
#define F(s) freopen(#s".in","r",stdin),freopen(#s".out","w",stdout);
using namespace std;
int n,m;
bool a[2005][1005];// 不用bitset优化也能过
char c[1005];
il int Gauss(int n,int m){
    int res=0,mx=n;
    for(ri int i=1;i<=n;++i){
        res=i;
        for(ri int j=i;j<=m;++j){
            if(a[j][i]){
                res=j;
                mx=max(res,mx);
                break;
            }
        }
        if(!a[res][i]) return -1;
        if(res!=i) swap(a[res],a[i]);
        for(ri int j=1;j<=m;++j){
            if(j!=i&&a[j][i]){
                for(ri int k=i;k<=n+1;++k){
                    a[j][k]^=a[i][k];
                }
            }
        }
    }
    return mx;
}
int main(){
    ios::sync_with_stdio(0);
    cin>>n>>m;
    for(ri int i=1;i<=m;++i){
        cin>>c;
        for(ri int j=0;j<strlen(c);++j){
            a[i][j+1]=(c[j]=='1');
        }
        cin>>a[i][n+1];
    }
    int f=Gauss(n,m);
    if(~f){
        cout<<f<<'\n';
        for(ri int i=1;i<=n;++i){
            if(a[i][n+1]) cout<<"?y7M#\n";
            else cout<<"Earth\n";
        }
    }
    else cout<<"Cannot Determine\n";
    return 0;
}

行列式求值

il int getdet(int n,int p,int a[][N]){
    ri int as=1;
    for(ri int i=1;i<=n;++i){
        for(ri int j=i+1;j<=n;++j){
            while(a[i][i]){//辗转相减
                ri int l=a[j][i]/a[i][i];
                a[j][i]%=a[i][i];
                swap(a[j][i],a[i][i]);
                for(ri int k=i+1;k<=n;++k){
                    a[j][k]=(a[j][k]-l*a[i][k]%p+p)%p;//!!!!!
                    swap(a[j][k],a[i][k]);
                }
                as*=-1;
            }
            swap(a[i],a[j]),as*=-1;
        }
        as=(as*a[i][i]%p+p)%p;
        if(!as) return 0;
    } 
    return as;
}
AC·code 
#include<bits/stdc++.h>
#define cs const
#define il inline
#define ri register
#define int long long
using namespace std;
il int rd(){
    ri int x=0,f=1;ri char c=getchar();
    while(c<'0'||c>'9') f=(c=='-')?-f:f,c=getchar();
    while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
    return x*f;
}
il void wt(int x){
    if(x<0) x=-x,putchar('-');if(x>9) wt(x/10);
    return putchar(x%10+48),void();
}
cs int N=605;
int n,p,a[N][N];
il int getdet(int n,int p,int a[][N]){
    ri int as=1;
    for(ri int i=1;i<=n;++i){
        for(ri int j=i+1;j<=n;++j){
            while(a[i][i]){
                ri int l=a[j][i]/a[i][i];
                a[j][i]%=a[i][i];
                swap(a[j][i],a[i][i]);
                for(ri int k=i+1;k<=n;++k){
                    a[j][k]=(a[j][k]-l*a[i][k]%p+p)%p;//!!!!!
                    swap(a[j][k],a[i][k]);
                }
                as*=-1;
            }
            swap(a[i],a[j]),as*=-1;
        }
        as=(as*a[i][i]%p+p)%p;
        if(!as) return 0;
    } 
    return as;
}
signed main(){
    n=rd(),p=rd();
    for(ri int i=1;i<=n;++i){
        for(ri int j=1;j<=n;++j){
            a[i][j]=rd();
        }
    }
    wt(getdet(n,p,a));
    return 0;
}

矩阵求逆

求\(A\)的逆矩阵，把\(A\)和单位矩阵\(I\)放在一个矩阵里
对\(A\)进行加减消元使\(A\)化成单位矩阵
此时原来单位矩阵转化成逆矩阵
如果无法将\(A\)化为单位矩阵，则不存在逆矩阵

il bool Gauss(cs int n,cs int m){
    for(ri int i=1;i<=n;++i){
        ri int l=i,inv=1;
        for(ri int j=i;j<=n;++j){
            if(a[j][i]){l=j;break;}
        }
        if(!a[l][i]) return 0;
        if(l!=i) swap(a[i],a[l]);
        inv=getinv(a[i][i]);
        for(ri int j=i;j<=m;++j){
            a[i][j]=(1ll*a[i][j]*inv)%mod; 
        } 
        for(ri int j=1;j<=n;++j){
            if(j==i) continue;
            for(ri int k=i+1;k<=m;++k){
                a[j][k]=(a[j][k]-1ll*a[j][i]*a[i][k]%mod+mod)%mod;
            }
            a[j][i]=0;
        }
    }
    return 1;
}
AC·code 
#include<bits/stdc++.h>
#define il inline
#define ri register
#define cs const
using namespace std;
il int rd(){
    ri int x=0,f=1;ri char c=getchar();
    while(c<'0'||c>'9') f=(c=='-')?-1:1,c=getchar();
    while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
    return x*f;
}
il void wt(int x){
    if(x<0) x=-x,putchar('-');
    if(x>=10) wt(x/10);
    return putchar(x%10+48),void();
}
cs int mod=1e9+7,N=405;
int n,a[N][N*2];
il void exgcd(int a,int b,long long &x,long long &y){
    if(!b) return x=1,y=0,void();
    return exgcd(b,a%b,y,x),y-=x*(a/b),void();
}
il int getinv(int a){
    ri long long x,y;exgcd(a,mod,x,y);
    return (x%mod+mod)%mod;
}
il bool Gauss(cs int n,cs int m){
    for(ri int i=1;i<=n;++i){
        ri int l=i,inv=1;
        for(ri int j=i;j<=n;++j){
            if(a[j][i]){l=j;break;}
        }
        if(!a[l][i]) return 0;
        if(l!=i) swap(a[i],a[l]);
        inv=getinv(a[i][i]);
        for(ri int j=i;j<=m;++j){
            a[i][j]=(1ll*a[i][j]*inv)%mod; 
        } 
        for(ri int j=1;j<=n;++j){
            if(j==i) continue;
            for(ri int k=i+1;k<=m;++k){
                a[j][k]=(a[j][k]-1ll*a[j][i]*a[i][k]%mod+mod)%mod;
            }
            a[j][i]=0;
        }
    }
    return 1;
} 
int main(){
    n=rd();
    for(ri int i=1;i<=n;++i){
        for(ri int j=1;j<=n;++j){
            a[i][j]=rd()%mod;
        }
        a[i][n+i]=1;
    }
    ri bool fl=Gauss(n,n*2);
    if(fl){
        for(ri int i=1;i<=n;++i){
            for(ri int j=n+1;j<=n*2;++j){
                wt(a[i][j]),putchar(' ');
            }
            putchar('\n');
        }
    }
    else printf("No Solution");
    return 0;
}

标签：return,int,il,cs,ri,高斯消,define
From： https://www.cnblogs.com/windymoon/p/17073064.html

相关文章

  • C++ 数学与算法系列之高斯消元法求解线性方程组
    1.前言什么是消元法？消元法是指将多个方程式组成的方程组中的若干个变量通过有限次地变换，消去方程式中的变量，通过简化方程式，从而获取结果的一种解题方法。消元法主要有代......
  • 数值分析·学习 | 解线性方程组的直接方法（高斯消去法以及LU求解）matlab实现
    ​ 目录一、前言：二、算法描述：三、实现代码：1、高斯消去法：2、高斯消去法-列主元消去法：3、LU分解：4、求逆矩阵：四、总结：一、前言：个人学习内容分享二、算法描述：1......
  • 高斯消元+组合数+卡特兰数
    高斯消元+组合数+卡特兰数高斯消元\(O(n^3)\)的线性时间内求解n元线性方程组\[\\\begin{cases}\a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1\\\a_{21}x_1+a_{22}x_2+.......
  • 算法学习笔记（39）——高斯消元
    高斯消元高斯消元高斯消元解线性方程组高斯消元解异或线性方程组高斯消元解线性方程组通过初等行变换把增广矩阵化为阶梯型矩阵，并回代得到方程的解适用于求解......
  • 高斯消元与线性基
    高斯消元与线性基Guass—约旦消元消元算法简介：这是求解线性方程组（也就是M个N元一次方程组）的方法思想：我们可以把方程组看作一个系数矩阵例如：\[\left\{\begin{aligned......
  • 高斯消元&高斯约旦消元
    高斯消元就是上三角，然后再回代。高斯约旦消元就是消的时候直接变成对角线了，你选取当前主元，然后把其他的都消去这个元。一般来说就写后者。注意二者都要特判自由元，但常数......
  • 高斯消元
    0x00背景高斯消元是求线性方程组的标准方法,原理和代码都不难0x01基本操作一个线性方程组有\(m\)个一次方程,\(n\)个变量,把所有系数都写成一个\(m\)行\(n\)列的矩阵,......
  • 第45届国际大学生程序设计竞赛 亚洲区域赛（济南）A // 高斯消元
    题目来源：第45届国际大学生程序设计竞赛（ICPC）亚洲区域赛（济南）A题目链接：A-MatrixEquation题意定义01矩阵的两种运算：对于大小为\(N\timesN\)的01矩阵\(X\)、\(Y\)，记......
  • 高斯消元
    高斯消元高斯消元适用于求解线性方程。线性方程形如\[\begin{cases}a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n=b_1\\a_{2,1}x_1+a_{2,2}x_2+...+a_{2,n}x_n=b_2\\...\\a_{......
  • 【XSY4180】串串游走（AC自动机，期望DP，高斯消元）
    假瑞出搬的神仙题。原题CFgym103119B。先把\(T\)去重。考虑先用\(O(nm\logk)\)建出所有串的AC自动机。注意建AC自动机的时候，为了保证空间，假设当前点\(u\)没......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99