学习笔记——高斯消元

简介

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

算法理解

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

\[\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}\\\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{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{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{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\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