欧几里德算法及其扩展算法

 一、欧几里德算法

1、欧几里德算法概述：

　　欧几里德算法又称辗转相除法，用于计算两个整数a,b的最大公约数。其计算原理依赖于下面的定理：

　　gcd函数就是用来求(a,b)的最大公约数的。

　　gcd函数的基本性质：

　　gcd(a,b)=gcd(b,a)=gcd(-a,b)=gcd(|a|,|b|)

2、欧几里得算法的公式表述

　　gcd(a,b)=gcd(b,a mod 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)的公约数是一样的，其最大公约数也必然相等，得证

3、欧几里德算法的C++语言描述

    int Gcd(int a, int b) 　　{ 　　  if(b == 0)　　  return a; 　　  return Gcd(b, a % b); 　　} 

当然你也可以写成迭代形式：

    int Gcd(int a, int b) 　　{ 　　  while(b != 0) 　　  { 　　    int r = b; 　　    b = a % b; 　　    a = r; 　　  } 　　  return a; 　　} 

二、扩展欧几里德算法

1、扩展欧几里德定理

　　对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然存在整

　　数对 x，y ，使得 gcd（a，b）=ax+by。

2、c++语言实现

   #include <iostream> 　　using namespace std; 　　int x,y,q; 　　void extend_Eulid(int a,int b) 　　{ 　　  if(b == 0)      {         x = 1;y = 0;q = a; 　　  } 　　  else 　　  { 　　    extend_Eulid(b,a%b);        int temp = x; 　　    x = y; 　　    y = temp - a/b*y; 　　  } 　　} 　　int main() 　　{ 　　  int a,b;     cin>>a>>b; 　　  if(a < b) swap(a,b); 　　  extend_Eulid(a,b); 　　  printf("%d=(%d)*%d+(%d)*%d\n",q,x,a,y,b); 　　  return 0; 　　} 

3、求解 x，y的方法的理解

　　设 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，所以递归可以

　　结束。

4、扩展欧几里德算法

　　扩展欧几里德算法是用来在已知a, b求解一组p，q使得ax+by = Gcd(a, b) =d(解一定存在，根据数论中的相关定理)。扩展欧几里德常用在求解模线性方程及方程组中。下面是一个使用C++的实现：

   int exGcd(int a, int b, int &x, int &y) 　　{ 　　  if(b == 0) 　　  { 　　  x = 1; 　　  y = 0; 　　  return a; ---很难找出一个这么实现的价值，因为扩展欧几里得还有更大的用途；个人认为定义全局数组更好,不用return r。 　　  } 　　  int r = exGcd(b, a % b, x, y); 　　  int t = x; 　　  x = y; 　　  y = t - a / b * y; 　　  return r; 　　} 

　　把这个实现和Gcd的递归实现相比，发现多了下面的x,y赋值过程，这就是扩展欧几里德算法的精髓。

　　可以这样思考:

　　对于a' = b, b' = a % b 而言，我们求得 x, y使得 a'x + b'y = Gcd(a', b')

　　由于b' = a % b = a - a / b * b (注：这里的/是程序设计语言中的除法)

　　那么可以得到:

　　a'x + b'y = Gcd(a', b') ===>

　　bx + (a - a / b * b)y = Gcd(a', b') = Gcd(a, b) ===>

　　ay +b(x - a / b*y) = Gcd(a, b)

　　因此对于a和b而言，他们的相对应的p，q分别是 y和(x-a/b*y)

（1）使用扩展欧几里德算法解决不定方程的办法

　　对于不定整数方程pa+qb=c，若 c mod Gcd(a, b)=0,则该方程存在整数解，否则不存在整数解。

　　上面已经列出找一个整数解的方法，在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后， /* p * a+q * b = Gcd(a, b)的其他整数解满足：

　　p = p0 + b/Gcd(a, b) * t

　　q = q0 - a/Gcd(a, b) * t(其中t为任意整数)

　　至于pa+qb=c的整数解，只需将p * a+q * b = Gcd(a, b)的每个解乘上 c/Gcd(a, b) 即可。*/有问题

　　在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后，应该是

　　得到p * a+q * b = c的一组解p1 = p0*(c/Gcd(a,b)),q1 = q0*(c/Gcd(a,b))，p * a+q * b = c的其他整数解满足：

　　p = p1 + b/Gcd(a, b) * t

　　q = q1 - a/Gcd(a, b) * t(其中t为任意整数)

　　p 、q就是p * a+q * b = c的所有整数解。

　　“有问题”的一段已经标出

　　也就是在得出p * a+q * b = Gcd(a, b)的解后是先乘上c/Gcd(a, b)从而导出p * a+q * b = c的其他整数解，

　　还是先导出p * a+q * b = Gcd(a, b)的其他整数解再乘上c/Gcd(a, b)。

　　编程时 exgcd 更多用于求解“中国余数定理”相关知识举个例子 比如n除以5余2 除以13余3 那么n最小是多少，所有的n满足什么条件？

　　n（min)=42

　　n=42+k*65

（2）欧几里德算法的扩展

　　扩展欧几里德算法不但能计算(a，b)的最大公约数，而且能计算a模b及b模a的乘法逆元，用C语言描述如下：

   int gcd(int a,int b,int &ar,int &br) 　　{ 　　  int x1,x2,x3; 　　  int y1,y2,y3; 　　  int t1,t2,t3; 　　  if(0 == a)      {        //有一个数为0,就不存在乘法逆元  　　   ar = 0;     　　br = 0 ;    　　 return b;   　　} 　　  if(0 == b) 　　  { 　　    ar = 0; 　　    br = 0 ; 　　    return a; 　　  }   　　x1 = 1;   　　x2 = 0;   　　x3 = a;   　　y1 = 0;   　　y2 = 1;   　　y3 = b;   　　int k;   　　for( t3 = x3 % y3 ; t3 != 0 ; t3 = x3 % y3)   　　{     　　k = x3 / y3;     　　t2 = x2 - k * y2;     　　t1 = x1 - k * y1;     　　x1 = y1;     　　x1 = y2;     　　x3 = y3;     　　y1 = t1;     　　y2 = t2;     　　y3 = t3;       }  　　if( y3 == 1)  　　{    　　//有乘法逆元    　　ar = y2;    　　br = x1;    　　return 1;  　　}  　　else 　　 {   　　//公约数不为1,无乘法逆元   　　ar = 0;   　　br = 0;   　　return y3; 　　  } 　 } 

　　扩展欧几里德算法对于最大公约数的计算和普通欧几里德算法是一致的。计算乘法逆元则显得很难明白。我想了半个小时才想出证明他的方法。

　　首先重复拙作整除中的一个论断：

　　如果gcd(a，b)=d，则存在m，n，使得d = ma + nb，称呼这种关系为a、b组合整数d，m，n称为组合系数。当d=1时，有 ma + nb = 1 ，此时可以看出m是a模b的乘法逆元，n是b模a的乘法逆元。

　　为了证明上面的结论，我们把上述计算中xi、yi看成ti的迭代初始值，考察一组数(t1，t2，t3)，用归纳法证明：当通过扩展欧几里德算法计算后，每一行都满足a×t1 + b×t2 = t3

　　第一行：1 × a + 0 × b = a成立

　　第二行：0 × a + 1 × b = b成立

　　假设前k行都成立，考察第k+1行

　　对于k-1行和k行有

　　t1(k-1) t2(k-1) t3(k-1)

　　t1(k) t2(k) t3(k)

　　分别满足：

　　t1(k-1) × a + t2(k-1) × b = t3(k-1)

　　t1(k) × a + t2(k) × b = t3(k)

　　根据扩展欧几里德算法，假设t3(k-1) = j t3(k) + r

　　则：

　　t3(k+1) = r

　　t2(k+1) = t2(k-1) - j × t2(k)

　　t1(k+1) = t1(k-1) - j × t1(k)

　　则

　　t1(k+1) × a + t2(k+1) × b

　　=t1(k-1) × a - j × t1(k) × a +

　　t2(k-1) × b - j × t2(k) × b

　　= t3(k-1) - j t3(k) = r

　　= t3(k+1)

　　得证

　　因此，当最终t3迭代计算到1时，有t1× a + t2 × b = 1，显然，t1是a模b的乘法逆元，t2是b模a的乘法逆元。

（3）C语言实现

   //扩展的欧几里德算法求乘法逆元 　　#include <stdio.h> 　　int ExtendedEuclid( int f,int d ,int *result); 　　int main() 　　{ 　　  int x,y,z; 　　  z = 0; 　　  printf("输入两个数：\n"); 　　  scanf("%d%d",&x,&y); 　　  if(ExtendedEuclid(x,y,&z)) 　　  printf("%d和%d互素，乘法的逆元是：%d\n",x,y,z); 　　  else 　　  printf("%d和%d不互素，最大公约数为：%d\n",x,y,z); 　　  return 0; 　　} 　　int ExtendedEuclid( int f,int d ,int *result) 　　{ 　    int x1,x2,x3,y1,y2,y3,t1,t2,t3,q; 　    x1 = y2 = 1;      x2 = y1 = 0; 　　  x3 = ( f>=d )?f:d; 　　  y3 = ( f>=d )?d:f; 　　  while( 1 ) 　　  { 　　    if ( y3 == 0 ) 　　    { 　　      *result = x3; /* 两个数不互素则result为两个数的最大公约数，此时返回值为零 */ 　　      return 0; 　　    } 　　    if ( y3 == 1 ) 　　    { 　　      *result = y2; /* 两个数互素则resutl为其乘法逆元，此时返回值为1 */ 　　      return 1; 　　    } 　　    q = x3/y3; 　　    t1 = x1 - q*y1; 　　    t2 = x2 - q*y2; 　　    t3 = x3 - q*y3; 　　    x1 = y1; 　　    x2 = y2; 　　    x3 = y3; 　　    y1 = t1; 　　    y2 = t2; 　　    y3 = t3; 　　  } 　　}