gausstmp
title: gausstmp
categories:
- ICPC
tags:
- null
abbrlink: dd749cdf
date: 2023-08-04 00:00:00
当我们枚举到第 $i$ 个变量，且前 $i$ 个变量都不是自由元时，如果第 $i\sim n$ 行中这个变量的系数均为 $0$，那么第 $i$ 个变量是自由元。
$$
\left[\begin{matrix}
a_{1,1} & a_{1,2} & a_{1, 3} & a_{1,4} & a_{1,5} \
0 & a_{2,2} & a_{2, 3} & a_{2,4} & a_{2,5} \
0 & 0 & \color{red}{0}& a_{3,4} & a_{3,5} \
0 & 0 & 0 & a_{4,4} & a_{4,5} \
0 & 0 & 0 & 0 & a_{5,5}
\end{matrix} \middle |
\begin{matrix}
c_{1} \
c_{2} \
c_{3} \
c_{4} \
c_{5}
\end{matrix} \right]
$$
上面的例子中，$x_3$ 即为自由元。可以发现，当 $x_4, x_5$ 被确定之后，第三行的方程可以记为：$$
0\cdot x_3=c_3-x_4\cdot a_{3, 4}-x_5\cdot a_{3, 5}=C
$$
若 $C=0$，则 $x_3$ 可以在 $\Bbb {R}$ 中任取。
若 $C\neq 0$，则 $x_3$ 一定无解，方程组也一定无解。
另外，存在一个自由元，其值可以任取，不代表方程组有无穷组解。因为若有其它变量是无解的，则方程组无解。
出现自由元这种情况，说明有些方程之间是线性相关的，所以加减消元之后就消没了。
解决 $n\neq m$ 的情况
在枚举变量的过程当中，会发现有些变量是自由元。这时我们应当跳过这一列，并且保持行不变，处理下一个变元。
因此我们用指针 $\operatorname {row}$ 和 $\operatorname {col}$ 表示当前所在的行和列，最后 $\operatorname {row}$ 的位置即为非自由元的个数。
由于遇到所有自由元时，我们都保持行不变，所以可以理解为：我们把所有自由元都挪到了最后面没有被 $\operatorname {row}$ 遍历到的几行当中去。
这可以结合上面 $5\times 5$ 的例子进行理解。
在遇到自由元 $x_3$ 后，跳过第 $3$ 列，保持 $\operatorname {row}$ 不变。所以第 $3$ 行保留的主元是 $x_4$，同理第 $4$ 行保留的主元是 $x_5$.
所以第 $5$ 行即代表第 $x_3$，第 $5$ 行的方程为：$0\times x_3=c_5$.
代码
下面代码会将所有自由元保留到最后的连续几行。
因此，要判断解的情况，只需要去判断第 $\operatorname {row+1}\sim n$ 行的系数 $c$ 是否为 $0$ 即可。
//高斯消元法，返回-1为无解，0为无穷解，1为唯一解
double A[MAXN][MAXN], ans[MAXN];
int Gauss() {
int row = 1, max_row;
for(int col = 1; col <= m; col++) {
max_row = row;
for(int i = row + 1; i <= n; i++) { //提高精度 找到一列中值最大的
if(fabs(A[i][col]) > fabs(A[max_row][col])) max_row = i;
}
if(fabs(A[max_row][col]) < eps)
continue; //如果col是自由元，那么跳过这一列，但行不变
if(max_row != row) { //swap(当前行，最大值所在的行)
for(int i = 1; i <= m + 1; i++) swap(A[row][i], A[max_row][i]);
}
for(int i = row + 1; i <= n; i++) {
double t = A[i][col] / A[row][col];
for(int j = 1; j <= m + 1; j++) {
A[i][j] -= A[row][j] * t;
}
}
row++;
}
row--;
if(row < n) {
for(int i = row + 1; i <= n; i++) {
//如果左边所有系数都是0，等号右边不是0，则无解
if(abs(A[i][m + 1]) > eps) return -1;
} //否则有无数解
return 0;
}
for(int i = m; i >= 1; i--) { //回代
for(int j = i + 1; j <= m; j++) A[i][m + 1] -= A[i][j] * ans[j];
ans[i] = A[i][m + 1] / A[i][i];
}
return 1;
}