高斯消元

高斯消元

网上的题解一个比一个难看，还一人一个码风。

如果我知道谁还要学高斯消元，我会让他先做 [SDOI2006] 线性方程组 再做 【模板】高斯消元法，因为前者需要判出无解或无穷解，而后者只需要输出 No Solution 同时代表这两种情况，前者更通用一些，而写法却差得不少，从前者往后者转就是薄纱。

首先这玩意的思路就是代入消元和加减消元，只不过元的个数来到了 \(100\)，复杂度是 \(O(n^3)\) 的。注意最好别开中间变量，因为这玩意全程要用 double 运算，我每次开变量都习惯开成 int，不知道啥时候就开始错了。还有判断零，由于精度问题，我们设一个 \(eps = 10^{-7}\)，只要一个 double 的 fabs 比 \(eps\) 小我们就认为它是 \(0\) 了。

先放完整模板，再慢慢解释。

const int N = 1e2+5;
double eps = 1e-8;
double a[N][N];
double ans[N];
int n;
void gauss()
{
	int now{1};
	for(int i{1};i<=n;i++) // main 
	{
		for(int j{now};j<=n;j++)
		if(fabs(a[j][i]) > eps) // turn to now top
		{
			swap(a[j],a[now]);
			break;
		}
		if(fabs(a[now][i]) > eps)
		{
			for(int j{n+1};j>=i;j--) a[now][j] /= a[now][i];// to 1
			for(int j{now+1};j<=n;j++)
			{
				for(int k{n+1};k>=i;k--) a[j][k] -= a[now][k] * a[j][i];//+- delete
			}
			now++;
		}
	}
	if(now<n+1)
	{
		for(int i{now};i<=n;i++)
		{
			if(fabs(a[i][n+1]) > eps) puts("-1"),exit(0); //no solution
		}
		puts("0"),exit(0); // infinity solution
	}
	ans[n] = a[n][n+1];
	for(int i{n-1};i>=1;i--)
	{
		ans[i] = a[i][n+1];
		for(int j{i+1};j<=n;j++)
			ans[i] -= a[i][j] * ans[j];
	}
}
signed main()
{
	#ifdef LOCAL
	freopen("in.in","r",stdin);
	#endif
	n = read();
	for(int i{1};i<=n;i++)
		for(int j{1};j<=n+1;j++)
		a[i][j] = 1.0 * read();
	gauss();
	for(int i{1};i<=n;i++)
	{
		printf("x%lld=%.2Lf\n",i,ans[i]); 
	}
	exit(0);
}

\(\begin{Bmatrix} 3x & + & 2y &+& z &=&10 \\5x & + & y &+& 6z &=&25 \\2x & + & 3y &+& 4z &=&20\end{Bmatrix}\)

我们发现一个 \(n\) 元的线性同余方程组有 \(n\) 个方程，每个方程有 \(n+1\) 项。我们称

\(\begin{bmatrix} 3 & 2 & 1 \\5 & 1 & 6 \\2 & 3 & 4 \end{bmatrix}\)​ 为它的系数矩阵，

\(\begin{bmatrix} 3 & 2 & 1 &10 \\5 & 1 & 6 &25 \\2 & 3 & 4 &20\end{bmatrix}\) 为它的增广矩阵。

首先你需要读入一个增广矩阵。

for(int i{1};i<=n;i++)
	for(int j{1};j<=n+1;j++)
	a[i][j] = 1.0 * read();

然后就可以开始 gauss() 了。

int now{1};

首先我们定义一个 now，表示的是当前已经手动消掉了多少行。为什么是“手动”？因为有的行那一项天生系数就是 \(0\)，就不需要自己消，那么就可能会出现无解、无穷解的情况。为什么初始化为 \(1\)？因为我们要从第一行开始计算，这刚好可以是我们的初始行位置。

for(int i{1};i<=n;i++) // main 
	{
		for(int j{now};j<=n;j++)
		if(fabs(a[j][i]) > eps) // turn to now top
		{
			swap(a[j],a[now]);
			break;
		}
		if(fabs(a[now][i]) > eps)
		{
			for(int j{n+1};j>=i;j--) a[now][j] /= a[now][i];// to 1
			for(int j{now+1};j<=n;j++)
			{
				for(int k{n+1};k>=i;k--) a[j][k] -= a[now][k] * a[j][i];//+- delete
			}
			now++;
		}
	}

接下来这部分是加减消元。

我们外层枚举的 \(i\) 是主元，然后枚举 \(j\) 从 now 到 \(n\)，找到一个当前元系数非 \(0\) 的式子并把它提到 now 的这一行（当然如果没找到就 continue），然后把下面所有的行的 \(i\) 元全部消掉，并给 now 递增，这样就保证了 \(1\) 到 now-1 行都是被消过的。我们下一次枚举行只需要从递增过的 now 行枚举即可了。

怎么加减消元？先把当前行的当前元的系数置为 \(1\)，当前行的其他系数相应改变。接着我们就拿着这行去消下面所有的行的 \(i\) 元，因为当前行的 \(i\) 系数已经为 \(1\) 了，那么遇到哪行就相应扩那行的 \(i\) 系数倍，一项一项减就行了。注意行内是倒着枚举，因为要用到 \(i\) 元的系数，不能一上来就清零了。

if(now<n+1)
	{
		for(int i{now};i<=n;i++)
		{
			if(fabs(a[i][n+1]) > eps) puts("-1"),exit(0); //no solution
		}
		puts("0"),exit(0); // infinity solution
	}

这一部分是判断有无解的。如果你没有 continue 过，那么你的 now 经过 \(n\) 次递增肯定会变成 \(n+1\)，而现在 now 如果比 \(n+1\) 小，则说明有元没有手动消，就会出现无解无穷解的情况。

解的情况：

唯一解：\(i\) 可以枚举完 \(n\) 行 \(\begin{bmatrix} 1 & 1 & 1 & 6\\0 & 1 & 2 & 8 \\0 & 0 & 1 & 3 \end{bmatrix}\)

无解：\(a[i][i] = 0,b\neq 0\) \(\begin{bmatrix}1 &1 & 1 & 6 \\0 & 1 & 2 & 3 \\ 0 & 0 & 0&3\end{bmatrix}\)

无穷多解：\(a[i][i]=0,b = 0\) \(\begin{bmatrix}1&1&1&6\\0&1&2&3\\0&0&0&0\end{bmatrix}\)。

我们既然能进到这个 if 里，说明从 now 到 \(n\) 行之内的 \(a_{i,i}\) 都是 \(0\) 了，那么我们只需要 \(\text{check}\) 一下 \(a_{i,{n+1}}\) 的大小，注意无解的优先级是比无穷解高的。

ans[n] = a[n][n+1];
	for(int i{n-1};i>=1;i--)
	{
		ans[i] = a[i][n+1];
		for(int j{i+1};j<=n;j++)
			ans[i] -= a[i][j] * ans[j];
	}

之后从下往上回带就可以了。