扩展欧几里得

1、GCD

 　　辗转相除法：

　　

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

2、扩展欧几里得算法

　　由贝祖等式可知ax+by=gcd(a,b) 。（由此可推一般式： ax+by=d）

　　如果要解x , y 。由数学推导可知：

　　1、当b=0 时可知x=1，y=0。

　　2、当b！=0时，a*x1+b*y1=gcd(a,b)，  b*x2+(a%b)*y2=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 。故带入后，由恒等可知 x1=y2  ； y1= x2-(a/b)*y2 。

所以扩展欧几里得算法代码可以写为

 

long long extgcd(long long a, long long b, long long &x, long long &y)
{
    long long gcd, t;
    if (b == 0) {
         x = 1;
         y = 0;
         return a;
    }
    gcd = extgcd(b, a % b, x, y);
    t = x - a/b*y; 
    x = y; 
    y = t;
    return gcd;
}        

 

 

扩展欧几里得算法的题目有时候会由于欧几里得算法解的无穷性而变得难以嵌入。（本次是以poj1061（青蛙的约会）为样例讨论  ）

ax + by = c其实就等价于ax ≡ c (mod b)，当存在一个特解x*使得ax*+by*=c。那么x*加上若干倍b还是这个方程的解，即x=x*+k*b仍为方程的解。

因而方程在[0, b-1]上一定有整数解（假如小于0，你加上若干倍b啊，就可以让它保持在0~b-1中；如果大于b-1，你减去若干倍b啊，它也保持0~b-1）。

那么怎样求最小的非负整数x呢？

1、若gcd(a, b) == 1，则方程ax ≡ c (mod b)在[0, b-1]上有唯一解。

证明如下：设有两个方程 a*x1+b*y1=c   ，   a*x2+b*y2=c 。两式相减得  a*(x1-x2)=b*(y2-y1) 。因为 gcd(a,b) ==1 。故x1,x2 在 [0, b-1] 区间内当且仅当 x1== x2 。原命题得证。

2、若gcd(a, b) ==d，则方程ax ≡ c (mod b)在[0, b/d - 1]上有唯一解。

 

证明：ax + by = c，如果有解，两边同除以d，就有a/d * x + b/d * y = c/d，即a/d * x ≡ c/d (mod b/d)，显然gcd(a/d, b/d) = 1，所以由定理二知道x在[0, b/d - 1]上有唯一解。

故可知poj1061题 在判断c%gcd ==0 后我们需要讨论：x 是否为正整数的最小解！

由上述两个结论我们可以知道只需要将所解的 x 拖回 [0,b/gcd-1] 当中就可以求得最小正整数值。即下面代码中的第40~44行。

（思路来源：https://www.cnblogs.com/comeon4mydream/archive/2011/07/18/2109060.html）

 

#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#define ll long long int
using namespace std;
long long exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    ll r=exgcd(b,a%b,y,x);
    y-=x*(a/b);
    return r;
}
int main()
{
    ll x,y,m,n,L;
    ll a,b,c,gcd;
    while(scanf("%lld%lld%lld%lld%lld",&x,&y,&m,&n,&L)!=EOF)
    {
        a=m-n;
        b=L;
        c=y-x;
        if(a<0)
        {
            a=-a;
            c=-c;
        }
        gcd=exgcd(a,b,x,y);
        if(c%gcd!=0)
            printf("Impossible");
        else
        {
            x=x*c/gcd;
            int t=b/gcd;
            if(x>=0)
                x=x%t;
            else
                x=x%t+t;
            printf("%lld\n",x);
        }
    }
    return 0;
}