扩展欧几里得算法及其相关应用

最新推荐文章于 2023-09-08 16:12:16 发布

猪猪奋斗记

最新推荐文章于 2023-09-08 16:12:16 发布

阅读量943

点赞数 2

分类专栏：数学

本文链接：https://blog.csdn.net/bigbigship/article/details/25073451

版权

数学专栏收录该内容

119 篇文章 0 订阅

订阅专栏

欧几里得算法：

欧几里得算法又称辗转相除法，用于求两个数的最大公约数

设 a=kb+r; 则gcd(a,b)=gcd(b,r)即gcd(a,b)=gcd(b,a%b);

a可以表示成a = kb + r，则r = a mod b

假设d是a,b的一个公约数，则有

d|a, d|b，而r = a - kb，因此d|r

因此d是(b,a mod b)的公约数

假设d 是(b,a mod b)的公约数，则

d | b , d |r ，但是a = kb +r

因此d也是(a,b)的公约数

因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证

实现代码：

1）递归：

LL gcd(LL a,LL b)
{
    if(b)
        return gcd(b,a%b);
    return a;
}

2）非递归

LL gcd( LL a,LL b)
{
    while(b!=0)
    {
        LL r=b;
        b=a%b;
        a=r;
    }
    return a;
}

扩展欧几里得算法：

基本算法：对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然存在整数对 x，y ，使得 gcd（a，b）=ax+by。

证明：设 a>b。

　　1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；

　　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;

这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.

　上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

代码：

void ex_gcd(LL a,LL b,LL &x,LL &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return ;
    }
    else ex_gcd(b,a%b,x,y);
    LL ans=x;
    x=y;
    y=ans-a/b*y;
}

题目链接：http://acm.nefu.edu.cn/JudgeOnline/problemshow.php?problem_id=630

首先列出方程 ax+by=1;

然后分析只有当a,b互质的时候才能得到1

因此经典扩展欧几里得算法

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;

typedef long long LL;

void ex_gcd(LL a,LL b,LL &x,LL &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return ;
    }
    else ex_gcd(b,a%b,x,y);
    LL ans=x;
    x=y;
    y=ans-a/b*y;
}
LL gcd(LL a,LL b)
{
    if(b)
        return gcd(b,a%b);
    return a;
}
int main()
{
    int T;
    cin>>T;
    while(T--)
    {
        LL a,b;
        cin>>a>>b;
        if(a < b) swap(a,b);
        LL g = gcd(a,b);
        if(g != 1)
        {
            puts("-1");
            continue;
        }
        if(b == 1 && a == 2)
        {
            puts("1");
            continue;
        }
        if(b == 1)
        {
            puts("2");
            continue;
        }
        if(b == 0 && a != 1)
        {
            puts("-1");
            continue;
        }
        if(b == 0 && a == 1)
        {
            puts("1");
            continue;
        }
        LL x,y;
        ex_gcd(a,b,x,y);
        x=abs(x);
        y=abs(y);
        cout<<x+y-1<<endl;
    }
    return 0;
}

扩展欧几里得的几个应用

1.求解不定方程

给定a,b,c 求a*x+b*y=c的一组解

step 1:应用扩展欧几里得得出（a, b）= t

step 2:判断。如果不满足 t | c 那么无解。否则得到一组解 x1 = x * c / t ,x2 = y * c / t;

其他的解满足 x = x1+ b / t * n;

y = y1 - a / t * n;(n为任意整数)

代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;

typedef long long LL;

LL ex_gcd(LL a,LL b,LL &x,LL &y){
    if(!b){
        x=1,y=0;
        return a;
    }
    else{
        int gcd = ex_gcd(b,a%b,x,y);
        int tmp = x;
        x = y;
        y = tmp - a/b*x;
        return gcd;
    }
}

int main()
{
    LL a,b,c,x,y;
    while(~scanf("%lld%lld%lld",&a,&b,&c)){
        LL t = ex_gcd(a,b,x,y);
        if(c % t) puts("No solution!");
        else{
            printf("ans: %lld %lld\n",c/t*x,c/t*y);
        }
    }
    return 0;
}

2.解线性同余数方程组

1）解一元线性同余方程

同余方程 ax≡b (mod n)对于未知数 x 有解，当且仅当 gcd(a,n) | b。且方程有解时，方程有 gcd(a,n) 个解。

求解方程 ax≡b (mod n) 相当于求解方程 ax+ ny= b, (x, y为整数)

然后得出一个解为( x * b / d ) % n, 全部解为 xi = ( x * b / d +i * (n / d) ) % n (o<=i<=d-1)

一个结论：设 A = x * b / d , B =n / d 最小的解 min = ( A % S + S ) % S;

代码如下：

#include <iostream>

using namespace std;

typedef long long LL;

LL ex_gcd(LL a,LL b,LL &x,LL &y){
    if(!b){
        x=1,y=0;
        return a;
    }
    else{
        LL gcd=ex_gcd(b,a%b,x,y);
        LL tmp = x;
        x = y;
        y = tmp - a / b * x;
        return gcd;
    }
}

int main()
{
    LL a,n,b,x,y;
    while(cin>>a>>n>>b){ //求ax ≡ b (mod n)
        LL d = ex_gcd(a,n,x,y);
        if(b%d) cout<<"No solution"<<endl;
        else{
            LL s = n / d;
            cout<<((x*b/d)%s + s)%s<<endl;
        }
    }
    return 0;
}

2）解一元线性同余方程组

理论请见：http://wenku.baidu.com/link?url=Wgvl2krKc6x5Yz8HVx4x2BHVPrPB0F0ByxeudHSLA_whQXMa9nCCPF6TM-UPoO8pZzchl2Pm0x-vmQj1ycLqbET5UOndFLr3x3lTke5BwVi

例如： x≡b1 (mod a1)

x≡b2 (mod a2)

b1+k1*a1 = b2 + k2*a2 然后可以求出k2 ,然后两两一求。

例题： POJ 2891 Strange Way to Express Integers

就是求一个同余方程组的最小解

代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;

LL ex_gcd(LL a,LL b,LL &x,LL &y){
    if(!b){
        x=1,y=0;
        return a;
    }
    else{
        LL gcd = ex_gcd(b,a%b,x,y);
        LL tmp = x;
        x = y;
        y = tmp-a/b*y;
        return gcd;
    }
}

int n;
LL r1,r2,a1,a2,x,y;

LL solve(){
    LL r1,r2,a1,a2,x,y;
    bool tag = 1;
    cin>>a1>>r1;
    for(int i=1;i<n;i++){
        cin>>a2>>r2;
        LL a = a1,b=a2,c = r2-r1;
        LL d = ex_gcd(a,b,x,y);
        if(c%d) tag = 0;
        LL t = b/d;
        x = (x*(c/d)%t+t)%t;
        r1 = a1*x+r1;
        a1 = a1*(a2/d);
    }
    if(tag) return r1;
    else return -1;
}

int main()
{
    while(~scanf("%d",&n)){
        printf("%lld\n",solve());
    }
    return 0;
}

3) 求乘法逆元

对于整数 a,p 如果存在一个x使得 ax ≡ 1 (mod p) 则说 x 是 a 的模 p 乘法逆元

同余方程ax≡b (mod n)，如果 gcd(a,n)== 1，则方程只有唯一解。

在这种情况下，如果 b== 1，同余方程就是 ax=1 (mod n ),gcd(a,n)= 1。

这时称求出的 x 为 a 的对模 n 乘法的逆元。

对于同余方程 ax= 1(mod n )， gcd(a,n)= 1 的求解就是求解方程

ax+ ny= 1，x, y 为整数。这个可用扩展欧几里德算法求出，原同余方程的唯一解就是用扩展欧几里德算法得出的 x 。

定理：a存在模p的乘法逆元的充要条件是gcd(a,p) = 1

代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;

typedef long long LL;

LL ex_gcd(LL a,LL b,LL &x,LL &y){
    if(!b){
        x=1,y=0;
        return a;
    }
    else{
        int gcd = ex_gcd(b,a%b,x,y);
        int tmp = x;
        x = y;
        y = tmp - a/b*x;
        return gcd;
    }
}

int main()
{
    LL p = 13,x,y;
    for(LL i=1;i<13;i++){
        ex_gcd(i,p,x,y);
        while(x<0) x+=p;
        printf("the inv of %lld mod %lld is %lld\n",i,p,x);
    }
    return 0;
}