一. 算法描述
欧几里得算法
我们先回忆一下欧几里得算法(辗转相除法):
这个很好证明:
首先,
,
。
1.现假设
,那么
,取
,这样
,所以a和b的约数是b和a mod b的约数。
2。现假设
,那么
,取
,这样
,所以b和 a mod b 的约数是a和b的约数。
因此两者的最大公约数是一样的。
扩展欧几里得算法
裴蜀定理
对于任意的正整数a,b,一定存在整数x,y,使得:
此处要讨论的即为,对于给定的a,b,如何求出这里的x,y?
由欧几里得算法可得:
由裴蜀定理可得:
由于是对任意的a与b都要成立,那么:
也就是说,我们想求x与y,我们可以先求
中的x‘和y’,假设我们求到了x‘和y’,我们就可以利用(*)式求出x和y。
由次我们发现,这里是一种
递归
的思想。那我们就要找到递归结束的条件。
对于
,我们可以发现,
最后b一定会为0
,也就是
,此时很容易求出
,这样,我们可以根据这一组x和y的值,不断返回前面所有组的(x,y),最后求出给顶我们的正整数a,b的那组解。
二.代码实现
由于我们要返回一个二元组(x,y),这样我们定义pair就可以。
那对于代码实现,就很简单了,在gcb函数中,我们先定义递归结束的条件:
if(b==0){
return {0,1};
}
否则,我们就去求
和
的对应的
:
II tmp=gcb(b,a%b);
我们再用(*)式,返回我输入的
对应的
,完整代码:
#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;
}
}
三.扩展欧几里得算法应用
线性同余方程
对于给定的
,求出一个
,使得
.
那我们考虑用扩展欧几里得算法求出这个
:
首先,
可以写成
这样,
现令
,若
不是
的倍数,那么此时不存在
,这是因为,
都是
的倍数,这样
一定是
的倍数。
所以上式我们可以写成:
这样,我们只需要把求到的
扩大
倍即可。
代码实现
这里代码实现和扩展欧几里得算法基本一样,但是需要注意一个问题是:经过exgcb函数传回来的
是int类型,
也是int类型,在最后得到结果时,如果直接写:
,在数字非常大的时候,两个int相乘是有可能爆掉的!!!因而我们修改成:
。但是题目有要求说,一定要输出一个int范围内的答案,但事实上我们观察
,
是可以对
取余的
!因而我们最终输出写成:
具体代码:
#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?
答案是显然的
。扩展欧几里得算法就是求出
修改后的代码为:
#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平台