学习笔记——高斯消元

简介

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

算法理解

考虑一个二元一次方程组怎么解。
{ a 1 x 1 + a 2 x 2 = a 3 b 1 x 1 + b 2 x 2 = b 3 \begin{cases}a_1x_1+a_2x_2=a_3\\b_1x_1+b_2x_2=b_3\end{cases} {a1​x1​+a2​x2​=a3​b1​x1​+b2​x2​=b3​​
首先将第一个方程化为 x 1 + a 2 a 1 x 2 = a 3 a 1 x_1+\frac{a_2}{a_1}x_2=\frac{a_3}{a_1} x1​+a1​a2​​x2​=a1​a3​​，然后两边再同乘 b 1 b_1 b1​，这样在两个方程中， x 1 x_1 x1​的系数就都化为 b 1 b_1 b1​，然后两式相减，就可以消去 x 1 x_1 x1​，然后就可以解出 x 2 x_2 x2​，再回代，求解 x 1 x_1 x1​。

那么多元方程组也是一样的思路。
{ a 1 , 1 x 1 + a 1 , 2 x 2 + ⋯ = a 1 , n + 1 a 2 , 1 x 1 + a 2 , 2 x 2 + ⋯ = a 2 , n + 1 ⋮ a n , 1 x 1 + a n , 2 x 2 + ⋯ = a n , n + 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} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​a1,1​x1​+a1,2​x2​+⋯=a1,n+1​a2,1​x1​+a2,2​x2​+⋯=a2,n+1​⋮an,1​x1​+an,2​x2​+⋯=an,n+1​​
首先把第一个方程化为 x 1 + a 1 , 2 a 1 , 1 x 2 + ⋯ = a 1 , n + 1 a 1 , 1 x_1+\frac{a_{1,2}}{a_{1,1}}x_2+\dots=\frac{a_{1,n+1}}{a_{1,1}} x1​+a1,1​a1,2​​x2​+⋯=a1,1​a1,n+1​​，然后对于之后的每一个方程，都把第一个方程乘上当前这个方程的 x 1 x_1 x1​的系数，然后相减，就可以得到$n-1$个不含 x 1 x_1 x1​的方程，然后再对这剩下的$n-1$个方程，重复上述操作即可。最后先求解出 x n x_n xn​，然后不停地代到上一个方程里去解，最终得到整个方程组的解。

那么在代码实现的时候，是先把每个未知数的系数拎出来，形成一个矩阵(一般是$n\ast (n+1)$)。
为了代码美观方便起见，下标从0开始。
[ a 0 , 0 a 0 , 1 ⋯ a 0 , n a 1 , 0 a 1 , 1 ⋯ a 1 , n ⋮ ⋮ ⋮ ⋮ a n − 1 , 0 a n − 1 , 1 ⋯ a n − 1 , n ] \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} ⎣⎢⎢⎢⎡​a0,0​a1,0​⋮an−1,0​​a0,1​a1,1​⋮an−1,1​​⋯⋯⋮⋯​a0,n​a1,n​⋮an−1,n​​⎦⎥⎥⎥⎤​
美其名曰：增广矩阵。注意，最后一列是等号右边的常数项。

算法实现

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

• $step1.1$
  选取矩阵中的一行$row$，将含这个元的系数化为$1$。即
  [ 1 a 0 , 1 a 0 , 0 ⋯ a 0 , n a 0 , 0 a 1 , 0 a 1 , 1 ⋯ a 1 , n ⋮ ⋮ ⋮ ⋮ a n − 1 , 0 a n − 1 , 1 ⋯ a n − 1 , n ] \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} ⎣⎢⎢⎢⎡​1a1,0​⋮an−1,0​​a0,0​a0,1​​a1,1​⋮an−1,1​​⋯⋯⋮⋯​a0,0​a0,n​​a1,n​⋮an−1,n​​⎦⎥⎥⎥⎤​
  特别地，为了减小精度损失和代码方便，一般是找系数最大的一行。而如果系数最大的系数也是$0$，那么说明当前这个元的系数不用消，则直接进入下一个消元，并对此元进行标记(至于标记什么后面会提到)。
• $step1.2$
  对于后面的每一行$i$，都将第$row$的$col$(刚刚已经化为$1$)乘上第$i$行第$col$的系数。然后两行相减，即可以将$i$行的第$col$这一元消去。

枚举完后，增广矩阵应该是变成了上三角矩阵，即
[ ⋯ ⋯ ⋯ ⋯ 0 ⋯ ⋯ ⋯ ⋮ ⋮ ⋮ ⋮ 0 0 ⋯ x n − 1 ] \begin{bmatrix}\cdots&\cdots&\cdots&\cdots\\0&\cdots&\cdots&\cdots\\\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&x_{n-1}\end{bmatrix} ⎣⎢⎢⎢⎡​⋯0⋮0​⋯⋯⋮0​⋯⋯⋮⋯​⋯⋯⋮xn−1​​⎦⎥⎥⎥⎤​
有重要性质：一个位置为$0$，则其左下所有位置都为$0$。
然后求解

• $step2.1$
  先判断是否有解。
  那么如果无解，说明有某一次消元时得到了$0=x(x\ne 0)$的情况。所以在最后的增广矩阵中找未知数系数全为$0$但是等号右边不为$0$的行，即$000\cdots 0x(x\ne 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$取 1 0 − 12 10^{-12} 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.1int 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.2for(int i=row+1;i<equ;i++)if(fabs(a[i][col])>eps){double tmp=a[i][col]/a[row][col];//把当前方程的col元系数化为1for(int j=col;j<=var;j++)a[i][j]-=a[row][j]*tmp;//乘上col元的系数，消元a[i][col]=0;//避免精度影响，当前方程该元赋为0}}//step2for(int i=row;i<equ;i++)if(fabs(a[i][col])>eps)return -1;//step2.1if(var>row) return var-row;//step2.2//step2.3for(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
计算电压

异或方程组

此种题目是形如
{ a 0 , 0 x 0 ⊗ a 0 , 1 x 0 ⊗ ⋯ = a 0 , n a 1 , 0 x 1 ⊗ a 1 , 1 x 1 ⊗ ⋯ = a 1 , n ⋮ a n − 1 , 0 x n − 1 ⊗ a n − 1 , 1 x n − 1 ⊗ ⋯ = a n − 1 , n \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} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​a0,0​x0​⊗a0,1​x0​⊗⋯=a0,n​a1,0​x1​⊗a1,1​x1​⊗⋯=a1,n​⋮an−1,0​xn−1​⊗an−1,1​xn−1​⊗⋯=an−1,n​​
众所周知，异或相当于模二，所以消元的时候就是系数之间异或就可以了。
一般是关联开关问题，对于每个开关，比较前后状态和相邻的开关，可以得出方程。

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;//step1int col=0,fnum=0,row;for(row=0;row<equ&&col<var;row++,col++){//step1.1int 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.2for(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];}//step2for(int i=row;i<equ;i++)if(a[i][col]!=0)return -1;//step2.1int temp=var-row;if(row<var) return temp;//step2.2//step2.3for(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;//step1int col=0,row;for(row=0;row<equ&&col<var;row++,col++){//step1.1int 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.2for(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]);//为了避免精度损失，要用lcmif(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;}}//step2for(int i=row;i<equ;i++)if(a[i][col]!=0)return -1;//step2.1if(var>row) return var-row;//step2.2//step2.3for(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;//step1int col=0,row;for(row=0;row<equ&&col<var;row++,col++){//step1.1int 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.2for(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;}}//step2for(int i=row;i<equ;i++)if(a[i][col]!=0)return -1;//step2.1if(row<var) return var-row;//step2.2//step2.3for(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

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！