前言

之前写过一篇关于求素数算法的博文，其中的算法主要是判断出某个数是不是一个素数。本文就来探讨一下有关素数算法的另一个问题——求最大公约数。本文的内容，来自于笔者最近看的两本书里面提及的算法，一本是《计算机是怎样跑起来的》，还有一本是《什么是数学 · 对思想和方法的基本研究》，下面就用代码来实现一遍书中的算法。

更相减损法

更相减损法出自中国古代的数学专著，其原文是

可半者半之，不可半者，副置分母、子之数，以少减多，更相减损，求其等也。以等数约之。

—— 《九章算术》

这段文言文翻译成白话文，包含了下面四个算法步骤

1. 如果两个数都是偶数，那就先将它们除以二

2. 若两个数都是奇数，那么就比较两者大小，用大的数字作被减数，小的数字作减数

3. 将减法得出的结果与减数作比较，重复上面的步骤

4. 直到最后两个数相等，那么最大公约数就是最后相等的那个数，或最后相等的数乘以每次约去的2最后的乘积

C语言实现

int gcd(int x, int y)
{
    // 任意两个自然数都能同时被1整除
	int gcdNum = 1;
	while (x != y) {
        // 如果两个数同时能被2整除，那么都先除以2
		if (x % 2 == 0 && y % 2 == 0) {
			x /= 2;
			y /= 2;
			gcdNum *= 2;
		}
		else {
			if (x > y)
				x -= y;
			else
				y -= x;
		}
	}
	return gcdNum * x;
}

上述的代码，其实并没有什么优化之处。只要稍作思考就知道，由奇数 - 偶数 = 奇数、偶数 - 奇数 = 奇数这一简单算术规律可知，一旦 x 与 y 其中一个不是偶数了，那么接下来的出现的所有数对都是一奇一偶的形式了，所以也就没有必要每次都去判断 x 和 y 是否同时满足偶数这一条件。

再加上更相减损法中，执行最多的操作是步骤2的减法操作，计算机对于这类计算可以计算地非常快，所以算法可以直接使用一个步骤——减法实现

Python代码实现

def gcd(x, y):
	while x != y:
		if x > y:
			x -= y
		else:
			y -= x
	else:
		return x

欧几里得辗转相除法

两千多年前，欧几里得在其著作《几何原本》中提出了辗转相除法，当时没有那么现代化的数学语言描述，所以他当时用的是几何线段来演示这个算法。

欧几里得的做法是：取两根不同长度的线段，视作两个不同的整数，不断地用其中较短的那条线段，去尽可能地截取较长的那条线段，直到较长的线段被截取成比原来较短线段更短，然后用新的较短的线段去尽可能截取较长的线段，如此重复这一过程，直到较短的线段正好将较长线段截取成整数倍，此时较短的线段长度，就是两者的最大公约数了。

我在维基百科上找到了一张描述这一过程的动态图，很生动直观地体现了这一过程

关于欧几里得辗转相除法的代数证明如下

• 前提一：如果 $a$ 是任一整数，而 $b$ 是任一大于 $0$ 的整数，那么总能找到一个整数 $q$ ，使得等式：$a=bq+r$ 恒成立，其中的整数 $r$ 满足： $0\le r<b$
• 前提二：如果整数 $d$ 是整数 $a$ 和 $b$ 的最大公约数，那么就将这样的关系记作：$d=(a,b)$
• 由前提一可知，形如：$a=bq+r$ 的等式中，可以推出关系：$d=(a，b)=(b，r)$
• 因为对于任意能同时整除 $a$ 和 $b（a>b）$的数 $u$ 来说，有 $a=s\cdot u$ 和 $b=t\cdot u$，此时：$a=bq+r$ $\iff$ $s\cdot u=t\cdot u\cdot q+r$ $\Rightarrow$ $r=(s-qt)\cdot u$ ，因为 $s-qt$ 是一个整数，所以证明 $u$也能整除$r$ ；即：能同时整除 $a$ 和 $b$ 的公因子 $u$ ，也是 $r$ 的因子。反之可证，$b$ 和 $r$ 的公因子，也是 $a$ 的因子。因此，若想求得 $a$ 和 $b$ 的最大公约数，可以转变成求 $a$ 和 $b$ 的余数与 $b$ 的最大公约数。

C语言实现

int gcd(int x, int y)
{
	while (x % y) {
		if (x > y)
			x %= y;
		else
			y %= x;
	}
	return x > y ? y : x;   //返回两个数中最小的那个
}

Python代码实现

def gcd(x, y):
    while x % y:
        if x > y:
            x %= y
        else:
            y %= x
    return min(x, y)   # 返回最后两个数中最小的那个

实际上，辗转相除法是具有递归结构的。其递归结构如下

• 基线条件：若 $a$ 除以 $b$ 没有余数，则 $d=(a,b)=b$，此时 $a$ % $b=0$
• 递归条件：若 $a$ 除以 $b$ 有余数，则 $d=(a,b)=(b,r)=(r,r_2)=(r_2,r_3)=(r_3,r_4)=......$

因此可以给出相应的递归算法的形式

C语言代码实现

int gcd(int x, int y)
{
    return y == 0 ? x : gcd(y, x % y)
}

Python代码实现

def gcd(x, y):
    return x if y == 0 else gcd(y, x%y)

可以看出，使用了递归算法的代码非常简洁优雅！

算法效率比较

两者的比较

• 更相减损法执行过程中，操作最多的运算就是“$-$” ；辗转相除法执行过程中，操作最多的运算是 “%”。对于计算机而言，减法操作要比取余操作更快。
• 更相减损法与辗转相除法都具有递归性，区别在于两者的递归算法中的基线条件不同。更相减损法的基线条件是：两个数相等；而辗转相除法的基线条件是：较大的数正好被较小的数整除。
• 更相减损法每一步的操作，都是在减去一倍的更小的数字，很多时候两个数相差较大时，会不断重复地减去较小的数字，这样的过程其实可以被加快，直接减去一个满足条件的n倍的小数字；辗转相除法每一步的操作会更快，本质上就是去找一个n倍数，然后大数字减去n倍小数字，缩小两数之间的差距，所以辗转相除法的速度瓶颈取决于较小的那个数字。
• 当处理较大整数时，辗转相除法往往比更相减损法更快，原因在于执行的步骤少很多。

若要给辗转相除法给出一个体现算法复杂度的大$O$表示法，说实话，这是一件挺复杂的事情。但与上述分析的一样，辗转相除法的时间复杂度关键在于较小的数字

辗转相除法的计算效率已经被彻底研究过了。一个算法的效率可以用计算所需步数乘以每步计算的开销表示。加百利·拉梅于1884年指出，用辗转相除法计算两个数的最大公约数所需的步数不会超过其中较小数十进制下的位数 h 的5倍。因为每一步的计算开销通常也是h数量级的，所以辗转相除法的复杂度是 h^2 (h的平方)。

—— 维基百科

在维基百科关于辗转相除法的大$O$表示法详细地记载了，其中算法最差情况下的复杂度，平均情况下的复杂度，公式挺复杂的，需要一定的数学水平才能看懂。（我看不懂。。。）

至于，最简单最暴力穷尽的因式分解法求最大公约数，算法太繁琐了，而且效率非常糟糕，所以也就不再讨论了。

相关系列文章：

求解素数的算法总结