数论

高斯消元

Posted by Monad on July 17, 2018

引入

我们把含有 n 个未知数且最高次数为 1 的方程称为“n 元线性方程”。

如 $ 3x_1 + 5x_2 - 2x_3 = 5 $ 就是一个 3 元线性方程。

然后我们把若干个 n 元线性方程列在一起，就构成了一个 n 元线性方程组。

如 $ \begin{cases} x + 2y = 5 \\ 2x - y = 1 \end{cases} $ 就是一个 2 元线性方程组。

假设我们现在有 m 条方程以及 n 个未知数，我们要把这 n 个未知数的值求出来。
我们可以简单地先对这 m 条线性方程进行简单的化简，让每一条线性方程都化简成 $a_1x_1 + a_2x_2 + a_3x_3 + … + a_nx_n = c$ 的形式。
这样，这个方程组就写成这种形式：

$$ \begin{cases} a_{1,1}x_1 + a_{1,2}x_2 + a_{1,3}x_3 + \cdots + a_{1,n}x_n &= c \\ a_{2,1}x_1 + a_{2,2}x_2 + a_{2,3}x_3 + \cdots + a_{2,n}x_n &= c \\ a_{3,1}x_1 + a_{3,2}x_2 + a_{3,3}x_3 + \cdots + a_{3,n}x_n &= c \end{cases} $$

但是，有意义的部分就只有系数 a 和 c。
我们可以把它们单独取出来，构成一个矩阵：

$$ \left [ \begin{array}{ccc|c} a_{1,1} & \cdots & a_{1,n} & c_1 \\ a_{2,1} & \cdots & a_{2,n} & c_2 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m,1} & \cdots & a_{m,n} & c_m \end{array} \right ] $$

高斯消元

高斯消元就是在矩阵上进行运算，从而得出原方程的解的算法。

矩阵的变换

要解出方程，我们要对这个矩阵进行 3 个基本操作，这 3 个基本操作是不会改变原方程的解的。

  1. 让某一行乘一个非零的常数 k
    $ a_1x_1 + a_2x_2 + \cdots + a_nx_n = c \Leftrightarrow ka_1x_1 + ka_2x_2 + \cdots + ka_nx_n = kc $
    矩阵形式就是 $ \left [ \begin{array}{ccc|c} a_{1,1} & \cdots & a_{1,n} & c_1 \\ a_{2,1} & \cdots & a_{2,n} & c_2 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m,1} & \cdots & a_{m,n} & c_m \end{array} \right ] \Leftrightarrow \left [ \begin{array}{ccc|c} ka_{1,1} & \cdots & ka_{1,n} & kc_1 \\ a_{2,1} & \cdots & a_{2,n} & c_2 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m,1} & \cdots & a_{m,n} & c_m \end{array} \right ] $

  2. 交换两行
    即 $ \left [ \begin{array}{ccc|c} a_{1,1} & \cdots & a_{1,n} & c_1 \\ a_{2,1} & \cdots & a_{2,n} & c_2 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m,1} & \cdots & a_{m,n} & c_m \end{array} \right ] \Leftrightarrow \left [ \begin{array}{ccc|c} a_{2,1} & \cdots & a_{2,n} & c_2 \\ a_{1,1} & \cdots & a_{1,n} & c_1 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m,1} & \cdots & a_{m,n} & c_m \end{array} \right ] $

  3. 让一行加上另外一行
    即 $ (a_{1,1} \pm a_{2,1})x_1 + (a_{1,2} \pm a_{2,2})x_2 + \cdots + (a_{1,n} \pm a_{2,n})x_n = c_1 + c_2 $。
    用矩阵表示，即 $ \left [ \begin{array}{ccc|c} a_{1,1} & \cdots & a_{1,n} & c_1 \\ a_{2,1} & \cdots & a_{2,n} & c_2 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m,1} & \cdots & a_{m,n} & c_m \end{array} \right ] \Leftrightarrow \left [ \begin{array}{ccc|c} a_{1,1} + a_{2,1} & \cdots & a_{1,n} + a_{2,n} & c_1 + c_2 \\ a_{2,1} & \cdots & a_{2,n} & c_2 \\ \vdots & \ddots & \vdots & \vdots \\ a_{m,1} & \cdots & a_{m,n} & c_m \end{array} \right ] $

消元

然后，解决这个多元方程的重要思想就是消元，就是要让方程中的未知数逐步减少。
其主要方法如下：首先选一个要消去的元，我们称之为主元。如果以 $ x_1 $ 为例，把第 2 ~ m 个方程中 $ x_1 $ 的系数都消成 0。要让 $ a_{i,1} $ 变成 0，我们可以让第 1 行乘上 $ \frac{a_{i,1}}{a_{1,1}} $，这样第 i 行减去第 1 行，$a_{i,1}$就会变成 0。当然，其它项也要相减。合并上述操作，可得要消除第 i 个元，对于第j 行第 k 列（i < j），我们要把 $a_{j,k}$ 减去 $ a_{i,k} \times \frac{a_{j,i}}{a_{i,i}} $；当然，$ c_j $ 也要减哦。
这样的话，我们一直这样消元下去，直到把每一行都处理完之后，就可以得到一个如下所示的矩阵，我们把它称为“上三角矩阵”。

$$ \left [ \begin{array}{ccccc|c} a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,n} & c_1 \\ 0 & a_{2,2} & a_{2,3} & \cdots & a_{2,n} & c_2 \\ 0 & 0 & a_{3,3} & \cdots & a_{2,n} & c_2 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & a_{m,n} & c_m \end{array} \right ] $$

求解

通过上面的过程，在一般情况下（特殊情况下面会提及），最后一个非零的方程就只有一个元，我们可以直接求出这个元的值。然后再代入倒数第二个方程，就可以算出倒数第二个元……这样以此类推，我们就可以求出这个方程组的解了。

高斯约旦消元

像上面，我们消元的结果是消成一个上三角矩阵，然后还要一个一个数代回去，有点小小的麻烦。其实只需要一点点的改进就可以做得更加好。
我们上面是把第 i+1 到第 m 行把系数 $ a_{j,i} $ 消成 0，然而我们可以更“暴力”一点：直接把第 1 到第 m 行把系数 $ a_{j,i} $ 消成 0！
这样的话，最后的结果就是 m 个一元一次方程，不需要回代就可以解除每一个数的解。

特殊情况

当然，也不是所有的方程都可以求出解，也不是所有的方程都有唯一解。 我们可以分成三种情况：

  1. 无解：有一个方程系数 a 都为 0，但是它的常数 c 不为 0 ，此时原方程组无解；
  2. 多解：有一个方程有两个或以上的系数不为 0，就像这个方程：$ 2x_1 + 3x_2 = 8 $，它有无数个解。以 $ x_1 $ 为例，$ x_1 $ 的解可以写成 $ x_1 = \frac{8 - 3x_2}{2} $。我们就把 $ x_1 $ 称为主元，把 $ x_2 $ 等其它的元称为自由元。这样的方程组有无数个解，一般出现在 m < n 的情况中；
  3. 唯一解：即正常情况，有且仅有一个解。

还有一个问题就是，我们刚刚选择主元的时候，第 i 行选择 $ x_i $ 作为主元。这里就会有一个问题，如果 $ a_{i,i} $ 为 0 怎么办？我们可以从 i+1 行找到第 m 行，找到一个 $ a_{i,i} $ 不为 0 的行，与第 i 行交换后继续消元即可。
从这里可以发现，理论上任意一行都是可以拿来消元的（只有 $ a_{i,i} $ 不为 0）。但是考虑到我们消元的时候我们要计算 $ \frac{a_{j,i}}{a_{i,i}} $，如果 $ a_{i,i} $ 较小，除出来的值就会很大。为了不爆炸精准起见，我们选一个绝对值最大的 $ a_{k,i} $ 代替 $ a_{i,i} $ 进行消元。如果绝对值最大的仍为 0，就说明发生了多解或无解的情况，可以跳过该次消元，到后面再来判断（有的题目可以直接当作错误或异常处理）。

代码

这里没有仔细判断多解的情况，有必要的话可以自己实现。

int n, m;
double a[MAXN][MAXN], ans[MAXN];

int Gauss() {
	for (int i=1, k; i<=n; i++) {
		k = i;
		for (int j=i+1; j<=n; j++)
			if (fabs(a[j][i]) > fabs(a[k][i]))
				k = j;
		for (int j=1; j<=n+1; j++)
			swap(a[i][j], a[k][j]);

		if (fabs(a[k][i]) <= EPS) // 无解或多解
			continue;

		for (int j=1; j<=m; j++) {
			if (j == i) continue;
			double multiple = a[j][i] / a[i][i];
			for (int k=1; k<=n+1; k++)
				a[j][k] = a[j][k] - a[i][k] * multiple;
		}
	}

	for (int i=n+1; i<=m; i++)
		if (abs(a[i][n+1]) > EPS) // 无解
			return 1;

	for (int i=1; i<=n; i++)
		ans[i] = a[i][n+1] / a[i][i];

	return 0;
}

时空复杂度

如代码所示，时间复杂度为 $ O(n^2m) $，空间复杂度为 $ O(nm) $。

其它

对于一些奇奇怪怪的题目的一些奇奇怪怪的优化或做法。

只有整数

在计算过程中，我们可以使用 int 来保证精度，当然前提是计算过程中不会爆 int。
我们在消元的时候，$ \frac{a_{j,i}}{a_{i,i}} $ 不一定是整数，我们只需要求出 $ {a_{j,i}} $ 和 $ {a_{i,i}} $ 的最小公倍数 Q，然后使 $ a_{j,k} = a_{j,k} \times \frac{Q}{a_{j,i}} - a_{i,k} \times \frac{Q}{a_{i,i}} $ 即可。其原理是先把第 i 行和第 j 行各自乘上一个整数，使 $ {a_{j,i}} $ 和 $ {a_{i,i}} $ 都变成了 Q，然后再相减。
如果答案可能为实数的话，只要再末尾 a[i][n+1] / a[i][i] 时化为实数除即可。

异或方程组

异或方程组与普通方程组差不多，但是它每个方程的形式为 $ a_1x_1 \verb!^! a_2x_2 \verb!^! \cdots \verb!^! a_nx_n = c $，其中所有系数都是 0 或 1。这个方程也可以表示为 $ a_1x_1 + a_2x_2 + \cdots + a_nx_n \equiv c \thickspace (mod \thickspace 2) $。
这时，上述的三个基本变换为：交换两行、取反一行、一个方程异或另一个方程。
在选择主元的时候，我们只要找到任意一个 $a_{k,i} = 1$ 即可。如果 $ a_{j,i} = 1 $，那么就拿第 i 行异或第 j 行即可。
最后既然所有数都是 0 或 1，所以我们可以用位运算或 bitset 来加速。加速后的时间复杂度为 $ O(\frac{n^2m}{64}) $。

CC 原创文章采用CC BY-NC-SA 4.0协议进行许可，转载请注明：“转载自：高斯消元”
CATALOG
    FRIENDS