欧几里得算法又称辗转相除法，是指用于计算两个非负整数a，b的最大公约数。应用领域有数学和计算机两个方面。计算公式gcd(a,b) = gcd(b,a mod b)1。
欧几里得算法和扩展欧几里得算法可使用多种编程语言实现。欧几里得算法是用来求两个正整数最大公约数的算法。古希腊数学家欧几里得在其著作《The Elements》中最早描述了这种算法,所以被命名为欧几里得算法。
欧几里得算法
Euclidean Algorithm
Euclid's algorithm
辗转相除法
数学、计算机
计算两个正整数a、b的最大公约数
gcd(a,b) = gcd(b,a mod b)
数学
算法简介
欧几里得算法是用来求两个正整数最大公约数的算法。古希腊数学家欧几里得在其著作《The Elements》中最早描述了这种算法,所以被命名为欧几里得算法。
扩展欧几里得算法可用于RSA加密等领域。
假如需要求 1997 和 615 两个正整数的最大公约数,用欧几里得算法，是这样进行的：
1997 / 615 = 3 (余 152)
615 / 152 = 4(余7)
152 / 7 = 21(余5)
7 / 5 = 1 (余2)
5 / 2 = 2 (余1)
2 / 1 = 2 (余0)
至此，最大公约数为1
以除数和余数反复做除法运算，当余数为 0 时，取当前算式除数为最大公约数，所以就得出了 1997 和 615 的最大公约数 1。
计算证明
其计算原理依赖于下面的定理：
定理：两个整数的最大公约数等于其中较小的那个数和两数相除余数的最大公约数。最大公约数（Greatest Common Divisor）缩写为GCD。
gcd(a,b) = gcd(b,a mod b) (不妨设a>b 且r=a mod b ,r不为0)
证法一
a可以表示成a = kb + r（a，b，k，r皆为正整数，且r
假设d是a,b的一个公约数，记作d|a,d|b，即a和b都可以被d整除。
而r = a - kb，两边同时除以d，r/d=a/d-kb/d，由等式右边可知m=r/d为整数，因此d|r
因此d也是b,a mod b的公约数。
因(a,b)和(b,a mod b)的公约数相等，则其最大公约数也相等，得证。
证法二
假设c = gcd(a,b),则存在m,n，使a = mc, b = nc;
令r = a mod b，即存在k，使r = a-kb = mc - knc = (m-kn)c;
故gcd(b,a mod b) = gcd(b,r) = gcd(nc,(m-kn)c) = gcd(n,m-kn)c;
则c为b与a mod b的公约数;
假设d = gcd(n,m-kn), 则存在x,y, 使n = xd, m-kn = yd; 故m = yd+kn = yd+kxd = (y+kx)d;
故有a = mc = (y+kx)dc, b = nc = xdc; 可得 gcd(a,b) = gcd((y+kx)dc,xdc) = dc;
由于gcd(a,b) = c, 故d = 1;
即gcd(n,m-kn) = 1, 故可得gcd(b,a mod b) = c;
故得证gcd(a,b) = gcd(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⑵ 以及 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 的最大公因数. 因为 r < 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，其中 r0 < 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⑴ 得 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 的最大公因子的：
⒈ 若 r 是 a ÷ b 的余数，且r不为0， 则
gcd(a,b) = gcd(b,r)
⒉ a 和其倍数之最大公因子为 a。
另一种写法是：
⒈ 令r为a/b所得余数（0≤r
若 r= 0，算法结束；b 即为答案。
⒉ 互换：置 a←b，b←r，并返回第一步。
算法版本
Swift 语言版本
import Foundation
func gcd(a: Int, b: Int) -> Int {
var x = a, y = b
while y != 0 {
(x, y) = (y, x % y)
}
return x
}
let x = 75, y = 100
let result = gcd(a: x, b: y)
print("/(x) 和 /(y) 的最大公约数是 /(result)")
Go语言版本
package main
import "fmt"
func main() {
var x, y int = 18, 12
result := gcd(x,y)
fmt.Printf("x, y 的最大公约数是 : %d",result)
}
func gcd(x,y int) int{
for y != 0 {
x, y = y, x%y
}
return x
}
Pascal语言版
var a,b,c:integer;
begin
readln(a,b);
c:=a mod b;
while c<>0 do
begin
a:=b;b:=c;c:=a mod b;
end;
write(b);
end.
C语言版
/*
欧几里得算法：辗转求余
原理： gcd(a,b)=gcd(b,a mod b)
当b为0时，两数的最大公约数即为a
getchar()会接受前一个scanf的回车符
*/
#include<stdio.h>
unsigned int MaxCommonFactor(int a,int b)
{
if(a % b == 0)
return b;
return MaxCommonFactor(b, a % b);
}
unsigned int Gcd(unsigned int M,unsigned int N)
{
unsigned int Rem;
while(N > 0)
{
Rem = M % N;
M = N;
N = Rem;
}
return M;
}
int main(void)
{
int a,b;
scanf("%d %d",&a,&b);
printf("the greatest common factor of %d and %d is ",a,b);
printf("%d/n",Gcd(a,b));
printf("recursion:%d/n",MaxCommonFactor(a,b));
return 0;
}
Ruby语言版
#用欧几里得算法计算最大公约数(排版略）
def gcd(x, y)
if y == 0
return x
else
return gcd(y, x % y)
end
end
C++版
#include <iostream>
using namespace std;
int gcd(int a, int b)
{
if (a % b==0) return b;
else return gcd(b, a % b);
}
int x, y;
int main(){
cin >> x >> y;
cout << gcd(x, y);
return 0;
}
Java版
int gcd(int m,int n)
{ if(n == 0){
return m;
}
int r = m%n;
return gcd(n,r);
}
JavaScript版
function gcd(a, b) {
if (a % b == 0) return b;
return gcd(b, a % b);
}
Python版
def gcd(a, b):
while a != 0:
a, b = b % a, a
return b
Erlang版
gcd(A, 0) -> A;
gcd(A, B) -> gcd(B, A rem B).
Rust版
pub fn gcd(x: u64, y: u64) -> u64 {
let remainder = x % y;
if remainder == 0 {
return y;
} else {
return gcd(y, remainder);
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn gcd_works() {
assert_eq!(gcd(2, 4), 2);
assert_eq!(gcd(6, 27), 3);
assert_eq!(gcd(4, 2), 2);
assert_eq!(gcd(27, 6), 3);
}
}
Bash版
function gcd() {
if [ ! -n "$2" ];then
return
fi
if [ "$2" == "0" ];then
echo "$1"
return
fi
gcd "$2" "$[$1%$2]"
}
模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(ab
{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,b);// a and b are odd
}
算法扩展
扩展欧几里得算法不但能计算(a，b）的最大公约数，而且能计算a模b及b模a的乘法逆元，用C语言描述如下：
#include <stdio.h>
unsigned int gcdExtended( int a, int b, int *x, int *y);
int main(void) {
int a, b,GCD;
int x, y;
a = 1232, b = 573;
/*
gcdExtended(1232, 573)时, x = 20 and y = –43
1232x + 573y = 1
24640-24639 = 1
或者gcdExtended( 573,1232) 时,x=-43, y=20
573x+1232y = 1
-43*573+1232*20 = -24639+57640 = 1
gcdExtended(9151, 5787) 时
x=2011, y=-3180
*/
GCD = gcdExtended(a, b,&x, &y);
printf("gcdExtended(%d, %d) = %d, x=%d, y=%d/n", a, b, GCD,x,y);
return 0;
}
// 欧几里得扩展算法的C语言实现
// ax+by=1
unsigned int gcdExtended(int a, int b, int *x, int *y){
if (a == 0){
*x = 0;
*y = 1;
return b;
}
int x1, y1;
int gcd = gcdExtended(b%a, a, &x1, &y1);
*x = y1 - (b/a) * x1;
*y = x1;
return gcd;
}
扩展欧几里得算法对于最大公约数的计算和普通欧几里得算法是一致的。计算乘法逆元则显得很难明白。我想了半个小时才想出证明他的方法。
首先重复拙作整除中的一个论断：
如果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的乘法逆元。
声明：本文部分文字与图片资源来自于网络，转载此文是出于传递更多信息之目的,若有来源标注错误或侵犯了您的合法权益，请立即【留言反馈】通知我们，情况属实，我们会第一时间予以删除，并同时向您表示歉意本文作者为北京精神，转载请注明。