拓展歐幾里得
歐幾里得演算法
歐幾里得演算法(也稱輾轉相除法)是目前已知求最大公約數的最快通用演算法,具有程式碼複雜度低、易理解、用途廣等諸多優點,也是OI中不可或缺的的一種演算法。(當然更相減損術也是)
若有兩個數a和b,需要求a和b的最大公約數,怎麼辦?
列舉它們的因子? 小資料可以,大資料的話,這個O(n)的演算法就圖森破了。
分解質因數?在大數面前也是拿衣服的。
總不能不搞吧,上帝喊了一聲“效率”,於是我們便有了O(log n)效率極高的歐幾里得演算法。
引理
歐幾里得有個非常強的定理,即gcd(a,b)=gcd(b,a mod b),讓我們來證明一下
約定:mod <=> 取餘,gcd <=> 最大公約數,| <=> 整除
假設a、b的公約數為k,a = bx + y
則 k | a,k | b,a mod b=y
因為 k | b,所以 k |bx,又因為 k | a,所以 k | (a - bx),即 k | y
而a mod b = y,所以 k | a mod b
再假設b、a mod b的公約數為kk,同理得 kk | a,所以(a,b)和(b,a mod b)的公約數是相同的,所以它們的最大公約數也是相同的.
所以gcd(a,b)=gcd(b,a mod b)
時間複雜度(討論情況為a>=b)
a mod b必然是小於a/2的,而上一次的b會變成下一次的a,上一次的a mod b會變成下一次的b,最壞情況也就是b在a/2附近,即a mod b在a/2附近。在最壞情況時每次的值也會減少一半,所以說時間複雜度是O(log n)的。(實際中往往會更低)
演算法實現
有了以上的引理演算法實現就非常簡單了,遞迴和非遞迴形式都可以。任何數與0的最大公約數是它本身,所以一直操作直到b = 0時的a就是原a,b的最大公約數了。
虛擬碼:
while(b){ r=a%b; a=b; b=r; } //最終的a即為原a,b最大公約數
擴充套件歐幾里得演算法
擴充套件歐幾里得演算法的功能就更強大了,它可以用來求二元一次方程的通解,還可以用來求乘法逆元。
在此順便簡介一下乘法逆元:
若有 a*x ≡ 1 (mod m),則稱 x 為a關於m的乘法逆元,等價式 a * x+m * y = 1 (y=a*x/m)
這就也是個二元一次方程了,ExGcd可搞。
引理
裴蜀定理:若ax+by = z,則 gcd(a,b)| z
再順手證明一下裴蜀定理:
設k = gcd(a,b),則 k | a, k | b,根據整除的性質,有 k | (ax+by)
設 s為ax+by的最小正數值
再設 q = [a / s](a整除s的值);r = a mod s = a-q(ax+by) = a(1 - qx)+b(-qy);
由此可見r也為a,b的線性組合;(ax+by稱為a,b的線性組合)
又因為s為a,b的線性組合的最小正數值,0<= r < s,所以r的值為0,即 a mod s = r =0;s | a;
同理可得 s | b,則 s | k;
又因為 k | (ax+by),s為ax+by的最小正數值,所以 k | s;
因為 s | k,k | s,所以s = k;
原命題得證。
如何求解(以下討論a>b)
顯然當 b=0,gcd(a,b)=a。此時 x=1,y=0;
當a>b>0 時
設 ax1+ by1= gcd(a,b);
bx2+ (a mod b)y2= gcd(b,a mod b);
根據歐幾里德原理有 gcd(a,b) = gcd(b,a mod b);
則:ax1+ by1= bx2+ (a mod b)y2;
即:ax1+ by1= bx2+ (a - [a / b] * b)y2 = ay2+ bx2- [a / b] * by2;
(a mod b = a - [a / b]*b;[a / b]代表a整除b)
也就是ax1+ by1 = ay2 + b(x2- [a / b] *y2);
根據恆等定理得:x1=y2;y1=x2- [a / b] *y2;
這樣我們就得到了求解 x1,y1 的方法:x1,y1 的值基於 x2,y2
由引理我們知道:ax+by = z,z為gcd(a,b)若干倍。
所以我們先求解ax+by = gcd(a,b),再將求出的解乘以 z/gcd(a,b)就好了。
演算法實現
我們依舊遞迴實現,因為上一次的x,y與下一次的x,y的值有關,所以我們從x=1,y=0(此時b=0)的情況開始遞迴上來,套公式就好了。(a,b下去,x,y上來)
虛擬碼:
int ExGcd(int a,int b,int &x,int &y){ //x,y可設為全域性
if(!b){
x=1;
y=0;
return a;
}
int r=ExGcd(b,a%b,x,y);
t=x;
x=y;
y=t-a/b*y;
return r;
}
本文轉載自@leader_one