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)