[Python] 体验用欧几里得算法计算最大公约数的过程

背景

我们在小学时，学过最大公约数（$\text{Greatest Common Divisor}$）的知识，但是当时并没有学习比较高效的求最大公约数的方法。使用欧几里得算法（$\text{Euclidean Algorithm}$），可以高效地计算出自然数的最大公约数（两个自然数不能同时为 $0$）。本文带您体验这一过程。

正文

根据定义暴力计算

根据最大公约数的定义，我们可以用暴力的方式进行计算自然数 $a$ 和自然数 $b$ 的最大公约数（$a=0$ 和 $b=0$ 不能同时成立）。对应的 $\text{Python}$ 代码如下

def brute_force_gcd(a, b):
    if a == 0 and b == 0:
        raise ValueError("a and b cannot be both 0")
    if a == 0:
        return b
    if b == 0:
        return a
    candidate = min(a, b)
    while True:
        if a % candidate == 0 and b % candidate == 0:
            return candidate
        candidate -= 1

由于所有的自然数都能被 $1$ 整除，所以上述代码求出的最大公约数 $gcd\ge 1$ 总是成立。

观察

如果固定 $b$ 的值，让 $a$ 变化（也可以固定 $a$，让 $b$ 变化，两种方式没有本质区别），我们看一看求出的结果是否有特殊的规律。

$b=1$ 时

当 $b=1$ 时，对任意自然数 $a$，$gcd(a,1)=1$ 显然成立。

$b=2$ 时

当 $b=2$ 时，我们可以列举比较小的 $a$ 值所对应的计算结果

看起来 $gcd(a,2)$ 的值是在 $1$ 和 $2$ 之间交替出现。解释如下 ⬇️

• 当 $a$ 是偶数时，$2\mid 2$ 并且 $2\mid a$，所以 $gcd(a,2)=2$
• 当 $a$ 是奇数时
  • 因为 $2\nmid a$，所以 $gcd(a,2)\ne 2$
  • 那么 $gcd(a,2)=1$

$b=3$ 时

当 $b=3$ 时，我们可以列举比较小的 $a$ 值所对应的计算结果

看起来 $gcd(a,3)$ 的值会以 $3,1,1$ 这样周期反复出现。解释如下 ⬇️

• 当 $a\equiv 0 \pmod 3$ 时，$3|3$ 并且 $3|a$，所以 $gcd(a,3)=3$
• 当 $a\equiv 1 \pmod 3$ 时
  • $3\nmid a$，所以 $gcd(a,3)\ne 3$
  • $2\nmid 3$，所以 $gcd(a,3)\ne 2$
  • 那么 $gcd(a,3)=1$
• 当 $a\equiv 2 \pmod 3$ 时
  • $3\nmid a$，所以 $gcd(a,3)\ne 3$
  • $2\nmid 3$，所以 $gcd(a,3)\ne 2$
  • 那么 $gcd(a,3)=1$

$b=4$

当 $b=4$ 时，我们可以列举比较小的 $a$ 值所对应的计算结果

看起来 $gcd(a,4)$ 的值会以 $4,1,2,1$ 这样周期反复出现。解释如下 ⬇️

• 当 $a\equiv 0 \pmod 4$ 时，$4|4$ 并且 $4|a$，所以 $gcd(a,4)=4$
• 当 $a\equiv 1 \pmod 4$ 时
  • $4\nmid a$，所以 $gcd(a,4)\ne 4$
  • $3\nmid 4$，所以 $gcd(a,4)\ne 3$
  • $2\nmid a$，所以 $gcd(a,4)\ne 2$
  • 那么 $gcd(a,3)=1$
• 当 $a\equiv 2 \pmod 3$ 时
  • $4\nmid a$，所以 $gcd(a,4)\ne 4$
  • $3\nmid 4$，所以 $gcd(a,4)\ne 3$
  • $2\mid 4$，并且 $2\mid a$，所以 $gcd(a,4)= 2$
• 当 $a\equiv 3 \pmod 4$ 时
  • $4\nmid a$，所以 $gcd(a,4)\ne 4$
  • $3\nmid 4$，所以 $gcd(a,4)\ne 3$
  • $2\nmid a$，所以 $gcd(a,4)\ne 2$
  • 那么 $gcd(a,4)=1$

一般的情形

分析了以上几个值比较小的 $b$ 之后，我们尝试从中找规律。假设 $a\ge b$ 成立，看起来以下等式是成立的

我们看看能否证明它。为了方便描述，我们记 $gcd(a,b)=g$, $gcd(a-b,b)=g'$。

我们的目标是证明 $g=g'$ (或者找到两者不相等的反例)。

$gcd(a,b)=g$，根据最大公约数的定义，以下两者成立

• $g\mid a$
• $g\mid b$

所以 $g\mid(a-b)$。$g$ 既是 $a-b$ 的约数，又是 $b$ 的约数，那么 $g$ 是 $a-b$ 和 $b$ 的公约数。按照最大公约数的定义，$g\le g'$ 成立（因为 $g'$ 是 $a-b$ 和 $b$ 的最大公约数）。

反过来看，$gcd(a-b,b)=g'$，根据最大公约数的定义，以下两者成立

• $g'\mid (a-b)$
• $g' \mid b$

所以 $g' \mid ((a-b)+b)$，即 $g' \mid a$。$g'$ 既是 $b$ 的约数，又是 $a$ 的约数，那么 $g'$ 是 $a$ 和 $b$ 的公约数。按照最大公约数的定义，$g'\le g$ 成立（因为 $g$ 是 $a$ 和 $b$ 的最大公约数）。

既然 $g\le g'$ 和 $g'\le g$ 同时成立，那么 $g=g'$ 成立。也就是说，对正整数 $a,b$ 而言，当 $a\ge b$ 成立时，以下等式总是成立。

在此基础上，可以推出 ⬇️ ($a,b$ 都是正整数)

所以对任意正整数 $a,b$ 而言，以下等式成立

这就是欧几里得算法 ($\text{Euclidean Algorithm}$) 的核心思想。

用 $\text{Python}$ 程序实现欧几里得算法

以 $55$ 和 $34$ 为例，我们利用欧几里得算法来计算两者的最大公约数，具体过程如下 ⬇️

用代码来实现，可以这样写 ⬇️

def gcd(a, b):
    if (a, b) == (0, 0):
        raise ValueError("a and b cannot be both 0")
    if b == 0:
        return a
    return gcd(b, a % b)

当我们用这段代码计算 $gcd(55,34)$ 时，函数的调用会形成一个栈。我想用比较直观的方式展示这个函数调用栈，于是在 trae 的帮助下，将其改造成了如下的代码（我觉得不必关心 show_stack 方法的细节，知道它可以展示函数调用栈的内容就够了，show_stack 方法的细节我自己也不懂 😂）

import inspect

def gcd(a, b):
    if (a, b) == (0, 0):
        raise ValueError("a and b cannot be both 0")
    if b == 0:
        show_stack()
        return a
    return gcd(b, a % b)

def show_stack():
    # show_stack 方法的原始内容由 trae 提供，我做了小调整
    for frame_info in inspect.stack()[:][1:]:
        frame = frame_info.frame
        args = []
        for name in frame.f_code.co_varnames[:frame.f_code.co_argcount]:
            if name in frame.f_locals:
                args.append(f"{name}={frame.f_locals[name]}")
        print(f"{frame_info.function}({', '.join(args)})")

if __name__ == "__main__":
    gcd(55, 34)

请将上述代码保存为 calc_gcd.py，使用如下命令可以运行 calc_gcd.py

python3 calc_gcd.py

运行结果如下

gcd(a=1, b=0)
gcd(a=2, b=1)
gcd(a=3, b=2)
gcd(a=5, b=3)
gcd(a=8, b=5)
gcd(a=13, b=8)
gcd(a=21, b=13)
gcd(a=34, b=21)
gcd(a=55, b=34)
<module>()

通过观察上述运行结果，我们可以想象出函数调用栈的结构。

验证

我们可以用如下的代码来简单验证欧几里得算法的计算结果

import matplotlib.pyplot as plt

def gcd(a, b):
    if (a, b) == (0, 0):
        raise ValueError("a and b cannot be both 0")
    if b == 0:
        return a
    return gcd(b, a % b)

def plot_gcd_results(nums, gcd_results, b):
    # plot_gcd_results 里的代码是 trae 帮我生成的，我做了些小调整
    plt.figure(figsize=(12, 8))
    plt.bar(nums, gcd_results, color='skyblue', edgecolor='black')
    plt.xlabel('a', fontsize=12)
    plt.ylabel('GCD(a, %d)' % b, fontsize=12)
    plt.title('GCD(a, %d) for a from %d to %d' % (b, nums[0], nums[-1]), fontsize=14)
    plt.xticks(nums)
    plt.yticks(range(max(gcd_results) + 1))
    plt.grid(axis='y', linestyle='--')
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    b = 4
    nums = range(5 * b + 1)
    gcd_results = [gcd(num, b) for num in nums]
    plot_gcd_results(nums, gcd_results, b)

运行结果如下图所示

可以看到，$gcd(a,4)$ 的值，以 $4$ 为周期的（值是 $4,1,2,1,\cdots$），这与我们前面分析的结果是一致的。

可以将第 $24$ 行的 $b$ 变量改为其他值来进行验证，这里就不赘述了。

参考资料

• A Friendly Introduction to Number Theory 中的
  • 第五章（在 Chapter 1~6 里可以看到第一章到第六章的内容）