（ 数论专题 ）【 扩展欧几里得 】

（ 数论专题 ）【 扩展欧几里得 】

扩展欧几里德算法

谁是欧几里德？自己百度去

先介绍什么叫做欧几里德算法

有两个数 a b，现在，我们要求 a b 的最大公约数，怎么求？枚举他们的因子？不现实，当 a b 很大的时候，枚举显得那么的naïve ，那怎么做？

欧几里德有个十分又用的定理： gcd(a, b) = gcd(b , a%b) ，这样，我们就可以在几乎是 log 的时间复杂度里求解出来 a 和 b 的最大公约数了，这就是欧几里德算法，用 C++ 语言描述如下：
 

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

由于是用递归写的，所以看起来很简洁，也很好记忆。那么什么是扩展欧几里德呢？

现在我们知道了 a 和 b 的最大公约数是 gcd ，那么，我们一定能够找到这样的 x 和 y ，使得: a*x + b*y = gcd 这是一个不定方程（其实是一种丢番图方程），有多解是一定的，但是只要我们找到一组特殊的解 x0 和 y0 那么，我们就可以用 x0 和 y0 表示出整个不定方程的通解（ 下面两个通解非常重要！！！ ）：

        x = x0 + (b/gcd)*t

        y = y0 – (a/gcd)*t

    为什么不是：

        x = x0 + b*t

        y = y0 – a*t

这个问题也是在今天早上想通的，想通之后忍不住喷了自己一句弱逼。那是因为：

b/gcd 是 b 的因子， a/gcd 是 a 的因子是吧？那么，由于 t的取值范围是整数，你说 (b/gcd)*t 取到的值多还是 b*t 取到的值多？同理，(a/gcd)*t 取到的值多还是 a*gcd 取到的值多？那肯定又要问了，那为什么不是更小的数，非得是 b/gcd 和a/gcd ？

注意到：我们令 B = b/gcd ， A = a、gcd ， 那么，A 和 B 一定是互素的吧？这不就证明了 最小的系数就是 A 和 B 了吗？要是实在还有什么不明白的，看看《基础数论》（哈尔滨工业大学出版社），这本书把关于不定方程的通解讲的很清楚

现在，我们知道了一定存在 x 和 y 使得 ： a*x + b*y = gcd ， 那么，怎么求出这个特解 x 和 y 呢？只需要在欧几里德算法的基础上加点改动就行了。

我们观察到：欧几里德算法停止的状态是： a= gcd ， b = 0 ，那么，这是否能给我们求解 x y 提供一种思路呢？因为，这时候，只要 a = gcd 的系数是 1 ，那么只要 b 的系数是 0 或者其他值（无所谓是多少，反正任何数乘以 0 都等于 0 但是a 的系数一定要是 1），这时，我们就会有： a*1 + b*0 = gcd

当然这是最终状态，但是我们是否可以从最终状态反推到最初的状态呢？

假设当前我们要处理的是求出 a 和 b的最大公约数，并求出 x 和 y 使得 a*x + b*y= gcd ，而我们已经求出了下一个状态：b 和 a%b 的最大公约数，并且求出了一组x1 和y1 使得： b*x1 + (a%b)*y1 = gcd ， 那么这两个相邻的状态之间是否存在一种关系呢？

我们知道： a%b = a - (a/b)*b（这里的 “/” 指的是整除，例如 5/2=2 , 1/3=0），那么，我们可以进一步得到：

        gcd = b*x1 + (a-(a/b)*b)*y1

               = b*x1 + a*y1 – (a/b)*b*y1

               = a*y1 + b*(x1 – a/b*y1)

    对比之前我们的状态：求一组 x 和 y 使得：a*x + b*y = gcd ，是否发现了什么？

    这里：

        x = y1

        y = x1 – a/b*y1

    以上就是扩展欧几里德算法的全部过程，依然用递归写：

int ex_gcd( int a, int b, int &x, int &y )
{
    if ( b==0 ) {
        x = 1;
        y = 0;
        return a;
    }
    int re = ex_gcd(b,a%b,y,x);
    y -= x*(a/b);
    return re;
}

依然很简短，相比欧几里德算法，只是多加了几个语句而已。

这就是理论部分，欧几里德算法部分我们好像只能用来求解最大公约数，但是扩展欧几里德算法就不同了，我们既可以求出最大公约数，还可以顺带求解出使得： a*x + b*y = gcd 的通解 x 和 y

 

---

 

例题1   青蛙的约会   POJ - 1061

题意： 两只青蛙在网上相识了，它们聊得很开心，于是觉得很有必要见一面。它们很高兴地发现它们住在同一条纬度线上，于是它们约定各自朝西跳，直到碰面为止。可是它们出发之前忘记了一件很重要的事情，既没有问清楚对方的特征，也没有约定见面的具体位置。不过青蛙们都是很乐观的，它们觉得只要一直朝着某个方向跳下去，总能碰到对方的。但是除非这两只青蛙在同一时间跳到同一点上，不然是永远都不可能碰面的。为了帮助这两只乐观的青蛙，你被要求写一个程序来判断这两只青蛙是否能够碰面，会在什么时候碰面。
我们把这两只青蛙分别叫做青蛙A和青蛙B，并且规定纬度线上东经0度处为原点，由东往西为正方向，单位长度1米，这样我们就得到了一条首尾相接的数轴。设青蛙A的出发点坐标是x，青蛙B的出发点坐标是y。青蛙A一次能跳m米，青蛙B一次能跳n米，两只青蛙跳一次所花费的时间相同。纬度线总长L米。现在要你求出它们跳了几次以后才会碰面。

思路： 设时间为 t 时两青蛙相遇，则两个青蛙的位置分别为(x+mt)mod L、(y+nt) mod L，相遇即是(x+mt)%L=(y+nt)%L，

即(m-n)*t+k*L=y-x   ， 式中t为经过的时间，k为多走了几圈。
OK，现在 t和k是未知量， 已经符合ax+by=c的方程了，设a=m-n，b=L，c=y-x，然后套用模板求出ax+by=gcd(a,b)的特解，然后方程两边乘以c/gcd， 这样解出的x也应该乘以c/gcd， 注意t>0，所以要用通解公式得出最小正整数。最后注意用long long~

其中用到了一个公式：

只要我们找到一组特殊的解 x0 和 y0 那么，我们就可以用 x0 和 y0 表示出整个不定方程的通解：

        x = x0 + (b/gcd)*t

        y = y0 – (a/gcd)*t

用这个通解来求出最小时间， 也就是x最小正整数的解。

代码：

#include <iostream>
 
using namespace std;
 
typedef long long ll;
 
ll exgcd( ll a, ll b, ll &x, ll &y ) // 扩展欧几里得模板
{
    if ( b==0 ) {
        x = 1;
        y = 0;
        return a;
    }
    ll re = exgcd(b,a%b,y,x);
    y -= x*(a/b);
    return re;
}
 
int main()
{
    ll x,y,m,n,L;
    cin >> x >> y >> m >> n >> L;
    ll a = m-n;
    ll b = L;
    ll c = y-x;
    if ( a<0 ) { // 最好把方程弄成正的，
        a = -a;  // 方程两边同乘-1
        c = -c;
    }
    ll gcd = exgcd(a,b,x,y);
    if ( c%gcd!=0 ) {  // c不是gcd的倍数，所以 ax+by=c 没有解
        cout << "Impossible" << endl;
    }
    else {
        x = x*c/gcd;   // 把 ax+by=gcd 化成 ( c/gcd )*(ax+by) = c ，特解也乘以倍数
        ll t = b/gcd;   // 通解需要加的倍数， 这个不用乘那个（c/gcd）
        if ( x>=0 ) {
            x = x%t ;
        }
        else {
            x = x%t+t;
        }
        cout << x << endl;
    }
 
    return 0;
}

 

例题2  Romantic   HDU - 2669

题意：给出a，b 求满足 ax+by=1 方程的解（ 要求x是最小正整数 ）。 根据扩展欧几里得，若gcd(a,b)不为1，该方程一定无解，否则求出特解x0，y0， 根据通解公式求得答案。

        x = x0 + (b/gcd)*t

        y = y0 – (a/gcd)*t

代码：

#include <iostream>
 
using namespace std;
 
int exgcd( int a, int b, int &x, int &y )
{
    if ( b==0 ) {
        x = 1;
        y = 0;
        return a;
    }
    int re = exgcd(b,a%b,y,x);
    y -= x*(a/b);
    return re;
}
 
int main()
{
    int a,b,x,y;
    while ( cin >> a >> b ) {
        int gcd = exgcd(a,b,x,y);
        if ( gcd!=1 ) {  // 如果 a b 的gcd不为1 一定无解
            cout << "sorry" << endl;
        }
        else {
            int t = b/gcd; // x通解需要加的那个倍数
            if ( x>=0 ) {
                int k = x/t; // 需要缩小多少个k
                x = x%t;
                y += k*(a/gcd);
                cout << x << " " << y << endl;
            }
            else {
                int k = x/t+1; // 需要扩大多少个k
                x = x%t + t;
                y -= k*(a/gcd);
                cout << x << " " << y << endl;
            }
        }
    }
 
    return 0;
}

附加一个我对青蛙走路的理解：关于青蛙走路问题的理解（ 扩展欧几里得的应用 ）