如何求解二元一次不定方程的整数解

这里讨论的二元一次不定方程专指ax+by=c(a*b≠0,a,b,c∈Z)-----①
定理一:
方程①有整数解的充分必要条件是(a,b)|c((a,b)即Gcd(a,b),下同)
定理二(裴蜀定理)
设(a,b)=1,则方程①的全部整数解可表为:
x=x0+bt,y=y0-at(x0,y0是方程的一组整数解,t为任意整数)
也可以表述为:
方程①的全部整数解可表为:
x=x0+bt/(a,b),y=y0-at/(a,b)(x0,y0是方程的一组整数解,t为任意整数)

我以前的一个想法是求出(a,b)并判断是否(a,b)|c,然后约简a,b,c(都除以(a,b)),这样(a,b)=1,如果方程①有整数解的话,根据定理二,在一个长度为|b|的整数范围内必存在一个x的解,在一个长度为|a|的整数范围内必存在一个y的解。那么在x的[0,|b|-1]或者y的[0,|a|-1]上试解吧,找到x0可以求出y0,找到y0可以求出x0。为了减少搜索次数,我们先比较|a|和|b|,如果|a|≤|b|就在y上试解,反之则在x上试解。我用这种办法ac了sgu106和poj1061,但是这种办法的代价并非是最低的,因为有一个著名的欧几里德扩展算法。
这个算法源自于求gcd的欧几里德算法,二者的算法复杂度基本上是一样的。
回忆一下欧几里德算法:
有两个数a,b,我们把它们写成u0和u1,求gcd的步骤如下:
u0=q0*u1+u2
u1=q1*u2+u3
u2=q2*u3+u4
...
uj=qj*u(j+1)
u(j+1)便是a和b的最大公约数。
要求出一组整数解,可以先求出ax+by=(a,b)的一组解x0,y0,然后x0=x0*c/(a,b),y0=x0*c/(a,b).
同样地,把a和b写成u0和u1,设ui=u0*xi+u1*yi(必然存在这样的xi和yi,因为uj是gcd(u0,u1)的倍数),
那么x0=1,y0=0,x1=0,y1=1;
u(i+1)=u(i-1)-q(i-1)*ui=u0*x(i-1)+u1*y(i-1)-q(i-1)*u0*x(i)-q(i-1)*u1*y(i)=u0*(x(i-1)-q(i-1)*x(i))+u1*(y(i-1)-q(i-1)*y(i)),而u(i+1)=u0*x(i+1)+u1*y(i+1),因此
x(i+1)=x(i-1)-q(i-1)*x(i)
y(i+1)=y(i-1)-q(i-1)*y(i)
伟大的递推式!
我们要求的是xj和yj(uj=u0*xj+u1*yj,uj为gcd),所以只要递推到xj和yj即可。
以下程序演示了欧几里德扩展算法:
#include<stdio.h>
int Ext_Gcd(int a,int b,int &x,int &y) //返回的Gcd带符号
{
int q,tmp,x0=0,x1=1,y0=1,y1=0;
while(a)
{
q=b/a;
tmp=a;
a=b%a;
b=tmp;
tmp=x1;
x1=x0-q*x1;
x0=tmp;
tmp=y1;
y1=y0-q*y1;
y0=tmp;
}
x=x0;
y=y0;
return b;
}
int main()
{
int a,b,c,x,y,gcd;
printf("*****欧几里德扩展算法演示程序*****/n/n");
for(;;)
{
printf("请输入方程的参量:/n");
scanf("%d %d %d",&a,&b,&c);
gcd=Ext_Gcd(a,b,x,y);
if (c%gcd) {printf("该方程无解。/n");continue;}
c/=gcd;
printf("a,b的公约数为:%d。/n方程的标准解为:%d %d/n/n",gcd,x*c,y*c);
}
}