歐幾里得（Euclid）與拓展的歐幾里得演算法
歐幾里得（Euclid）與拓展的歐幾里得演算法
- 歐幾里得（Euclid）與拓展的歐幾里得演算法
- 歐幾里得演算法
- 原理
- 實現
- 拓展的歐幾里得演算法
- 原理
- 遞迴求解
- 迭代求解
- 歐幾里得演算法
歐幾里得演算法
原理
歐幾里得演算法是一種快速計算最大公約數的演算法，對於任意的兩個數\((a,b)\)，其最大公約數表示為\(gcd(a,b)\)，根據歐幾里得演算法，\(gcd(a,b)=gcd(b,a\%b)\)。證明如下：
如果\(b>a\)，顯然成立；因此只需考慮\(b<a\)的情況。根據初等數學知識，可知\(a,b\)的關係可表示為\(a=qb+r\)，其中\(q\)為商，\(r\)為餘數。
對於\((a,b)\)的最大公約數\(g1=gcd(a,b)\)，當然\(g1|a,g1|b\)（\(g1|a\)表示g1整除a），所以易知對於\(r=a-qb\)，同樣滿足\(g1|r\)；
又因為\(a\%b=r\)，所以對於\(a,b\)的最大公約數g1，同樣滿足\(g1|a\%b,g1|b\)，即\((b,a\%b)\)的最大公約數至少為\(g1\)，即\(gcd(b,a\%b)>g1=gcd(a,b)\)。
反過來，對於\((b,a\%b)\)的最大公約數\(g2=gcd(b,a\%b)\)，同樣滿足\(g2|a, g2|b\)，即\(gcd(a,b)>g2=gcd(b,a\%b)\)。
因此\(gcd(a,b)=gcd(b,a\%b)\)證明成立。下面對該演算法進行實現。
實現
#include <iostream> using namespace std; int euclid(int a, int b) { if (b!=0) { return euclid(b, a%b); } else { return a; } } int main() { int a(0),b(0); cin >> a >> b; cout << euclid(a,b); return 0; }
拓展的歐幾里得演算法
原理
拓展的歐幾里得演算法在密碼學中有著重要的應用，現給出定理：
對正整數a，b；總是存在一組整數X,Y，使得\(Xa+Yb=gcd(a,b)\)成立，且\(gcd(a,b)\)為滿足這種條件的最小整數。
這裡不對該定理進行證明，歐幾里得演算法給出了在已知a，b的情況下求\(gcd(a,b)\)的方法，但是如果想要求得X，Y的值，就要求助於拓展的歐幾里得演算法。怎麼才能從歐幾里得演算法的計算過程當中得到我們想要求解的值呢？我們再次詳細回顧歐幾里得演算法的求解過程。
對於已知整數\(a,b\)，我們的演算法求解過程如下：
| \(a\) | \(b\) | 餘數\(r\) | 商\(q\) |
|---|---|---|---|
| a | b | \(r1=a\%b\) | \(q1=a/b\) |
| b | r1 | \(r2=b\%r1\) | \(q2=b/r1\) |
| r1 | r2 | \(r3=r1\%r2\) | \(q3=r1/r2\) |
| ... | ... | ... | ... |
| \(r_{n-1}\) | \(r_n\) | \(r_{n+1}=r_{n-1}\%r_n\) | \(q_{n+1}=r_{n-1}/r_n\) |
逐步計算，直到某一步出現\(r_{n-1}\%r_{n}=0\)的情況，這時候就找到了最大公約數，最大公約數即為\(r_n\)，以上就是歐幾里得演算法的全過程。通過這個過程當中多產生的一些中間結果我們能不能求得\(X,Y\)的值呢？下面進行兩種求解方法的推導。
遞迴求解
根據上面的表格我知道，\(Xa+Yb=gcd(a,b)\)，並且對於中間所求解的每一步我們所得到的\(r_i\)都滿足\(X_i\cdot r_i+Y_i\cdot r_{i+1}=gcd(a,b)\)，因為很明顯每一對\(r_i,r_{i+1}\)都滿足最大公約數為\(gcd(a,b)\)，這也是歐幾里得演算法的原理。
我們試著尋找\((X_i,Y_i),(X_{i+1},Y_{i+1})\)之間的遞推關係，由以上闡述可知:
\[
\begin{cases}
X_i\cdot r_i+Y_i\cdot r_{i+1}=gcd(a,b), &\text{(1)}\\
X_{i+1}\cdot r_{i+1}+Y_{i+1}\cdot r_{i+2}=gcd(a,b),&\text{(2)}
\end{cases}
\]
為了將上式轉換成\(r_i\)的方程組，我們使用\(r_{i+1},r_{i+2}來表示r_i\)，通過以上可知\(r_i=r_i/r_{i+1}\cdot r_{i+1}+r_{i+2}\)，將該式帶入上式（1），並將兩式合併可得：
\[
X_i\cdot (r_i/r_{i+1}\cdot r_{i+1}+r_{i+2})+Y_i\cdot r_{i+1}=X_{i+1}\cdot r_{i+1}+Y_{i+1}\cdot r_{i+2}
\]
進一步化簡可得：
\[
(X_i\cdot r_i/r_{i+1}+Y_i)\cdot r_{i+1}+X_i\cdot r_{i+2}=X_{i+1}\cdot r_{i+1}+Y_{i+1}\cdot r_{i+2}
\]
根據係數相等的原則可得：
\[
\begin{cases}
X_i\cdot r_i/r_{i+1}+Y_i=X_{i+1}, &\text{(1)}\\
X_i=Y_{i+1},&\text{(2)}
\end{cases}
\]
以上就得到了\((X_i,Y_i),(X_{i+1},Y_{i+1})\)之間的遞推關係，那麼我們接下來的工作就是找到一對可以求出其值的\((X_i,Y_i)\)，通過以上可知當出現某一次計算使得\(r_{i}\%r_{i+1}=0\)時，我們可知對於\(X_i\cdot r_i+Y_i\cdot r_{i+1}=gcd(a,b)\)，滿足\(gcd(a,b)=r_{i+1}\)，那麼很顯然\(X_i=0,Y_i=1\)。於是我們就得到了一對\((X_i,Y_i)\)的值，我們已經知道了最後一對\(r_{i},r_{i+1}\)所對應的\((X_i,Y_i)\)才能夠推知前面的值，所以我們的推導時從後往前推的，因此我們將上面的遞推關係稍微變換一下形式：
\[
\begin{cases}
X_i=Y_{i+1},&\text{(1)}\\
Y_i=X_{i+1}-Y_{i+1}\cdot r_i/r_{i+1}, &\text{(2)}
\end{cases}
\]
此時我們就得到了推導關係和初值，通過計算我們就可以求得滿足\(Xa+Yb=gcd(a,b)\)的\(X,Y\)值。下面通過程式碼對其進行實現：
#include <iostream>
using namespace std;
void extEuc(int a, int b, int&x, int&y)
{
int rx,ry;
int r(0);
if (a%b==0)
{
x=0;y=1;
return;
}
else
{
r = a%b;
extEuc(b, r, rx, ry);
x = ry;
y = rx - ry*a/b;
return;
}
}
int main()
{
int a,b;
int x,y;
cin >> a >> b;
extEuc(a, b, x, y);
cout << x <<' '<< y << endl;
return 0;
} 輸入42 2017可求得輸出為-48，1。
迭代求解
較多的遞迴呼叫可能會影響計算速度，所以我們接下來推一下迭代的計算方式，已知上面表格中所列歐幾里得演算法的計算步驟。已知\(gcd(a,b)\)是滿足該集合的最小值\(\lbrace Xa+Yb | X,Y\in Z \rbrace\)，已知對於每一步所產生的餘數均能被\(gcd(a,b)\)整除，現在考慮每一步迭代所產生的餘數滿足的等式：
\[
\begin{cases}
r_i=X_i\cdot a+Y_i\cdot b \\
r_{i+1}=X_{i+1}\cdot a + Y_{i+1}\cdot b
\end{cases}
\]
且已知\(r_i, r_{i+1}\)滿足\(r_{i-1}=r_{i-1}/r_i\cdot r_i+r_{i+1}\)，將上面兩式代入到該式，可得：
\[
r_{i-1}=(r_{i-1}/r_i\cdot X_i+X_{i+1})\cdot a+(r_{i-1}/r_i\cdot Y_i +Y_{i+1})\cdot b.
\]
值得注意的是此處的\(X_i,Y_i\)與遞迴方法中的值含義不同。根據上式可推知以下遞推關係：
\[\begin{cases}
X_{i+1}=X_{i-1}-r_{i-1}/r_i\cdot X_i \\
Y_{i+1}=Y_{i-1}-r_{i-1}/r_i\cdot Y_i
\end{cases}\]
已知中間的遞推關係，關鍵是考慮如何判斷迴圈的起始值和結束條件，對於\(a,b\)也可看做是餘數\(r_i\)，那麼對於\(a,b\)來說，其滿足的值為：
\[\begin{cases}
a = a\cdot 1 + b\cdot 0\\
b = a\cdot 0 + b\cdot 1
\end{cases}\]
所以就得到了兩對\((X,Y)\)的值，分別為\((1,0),(0,1)\)，並且已知\(r_i\)之間的遞推關係為\(r_{i+1}=r_{i-1}\%r_i\)。我們也知道迴圈結束的條件為\(r_i=gcd(a,b)\)，其最後的形式為\(Xa+Yb=gcd(a,b)\)，其直接判斷方式為\(r_{i-1}\%r_i=0\)，然後我們就得到了最終的\(X,Y\)值，根據以上遞推形式，我們有以下實現：
#include <iostream>
using namespace std;
int main()
{
int a,b;
int x1(1),y1(0),x2(0),y2(1);
int temp;
cin >> a >> b;
while (a%b!=0)
{
temp = x2;
x2 = x1 - a/b*x2;
x1 = temp;
temp = y2;
y2 = y1 - a/b*y2;
y1 = temp;
temp = a%b;
a = b;
b = temp;
}
cout << x2 <<' '<<y2<<endl;
return 0;
} 輸入42 2017可求得輸出為-48，1。
這裡有一個用盡可能多的程式語言實現求逆元的網站，大家也可以參考這裡的不同實現。
參考文獻
[1] Katz J，Lindel Y．Introduction to Modern Cryptography—Principle and Protocol現代密碼學——原理與協議【M】任偉．北京：國防工業出版社．2010：10-15．