辗转相除法

自从转载了来自圣经的算法一文后，就想把这些算法比较详细地搞清楚，先拿辗转相除法开刀了，谁让她最简单呢。呵呵。

　　下面的大部分内容来自维基百科。

　　辗转相除法，又被称为欧几里德（Euclidean）算法， 是求最大公约数的算法。辗转相除法首次出现于欧几里得的《几何原本》（第VII卷，命题i和ii）中，而在中国则可以追溯至东汉出现的《九章算术》。

　　两个数的最大公约数是指能同时整除它们的最大正整数。辗转相除法的基本原理是：两个数的最大公约数等于它们中较小的数和两数之差的最大公约数。例如，252和105的最大公约数是21（252 = 21 × 12；105 = 21 × 5）；因为252 − 105 = 147，所以147和105的最大公约数也是21。在这个过程中，较大的数缩小了，所以继续进行同样的计算可以不断缩小这两个数直至其中一个变成零。这时，所剩下的还没有变成零的数就是两数的最大公约数。由辗转相除法也可以推出，两数的最大公约数可以用两数的整数倍相加来表示，如21 = 5 × 105 + (−2) × 252。这个重要的等式叫做贝祖等式（Bézout's identity）。

　　辗转相除法法原先只用来处理自然数，但在19世纪，辗转相除法被推广至其他类型的数，如高斯整数和一元多项式。另外，还被用来解决丢番图方程（Diophantine equations）和构造连分数等。

算法描述

　　两个数a，b的最大公约数记为GCD(a,b)。a，b的最大公约数是两个数的公共素因子的乘积。如462可以分解成2 × 3 × 7 × 11；1071可以分解成3 × 3 × 7 × 17。462和1071的最大公约数等于它们共有的素因数的乘积3 × 7 = 21。如果两数没有公共的素因数，那么它们的最大公约数是1，也即这两个数互素，即GCD(a,b)=1。另g=GCD(a,b)，则a=mg, b=ng，其中m，n均为正整数。由上述分析可知，m，n互素。因为m，n没有公共素因子，GCD(m,n)=1。

　　辗转相除法是一种递归算法。设k表示步骤数（从0开始计数），算法的计算过程如下。每一步的输入是都是前两次计算的余数r_k−1和r_k−2。因为每一步计算出的余数都在不断减小，所以，r_k−1小于r_k−2。在第k步中，算法计算出满足以下等式的商q_k和余数 r_k：

r_k−2 = q_k r_k−1 + r_k

其中r_k < r_k−1。也就是r_k−2要不断减去r_k−1直到比r_k−1小。

　　在第一步计算时（k = 0），设r₋₂和r₋₁分别等于a和b，第2步（此时k = 1）时计算r₋₁（即b）和r₀（第一步计算产生的余数）相除产生的商和余数，以此类推。整个算法可以用如下等式表示：

a = q₀ b + r₀

b = q₁ r₀ + r₁

r₀ = q₂ r₁ + r₂

r₁ = q₃ r₂ + r₃

…

　　如果输入参数a小于b，则第一步计算的结果是交换两个变量的值：因为a < b，所以a和b相除得到的商q₀等于0，余数r₀等于a。所以在运算的每一步中得出的余数一定小于上一步计算的余数（r_k一定小于r_k−1）。由于每一步的余数都在减小并且不为负数，必然存在第N步时r_N等于0，使算法终止，r_N−1 就是a和b的最大公约数。其中N不可能无穷大，因为在r₀和0之间只有有限个自然数。

　　辗转相除法的正确性可以用两步来证明。首先，算法的最终结果r_N−1同时整除a和b。因为它是一个公约数，所以必然小于或者等于最大公约数g。然后，任何a和b的公约数都能整除r_N−1。所以g一定小于或等于r_N−1。两个不等式只在r_N−1 = g是同时成立。

算法实现

　　递归版本

function　gcd(a, b)
2 　　if a<b
3 　　　　swap(a,b);　　
4 　　if b==0
5 　　　　then return a;
6 　　else
7 　　　　return gcd(b, a mod b);
8 end

　　循环版本

function gcd(a,b)
2 　　if a<b
3 　　　　then swap(a,b);
4 　　while(b!=0)
5 　　{
6 　　　　c = a mod b;
7 　　　　a = b;
8 　　　　b = c;
9 　　}
10 　　return a;
11 end

扩展欧几里德算法

　　贝祖等式说明，两个数a和b的最大公约数g可以表示为a和b的线性和。也就是说，存在整数s和t使g = sa + tb。^{整数s和t可以从辗转相除法算出的商q₀、q₁……计算出。贝祖等式的整数s和t可以通过扩展欧几里得算法算出。这个扩展算法在原有辗转相除法的基础上增加了两个递归等式：}

s_k = s_k−2 − q_k−1s_k−1

t_k = t_k−2 − q_k−1t_k−1

算法开始时：

s₋₂ = 1, t₋₂ = 0

s₋₁ = 0, t₋₁ = 1

加上这个递归式后，当算法终止于r_N = 0，贝祖等式的整数s和t分别由s_N和t_N给出。

这个算法的正确性可以用数学归纳法来证明。假设递归至第k−1步是正确的，也就是假设：

r_j = s_j a + t_j b

对所有j小于k成立。则第k步运算得出以下等式：

r_k = r_k−2 − q_k−1r_k−1

因为r_k−2和r_k−1被假定是正确的，所以可以用s和t表示：

r_k = (s_k−2 a + t_k−2 b) − q_k−1(s_k−1 a + t_k−1 b)

整理后得到第k步的结果，和我们期望得到的结果一致：

r_k = s_k a + t_k b = (s_k−2 − q_k−1s_k−1) a + (t_k−2 − q_k−1t_k−1) b

算法分析

　　加百利·拉梅于1884年指出，一个算法的效率可以用计算所需步数乘以每步计算的开销表示。用辗转相除法计算两个数的最大公约数所需的步数不会超过其中较小数的位数h的5倍。因为每一步的计算开销通常也是h数量级的，所以辗转相除法的复杂度是h²。

计算步数

计算两个自然数a和b的最大公约数所需的步数可以表示为T(a, b)。如果a和b的最大公约数是g，m和n是两个互素整数，使a = mg，b = ng，那么：

T(a, b) = T(m, n)

这可以通过在辗转相除法的计算过程中的每一步都除以g来证明。同样，当a和b同时乘以w时，计算步数不变：T(a, b) = T(wa, wb)。所以，对于数值上相近的数，如T(a, b)和T(a, b + 1)，计算步数可能相差很大。

根据辗转相除法的递归性质可以得出另一个公式：

T(a, b) = 1 + T(b, r₀) = 2 + T(r₀, r₁) = … = N + T(r_N−2, r_N−1) = N + 1

并且定义T(x, 0) = 0。

最差情况

假设用辗转相除法求自然数a和b（a > b > 0）的最大公约数需要N步，那么满足这一条件的a和b的最小值是分别是斐波那契数F_N+2和F_N+1。这可以用数学归纳法证明。假设N = 1，b整除a，满足这一条件的a和b最小是b = 1、a = 2，正是F₂和F₃。现在假设这一规律对M − 1有效。一个需要M步的算法的第一步是a = q₀b + r₀，第二步是b = q₁r₀ + r₁。因为算法是递归的，它需要M − 1 步才能算出 GCD(b, r₀)，其中b和r₀的最小值是 F_M+1 和 F_M。所以a的最小值是当q₀ = 1 的时候，此时 a = b + r₀ = F_M+1 + F_M = F_M+2。1844年，加百利·拉梅发现这个证明标志着计算复杂性理论的诞生。这也是斐波那契数列的第一个实际应用。

这个结果也证明了辗转相除法的运算步骤不会超过较小数十进制下的位数的五倍。因为如果算法需要N步，那么b一定大于或等于F_N+1，也就是一定大于或等于φ^N−1，其中φ是黄金分割比。因为b ≥ φ^N−1，所以N − 1 ≤ log_φb。因为log₁₀φ > 1/5，(N − 1)/5 < log₁₀φ log_φb = log₁₀b，所以N ≤ 5 log₁₀b。所以，辗转相除法不会进行超过O(h)次除法，其中h是较小数b在十进制下的位数。

平均情况

辗转相除法的平均步骤数有三种不同的定义。第一种定义是计算已知自然数a和从0到a − 1范围内随机选取的自然数b的最大公约数所需的时间T(a)：

但是因为T(a, b)在连续整数间变化非常剧烈，所以T(a)的平均值也会显得很杂乱。

为了解决这个问题，第二种定义规定τ(a)只要取遍其中所有和a互素的数即可:

在小于a的数中，有φ(a)个数与a互素，其中φ是欧拉函数。在这个定义中，τ(a)的函数值增长地平稳很多。

但是这个公式中仍然存在一个问题：ε是无穷小量。公式中的常熟C等于：

其中γ是欧拉-马歇罗尼常数，ζ′是黎曼ζ函数的导数。公式最左边的由两个独立的方法确定。

因为第一种定义可以通过用第二种定义的求和来完成：

所以也可以由以下公式近似：

其中Λ(d)是冯·曼戈尔特函数。

第三种定义Y(n)定义为从1到n间随机选取a和b时计算他们的最大公约数所需的平均步骤数：

将T(a)的近似公式代入，得到Y(n)的近似值：

每一步的计算开销

在辗转相除法的每一步中，商q_k和余数r_k都通过r_k−2和r_k−1求出：

r_k−2 = q_k r_k−1 + r_k.

所以每一步的计算开销主要与计算商q_k的算法有关，因为余数r_k可以很迅速地从r_k−2、r_k−1和q_k计算出来：

r_k = r_k−2 − q_k r_k−1.

而计算一个h为整数的除法的算法复杂度是O(h(ℓ+1))，其中ℓ是商的位数。

作为对比，辗转相除法原先的版本使用的是减法，因此效率要慢很多。进行一次除法等同于进行q次减法（其中q是商）。如果a和b的比很大，计算出的商也很大，也就需要进行很多次减法。但在另一方面，计算出来的商在大多数情况下都是非常小的，除法中得出一个确定的商q的概率大约是。其中u = (q + 1)²。比如，商是1、2、3、4的可能性分别是大约41.5%、17.0%、9.3%、5.9%。因为计算机计算减法要快于除法，所以对于很大的数字，基于减法的辗转相除法的性能可以比得上基于除法的算法。这也被运用于二进制最大公约数算法。

综合考虑算法需要的步数和每一步的计算开销，辗转相除法的随两个数字a和b的平均位数成平方级的速度增长(h²)。设h₀、h₁、…、h_N−1表示计算过程中的余数r₀、r₁、…、r_N−1的位数，因为算法的步数N随h线性增长，所以算法的运算时间为：