数论基础——扩展欧几里德算法解析

扩展欧几里德算法是用来在已知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