学习笔记——高斯消元

简介

高斯消元是一个用于求方程组的解的算法。在线性代数中非常重要。一般而言，其复杂度为 $O(n^3)$ ，可以承受 $200$ 及以下的数据范围(当然有的题目时限是 $10000ms$ 什么的，特殊情况特殊对待)。

算法理解

考虑一个二元一次方程组怎么解。

{\begin{cases} a_{1} x_{1} + a_{2} x_{2} = a_{3} \\ b_{1} x_{1} + b_{2} x_{2} = b_{3} \end{cases}

$\begin{cases}a_1x_1+a_2x_2=a_3\\b_1x_1+b_2x_2=b_3\end{cases}$

首先将第一个方程化为 $x_1+\frac{a_2}{a_1}x_2=\frac{a_3}{a_1}$ ，然后两边再同乘 $b_1$ ，这样在两个方程中， $x_1$ 的系数就都化为 $b_1$ ，然后两式相减，就可以消去 $x_1$ ，然后就可以解出 $x_2$ ，再回代，求解 $x_1$ 。

那么多元方程组也是一样的思路。

{\begin{cases} a_{1, 1} x_{1} + a_{1, 2} x_{2} + \dots = a_{1, n + 1} \\ a_{2, 1} x_{1} + a_{2, 2} x_{2} + \dots = a_{2, n + 1} \\ ⋮ \\ a_{n, 1} x_{1} + a_{n, 2} x_{2} + \dots = a_{n, n + 1} \end{cases}

$\begin{cases}a_{1,1}x_1+a_{1,2}x_2+\dots=a_{1,n+1}\\a_{2,1}x_1+a_{2,2}x_2+\dots=a_{2,n+1}\\\vdots\\a_{n,1}x_1+a_{n,2}x_2+\dots=a_{n,n+1}\end{cases}$

首先把第一个方程化为 $x_1+\frac{a_{1,2}}{a_{1,1}}x_2+\dots=\frac{a_{1,n+1}}{a_{1,1}}$ ，然后对于之后的每一个方程，都把第一个方程乘上当前这个方程的 $x_1$ 的系数，然后相减，就可以得到 $n-1$ 个不含 $x_1$ 的方程，然后再对这剩下的 $n-1$ 个方程，重复上述操作即可。最后先求解出 $x_n$ ，然后不停地代到上一个方程里去解，最终得到整个方程组的解。

那么在代码实现的时候，是先把每个未知数的系数拎出来，形成一个矩阵(一般是 $n*(n+1)$ )。
为了代码美观方便起见，下标从0开始。

[\begin{matrix} a_{0, 0} & a_{0, 1} & \dots & a_{0, n} \\ a_{1, 0} & a_{1, 1} & \dots & a_{1, n} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ a_{n - 1, 0} & a_{n - 1, 1} & \dots & a_{n - 1, n} \end{matrix}]

$\begin{bmatrix}a_{0,0}&a_{0,1}&\cdots&a_{0,n}\\a_{1,0}&a_{1,1}&\cdots&a_{1,n}\\\vdots&\vdots&\vdots&\vdots\\a_{n-1,0}&a_{n-1,1}&\cdots&a_{n-1,n}\end{bmatrix}$

美其名曰：增广矩阵。注意，最后一列是等号右边的常数项。

算法实现

枚举当前要消的元，在增广矩阵中就是枚举第几列 $col$ 。

$step1.1$
选取矩阵中的一行 $row$ ，将含这个元的系数化为 $1$ 。即

[\begin{matrix} 1 & \frac{a_{0, 1}}{a_{0, 0}} & \dots & \frac{a_{0, n}}{a_{0, 0}} \\ a_{1, 0} & a_{1, 1} & \dots & a_{1, n} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ a_{n - 1, 0} & a_{n - 1, 1} & \dots & a_{n - 1, n} \end{matrix}]

$\begin{bmatrix}1&\frac{a_{0,1}}{a_{0,0}}&\cdots&\frac{a_{0,n}}{a_{0,0}}\\a_{1,0}&a_{1,1}&\cdots&a_{1,n}\\\vdots&\vdots&\vdots&\vdots\\a_{n-1,0}&a_{n-1,1}&\cdots&a_{n-1,n}\end{bmatrix}$

特别地，为了减小精度损失和代码方便，一般是找系数最大的一行。而如果系数最大的系数也是 $0$ ，那么说明当前这个元的系数不用消，则直接进入下一个消元，并对此元进行标记(至于标记什么后面会提到)。

$step1.2$
对于后面的每一行 $i$ ，都将第 $row$ 的 $col$ (刚刚已经化为 $1$ )乘上第 $i$ 行第 $col$ 的系数。然后两行相减，即可以将 $i$ 行的第 $col$ 这一元消去。

枚举完后，增广矩阵应该是变成了上三角矩阵，即

[\begin{matrix} \dots & \dots & \dots & \dots \\ 0 & \dots & \dots & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & \dots & x_{n - 1} \end{matrix}]

$\begin{bmatrix}\cdots&\cdots&\cdots&\cdots\\0&\cdots&\cdots&\cdots\\\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&x_{n-1}\end{bmatrix}$

有重要性质：一个位置为 $0$ ，则其左下所有位置都为 $0$ 。
然后求解

$step2.1$
先判断是否有解。
那么如果无解，说明有某一次消元时得到了 $0=x(x\not=0)$ 的情况。所以在最后的增广矩阵中找未知数系数全为 $0$ 但是等号右边不为 $0$ 的行，即 $0\;0\;0\cdots0\;x(x\not=0)$ 的行，则无解。
在代码中，运用上面增广矩阵的性质，可以从最后的 $row$ 开始直接遍历最后一个常数项看是否有不等于 $0$ 的，则其左边的所有系数一定都为 $0$ ，此时可以判断无解。
$step2.2$
如果多解，那么说明消元的时候出现了 $0=0$ 的情况。所以在增广矩阵中找全 $0$ 行，如果有，则说明有多个解。
注意：对于无解或多解情况，增广矩阵还是遵照类三角矩阵，因此，如果按照上面所述的方式进行消元，那么最后 $row$ 比方程个数小，就会多解，而之前标记的就是可以取任意数的元，我们称之为自由元，一般而言，最后自由元个数 $freenum=equ-row$ ，其中 $equ$ 为方程个数。
$step2.3$
如果有唯一解，那么就回代求解。此时增广矩阵一定是严格的上三角矩阵，先求出最后一个元，然后将最后一个元代入倒数第二个消元后的方程中，一定可以化为 $ax=b$ ，就可以依次求出每一个元，从而解得方程。

几种常见方程类型

算法放在 $Gauss$ 函数中，你需要先将方程组构造在 $a[][]$ 中，然后调用：

int res=Gauss(equ,var);
if(res==-1) 无解
//else if(res==-2) 浮点解，整数线性方程组中用
else if(res>0) 多解，此时res为自由元个数
else 有唯一解，输出x[]即可

浮点线性方程组#

此种题目要特别注意精度问题，一般 $eps$ 取 $10^{-12}$ 。
这种方程组一般放在电学题目中，此时运用基尔霍夫电流定律，对于每个节点列出方程。

double a[MAXN][MAXN],x[MAXN];//a:增广矩阵	x:解
bool frex[MAXN];//标记是否为自由元
int Gauss(int equ,int var){//equ:方程个数	var:未知数个数
	for(int i=0;i<=var;i++) x[i]=0,frex[i]=1;//初始化
	//step1	消元
	int col=0,row;//枚举第row行，第col列，即第row个方程，消去第col个元
	for(row=0;row<equ&&col<var;row++,col++){
		//step1.1
		int mx_r=row;//找出系数最大的一行，减小精度，便于判断
		for(int i=row+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[mx_r][col]))
				mx_r=i;
		if(mx_r!=row)
			for(int j=col;j<=var;j++)
				swap(a[row][j],a[mx_r][j]);//把最大的那一行与当前处理的交换
		if(fabs(a[row][col])<eps){row--;continue;}//如果最大的一行都是0，那么该元已经是0，跳过
		//step1.2
		for(int i=row+1;i<equ;i++)
			if(fabs(a[i][col])>eps){
				double tmp=a[i][col]/a[row][col];//把当前方程的col元系数化为1
				for(int j=col;j<=var;j++)
					a[i][j]-=a[row][j]*tmp;//乘上col元的系数，消元
				a[i][col]=0;//避免精度影响，当前方程该元赋为0
			}
	}
	//step2
	for(int i=row;i<equ;i++)
		if(fabs(a[i][col])>eps)
			return -1;//step2.1
	if(var>row) return var-row;//step2.2
	//step2.3
	for(int i=var-1;i>=0;i--){
		double tmp=a[i][var];
		for(int j=i+1;j<var;j++)
			tmp-=x[j]*a[i][j];
		x[i]=tmp/a[i][i];//回代求解
	}return 0;
}

例题
模板题
 Electric resistance
计算电压

异或方程组#

此种题目是形如

{\begin{cases} a_{0, 0} x_{0} \otimes a_{0, 1} x_{0} \otimes \dots = a_{0, n} \\ a_{1, 0} x_{1} \otimes a_{1, 1} x_{1} \otimes \dots = a_{1, n} \\ ⋮ \\ a_{n - 1, 0} x_{n - 1} \otimes a_{n - 1, 1} x_{n - 1} \otimes \dots = a_{n - 1, n} \end{cases}

$\begin{cases}a_{0,0}x_0\otimes a_{0,1}x_0\otimes\cdots=a_{0,n}\\a_{1,0}x_1\otimes a_{1,1}x_1\otimes\cdots=a_{1,n}\\\vdots\\a_{n-1,0}x_{n-1}\otimes a_{n-1,1}x_{n-1}\otimes\cdots=a_{n-1,n}\end{cases}$

众所周知，异或相当于模二，所以消元的时候就是系数之间异或就可以了。
一般是关联开关问题，对于每个开关，比较前后状态和相邻的开关，可以得出方程。

int a[MAXN][MAXN],x[MAXN],frX[MAXN];//记录自由元
int Gauss(int equ,int var){
	for(int i=0;i<=var;i++)
		x[i]=frX[i]=0;
	//step1
	int col=0,fnum=0,row;
	for(row=0;row<equ&&col<var;row++,col++){
		//step1.1
		int mx_r=row;
		for(int i=row+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[mx_r][col]))
				mx_r=i;
		if(mx_r!=row)
			for(int j=row;j<=var;j++)
				swap(a[row][j],a[mx_r][j]);
		if(a[row][col]==0){
			frX[fnum++]=col;//如果当前的元不用消，说明是自由元
			row--;continue;
		}
		//step1.2
		for(int i=row+1;i<equ;i++)
			if(a[i][col]!=0)
				for(int j=col;j<=var;j++)
					a[i][j]^=a[row][j];
	}
	//step2
	for(int i=row;i<equ;i++)
		if(a[i][col]!=0)
			return -1;//step2.1
	int temp=var-row;
	if(row<var) return temp;//step2.2
	//step2.3
	for(int i=var-1;i>=0;i--){
		x[i]=a[i][var];
		for(int j=i+1;j<var;j++)
			x[i]^=(a[i][j]&&x[j]);
	}return 0;
}

例题
EXTENDED LIGHTS OUT
开关问题

特别地，对于异或方程组，还有对自由元的处理，可能出现问你有多少自由元，或者解为 $1$ 的最少个数等等。其中问 $1$ 的最少个数，只能 $dfs$ 枚举算一下，目前没有更好的办法。
$dfs$ 时，注意上面 $Gauss$ 的板子里进行了行交换，所以算的时候不能直接用 $a[i][j]$ ，要记录一下当前的行，如果是自由元，则行不变，如果不是，则行上升。由于交换后，有效的 $a[][]$ 都在增广矩阵的上部，所以行直接从总方程数减去自由元开始，这样可以使每次都能找到相对应的 $a[][]$ 进行计算。

int ans,var;
void dfs(int row,int i,int sum){
	if(sum>=ans) return;//剪枝，很快啊
	if(i<0){ans=sum;return;}
	if(!frX[i]){
		x[i]=a[row][var];
		for(int j=i+1;j<var;j++)
			x[i]^=(a[row][j]&&x[j]);
		dfs(row-1,i-1,sum+x[i]);
	}else{
		x[i]=0;dfs(row,i-1,sum);
		x[i]=1;dfs(row,i-1,sum+1);
	}
}
int main()
{
	...dfs(equ-freenum-1,var-1,0);
	printf("%d",ans);//即最少的1的个数
}

例题
Lights G
The Water Bowls

整数线性方程组#

此种题目不允许浮点数解。一般会判断或者输出无整数解等。~~一般都是裸题。~~

int GCD(int x,int y){return x%y==0?y:GCD(y,x%y);}
int LCM(int x,int y){return x/GCD(x,y)*y;}
int a[MAXN][MAXN],x[MAXN];
bool frex[MAXN];
int Gauss(int equ,int var){
	for(int i=0;i<=var;i++) x[i]=0,frex[i]=1;
	//step1
	int col=0,row;
	for(row=0;row<equ&&col<var;row++,col++){
		//step1.1
		int mx_r=row;
		for(int i=row+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[mx_r][col]))
				mx_r=i;
		if(mx_r!=row)
			for(int j=col;j<=var;j++)
				swap(a[row][j],a[mx_r][j]);
		if(a[row][col]==0){row--;continue;}
		//step1.2
		for(int i=row+1;i<equ;i++)
			if(a[i][col]!=0){
				int lcm=LCM(abs(a[i][col]),abs(a[row][col]));
				int ta=lcm/abs(a[i][col]),tb=lcm/abs(a[row][col]);//为了避免精度损失，要用lcm
				if(a[i][col]*a[row][col]<0) tb=-tb;
				for(int j=col;j<=var;j++)
					a[i][j]=a[i][j]*ta-a[row][j]*tb;
			}
	}
	//step2
	for(int i=row;i<equ;i++)
		if(a[i][col]!=0)
			return -1;//step2.1
	if(var>row) return var-row;//step2.2
	//step2.3
	for(int i=var-1;i>=0;i--){
		int tmp=a[i][var];
		for(int j=i+1;j<var;j++)
			tmp-=x[j]*a[i][j];
		if(tmp%a[i][i]!=0) return -2;//如果无法整除，则无整数解
		x[i]=tmp/a[i][i];
	}return 0;
}

例题
模板题
暂时没有找到别的。

模线性方程组#

与整数方程组类似，但是当无整数解的时候可以加模数进行更正。比较麻烦，所以放在最后。

int a[MAXN][MAXN],x[MAXN],frX[MAXN];
int GCD(int x,int y){return x%y==0?y:GCD(y,x%y);}
int LCM(int x,int y){return x/GCD(x,y)*y;}
int Gauss(int equ,int var){
	for(int i=0;i<=var;i++)
		x[i]=0,frX[i]=1;
	//step1
	int col=0,row;
	for(row=0;row<equ&&col<var;row++,col++){
		//step1.1
		int mx_r=row;
		for(int i=row+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[mx_r][col]))
				mx_r=i;
		if(mx_r!=row)
			for(int j=row;j<=var;j++)
				swap(a[row][j],a[mx_r][j]);
		if(a[row][col]==0){row--;continue;}
		//step1.2
		for(int i=row+1;i<equ;i++)
			if(a[i][col]!=0){
				int lcm=LCM(abs(a[i][col]),abs(a[row][col]));
				int ta=lcm/abs(a[i][col]);
				int tb=lcm/abs(a[row][col]);
				if(a[i][col]*a[row][col]<0) tb=-tb;
				for(int j=col;j<=var;j++)
					a[i][j]=((a[i][j]*ta-a[row][j]*tb)%MOD+MOD)%MOD;
			}
	}
	//step2
	for(int i=row;i<equ;i++)
		if(a[i][col]!=0)
			return -1;//step2.1
	if(row<var) return var-row;//step2.2
	//step2.3
	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]!=0)
				temp-=a[i][j]*x[j];
			temp=(temp%MOD+MOD)%MOD;
		}
		while(temp%a[i][i]!=0)
			temp+=MOD;//如果不能整除，则不断加模数来更正
		x[i]=(temp/a[i][i])%MOD;
	}return 0;
}

可能与星期月份等相结合。模数一般较小。

但是也不排除个别毒瘤题目，模数很大，导致加了很久还无法整除，会使得时间复杂度大大上升。因此，可以用扩欧来优化这一过程。

int exGCD(int a,int b,int &x,int &y){
	if(b==0){x=1;y=0;return a;}
	int c=exGCD(b,a%b,y,x);
	y-=a/b*x;return c;
}//扩欧
int solve(int a,int b,int c){
	int x,y;
	int r=exGCD(a,b,x,y);
	x=x*c/r;b/=r;
	x=(x%b+b)%b;
	return x;
}//解ax+by=c的不定方程，返回x的最小正整数解
const int MAXN=75;
int a[MAXN][MAXN],x[MAXN],frX[MAXN],MOD;
int GCD(int x,int y){return x%y==0?y:GCD(y,x%y);}
int LCM(int x,int y){return x/GCD(x,y)*y;}
int Gauss(int equ,int var){
	for(int i=0;i<=var;i++)
		x[i]=0,frX[i]=1;
	int col=0,row;
	for(row=0;row<equ&&col<var;row++,col++){
		int mx_r=row;
		for(int i=row+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[mx_r][col]))
				mx_r=i;
		if(mx_r!=row)
			for(int j=row;j<=var;j++)
				swap(a[row][j],a[mx_r][j]);
		if(a[row][col]==0){row--;continue;}
		for(int i=row+1;i<equ;i++)
			if(a[i][col]!=0){
				int lcm=LCM(abs(a[i][col]),abs(a[row][col]));
				int ta=lcm/abs(a[i][col]);
				int tb=lcm/abs(a[row][col]);
				if(a[i][col]*a[row][col]<0) tb=-tb;
				for(int j=col;j<=var;j++)
					a[i][j]=((a[i][j]*ta-a[row][j]*tb)%MOD+MOD)%MOD;
			}
	}
	for(int i=row;i<equ;i++)
		if(a[i][col]!=0)
			return -1;
	if(row<var) return 1;
	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]!=0)
				temp-=a[i][j]*x[j];
			temp=(temp%MOD+MOD)%MOD;
		}
		temp+=MOD*solve(MOD,a[i][i],-temp);//直接扩欧得出
		//其他部分和暴力枚举是一样的
		x[i]=(temp/a[i][i])%MOD;
	}return 0;
}

例题
Widget Factory
SETI

The End

posted @ 2021-04-12 17:32 ZCETHAN 阅读(195) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· TypeScript + Deepseek 打造卜卦网站：技术与玄学的结合
· Manus的开源复刻OpenManus初探
· AI 智能体引爆开源社区「GitHub 热点速览」
· 三行代码完成国际化适配，妙~啊~
· .NET Core 中如何实现缓存的预热？

公告

昵称： ZCETHAN
园龄： 4年2个月
粉丝： 11
关注： 10

+加关注

2025年3月

日

一

二

三

四

五

六

Loading

ZCETHAN

学习笔记——高斯消元

简介

算法理解

算法实现

几种常见方程类型

浮点线性方程组#

异或方程组#

整数线性方程组#

模线性方程组#

The End

公告

搜索

常用链接

最新随笔

我的标签

积分与排名

随笔分类 (139)

随笔档案 (107)

阅读排行榜

评论排行榜

推荐排行榜

最新评论