扩展欧几里得算法

本文是我上篇博客《最大公约数》的延续，见 。

扩展欧几里得算法（简称扩欧）是欧几里得算法（又称辗转相除法）的扩展。扩欧可以在求得 a、b的最大公约数的同时，能找到整数 x、y，使它们满足贝祖等式：ax+by = gcd(a, b) 。如果 a 是负数，可以把问题转化成 |a| (-x) + by = gcd(|a|, b) ，然后令 x' = (-x) 。

例如 a=6，b=9，使用欧几里得算法我们可以得到 gcd(6, 9) = 3，而使用扩欧算法，我不仅可以计算得到最大公约数 3 ，而且可以得到方程 6x+9y = 3 的整数解 x、y 。贝祖定理告诉我们该方程必有整数解，且可以简单证明：

方程 ax+by = c 有整数解  <=>  c 是 a 和 b 的最大公约数的倍数，即 c % gcd(a, b) == 0 成立。

证明：

方程 ax+by = c 有整数解，即 x、y为整数

则 ax % gcd(a, b) == 0 成立

by % gcd(a, b) == 0 成立

=> (ax+by) % gcd(a, b) == 0

=> c % gcd(a, b) == 0

证明完毕。

前面我们一直围绕 ax+by = gcd(a, b) 这个方程做思考，其实这个方程就是线性同余方程。线性同余方程的表述为 ax ≡ b ( mod n ) ，三横表示恒等于，此方程有解当且仅当 b 能被 a 与 n 的最大公约数整除，此时，如果 x0 是方程的一个解，那么方程的通解可以表示为 { x0 + kn / gcd(a, n) | (k∈Z) } 。有两个问题需要解决：为什么说贝祖等式与线性同余方程等价；方程通解的推导。

1）贝祖等式与线性同余方程等价：

ax ≡ b ( mod n )  <=>  ax mod n = b mod n

b mod n = b - n ?b/n?（ 向下取整 ）

ax mod n = ax - n ?ax/n?

则 b - n ?b/n? = ax - n ?ax/n?

令 y = ?b/n? - ?ax/n?（ y 值域为整数 ）

=> ax + ny = b

转化成了贝祖等式。

2）方程通解的推导：

在推导前，我们先验证下公式，例如方程 6x ≡ 3 ( mod 9 ) 的一个特解为 5（ 把方程看成 6x mod 9 = 3 mod 9 ），则其通解可以表示为 5 + 9k / 3 = 5 + 3k ，k∈Z 。代入 k = -1，1，2 分别得到 x1 = 2，x2 = 8，x3 = 11 。可以代入原方程进行验算，都是正确的。

我们从贝祖等式 ax+by = gcd(a, b) 入手推导其通解：

设 x0、y0 为方程的一组特解

则 a*x0 + b*y0 = gcd(a, b)

用原等式减上式得：

a (x - x0) + b (y - y0) = 0

两边同时除以 gcd(a, b)

(x - x0) a / gcd(a, b) = (y0 - y) b / gcd(a, b)

容易知道 a / gcd(a, b) 与 b / gcd(a, b) 互质

要使方程左右两边相等，则 (x - x0) 必须是 b / gcd(a, b) 的倍数

同理， (y - y0) 必须是 a / gcd(a, b) 的倍数

所以 x - x0 = b / gcd(a, b) * t，y0 - y = a / gcd(a, b) * t ，( t∈Z )

所以 x = x0 + b / gcd(a, b) * t，y = y0 - a / gcd(a, b) * t ，( t∈Z )

贝祖等式通解推导完毕。

因为贝祖等式可以等价成线性同余方程，所以可得线性同余方程 ax ≡ b ( mod n ) 的通解为 { x0 + kn / gcd(a, n) | (k∈Z) } 。

使用扩欧算法可以得到线性同余方程的特解，根据这个特解可以得到方程的通解。

扩展欧几里得算法的C语言代码：

#include 

// 返回最大公约数 
int e_gcd(int a, int b, int &x, int &y) {
    if(b == 0) {
        x = 1; y = 0;
        return a;
    }
    int x1;
    int d = e_gcd(b, a % b, x1, x);
    y = x1 - a / b * x;
    return d;
}

int main() {
    int a=6, b=9;
    int x0, y0, gcd;
    gcd = e_gcd(a, b, x0, y0);
    
    printf("x0:%d y0:%d gcd:%d", x0, y0, gcd);
    // x = x0 + (b/gcd)*t
    // y = y0 – (a/gcd)*t
    
    // 6x + 9y = 3
    // x = -1 + 3*t
    // y =  1 - 2*t
    
    return 0;
}

可以看到是用递归的方式来写的，要理解代码的原理，首先我们看看欧几里得算法的递归写法：

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

欧几里得算法的原理是 gcd(a, b) = gcd(b, a mod b) ，通过连续计算余数，直到余数等于 0 为止，最后得到的非 0 余数就是最大公约数。

那么我们如何从欧几里得算法得到扩展欧几里得算法呢？

我们知道方程 ax+by = gcd(a, b) 一定有整数解，我们要做的就是求出一个特解出来。

我们注意到算法的终止条件是 b = 0 ，a = gcd 。我们根据方程 ax+by = gcd(a, b) 构造这样一个等式 gcd * x1 + 0 * y1 = gcd ，可得 x1 = 1 ，y1 = 0 为一组整数解（ y 有无穷个整数解 ），作为算法的最终状态。

那么我们返回算法的前一个状态，此时我们已经求得 b 和 a%b 的最大公约数为 gcd ，并且求得了一组整数解 x1 和 y1 使等式 b * x1 + (a%b) * y1 = gcd 成立。我们知道 a%b = a - b ?a/b? ，代入前式可以得到：

b * x1 + ( a - b ?a/b? ) * y1 = gcd

=> a * y1 + b * ( x1 - y1 ?a/b? ) = gcd

则有 x2 = y1 ，y2 = x1 - y1 ?a/b? 为该状态的一组整数解。

我们通过递归不断返回上一个状态，直到回到最初状态，就能求出最初方程的一个特解了。

由此我们可以写出扩欧的递归算法：

// 返回最大公约数 
int e_gcd(int a, int b, int &x, int &y) {
    if(b == 0) {
        x = 1; y = 0;
        return a;
    }
    int x1;
    int d = e_gcd(b, a % b, x1, x);
    y = x1 - a / b * x;
    return d;
}

我们得到方程 ax+by = gcd(a, b) 的一组特解 x0、y0 后，就可以得到方程的通解：x = x0 + b / gcd(a, b) * t ，y = y0 - a / gcd(a, b) * t ，( t∈Z ) 。