扩展欧几里得算法总结

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)
原文地址:https://www.cnblogs.com/byronsh/p/extended-euclidean-algorithm.html