扩展欧几里德算法是用来在已知a, b求解一组x,y,使它们满足贝祖等式: ax+by = gcd(a, b) =d(解一定存在,根据数论中的相关定理)。扩展欧几里德常用在求解模线性方程及方程组中。
我们先来看欧几里德算法(即gcd算法)的一个重要性质:gcd(a,b)=gcd(b,a%b)。这一性质的证明可以见百度百科,在此不再赘述。由此可以写出gcd算法的递归形式:
int gcd(int a,int b)
{
return b ? gcd(b,a%b) : a;
}
只有一行语句,十分简洁。
而扩展欧几里德算法是为了求方程组ax+by=gcd(a,b)中的x和y,且x和y为整数。它是欧几里德算法的扩展,我们先来看代码:
int exgcd(int a,int b,int &x,int &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
int r=exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
return r;
}
第一眼看上去有些难以理解,我第一眼看上去总是被这个递归搞混。。。而且对于x和y为什么在递归过程中要这样变化弄得不知所云。不要紧,我们先来看解析:
对于以上代码,因为在每一次gcd的递归中,我们令现阶段的a`=b,b`=a mod b,其中a和b是上一阶段的数值。对于现阶段我们的问题转化成了:求x,y使得a`x+b`y=gcd(a`,b`)=gcd(a,b)。由于b`=a mod b=a-a/b*b,我们带入式中有:a`x+b`y=gcd(a,b)--->bx+(a-a/b*b)y=gcd(a,b)--->ay+b(x-a/b*y)=gcd(a,b)
比较始末式子,观察到转换过程为:
- x=y;
- y=x-a/b*y
这也是我们程序里后面几句代码的意思。
那么为什么b==0时要把x设为1,y设为0呢?我们用一个实例来看这个函数递归的过程:
#include<stdio.h>
#include<string.h>
#include<iostream>
#include<algorithm>
using namespace std;
int exgcd(int a,int b,int &x,int &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
int r=exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
cout<<"r:"<<r<<endl;
cout<<"x:"<<x<<endl;
cout<<"y:"<<y<<endl;
cout<<"a:"<<a<<endl;
cout<<"b:"<<b<<endl;
cout<<"a/b:"<<a/b<<endl;
system("pause");
return r;
}
int main()
{
int a,b,x,y;
while(cin>>a>>b)
{
x=y=-1;
int r=exgcd(a,b,x,y);
}
return 0;
}
我们输入a,b的值,来看递归中exgcd里的值怎么变化的。我取a=36,b=21
一开始我们递归到b=0时,此时值因为程序写的问题没有显示出来,不过可以推断出:
r:3
x:1
y:0
a:3
b:0
我们发现此时x,y,a,b满足方程ax+by=gcd(a,b)=r=3。也就是说,在初始状态,因为a=r,b=0,所以满足这个方程x只能取1,y取0,也就是代码里赋值的含义。
接着递归,我们有:
此时是初始状态的上一层递归状态,其中x和y的值由语句
int t=x;
x=y;
y=t-a/b*y;
变化而来,仍然满足方程ax+by=gcd(a,b)=r=3。
后面几个递归都一样,一直到最后一层,我们计算一下:36*3-5*21=3,嗯,得到了我们要求的x和y。
现在我们知道了exgcd函数的计算流程:先由
int r=exgcd(b,a%b,x,y);
递归到最后一层,然后在每一层返回中逐步更新当前状态x,y的值,最终得到解答。
说了这么多,不如来做个练习吧?
练习题:http://codeforces.com/problemset/problem/7/C