密码学数学基础 笔记 2：欧几里得算法及其变形、Bézout 等式，以及线性丢番图方程的求解

欧几里得算法

最大公因数与辗转相除法

设整数 $a,b$ 不同时为 $0$。它们的最大公因数记为 $\gcd(a,b)$，是二者最大的正公因子

欧几里得算法还会证明一个更强的性质：$a,b$ 的每一个公因子都整除 $\gcd(a,b)$

一步变换为什么成立

若

$$a=qb+r,$$

则 $a,b$ 与 $b,r$ 有完全相同的公因子：

• 若 $d\mid a$ 且 $d\mid b$，则 $d\mid (a-qb=r)$
• 若 $d\mid b$ 且 $d\mid r$，则 $d\mid (qb+r=a)$

因此

$$\boxed{\gcd(a,b)=\gcd(b,r)}$$

每次做带余除法，并用“除数、余数”替换原来的两个数，最大公因数就保持不变

算法与例子

以 $252,198$ 为例：

$$\begin{aligned} 252&=1\times198+54,\\ 198&=3\times54+36,\\ 54&=1\times36+18,\\ 36&=2\times18+0 \end{aligned}$$

因此

$$\gcd(252,198)=\gcd(198,54)=\gcd(54,36) =\gcd(36,18)=18$$

每一步的非零余数都比除数小，所以过程必然终止。到达 $(d,0)$ 时，最大公因数就是 $d$，即最后一个非零余数。若初始就有一个数为 $0$，直接返回另一个数的绝对值

更进一步，因为整个过程中公因子的集合保持不变，而 $(d,0)$ 的公因子恰好是 $d$ 的因子，所以原来每个公因子都整除 $d$

1
2
3
4
5

def gcd(a, b):
    a, b = abs(a), abs(b)
    while b:
        a, b = b, a % b
    return a

这段实现也采用常见的约定 $\gcd(0,0)=0$

算法需要多少次除法

先假设 $a>b>0$。把余数依次记为

$$r_{-1}=a,\qquad r_0=b,\qquad r_{i+1}=r_{i-1}\bmod r_i$$

仅由“余数严格递减”，可以得到至多 $b$ 次除法的粗略上界。更好的观察是：每两步，余数至少缩小一半

考虑相邻的 $r_i,r_{i+1}$，在下一步仍能执行时：

• 若 $r_{i+1}\leq r_i/2$，则 $r_{i+2}<r_{i+1}\leq r_i/2$
• 若 $r_{i+1}>r_i/2$，则 $r_i$ 除以 $r_{i+1}$ 的商为 $1$，于是 $r_{i+2}=r_i-r_{i+1}<r_i/2$

因此

$$r_{i+2}<\frac{r_i}{2}$$

从 $b$ 开始，连续减半 $O(\log b)$ 次后就会降到 $1$ 以下，算法必须终止。一个宽松的除法次数上界为

$$2\log_2 b+2$$

这里统计的是整数除法的次数。若输入有 $L$ 个二进制位，那么除法次数为 $O(L)$；大整数除法本身还需要计算成本，不能把它当成一次固定时间的操作

两种变形

最小余数版本

等式 $\gcd(a,b)=\gcd(b,r)$ 只要求 $a=qb+r$，对余数的正负没有要求。又因为

$$\gcd(b,r)=\gcd(b,-r),$$

所以可以选择绝对值不超过 $b/2$ 的余数，再取其绝对值继续计算

例如

$$\begin{aligned} 123987&=4\times34762-15061,\\ 34762&=2\times15061+4640,\\ 15061&=3\times4640+1141,\\ 4640&=4\times1141+76,\\ 1141&=15\times76+1,\\ 76&=76\times1+0 \end{aligned}$$

因此 $\gcd(123987,34762)=1$

实现时，先求普通余数 $r=a\bmod b$。如果 $2r>b$，就用 $b-r$ 替换它。这相当于选择负余数 $r-b$ 后取绝对值

1
2
3
4
5
6
7
8

def balanced_gcd(a, b):
    a, b = abs(a), abs(b)
    while b:
        r = a % b
        if 2 * r > b:
            r = b - r
        a, b = b, r
    return a

此时每一步的新余数都不超过当前除数的一半。对 $a>b>0$，除法次数至多为

$$\operatorname{bits}(b)=\lfloor\log_2 b\rfloor+1$$

这个上界比普通版本更小，但实际运行时间还取决于额外判断和具体的大整数实现

二进制 GCD

二进制 GCD 根据奇偶性，用除以 $2$ 和减法缩小输入。对非负整数 $u,v$，基本规则为：

例如，两个数都为奇数时，先由 $\gcd(u,v)=\gcd(u-v,v)$ 做减法；由于 $u-v$ 为偶数、$v$ 为奇数，再去掉前者的因子 $2$

这些版本共同依赖两点：变换保持最大公因数不变，并且不断减小待处理的整数

Bézout 等式与扩展欧几里得算法

普通算法只返回最大公因数。如果在计算过程中同时记录每个余数如何由原来的 $a,b$ 线性组合得到，就能找到整数 $u,v$，使得

$$\boxed{au+bv=\gcd(a,b)}$$

这就是 Bézout 等式

从余数反向代入

继续使用 $252,198$ 的例子。从最后一个非零余数出发：

$$\begin{aligned} 18&=54-36\\ &=54-(198-3\times54)\\ &=4\times54-198\\ &=4(252-198)-198\\ &=4\times252-5\times198 \end{aligned}$$

因此一组 Bézout 系数为 $(u,v)=(4,-5)$

一般地，若 $a=qb+r$，递归求得

$$d=u'b+v'r,$$

将 $r=a-qb$ 代回，就有

$$d=v'a+(u'-qv')b$$

因此，扩展算法的系数递推为

$$(u,v)=(v',u'-qv')$$

正向维护系数

也可以在求余数时同步更新系数，避免最后反向代入。始终维护

$$r_i=u_i a+v_i b$$

因为

$$r_{i+1}=r_{i-1}-q_{i+1}r_i,$$

对应的系数满足同样的递推：

$$\begin{aligned} u_{i+1}&=u_{i-1}-q_{i+1}u_i,\\ v_{i+1}&=v_{i-1}-q_{i+1}v_i \end{aligned}$$

初始时，$a=1a+0b$，$b=0a+1b$。前面的例子可以整理为：

对应的迭代实现如下，返回 $(d,u,v)$，满足 $d\geq0$ 且 $au+bv=d$：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

def extended_gcd(a, b):
    old_r, r = abs(a), abs(b)
    old_u, u = 1, 0
    old_v, v = 0, 1

    while r:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_u, u = u, old_u - q * u
        old_v, v = v, old_v - q * v

    if a < 0:
        old_u = -old_u
    if b < 0:
        old_v = -old_v

    return old_r, old_u, old_v

Python 的多重赋值先计算右侧所有表达式，再更新左侧变量，因此每次递推使用的都是更新前的系数

矩阵形式

一次欧几里得变换也可以写成

$$(b,r)=(a,b) \begin{pmatrix} 0&1\\ 1&-q \end{pmatrix}$$

设每一步的变换矩阵为 $M_i$，算法结束时有

$$(d,0)=(a,b)M_1M_2\cdots M_t$$

如果把矩阵乘积写成

$$M_1M_2\cdots M_t= \begin{pmatrix} u&e\\ v&f \end{pmatrix},$$

就得到 $d=au+bv$，所以第一列正是 Bézout 系数

例如，对 $15,6$，两次除法的商都是 $2$，因此

$$\begin{pmatrix}0&1\\1&-2\end{pmatrix}^2 =\begin{pmatrix}1&-2\\-2&5\end{pmatrix}, \qquad 3=15-2\times6$$

矩阵形式与前面的系数递推描述的是同一个计算过程

线性丢番图方程

丢番图方程要求未知数取整数值。考虑最基本的情形

$$ax+by=c,$$

下面设 $a,b$ 均非零，并记 $d=\gcd(a,b)$

有解的条件

方程有整数解，当且仅当

$$\boxed{d\mid c}$$

必要性来自整除的性质：$d$ 整除 $a,b$，因此整除任何 $ax+by$

充分性来自 Bézout 等式。若 $au+bv=d$，且 $c=kd$，则

$$a(ku)+b(kv)=c$$

所以一组特解为

$$x_0=\frac cd u,\qquad y_0=\frac cd v$$

例如，$\gcd(713,851)=23$，但 $23\nmid10$，所以 $713x+851y=10$ 无整数解

全部整数解

若 $(x_0,y_0)$ 是一组特解，则全部整数解为

$$\boxed{ x=x_0+\frac bd t,\qquad y=y_0-\frac ad t,\qquad t\in\mathbb{Z} }$$

代回原式，新增的两项相互抵消，因此这些确实都是解

再说明没有遗漏。记 $A=a/d$、$B=b/d$，由 Bézout 等式可得某些整数 $u,v$ 满足 $Au+Bv=1$，因此 $A,B$ 互素

任意另一组解与特解相减，得到

$$A(x-x_0)+B(y-y_0)=0$$

于是 $B\mid A(x-x_0)$。将 $Au+Bv=1$ 乘以 $x-x_0$，可知 $B\mid x-x_0$，所以 $x-x_0=Bt$。代回上式就得到 $y-y_0=-At$

例如，求解

$$252x+198y=36$$

将 $18=4\times252-5\times198$ 乘以 $2$，得到特解 $(8,-10)$。全部解为

$$x=8+11t,\qquad y=-10-14t,\qquad t\in\mathbb{Z}$$

应用：模逆元

设 $n>1$。若

$$ax\equiv1\pmod n,$$

也就是 $n\mid ax-1$，则称 $x$ 是 $a$ 模 $n$ 的乘法逆元

这等价于存在整数 $y$，使得

$$ax+ny=1$$

因此逆元存在的充要条件是 $\gcd(a,n)=1$，扩展欧几里得算法可以直接给出逆元

例如

$$1=2\times26-3\times17,$$

所以 $17$ 模 $26$ 的逆元为 $-3\equiv23\pmod{26}$。检查可得 $17\times23=391=15\times26+1$

在 RSA 的一种常见构造中，给定不同素数 $p,q$，记 $\varphi(n)=(p-1)(q-1)$。选取与 $\varphi(n)$ 互素的公钥指数 $e$ 后，私钥指数 $d$ 满足

$$ed\equiv1\pmod{\varphi(n)}$$

这一步就是用扩展欧几里得算法求模逆元