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)