数学知识——扩展欧几里得算法

  • Post author:
  • Post category:其他


一. 算法描述

欧几里得算法

我们先回忆一下欧几里得算法(辗转相除法):

gcb(a,b)=gcb(b,a\, mod\, b)

这个很好证明:

首先,
d\mid a,d\mid b\rightarrow d\mid ax+by

a\, mod \, b=a-\left \lfloor \frac{a}{b} \right \rfloor b=a-cb

1.现假设
gcb(a,b)=x
,那么
x\mid a,x\mid b
,取
x=1,y=-c
,这样
x\mid a-cb,x\mid b
,所以a和b的约数是b和a mod b的约数。

2。现假设
gcb(b,a\, mod\, b)=x
,那么
x\mid b,x\mid a-cb
,取
x=c,y=1
,这样
x\mid a,x\mid b
,所以b和 a mod b 的约数是a和b的约数。

因此两者的最大公约数是一样的。

扩展欧几里得算法

裴蜀定理

对于任意的正整数a,b,一定存在整数x,y,使得:

ax+by=gcb(a,b)

此处要讨论的即为,对于给定的a,b,如何求出这里的x,y?

由欧几里得算法可得:

ax+by=gcb(a,b)=gcb(b,a\, mod \, b)=gcb(b,a-\left \lfloor \frac{a}{b} \right \rfloor b)

由裴蜀定理可得:

bx'+(a-\left \lfloor \frac{a}{b} \right \rfloor b)y'=gcb(b,a-\left \lfloor \frac{a}{b} \right \rfloor b)=ay'+b(x'-\left \lfloor \frac{a}{b} \right \rfloor y')

由于是对任意的a与b都要成立,那么:

x=y',y=x'-\left \lfloor \frac{a}{b} \right \rfloor y'\cdot \cdot \cdot\cdot \cdot \cdot (*)

也就是说,我们想求x与y,我们可以先求
bx'+(a-\left \lfloor \frac{a}{b} \right \rfloor b)y'=gcb(b,a-\left \lfloor \frac{a}{b} \right \rfloor b)
中的x‘和y’,假设我们求到了x‘和y’,我们就可以利用(*)式求出x和y。

由次我们发现,这里是一种

递归

的思想。那我们就要找到递归结束的条件。

对于
gcb(a,b)=gcb(b,a%b)
,我们可以发现,

最后b一定会为0

,也就是
gcb(a,0)
,此时很容易求出
x=1,y=0
,这样,我们可以根据这一组x和y的值,不断返回前面所有组的(x,y),最后求出给顶我们的正整数a,b的那组解。

二.代码实现

由于我们要返回一个二元组(x,y),这样我们定义pair就可以。

那对于代码实现,就很简单了,在gcb函数中,我们先定义递归结束的条件:

if(b==0){
    return {0,1};
}

否则,我们就去求
b

a%b
的对应的
x',y'
:

II tmp=gcb(b,a%b);

我们再用(*)式,返回我输入的
a,b
对应的
x,y
,完整代码:

#include<iostream>
using namespace std;
typedef pair<int,int> II;
int n;
II gcb(int a,int b){
    if(b==0){
        return {1,0};
    }
    else{
        II tmp=gcb(b,a%b);
        return {tmp.second,tmp.first-a/b*(tmp.second)};
    }
}
int main(){
    cin>>n;
    while(n--){
        int a,b;
        cin>>a>>b;
        II res=gcb(a,b);
        cout<<res.first<<" "<<res.second<<endl;
    }
}

三.扩展欧几里得算法应用

线性同余方程

对于给定的
a,b,m
,求出一个
x
,使得
ax\equiv b\left ( mod\, \, m \right )
.

那我们考虑用扩展欧几里得算法求出这个
x

首先,
ax\equiv b\left ( mod\, \, m \right )
可以写成

ax=my+b

这样,
ax-my=b\rightarrow ax+my=b

现令
d=gcb(a,m)
,若
b
不是
d
的倍数,那么此时不存在
x
,这是因为,
ax,my
都是
d
的倍数,这样
b
一定是
d
的倍数。

所以上式我们可以写成:
\frac{d}{b}(ax+my)=d\rightarrow ax'+my'=d

这样,我们只需要把求到的
x'
扩大
\frac{b}{d}
倍即可。

代码实现

这里代码实现和扩展欧几里得算法基本一样,但是需要注意一个问题是:经过exgcb函数传回来的
res.first
是int类型,
\frac{b}{d}
也是int类型,在最后得到结果时,如果直接写:
int \, x=res.first*\frac{b}{d}
,在数字非常大的时候,两个int相乘是有可能爆掉的!!!因而我们修改成:
long\,\, long \, x=res.first*\frac{b}{d}
。但是题目有要求说,一定要输出一个int范围内的答案,但事实上我们观察
ax\equiv b\left ( mod\, \, m \right )


x
是可以对
m
取余的

!因而我们最终输出写成:

int \, x=((long\, \, \, \, long)res.first*\frac{b}{d})%m

具体代码:

#include<iostream>
using namespace std;
typedef pair<int,int> II;
int gcb(int a,int b){
    if(b==0){
        return a;
    }
    else{
        gcb(b,a%b);
    }
}
II exgcb(int a,int b){
    if(b==0) return {1,0};
    else{
        II tmp=exgcb(b,a%b);
        return {tmp.second,tmp.first-a/b*(tmp.second)};
    }
}
int main(){
    int n;
    scanf("%d",&n);
    while(n--){
        int a,b,m;
        scanf("%d%d%d",&a,&b,&m);
        int d=gcb(a,m);
        if(b%d!=0) puts("impossible");
        else{
            II res=exgcb(a,m);
            int x=((long long)b/d*res.first)%m;
            printf("%d\n",x);
        }
    }
}

代码优化

那我们现在在想另外一个问题,函数exgcb能不能代替gcb?


答案是显然的

。扩展欧几里得算法就是求出
xa+by=gcb(a,b)

d=res.first*a+res.second*b

修改后的代码为:

#include<iostream>
using namespace std;
typedef pair<int,int> II;
II exgcb(int a,int b){
    if(b==0) return {1,0};
    else{
        II tmp=exgcb(b,a%b);
        return {tmp.second,tmp.first-a/b*(tmp.second)};
    }
}
int main(){
    int n;
    scanf("%d",&n);
    while(n--){
        int a,b,m;
        scanf("%d%d%d",&a,&b,&m);
        II tmp=exgcb(a,m);
        int d=a*tmp.first+m*tmp.second;
        if(b%d!=0) puts("impossible");
        else{
            II res=exgcb(a,m);
            int x=((long long)b/d*res.first)%m;
            printf("%d\n",x);
        }
    }
}

感谢AcWing平台



版权声明:本文为qq_64637770原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。