扩展欧几里得算法总结

Extended Euclidean algorithm 扩展欧几里得算法

扩展欧几里得算法是欧几里得算法的扩展。欧几里得算法又称辗转相除法，看名字就知道是大数学家欧几里得发现的，它可以用于计算两个数 (a) 和 (b) 的最大公约数（greatest common divisor, 简称 (gcd)）。本文把 (a) 和 (b) 的最大公约数记作 (gcd(a, b))。

欧几里得扩展算法不仅计算两个整数 a 和 b 的最大公约数 (gcd(a, b))（记作(d)），也计算满足下面关系的 (x) 和 (y)
(a×x+b×y=d= gcd⁡(a, b)) --(1)

给定两个整数 a 和 b，欧几里得算法只计算(1)式中的 (d)，欧几里得扩展算法要计算 (d), (x) 和 (y)。

计算方法

下面是一种计算方法。

欧几里得扩展算法要计算 (x), (y) 和 (d)，它们分别是下面算法里面 (s), (t) 和 (r) 变量的最终值，算法里 (s), (t) 和 (r) 关系为
(r=a×s+b×t)

为什么构造这么一个式子？ 你不觉得和 式子1很类似吗？这就是数学家的灵光一现。我们能证明它正确，但却不能证明他怎么会想到这个式子。还好我们只需要记住这个式子和笃信它正确即可。

下面就演示如何计算得到 (d), (x) 和 (y)。我自己觉得是很直观。如果你觉得很难理解，可以下面留言交流。

按照下面的序列计算
序列 0：(r_0=a=a×s_0+b×t_0)，其中 (s_0=1), (t_0=0)
序列 1：(r_1=b=a×s_1+b×t_1)，其中 (s_1=0), (t_1=1)
序列 2：(r_2=r_0 −q_1× r_1)， 其中 (q_1)= (r_0)除 (r_1) 的商（quotients)。所以 (s_2=s_0−q_1×s_1), (t_2=t_0 −q_1×t_1)
…
序列 i+1：(r_{i+1}=r_{i−1} −q_i×r_i)， 其中 (q_i)= (r_{i−1})除 (r_i) 的商。所以 (s_{i+1}=s_{i−1}−q_i×s_i), (t_{i+1}=t_i −q_i×t_i)
...

按照上面的序列计算 (r_i), (s_i) 和 (t_i)。如果你看过欧几里得算法的计算序列，一定觉得上面的计算序列很亲切，几乎一样的套路处理 (r_i), (s_i) 和 (t_i)。

理解的关键点： 上面的计算序列中，(r_i), (s_i) 和 (t_i) 之间一直满足关系 (r_i=a×s_i+b×t_i)。比如在 (i=0) 时， (r_0 = a imes s_0 + b imes t_0), 因为其中 (r_0 = a), (s_0 = 1) 和 (t_0 = 0)。如果(r_i) 是 (gcd(a, b))，那自然就得到期望的 (x) 和 (y)。

如何得到计算结果

假设在i+1次计算时 (r_{i+1}) 为 0 ，那么 (r_i) 就是 (gcd⁡(a, b))，这是因为 (r_{i-1}) 除 (r_i) 余 0， 那么 (r_i) 就是 (r_{i-1}) 和 (r_i) 的最大公约数，也就是 (a) 和 (b) 的最大公约数。具体证明网上搜欧几里得算法或辗转相除法。又因为 (r_i), (s_i) 和 (t_i) 之间一直满足关系 (r_i=a×s_i+b×t_i)，所以 (s_i) 和 (t_i) 就满足 (r_i=gcd⁡(a, b)= a×s_i+b×t_i)，即 (s_i) 和 (t_i) 就是我们要找的 (x) 和 (y)。

Python 实现

下面是扩展欧几里德算法的python函数。非常简单，可能看代码比看上面解释更容易懂。

def exgcd(a, b):
    """Calculate x, y and d of below expression based on input a and b:
    a*x + b*y = gcd(a, b)  -- (1)
    Where gcd(a, b) is the greatest common divisor（最大公约数） of a and b.
    Let's call gcd(a, b) as *g* from now on.
    
    Return:
    (x, y, g) - one tuple contains solution x, y and g of formula (1)
    
    Exception:
    Input a and b are not integers.
    """
    if False==isinstance(a, int) or False==isinstance(b, int):
        raise Exception("a, b must be integers!")
    r0, r1 = a, b
    s0, s1 = 1, 0
    t0, t1 = 0, 1
    while True:
        q, r = divmod(r0, r1)
        if 0==r:
            break # r1, s1 and t1 stores g, x, y
        r0, r1 = r1, r
        s0, s1 =s1, s0 - q*s1
        t0, t1 =t1, t0 - q*t1
     
    return (s1, t1, r1)