欧几里得算法（展转相除法）

<!-- 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个整数线性组合的最小正整数。

<!-- 一些废话 -->

数论的内容晦涩难懂，机关重重，解题方式也是变化无穷，对数学思惟要求很高。也是我本身又爱又恨的一部份内容（爱是以为很巧妙，恨是想不出来）。

我尽可能把这些内容写得清楚明白，但愿跟我同样的算法爱好者在入门的时候少走弯路，不被抹杀掉搞算法的兴趣和激情。

最后吐槽一下开源中国的编辑器，实在是用着难受。。。多是我不怎么会用吧，各类符号不能使用，代码也莫名奇妙就给我注释掉了。