自动搬运

查看原文

来自洛谷，原作者为

	Piwry 空想上の人格保持者

搬运于2025-08-24 21:37:52，当前版本为作者最后更新于2020-01-12 18:54:11，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

作为刚刚刷完模板题想来练练手的蒟蒻，还没摆脱模板题top1题解的影响，于是就想拿 高斯-约旦消元法 来过这道题。

然而当我读完题后，才发现 高斯-约旦消元法 那美妙的"\"字形在这题并不能体现，于是菜的不行的窝就去翻题解...

遗憾的是，当时我只找到了一篇来自小诗乃的题解，可这篇题解的代码也有很多冗余的计算。

最后，我就只好自力更生了

• 我们先来回顾一下用于消元的工具 增广矩阵 ；以及它的两种变换：交换行，对两行数字做加减。这两种变换后得到的矩阵从求解方程的意义上是没有变化的

• 现在再来回顾 高斯-约旦消元法 （假设有唯一解）的步骤：

1. 以列遍历，从一侧到另一侧；为了方便操作，通常从上至下遍历行，将操作行转移至当前行
2. 先从当前操作行及尚未操作的行找到一个主元（通常取当前列系数最大的，可以证明这样精度误差最小（虽然我不会证 简略解释可见文章结尾）），再从除该行的所有行里消去这个元
3. 重复 2.，直至遍历完毕

如果你手推过或仔细观察过它的过程，你会发现这个算法有两点重要的性质（为了描述方便，以下均默认用上文描述为“通常”的操作方式）：

• 每次消元完后，当前列只有当前行不为 $0$
• 当前列左边的数均为 $0$

现在我们回到这道题；因为有唯一解的情况可以直接用裸的 高斯-约旦消元法 做出，显然其难点就在于无解和无穷解的判定

而无解和无穷解情况的不同就在于，做到某一行你会找不到主元（即当前及下面的行的该列均为 $0$）

假设当我们第一次遇到找不到主元的情况时

我们可以发现：

• 当前列左边的数仍均为 $0$
• 当前操作列及其下面的数均为 $0$

于是便想到转而处理下一列，操作行不变，以维持这样的性质

现在我们对所有找不到主元的情况都按上述方法进行处理，当我们遍历完了所有列，当前矩阵也仍满足性质：

• 每次操作列左边的数均为 $0$
• 每次操作列从那次操作行下面一行开始全为 $0$

于是发现它就是一个“坡度不一”的“倒三角”（不考虑等式右边常数项），对于未操作的行（即“三角”下面的行），其中的系数全为 $0$！

这样的“倒三角”中可以确定最下面的一个未知数（或可能已经给出解），往上一一带回一定可以求解；
也就是说只要判定下面未操作行是否会导出矛盾就可以对整个方程组做出判定了！

那么就好办了，反正系数全为 $0$，常数项只要都为 $0$ 就是无穷解，否则就是无解

（至于不存在找不到主元的情况，直接按通常的 高斯-约旦消元法 做就可以了，一定是存在唯一解的）

$\color{#ffb464}\texttt{顺便推下我的}$$\color{#32dddd}\texttt{博客}$

参考代码

#include <cstdio>
#include <cmath>
#include <iostream>
using namespace std;

double A[50][51];
int N;

double Abs(double x){
	if(x < 0) return -x;
	else return x;
}

bool eq(double x, double y){// 因为存在微小的精度误差，两数的相差值在一定范围内即可视为相等 
	return (Abs(x-y) < 1e-9);
}

int main(){
	scanf("%d", &N);
	for(int i =0; i < N; ++i) for(int j =0; j < N+1; ++j)
		scanf("%lf", &A[i][j]);
	int nwline =0;
	// k 指主元序号（列）
	for(int k =0; k < N; ++k){// 需要考虑无穷解，循环到N 
		int maxi =nwline;
		for(int i =nwline+1; i < N; ++i)
			if(Abs(A[i][k]) > Abs(A[maxi][k])) maxi =i;
		if(eq(A[maxi][k], 0))
			continue;
		for(int j =0; j < N+1; ++j)
			swap(A[nwline][j], A[maxi][j]);
		for(int i =0; i < N; ++i){
			if(i == nwline)
				continue;
			double mul =A[i][k]/A[nwline][k];
			for(int j =k; j < N+1; ++j)
				A[i][j] -=A[nwline][j]*mul;
		}
		++nwline;
	}
	if(nwline < N){// 存在找不到主元的情况
		while(nwline < N)
			if(!eq(A[nwline++][N], 0))
				return puts("-1") && 0;
		putchar('0');
	}
	else
		for(int i =0; i < N; ++i)
			printf("x%d=%.2lf\n", i+1, A[i][N]/A[i][i]);
}

upd at 20240227_0113

诈尸修了下判断浮点数相等处的万年老锅（具体见代码中 eq 函数），以及重新润写了下文章正文；

另外上文提到的：

2. 先从当前操作行及尚未操作的行找到一个主元（通常取当前列系数最大的，可以证明这样精度误差最小（虽然我不会证 简略说明可见文章结尾）），再从除该行的所有行里消去这个元

这里的原理和误差传播有关；简单来说，输入误差相同时，大数字除小数字比起小数字除大数字结果误差更大。于是“取当前列系数最大的”可以有效减少后面代码中 mul 变量的精度误差

想更深入了解的读者可以以“误差传播”为关键词在互联网上搜索