扩展欧几里得算法总结

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 \times s_0 + b \times 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)