之前一篇随笔"算法学习(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为任意整数,故此一次不定方程有无限多解。请参见贝祖等式。
故得到代码如下:
// diophantine.cpp
#include <iostream>
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;
}
编译运行测试:
root@javis:~/Documents/Code$ g++ diophantine.c -o diophantine.out
root@javis:~/Documents/Code$ ./diophantine.out
Give me 3 interger(a,b,c): 3 4 5
a*x+b*y=c:3*-5+4*5=5
Give me 3 interger(a,b,c): 6 3 5
There is no answers of x,y for 6 and 3
Give me 3 interger(a,b,c): 9 5 4
a*x+b*y=c:9*-4+5*8=4
Give me 3 interger(a,b,c): 8 4 2
There is no answers of x,y for 8 and 4
Give me 3 interger(a,b,c): 100 44 88
a*x+b*y=c:100*88+44*-198=88
Give me 3 interger(a,b,c):
成功!