本文是我上篇博客《最大公約數》的延續，見 https://www.cnblogs.com/hanselhuang/p/gcd.html 。
擴展歐幾里得算法（簡稱擴歐）是歐幾里得算法（又稱輾轉相除法）的擴展。擴歐可以在求得 a、b的最大公約數的同時，能找到整數 x、y，使它們滿足貝祖等式：ax+by = gcd(a, b) 。如果 a 是負數，可以把問題轉化成 |a| (-x) + by = gcd(|a|, b) ，然后令 x' = (-x) 。
例如 a=6，b=9，使用歐幾里得算法我們可以得到 gcd(6, 9) = 3，而使用擴歐算法，我不僅可以計算得到最大公約數 3 ，而且可以得到方程 6x+9y = 3 的整數解 x、y 。貝祖定理告訴我們該方程必有整數解，且可以簡單證明：
方程 ax+by = c 有整數解 <=> c 是 a 和 b 的最大公約數的倍數，即 c % gcd(a, b) == 0 成立。
證明：
方程 ax+by = c 有整數解，即 x、y為整數
則 ax % gcd(a, b) == 0 成立
by % gcd(a, b) == 0 成立
=> (ax+by) % gcd(a, b) == 0
=> c % gcd(a, b) == 0
證明完畢。
前面我們一直圍繞 ax+by = gcd(a, b) 這個方程做思考，其實這個方程就是線性同余方程。線性同余方程的表述為 ax ≡ b ( mod n ) ，三橫表示恆等於，此方程有解當且僅當 b 能被 a 與 n 的最大公約數整除，此時，如果 x0 是方程的一個解，那么方程的通解可以表示為 { x0 + kn / gcd(a, n) | (k∈Z) } 。有兩個問題需要解決：為什么說貝祖等式與線性同余方程等價；方程通解的推導。
1）貝祖等式與線性同余方程等價：
ax ≡ b ( mod n ) <=> ax mod n = b mod n
b mod n = b - n ⌊b/n⌋（ 向下取整 ）
ax mod n = ax - n ⌊ax/n⌋
則 b - n ⌊b/n⌋ = ax - n ⌊ax/n⌋
令 y = ⌊b/n⌋ - ⌊ax/n⌋（ y 值域為整數 ）
=> ax + ny = b
轉化成了貝祖等式。
2）方程通解的推導：
在推導前，我們先驗證下公式，例如方程 6x ≡ 3 ( mod 9 ) 的一個特解為 5（ 把方程看成 6x mod 9 = 3 mod 9 ），則其通解可以表示為 5 + 9k / 3 = 5 + 3k ，k∈Z 。代入 k = -1，1，2 分別得到 x1 = 2，x2 = 8，x3 = 11 。可以代入原方程進行驗算，都是正確的。
我們從貝祖等式 ax+by = gcd(a, b) 入手推導其通解：
設 x0、y0 為方程的一組特解
則 a*x0 + b*y0 = gcd(a, b)
用原等式減上式得：
a (x - x0) + b (y - y0) = 0
兩邊同時除以 gcd(a, b)
(x - x0) a / gcd(a, b) = (y0 - y) b / gcd(a, b)
容易知道 a / gcd(a, b) 與 b / gcd(a, b) 互質
要使方程左右兩邊相等，則 (x - x0) 必須是 b / gcd(a, b) 的倍數
同理， (y - y0) 必須是 a / gcd(a, b) 的倍數
所以 x - x0 = b / gcd(a, b) * t，y0 - y = a / gcd(a, b) * t ，( t∈Z )
所以 x = x0 + b / gcd(a, b) * t，y = y0 - a / gcd(a, b) * t ，( t∈Z )
貝祖等式通解推導完畢。
因為貝祖等式可以等價成線性同余方程，所以可得線性同余方程 ax ≡ b ( mod n ) 的通解為 { x0 + kn / gcd(a, n) | (k∈Z) } 。
使用擴歐算法可以得到線性同余方程的特解，根據這個特解可以得到方程的通解。
擴展歐幾里得算法的C語言代碼：
#include <stdio.h> // 返回最大公約數 int e_gcd(int a, int b, int &x, int &y) { if(b == 0) { x = 1; y = 0; return a; } int x1; int d = e_gcd(b, a % b, x1, x); y = x1 - a / b * x; return d; } int main() { int a=6, b=9; int x0, y0, gcd; gcd = e_gcd(a, b, x0, y0); printf("x0:%d y0:%d gcd:%d", x0, y0, gcd); // x = x0 + (b/gcd)*t // y = y0 – (a/gcd)*t // 6x + 9y = 3 // x = -1 + 3*t // y = 1 - 2*t return 0; }
可以看到是用遞歸的方式來寫的，要理解代碼的原理，首先我們看看歐幾里得算法的遞歸寫法：
int gcd(int a, int b) { if(b==0) return a; else return gcd(b, a%b); }
歐幾里得算法的原理是 gcd(a, b) = gcd(b, a mod b) ，通過連續計算余數，直到余數等於 0 為止，最后得到的非 0 余數就是最大公約數。
那么我們如何從歐幾里得算法得到擴展歐幾里得算法呢？
我們知道方程 ax+by = gcd(a, b) 一定有整數解，我們要做的就是求出一個特解出來。
我們注意到算法的終止條件是 b = 0 ，a = gcd 。我們根據方程 ax+by = gcd(a, b) 構造這樣一個等式 gcd * x1 + 0 * y1 = gcd ，可得 x1 = 1 ，y1 = 0 為一組整數解（ y 有無窮個整數解 ），作為算法的最終狀態。
那么我們返回算法的前一個狀態，此時我們已經求得 b 和 a%b 的最大公約數為 gcd ，並且求得了一組整數解 x1 和 y1 使等式 b * x1 + (a%b) * y1 = gcd 成立。我們知道 a%b = a - b ⌊a/b⌋ ，代入前式可以得到：
b * x1 + ( a - b ⌊a/b⌋ ) * y1 = gcd
=> a * y1 + b * ( x1 - y1 ⌊a/b⌋ ) = gcd
則有 x2 = y1 ，y2 = x1 - y1 ⌊a/b⌋ 為該狀態的一組整數解。
我們通過遞歸不斷返回上一個狀態，直到回到最初狀態，就能求出最初方程的一個特解了。
由此我們可以寫出擴歐的遞歸算法：
// 返回最大公約數 int e_gcd(int a, int b, int &x, int &y) { if(b == 0) { x = 1; y = 0; return a; } int x1; int d = e_gcd(b, a % b, x1, x); y = x1 - a / b * x; return d; }
我們得到方程 ax+by = gcd(a, b) 的一組特解 x0、y0 后，就可以得到方程的通解：x = x0 + b / gcd(a, b) * t ，y = y0 - a / gcd(a, b) * t ，( t∈Z ) 。