歐幾里得算法（展轉相除法）

<!-- Euclid算法 -->

歐幾里得算法即「展轉相除法」，用於求兩個正整數的最大公約數（gcd）。

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

這是遞歸形式的代碼，遞歸層數爲4.785lgN + 1.6723，N是max(a, b)，可見它遞歸不了多少層，不須要擔憂棧溢出的問題。

更精檢的形式以下：

int gcd(int a, int b)
{
    return (b? gcd(b, a%b): a);
}

咱們使用這段代碼的時候，須要保證a和b是非負整數，a、b同爲0的狀況要預先處理掉。事實上咱們不須要管a和b的大小關係，代碼在遞歸中，始終將大數放在前，咱們即便將小數放在前也就是多一層遞歸而已。

另外還有一種僅經過位運算就能夠達到目的的gcd：

int gcd(int a, int b)
{
    int t = 1, c, d;
    while(a != b)
    {
        if(a < b) swap(a, b);
        if(! (a & 1)) { a >>= 1; c = 1; } else c = 0;
        if(! (b & 1)) { b >>= 1; d = 1; } else d = 0;
        if(c && d) t <<= 1;
        else if(! c && ! d) a -= b;
    }
    return t * a;
}

代碼的思想是： 若是a和b均爲偶數，則將公因子2提出來累加到因數t上，最終做爲結果的係數一併返回； 若是a和b一個爲偶，一個爲奇，那麼不可能同時含有公因子2，將爲偶的數剔除2這個因子； 若是a和b都爲奇數，則使用「更相減損術」計算其gcd。

同時注意代碼始終保證a大於等於b，不然交換位置。其可行性在於gcd(a, b) = gcd(b, a)。 代碼儘管是循環形式，可是卻利用到了遞推關係，如上述三句話用公式表示以下： a、b爲偶數：gcd(a, b) = 2gcd(a/2, b/2); a、b其中之一爲偶數（不妨設a）：gcd(a, b) = gcd(a/2, b); a、b均爲奇數：gcd(a, b) = gcd(a-b, b); 這個關係是遞推的，因此能夠利用循環去迭代。 另外注意此次的終止條件是a == b，不是b == 0。

這個算法的複雜度可能高於使用「展轉相除」的歐幾里得算法，可是運行時間上可能相差不大。 由於代碼中使用的都是位運算（好比咱們斷定某數位奇數或偶數的時候用到的是(x & 1)這個條件，除以2乘以2用到的是>>=1和<<=1），速度極快，不是除法運算和取模運算比得了的。另外代碼很是巧妙的用c和d兩個標記記錄當前a、b的奇偶性。

<!-- 擴展Euclid算法 -->

歐幾里得算法很好理解，代碼也特別好寫。下面說的這種算法是歐幾里得算法的擴展，代碼同樣簡練，可是理解起來卻特別費勁。

擴展歐幾里得算法的用處之一，是用來解決「線性不定方程求整數根」問題（數論是解決整數問題的理論）。形如ax+by=c這樣的方程叫作線性不定方程，也就是咱們初高中學的直線方程。

咱們當前的目的是求ax+by=c的一組解，再考慮其餘的解能不能推出來。

設d = gcd(a, b)，那麼ax+by=d必定有解，數論基礎裏面有一條重要的結論：兩個數a、b（在這裏咱們都將其當作正整數來討論）的最大公約數，必定等於a與b的線性組合（ax+by）中最小的正整數。

如今假設咱們獲得了ax+by=d的一組解(x0，y0)，若是d是c的約數，那麼ax+by=c對應的一組解只要將x0、y0乘以相應的倍數就好。

反之若是d不是c的約數，那麼ax+by=c必定無解。爲何呢？由於ax+by必定是d的倍數，而c不是d的倍數，等號條件不可能成立。

下面咱們用拓展歐幾里得算法去求ax+by=d的解，這是咱們以前沒作的工做。

void gcd(int a, int b, int& d, int& x, int& y)
{
    if(!b) { d = a; x = 1; y = 0; }
    else
    {
        gcd(b, a%b, d, y, x);    // ( *  * )
        y -= x * (a/b);
    }
}

注意咱們在調用gcd(a, b, d, x, y)的時候，最終得出的d、x、y是ax+by=d的解。注意d、x、y傳遞的是引用，剛聲明的變量，直接傳進去就好，函數會改變它們的值。 如今解釋下代碼的含義：

很容易看出來這段代碼其中的一部分來自於咱們剛纔講過的普通歐幾里得算法，最終的d便是gcd(a, b)，它自從在遞歸底層被賦值之後就沒再變過。

咱們展轉相除去求得了d之後，將x賦值爲1，y賦值爲0，獲得了a x 1 + b x 0 = d這個等式，注意這裏b也是0，因此事實上y賦值爲何整數均可以。

可是最須要注意的是，此時的x和y並非ax+by=d的解，爲何？由於這時候的a和b已經不是原先的a和b了，通過數層的遞歸，a和b的值已經改變了。接下來要作的就是在回溯階段算出上一層的x和y值，若是算得出來，在遞歸結束時，咱們就能獲得原方程的x和y（只有遞歸的第一層的a和b是原方程ax+by=d的a和b）。

注意加星號的那一行，咱們在向下層遞歸的時候，傳遞的值分別是：b、a%b（這是爲了展轉相除算gcd），d（最終保存gcd）、y、x（注意不是x、y）。 當回溯到當前層時，x和y的值已經算出來了，它知足下一層遞歸中式子的關係，咱們須要找當前層的x、y和下層x、y的關係。

看這段代碼：gcd(b, a%b, d, y, x)，下層的參數列表是(int a, int b, int& d, int& x, int& y)。

這說明對這一層，有以下關係： by + (a % b)x = d。下面一層保證它成立，想一想最下面一層是a x 1 + b x 0 = d。

咱們但願對這一層有ax + by = d，顯然從上面那個式子得不出這個式子，咱們須要改變一下x或y。

正確的作法是將y減去x(a/b)，咱們來驗證一下： a x + b (y - x(a/b)) = a x + b (y - x(a - a % b) / b) = b y - (a % b) x = d

可見對本層的a、b，將y減去x(a/b)，x不變，就會獲得知足ax + by = d的x和y。

可能有人會好奇這個x(a / b)是怎麼來的，這個我以爲咱們不用糾結（或者數學特別厲害的人能想辦法找到這個差值），這個差值是算法的發明者找出來的，咱們記下之後能寫出算法來就好，畢竟咱們只是利用擴展歐幾里得算法解決問題。

以後的問題：

咱們費了一番功夫，獲得了ax + by = d的一組(x, y)，若是d是c的約數，那麼知足ax + by = c的x、y值咱們也拿到了。

如今咱們來找通解，將獲得的這一組x和y記做x0、y0。

另外一組x、y若是也知足ax + by = c，那麼會有等式 a x + b y = a x0 + b y0。 進而有 a(x - x0) = b (y0 - y)。 咱們將a、b約去d後，a'(x - x0) = b'(y0 - y)，這個式子裏因爲a'和b'互素，式子要想成立，必須使x - x0是b'倍數，y - y0是a'倍數。 即： x = k b' + x0, y = k a' + y0。 注意咱們取解的時候，這個k一旦定了，新的x和y要同時算出來，一個k值對應一組解。

到這裏，咱們能求出ax+by = c的通解了。

<!-- 更多變元 -->

如今擴展一下問題，咱們來求ax + by + cz = d這個方程是否有解。

利用gcd(a, b, c) = gcd(gcd(a, b), c)就能夠求得a、b、c的最大公約數，若是它也是d的約數，那麼就有解，不然無解。

這麼看來還能夠推廣一下結論：n個整數（這裏指正整數）的gcd，是這n個整數線性組合的最小正整數。

<!-- 一些廢話 -->

數論的內容晦澀難懂，機關重重，解題方式也是變化無窮，對數學思惟要求很高。也是我本身又愛又恨的一部份內容（愛是以爲很巧妙，恨是想不出來）。

我儘可能把這些內容寫得清楚明白，但願跟我同樣的算法愛好者在入門的時候少走彎路，不被抹殺掉搞算法的興趣和激情。

最後吐槽一下開源中國的編輯器，實在是用着難受。。。多是我不怎麼會用吧，各類符號不能使用，代碼也莫名奇妙就給我註釋掉了。