算法学习(2)----丢番图方程

　　之前一篇随笔"算法学习(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<<"
";
        }
        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):

成功！