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

轉載自：http://www.cnblogs.com/drizzlecrj/archive/2007/09/14/892340.html

很老的東東了，其實也沒啥好整理的，網上很多資料了，就當備用把:-)

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)

2. 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 
這個題目要發現規律，不錯的題目