高斯消元简单笔记

最新推荐文章于 2022-11-21 20:46:06 发布

佯设

最新推荐文章于 2022-11-21 20:46:06 发布

阅读量216

点赞数

分类专栏： acm算法学习笔记文章标签：算法 c++

本文链接：https://blog.csdn.net/m0_63737271/article/details/125639483

版权

acm算法学习笔记专栏收录该内容

36 篇文章 0 订阅

订阅专栏

萌新一枚，借鉴了许多大佬而总结下的笔记，可能有瑕疵或者错误，望指正

1将方程组记录成矩阵形式，通过矩阵变换，得到一个上三角的矩阵，这样就方便求解

2本文记录了异或，整数，浮点数，和模意义的方程组

3主要步骤为：

（1）寻找当前列最大的一行，换到当前行上去，再将当前行以下同列的元素全都变为0

（2）判断消元后的结果，如果方程组右边的数为0，而左边的数为0，说明无解。

如果有效的方程的数量和小于未知数的数量，说明有多解，有自由变量。

如果有效方程的数量等于未知数的数量，方程有唯一解

代码实现

1数据存储

double a[maxn][maxn],ans[maxn];//增广矩阵的数据，以及每个方程的解 
const double eps=1e-6;//浮点数误差大小 
int Free[maxn];//记录自由变量

一，整数方程组

1求最大公约数方程

int gcd(int a,int b)//求最大公因数 
{
	while(b)
	{
		int t=b;
		b=a%b;
		a=t;
	}
	return a;
}

2求最小公倍数方程

int lcm(int a,int b)//求最小公倍数 
{
	return a/gcd(a,b)*b;
}

3寻找当前列最大元素的行数，并交换，消去当前行以下当前列的元素

int row,col,maxrow;//当前行数，当前列数，当前列中元素最大的一行 
	col=0;
	for(row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row;
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		if(abs(a[i][col]>a[maxrow][col]))
		maxrow=i;
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		for(int i=row;i<=var;i++)
		swap(a[row][i],a[maxrow][i]);
		if(a[row][col]==0)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变
		{
			row--;
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		{
			if(a[i][col])
			{
				//先寻找当前列第i行和第row行元素的最小公倍数，并分别找到需要乘以的倍数 
				int temp1=lcm(abs(a[i][col]),abs(a[row][col]))/abs(a[i][col]);
				int temp2=lcm(abs(a[i][col]),abs(a[row][col]))/abs(a[row][col]);
				//第i行整行乘以temp2再减去temp1倍的row行，就可以消去当前行元素 
				if(a[i][col]*a[row][col]<0)
				temp2=-temp2;
				for(int j=col;j<=var;j++)
				a[i][j]=a[i][j]*temp1-a[row][j]*temp2;
			}
		}
	}

4判断有无解，并求解

	for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解
	if(a[i][col])
	return -1;
	if(var>row)//如果有效行数小于总行数，说明有自由变量 
	return var-row;
	//如果方程数等于未知数的个数，有唯一解 
	for(int i=var-1;i>=0;i--)//方程常规求解 
	{
		int temp=a[i][var];
		for(int j=i+1;j<var;j++)
		{
			if(a[i][j])
			temp-=ans[j]*a[i][j];
		}
		ans[i]=temp/a[i][i];
	}

5完整代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=100;
int a[maxn][maxn],ans[maxn];//记录矩阵，记录未知数的解 
int gcd(int a,int b)//求最大公因数 
{
	while(b)
	{
		int t=b;
		b=a%b;
		a=t;
	}
	return a;
}
int lcm(int a,int b)//求最小公倍数 
{
	return a/gcd(a,b)*b;
}
int Gauss(int equ,int var)//返回自由变量的数量   equ为行数,var为列数 
{
	for(int i=0;i<=var;i++)//初始化 
	ans[i]=0;
	int row,col,maxrow;//当前行数，当前列数，当前列中元素最大的一行 
	col=0;
	for(row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row;
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		if(abs(a[i][col]>a[maxrow][col]))
		maxrow=i;
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		for(int i=row;i<=var;i++)
		swap(a[row][i],a[maxrow][i]);
		if(a[row][col]==0)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变
		{
			row--;
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		{
			if(a[i][col])
			{
				//先寻找当前列第i行和第row行元素的最小公倍数，并分别找到需要乘以的倍数 
				int temp1=lcm(abs(a[i][col]),abs(a[row][col]))/abs(a[i][col]);
				int temp2=lcm(abs(a[i][col]),abs(a[row][col]))/abs(a[row][col]);
				//第i行整行乘以temp2再减去temp1倍的row行，就可以消去当前行元素 
				if(a[i][col]*a[row][col]<0)
				temp2=-temp2;
				for(int j=col;j<=var;j++)
				a[i][j]=a[i][j]*temp1-a[row][j]*temp2;
			}
		}
	}
	for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解
	if(a[i][col])
	return -1;
	if(var>row)//如果有效行数小于总行数，说明有自由变量 
	return var-row;
	//如果方程数等于未知数的个数，有唯一解 
	for(int i=var-1;i>=0;i--)//方程常规求解 
	{
		int temp=a[i][var];
		for(int j=i+1;j<var;j++)
		{
			if(a[i][j])
			temp-=ans[j]*a[i][j];
		}
		ans[i]=temp/a[i][i];
	}
	return 0;
 } 
int main()
{
	return 0;
}

二，浮点数方程

1寻找当前列最大元素的行数，并交换，消去当前行以下当前列的元素

for(int row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row; 
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		{
			if(fabs(a[i][col])>fabs(a[maxrow][col]))
			maxrow=i;
		}
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		{
			for(int i=row;i<=var;i++)
			swap(a[maxrow][i],a[row][i]);
		}
		if(fabs(a[row][col])<eps)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变 
		{
			row--;
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		{
			if(fabs(a[i][col])>eps)
			{
				double temp=a[i][col]/a[row][col];//第i行要减去几倍的row行，使得col列变为0。 
				for(int j=col;j<=var;j++)//当前列后面的元素也随之发生变化 
				a[i][j]-=a[i][row]*temp;
				a[i][col]=0;
			}
		}
	}

2判断有无解，并求解

for(int i=0;i<=var;i++)//初始化 
		{
		ans[i]=0;
		Free[i]=1;//开始假定所由元素都是自由变量 
		}
		for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解 
		if(fabs(a[i][col])>eps)
		return -1;
		if(row<var)//如果有效行数小于总行数，说明有自由变量 
		{
			double temp;
			for(int i=row-1;i>=0;i--)//从最下往上依次求解 
			{
				int cnt=0,idx;
				for(int j=0;j<var;j++)
				if(a[i][j]&&Free[j])//如果xj有系数并且是自由变量（暂时是） 
				{
					cnt++;//记录自由变量的值 
					idx=j;//计算当前xj的值 
				}
				if(cnt>1)//如果只有xj（暂时是）是自由变量，那么也就是说其他x都已经确定了 
				continue;
				temp=a[i][var];//取右边的数 
				for(int j=0;j<var;j++)
				if(a[i][j]&&j!=idx)
				temp-=ans[j]*a[i][j];//常规的方程求解 
				ans[idx]=temp/=a[i][idx];//xj的值确定了，所以xj不再是自由变量了 
				Free[idx]=0;
			}
			return var-row;//返回自由变量的值 
		}
		//如果方程数等于未知数的个数，有唯一解 
		for(int i=var-1;i>=0;i--)//常规的方程求解 
		{
			double temp=a[i][var];
			for(int j=i+1;j<var;j++)
			if(a[i][j])
			temp-=ans[j]*a[i][j];
			ans[i]=temp/=a[i][i];
		}

3完整代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=200;
double a[maxn][maxn],ans[maxn];//增广矩阵的数据，以及每个方程的解 
const double eps=1e-6;//浮点数误差大小 
int Free[maxn];//记录自由变量 
int Gauss(int equ,int var)//返回自由变量的数量   equ为行数,var为列数 
{
	
	int row,col,maxrow;//当前行数，当前列数，当前列中元素最大的一行 
	col=0;
	for(int row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row; 
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		{
			if(fabs(a[i][col])>fabs(a[maxrow][col]))
			maxrow=i;
		}
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		{
			for(int i=row;i<=var;i++)
			swap(a[maxrow][i],a[row][i]);
		}
		if(fabs(a[row][col])<eps)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变 
		{
			row--;
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		{
			if(fabs(a[i][col])>eps)
			{
				double temp=a[i][col]/a[row][col];//第i行要减去几倍的row行，使得col列变为0。 
				for(int j=col;j<=var;j++)//当前列后面的元素也随之发生变化 
				a[i][j]-=a[i][row]*temp;
				a[i][col]=0;
			}
		}
	}	
		for(int i=0;i<=var;i++)//初始化 
		{
		ans[i]=0;
		Free[i]=1;//开始假定所由元素都是自由变量 
		}
		for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解 
		if(fabs(a[i][col])>eps)
		return -1;
		if(row<var)//如果有效行数小于总行数，说明有自由变量 
		{
			double temp;
			for(int i=row-1;i>=0;i--)//从最下往上依次求解 
			{
				int cnt=0,idx;
				for(int j=0;j<var;j++)
				if(a[i][j]&&Free[j])//如果xj有系数并且是自由变量（暂时是） 
				{
					cnt++;//记录自由变量的值 
					idx=j;//计算当前xj的值 
				}
				if(cnt>1)//如果只有xj（暂时是）是自由变量，那么也就是说其他x都已经确定了 
				continue;
				temp=a[i][var];//取右边的数 
				for(int j=0;j<var;j++)
				if(a[i][j]&&j!=idx)
				temp-=ans[j]*a[i][j];//常规的方程求解 
				ans[idx]=temp/=a[i][idx];//xj的值确定了，所以xj不再是自由变量了 
				Free[idx]=0;
			}
			return var-row;//返回自由变量的值 
		}
		//如果方程数等于未知数的个数，有唯一解 
		for(int i=var-1;i>=0;i--)//常规的方程求解 
		{
			double temp=a[i][var];
			for(int j=i+1;j<var;j++)
			if(a[i][j])
			temp-=ans[j]*a[i][j];
			ans[i]=temp/=a[i][i];
		}
	return 0;
}
int main()
{
	
	return 0;
}

三，异或方程组

1数据存储

bitset<maxn>a[maxn];//记录矩阵数据，可以将整个a[i]与a[x]异或 ，比较方便 
int ans[maxn];//记录未知数的解 
int Free[maxn];//记录自由变量

2寻找当前列最大元素的行数，并交换，消去当前行以下当前列的元素

int row,col,maxrow;//当前行，当前列，当前列元素最大的一行 
	col=0;
	for(row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row;
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		if(abs(a[i][col])>abs(a[maxrow][col]))
		maxrow=i;
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		swap(maxrow,row);
		if(a[row][col]==0)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变 
		{
			row--;
			Free[++cnt]=col;//当前列的未知数就是自由变量 
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		if(a[i][col])
		a[i]^=a[maxrow];	
	}

3判断有无解，并求解

for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解
if(a[i][col])
return -1;
if(var>row)//如果有效行数小于总行数，说明有自由变量 
	return var-row;
	//如果方程数等于未知数的个数，有唯一解 
	for(int i=var-1;i>=0;i--)//方程常规求解 
	{
		int temp=a[i][var];
		for(int j=i+1;j<var;j++)
		{
			if(a[i][j])
			temp-=ans[j]*a[i][j];
		}
		ans[i]=temp/a[i][i];
	}

4完整代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=100;
int a[maxn][maxn],ans[maxn];//记录矩阵，记录未知数的解 
int gcd(int a,int b)//求最大公因数 
{
	while(b)
	{
		int t=b;
		b=a%b;
		a=t;
	}
	return a;
}
int lcm(int a,int b)//求最小公倍数 
{
	return a/gcd(a,b)*b;
}
int Gauss(int equ,int var)//返回自由变量的数量   equ为行数,var为列数 
{
	for(int i=0;i<=var;i++)//初始化 
	ans[i]=0;
	int row,col,maxrow;//当前行数，当前列数，当前列中元素最大的一行 
	col=0;
	for(row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row;
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		if(abs(a[i][col]>a[maxrow][col]))
		maxrow=i;
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		for(int i=row;i<=var;i++)
		swap(a[row][i],a[maxrow][i]);
		if(a[row][col]==0)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变
		{
			row--;
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		{
			if(a[i][col])
			{
				//先寻找当前列第i行和第row行元素的最小公倍数，并分别找到需要乘以的倍数 
				int temp1=lcm(abs(a[i][col]),abs(a[row][col]))/abs(a[i][col]);
				int temp2=lcm(abs(a[i][col]),abs(a[row][col]))/abs(a[row][col]);
				//第i行整行乘以temp2再减去temp1倍的row行，就可以消去当前行元素 
				if(a[i][col]*a[row][col]<0)
				temp2=-temp2;
				for(int j=col;j<=var;j++)
				a[i][j]=a[i][j]*temp1-a[row][j]*temp2;
			}
		}
	}
	for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解
	if(a[i][col])
	return -1;
	if(var>row)//如果有效行数小于总行数，说明有自由变量 
	return var-row;
	//如果方程数等于未知数的个数，有唯一解 
	for(int i=var-1;i>=0;i--)//方程常规求解 
	{
		int temp=a[i][var];
		for(int j=i+1;j<var;j++)
		{
			if(a[i][j])
			temp-=ans[j]*a[i][j];
		}
		ans[i]=temp/a[i][i];
	}
	return 0;
 } 
int main()
{
	return 0;
}

四，模意义方程组
1快速幂，用于求逆元

int qpow(int base,int n)//快速幂，这里为了求逆元 
{
	int ret=1;
	while(n)
	{
		if(n&1)
		ret=ret*base%mod;
		base=base*base%mod;
		n>>=1;
	}
	return ret;
}

2寻找当前列最大元素的行数，并交换，消去当前行以下当前列的元素

int maxrow,row,col;//当前行数，当前列数，当前列中元素最大的一行 
	col=0;
	for(row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row;
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		if(abs(a[i][col])>abs(a[maxrow][col]))
		maxrow=i;
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		{
			for(int i=col;i<=var;i++)
			swap(a[row][i],a[maxrow][i]);
		}
		if(a[row][col]==0)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变 
		{
			row--;
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		{
			if(a[i][col])
			{
				int temp=a[i][col]*qpow(a[row][col],mod-2)%mod;//倍数等于a[i][col]/a[row][col],相当于乘以a[row][col]的逆元 
				for(int j=col;j<=var;j++)
				a[i][j]=(a[i][j]-a[row][j]*temp%mod+mod)%mod;
				//第i行减去temp倍的第row行 
			}
		}
	}

3判断有无解，并求解

for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解 
	if(a[i][col])
	return -1;
	if(var>row)//如果有效行数小于总行数，说明有自由变量 
	return var-row;
		//如果方程数等于未知数的个数，有唯一解 
	for(int i=var-1;i>=0;i--)//常规的方程求解,这里除法变为乘以逆元 
	{
		int temp=a[i][var];
		for(int j=i+1;j<var;j++)
		{
			if(a[i][j])
			{
				temp-=a[i][j]*ans[j];
				temp=(temp%mod+mod)%mod;
			}
		}
		ans[i]=temp*qpow(a[i][i],mod-2)%mod;
	}

4完整代码

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int maxn=2e2+10;
int a[maxn][maxn];
int ans[maxn];//增广矩阵的数据，以及每个方程的解 
int qpow(int base,int n)//快速幂，这里为了求逆元 
{
	int ret=1;
	while(n)
	{
		if(n&1)
		ret=ret*base%mod;
		base=base*base%mod;
		n>>=1;
	}
	return ret;
}
int Gauss(int equ,int var)//返回自由变量的数量   equ为行数,var为列数 
{
	int maxrow,row,col;//当前行数，当前列数，当前列中元素最大的一行 
	col=0;
	for(row=0;row<equ&&col<var;row++,col++)
	{
		maxrow=row;
		for(int i=row+1;i<equ;i++)//寻找当前列中元素绝对值最大的一行 
		if(abs(a[i][col])>abs(a[maxrow][col]))
		maxrow=i;
		if(maxrow!=row)//如果当前行不是最大的一行，就交换 
		{
			for(int i=col;i<=var;i++)
			swap(a[row][i],a[maxrow][i]);
		}
		if(a[row][col]==0)//如果当前一整列全为0，那么让row减一，这样进入下一循环就能保持不变 
		{
			row--;
			continue;
		}
		for(int i=row+1;i<equ;i++)//将当前列的当前行一下的元素全都消为0 
		{
			if(a[i][col])
			{
				int temp=a[i][col]*qpow(a[row][col],mod-2)%mod;//倍数等于a[i][col]/a[row][col],相当于乘以a[row][col]的逆元 
				for(int j=col;j<=var;j++)
				a[i][j]=(a[i][j]-a[row][j]*temp%mod+mod)%mod;
				//第i行减去temp倍的第row行 
			}
		}
	}
	for(int i=row;i<equ;i++)//如果进行高斯消元后，增广矩阵右边的b不为0，而左边为0，说明无解 
	if(a[i][col])
	return -1;
	if(var>row)//如果有效行数小于总行数，说明有自由变量 
	return var-row;
		//如果方程数等于未知数的个数，有唯一解 
	for(int i=var-1;i>=0;i--)//常规的方程求解,这里除法变为乘以逆元 
	{
		int temp=a[i][var];
		for(int j=i+1;j<var;j++)
		{
			if(a[i][j])
			{
				temp-=a[i][j]*ans[j];
				temp=(temp%mod+mod)%mod;
			}
		}
		ans[i]=temp*qpow(a[i][i],mod-2)%mod;
	}
	return 0;
}
int main()
{
	return 0;
}