高斯消元

upd：2024/7/14 几乎是重构本文章。

upd：2024/7/15 修改了本文章的部分描述。

前言：参考自《高等代数（上册）》。本文适用于中国宝宝体质食用。

一、算法介绍

\[\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases} \]

定义1： 一个 \(n \times m\) 的矩阵 \(A\) 如果满足以下条件就被称为阶梯形矩阵。

• 零行（元素全为 \(0\) 的行）在最下方。
• 主元（非零行的第一个不为 \(0\) 的元素）的列随行的增大而增大。

定义2： 将以下三种变换称为矩阵的初等行变换。

• 1° 型：某行加上另外一行的某个倍数
• 2° 型：交换任意两行
• 3° 型：用一个非零数乘上某一行

命题1： 将初等行变换把线性方程组变成同解的方程组。以 1° 型为例，设原方程组为 \(A\)，变换一次后为 \(B\)，容易证明 \(A\) 的每一个解都是 \(B\) 的解，将 \(B\) 变换回去又可以证明 \(B\) 的每一个解都是 \(A\) 的解，因此 \(A\) 与 \(B\) 同解。

命题2： 任意一个矩阵都可以通过初等行变换化成阶梯形矩阵。这个可以通过数学归纳法证明。

当我们把方程组（设有 \(m\) 个方程，\(n\) 个未知数）化为阶梯形矩阵后，就不难得出解了。如果存在 \(0=d \ (d \ne 0)\) 的行，原方程组无解。否则设 \(r\) 为阶梯形矩阵中不为 \(0=0\) 型的行的个数。分为以下两种情况：

• \(r=n\) ：此时原方程组有唯一解，且前 \(n\) 行一定为上三角的形式，从下往上依次带入计算即可。
• \(r<n\) ：我们称不为 \(0=0\) 型的行的主元为系数的未知量为主变量，其余未知量称为自由未知量。可以发现对自由未知量任意赋一组值后，所有主变量的值也就唯一确定了。因此这样的方程组具有无穷个解。

除了通过将矩阵变换为阶梯形矩阵来求解外，还可以变换为简化阶梯形矩阵。

定义3：如果一个阶梯形矩阵满足其所有主元全为 \(1\)，且每个主元所在的列的其余元素均为 \(0\)，那么称它为简化阶梯形矩阵。

这样做的优点是当 \(r=n\) 时不需要回带。以上便是高斯消元的两种方法。该算法的时间复杂度为 \(O(n^2m)\) 的。

以下就是 P2455 [SDOI2006] 线性方程组 的代码（注释里附有实现细节）：

#include<bits/stdc++.h>
#define int long long
#define N 110
#define lb long double
using namespace std;
int n;
lb a[N][N], b[N]; //a[i][j]表示第i个方程,第j个系数 
lb x[N];
lb Abs(lb x) {
	return x > 0 ? x : -x;
}
void My_work() { //第二种消元，也被称为约旦消元QAQ 

	int r = 1; //r表示0=0和0=d(d!=0)的行的个数（具体实现时考虑到有无解的情况与文章定义略有不同）。
	for(int i = 1; i <= n; i++) { //每个未知数 
		lb Mx = 0;
		int id = 0;
		for(int j = r; j <= n; j++)	{ //第j个方程 
			if(Abs(a[j][i]) > Mx) { //选一个绝对值最大的（据说这样做能降低精度）
				Mx = a[j][i];
				id = j;
			}
		}
		if(id == 0) continue;
		for(int j = i; j <= n + 1; j++) swap(a[r][j], a[id][j]); //把找到的行放在第r行

		for(int j = 1; j <= n; j++) { //将第j个方程的第i个方程的系数消成0 
			if(a[j][i] == 0 || j == r) continue;
			lb Num = 1.0 * a[j][i] / a[r][i];
			for(int k = 1; k <= n + 1; k++) a[j][k] -= a[r][k] * Num;
		}
		/*for(int i = 1; i <= n; i++) {
			for(int j = 1; j <= n + 1; j++) printf("%.2Lf ", a[i][j]);
			printf("\n");
		}
		printf("\n");*/
		r++;
	}r--;
	//printf("%lld\n",r);
	if(r < n) {
		for(int i = r + 1; i <= n; i++) {
			if(Abs(a[i][n + 1]) > 1e-9) { //无解
				cout << -1 << endl;
				return;
			}
		}
		cout << 0 << endl; //无数解
		return;
	}
	for(int i = 1; i <= n; i++) 
		printf("x%lld=%.2Lf \n", i, a[i][n + 1] / a[i][i]);
}

signed main() {
	//ios::sync_with_stdio(0); 
	//cin.tie(0); cout.tie(0);
	scanf("%lld", &n);
	for(int i = 1; i <= n; i++) {
		for(int j = 1; j <= n + 1; j++) scanf("%Lf", &a[i][j]);

	}
	My_work();
	return 0;

}

二、例题讲解

1、P3164 [CQOI2014] 和谐矩阵

题意：

我们称一个由 \(0\) 和 \(1\) 组成的矩阵是和谐的，当且仅当每个元素都有偶数个相邻的 \(1\)。一个元素相邻的元素包括它本身，及他上下左右的 \(4\) 个元素（如果存在）。给定矩阵的行数和列数，请计算并输出一个和谐的矩阵。注意：所有元素为 \(0\) 的矩阵是不允许的。

\(1\le n,m\le 40\)。

分析：

构造异或线性方程组，对高斯消元稍加修改就能支持了。利用 bitset 优化，时间复杂度 \(O(\frac{n^2m^2}{w})\)，但暴力消元也能通过本题。

2、P2447 [SDOI2010] 外星千足虫

题意：

现在你面前摆有 \(1\ldots N\) 编号的 \(N\) 只千足虫，你的任务是鉴定每只虫子所属的星球，但不允许亲自去数它们的足。

Charles 每次会在这 \(N\) 只千足虫中选定若干只放入“昆虫点足机”（the Insect Feet Counter, IFC）中，“点足机”会自动统计出其内所有昆虫足数之和。Charles 会将这个和数 \(\bmod\) \(2\) 的结果反馈给你，同时告诉你一开始放入机器中的是哪几只虫子。

他的这种统计操作总共进行 \(M\) 次，而你应当尽早得出鉴定结果。

假如在第 \(K\) 次统计结束后，现有数据就足以确定每只虫子的身份，你就还应将这个 \(K\) 反馈给 Charles，此时若 \(K<M\)，则表明那后 \(M-K\) 次统计并非必须的。

如果根据所有 \(M\) 次统计数据还是无法确定每只虫子身份，你也要跟 Charles 讲明：就目前数据会存在多个解。

对于 \(100\%\) 的数据，满足 \(1\leq N\leq 10^3\)，\(1\leq M\leq 2\times 10^3\)。

分析：

\(\mod2\) 这个条件同样可以转异或方程组。保证使用的方程最少，只需要选主元时贪心取编号最小的方程即可。同样也能用 bitset 优化。