之前一篇随笔"算法学习(1)----扩展欧几里得算法"记录了对朴素欧几里得算法和扩展欧几里得算法的学习和认识。学习所用书籍为 [美]Anany Levitin 所著《算法设计与分析基础》。
书后习题1.1 第10题 b 题如下：
改写扩展欧几里得算法对丢番图方程 ax+by=c 求解，系数 a,b,c 为任意整数。
看完题目开始思考，该题与扩展欧几里得算法有何联系？如果 c 为 a 和 b 的最大公约数，那么直接用扩展欧几里得算法就可以求出 x 和 y 了，但这里 a,b,c 为任意整数啊。
苦思片刻，终于顿悟：假设 a 和 b 的最大公约数为 d ，如果 c= k*d，那么先令 ax+by=d，再取 x=x*k, y=y*k 不就行了吗？如果 c 不是 d 的整数倍，那么是否就说明该方程无解呢？这个就不愿意再去费脑筋证明了，直接百度“丢番图方程”，果然验证了我的猜想：
不定方程一次不定方程是形式如a1x1 + a2x2 + ... + anxn= c的方程，一次不定方程有整数解的充要条件为：(a1,...,an)须是c的因子，其中(a1,...,an)表示a1,...,an的最大公因子。若有二元一次不定方程ax + by= c，且(a,b) | c，则其必有一组整数解x1,y1，并且还有以下关系式：
* x= x1 +[b / (a,b)]t
* y= y1 −[a / (a,b)]t
t为任意整数，故此一次不定方程有无限多解。请参见贝祖等式。
故得到代码如下：
using namespace std;
int diophantine(int a,int b,int c,int &x,int &y);
int exgcd(int m,int n,int &x,int &y);
int main(int argc,char* argv[])
{
int a,b,c,x,y,r;
while(true)
{
cout<<"Give me 3 interger(a,b,c): ";
cin>>a>>b>>c;
r=diophantine(a,b,c,x,y);
if(r==-1)
{
cout<<"There is no answers of x,y for "<<a<<" and "<<b<<"\n";
}
else
{
cout<<"a*x+b*y=c:"<<a<<"*"<<x<<"+"<<b<<"*"<<y<<"="<<c<<endl;
}
}
return 0;
}
int exgcd(int m,int n,int &x,int &y)
{
if(n==0){
x=1,y=0;
return m;
}
int r=exgcd(n,m%n,y,x);
y -= m/n*x;
return r;
}
int diophantine(int a,int b,int c,int &x,int &y)
{
int r=exgcd(a,b,x,y);
int k=c/r;
if(k*r==c){
x*=k,y*=k;
return r;
}
return -1;
}
在Linux控制台编译运行验证：
成功！
