首页 > Python算法 > 数论算法 阅读数:50

欧几里得算法

很多数学算法,尤其是离散数学的算法,在计算机领域中拥有重要的位置。离散数学中包括数论与图论,我们已经在上节了解了一些基本的图论算法,接下来我们学习数论算法中的欧几里得算法。

利用欧几里得算法(Euclidean Algorithm,又称辗转相除法),我们可以得出任何两个自然数的最大公约数(Greatest Common Divisor,GCD)。最大公约数指最大的更够被两个数字都整除的数字。比如,5 与 15 的最大公约数为 5,6 与 8 的最大公约数为 2,3 与 2 的最大公约数为 1。

当遇到较大的数字时,如 409 和 1078、1030498 和 45091234,我们怎样得出它们的最大公约数呢?本节教程,我们学习欧几里得在公元前 300 年想出来的解题方法。

欧几里得算法是数论中一个非常重要的基本工具,在密码学中尤其重要。它的应用十分广泛,包括解贝祖等式,解线性丢番图方程,在下节的中国余数定理中我们还会用到它来寻找满足条件的解。分析欧几里得算法后我们还会了解一下它的两个应用。

1. 算法分析与证明

两个自然数 a 与 b 的最大公约数的正式的书写方式为 gcd(a,b)。举例,gcd(5,15)=5。在学习算法之前,我们需要首先理解欧几里得证明的以下定理:

定理 1:

gcd(a,b)=gcd(b,r)  其中 a,b,r 为自然数,a>b 且 a=kb+r,k 为正整数。定理表明两个数字的最大公约数同时也是较小数和两个数字的余数的最大公约数。


例如,gcd(36,15)=gcd(15,6),因为 36=2×15+6。

定理证明如下:
  1. 令 g 为 a 与 b 的公约数,g 不一定是最大公约数;
  2. 因为 a 与 b 为 g 的倍数,所以可以表示为 a=mg,b=ng,m 与 n 为自然数;
  3. 由 a=kb+r 得 r=a−kb。又由(3)得 r=mg−kng=(m−kn)g;
  4. 从(3)得知 g 为 r 的因子。但同时 g 又为 b 的因子,所以 g 肯定也是 b 与 r 的公约数;
  5. 从(4)得知所有 a 与 b 的公约数都是 b 与 r 的公约数,所以 a 与 b 的最大公约数肯定也是 b 与 r 的最大公约数,即 gcd(a,b)=gcd(b,r)。

欧几里得算法如下:
  1. 已知两个自然数 a 与 b,a>b;
  2. 计算 r,r=a%b,即 a 除 b 的余数;
  3. 如果 r 为 0,代表 a 可以被 b 整除,输出 b。如果 r 不为 0,返回步骤(1),设 a 为 b,b 为 r。

我们通过一个例子理解以上过程:找 54 与 20 的最大公约数,如下所示。

gcd(54, 20)

54 % 20 = 14

gcd(20,14)

20 % 14 = 6

gcd(14, 6)

14 % 6 = 2

gcd(6, 2)

6%2 = 0

gcd(2, 0)

2


欧几里得算法就是一个重复使用定理 1 的过程,算法在 b=0 时结束。欧几里得算法的优点在于它不断缩小两个数,因此在处理大数时十分快速。欧几里得算法的时间复杂度为 O(lgn),n 为 min(a,b)。

2. 算法代码

欧几里得算法代码如下,代码具有典型递归结构:
#输出a与b的最大公约数
def euclidean(a, b):
    if (b==0):            #边界条件
        return a
    else:
        return euclidean(b,a%b)    #定理14.1
测试代码:

print(euclidean(54,20)

结果如下:

2

3. 算法应用

接下来我们看一下欧几里得算法的几个应用。

1) 解贝祖等式

解贝祖等式(Bezout’s Identity)需要用到扩展欧几里得算法(Extended Euclidean Algorithm)。顾名思义,扩展欧几里得算法由欧几里得算法扩展得来,在原算法的基础上增加了找到整数 m,n,使得贝祖等式 gcd(a, b)=ma+nb 成立的作用。

定理 2:

设 g 为自然数 a 与 b 的最大公约数,m 与 n 为整数常数,a,b 的贝祖等式定义为 g=ma+nb

通过变化 m 和 n 的值,任何正整数组合 (a,b) 都有一个或多个对应的贝祖等式。利用扩展欧几里得算法,我们可以得到任意两个自然数的贝祖等式。来看以下例子,我们将计算 a=54,b=20 时的贝祖等式。

1) 首先我们用欧几里得算法找出 54 与 20 的最大公约数。
gcd(54,20)
=gcd(20,14)
=gcd(14,6)
=gcd(6,2)
=gcd(2,0)
=2

2) 接着我们从后往前推出贝祖等式,永远保持 2 在等式左边。
欧几里得算法

3) 层层递进,保持式子的成立性,并且用上一层的数字代替当前层的数字。
欧几里得算法

4) 最终得到的等式为 3×54−8×20=2,所以 m=3,n=-8。扩展欧几里得算法的规律是用上一行的两个数字代替当前行的第二个数字,比如上例中,用 14 与 6 代替 2,用 20 与 14 代替 6,依此类推。如果再仔细观察一下,会发现当前行的第一个常数是下一行的第二个常数,如图 1 所示。
每一行的第一个常数是下一行的第二个常数
图 1:每一行的第一个常数是下一行的第二个常数

另外,当前行的第二个常数是下一行的第一个常数减去下一行的第二个常数与当前两个数字的整数商,如图 2 所示。
每一行的第二个常数计算方式的规律
图 2:每一行的第二个常数计算方式的规律

读者能写出解扩展欧几里得算法的伪代码吗?可以先自己尝试一下,然后再看以下代码。

扩展欧几里得算法代码如下:
#输出贝祖等式gcd(a,b)=ma+nb的整数常数m,n
def extendedEuclidean(a, b):
    if (b == 0):             #边界条件
        return 1, 0            #返回(1,0)因为g = 1×a +0× 0
    else:
        m,n = extendedEuclidean(b, a % b)   #找到b和r的贝祖等式常数(m,n)
        quotient = a//b           #计算a与b的整数商
        return n, m - (n × quotient)

2) 解线性丢番图方程

欧几里得算法的第二个应用是解线性丢番图方程(Linear Diophantine Equations)。解线性丢番图方程时,我们永远都不会和小数点打交道,方程的系数与解全部限制在整数范围。

定理 3:

a ,b,c 为已知整数,常见的线性丢番图方程的形态为 ax+by=c

例如 120x+16y=4。我们的目的是得出满足条件的所有 (x,y) 整数对。

定理 4:

线性丢番图方程只有两种结果:无解或者无穷解。当 c 是 gcd(a,b) 的倍数时线性丢番图方程有解,否则无解。令 g 为 gcd(a, b),k 为整数常数,我们表达有解条件为 c=kg

证明定理 4:已知 g 是(a, b) 的最大公约数,因此我们肯定 ax+by 可以被 g 整除。又因为 ax+by=c,所以 c 必须也可以被 g 整除,不然方程就无解。比如 5x+6y=2 有解,因为 gcd(5,6)=1,而且 2 是 1 的倍数;8x+6y=3 无解,因为 gcd(8,6)=2,但 2 不是 3 的倍数。

下面我们来解线性丢番图方程。观察例子 8x+15y=2,我们判断方程有解,因为 gcd(8,15)=1,而 2 是 1 的倍数。暂时直接通过观察,我们得出其中一个解 (x1,y1) 为 (4,−2)。

剩余解我们可以直接通过线性丢番图方程解的表达方式,t 为整数常数,直接得出。

回到例子代入公式,得到解为:x=4+15t,y=−2−8t,t∈Z。

然而大多数时候我们不容易直接观察出答案,所以需要采用以下方法:首先解贝祖等式 g=ma+nb,得出 m 与 n,令

例如,在解 8x+15y=2 时,先解 8x+15y=1,得到 (m,n)=(2,−1),然后再接着推算 (x1,y1)=(2× 2,−1× 2)。

总结一下,在解线性丢番图方程时我们有两个地方可能需要用到欧几里得算法。
  1. 判断方程是否有解,即 c 是否为 gcd(a, b) 的倍数,我们需要计算 gcd(a, )b;
  2. 计算 (x1,y1),用扩展欧几里得算法解贝祖等式 g=ma+nb。