【笔记】高斯消元
用于求解 \(n\) 元一次方程组解的算法。
给定 \(n\) 个未知数 \(x_1,x_2,\dots,x_n\) 及若干个方程。
将其填入矩阵 \(A\)，\(A_{i,j}\) 表示第 \(i\) 个方程中 \(x_j\) 的系数，并将每个方程的常数项填入 \(A_{i,n+1}\)。
高斯消元
我们要一个一个加减消元，直到出现一个方程只有一个元，再依次回代求解。
假如我们现在要消 \(x_i\)，我们可以找一个 \(x_i\) 系数不为 \(0\) 的方程，左右两边同除这个系数，这样就得到了一个 \(x_i\) 系数为 \(1\) 的方程。
然后拿这个方程左右两边同乘其他方程的 \(x_i\) 系数，然后再与其相减，这样就能使其他方程的 \(x_i\) 系数变成 \(0\)。
最后只有选定的方程 \(x_i\) 系数是 \(1\) 了。重复此操作直到消到 \(x_n\)，期间我们保证消每个元时选的式子互不相同。最后如果有解，每个元都对应一个式子。
我们可以从上往下求解，按顺序选择算式，每次只消他下面的式子。这样，完成消元后就会得到一个上三角矩阵，第一行有所有未知数的系数，最后一行只有 \(x_n\) 的系数。
然后 \(x_n\) 就可以直接求了，求出来之后代入第 \(n-1\) 个式子中求得 \(x_{n-1}\)，以此类推直到全部求解。
这时候有个问题：如果从上往下求解过程中，下一个式子的 \(x_i\) 系数是 \(0\) 怎么办？
可以去找 \(x_i\) 系数绝对值最大的式子，与下一个式子交换，这样既可以规避 \(0\)，又可以尽量减小精度损失。如果整个方程都找不到非 \(0\) 系数的 \(x_i\)，那么方程无解。
int Gauss(){
for(int i=1;i<=n;i++){//当前消 x_i
int res=i;
for(int j=i+1;j<=n;j++){//寻找一个最大的系数
if(fabs(a[res][i])<fabs(a[j][i])) res=j;
}
swap(a[i],a[res]);//交换
if(fabs(a[i][i])<eps) return 0;//如果最大的都是 0，无解
double div=a[i][i];
for(int j=i;j<=n+1;j++) a[i][j]/=div;//同除 x_i 的系数，使其为 1
for(int j=i+1;j<=n;j++){
double tim=a[j][i];//同乘其他方程 x_i 的系数，使其相减得零
for(int k=i;k<=n+1;k++) a[j][k]-=tim*a[i][k];//加减消元
}
}
ans[n]=a[n][n+1];//求出 x_n
for(int i=n-1;i>=1;i--){//回代求解
ans[i]=a[i][n+1];
for(int j=i+1;j<=n;j++) ans[i]-=(a[i][j]*ans[j]);
}
return 1;
}
时间复杂度 \(O(n^3)\)。
高斯 - 约旦消元
和高斯消元差不多。区别是我们不去只消每个式子下面式子的元，而是消所有。这样最终就得到一个每一行只有一个未知数有系数的矩阵，\(x_i=\frac{A_{i,n+1}}{A_{i,i}}\)。
复杂度会更劣一些但是更容易实现。
void Gauss_Jordan(){
for(int i=1;i<=n;i++){
int res=i;
for(int j=i+1;j<=n;j++){
if(fabs(a[res][i])<fabs(a[j][i])) res=j;
}
swap(a[i],a[res]);
if(fabs(a[i][i])<eps) return 0;
double div=a[i][i];
for(int j=i;j<=n+1;j++) a[i][j]/=div;
for(int j=1;j<=n;j++){
if(j==i) continue;
double tim=a[j][i];
for(int k=i;k<=n+1;k++) a[j][k]-=tim*a[i][k];
}
}
for(int i=1;i<=n;i++) ans[i]=a[i][n+1]/a[i][i];
return 1;
}
例题
高斯消元
Luogu P2455 [SDOI2006] 线性方程组
高斯消元模板，区别在于需要区分无解和无穷解。
正常来说消到哪个判哪个没有什么问题，但是当方程组同时出现 \(0x=0\) 和 \(0x=y\) 时，需要输出无解，但如果先处理到 \(0x=0\)，会输出无穷解。
如果只加一个优先级判断，依然会是有问题的。
考虑如下一组样例：
答案应是无穷解。但我们在处理第一个式子的时候，就会认定他是最大的 \(x_1\) 系数，而不会在求解 \(x_2\) 最大系数时考虑。我们从下往上代入，发现求不出 \(x_2\)，上面就只能变成 \(0x_1=3\)，判定为无解。
修改的方式是，在找最大系数时不能只往后找，应该找所有的，还没被确定最大系数或者最大系数被确定为 \(0\) 的行。
int Gauss(){
for(int i=1;i<=n;i++){
int res=i;
for(int j=1;j<=n;j++){//此处即为修改
if(fabs(a[j][j])>eps&&j<i) continue;
if(fabs(a[j][i])>fabs(a[res][i])) res=j;
}
if(res!=i) swap(a[i],a[res]);
if(fabs(a[i][i])<eps) continue;
double div=a[i][i];
for(int j=i;j<=n+1;j++) a[i][j]/=div;
for(int j=i+1;j<=n;j++){
double tim=a[j][i];
for(int k=i;k<=n+1;k++) a[j][k]-=tim*a[i][k];
}
}
int flag1=0,flag2=0;
for(int i=n;i>=1;i--){
for(int j=i+1;j<=n;j++) a[i][n+1]-=a[i][j]*ans[j];
if(fabs(a[i][i])<eps){
if(fabs(a[i][n+1])<eps) flag1=1;
else flag2=1;
}
ans[i]=a[i][n+1];
}
if(flag2) return -1;
else if(flag1) return 0;
else return 1;
}
Luogu P2447 [SDOI2010] 外星千足虫
阅读题目发现即让求解模 \(2\) 下的方程组，且只关心解的奇偶。
因为只关心奇偶我们可以奇数为 \(1\) 偶数为 \(0\)，然后发现题目等价于异或。
异或的消元更为简单：系数为仅为 \(0\) 和 \(1\)，找一个第 \(i\) 列为 \(1\) 的方程跟其他所有第 \(i\) 列也为 \(1\) 的方程异或，就可以消去所有第 \(i\) 位。回代比较麻烦所以直接使用高斯约旦。
对于 \(k\) 的求解，每次尽量找第一个，然后取所有用过的行数 \(\max\)。
因为数据范围比较大我们使用 bitset 优化这个过程。
#include<bits/stdc++.h>
using namespace std;
const int N=1e3+10;
int n,m,ans;
bitset<N> a[2*N];
int gaus(){
for(int i=1;i<=n;i++){
int res=i;
for(int j=i;j<=m;j++){
if(a[j][i]==1){
res=j;
break;
}
}
if(!a[res][i]) return 0;
ans=max(ans,res);
if(res!=i) swap(a[res],a[i]);
for(int j=1;j<=m;j++){
if(i!=j&&a[j][i]) a[j]^=a[i];
}
}
for(int i=1;i<=n;i++){
if(!a[i][i]&&!a[i][n+1]) return 0;
}
return 1;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin>>n>>m;
for(int i=1;i<=m;i++){
for(int j=1;j<=n+1;j++){
char c;
cin>>c;
a[i][j]=bool(c-'0');
}
}
int d=gaus();
if(!d) cout<<"Cannot Determine\n";
else{
cout<<ans<<'\n';
for(int i=1;i<=n;i++){
if(a[i][n+1]) cout<<"?y7M#\n";
else cout<<"Earth\n";
}
}
return 0;
}
高斯消元求解 DP
一些 DP 的转移有后效性，可以使用高斯消元解方程组求解。
浙公网安备 33010602011771号