4、欧几里德及其扩展算法，Stein算法（高效有用）

1、欧几里德算法

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

定理：gcd(a,b) = gcd(b,a mod b)，其中，gcd(a,b)就表示求最大公约数。

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

假设d是a,b的一个公约数，则有a=d*k1,b=d*k2,r=a-kb=d*k1-k*d*k2=d*(k1-k*k2),也就是说d也是r的公约数。

假设d 是(b,a mod b)的公约数，则b=d*k1,r=d*k2,而因为r=a-kb,所以a=kb+r=k*d*k1+d*k2=d*(k*k1+k2),因此d也是(a,b)的公约数。因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证

其算法用描述为：

E1： r=m%n

E2   if r= 0

        return n  //(即为最大公因数)

       else

       ｛ m = n; n = r; ｝
       返回E1

     算法用程序实现为：

int gcb(int m,int n)

{

int r;

r=m%n;

if(r==0) return n;

while(r)

 {

  m=n;

  n=r;

  r=m%n; 

 }//while

return n;

}

这两个数的最小公倍倍数也就是他们的乘积除以最大公约数了。

2、扩展的欧几里德算法

A）给定的两个正整数a,b，计算它们的最大公因子d和两个整数a和b，使得ax+by=d，即gcd(a,b)=ax+by。

要证明这个定理，首先需给出一个定理带余除法：

若a,b是两个整数，其中b>0,则存在两个整数q,r，使得

               a = b*q + r  0<=r<b

其中q和r是唯一的。

这个定理证明是简单的，有了这个定义，我们就来证明欧几里德算法          

a = q0*b + r0

b = q1*r0 + r1

r0= q2*r1 + r2

 ...

rn-3=qn-1*rn-2+ r n-1

rn-2 = qn* rn-1 + r n(r n = 0)

这里r n-1就是最大公约数。

现在就需要反推回去，可以得到r可以用a,b线性表示

ri=(Xi-2-Qi*Xi-1)*a-(Yi-2-Qi*Yi-1)*b

下一步就是求x,y的值。

Xi=Xi-2-Qi*Xi-1

Yi=Yi-2-Qi*Yi-1

B）求解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%b)y2=gcd(b,a%b); 根据朴素欧几里德原理有
gcd(a,b)=gcd(b,a%b); 则:ax1+by1=bx2+(a%b)y2;

即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+b(x2-(a/b)*y2); 根据恒等定理得：x1=y2;y1=x2-(a/b)*y2; 这样我们就得到了求解x1,y1的方法：x1，y1的值基于x2，y2,依此类推，这就是递归了。

C)程序描述

#include<stdio.h> 

int x,y,q; 

void swap(int&m,int &n)

{

int temp;

temp=m;

m=n;

n=temp;

}

 

void extend_Eulid(inta,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;

 } //else

} //extend_Euild

 

int main() 

{ int a,b; 

 scanf("%d",&a); 

 scanf("%d",&b); 

if (a<b)  swap(a,b); 

extend_Eulid(a,b); 

printf("%d=(%d)*%d+(%d)*%d\n",q,x,a,y,b);

} //main

3、Stein算法

欧几里德算法是计算两个数最大公约数的传统算法，他无论从理论还是从效率上都是很好的。但是他有一个致命的缺陷，这个缺陷只有在大素数时才会显现出来。

考虑现在的硬件平台，一般整数最多也就是64位，对于这样的整数，计算两个数之间的模是很简单的。对于字长为32位的平台，计算两个不超过32位的整数的模，只需要一个指令周期，而计算64位以下的整数模，也不过几个周期而已。但是对于更大的素数，这样的计算过程就不得不由用户来设计，为了计算两个超过 64位的整数的模，用户也许不得不采用类似于多位数除法手算过程中的试商法，这个过程不但复杂，而且消耗了很多CPU时间。对于现代密码算法，要求计算 128位以上的素数的情况比比皆是，设计这样的程序迫切希望能够抛弃除法和取模。

Stein算法由J. Stein 1961年提出，这个方法也是计算两个数的最大公约数。和欧几里德算法算法不同的是，Stein算法只有整数的移位和加减法。

首先必须注意到以下结论：

（1）gcd(a,a) = a，也就是一个数和他自身的公约数是其自身

（2）gcd(ka,kb) = k gcd(a,b)，也就是最大公约数运算和倍乘运算可以交换，特殊的，当k=2时，说明两个偶数的最大公约数必然能被2整除。

Stein算法描述如下：

（1）     如果A=0，B是最大公约数，算法结束

（2）     如果B=0，A是最大公约数，算法结束

（3）     设置A1 = A、B1=B和C1 = 1

（4）     如果An和Bn都是偶数，则An+1 =An /2，Bn+1=Bn /2，Cn+1 =Cn *2(注意，乘2只要把整数左移一位即可，除2只要把整数右移一位即可)

（5）     如果An是偶数，Bn不是偶数，则An+1 =An /2，Bn+1 =Bn ，Cn+1 =Cn (很显然啦，2不是奇数的约数)

（6）     如果Bn是偶数，An不是偶数，则Bn+1 =Bn /2，An+1 =An ，Cn+1 =Cn (很显然啦，2不是奇数的约数)

（7）     如果An和Bn都不是偶数，则An+1 =|An -Bn|，Bn+1 =min(An,Bn)，Cn+1 =Cn

n++，转（4）.

     注：上述中n+1,n等都是下标。

程序描述如下：

（1）int gcd(inta,int b)

{

if (a == 0) return b;

if (b == 0) return a;

if (a % 2 == 0 && b% 2 == 0) return 2 * gcd(a/2,b/2);

else if (a % 2 == 0) returngcd(a/2,b);

else if (b % 2 == 0) returngcd(a,b/2);

else returngcd(abs(a-b),min(a,b));

}

(2)intgcd(int a,int b)

{

if(a == 0) return b;

if(b == 0) return a;

if(a % 2 == 0 && b % 2 == 0) return gcd(a/2,b/2) << 1;

elseif (a % 2 == 0) return gcd(a>>1,b);

elseif (b % 2 == 0) return gcd(a,b>>1);

elsereturn gcd(abs(a-b),min(a,b));

}