利用基本列運算實現擴展歐幾里得演算法

本文的閱讀等級:初級

美國計算機科學家高德納 (Donald Knuth) 說[1]

我們可以稱它是所有演算法的始祖,因為它是至今尚存最古老的不平凡演算法。

高德納所說的「所有演算法的始祖」即為歐幾里得演算法 (Euclidean algorithm),又稱為輾轉相除法,它是求最大公約數的一種算法。最大公約數是指能夠同時整除兩個整數的最大正整數。

 
a, b 為二整數,它們的最大公約數記為 \gcd(a,b)。因為 \gcd(\vert a\vert, \vert b\vert)=\gcd(a,b),為了方便,我們可假設 a\ge b>0。歐幾里得演算法的基本原理在於這個性質:若 a=qb+r,其中 q 是商,r 是餘數且 0\le r<b,則 \gcd(a,b)=\gcd(b,r)。歐幾里得演算法重複運用此性質輾轉相除以得到最大公約數,運算過程可表示為下列方程組:

\begin{aligned} a&=q_1b+r_1~~~~~0<r_1<b\\ b&=q_2r_1+r_2~~~~0<r_2<r_1\\ r_1&=q_3r_2+r_3~~~~0<r_3<r_2\\ &~\vdots\\ r_{n-2}&=q_nr_{n-1}+r_n~~~0<r_n<r_{n-1}\\ r_{n-1}&=q_{n+1}r_n+0 \end{aligned}

因為餘數不為負,且每一步驟的餘數必小於前一步驟的餘數,故存在 r_{n+1}=0,最後一個非零餘數 r_n 即為 \gcd(a,b)。例如,欲求 \gcd(5892,1902),運用輾轉相除可得

\begin{aligned} 5892&=3\cdot 1902+186\\ 1902&=10\cdot 186+42\\ 186&=4\cdot 42+18\\ 42&=2\cdot 18+6\\ 18&=3\cdot 6+0 \end{aligned}

所以,

\begin{aligned} \gcd(5892,1902)&=\gcd(1902,186)=\gcd(186,42)=\gcd(42,18)\\ &=\gcd(18,6)=\gcd(6,0)=6\end{aligned}

 
貝祖等式 (Bézout’s identity) 說 [2]:給定兩個不全為零的整數 ab,存在整數 xy,稱為貝祖係數,使得

ax+by=\gcd(a,b)

上例中,為了將 6 寫成 58921902 的線性組合,我們從先前所示倒數第二個方程式開始,逐次迭代消去餘數 18, 42186,過程如下:

\begin{aligned} 6&=42-2\cdot 18\\ &=42-2(186-4\cdot 42)\\ &=(-2)186+9\cdot 42\\ &=(-2)186+9(1902-10\cdot 186)\\ &=9\cdot 1902-92\cdot 186\\ &=9\cdot 1902-92(5892-3\cdot 1902)\\ &=(-92)5892+285\cdot 1902, \end{aligned}

可知 x=-92y=285 滿足貝祖等式 5892x+1902y=\gcd(5892,1902)=6。注意,貝祖係數並不唯一,譬如,

\begin{aligned} 6&=(-92+1902)5892+(285-5892)1902\\ &=1810\cdot 5892-5607\cdot 1902. \end{aligned}

 
事實上,欲解出貝祖等式並不需要如上例先以歐幾里得演算法得出 \gcd(a,b),再用反向迭代求出貝祖係數 x, y。貝祖等式可由擴展歐幾里得演算法 (extended Euclidean algorithm) 算出[2],即在原有的歐幾里得演算法上增加兩個遞歸等式:對於 k=1,2,\ldots

\begin{aligned} r_k&=r_{k-2}-q_kr_{k-1}\\ s_k&=s_{k-2}-q_ks_{k-1}\\ t_k&=t_{k-2}-q_kt_{k-1} \end{aligned}

初始值設為 r_{-1}=ar_0=bs_{-1}=1s_0=0t_{-1}=0t_{0}=1。當 r_{n+1}=0,演算法終止,最大公約數是 r_n,貝祖係數分別是 x=s_ny=t_n

 
接著我們解說如何利用矩陣方法實現擴展歐幾里得演算法並驗證其正確性。擴展歐幾里得演算法的遞歸算式可用矩陣的列取代運算表示如下:

\begin{bmatrix} s_{-1}&t_{-1}&r_{-1}\\ s_{0}&t_{0}&r_{0} \end{bmatrix}\to\begin{bmatrix} s_{1}&t_{1}&r_{1}\\ s_{0}&t_{0}&r_{0} \end{bmatrix}\to\begin{bmatrix} s_{1}&t_{1}&r_{1}\\ s_{2}&t_{2}&r_{2} \end{bmatrix}\to\cdots

注意,每一步驟總是取代指標較小的那一列,因此交互改變列1和列2。遞歸等式定義的列取代運算可表示成矩陣乘法

\left[\!\!\begin{array}{cr} 1&-q_k\\ 0&1 \end{array}\!\!\right]\begin{bmatrix} s_{k-2}&t_{k-2}&r_{k-2}\\ s_{k-1}&t_{k-1}&r_{k-1} \end{bmatrix}=\begin{bmatrix} s_{k}&t_{k}&r_{k}\\ s_{k-1}&t_{k-1}&r_{k-1} \end{bmatrix}

\left[\!\!\begin{array}{rc} 1&0\\ -q_k&1 \end{array}\!\!\right]\begin{bmatrix} s_{k-1}&t_{k-1}&r_{k-1}\\ s_{k-2}&t_{k-2}&r_{k-2} \end{bmatrix}=\begin{bmatrix} s_{k-1}&t_{k-1}&r_{k-1}\\ s_{k}&t_{k}&r_{k} \end{bmatrix}

上式中,左乘矩陣即為對應基本列運算的基本矩陣 (見“特殊矩陣 (10):基本矩陣”)

E_{12}(-c)=\left[\!\!\begin{array}{cr} 1&-c\\ 0&1 \end{array}\!\!\right],~~E_{21}(-c)=\left[\!\!\begin{array}{rc} 1&0\\ -c&1 \end{array}\!\!\right]

基本矩陣 E_{ij}(-c) 的作用是列 i 減去列 j 通乘一數 c,即 \mathrm{row}_i\leftarrow\mathrm{row}_i-c\cdot\mathrm{row}_j

 
下面展示以基本列運算求解貝祖等式 5892x+1092y=\gcd(5892,1902) 的詳細計算步驟:

\begin{aligned} \begin{bmatrix} 1&0&5892\\ 0&1&1902 \end{bmatrix}&\xrightarrow[]{E_{12}(-3)}\left[\!\!\begin{array}{crr} 1&-3&186\\ 0&1&1902 \end{array}\!\!\right]\xrightarrow[]{E_{21}(-10)}\left[\!\!\begin{array}{rrr} 1&-3&186\\ -10&31&42 \end{array}\!\!\right]\\ &\xrightarrow[]{E_{12}(-4)}\left[\!\!\begin{array}{rrr} 41&-127&18\\ -10&31&42 \end{array}\!\!\right]\xrightarrow[]{E_{21}(-2)}\left[\!\!\begin{array}{rrr} 41&-127&18\\ -92&285&6 \end{array}\!\!\right]\\ &\xrightarrow[]{E_{12}(-3)}\left[\!\!\begin{array}{rrc} 317&-982&0\\ -92&285&6 \end{array}\!\!\right] \end{aligned}

從最後得到的矩陣可知 (-92)5892+285\cdot 1902=6。為甚麼如此?設想第一個增廣矩陣對應線性方程

\begin{bmatrix} 1&0\\ 0&1 \end{bmatrix} \begin{bmatrix} u\\ v \end{bmatrix}=\begin{bmatrix} 5892\\ 1902 \end{bmatrix}

因為基本列運算不改變線性方程的解,故知最後一個增廣矩陣所對應的線性方程

\left[\!\!\begin{array}{rr} 317&-982\\ -92&285 \end{array}\!\!\right] \begin{bmatrix} u\\ v \end{bmatrix}=\begin{bmatrix} 0\\ 6 \end{bmatrix}

也有相同解 u=5892v=1902,即

\begin{aligned} 317\cdot 5892-982\cdot 1902&=0\\ (-92)5892+285\cdot 1902&=6. \end{aligned}

 
法國數學家拉米 (Gabriel Lamé) 於公元1844年證明歐幾里得演算法需要的計算步驟不超過較小數 (即 b) 的位數的5倍[2]。上例中,b=1902,不論 a\ge b 為何數,歐幾里得演算法至多執行 5\times 4=20 個步驟 (對於 a=5892,我們實際執行了5個步驟)。歐幾里得演算法對於處理大數字問題具有很高的運算效率,這說明了何以高德納說它是流傳至今最古老的非凡演算法。

 
參考來源:
[1] Donald E. Knuth, The Art of Computer Programming, Vol. 2, Seminumerical Algorithms, Second Edition, Addison-Wesley, 1981, pp 318. 原文是 “We might call it the granddaddy of all algorithms, because it is the oldest nontrivial algorithm that has survived to the present day.”
[2] 維基百科:輾轉相除法

Advertisements
本篇發表於 線性代數專欄, 應用之道 並標籤為 , , , 。將永久鏈結加入書籤。

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s