拓展欧几里得算法

普通的欧几里得算法

这大概是最为人熟知的数论算法了吧，作用是求出gcd(a,b)，代码如下: 递归版本:

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

循环版本:

int gcd(int a,int b)
{
    while(b)
    {
        int tmp=a%b;
        a=b;
        b=tmp;
    }
    return a;
}

另外，求最大公约数还可以用库函数__gcd(),这是GNU里的库函数，不是标准库函数，头文件为 < algorithm >

拓展欧几里得算法

而拓展欧几里得算法则是在算出a,b的最大公约数的同时，顺带计算贝祖公式的一组解，即ax+by=gcd(a,b)的解( 贝祖等式 )

那么，如果要手算ax+by=gcd(a,b)，我们要怎么算呢?

首先，因为gcd(a,b)=gcd(b,a%b),

所以a x + b y = gcd ( a,b ) = gcd ( b,a%b )= c x' + d y' ( 其中，c = b，d = a % b )

又因为a % b = a - (a/b) * b

所以有a y' + b ( x' - (a/b) * y) = a x + b y

比较两式，可得x = y' , y = x' - ( a / b ) * y’

又因为当b = 0时，解为x = 1 , y = 0，所以只要递归求解至b = 0即可

代码如下:

int exgcd(int a,int b,int& x,int& y)
{
    if(!b)
    {
        x=1;y=0;
        return a;
    }
    int r=exgcd(b,a%b,x,y);
    int tmp=x;
    x=y;
    y=tmp-(a/b)*y;
    return r;
}

另一种写法为:

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

拓展欧几里得的一个作用，是求出满足形如ax+by=c的方程的解。根据贝祖等式可以看出，若c % gcd( a,b ) != 0，则方程无解。代码如下:

ll solve(ll a,ll b,ll b)
{
    ll x,y;
    ll g=exgcd(a,b,x,y);
    if(c%g)
        return -1;
    x*=c/gcd;
    b/=gcd;
    if(b<0)
        b=-b;        //如果b < 0，就取其绝对值
    ll ans=x%b;        //如果求出的解小于0，就将其加上b
    if(ans<=0)
        ans+=b;
    return ans; 
}

这个算法的另一个作用是求最小逆元。根据a x = 1 ( mod p) ==> ax = py + 1 ==> ax + by = 1，不难看出，求最小逆元其实就是求ax + by = c的解的特殊情况。因此，利用拓展欧几里得求最小逆元的代码如下：

int inv(int a,int b)
{
    int x,y;
    int g=exgcd(a,bx,y);
    if(1%gcd)        
        return -1;
    x*=1/g;
    b=abs(b);        //假如b是一个负数，就取b的绝对值
    int ans=x%b;    //假如求得的解不是正数，就将解对b取模再加上b
    if(ans<=0)        
        ans+=b;
    return ans;
}

一道求ax + by = c 的模板题：青蛙约会

我的代码:

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
using namespace std;
typedef long long ll;
void exgcd(ll a,ll b,ll& d,ll& x,ll& y)
{
    if(!b){
        d=a;x=1;y=0;
    }
    else{
        exgcd(b,a%b,d,y,x);
        y-=(a/b)*x;
    }
}
ll solve(ll a,ll b,ll c)
{
    ll x;ll y;ll g;
    exgcd(a,b,g,x,y);
    if(c%g)        return -1;
    x=(x*c)/g;
    b/=g;
    if(b<0)        b=-b;
    ll ans=x%b;
    if(ans<=0)        ans+=b;
    return ans;
}
int main()
{
    ll x,y,m,n,L;
    scanf("%lld %lld %lld %lld %lld",&x,&y,&m,&n,&L);
    ll ans=solve(n-m,L,x-y);
    if(ans==-1)
        printf("Impossible\n");
    else
        printf("%lld\n",ans);
    return 0;
}

~~为什么求ax + by = c的代码可以这样写呢？我也不知道啊OTZ(发出了数学太差的声音)~~

NeoRuTayE's blog