(扩展)欧几里得

2017-12-13

首先是欧几里得定理，即gcd(a,b)=gcd(b,a%b)

证明一下吧
c=gcd(a,b) d=gcd(b,a%b)
假设a=b*k+t;
k是商，t是余数
那么d=gcd(b,a%b)=gcd(b,t)
因为d|b,d|t，并且a=b*k+t
所以d|a，即d是a，b的公因数，即d小于等于c
又因为c|a,c|b，并且t=a-b*k
所以c|t，即c是b，t的公因数，即c小于等于d
所以说c==d的

或许我们大家都知道gcd(a,b)=gcd(|a|,|b|)，所以我们在求最大公约数
的时候a和b通常都是取正数的。
对于gcd(a,b)=gcd(b,a-b)(a>b)，这里可以用上述方法进行证明的。

我之前在《编程之美》上求解a与b的最大公约数中看到了三种解法，在这里简单
的介绍一下吧。
（1）直接利用__gcd(a,b)，这是c++里面的库函数，包含在头文件<algorithm>
中，我查看了一下他的源码，是用迭代的方法求解的。
类似于这个吧！
int gcd1(int x,int y){
    while(y){
        int r=x%y;
        x=y;
        y=r;
    }
    return x;
}
当然，我们也可以用递归的方法
int dg1(int x,int y){
    return y==0?x:dg1(y,x%y);
}

（2）学过计算机组成原理的人都知道，求余操作耗时是比较长的，那么我们可以用
上面提到的减法操作，并且这种方法对于大整数也是可以求得的，代码我就不赘述了。

（3）不知道大家还记得这样一个定理吗？即gcd(a1*k,b1*k)=k*gcd(a1,b1)；
还有一个，如果说我们的b%k！=0，那么gcd(a*k,b)=gcd(a,b)。这里的证明我也
不是很清楚。
接下来我要介绍的方法就是和这个是有关系的，我们大家都知道，在计算机里位操作
是非常快的，左移<<相当于*2，右移>>相当于/2，是吧！那么我们可以利用移位操作
来加快我们的执行速度。直接看代码吧！
int gcd(int x,int y){
    if (x<y) return gcd(y,x);
    if (y==0) return x;
    if (iseven(x)&&iseven(y)){
        return gcd(x>>1,y>>1)<<1;
    }if (iseven(x)){
        return gcd(x>>1,y);
    }if (iseven(y)){
        return gcd(x,y>>1);
    }else{
        return gcd(y,x-y);
    }
}
对于最后一种情况x与y都是奇数的话，它们俩进行减操作，最后的结果一定是
一个偶数吧！还有，我们在判断奇偶性的时候，只要看二进制表示的最后一位是
0还是1即可，所以就有了下面的代码：
bool iseven(int n){
    return !(n&1);
}
这就是我介绍的三种方法。

接下来就是扩展欧几里得了 ，说实话，我对这个模板并没有很好的理解，我仅仅能够按照他的代码自己模拟，但是至于代码怎么得来的以及它的证明过程我并不是很清楚。我参考的是我的一个同学的博客。
http://blog.csdn.net/yoer77/article/details/69568676

在数论中，裴蜀等式（英语：Bézout’s identity）或贝祖定理（Bézout’s
 lemma）是一个关于最大公约数（或最大公约式）的定理。裴蜀定理得名于法国
数x和y的线性丢番图方程（称为裴蜀等式）： 

ax + by = m 有整数解时当且仅当m是d的倍数。
（这句话很关键，注意这里的'''当且仅当'''）

裴蜀等式有解时必然有无穷多个整数解，每组解x、y都称为裴蜀数，可用扩展
欧几里得算法(Extended Euclidean algorithm)求得。
例如，12和42的最大公因数是6，则方程12x+42y=6有解。事实上有
(-3)×12 + 1×42 = 6及4×12 + (-1)×42 = 6。
特别来说，方程 ax+by=1 有整数解当且仅当整数a和b互素。 
裴蜀等式也可以用来给最大公约数定义： d其实就是最小的可以写成
ax+by形式的正整数。(这个定义的本质是整环中“理想”的概念。因此
对于多项式整环也有相应的裴蜀定理。)

扩展欧几里得就是在求出gcd(a,b)的同时求出a*x+b*y=gcd(a,b)的一个解。

用类似辗转相除法，求二元一次不定方程63x+22y=1的整数解。
首先
63=22*2+19
22=19*1+3
19=3*6+1
这里的1就是我们的最大公约数了
然后我们左右换个表示方式
19=63+22*(-2)
3=22+19*(-1)
1=19+3*(-6)
最后我们一步步的回带进去
1=19+3*(-6)
1=19+[22+19*(-1)]*(-6)
1=19*7+22*(-6)
1=[63+22*(-2)]*7+22*(-6)
1=63*7+22*(-20)
这样我们就求解出来了
x=7,y=-20
这就是我们所要求得的答案啊！

设：a>b。
显然当 b=0，gcd(a,b)=a。此时 x=1，y=0;
ab!=0 时
对于一般的情况而言
a*x1+b*y1=gcd(a,b)
同时我们也可以得出
b*x2+(a%b)*y1=gcd(b,a%b)
由于
gcd(a,b)=gcd(b,a%b)
那么我们可以进一步得到
a*x1+b*y1=b*x2+(a%b)*y2，a%b=(a-a/b*b)
我们要知道，这里的a/b*b并不一定会等于a，因为我们并不一定会整除
啊！进一步得到
a*x1+b*y1=a*y2+b*(x2-（a/b）*y2)
这个等式要保证恒成立，那么
x1=y2，y1=x2-(a/b)*y2
所以我们就得到了递推公式
但是我本人并不知道这个所以是如何的来的

递归的代码
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 = y;
    y = x - (a/b) * y;
    x = t;
    return r;
}
如有不理解，可以手动模拟一下，发现的确是这样，但是还是觉得自己
对这个的理解不是很清晰。

扩展欧几里得的非递归，据说是参考的这两个链接
http://anh.cs.luc.edu/331/notes/xgcd.pdf 
http://math.cmu.edu/~bkell/21110-2010s/extended-euclidean.html 

非递归的代码
int exgcd(int m, int n, int &x, int &y) {
    if (n == 0) {
        x = 1; y = 0;
        return m;
    }
    int a, a1, b, b1, c, d, q, r, t;
    a1 = b = 1;
    a = b1 = 0;
    c = m; d = n;

    q = c/d; r = c%d;
    while (r) {
        c = d;
        d = r;
        t = a1;
        a1 = a;
        a = t - q * a;
        t = b1;
        b1 = b;
        b = t - q * b;

        q = c/d;
        r = c%d;
    }
    x = a; y = b;
    return d;
}

我自己觉得这个的确很有道理，但是思维怎么就转变不过来，还是看一个题目吧！

题目链接
http://hihocoder.com/problemset/problem/1297

首先我们这道题目用的是扩展欧几里得的方法求解的
题目输入s1,s2,v1,v2以及m，其中m是石板的总数目，编号是0到
m-1，s1是第一个人的起始位置，v1是第一个人的步长，s2是第二个
人的起始位置，v2是第二个人的步长，我们假设第二个人比第一个人
多走了t圈然后和他相遇，假设第一个人走了k步，那么我们可以列出
等式
s1+v1*k=s2+v2*k-m*t，
其中我们只有k与t是变量
(v1-v2)*k+m*t=s2-s1
这样我们是不是就把这个转换为A*x+B*y=C了，我们令
A=v1-v2
B=m
C=s2-s1
首先我们要明白一点，我们要保证A>0，那么如果A<0的话，那么我
们令A+m，在这里其实就是等同于v1-v2的值加上了m，其实这就等
同于每次多跳了一圈，对最后的结果是不影响的。
首先我们由上述定理可以得出，如果我们这里的C%gcd(a,b)不为
0的话，那么这个放方程是没有解的。直接输出-1。
否则的话，我们对A*x+B*y=C等式两边都除以gcd(a,b)，那么
我们就得到了A'*x+B'*y=C'，其中
A'=A/gcd(a,b)，B'=B/gcd(a,b)，C'=C/gcd(a,b)
其中A'与B'是互质的
我们根据扩展欧几里得可以算出A'*x+B'*y=1的所有解，我们的
目标是A'*x+B'*y=C'，那么我们可以把得出的结果*C'，但是我
们尽管通过欧几里得求得x的值，我们也不能保证x的值是我们想要
的最小的正整数，所以我们构建了x的解集
A'*x+B'*y=1
A'*x+B'*y+[u+(-u)]*A'*B'=1
转化一下
A'*(x+u*B')+B'*(y-u*A')=1
那么
x=x+u*B'
y=y-u*A'
那么我们把x与y进行相同的操作同样也是这个方程的解，这就是
代码中为什么会对求得的结果进行+B'的原因

给出代码
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long ll;

ll exgcd(ll m, ll n, ll &x, ll &y) {
    if (n == 0) {
        x = 1; y = 0;
        return m;
    }
    int a, a1, b, b1, c, d, q, r, t;
    a1 = b = 1;
    a = b1 = 0;
    c = m; d = n;

    q = c/d; r = c%d;
    while (r) {
        c = d;
        d = r;
        t = a1;
        a1 = a;
        a = t - q * a;
        t = b1;
        b1 = b;
        b = t - q * b;

        q = c/d;
        r = c%d;
    }
    x = a; y = b;
    return d;
}

int main() {
    ll s1, s2, v1, v2, m;
    ll A, B, C;
    ll k, t;

    while(cin >> s1 >> s2 >> v1 >> v2 >> m){

        A = v1-v2;
        B = m;
        C = s2-s1;

        if (A < 0) A += B;
        ll d = __gcd(A, B);

        if (C % d) cout << -1 << endl;
        else {
            A = A/d;
            B = B/d;
            C = C/d;
            exgcd(A, B, k, t);
            k = (k * C) % B;
            while (k < 0) {
                k += B;
            }
            cout << k << endl;
        }
    }

    return 0;
}