最大公約數(Gcd)兩種算法(Euclid && Stein)轉載

寫的不錯,代碼看起來很舒服。

  1. 歐幾里德算法和擴展歐幾里德算法
    歐幾里德算法
    歐幾里德算法又稱輾轉相除法,用於計算兩個整數a,b的最大公約數。其計算原理依賴於下面的定理:

定理:gcd(a,b) = gcd(b,a mod b)

證明:a可以表示成a = kb + r,則r = a mod b
假設d是a,b的一個公約數,則有
d|a, d|b,而r = a – kb,因此d|r
因此d是(b,a mod b)的公約數

假設d 是(b,a mod b)的公約數,則
d | b , d |r ,但是a = kb +r
因此d也是(a,b)的公約數

因此(a,b)和(b,a mod b)的公約數是一樣的,其最大公約數也必然相等,得證

歐幾里德算法就是根據這個原理來做的,其算法用C++語言描述爲:

int Gcd(int a, int b)
{
if(b == 0)
return a;
return Gcd(b, a % b);
}

當然你也可以寫成迭代形式:
int Gcd(int a, int b)
{
while(b != 0)
{
int r = b;
b = a % b;
a = r;
}
return a;
}

本質上都是用的上面那個原理。

補充: 擴展歐幾里德算法是用來在已知a, b求解一組p,q使得p * a+q * b = Gcd(a, b) (解一定存在,根據數論中的相關定理)。擴展歐幾里德常用在求解模線性方程及方程組中。下面是一個使用C++的實現:
int exGcd(int a, int b, int &x, int &y)
{
if(b == 0)
{
x = 1;
y = 0;
return a;
}
int r = exGcd(b, a % b, x, y);
int t = x;
x = y;
y = t – a / b * y;

return r;

}

把這個實現和Gcd的遞歸實現相比,發現多了下面的x,y賦值過程,這就是擴展歐幾里德算法的精髓。
可以這樣思考:
對於a’ = b, b’ = a % b 而言,我們求得 x, y使得 a’x + b’y = Gcd(a’, b’)
由於b’ = a % b = a – a / b * b (注:這裏的/是程序設計語言中的除法)
那麼可以得到:
a’x + b’y = Gcd(a’, b’) ===>
bx + (a – a / b * b)y = Gcd(a’, b’) = Gcd(a, b) ===>
ay +b(x – a / b*y) = Gcd(a, b)
因此對於a和b而言,他們的相對應的p,q分別是 y和(x-a/b*y)

  1. Stein算法
    歐幾里德算法是計算兩個數最大公約數的傳統算法,他無論從理論還是從效率上都是很好的。但是他有一個致命的缺陷,這個缺陷只有在大素數時纔會顯現出來。

考慮現在的硬件平臺,一般整數最多也就是64位,對於這樣的整數,計算兩個數之間的模是很簡單的。對於字長爲32位的平臺,計算兩個不超過32位的整數的模,只需要一個指令週期,而計算64位以下的整數模,也不過幾個週期而已。但是對於更大的素數,這樣的計算過程就不得不由用戶來設計,爲了計算兩個超過64位的整數的模,用戶也許不得不採用類似於多位數除法手算過程中的試商法,這個過程不但複雜,而且消耗了很多CPU時間。對於現代密碼算法,要求計算128位以上的素數的情況比比皆是,設計這樣的程序迫切希望能夠拋棄除法和取模。 (注:說到拋棄除法和取模,其實輾轉相除法可以寫成減法的形式)

Stein算法由J. Stein 1961年提出,這個方法也是計算兩個數的最大公約數。和歐幾里德算法 算法不同的是,Stein算法只有整數的移位和加減法,這對於程序設計者是一個福音。

爲了說明Stein算法的正確性,首先必須注意到以下結論:

gcd(a,a) = a,也就是一個數和他自身的公約數是其自身
gcd(ka,kb) = k gcd(a,b),也就是最大公約數運算和倍乘運算可以交換,特殊的,當k=2時,說明兩個偶數的最大公約數必然能被2整除。

有了上述規律就可以給出Stein算法如下:

如果A=0,B是最大公約數,算法結束
如果B=0,A是最大公約數,算法結束
設置A1 = A、B1=B和C1 = 1
如果An和Bn都是偶數,則An+1 =An /2,Bn+1 =Bn /2,Cn+1 =Cn *2(注意,乘2只要把整數左移一位即可,除2只要把整數右移一位即可)
如果An是偶數,Bn不是偶數,則An+1 =An /2,Bn+1 =Bn ,Cn+1 =Cn (很顯然啦,2不是奇數的約數)
如果Bn是偶數,An不是偶數,則Bn+1 =Bn /2,An+1 =An ,Cn+1 =Cn (很顯然啦,2不是奇數的約數)
如果An和Bn都不是偶數,則An+1 =|An -Bn|,Bn+1 =min(An,Bn),Cn+1 =Cn
n++,轉4
這個算法的原理很顯然,所以就不再證明了。現在考察一下該算法和歐幾里德方法效率上的差別。
給出一個C++的實現:
int Gcd(int a, int b)
{
if(a == 0) return b;
if(b == 0) return a;
if(a % 2 == 0 && b % 2 == 0) return 2 * gcd(a >> 1, b >> 1);
else if(a % 2 == 0) return gcd(a >> 1, b);
else if(b % 2 == 0) return gcd(a, b >> 1);
else return gcd(abs(a – b), Min(a, b));
}

考慮歐幾里德算法,最惡劣的情況是,每次迭代a = 2b -1,這樣,迭代後,r= b-1。如果a小於2N,這樣大約需要 4N次迭代。而考慮Stein算法,每次迭代後,顯然AN+1BN+1≤ ANBN/2,最大迭代次數也不超過4N次。也就是說,迭代次數幾乎是相等的。但是,需要注意的是,對於大素數,試商法將使每次迭代都更復雜,因此對於大素數Stein將更有優勢

練習:
OJ上面的赤裸裸的Gcd算法的題不多,大多都是套了一個外殼。
找了兩道,可以試試看
HDOJ 2028 Lowest Common Multiple Plus 這個是求n個數的最小公倍數(有了最大公約數,最小公倍數應該很容易了)
ZJU 2678 Bishops on a Toral Board 這個題目要發現規律,不錯的題目

点赞