欧几里德算法

欧几里得算法的概述

欧几里德算法又称辗转相除法，用于计算两个正整数a，b的最大公约数。其计算原理依赖于下面的定理： 　　定理：gcd(a,b) = gcd(b,a mod b) (a>b 且a mod b 不为0) 　　证明：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)的公约数 　　因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证

编辑本段欧几里得算法原理

Lemma 1.3.1 若 a, b 且 a = bh + r, 其中 h, r , 则 gcd(a, b) = gcd(b, r). 　　

证 明. 假设 d1 = gcd(a, b) 且 d2 = gcd(b, r).

我们证明 d1| d2 且 d2| d1, 因而可利用 Proposition 1.1.3(2) 以及 d1, d2 皆为正数得证 d1 = d2. 　　

因 d1| a 且 d1| b 利用 Corollary 1.1.2 我们知 d1| a - bh = r.

因为 d1| b, d1| r 且 d2 = gcd(b, r) 故由 Proposition 1.2.5 知 d1| d2. 另一方面, 因为 d2| b 且 d2| r 故 d2| bh + r = a. 因此可得 d2| d1. 　　

Lemma 1.3.1 告诉我们当 a > b > 0 时, 要求 a, b 的最大公因数我们可以先将 a 除以 b 所得馀数若为 r, 则 a, b 的最大公因数等于 b 和 r 的最大公因数. 因为 0r < b < a, 所以当然把计算简化了. 接着我们就来看看辗转相除法. 由于 gcd(a, b) = gcd(- a, b) 所以我们只要考虑 a, b 都是正整数的情况. 　　

Theorem 1.3.2 (The Euclidean Algorithm) 假设 a, b 且 a > b. 由除法原理我们知存在 h0, r0 使得　　

a = bh0 + r0, 其中 0r0 < b. 　　若 r0 > 0, 则存在 h1, r1 使得　　b = r0h1 + r1, 其中 0r1 < r0. 　　

若 r1 > 0, 则存在 h2, r2 使得　　r0 = r1h2 + r2, 其中 0r2 < r1. 　　

如此继续下去直到 rn = 0 为止. 若 n = 0 (即 r0 = 0), 则 gcd(a, b) = b.

若 n1, 则 gcd(a, b) = rn - 1. 　　

证 明. 首先注意若 r0 0, 由于 r0 > r1 > r2 > ... 是严格递减的, 因为 r0 和 0 之间最多仅能插入 r0 - 1 个正整数, 所以我们知道一定会有 nr0 使得 rn = 0. 　　若 r0 = 0, 即 a = bh0, 故知 b 为 a 之因数, 得证 b 为 a, b 的最大公因数.

若 r0 > 0, 则由 Lemma 1.3.1 知　　gcd(a, b) = gcd(b, r0) = gcd(r0, r1) = ... = gcd(rn - 1, rn) = gcd(rn - 1, 0) = rn - 1. 　　

现在我们来看用辗转相除法求最大公因数的例子　　

Example 1.3.3 我们求 a = 481 和 b = 221 的最大公因数.

首先由除法原理得 481 = 2 . 221 + 39, 知 r0 = 39. 因此再考虑 b = 221 除以 r0 = 39 得 221 = 5 . 39 + 26, 知 r1 = 26. 再以 r0 = 39 除以 r1 = 26 得 39 = 1 . 26 + 13, 知 r2 = 13. 最后因为 r2 = 13 整除 r1 = 26 知 r3 = 0, 故由 Theorem 1.3.2 知 gcd(481, 221) = r2 = 13. 　　

在利用辗转相除法求最大公因数时, 大家不必真的求到 rn = 0. 例如在上例中可看出 r0 = 39 和 r1 = 26 的最大公因数是 13, 利用 Lemma 1.3.1 马上得知 gcd(a, b) = 13. 　　

在上一节 Corollary 1.2.5 告诉我们若 gcd(a, b) = d, 则存在 m, n 使得 d = ma + nb. 当时我们没有提到如何找到此 m, n. 现在我们利用辗转相除法来介绍一个找到 m, n 的方法.

我们沿用 Theorem 1.3.2 的符号. 首先看 r0 = 0 的情形, 此时 d = gcd(a, b) = b 所以若令 m = 0, n = 1, 则我们有 d = b = ma + nb. 当 r0 0 但 r1 = 0 时, 我们知 d = gcd(a, b) = r0. 故利用 a = bh0 + r0 知, 若令 m = 1, n = - h0, 则 d = r0 = ma + nb. 同理若 r0 0, r1 0 但 r2 = 0, 则知 d = gcd(a, b) = r1. 故利用 a = bh0 + r0 以及 b = r0h1 + r1 知　　r1 = b - r0h1 = b - (a - bh0)h1 = - h1a + (1 + h0h1)b. 　　

因此若令 m = - h1 且 n = 1 + h0h1, 则 d = r1 = ma + nb. 依照此法, 当 r0, r1 和 r2 皆不为 0 时, 由于 d = gcd(a, b) = rn - 1 故由 rn - 3 = rn - 2hn - 1 + rn - 1 知 d = rn - 3 - hn - 1rn - 2. 利用前面推导方式我们知存在 m1, m2, n1, n2 使得 rn - 3 = m1a + n1b 且 rn - 2 = m2a + n2b 故代入得　　

d = (m1a + n1b) - hn - 1(m2a + n2b) = (m1 - hn - 1m2)a + (n1 - hn - 1n2)b. 　　

因此若令 m = m1 - hn - 1m2 且 n = n1 - hn - 1n2, 则 d = ma + nb. 　　

上面的说明看似好像当 r0 0 时对每一个 i {0, 1,..., n - 2} 要先将 ri 写成 ri = mia + nib, 最后才可将 d = rn - 1 写成 ma + nb 的形式. 其实这只是论证时的方便, 在实际操作时我们其实是将每个 ri 写成 mi'ri - 2 + ni'ri - 1 的形式慢慢逆推回 d = ma + nb. 请看以下的例子. 　　

Example 1.3.4 我们试着利用 Example 1.3.3 所得结果找到 m, n 使得 13 = gcd(481, 221) = 481m + 221n. 首先我们有 13 = r2 = 39 - 26 = r0 - r1. 而 r1 = 221 - 5 . 39 = b - 5r0, 故得 13 = r0 - (b - 5r0) = 6r0 - b. 再由 r0 = 481 - 2 . 221 = a - 2b, 得知 13 = 6(a - 2b) - b = 6a - 13b. 故得 m = 6 且 n = - 13 会满足 13 = 481m + 221n. 　　

要注意这里找到的 m, n 并不会是唯一满足 d = ma + nb 的一组解. 虽然上面的推演过程好像会只有一组解, 不过只能说是用上面的方法会得到一组解, 并不能担保可找到所有的解. 比方说若令 m' = m + b, n' = n - a, 则 m'a + n'b = (m + b)a + (n - a)b = ma + nb = d. 所以 m', n' 也会是另一组解. 所以以后当要探讨唯一性时, 若没有充分的理由千万不能说由前面的推导过程看出是唯一的就断言是唯一. 一般的作法是假设你有两组解, 再利用这两组解所共同满足的式子找到两者之间的关系. 我们看看以下的作法. 　　

Proposition 1.3.5 假设 a, b 且 d = gcd(a, b). 若 x = m0, y = n0 是 d = ax + by 的一组整数解, 则对任意 t , x = m0 + bt/d, y = n0 - at/d 皆为 d = ax + by 的一组整数解, 而且 d = ax + by 的所有整数解必为 x = m0 + bt/d, y = n0 - at/d 其中 t 这样的形式. 　　

证 明. 假设 x = m, y = n 是 d = ax + by 的一组解. 由于已假设 x = m0, y = n0 也是一组解, 故得 am + bn = am0 + bn0. 也就是说 a(m - m0) = b(n0 - n). 由于 d = gcd(a, b), 我们可以假设 a = a'd, b = b'd 其中 a', b' 且 gcd(a', b') = 1 (参见 Corollary 1.2.3). 因此得 a'(m - m0) = b'(n0 - n). 利用 b'| a'(m - m0), gcd(a', b') = 1 以及 Proposition 1.2.7(1) 得 b'| m - m0. 也就是说存在 t 使得 m - m0 = b't. 故知 m = m0 + b't = m0 + bt/d. 将 m = m0 + bt/d 代回 am + bn = am0 + bn0 可得 n = n0 - at/d, 因此得证 d = ax + by 的整数解都是 x = m0 + bt/d, y = n0 - at/d 其中 t 这样的形式. 最后我们仅要确认对任意 t , x = m0 + bt/d, y = n0 - at/d 皆为 d = ax + by 的一组整数解. 然而将 x = m0 + bt/d, y = n0 - at/d 代入 ax + by 得 a(m0 + bt/d )+ b(n0 - at/d )= am0 + bn0 = d, 故得证本定理. 　　利用 Proposition 1.3.5 我们就可利用 Example 1.3.4 找到 13 = 481x + 221y 的一组整数解 x = 6, y = - 13 得到 x = 6 + 17t, y = - 13 - 37t 其中 t 是 13 = 481x + 221y 所有的整数解.

编辑本段欧几里得算法设计

辗转相除法是利用以下性质来确定两个正整数 a 和 b 的最大公因子的：　　

1. 若 r 是 a ÷ b 的余数,且r不为0， 则　　gcd(a,b) = gcd(b,r) 　　

2. a 和其倍数之最大公因子为 a。　　

另一种写法是：　　

1. 令r为a/b所得余数（0≤r＜b）　　若 b= 0，算法结束；a 即为答案。　　

2. 互换：置 a←b，b←r，并返回第一步。

编辑本段欧几里德算法的C语言版

/*欧几里德算法：辗转求余　　

原理: gcd(a,b)=gcd(b,a mod b) 　　

当b为0时,两数的最大公约数即为a 　　

getchar()会接受前一个scanf的回车符　　

*/ 　　

#include<stdio.h> 　　

unsigned int Gcd( unsigned int M, unsigned int N ) 　　

{ 　　

　　unsigned int Rem; 　　

　　while( N > 0 ) 　　

　　{ 　　

　　　　Rem = M % N; 　　

　　　　M = N; 　　

　　　　N = Rem; 　　

　　} 　　

　　return M; 　　

} 　　

void main() 　　

{ 　　

　　int temp; 　　

　　int a,b; 　　

　　scanf("%d",&a); 　　

　　scanf("%d",&b); 　　

　　printf("the greatest common factor of %d and %d is ",a,b); 　　

　　printf("%d\n",Gcd(a,b)); 　　

}

编辑本段欧几里德算法的C++/java语言版

欧几里德算法就是根据这个原理来做的，其算法用C++语言描述为： 　　void swap(int &； a， int &； b) 　　{ 　　int c = a；　　a = b；　　b = c；　　} 　　int gcd(int a，int b) 　　{ 　　if(0 == a ) 　　{ 　　return b；　　} 　　if( 0 == b) 　　{ 　　return a；　　} 　　if(a > b) 　　{ 　　swap(a，b)；　　} 　　int c；　　for(c = a % b ； c >0 ； c = a % b) 　　{ 　　a = b；　　b = c；　　} 　　return b；　　} 　　用java描述为: 　　int divisor(int m, int n) { 　　if (m % n == 0) { 　　return n; 　　} else { 　　return divisor(n, m % n); 　　} 　　} 　　模P乘法逆元 　　对于整数a、p，如果存在整数b，满足ab mod p =1，则说，b是a的模p乘法逆元。　　定理：a存在模p的乘法逆元的充要条件是gcd(a，p) = 1 　　证明： 　　首先证明充分性 　　如果gcd(a，p) = 1，根据欧拉定理，aφ(p) ≡ 1 mod p，因此　　显然aφ(p)-1 mod p是a的模p乘法逆元。　　再证明必要性 　　假设存在a模p的乘法逆元为b 　　ab ≡ 1 mod p 　　则ab = kp +1 ，所以1 = ab - kp 　　因为gcd(a，p) = d 　　所以d | 1 　　所以d只能为1

编辑本段求最大公约数的Stein算法

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

Stein算法由J. Stein于1961年提出，这个方法也是计算两个数的最大公约数。和欧几里德算法不同的是，Stein算法只有整数的移位和加减法，这对于程序设计者是一个福音。为了说明Stein算法的正确性，首先必须注意到以下结论：　　

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

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

C++/java 实现　　// c++/java stein 算法　　

int gcd(int a,int b)

{ 　　

　　if(a<b){//arrange so that a>b 　　

　　int temp = a; 　　

　　a = b; 　　

　　b=temp; 　　

　　} 　　

　　if(0==b)//the base case 　　

　　return a; 　　

　　if(a%2==0 && b%2 ==0)//a and b are even 　　

　　return 2*gcd(a/2,b/2); 　　

　　if ( a%2 == 0)// only a is even 　　

　　return gcd(a/2,b); 　　

　　if ( b%2==0 )// only b is even 　　

　　return gcd(a,b/2); 　　

　　return gcd((a+b)/2,(a-b)/2);// a and b are odd 　　

}

编辑本段欧几里德算法的扩展

扩展欧几里德算法不但能计算(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的乘法逆元。