[笔记]高斯消元

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

内容

求解如下的线性方程组（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];//对角线已经是 1
	for(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;
}
posted @ 2025-10-20 21:55  Sinktank  阅读(36)  评论(0)    收藏  举报
★CLICK FOR MORE INFO★ TOP-BOTTOM-THEME
Enable/Disable Transition
Copyright © 2023 ~ 2025 Sinktank - 1328312655@qq.com
Illustration from 稲葉曇『リレイアウター/Relayouter/中继输出者』,by ぬくぬくにぎりめし.