高斯消元法是求解线性方程组的经典算法。
内容
求解如下的线性方程组(P3389 【模板】高斯消元法):
我们考虑将所有变量丢到一个 \(n\times (n+1)\) 的矩阵中,其中最后一列存储常数项 \(b_1,\dots,b_n\):
我们称这个矩阵为增广矩阵。
高斯消元的理论依据:
- 交换律:将两行交换,方程的解不变。
- 加法律:将某行加上另外一行,方程的解不变。
- 乘法律:将某行乘一个非零常数,方程的解不变。
- 乘加律:将某行乘上一个常数,加到另一行上,方程的解不变。可由 \(2,3\) 推出。
我们的目标是利用这些性质,对增广矩阵做一些变换,使它变成下面的形态,一个下三角矩阵(Step 1):
最后再变成这样(Step 2):
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\) 个开关按下 / 不按下。
则所有约束条件可以用下面的异或方程组来描述:
异或可以理解为模 \(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;
}