欧几里得算法
欧几里得算法是一个用来计算两个整数的最大公因数的算法，最早出现在欧几里得的《几何原本》（wiki语）。
辗转相除法
在天朝，该算法又叫做“辗转相除法”。为了求得两个正整数 $a$ 和 $b$ 的最大公因数，首先用 $b$ 除 $a$ 得到余数 $d_1$，接着将 $b$ 作为被除数、$d_1$ 作为除数，用 $d_1$ 除 $b$ 得到余数 $d_2$，然后分别将 $d_1$ 和 $d_2$ 作为被除数和除数重复此流程，直到得到 $d_{n-1}=kd_n$ 且余数为 $0$，这时候的 $d_n$ 就是整数 $a$ 和 $b$ 的最大公因数。
算法
用数学符号表示辗转相除法的过程，$d_0=b$，$d_i$ 非负且 $d_i\lt d_{i-1}$
\[\begin{aligned} d_1&=a-bq_1\\ d_2&=b-d_1q_2\\ &\cdots\\ d_{i}&=d_{i-2}-d_{i-1}q_n\\ &\cdots\\ d_{n}&=d_{n-2}-d_{n-1}q_n\\ d_{n+1}&=d_{n-1}-d_nq_{n+1} \end{aligned}\]当出现 $d_{n+1}=0$ 时，正整数 $d_n$ 就是 $a$ 和 $b$ 的最大公因数。
证明
首先，因为 $d_i$ 非负且 $d_i\lt d_{i-1}$，
\[0\le d_{n+1}\lt d_{n} \lt d_{n-1} \lt\cdots\lt d_2\lt d_1\]所以 $d_i$ 是一个不断减小的非负整数的序列，最多在 $n+1$ 个步骤之后出现 $d_{n+1}=0$ 而且 $d_{n}\gt 0$。因此该算法总能在有限个步骤内结束。
自底向上看，在第 $n+1$ 步可知 $d_n\mid d_{n-1}$，结合 $d_{n-2}-d_{n-1}q_n=d_{n}$ 可以推出 $d_n\mid d_{n-2}$。同理，可以相继推断出
\[d_n\mid d_{n-3}\text{, }d_n\mid d_{n-4}\text{, }\cdots\text{, }d_n\mid d_{2}\text{, }d_n\mid d_{1}\]所以 $d_n\mid b$ 且 $d_n\mid a$，因此 $d_n$ 是 $a$ 和 $b$ 的公因数，即
\[d_n\le\gcd(a,b)\tag{1}\]自顶向下看，在第一步可知 $\gcd(a,b)\mid d_1$，结合第二步 $b-d_1q_2=d_2$ 可知 $\gcd(a,b)\mid d_2$。同理，可以相继推断出
\[\gcd(a,b)\mid d_3\text{, }\cdots\text{, }\gcd(a,b)\mid d_{n-2}\text{, }\gcd(a,b)\mid d_{n-1}\]所以 $\gcd(a,b)\mid d_{n}$；加上 $d_{n}\gt0$，因此 $\gcd(a,b)\le d_{n}$。结合 $(1)$ 式，得到
\[d_n=\gcd(a,b)\tag{Q.E.D}\]这个算法可以用非常简单，可以用以下代码实现：
def gcd(a, b):
while True:
d = a % b
if d == 0:
return b
a, b = b, d
顺便提一个点就是代码无需比较 $a$ 和 $b$ 的大小并且调换被除数和除数的顺序。即使 $a\lt b$，在第二轮里数值较大的 $b$ 就会变成被除数、上一轮较小的 $a$ 会变成除数，跟调换顺序之后相比只多了一次除法而已。
扩展欧几里得算法
在欧几里得算法的基础上还可以做进一步推断：存在整数 $s$ 和 $t$ 满足
\[as+bt=\gcd(a,b)\tag{2}\]整数 $s$ 和 $t$ 可以使用改进之后的欧几里得算法计算得到，这里先从证明入手。首先，欧几里得算法中每一个步骤中的 $d_i$ 都可以写成 $a$ 和 $b$ 的多项式
\[as_i+bt_i=d_i\tag{3}\]这个结果不是瞎猜的：观察辗转相除法的前两条等式
\[\begin{aligned} d_1&=a\times1-b\times q_1\\ d_2&=b-d_1q_2=a\times(-1)+b\times(1+q_1q_2)\\ \end{aligned}\]将上述两个等式代入 $d_3=d_1-d_2q_3$
\[\begin{aligned} d_3&=(a-bq_1)-[a\times(-1)+b\times(1+q_1q_2)]q_3\\ &=a\times(1+q_3)+b\times[-q_1-(1+q_1q_2)q_3] \end{aligned}\]所以 $d_3$ 也可以用关于 $a$ 和 $b$ 的多项式表示。实际上
证明
归纳法证明 先来看trivial case。当 $i=1$ 时，$s_1=1$ 且 $t_1=-q_1$
\[a-bq_1=a\times1+b\times(-q_i)=d_1\]当 $i=2$ 时，将 $d_1=a-bq_1$ 代入 $b-d_1q_2=d_2$
\[b-(a-bq_1)q_2=a(-q_2)+b(1+q_1q_2)=d_2\]所以 $s_2=-q_2$ 且 $t_2=1+q_1q_2$，当 $i=1$ 或者 $i=2$ 时，$(3)$ 式成立。
假设当 $i\lt k$ 时，$(3)$ 式成立，换言之 $i=k-2$ 或 $i=k-1$ 时，存在以下等式
\[\begin{aligned} as_{k-2}+bt_{k-2}&=d_{k-2}\\ as_{k-1}+bt_{k-1}&=d_{k-1}\\ \end{aligned}\]将上述等式代入 $b_{k-2}-b_{k-1}q_{k}=d_{k}$ 可以得到
\[\begin{aligned} d_{k}&=(as_{k-2}+bt_{k-2})-(as_{k-1}+bt_{k-1})q_{k}\\ &=a(s_{k-2}-s_{k-1}q_k)+b(t_{k-2}-t_{k-1}q_k) \end{aligned}\]令 $s_k=s_{k-2}-s_{k-1}q_k$，$t_k=t_{k-2}-t_{k-1}q_k$，因此 $d_k=as_k+bt_k$。这说明 $d_k$ 也可以写成关于 $a$ 和 $b$ 的多项式的形式，原假设 $(3)$ 对于 $i=k$ 也成立。换句话说，最大公因数 $d_n=\gcd(a,b)$ 也可以用 $a$ 和 $b$ 的多项式来表达
\[as_n+bt_n=d_n\tag{Q.E.D}\]代码实现
上述归纳证明还给出了 $s_i$ 和 $t_i$ 的递推计算方式
\[\begin{aligned} s_i&=s_{i-2}-s_{i-1}q_i\\ t_i&=t_{i-2}-t_{i-1}q_i\\ d_i&=d_{i-2}-d_{i-1}q_i\\ \end{aligned}\]且 $i\ge2$。这里为了实现方便可以对 $i\lt2$ 的情况可以稍加修改，根据 $s_1$ 和 $s_2$ 逆推 $s_0$ 的值，${t_0}$ 同理
\[\begin{cases} s_0=0\\ t_0=1 \end{cases} \text{ and } \begin{cases} s_1=1\\ t_1=-q_1&(a=bq_1+d_1)\\ \end{cases}\]实际上用代码实现的时候只需要计算 $t$ 或 $s$ 当中的其中之一就可以了，另一个可以通过 $a$，$b$ 以及 $\gcd(a,b)$ 计算出来。
def euclid_ex(a,b):
to, tn, ao, bo = 0, 1, a, b
while True:
d = a % b
if d == 0:
return (bo - bo * tn) // ao, tn, b
q = a // b
to, tn = tn, to - tn * q
a, b = b, d
贝祖定理
前面的扩展欧几里得算法提到的 $(2)$ 式：存在整数 $s$ 和 $t$ 满足
\[as+bt=\gcd(a,b)\tag{2}\]实际上，能满足 $as+bt$ 形式的最小整数就是 $\gcd(a,b)$，是为贝祖定理。
这是显而易见的。假设存在 $d\lt\gcd(a,b)$ 使得 $as’+bt’=d$，由于 $\gcd(a,b)\mid a$ 且 $\gcd(a,b)\mid b$，所以 $\gcd(a,b)\mid d\implies d\ge\gcd(a,b)$。这个结果和假设矛盾，所以原假设不能成立，也就是说不存在小于 $\gcd(a,b)$ 且满足 $(2)$ 的整数。
特别地，如果 $a$ 和 $b$ 互质，那么存在 $s$ 和 $t$ 使得
\[ax+by=1\]这个推论反过来也是成立的：对于 $a$ 和 $b$，如果存在整数 $s$ 和 $t$ 满足 $as+bt=1$，那么 $a$ 和 $b$ 互质。这是两个整数互质的充分必要条件。
算法效率
这部分内容参考了 A Concrete Introduction to Higher Algebra 第一部分第三章F小节 The Efficiency of Euclid’s Algorithm
假如用函数 $N(a, b)$ $(a\gt b)$ 表示欧几里得算法消耗的时间，而且使用一次除法同时计算出商和余数要消耗的时间为常数 $1$，那么
\[\begin{aligned} N(a,b)&=N(b,d_1)+1\\ &=[N(d_1,d_2)+1]+1\\ &=N(d_{n-1},d_n)+n \end{aligned}\]最后一步计算出余数 $d_{n+1}$，由于下一步 $d_{n+1}$ 整除 $d_n$ 没有余数，这里就不再视为消耗 $1$ 个时间。于是 $N(a,b)=n$，$n$ 的大小跟余数的序列有关
\[a,b,d_1,d_2,\dots,d_{n-1},d_{n}\]因为余数总是小于除数，所以这个序列本身是递减的。$n$ 的值和 $d_i$ 下降的“速度”有关：如果两个余数的差 $\Delta=d_{i-1}-d_i$ 较大，下一轮除法的商
\[\begin{aligned} q_{i+1}&=\frac{d_{i-1}-d_{i+1}}{d_i}\\ &=\frac{d_{i-1}}{d_i}-\frac{d_{i+1}}{d_i}\gt\frac{d_{i-1}}{d_i}-1 \end{aligned}\]也会是一个相对较大的值，所以，商的数值越大意味着余数下降得越“快”，$N(a,b)=n$ 的值越小；反之，商越小则 $N(a,b)=n$ 越小。以 $a=89$ 和 $b=55$ 为例子：
\[\begin{aligned} 89&=55\times1+34\\ 55&=34\times1+21\\ 34&=21\times1+13\\ 21&=13\times1+8\\ 13&=8\times1+5\\ 8&=5\times1+3\\ 5&=3\times1+2\\ 3&=2\times1+1\\ 2&=1\times2 \end{aligned}\]除去最后一步，这里使用了 $8$ 次计算余数的除法，每一步的商都不超过 $1$。仔细看每一步的余数，按照反过来的顺序其实是斐波那契数列的一部分子序列：
\[1,1,2,3,5,8,13,21,34,55,89,\dots\] \[F_{n}=\begin{cases} F_{n-1}+F_{n-2}&(n\ge3)\\ 1&(n=1\text{ or }n=2) \end{cases}\]如果对两个相邻的斐波那契数 $F_{n+1}$ 和 $F_{n}$ 使用欧几里得算法，那么除开最后一步以外，每一步除法的商都是 $1$，余数将会是 $F_{n-1}$，一共需要 $(n-1)-1=n-2$ 次除法。例子 $a=F_{11}=89$ 和 $b=F_{10}=55$ 使用了 $8$ 次除法。再来看另一个例子 $a’=F_{11}=89$ 和 $b’=53$
\[\begin{aligned} 89&=53\times1+36\\ 53&=36\times1+17\\ 36&=17\times2+2\\ 17&=2\times8+1\\ 2&=1\times2 \end{aligned}\]这个例子 $N(89, 53)=5$，如果商越大，算法使用的除法次数越少；欧几里得算法对任意两个连续的斐波那契数使用的除法次数，会比其他同样规模的两个整数要多。
Lamé’s Theorem
定理 对于任意两个整数 $a\gt b$，如果 $b\lt F_n$，那么
\[N(a,b)\le n-3\lt n-2=N(F_{n+1},F_{n})\]这个定理以法国人 Gabriel Lamé 的名字命名，他以针对欧几里得算法的运行时间分析而闻名，这个分析被认为是计算复杂性理论的开端。奇怪的是，英文wiki有记录 Lamé 本人生平的页面，却没有这个定理的页面😅。作为一道练习题(Exercise 69)，书里并没有证明，但是给出了可使用归纳法证明的提示：
归纳法证明 先来看trivial case。当 $n=3$ ，$F(4)=3$，$F(3)=2$，$N(F_4,F_3)=1$。对于任意 $b\lt F_3=2\implies b=1$，因此 $N(a,1)=0$，命题对于 $n=3$ 成立。
假设对于 $n\le k$ 时成立，也就是说对于 $b\lt F_k$
\[N(a,b)\le k-3\lt k-2=N(F_{k+1},F_{k})\tag{4}\]然后考虑 $n={k+1}$ 时假设是否能成立。用 $b$ 除 $a$ 可以得到 \(a=bq+r\text{ and }r\lt b\)
所以 $N(a,b)=N(b,r)+1$；其次又因为 $b\lt F_{k+1}$，所以 $r\lt F_{k+1}$。这里需要对 $r$ 分两种情况考虑：$r\lt F_k$ 或 $F_k\le r\lt F_{k+1}$。
首先假如 $r\lt F_k$，可以应用归纳假设 $(4)$ 得到
\[N(b,r)\le k-3\lt k-2=N(F_{k+1},F_k)\]所以
\[\begin{aligned} N(a,b)&=N(b,r)+1\\ &\le (k+1)-3\\ &\lt (k+1)-2 = N(F_{k+2}, F_{k+1}) \end{aligned}\]其次针对第二种情况 $F_k\le r\lt b\lt F_{k+1}$。用 $r$ 除 $b$ 得到
\[b=rq'+r'\text{ and }r'\lt r\]所以
\[N(a,b)=N(b,r)+1=N(r,r')+2\tag{5}\]考虑余数 $r’$ 的大小。因为 $r\lt b\implies q’\ge1$，所以
\[r'=b-rq'\le b-r\lt F_{k+1}-F_{k}=F_{k-1}\]所以对 $r’\lt F_{k-1}$ 应用归纳假设 $(4)$ 可知
\[N(r,r')\le (k-1)-3\lt (k-1)-2=N(F_{k},F_{k-1})\]又因为 $(5)$，所以
\[\begin{aligned} N(a,b)&=N(r,r')+2\\ &\le (k-1)-3+2=(k+1)-3\\ &\lt (k+1)-2\\ &=N(F_{k+2},F_{k+1}) \end{aligned}\]所以归纳假设 $(4)$ 对于 $n=k+1$ 也成立。综上所述，对于任意两个整数 $a\gt b$ 和 $b\lt F_n$，都有
\[N(a,b)\le n-3\lt n-2=N(F_{n+1},F_{n})\tag{Q.E.D}\]最坏情况
上面的 Lamé’s Theorem 证明了算法的最坏情况出现在两个斐波那契数之间。为了能用 $O$ 记号描述最坏情况，需要搞清楚输入规模和执行时间的上界 $N(F_{n+1}, F_{n})$ 的关系。先来看下面两个推论：
推论 $1$：$F_{n+5}\gt10F_{n}$
这个推论表示每相隔 $5$ 个斐波那契数、十进制数位就会增加 $1$ 。利用递推关系展开 $F_{n+5}$
\[\begin{aligned} F_{n+5}&=F_{n+4}+F_{n+3}\\ &=2F_{n+3}+F_{n+2}\\ &=3F_{n+2}+2F_{n+1}\\ &=5F_{n+1}+3F_{n}\\ &=8F_{n}+5F_{n-1}\\ \end{aligned}\]因此 $F_{n+5}/F_n=8+5F_{n-1}/F_n$。假如 $5F_{n-1}/F_n\gt2$ 为真，推论 $F_{n+5}\gt10F_{n}$ 就可以得到证明了。实际上，用归纳法证明$5F_{n-1}\gt2F_n$的过程是十分简单的：
归纳法证明 当 $n=2$ 时，$F_{n-1}=F_1=1$，$F_{n}=F_2=2$，因此 $5F_{n-1}=5F_1=5\gt2=2F_2=2F_{n}$，假设成立。
假设 $5F_{n-1}\gt2F_n$ 对 $n=k$ 成立，即
\[5F_{k-1}\gt2F_k\]那么当 $n=k+1$ 时，
\[\begin{aligned} 5F_k&=5F_{k-1}+5F_{k-2}&(\text{定义})\\ &\gt2F_k+2F_{k-1}&(\text{归纳假设})\\ &=2F_{k+1}&(\text{定义}) \end{aligned}\]所以归纳假设对 $n=k+1$ 也成立，这表示 $5F_{n-1}\gt2F_n$，因此 $F_{n+5}/F_n=8+5F_{n-1}/F_n\gt10$，推论得证。
推论 $2$：$F_{5d+2}$ 的十进制数位至少为 $d+1$
由推论 $1$ 可以知道 $F_{5d+2}/F_{5(d-1)+2}\gt10^1$，$F_{5d+2}/F_{2}\gt10^d$。又因为 $F_2=1$，所以 $F_{5d+2}\gt10^d$，这表示 $F_{5d+2}$ 的十进制数位至少为 $d+1$。
结论
推论 $2$ 可以换个说法：所有整数 $b\lt F_{5d+2}$ 的十进制数位不超过 $d+1$。加上前面提到的Lamé定理的结果，可以知道：
对于任何十进制数位不超过 $d$ 的整数 $b$，也就是 $b\lt 10^d\lt F_{5d+2}$，欧几里得算法所需要的除法次数（或者说消耗时间）$N(a,b)$ 的上界 $N(a, F_{5d+2})=5d$。
这个上界是 $b$ 的十进制数位 $d$ 的 $5$ 倍。用一个例子来说明，对 $b=10^2$ 以内的整数使用欧几里得算法，除法次数的上界是 $N(F_{13}, F_{12})=10$。换句话说，如果输入的其中一个整数是 $n$ ，算法的最坏情况的时间复杂度是 $O(\log n)$。
平均情况
书里提到欧几里得算法在平均情况下使用的除法次数是 $2d$，据说这个结论在 Knuth 的大作 The Art of Computer Programming (4.5.3) 也有，后面有时间的话也研究下 🙏。
