擴展歐幾里得算法

算法介紹

歐幾里得算法（Euclid's Algorithm）又稱輾轉相除法。古希臘數學家歐幾里得在其著作 The Elements 中最早描述了這種算法，所以該算法被命名為歐幾里得算法。算法利用公式 gcd(a,b) = gcd(b, a mod b)，求兩個非負整數 a 和 b 的最大公約數。歐幾里得算法的具體步驟如下：

1. 若 b ≠ 0，使用帶余除法，用 b 除以 a 得到余數 r；否則轉到第（3）步
2. 用 b 代替 a, 用 r 代替 b，重復第（1）步
3. a 的值就是最大公約數 d

擴展歐幾里得算法（Extend Euclid's Algorithm）是歐幾里得算法的擴展，基於以下定理：對於任意兩個整數 a, b，必存在整數 x, y，使得 xa + yb = gcd(a,b) 成立，求出整數解 x, y，可以得到 gcd(a,b)。擴展歐幾里得算法具體步驟如下：

1. 定義 x1 = 1, y1 = 0, x2 = 0, y2 = 1
2. 若 b ≠ 0，使用帶余除法，用 b 除以 a 得到余數 r；否則轉到第（4）步
3. 用 b 代替 a, 用 r 代替 b， 令 x1, y1, x2, y2 = x2, y2, x1 - q·x2, y1 - q·y2，重復第（2）步
4. a 的值就是最大公約數 d，x1 和 y1 的值就是所求 x 和 y

擴展歐幾里得原理

首先，有 \(\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} a \\ b \end{bmatrix}\)

令 \(q=[\frac{a}{b}]\)，即 q 為 \(a \div b\) 的整數部分，構建矩陣 \(\begin{bmatrix} 0 & 1 \\ 1 & -q \end{bmatrix}\)

將等式兩端同時左乘上述矩陣，則有 \(\begin{bmatrix} 0 & 1 \\ 1 & -q \end{bmatrix} \cdot \begin{bmatrix} a\\b \end{bmatrix} = \begin{bmatrix} b \\ a-q \cdot b \end{bmatrix} = \begin{bmatrix} b \\ a \; mod \; b \end{bmatrix}\)

再令 \(q'=[\frac{b}{a \; mod \; b}]\)，重復上述步驟

最終有 \(\begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \end{bmatrix} \cdot \begin{bmatrix} a\\b \end{bmatrix} = \begin{bmatrix} gcd(a,b) \\ 0 \end{bmatrix}\)

即求出 \(x_1 \cdot a + y_1 \cdot b = gcd(a,b)\)

算法流程圖

具體代碼(python)

def ex_gcd(a: int, b: int) -> (int, int, int):
    """擴展歐幾里得算法

    :return: gcd x y
    """
    # 0與0沒有最大公約數，因為任何非零數都能整除 0
    if a == 0 and b == 0:
        return None
    else:
        x1, y1, x2, y2 = 1, 0, 0, 1  # 初始化x1,y1,x2,y2
        while b:
            q, r = divmod(a, b)
            a, b = b, r  # gcd(a,b)=gcd(b,a%b)
            # 下面操作模仿矩陣乘法，即
            # [[0,  1]    [[x1, y1]
            #  [1, -q]] x  [x2, y2]]
            x1, y1, x2, y2 = x2, y2, x1 - q * x2, y1 - q * y2
        # 返回一個三元組，依次是(gcd,x,y)，使得 xa+yb=gcd
        # 若傳入的 a,b 是負數，那么最后結果有可能為負，需要取負
        return (a, x1, y1) if a > 0 else (-a, -x1, -y1)

測試樣例和結果

# 樣例一
a = 31
b = -13
# 結果 gcd = xa + yb
gcd = 1
x = -5
y = -32

# 樣例二
a = 2461502723515673086658704256944912426065172925575
b = 1720876577542770214811199308823476528929542231719
# 結果 gcd = xa + yb
gcd = 1
x = -279503883787274514253902210297369853780943628056
y = 399795999407450264570090854758429065117844005679

# 樣例三
a = 13709616469144948883512229123502305176385931810284088906755090238431898972708904439178898468021710798401875986657125211084472621499595371254346390738382042
b = 19235039994987625167590963480899777255933775238312044097122773255647530276806317636026727679800825370459321617724871515442147432420951257037823141069640181
# 結果 gcd = xa + yb
gcd = 1
x = -1076045803680437575165069317517720056816012968550552297497847522201869587140865136718292159497811949579872126471694232065112371318866682108000433819750738
y = 766942791672689226450200919445856853258490754050379704483194087100115700225755919066211027540484178869823084929299459228617007521725575705884841982026337

用途

在密碼學中，互素是一個非常重要的概念，判斷兩個數互素就是計算其最大公約數是否為 1

同時，密碼學中的逆元也是非常重要的概念，整數 e 在模數 p 下的逆元 d 需要滿足 e·d mod p = 1。也就是 e·d + k·p = 1。使用擴展歐幾里得算法計算 ex_gcd(e,p) 即可算出：如果結果的 gcd = 1，那么 x 就是 e 在模 p 下的逆元；如果結果的 gcd ≠ 1，那么表明 e 在模 p 下沒有逆元

利用歐幾里得求逆元

def inverse(n: int, p: int):
    """
    :return: n 在 mod p 下的逆元
    """
    gcd_, iv, _ = ex_gcd(n, p)
    if gcd_ != 1:
        raise ValueError(str(n) + " has no inverse in mod " + str(p))
    return iv % p

參考資料：《密碼學實驗教程》