题解 P2455 【[SDOI2006]线性方程组】
Piwry
2020-01-12 18:54:11
作为刚刚刷完模板题想来练练手的蒟蒻，还没摆脱模板题top1题解的影响，于是就想拿 _高斯-约旦消元法_ 来过这道题。
然而当我读完题后，才发现 _高斯-约旦消元法_ 那美妙的"\\"字形在这题并不能体现，于是~~菜的不行的~~窝就去翻题解...
遗憾的是，当时我只找到了一篇来自~~小~~**诗乃**的[**题解**](https://www.luogu.com.cn/blog/asadashino/solution-p2455)，可这篇题解的代码也有很多**冗余的计算**。
最后，我就只好自力更生了
---
- 我们先来回顾一下用于消元的工具 **增广矩阵** ；以及它的两种变换：交换行，对两行数字做加减。这两种变换后得到的矩阵从求解方程的意义上**是没有变化的**
- 现在再来回顾 _高斯-约旦消元法_ （**假设有唯一解**）的步骤：
1. 以列遍历，从一侧到另一侧；为了方便操作，通常从上至下遍历行，将操作行转移至当前行
2. 先从当前操作行及尚未操作的行找到一个主元（通常取当前列系数最大的，可以证明这样精度误差最小（~~虽然我不会证~~ 简略解释可见文章结尾）），再从除该行的所有行里消去这个元
3. 重复 `2.`，直至遍历完毕
如果你手推过或仔细观察过它的过程，你会发现这个算法**有两点重要的性质**（为了描述方便，以下均默认用上文描述为“_通常_”的操作方式）：
- 每次消元完后，当前列只有当前行不为 $0$
- 当前列左边的数均为 $0$
---
现在我们回到这道题；因为有唯一解的情况可以直接用裸的 _高斯-约旦消元法_ 做出，显然其难点就在于**无解**和**无穷解**的判定
而无解和无穷解情况的**不同**就在于，做到某一行你会找不到主元（即当前及下面的行的该列**均为** $0$）
假设当我们**第一次**遇到找不到主元的情况时
我们可以发现：
- 当前列左边的数仍均为 $0$
- 当前操作列及其下面的数均为 $0$
于是便想到**转而处理下一列，操作行不变**，以**维持**这样的性质
---
现在我们对所有找不到主元的情况**都**按上述方法进行处理，当我们遍历完了所有列，当前矩阵也仍满足性质：
- 每次操作列左边的数均为 $0$
- 每次操作列从那次操作行下面一行开始全为 $0$
于是发现它就是一个“**坡度不一”的“倒三角**”（不考虑等式右边常数项），对于**未操作的行**（即“三角”下面的行），**其中的系数全为 $0$！**
这样的“**倒三角**”中可以确定最下面的一个未知数（或可能已经给出解），**往上一一带回一定可以求解**；
\
也就是说**只要判定下面未操作行是否会导出矛盾**就可以对整个方程组做出判定了！
那么就好办了，反正**系数全为** $0$，常数项只要都为 $0$ 就是**无穷解**，否则就是**无解**
（至于**不**存在找不到主元的情况，直接按通常的 _高斯-约旦消元法_ 做就可以了，一定是**存在唯一解**的）
---
$\color{#ffb464}\texttt{顺便推下我的}$[$\color{#32dddd}\texttt{博客}$](https://www.luogu.com.cn/blog/105254/)
---
### 参考代码
```cpp
#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` 变量的精度误差
想更深入了解的读者可以以“误差传播”为关键词在互联网上搜索