欧几里得算法&扩展欧几里得

本文包括：
1.几里德算法 递归实现
2.扩展欧几里德算法 递归实现
…实际上就是把别人总结的，我认为有助于自己理解的内容copy过来，再加上几句自己的理解。

欧几里得算法

欧几里德算法又称辗转相除法，用于计算两个整数a,b的最大公约数 gcd(a,b)。基本算法：设 a = qb + r，其中a，b，q，r都是整数，则 gcd(a,b) = gcd(b,r)，即 gcd(a,b) = gcd(b,a%b)。

证明：
a = qb + r
如果 r = 0，那么 a 是 b 的倍数，此时显然 b 是 a 和 b 的最大公约数。
如果 r ≠ 0，任何整除 a 和 b 的数必定整除 a - qb = r，而且任何同时整除 b 和 r 的数必定整除 qb + r = a，所以 a 和 b 的公约数集合与 b 和r 的公约数集合是相同的。特别的，a 和 b 的最大公约数是相同的。

递归实现：

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

扩展欧几里得算法

有两个数 a,ba,b，对它们进行辗转相除法，可得它们的最大公约数——这是众所周知的。然后，收集辗转相除法中产生的式子，倒回去，可以得到 ax+by=gcd(a,b)ax+by=gcd(a,b)的整数解。

先来看下这个几乎所有总结扩展欧几里得算法的帖子中都会用到的例子

用类似辗转相除法，求二元一次不定方程47x+30y=1的整数解。
47=30*1+17
30=17*1+13
17=13*1+4
13=4*3+1
然后把它们改写成“余数等于”的形式

17=47*1+30*(-1) //式1
13=30*1+17*(-1) //式2
4=17*1+13*(-1) //式3
1=13*1+4*(-3)
然后把它们“倒回去”

1=13*1+4*(-3) //应用式3
1=13*1+[17*1+13*(-1)]*(-3)
1=13*4+17*(-3) //应用式2
1=[30*1+17*(-1)]*4+17*(-3)
1=30*4+17*(-7) //应用式1
1=30*4+[47*1+30*(-1)]*(-7)
1=30*11+47*(-7)
得解x=-7, y=11。

这些式子自己手写一遍才能理解的更清楚。
这个“倒回去”的过程对应的就是代码中递归向上返回的过程。现在的问题就是要找出这个可以倒回去的递推关系。

   设：a>b。
　　推理1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0;//推理1
　　推理2，ab!=0 时
　　设 ax1+by1=gcd(a,b);
　　bx2+(a mod b)y2=gcd(b,a mod b);
　　根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);
　　则:ax1+by1=bx2+(a mod b)y2;
　　即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;
　　根据恒等定理得：x1=y2 ,y1=x2-(a/b)*y2;//推理2
　　这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.
   上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

这样我们就找到了递推关系：

x1=y2;x1=y2;
y1=x2−(a/b)∗y2;y1=x2−(a/b)∗y2;

递推的终止条件再上面也已经给出：

当b=0，gcd（a，b）=a。
此时x=1，y=0;
当b=0，gcd（a，b）=a。
此时x=1，y=0;
由此我们可以得到递归的扩展欧几里得算法的代码：

int exgcd(int a, int b, int &x, int &y)
{
    if(b == 0)
    {//推理1，终止条件
        x = 1;
        y = 0;
        return a;
    }
    int r = exgcd(b, a%b, x, y);
    //先得到更底层的x2,y2,再根据计算好的x2,y2计算x1，y1。
    //推理2，递推关系
    int t = y;
    y = x - (a/b) * y;
    x = t;
    return r;
}

比如上面的 47x+30y=1 的例子，exgcd() 要求出 gcd(47, 30)，同时得到一组 （x，y）的解。

 用一个方便我自己理解的方式，把求解的过程写出来好了。这里的大括号"{}"里面不是函数体，"{"表示一层递归调用的开始,"}"表示该层递归的结束。 
 exgcd(47, 30, x, y)
 {
     r = exgcd(30, 17，x, y)
     {
         r = exgcd(17, 13, x, y)
         {
             r = exgcd(13, 4, x, y)
             {
                 r = exgcd(4, 1, x, y)
                 {
                     r = exgcd(1, 0, x, y)
                     {
                         x = 1; 
                         y = 0;
                         return 1;
                     }
                     t = y = 0;
                     y = x - (4/1) * y = 1;
                     x = t = 0;
                     return r = 1; 
                 }
                 t = 1;
                 y = 0 - (13/4) * 1 = -3;
                 x = 1;
                 return 1;
             }
             t = -3;
             y = 1 - (17/13) * (-3) = 4;
             x = -3;
             return 1;
         }
         t = 4;
         y = -3 - (30/17) * 4 = -7;
         x = 4;
         return 1;
     }
     t = -7;
     y = 4 - (47/30) * (-7) = 11;
     x = -7;
     return 1;
 }
 最后的结果:
 r = exgcd(47,30,x,y) = 1;
 x = -7;
 y = 11;

当gcd(a , m) != 1 的时候是没有解的这也是 ax + by = c 有解的充要条件： c % gcd(a , b) == 0

接着乘法逆元讲，一般，我们能够找到无数组解满足条件，但是一般是让你求解出最小的那组解，怎么做？我们求解出来了一个特殊的解 x0 那么，我们用 x0 % m其实就得到了最小的解了。为什么？

可以这样思考：

x 的通解不是 x0 + m*t 吗？

那么，也就是说， a 关于 m 的逆元是一个关于 m 同余的，那么根据最小整数原理，一定存在一个最小的正整数，它是 a 关于m 的逆元，而最小的肯定是在（0 , m）之间的，而且只有一个，这就好解释了。

可能有人注意到了，这里，我写通解的时候并不是 x0 + (m/gcd)*t ，但是想想一下就明白了，gcd = 1，所以写了跟没写是一样的，但是，由于问题的特殊性，有时候我们得到的特解 x0 是一个负数，还有的时候我们的 m 也是一个负数这怎么办？

当 m 是负数的时候，我们取 m 的绝对值就行了，当 x0 是负数的时候，他模上 m 的结果仍然是负数（在计算机计算的结果上是这样的，虽然定义的时候不是这样的），这时候，我们仍然让 x0 对abs(m) 取模，然后结果再加上abs(m) 就行了，于是，我们不难写出下面的代码求解一个数 a 对于另一个数 m 的乘法逆元：
便有了如下代码：

ll cal(ll a,ll b,ll c){
    ll x,y;
    ll gcd=exgcd(a,b,x,y);
    if(c%gcd!=0) return -1; //a*x+b*y=c 有解的充要条件：c%gcd(a,b)==0
    x*=c/gcd;    //a*x+b*y=gcd
    b/=gcd;
    if(b<0) b=-b; //得mod正数
    ll ans=x%b;   //求最小的解
    if(ans<=0) ans+=b;  //把解换成正值
    return ans;
}

例题1：http://acm.hdu.edu.cn/showproblem.php?pid=2669
裸的模拟 Xa + Yb = 1.求X，Y。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <stdio.h>
using namespace std;
#define ll long long

ll exgcd(ll a,ll b,ll &x,ll &y) {
 if(b==0){
   x=1;
   y=0;
   return a;
 }
  ll ans=exgcd(b,a%b,x,y);
  ll temp=x;
  x=y;
  y=temp-a/b*y;
  return ans;
}

ll cal(ll a,ll b,ll c){
   ll x,y;
   ll gcd=exgcd(a,b,x,y);
   if(c%gcd!=0) return -1; //a*x+b*y=c 有解的充要条件：c%gcd(a,b)==0
   x*=c/gcd;    //a*x+b*y=gcd
   b/=gcd;
   if(b<0) b=-b; //得mod正数
   ll ans=x%b;   //求最小的解
   if(ans<=0) ans+=b;  //把解换成正值
   return ans;
}

int main(){
 ll a,b;
 while(scanf("%lld %lld",&a,&b)!=EOF){
   ll ans=cal(a,b,1);
   if(ans==-1) printf("sorry\n");
   else printf("%lld %lld\n",ans,(1-ans*a)/b);
 }
 return 0;
}

例题二：
http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=4712
求最小逆元
乘法逆元 ax ≡1(mod m)
这里，我们称 x 是 a 关于 m 的乘法逆元
可以等价于这样的表达式： ax + m*y = 1

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
 
#define INF 0x7fffffff
#define EPS 1e-12
#define MOD 1000000007
#define PI 3.141592653579798
#define N 100000
 
using namespace std;
 
typedef long long LL;
typedef double DB;
 
LL e_gcd(LL a,LL b,LL &x,LL &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    LL ans=e_gcd(b,a%b,x,y);
    LL temp=x;
    x=y;
    y=temp-a/b*y;
    return ans;
}
 
LL cal(LL a,LL b,LL c){
    LL x,y;
    LL gcd=e_gcd(a,b,x,y);
    if(c%gcd!=0) return -1;
    x*=c/gcd;
    b/=gcd;
    if(b<0) b=-b;
    LL ans=x%b;
    if(ans<=0) ans+=b;
    return ans;
}
 
int main(){
    LL a,b,t;
    scanf("%lld",&t);
    while(t--){
        scanf("%lld%lld",&a,&b);
        LL ans=cal(a,b,1);
        if(ans==-1) printf("Not Exist\n");
        else printf("%lld\n",ans);
    }
    return 0;
}

例题三：
http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3593
求最小的步数，处理特殊一点就过去了
a * x+b * y=A-B;
x,y必须大于0，min(x,y)=C;
ans=C+abs(x-y);

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <algorithm>
 
#define INF 0x7fffffff
#define EPS 1e-12
#define MOD 100000007
#define PI 3.14159265357979823846
#define N 100005
 
using namespace std;
 
typedef long long LL;
 
LL e_gcd(LL a,LL b,LL &x,LL &y){
    if(b==0) {
        x=1;
        y=0;
        return a;
    }
    LL ans=e_gcd(b,a%b,x,y);
    LL temp=x;
    x=y;
    y=temp-a/b*y;
    return ans;
}
LL cal(LL a,LL b,LL L){
    LL x,y;
    LL gcd=e_gcd(a,b,x,y);
    if(L%gcd!=0) return -1;
    x*=L/gcd;
    y*=L/gcd;
    a/=gcd;
    b/=gcd;
    LL ans=((LL)INF)*((LL)INF), f;
    LL mid=(y-x)/(a+b);
    for(LL T=mid-1;T<=mid+1;T++){
        if(abs(x+b*T)+abs(y-a*T)==abs(x+b*T+y-a*T))
            f=max(abs(x+b*T),abs(y-a*T));
        else
            f=fabs(x-y+(a+b)*T);
        ans=min(ans,f);
    }
    return ans;
}
 
int main()
{
    //freopen("in.in","r",stdin);
    //freopen("out.out","w",stdout);
    LL A,B,a,b,x,y;
    int t; scanf("%d",&t);
    while(t--){
        scanf("%lld%lld%lld%lld",&A,&B,&a,&b);
        LL L=B-A;
        LL ans=cal(a,b,L);
        if(ans==-1) printf("-1\n");
        else printf("%lld\n",ans);
    }
    return 0;
}

例题四：
http://poj.org/problem?id=1061
青蛙的约会，裸的扩展欧几里得
假设同一起点出发：
L1=(m-n) * x1;
L2=L * y1;
L1-L2=y-x;
(m-n) * x1+L * y1=y-x;

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
 
#define INF 0x7fffffff
#define EPS 1e-12
#define MOD 1000000007
#define PI 3.141592653579798
#define N 100000
 
using namespace std;
 
typedef long long LL;
typedef double DB;
 
LL e_gcd(LL a,LL b,LL &x,LL &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    LL ans=e_gcd(b,a%b,x,y);
    LL temp=x;
    x=y;
    y=temp-a/b*y;
    return ans;
}
 
LL cal(LL a,LL b,LL c){
    LL x,y;
    LL gcd=e_gcd(a,b,x,y);
    if(c%gcd!=0) return -1;
    x*=c/gcd;
    b/=gcd;
    if(b<0) b=-b;
    LL ans=x%b;
    if(ans<=0) ans+=b;
    return ans;
}
 
int main(){
    LL x,y,m,n,L;
    while(scanf("%lld%lld%lld%lld%lld",&x,&y,&m,&n,&L)!=EOF){
        LL ans=cal(m-n,L,y-x);
        if(ans==-1) printf("Impossible\n");
        else printf("%lld\n",ans);
    }
    return 0;
}

例题五：
http://acm.hdu.edu.cn/showproblem.php?pid=1576
做点处理即可。
(A/B)%9973，
A=x * B
A= n+y * 9973
x * B + y * 9973 = n

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
 
#define INF 0x7fffffff
#define EPS 1e-12
#define MOD 1000000007
#define PI 3.141592653579798
#define N 100000
 
using namespace std;
 
typedef long long LL;
typedef double DB;
 
LL e_gcd(LL a,LL b,LL &x,LL &y){
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    LL ans=e_gcd(b,a%b,x,y);
    LL temp=x;
    x=y;
    y=temp-a/b*y;
    return ans;
}
 
LL cal(LL a,LL b,LL c){
    LL x,y;
    LL gcd=e_gcd(a,b,x,y);
    if(c%gcd!=0) return -1;
    x*=c/gcd;
    b/=gcd;
    if(b<0) b=-b;
    LL ans=x%b;
    if(ans<=0) ans+=b;
    return ans;
}
 
int main(){
    LL n,b,t;
    scanf("%I64d",&t);
    while(t--){
        scanf("%I64d%I64d",&n,&b);
        LL ans=cal(b,9973,n);
        if(ans==-1) printf("Impossible\n");
        else printf("%lld\n",ans);
    }
    return 0;
}

暂时就写这么多啦，等以后遇到后再进行补充~