[笔记]高斯消元

高斯消元法是求解线性方程组的经典算法。

内容

求解如下的线性方程组（P3389 【模板】高斯消元法）：

\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+\dots+a_{1,n}x_n=b_1\\ a_{2,1}x_1+a_{2,2}x_2+\dots+a_{2,n}x_n=b_2\\ \dots\\ a_{n,1}x_1+a_{n,2}x_2+\dots+a_{n,n}x_n=b_n \end{cases} \]

我们考虑将所有变量丢到一个 \(n\times (n+1)\) 的矩阵中，其中最后一列存储常数项 \(b_1,\dots,b_n\)：

\[\begin{bmatrix} a_{1,1}&a_{1,2}&\cdots&a_{1,n}&b_1\\ a_{2,1}&a_{2,2}&\cdots&a_{2,n}&b_2\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ a_{n,1}&a_{n,2}&\cdots&a_{n,n}&b_n \end{bmatrix} \]

我们称这个矩阵为增广矩阵。

高斯消元的理论依据：

1. 交换律：将两行交换，方程的解不变。
2. 加法律：将某行加上另外一行，方程的解不变。
3. 乘法律：将某行乘一个非零常数，方程的解不变。
4. 乘加律：将某行乘上一个常数，加到另一行上，方程的解不变。可由 \(2,3\) 推出。

我们的目标是利用这些性质，对增广矩阵做一些变换，使它变成下面的形态，一个下三角矩阵（Step 1）：

\[\begin{bmatrix} a'_{1,1}&a'_{1,2}&a'_{1,3}&\cdots&a'_{1,n}&b'_1\\ 0&a'_{2,2}&a'_{2,3}&\cdots&a'_{2,n}&b'_2\\ 0&0&a'_{3,3}&\cdots&a'_{3,n}&b'_3\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&a'_{n,n}&b'_n \end{bmatrix} \]

最后再变成这样（Step 2）：

\[\begin{bmatrix} 1&0&0&\cdots&0&b''_1\\ 0&1&0&\cdots&0&b''_2\\ 0&0&1&\cdots&0&b''_3\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&1&b''_n \end{bmatrix} \]

Step 1

我们依次遍历 \(i=1,2,\dots,n\)，对于第 \(i\) 行，我们选定 \(a_{i,i}\) 为主元（\(a_{i,i}\ne 0\)）用乘加律将 \(i+1,\dots,n\) 这些行的第 \(i\) 列全部消成 \(0\)。

这样进行下去就能获得下三角矩阵了。

需要注意，若 \(a_{i,i}=0\)，则上面的步骤无法进行。需要向下找一个 \(j\) 使得 \(a_{j,i}\ne 0\)，然后利用交换律将 \(i,j\) 行交换，再继续进行。

若无法找到这样的 \(j\)，则说明无解，或有无穷多解。

Step 2

从最后一行开始，逐个回带即可。相当于用第 \(n\) 行将第 \(n-1\) 行的第 \(n\) 列消成 \(0\)……以此类推。

• 若回带过程中遇到 \(0x=k(k\ne 0)\) 则说明无解。
• 在有解得情况下，若遇到 \(0x=0\) 则说明有无穷多解。
• 否则有唯一解。

---

通常高斯消元会用 double 存储结果，为了减小精度误差，两个值差距在一定范围内（如 \(10^{-7}\)）即视作相等。

另外为了减小误差，通常选取 \(a_{j,i}\) 最大的 \(j\) 和第 \(i\) 行进行交换（这样做的可行性在这里有更详细的说明）。

代码实现的 Step 1，先将第 \(i\) 行整体乘法，使得 \(a_{i,i}=1\) 再进行其他操作，这样乘加时就不需要考虑 \(a_{i,i}\) 的取值，且 Step 2 也不需要再考虑系数问题了。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
const int N=105;
double a[N][N],ans[N];
int n;
signed main(){cin>>n;for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) cin>>a[i][j];for(int i=1;i<=n;i++){//枚举主元序号（即列的编号）int r=i;for(int j=i+1;j<=n;j++){if(abs(a[j][i])>abs(a[r][i])) r=j;}if(abs(a[r][i])<eps) cout<<"No Solution\n",exit(0);//找不到非 0 的行，无解或无穷多解if(i^r) swap(a[i],a[r]);double div=a[i][i];for(int j=i;j<=n+1;j++) a[i][j]/=div;//使 a[i][i]=1，更方便一些for(int j=i+1;j<=n;j++){div=a[j][i];for(int k=i;k<=n+1;k++) a[j][k]-=a[i][k]*div;}}ans[n]=a[n][n+1];//对角线已经是 1for(int i=n-1;i;i--){ans[i]=a[i][n+1];for(int j=i+1;j<=n;j++) ans[i]-=a[i][j]*ans[j];}for(int i=1;i<=n;i++)cout<<fixed<<setprecision(2)<<ans[i]<<"\n";return 0;
}

高斯-约旦消元

高斯-约旦消元相比高斯消元，通过对 Step 1 进行少量修改，可以直接使原矩阵变换成第二个矩阵，从而省掉 Step 2。

具体地，对于第 \(i\) 行，除了将 \(i+1,\dots,n\) 行的第 \(i\) 列消成 \(0\)，我们把 \(1,\dots,i-1\) 行的第 \(i\) 列也消成 \(0\)。这样我们直接就得出了第二个矩阵，直接访问 \(f_{i,n+1}\) 即可获得 \(x_i\) 的值了。

码量小了不少。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
const int N=105;
double a[N][N];
int n;
signed main(){cin>>n;for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) cin>>a[i][j];for(int i=1;i<=n;i++){//枚举主元序号（即列的编号）int r=i;for(int j=i+1;j<=n;j++){if(abs(a[j][i])>abs(a[r][i])) r=j;}if(abs(a[r][i])<eps) cout<<"No Solution\n",exit(0);if(i^r) swap(a[i],a[r]);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;div=a[j][i];for(int k=i;k<=n+1;k++) a[j][k]-=a[i][k]*div;}}for(int i=1;i<=n;i++)cout<<fixed<<setprecision(2)<<a[i][n+1]<<"\n";return 0;
}

题

P2455 [SDOI2006] 线性方程组

板子题，和刚才的 P3389 区别仅在与无解和无穷多解需要不同的输出。

然而我们刚才的代码会 WA 90pts。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
const int N=55;
double a[N][N];
int n;
bool infi; 
signed main(){cin>>n;for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) cin>>a[i][j];for(int i=1;i<=n;i++){//枚举主元序号（即列的编号）int r=i;for(int j=i+1;j<=n;j++){if(abs(a[j][i])>abs(a[r][i])) r=j;}if(abs(a[r][i])<eps) continue;//无解或无穷多解 if(i^r) swap(a[i],a[r]);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;div=a[j][i];for(int k=i;k<=n+1;k++) a[j][k]-=a[i][k]*div;}}for(int i=1;i<=n;i++){if(abs(a[i][i])<eps){if(abs(a[i][n+1])>eps) cout<<"-1\n",exit(0);infi=1;}}if(infi) cout<<"0\n",exit(0);for(int i=1;i<=n;i++)cout<<"x"<<i<<"="<<fixed<<setprecision(2)<<a[i][n+1]<<"\n";return 0;
}

其原因在于，当我们发现第 \(i\) 列不存在非零的元素后，我们会 continue 掉，继续对第 \(i+1\) 行进行处理。

前 \(i\) 行中只是第 \(i\) 列不存在非零的元素，第 \(i+1\) 列是可能存在非零元素的，但是它被我们 continue 掉了，不会再处理。这就可能导致无穷多解被判成无解。

可以参考这组 Hack。

2
0 2 3
0 0 0

我们的对策也很简单，只需要在找第 \(i\) 列非零元素时，一并将第 \(1,\dots,i-1\) 行也遍历一遍就可以了。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
const int N=55;
double a[N][N];
int n;
bool infi; 
signed main(){cin>>n;for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) cin>>a[i][j];for(int i=1;i<=n;i++){//枚举主元序号（即列的编号）int r=i;for(int j=1;j<=n;j++){if(abs(a[j][j])>eps&&j<i) continue;//行 j 的在第 j 列已经存在主元，则不可用if(abs(a[j][i])>abs(a[r][i])) r=j;}if(abs(a[r][i])<eps) continue;//无解或无穷多解if(i^r) swap(a[i],a[r]);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;div=a[j][i];for(int k=i;k<=n+1;k++) a[j][k]-=a[i][k]*div;}}for(int i=1;i<=n;i++){if(abs(a[i][i])<eps){if(abs(a[i][n+1])>eps) cout<<"-1\n",exit(0);infi=1;}}if(infi) cout<<"0\n",exit(0);for(int i=1;i<=n;i++)cout<<"x"<<i<<"="<<fixed<<setprecision(2)<<a[i][n+1]<<"\n";return 0;
}

另一种方法是将构成阶梯的行都 swap 到上面去，这样每次向下遍历即可保证正确性。更简洁一些，我更中意这个。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
const int N=55;
double a[N][N];
int n,idx=1;
signed main(){cin>>n;for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) cin>>a[i][j];for(int i=1;i<=n;i++){//枚举主元序号（即列的编号）int r=idx;for(int j=r+1;j<=n;j++){if(abs(a[j][i])>abs(a[r][i])) r=j;}if(abs(a[r][i])<eps) continue;//无解或无穷多解if(idx^r) swap(a[idx],a[r]);double div=a[idx][i];for(int j=i;j<=n+1;j++) a[idx][j]/=div;for(int j=1;j<=n;j++){if(j==idx) continue;div=a[j][i];for(int k=i;k<=n+1;k++) a[j][k]-=a[idx][k]*div;}idx++;}if(idx<=n){while(idx<=n) if(abs(a[idx++][n+1])>eps) cout<<"-1\n",exit(0);cout<<"0\n";}else{for(int i=1;i<=n;i++)cout<<"x"<<i<<"="<<fixed<<setprecision(2)<<a[i][n+1]<<"\n";}return 0;
}

P10499 开关问题

若按下开关 \(j\) 会改变开关 \(i\) 的状态，则令 \(a_{i,j}=1\)，否则为 \(0\)。

同时，令 \(x_i=1/0\) 分别为 \(i\) 个开关按下 / 不按下。

则所有约束条件可以用下面的异或方程组来描述：

\[\begin{cases} a_{1,1}x_1\oplus a_{1,2}x_2\oplus \dots\oplus a_{1,n}x_n=s_1\oplus t_1\\ a_{2,1}x_1\oplus a_{2,2}x_2\oplus \dots\oplus a_{2,n}x_n=s_2\oplus t_2\\ \dots\\ a_{n,1}x_1\oplus a_{n,2}x_2\oplus \dots\oplus a_{n,n}x_n=s_n\oplus t_n \end{cases} \]

异或可以理解为模 \(2\) 意义下的加法，所以可以用高斯消元做。

考虑到两个 bitset 可以直接异或，所以写起来会方便很多。

或者利用 \(n\) 很小，将每行的 \(a\) 压缩到一个整数里也行。

时间复杂度为 \(O(\dfrac{n^3}{\omega})\) 或者 \(O(n^2)\)。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const int N=32;
int k,n,a[N],g[N];
inline int gauss(){for(int i=0;i<n;i++){for(int j=i+1;j<n;j++){if(a[j]&(1<<i)){swap(a[i],a[j]),swap(g[i],g[j]);break;}}if(a[i]&(1<<i))for(int j=0;j<n;j++){if((i^j)&&(a[j]&(1<<i))) a[j]^=a[i],g[j]^=g[i];}}int ans=0;for(int i=0;i<n;i++){if(!a[i]){if(g[i]) return -1;else ans++;}}return 1<<ans;
}
signed main(){ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);cin>>k;while(k--){cin>>n;for(int i=0;i<n;i++) cin>>g[i],a[i]=(1<<i);for(int i=0,x;i<n;i++) cin>>x,g[i]^=x;int u,v;while(cin>>u>>v){if(!u) break;u--,v--;a[v]|=(1<<u);}int ans=gauss();if(~ans) cout<<ans<<"\n";else cout<<"Oh,it's impossible~!!\n";}return 0;
}

这就是高斯消元求解异或方程组。

P5027 Barracuda

我们可以枚举不合法的方程，然后判断解的情况即可。

具体来说，必须恰有一种方案有合法解。其中合法解要满足：

• 所有解均为正整数。
• 最大值唯一。

时间复杂度 \(O(n^4)\)。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const int N=105;
const double eps=1e-7;
int n,idx,b[N][N],ans;
double a[N][N];
inline int gauss(){idx=1;for(int i=1;i<=n;i++){int r=idx;for(int j=r+1;j<=n;j++){if(abs(a[j][i])>abs(a[r][i])) r=j;}if(abs(a[r][i])<eps) return 0;if(idx^r) swap(a[idx],a[r]);double div=a[idx][i];for(int j=i;j<=n+1;j++) a[idx][j]/=div;for(int j=1;j<=n;j++){if(j==idx) continue;div=a[j][i];for(int k=i;k<=n+1;k++) a[j][k]-=a[idx][k]*div;}idx++;}int mx=0,f=0,p=0;for(int i=1,t;i<=n;i++){t=a[i][n+1]+0.5;if(a[i][n+1]<eps||abs(a[i][n+1]-t)>eps){return 0;}if(t>mx) mx=t,p=i;}for(int i=1;i<=n;i++) if(abs(a[i][n+1]-mx)<eps){if(f) return 0;f=1;}return p;
}
signed main(){cin>>n;for(int i=1,m,x;i<=n+1;i++){cin>>m;while(m--) cin>>x,b[i][x]=1;cin>>x,b[i][n+1]=x;}for(int x=1;x<=n+1;x++){for(int i=1,idx=0;i<=n+1;i++){if(i==x) continue;++idx;for(int j=1;j<=n+1;j++) a[idx][j]=b[i][j];}int t=gauss();if(t){if(ans) cout<<"illegal\n",exit(0);ans=t;}}if(ans) cout<<ans<<"\n";else cout<<"illegal\n";return 0;
}