首页 > > 程序设计 > C/C++ >

扩展欧几里得算法

2018-06-17 23:43:25来源：未知阅读 ()

【转载】http://blog.csdn.net/qq_34494458/article/details/52637193

一：欧几里得算法（辗转相除法）

基本算法：设a=qb+r，其中a，b，q，r都是整数，则gcd(a,b)=gcd(b,r)，即gcd(a,b)=gcd(b,a%b)。

证明：

a可以表示成a = kb + r，则r = a mod b

　　假设d是a,b的一个公约数，则有

　　d|a, d|b，而r = a - kb，因此d|r

　　因此d是(b,a mod b)的公约数

　　假设d 是(b,a mod b)的公约数，则

　　d | b , d |r ，但是a = kb +r

　　因此d也是(a,b)的公约数

　　因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证

算法实现：

int gcd(int a,int b) {             //递归算法  
    return b ? gcd(b, a%b) : a;  
}  
  
int gcd(int a, int b) {      //迭代算法  
    while(b) {  
        int c = a%b;  
        a = b;  
        b = c;  
    }  
    return a;  
}

二扩展欧几里得算法：

基本算法：对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然存在整数对 x，y ，使得 gcd（a，b）=ax+by。

证明：设 a>b。

　　1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；

　　2，ab!=0 时

　　设 ax1+by1=gcd(a,b);

　　bx2+(a mod b)y2=gcd(b,a mod b);

　　根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);

　　则:ax1+by1=bx2+(a mod b)y2;

　　即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;

　　根据恒等定理得：x1=y2; y1=x2-(a/b)*y2;

这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.

　上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

//递归代码  
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, y, x);  
   int t = x;  
   y = y - a/b*t;  
   return r;  
  }
   // 非递归算法  
int exgcd(int m, int n, int &x, int &y) {  
   int x1, y1, x0, y0;  
   x0 = 1; y0 = 0;  
   x1 = 0; y1 = 1;  
   x = 0; y = 1;  
   int r = m%n;  
   int q = (m-r)/n;  
   while(r) {  
       x = x0 - q*x1;  
       y = y0 - q*y1;  
       x0 = x1; y0 = y1;  
       x1 = x; y1 = y;  
       m = n; n = r; r = m%n;  
       q = (m-r)/n;  
   }  
   return n;  
}

刘汝佳的：

void gcd(int a, int b, int& d, int& x, int& y) {  
    if (!b) {d = a; x = 1; y = 0; }  
    else {gcd(b, a%b, d, y, x); y -= x*(a/b); }  
}

上面求出的 x（当a>b时）都是最接近0的正整数。（不知道为什么）

扩展欧几里德算法的应用主要有以下两个方面：

（1）求解不定方程；

（2）求解模线性方程（线性同余方程）与逆元；

（1）求解不定方程；

1.对于不定整数方程ax+by=m，若 m mod Gcd(a, b)=0,则该方程存在整数解，否则不存在整数解。

证明：

首先扩展欧几里德主要是用来与求解线性方程相关的问题，所以我们从一个线性方程开始分析。现在假设这个线性方程为a*x+b*y=m，如果这个线性方程有解，那么一定有gcd(a,b) | m，即a，b的最大公约数能够整除m(m%gcd(a,b)==0)。证明很简单，由于a%gcd(a,b)==b%gcd(a,b)==0，所以a*x+b*y肯定能够整除gcd(a,b)，如果线性方程成立，那么就可以用m代替a*x+b*y，从而得到上面的结论，利用上面的结论就可以用来判断一个线性方程是否有解。
2.ax+by=Gcd（a， b）两边同时乘以 m/Gcd（a， b）得a（x*m/Gcd（a， b）)+b(y*m/Gcd（a， b）)=m;

上面已经列出找一个整数解的方法，在找到 a*x+ b*y = Gcd（a， b）的一组解x0,y0后，不定方程ax+by=m的一组解满足 x = m(x0)/gcd ， y = m*(y0)/gcd；
所以通解为 x = m*(x0)/gcd + k*lcm/a , y = m*(y0)/gcd + k*lcm/b (其中k为任意整数)；
证明：

令a1=a/gcd(a,b)，b1=b/gcd(a,b)，m1=m/gcd(a,b)。那么x，y的一组解就是x0*m1，y0*m1，但是由于满足方程的解无穷多个，在实际的解题中一般都会去求解x或是y的最小正数的值。以求x为例，又该如何求解呢？还是从方程入手，现在的x，y已经满足a*x+b*y=m，那么a*(x+n*b)+b*(y-n*a)=m显然也是成立的。可以得出x+n*b(n=…，-2，-1，0，1，2，…)就是方程的所有x解的集合，由于每一个x都肯定有一个y和其对应，所以在求解x的时候可以不考虑y的取值。取k使得x+k*b>0，x的最小正数值就应该是(x+k*b)%b，但是这个值真的是最小的吗？？如果我们将方程最有两边同时除以gcd(a,b)，则方程变为a1*x+b1*y=m1，同上面的分析可知，此时的最小值应该为(x+k*b1)%b1，由于b1<=b，所以这个值一定会小于等于之前的值。在实际的求解过程中一般都是用while(x<0)x+=b1来使得为正的条件满足，为了更快的退出循环，可以将b1改为b(b是b1的倍数)，并将b乘以一个倍数后再加到x上。

代码：

    //不定方程  
void linear_equation(int a, int b, int c, int &x, int &y) {  
    int d = exgcd(a, b, x, y);  
    if (c%d)  
        return ;  
    int k = c/d;  
    x *= k; y *= k; // 一组解  
    return ;  
}

相关题目：poj2115 poj2142 poj1061。

（2）求解模线性方程（线性同余方程）与逆元；

同余方程 ax≡b (mod n)对于未知数 x 有解，当且仅当 gcd(a,n) | b。且方程有解时，方程有 gcd(a,n) 个解。

求解方程 ax≡b (mod n) 相当于求解方程 ax+ （-ny）= b, (x, y为整数)。

　　//ax = b (mod n) 即 (ax) mod n = b mod n

　算法导论上有两个定理：

　　定理一：设d = gcd(a, n), 假定对整数x', y', 有d = ax' + ny'，如果d | b, 则方程ax = b(mod n)有一个解的值为x0, 满足：、

　　　　　　x0 = x'(b / d)(mod n）

　　定理二：假设方程ax = b(mod n)有解， x0是方程的任意一个解，则方程对模n恰有d个不同的解，分别为： xi = x0 + i * (n / d), 其中 i = 1,2,3......d - 1

　　有了这两个定理，解方程就不难了。

代码：

   // 模线性方程  
bool modular_linear_equation(int a, int b, int n) {  
    int x, y, x0, i;  
    int d = exgcd(a, n, x, y);  
    if (b%d)  
        return false;  
    x0 = x*(b/d)%n; //特解  
    for(i = 0; i < d; i++)  
        printf("%d\n", (x0 + i*(n/d))%n);  
    return true;  
}