高斯消元

posted on 2024-07-23 01:25:34 | under | source

其实介绍的是高斯-约旦消元法，即 \(\rm Gauss-Jordan\) 消元法，GJ。

思想：依次枚举每个主元 \(x_i\)，选取满足 \(x_i\) 系数非 \(0\) 的一行，将它交换至第 \(i\) 行，并消掉其它行 \(x_i\) 的系数。最终就可以得到一个对角矩阵啦。

注意一下，当我们枚举到 \(x_i\) 时，\(x_1\dots x_{i-1}\) 已经呈对角分布，且第 \(1\dots {i-1}\) 行必然存在 \(j<i\) 使得 \(x_{j}\ne 0\)。也就是说，在选取 \(x_i\) 的行时，不能选择前 \(i-1\) 行，否则，会导致 \(1\dots i-1\) 的系数又 \(\ne 0\) 了。

假如某一个主元找不到非 \(0\) 系数，那么一定是无正常解，对于 P3389 【模板】高斯消元法 而言，直接返回 No Solution 即可。

但是 P2455 [SDOI2006] 线性方程组 要求我们返回这到底是无解还是无穷解，一个 \(\rm naive\) 的想法是判断最终的对角矩阵，假如系数为 \(0\) 且常数为 \(0\)，无穷解；若系数为 \(0\) 且常数不为 \(0\)，无解。

可惜有一组 hack，上述算法会出现问题。

hack：

0 2 3
0 0 0

原程序返回无解，但显然是无穷解。错因：对于主元 \(x_1\)，我们直接跳过，并将原第一行视为“已被用过”，那么在枚举主元 \(x_2\) 时，也会跳过原第一行，而直接将原第二行视为 \(x_2\) 的。

也就是说，假如某个主元 \(x_p\) 不存在非 \(0\) 系数行，那么就不能将其视为“已被用过”，它应该要参与到之后主元的过程中。

仔细想想，这确实没错。还记得一开始说的为啥不能选取之前“已被用过”的行吗？因为之前行的前几个主元存在非 \(0\) 系数，用它来消会导致后面的行的这些系数也非 \(0\)。但是假如该主元根本不存在非 \(0\) 系数行，那么也就不会出现上述矛盾了~

所以特判一下即可。

模板代码

#include<bits/stdc++.h>
using namespace std;

const int N = 5e2 + 5;
const double eps = 1e-7;
int n, no_ans, many_ans;
double a[N][N];

signed main(){
	cin >> n;
	for(int i = 1; i <= n; ++i)
		for(int j = 1; j <= n + 1; ++j) scanf("%lf", &a[i][j]);
	for(int i = 1; i <= n; ++i){
		int p = 0; 
		for(int j = 1; j <= n; ++j)
			if(((j < i && fabs(a[j][j]) < eps) || j >= i) && fabs(a[j][i]) > eps){
				p = j;
				break;
			}
		if(!p) continue;
		swap(a[i], a[p]);
		for(int j = 1; j <= n; ++j){
			if(i == j || fabs(a[j][i]) < eps) continue;
			double x = a[j][i] / a[i][i];
			for(int k = i; k <= n + 1; ++k) a[j][k] -= x * a[i][k];
		}
	}
	for(int i = 1; i <= n; ++i){
		if(fabs(a[i][i]) < eps){
			if(fabs(a[i][n + 1]) < eps) many_ans = true;
			else no_ans = true;
		}
	}
	if(no_ans) {cout << -1; return 0;}
	if(many_ans) {cout << 0; return 0;}
	for(int i = 1; i <= n; ++i) printf("x%d=%.2lf\n", i, a[i][n + 1] / a[i][i]);
	return 0;
}