被迫学习之拓展欧几里得算法 & 高斯消元

    ICPC 数学   algorithm math  Comments

为了在我自己什么都不会的背景下给 OI 队学弟学妹写材料，“被迫”去翻出书本和课件，重新啃了数论的内容学这两个东西……（发出退役咸鱼的声音

拓展欧几里得算法 (exgcd)

对于不全为 \(0\) 的非负整数 \(a, b\)，\(gcd(a, b)\) 表示 \(a, b\) 的最大公约数，则必然存在整数对 \(x, y\) ，使得 \(gcd(a, b)= ax + by\)。

比如 \(a = 3, b = 5\) 的时候，\(gcd(a, b) = 1\), 那么有 \(x, y\) 的一组解为 \(x = 2, y = -1\) 使得 \(ax + by = 6 - 5 = 1 = gcd(a, b)\).

拓展欧几里得算法就是用来求解 x, y 的。那么它是如何做到的？现在我们已知一组非负整数 \(a, b (a \geq b)\):

  1. 显然当 \(b = 0\) 的时候，\(gcd(a, b) = a\)，此时容易知道当 \(x = 1, y = 0\) 时为一组解。

  2. 否则，当 \(a > b > 0\) 的时候，设 \(ax\_{1} + by\_{1} = gcd(a, b)\)

  3. 那么有 \(bx_{2} + (a \mod b) \cdot y_{2} = gcd(b, a \mod b)\)

  4. 根据欧几里得原理，容易知道 \(gcd(a, b) = gcd(b, a \mod b)\)；

  5. 从而 \(ax_{1} + by_{1} = bx_{2} + (a \mod b) \cdot y_{2}\),

  6. 有 \(ax_{1} + by_{1} = bx_{2} + (a - [a \div b] \cdot b) \cdot y_{2} = ay_{2} + bx_{2} - [a \div b] \times by_{2}\), 其中 \([]\) 表示向下取整符号。

  7. 消去系数，最后求得 \(x_{1} = y_{2}\), \(y_{1} = x_{2} - [a \div b]y_{2}\).

我们得到了 \(x1, y1\) 和 \(x2, y2\) 的关系，也就是 \(x1, y1\) 可以由 \(x2, y2\) 表示。其中 \(x2, y2\) 是 \(bx + (a \mod b)y = gcd(b, a \mod b)\) 的一组解。但我们两组解都不知道，此时我们可以倒回去推。

可以证明在求解 \(gcd(a, b)\) 的递归过程中最后 \(b\) 一定等于 0, 此时 \(x = 1, y = 0\), 然后我们可以据此倒着推回去。

如求解 \(exgcd(6, 4)\), 即求使 \(6x + 4y = gcd(6, 4) = 2\) 成立的 \(x, y\)，

由 \(gcd(6, 4) = gcd(4, 2)\) 得当 \(a’ = 4, b’ = 2\) 时，\(a’x’ + b’y’ = gcd(4, 2) = 2\)，根据上面的推断式不难得到 \(x = y’, y = x’ - [6 \div 4]y’\)

由 \(gcd(4, 2) = gcd(2, 0)\) 得当 \(a” = 2, b” = 0\) 时，\(a”x” + b”y” = gcd(2, 0) = 2\)，此时由于 \(b” = 0\)，可得 \(x” = 1, y” = 0\).

再由上面的推断式 \(x’ = y”, y’ = x” - [4 \div 2]y”\)，代入 \(x”, y”\) 得到 \(x’ = 0, y’ = 1 - 2 \times 0 = 1\)

最后把 \(x’ = 0, y’ = 1\) 代入 \(x = y’, y = x’ - [6 \div 4]y’\) 得到 \(x = 1, y = 0 - 1 \times 1 = -1\).

将 \(x = 1, y = -1\) 代入 \(6x + 4y = gcd(6, 4) = 2\)，此组解有效。

可以看到我们的求解过程是层层递归的，我们利用欧几里得定理，不断让 \(a’ = b, b’ = a \mod b\)，直到 \(b = 0\) 得到 \(x = 1, y = 0\)，再用这个 \(x, y\) 层层回溯根据推导式，计算出上一层的 \(x, y\) 直到得到结果。明白了这样的思想之后，我们就可以像下面这样实现拓展欧几里得算法：

#include <iostream>
int x, y, a, b;

void exgcd(int a, int b) {
    if (b == 0) {
        x = 1;
        y = 0;
        return;
    } else {
        // x = y', y = x' - [a / b] * y'
        exgcd(b, a % b);
        int tmp = y;
        y = x - (a / b) * y;
        x = tmp;
    }
}

int main() {
    a = 6, b = 4;
    exgcd(6, 4);
    std::cout << x << " " << y;
    return 0;
}

输出：

1 -1

有时候，exgcd 可以简化如下, 其中 d = gcd(a, b):

void exgcd (int a, int b, int &d, int &x, int &y) {
    if (!b) {
        d = a, x = 1, y = 0;
    } else {
        exgcd(b, a % b, d, y, x);       // 注意 x,y 的顺序
        y -= x * (a / b);               // 原理与上文相同
    }
}

显然，我们求出的这个 x, y 可能是不唯一的。出现多解原因在于我们规定当 \(b = 0\) 的时候 \(y = 0\)，实际上此时 \(y\) 可以取任何整数。

那么这个算法有什么卵用？它可以用来解二元一次的不定方程。对于不定方程 \(ax + by = c\)，设 \(gcd(a, b) = d\)，如果 \(ax + by = c\) 有解，那么 \(d\) 是 \(c\) 的约数，即 \(d|c, c = n \cdot d\) （充分必要条件）.

那么我们只要先求出 \(ax’ + by’ = gcd(a, b) = d\) 的 \(x’, y’\)，分别对 \(x”, y\) 乘 \(\frac{c}{d}\) 倍即可。或者，我们可以让 \(b = 0\) 时的 exgcd 结果为 \(\frac{c}/{d}\) ，这样一遍 exgcd 跑完之后的 x, y 就是方程的解了。

int a, b, x, y, c;      // ax+by = c
void exgcd(int a, int b) {
    if (b == 0) {
        x = c / a, y = 0;        // 当 b' = 0时，a' = gcd(a, b)
    } else {
        exgcd(b, a % b);
        int tmp = y;
        y = x - (a / b) * y;
        x = tmp;
    }
}

因为 exgcd 所求得的 x, y 可以不唯一，而有些题目要求我们输出符合某些条件的解，所以我们还可以对求得的最终结果加以限定条件。例如题目要求 x < 某一个特定的数 p 时：

y = x - (a / b) * y;
x = tmp;
// ------
if (x > p) {
    int t = x / p;      // 求 x 是 p 的多少倍
    x -= t * p;         // 从 x 中减去 tp
    y += t * p * a / b;     // 根据 ax+by=c, a(x-tp) + b(y+k) = c, 得 -atp+bk = 0, k = tpa/b
}

或者，求最小正整数解，参考 https://blog.csdn.net/qq_20200047/article/details/71159677：

if (c % d == 0) {
    x *= c / d;
    int t = b / d;
    x = (x % t + t) % t;        // x
    y = (c - (a * x)) / b;      // y
}

高斯消元

高斯消元法是用来解线性规划（线性方程组）的算法。

高斯消元的思想和我们解多元方程其实是差不多的。我们来看一个例子，用数学的方法解下列方程组：

\(\begin{cases}\ 2x + y + z = 1 (1)\ 6x + 2y + z = -1 (2)\ \ -2x + 2y + z = 7 (3)\ \end{cases}\)

  1. 消去 \(x\):

\((1) \times (-3) + (2)\) 得到 \(0x - y - 2z = -4\);

\((1) + (3)\) 得到 \(0x + 3y + 2z = 8\)，有：

\(\begin{cases}\ 2x + y + z = 1\ 0x - y - 2z = -4\ 0x + 3y + 2z = 8\end{cases}\)

  1. 消去 \(y\):

\((2) \times 3 + (3)\) 得到 \(0x + 0y - 4z = -4\)，有：

\(\begin{cases}\ 2x + y + z = 1\ 0x - y - 2z = -4\ 0x + 0y - 4z = -4\end{cases}\)

  1. 从而我们得到 \(z = 1\)，将 \(z = 1\) 回带到 (2) 得到 \(y = 2\).

  2. 再把 \(z = 1, y = 2\) 回带到 (1) 得到 \(x = -1\).

在我们具体的求解中，我们是把方程组的系数和值存储在数组中的：

index x (0) y (1) z (2) val (3)
0 2 1 1 1
1 6 2 1 -1
2 -2 2 1 7

因此，我们希望这个方程组最终变为如下的矩阵：

index x (0) y (1) z (2) val (3)
0 1 0 0 x
1 0 1 0 y
2 0 0 1 z

使用高斯消元法, 是以列为单位消元的。最终，第 i 行的 val 表示的就是第 i 个未知数的解。所以第一步我们需要枚举每一个式子 (i, 竖直行方向);

为了减少浮点误差，我们规定第 i 个式子往下的式子中，第 i 个未知数中系数最大的一个式子应当排在当前被减的这一行,假设其在第 k 行，那么我们就需要把 a[i][*] 与 a[k][*] 交换。以 i + 1 为起点枚举找到符合条件的 k，交换即可。例如当 i = 0 时执行上述交换，得到的新矩阵为：

x (0) y (1) z (2) val (3)
0 6 2 1 -1
1 2 1 1 1
2 -2 2 1 7

然后将正在处理的方程式化简，让正被处理的系数化 1:

x (0) y (1) z (2) val (3)
0 1 1⁄3 1⁄6 -1⁄6
1 2 1 1 1
2 -2 2 1 7

再使用加减法，将剩下式子的第一个系数化 0。如化简 1 式时，因为 x 的系数为 2，那么就让 \(a[1][k] -= a[0][k] * 2\); 化简 3 式时，因为 x 的系数为 -2，那么就让 \(a[2][k] -= a[0][k] * (-2)\)，结果如下：

x (0) y (1) z (2) val (3)
0 1 1⁄3 1⁄6 -1⁄6
1 0 1⁄3 2⁄3 4⁄3
2 0 8⁄3 4⁄3 20⁄3

同理，用这样的方法继续处理第 1,2 行，最后即可得到结果。例如处理第 1 行，把 \(y\) 项系数最大的式子提前：

x (0) y (1) z (2) val (3)
0 1 1⁄3 1⁄6 -1⁄6
1 0 8⁄3 4⁄3 20⁄3
2 0 1⁄3 2⁄3 4⁄3

让第二个未知数 \(y\) 的系数为 1：

x (0) y (1) z (2) val (3)
0 1 1⁄3 1⁄6 -1⁄6
1 0 1 1⁄2 5⁄2
2 0 1⁄3 2⁄3 4⁄3

然后使用加减法将剩余的式子化零：

x (0) y (1) z (2) val (3)
0 1 0 0 -1
1 0 1 1⁄2 5⁄2
2 0 0 1⁄2 1⁄2

最后按照上面的步骤处理第 2 行的式子。最终我们得到的矩阵如下：

x (0) y (1) z (2) val (3)
0 1 0 0 -1
1 0 1 0 2
2 0 0 1 1

要注意一下，高斯消元存在无解/多解的情况：如果有一个未知数在所有式子中，它的系数都是 0，那么该未知数是无法求解的，那么方程组也就不存在唯一解了。

luogu P3389 模板：高斯消元法

#include <cstdio>
#include <cmath> 
#include <algorithm>

const int MAXN = 205;
const double eps = 1e-8;        // 浮点误差

int n;
double del;
double a[MAXN][MAXN];

bool gauss() {
    // 对每个式子进行枚举 
    for (int i = 0; i < n; i++) {
        int k = i;
        // 从 i 的下一个等式开始枚举 
        for (int j = i + 1; j < n; j++) {
            // 比较第 k 个式子和第 j 个式子中未知数 i 系数的大小, fabs() 表示对浮点数取绝对值 
            if (fabs(a[j][i]) > fabs(a[k][i])) {
                k = j;      // 找到第 i 个未知数系数最大的 
            }
        }
        
        del = a[k][i];
        
        // 如果有一个未知数在所有式子中的系数都为 0（绝对值最大值小于浮点误差），则无解 
        if (fabs(del) < eps) {
            return false;
        }
        
        // 对于当前式子，与第 k 个式子交换所有项的系数，即当前式子实为未知数 i 系数最大的第 k 个式子 
        for (int j = i; j <= n; j++) {
            std::swap(a[i][j], a[k][j]);
        }
        
        // 然后对当前式子的所有项（包括值）消去第 k 个式子中第 i 个未知数的系数 
        for (int j = i; j <= n; j++) {
            a[i][j] /= del;
        }
        
        // 处理其它的式子，对各项加上/减去相同系数 
        for (k = 0; k < n; k++) {
            // 只对非当前式子操作 
            if (k != i) {
                del = a[k][i];
                for (int j = i; j <= n; j++) {
                    a[k][j] -= a[i][j] * del;
                }
            }
        }
    }
    return true;
}

int main() {
    // 未知数的数目 
    scanf("%d", &n);
    
    // n 个未知数，至少需要有 n 个等式才能求解，否则无解 
    for (int i = 0; i < n; i++) { 
        // 注意第二个循环体的条件是 j <= n，不要忘记了除了系数外还要输入一个多项式的值 
        for (int j = 0; j <= n; j++) {
            scanf("%lf", &a[i][j]);
        }
    }
    // a[i][j] 表示在第 i 个等式中第 j 个未知数的系数 
    
    // 高斯消元求解，如果该方程组有解，返回 true 
    bool flag = gauss(); 
    
    if (!flag) {
        // 无解
        printf("No Solution");
    } else {
        for (int i = 0; i < n; i++) {
            // 输出第 i 个式子的值，即为第 i 个未知数的答案
            printf("%.2lf\n", a[i][n]);
        }
    }
    
    return 0;
}

Yume Maruyama

まんまるお山に彩りを✨~