高斯消元
简述
高斯消元法（Gauss-Jordan elimination）是求解线性方程组的经典算法，它在当代数学中有着重要的地位和价值，是线性代数课程教学的重要组成部分。
高斯消元法除了用于线性方程组求解外，还可以用于行列式计算、求矩阵的逆，以及其他计算机和工程方面。
主要步骤
高斯消元法在将增广矩阵化为最简形后对于自由未知量的赋值，需要掌握线性相关知识，且赋值存在人工经验的因素，使得在学习过程中有一定的困难，将高斯消元法划分为五步骤，从而提出五步骤法，内容如下：
- 增广矩阵行初等行变换为行最简形；
- 还原线性方程组；
- 求解第一个变量；
- 补充自由未知量；
- 列表示方程组通解。
解线性方程组
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;
}
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行可以确定唯一解
}
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;
}
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;
}
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;
}
I went to the woods because I wanted to live deliberately, I wanted to live deep and suck out all the marrow of life, and not when I had come to die, discover that I had not live.
