如何利用扩展欧几里得算法求解不定方程_扩展欧几里得算法

乘法逆元

在中国剩余定理的计算里，需要求一个数字在一个模下的逆元，也就是对于给定的 a，b，找到方程

的一个整数解

。接下来我们分析一下这个方程背后隐藏着什么。

根据同余的定义，有

，也就是存在整数

使得

。移一下项，就得到了

。

这个形式恰好符合裴蜀定理

的形式，于是

，这表明

a，b互质是逆元存在的必要条件。同样可以证明： a，b互质是 a 在模 b 下存在逆元的充分条件。

---

欧几里得算法

从上面的论述可以看出，要求一个数的逆元，实际上就是找到裴蜀定理

中的 x 和 y。实际上这一步工作可以通过辗转相除法完成。

碾转相除法的依据是一个基本的事实：

。

这样，在 a > b时，令

，我们得到

。注意到

，所以这个操作可以继续进行下去，直到 b=0，此时

。

这种求最大公约数的算法叫做辗转相除法，又叫欧几里得算法。

根据上面的描述，很容易得出一个碾转相除法的递归实现：

int gcd(int a, int b) {
    if(a < b) return gcd(b, a);
    if(b == 0) return a;
    return gcd(b, a % b);
}

对这个函数进行尾递归内联优化，就得到了它的迭代版本：

int gcd(int a, int b) {
    if(a < b) return gcd(b, a);
    while(b != 0) {
        int t = a % b;
        a = b; b = t;
    }
    return a;
}

辗转相除法是计算最大公约数很快的一种方式，但是它的潜力远不止此。下面我们要用它完成更加神奇的事情。

---

扩展欧几里得算法

欧几里得算法的另一个用途就是求解

中的 x 和 y。

在欧几里得算法里，递归部分的核心是这样的：

，我们将由此得出与 x，y 相关的递归关系。

根据裴蜀定理，我们可以找到四个整数

，使得

且

。

这两个等式的右侧是相等的，于是我们得到：

，整理一下就是

。

我们希望这个等式对一切 a，b 都成立，于是

且

。

这样，要求解 x 和 y，只需要求解 x'，y'，由于 x'，y' 对应的问题规模更小，所以可以进行递归的运算。最后找一下边界条件：当 b=0 时，

对应

，

。

这种算法的思想和欧几里得算法是一致的，而且它在求出最大公约数的同时，求解了一组适于裴蜀定理的系数，所以叫做扩展欧几里得算法。

递归算法也是很好写的：

int exgcd(int a, int b, int& x, int& y) {
    if(a < b) return exgcd(b, a, y, x);
    if(b == 0) {
        x = 1; y = 0;
        return a;
    } else {
        int x1;
        int d = exgcd(b, a % b, x1, x);
        y = x1 - a / b * x;
        return d;
    }
}

x 和 y 用引用的形式传入，返回值是最大公约数。

---

扩展欧几里得算法的迭代形式

有一个不太令人高兴的事实：扩展欧几里得算法的递归是在函数中央的，不能进行尾递归内联。

有句古话说得好：一时迭代一时爽，一直迭代一直爽。

鲁迅：这锅我不背

这点困难怎么能阻挡我们寻找迭代算法的脚步！

注意到，递归的每一步都是关于 x'，y' 的线性组合，所以我们可以考虑用类似于线性递推数列的方式，用矩阵的形式改写递归式，从而利用矩阵乘法的结合律找到迭代算法。

这是我们找到的递归式：

然后变成矩阵形式：

这里

代表第一次迭代时的

。

然后就可以一直乘下去：

在每次迭代时，都可以从左向右计算一个矩阵，这样就达到了迭代的目的。

假设在一次迭代之前计算的矩阵之积为

，然后下一个矩阵

是右乘上去的：

把所有矩阵的积计算完成之后，得到一个最终的矩阵

，这个矩阵要满足

，展开整理，就得到

，

。

再看一看初始值，就是单位阵

。

所以这里可以有一个小(负)优化，在迭代时直接把矩阵写成

，这样迭代完成之后 x 和 y 的值就计算完了。

现在梳理一下思路：

1. 第一步是给矩阵赋初始值。
   
2. 然后开始迭代。每一次迭代时，计算
   。
3. 更新矩阵，把矩阵右乘迭代矩阵的结果赋值给原先的矩阵。
   
4. 进行一般欧几里得算法的迭代。
   
5. 当 b = 0 时，计算结束，此时 x 和 y 已经计算完成，并且最大公约数储存在 a 中。

给一个迭代版本的代码：

int exgcd(int a, int b, int& x, int& y) {
    if(a < b) return exgcd(b, a, y, x);
    int m = 0, n = 1;
    x = 1; y = 0;
    while(b != 0) {
        int d = a / b, t;
        t = m; m = x - d * t; x = t;
        t = n; n = y - d * t; y = t;
        t = a % b; a = b; b = t;
    }
    return a;
}

---

扩展欧几里得算法与乘法逆元

当我们用扩展欧几里得算法找到一组 x 和 y 满足

时，就可以得出

，也就是说 x 是 a 的乘法逆元。

显然，如果 x 是 a 的乘法逆元，那么所有的

也是 a 的乘法逆元。这表明，一定有 a 的一个乘法逆元在区间

内。

这是找到 n 在模 p 下的乘法逆元的一个算法，使用了上面的 exgcd：

int inv(int n, int p) {
    int x, y;
    if(exgcd(n, p, x, y) == 1) {
        x = x % p;
        return x >= 0 ? x : p + x;
    } else {
        return -1;
    }
}

如果返回 -1，表示乘法逆元不存在。

---

（封面pid：72665414）