扩展欧几里得与模逆元

公元前三世纪，欧几里得在《几何原本》第七卷中记录了一种求两个正整数最大公因数的方法。两千三百年后，这个方法的扩展形式每天在全球数十亿次 TLS 握手中被调用，支撑着互联网通信的安全基础。

这不是一个”教科书里的古老算法”的故事。这是一个关于数论如何从纯粹的智力游戏变成工程基石的故事。从 RSA 密钥生成到椭圆曲线运算，从中国剩余定理的加速到 Montgomery 乘法的底层机制——扩展欧几里得算法（Extended Euclidean Algorithm）及其衍生技术，构成了现代公钥密码学的计算核心。

本文将从最基本的辗转相除法出发，逐步推导扩展欧几里得算法、模逆元、中国剩余定理、Montgomery 模乘，直到 RSA 私钥运算的完整实现。每一步都配有完整的 C 代码，每一个数学结论都附带手工验算。

一、辗转相除法与最大公因数

基本原理

给定两个非负整数 \(a\) 和 \(b\)（\(b \neq 0\)），辗转相除法基于以下观察：

\[\gcd(a, b) = \gcd(b, a \bmod b)\]

证明很直观：设 \(d = \gcd(a, b)\)，则 \(d \mid a\) 且 \(d \mid b\)。因为 \(a \bmod b = a - \lfloor a/b \rfloor \cdot b\)，所以 \(d \mid (a \bmod b)\)。反过来，如果 \(d' \mid b\) 且 \(d' \mid (a \bmod b)\)，则 \(d' \mid (a \bmod b + \lfloor a/b \rfloor \cdot b) = a\)，所以 \(d' \mid \gcd(a, b)\)。两个方向的整除关系证明了 \(\gcd(a, b) = \gcd(b, a \bmod b)\)。

反复应用这个恒等式，\(b\) 严格递减直到为零，此时 \(a\) 就是最大公因数。

最简实现

/* 递归版本——简洁但栈深度为 O(log min(a,b)) */
int64_t gcd(int64_t a, int64_t b)
{
    if (b == 0) return a;
    return gcd(b, a % b);
}

/* 迭代版本——生产环境首选 */
int64_t gcd_iter(int64_t a, int64_t b)
{
    while (b != 0) {
        int64_t t = b;
        b = a % b;
        a = t;
    }
    return a;
}

复杂度分析

辗转相除法的迭代次数上界是 \(O(\log \min(a, b))\)。更精确地说，最坏情况发生在输入为相邻 Fibonacci 数时。Lame 定理（1844）给出：如果 \(a > b \geq 1\)，则辗转相除法的步数不超过 \(5 \cdot \lfloor \log_{10} b \rfloor\)。

这意味着对于 64 位整数，最多 93 步；对于 2048 位的 RSA 模数，最多约 3000 步。每步一次除法，在现代 CPU 上整数除法延迟约 20-40 个时钟周期，所以即使是大整数的 GCD，开销也很可控。

二、Bezout 恒等式与扩展欧几里得

Bezout 恒等式

法国数学家 Bezout 证明了：对于任意整数 \(a, b\)（不全为零），存在整数 \(x, y\) 使得：

\[ax + by = \gcd(a, b)\]

这个等式被称为 Bezout 恒等式。注意 \(x, y\) 不唯一——如果 \((x_0, y_0)\) 是一组解，那么 \((x_0 + kb/\gcd(a,b), y_0 - ka/\gcd(a,b))\) 对任意整数 \(k\) 也是解。

扩展欧几里得算法

扩展欧几里得算法在计算 \(\gcd(a, b)\) 的同时，找到满足 Bezout 恒等式的系数 \(x, y\)。核心思想是在辗转相除的回溯过程中逐步构造这两个系数。

设递归调用 \(\text{ext\_gcd}(b, a \bmod b)\) 返回了 \((g, x_1, y_1)\)，满足：

\[b \cdot x_1 + (a \bmod b) \cdot y_1 = g\]

将 \(a \bmod b = a - \lfloor a/b \rfloor \cdot b\) 代入：

\[b \cdot x_1 + (a - \lfloor a/b \rfloor \cdot b) \cdot y_1 = g\] \[a \cdot y_1 + b \cdot (x_1 - \lfloor a/b \rfloor \cdot y_1) = g\]

因此，对于原问题 \(ax + by = g\)，有：

\[x = y_1, \quad y = x_1 - \lfloor a/b \rfloor \cdot y_1\]

递归实现

typedef struct { int64_t g, x, y; } ext_gcd_result;

ext_gcd_result ext_gcd(int64_t a, int64_t b)
{
    if (b == 0)
        return (ext_gcd_result){ a, 1, 0 };

    ext_gcd_result r = ext_gcd(b, a % b);
    return (ext_gcd_result){ r.g, r.y, r.x - (a / b) * r.y };
}

迭代实现

递归版本在概念上更清晰，但迭代版本避免了栈溢出风险，在工程中更常用：

ext_gcd_result ext_gcd_iter(int64_t a, int64_t b)
{
    int64_t old_r = a, r = b;
    int64_t old_s = 1, s = 0;
    int64_t old_t = 0, t = 1;

    while (r != 0) {
        int64_t q = old_r / r;

        int64_t tmp = r;
        r = old_r - q * r;
        old_r = tmp;

        tmp = s;
        s = old_s - q * s;
        old_s = tmp;

        tmp = t;
        t = old_t - q * t;
        old_t = tmp;
    }

    return (ext_gcd_result){ old_r, old_s, old_t };
}

执行追踪

以 \(\text{ext\_gcd}(240, 46)\) 为例，下图展示了完整的前向除法和回溯替换过程：

手工验证：\(240 \times (-9) + 46 \times 47 = -2160 + 2162 = 2 = \gcd(240, 46)\)。

三、模逆元

定义与存在条件

给定整数 \(a\) 和正整数 \(m\)，如果存在整数 \(x\) 使得：

\[ax \equiv 1 \pmod{m}\]

则称 \(x\) 为 \(a\) 模 \(m\) 的逆元，记作 \(a^{-1} \pmod{m}\)。

模逆元存在的充要条件是 \(\gcd(a, m) = 1\)，即 \(a\) 和 \(m\) 互素。

证明：由 Bezout 恒等式，\(\gcd(a, m) = 1\) 当且仅当存在整数 \(x, y\) 使得 \(ax + my = 1\)。对两边取模 \(m\)，得 \(ax \equiv 1 \pmod{m}\)。

用扩展欧几里得求模逆元

这是最直接的方法：调用 \(\text{ext\_gcd}(a, m)\)，如果返回的 \(g = 1\)，则 \(x \bmod m\) 就是所求的逆元。

/* 返回 a 模 m 的逆元；若不存在返回 -1 */
int64_t mod_inverse(int64_t a, int64_t m)
{
    ext_gcd_result r = ext_gcd(a, m);
    if (r.g != 1) return -1;  /* 逆元不存在 */
    return ((r.x % m) + m) % m;  /* 确保结果为正 */
}

注意最后一行的 ((r.x % m) + m) % m——这是 C 语言中处理负余数的标准做法。C 标准规定 % 运算符的结果与被除数同号，所以当 r.x 为负数时，r.x % m 是负数，需要加上 m 再取模。

复杂度

扩展欧几里得法的时间复杂度为 \(O(\log m)\) 次除法运算，空间为 \(O(1)\)（迭代版本）。对于任意模数——包括合数——只要 \(\gcd(a, m) = 1\)，都能正确工作。

四、Fermat 小定理与两种求逆方法的对比

Fermat 小定理法

当模数 \(p\) 是素数时，Fermat 小定理给出另一种求逆的途径：

\[a^{p-1} \equiv 1 \pmod{p} \implies a^{-1} \equiv a^{p-2} \pmod{p}\]

用快速幂计算 \(a^{p-2} \bmod p\) 即可。

/* 快速幂求模逆元——仅当 m 为素数时正确 */
int64_t mod_inverse_fermat(int64_t a, int64_t m)
{
    int64_t result = 1;
    int64_t base = a % m;
    int64_t exp = m - 2;

    while (exp > 0) {
        if (exp & 1)
            result = (__int128)result * base % m;
        base = (__int128)base * base % m;
        exp >>= 1;
    }
    return result;
}

对比

在密码学实现中，选择哪种方法取决于具体场景：

• RSA 密钥生成：求 \(d = e^{-1} \bmod \phi(n)\)，其中 \(\phi(n)\) 不是素数，必须用扩展欧几里得。
• 椭圆曲线点运算：在素域 \(\mathbb{F}_p\) 上求逆，两种方法都可用。如果需要常数时间实现，Fermat 方法更容易做到（固定 \(\log_2 p\) 次迭代），但扩展欧几里得也可以通过条件移动（constant-time conditional swap）实现常数时间。
• 批量求逆：Montgomery 的批量求逆技巧（一次 ext-gcd 加 \(3(n-1)\) 次乘法求 \(n\) 个逆元）使得扩展欧几里得在批量场景下更具优势。

五、中国剩余定理

定理陈述

设 \(m_1, m_2, \ldots, m_k\) 是两两互素的正整数，\(M = m_1 m_2 \cdots m_k\)。对于任意整数 \(a_1, a_2, \ldots, a_k\)，同余方程组：

\[\begin{cases} x \equiv a_1 \pmod{m_1} \\ x \equiv a_2 \pmod{m_2} \\ \vdots \\ x \equiv a_k \pmod{m_k} \end{cases}\]

在模 \(M\) 意义下有唯一解。

构造性证明

令 \(M_i = M / m_i\)，即除了 \(m_i\) 之外所有模数的乘积。因为 \(m_1, \ldots, m_k\) 两两互素，所以 \(\gcd(M_i, m_i) = 1\)。

令 \(y_i = M_i^{-1} \bmod m_i\)（用扩展欧几里得计算），则：

\[x = \sum_{i=1}^{k} a_i M_i y_i \pmod{M}\]

验证：对于任意 \(j\)，当 \(i \neq j\) 时 \(M_i \equiv 0 \pmod{m_j}\)（因为 \(m_j \mid M_i\)），当 \(i = j\) 时 \(M_j y_j \equiv 1 \pmod{m_j}\)。所以 \(x \equiv a_j \pmod{m_j}\)。

实现

/*
 * CRT: 给定 k 个同余方程 x = a[i] (mod m[i])，
 *      其中 m[i] 两两互素，求 x mod (m[0]*m[1]*...*m[k-1])。
 *
 * 注意：此实现假设中间结果不溢出 int64_t。
 * 对于大模数，需要使用多精度整数或 __int128。
 */
int64_t crt(const int64_t a[], const int64_t m[], int k)
{
    int64_t M = 1;
    for (int i = 0; i < k; i++)
        M *= m[i];

    int64_t x = 0;
    for (int i = 0; i < k; i++) {
        int64_t Mi = M / m[i];
        int64_t yi = mod_inverse(Mi, m[i]);
        /* 使用 __int128 防止乘法溢出 */
        x = (x + (__int128)a[i] % M * Mi % M * yi % M) % M;
    }

    return (x % M + M) % M;
}

手工示例

求解：\(x \equiv 2 \pmod{3}\)，\(x \equiv 3 \pmod{5}\)，\(x \equiv 2 \pmod{7}\)。

• \(M = 3 \times 5 \times 7 = 105\)
• \(M_1 = 35\)，\(M_2 = 21\)，\(M_3 = 15\)
• \(y_1 = 35^{-1} \bmod 3\)：\(35 \equiv 2 \pmod{3}\)，\(2 \times 2 = 4 \equiv 1 \pmod{3}\)，所以 \(y_1 = 2\)
• \(y_2 = 21^{-1} \bmod 5\)：\(21 \equiv 1 \pmod{5}\)，所以 \(y_2 = 1\)
• \(y_3 = 15^{-1} \bmod 7\)：\(15 \equiv 1 \pmod{7}\)，所以 \(y_3 = 1\)
• \(x = 2 \times 35 \times 2 + 3 \times 21 \times 1 + 2 \times 15 \times 1 = 140 + 63 + 30 = 233\)
• \(233 \bmod 105 = 23\)

验证：\(23 = 7 \times 3 + 2\)，\(23 = 4 \times 5 + 3\)，\(23 = 3 \times 7 + 2\)。

六、CRT 在 RSA 解密中的加速：Garner 算法

RSA-CRT 的动机

标准 RSA 解密计算 \(c^d \bmod n\)，其中 \(n = pq\) 是两个大素数的乘积。直接计算需要对一个 2048 位的模数做模幂运算。

利用 CRT，可以将其分解为两个独立的、在更小模数上的运算：

\[m_p = c^{d_p} \bmod p, \quad m_q = c^{d_q} \bmod q\]

其中 \(d_p = d \bmod (p-1)\)，\(d_q = d \bmod (q-1)\)。

因为模幂运算的开销大致与模数位长的立方成正比，将一个 2048 位的模幂替换为两个 1024 位的模幂，速度提升约 \(2 \times (1/2)^3 = 4\) 倍。实际测试中，RSA-CRT 通常比标准 RSA 快 3-4 倍。

Garner 算法

标准 CRT 公式涉及大模数 \(M\) 的乘法，中间结果可能非常大。Garner 算法用一种增量式的计算方式避免了这个问题。

对于两个模数 \(p, q\) 的情况（即 RSA），Garner 公式为：

\[h = q^{-1} \cdot (m_p - m_q) \bmod p\] \[m = m_q + h \cdot q\]

其中 \(q^{-1} \bmod p\) 可以在密钥生成时预计算并存储（称为 \(q_{\text{inv}}\)）。

/*
 * RSA-CRT 解密。
 * 预计算参数：dp, dq, qinv (= q^{-1} mod p)
 * 此处用 64 位整数演示逻辑；实际 RSA 使用多精度整数。
 */
typedef struct {
    int64_t p, q, n;
    int64_t dp, dq;
    int64_t qinv;  /* q^{-1} mod p */
} rsa_crt_key;

int64_t mod_pow(int64_t base, int64_t exp, int64_t mod)
{
    int64_t result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1)
            result = (__int128)result * base % mod;
        base = (__int128)base * base % mod;
        exp >>= 1;
    }
    return result;
}

int64_t rsa_crt_decrypt(const rsa_crt_key *key, int64_t c)
{
    int64_t mp = mod_pow(c, key->dp, key->p);
    int64_t mq = mod_pow(c, key->dq, key->q);

    /* Garner 公式 */
    int64_t h = (__int128)(mp - mq + key->p) * key->qinv % key->p;
    int64_t m = mq + h * key->q;

    return m;
}

为什么存储五个参数

PKCS#1 标准要求 RSA 私钥存储五元组 \((p, q, d_p, d_q, q_{\text{inv}})\) 而不仅仅是 \((n, d)\)。原因是：

1. 性能：CRT 加速约 4 倍。
2. 预计算：\(q_{\text{inv}}\) 只需在密钥生成时计算一次。
3. 安全性：可以用 \(p, q\) 做 Bellcore 攻击防护（计算后验证 \(m^e \equiv c \pmod{n}\)）。

七、二进制 GCD：Stein 算法

动机：避免除法

在某些体系结构上（特别是嵌入式系统和早期处理器），整数除法比移位和减法昂贵得多。Stein 算法（1961）利用以下恒等式，仅用移位、减法和比较来计算 GCD：

1. 若 \(a, b\) 均为偶数：\(\gcd(a, b) = 2 \cdot \gcd(a/2, b/2)\)
2. 若 \(a\) 为偶数，\(b\) 为奇数：\(\gcd(a, b) = \gcd(a/2, b)\)
3. 若 \(a, b\) 均为奇数且 \(a \geq b\)：\(\gcd(a, b) = \gcd((a-b)/2, b)\)

实现

/*
 * 二进制 GCD（Stein 算法）。
 * 仅使用移位、减法、比较——无除法。
 */
uint64_t binary_gcd(uint64_t a, uint64_t b)
{
    if (a == 0) return b;
    if (b == 0) return a;

    /* 提取公因子 2 的幂次 */
    int shift = __builtin_ctzll(a | b);
    a >>= __builtin_ctzll(a);

    do {
        b >>= __builtin_ctzll(b);
        if (a > b) {
            uint64_t t = a;
            a = b;
            b = t;
        }
        b -= a;
    } while (b != 0);

    return a << shift;
}

复杂度与适用场景

Stein 算法的迭代次数同样是 \(O(\log \max(a,b))\)，但每步只需要移位和减法。在以下场景中优于经典辗转相除：

• 嵌入式系统：无硬件除法器，除法需要软件模拟。
• 多精度整数：大整数除法（\(O(n^2)\) 或更优）比移位（\(O(n)\)）昂贵得多。
• 常数时间实现：可以将条件分支替换为条件移动（CMOV），更容易做到常数时间。

不过，在现代 x86-64 处理器上，硬件除法已经很快（div r64 大约 20-40 个周期），经典辗转相除法通常更快，因为它每步消除更多位。

扩展二进制 GCD

Stein 算法也可以扩展来计算 Bezout 系数，但实现复杂度显著增加。Bernstein 和 Yang（2019）提出的 “divstep” 算法是目前已知最好的常数时间扩展 GCD 实现，被用在 libsecp256k1（Bitcoin 的椭圆曲线库）中。

八、Montgomery 模乘

问题：模乘中的除法

模乘 \(a \cdot b \bmod m\) 中的”取模”操作本质上是一次除法（或等价地，一次乘法加移位的 Barrett 约减）。对于多精度整数，这个除法开销很大。

Montgomery（1985）提出了一种巧妙的方法：选择一个 \(R = 2^k > m\) 且 \(\gcd(R, m) = 1\)（当 \(m\) 为奇数时自动满足），将数字 \(a\) 表示为”Montgomery 形式” \(\hat{a} = aR \bmod m\)，然后在 Montgomery 域中做乘法。

核心操作：Montgomery 约减（REDC）

给定 \(T < mR\)，Montgomery 约减计算 \(TR^{-1} \bmod m\)，过程中不需要做关于 \(m\) 的除法：

1. 预计算 \(m' = -m^{-1} \bmod R\)（用扩展欧几里得，只需做一次）
2. \(q = (T \bmod R) \cdot m' \bmod R\)（低位乘法，天然被 \(R\) 截断）
3. \(t = (T + q \cdot m) / R\)（此除法是精确的位移）
4. 若 \(t \geq m\)，返回 \(t - m\)；否则返回 \(t\)

关键洞察：步骤 2 保证 \(T + qm \equiv 0 \pmod{R}\)，所以步骤 3 的除法是精确的（无余数），可以用右移实现。

Montgomery 乘法

两个 Montgomery 形式的数相乘：

\[\text{MontMul}(\hat{a}, \hat{b}) = \text{REDC}(\hat{a} \cdot \hat{b}) = abR \cdot R^{-1} \bmod m = abR \bmod m = \widehat{ab}\]

结果仍然在 Montgomery 形式中。进入和离开 Montgomery 域各需一次 REDC：

• 进入：\(\hat{a} = \text{REDC}(a \cdot R^2 \bmod m)\)（\(R^2 \bmod m\) 预计算）
• 离开：\(a = \text{REDC}(\hat{a} \cdot 1) = \hat{a} R^{-1} \bmod m = a\)

实现

/*
 * Montgomery 模乘——64 位版本。
 * 适用于奇数模数 m < 2^63。
 */
typedef struct {
    uint64_t m;       /* 模数（奇数） */
    uint64_t m_inv;   /* -m^{-1} mod 2^64 */
    uint64_t r2;      /* R^2 mod m，其中 R = 2^64 */
} mont_ctx;

/* 计算 -m^{-1} mod 2^64，用 Newton 迭代（比 ext-gcd 快） */
static uint64_t mont_neg_inv(uint64_t m)
{
    /* m 必须为奇数 => m^{-1} mod 2 = 1 */
    uint64_t x = 1;
    for (int i = 0; i < 6; i++)  /* 6 次迭代对 64 位足够 */
        x *= 2 - m * x;
    return -x;  /* 返回 -m^{-1} mod 2^64 */
}

/* 计算 R^2 mod m */
static uint64_t mont_r2(uint64_t m)
{
    /* R mod m */
    uint64_t r = ((unsigned __int128)1 << 64) % m;
    /* R^2 mod m = (R mod m)^2 mod m */
    return (unsigned __int128)r * r % m;
}

void mont_init(mont_ctx *ctx, uint64_t m)
{
    ctx->m = m;
    ctx->m_inv = mont_neg_inv(m);
    ctx->r2 = mont_r2(m);
}

/* REDC: 计算 T * R^{-1} mod m */
uint64_t mont_redc(const mont_ctx *ctx, unsigned __int128 T)
{
    uint64_t t_lo = (uint64_t)T;
    uint64_t q = t_lo * ctx->m_inv;
    unsigned __int128 qm = (unsigned __int128)q * ctx->m;
    uint64_t t = (uint64_t)((T + qm) >> 64);
    return (t >= ctx->m) ? (t - ctx->m) : t;
}

/* Montgomery 乘法：计算 a*b*R^{-1} mod m */
uint64_t mont_mul(const mont_ctx *ctx, uint64_t a, uint64_t b)
{
    return mont_redc(ctx, (unsigned __int128)a * b);
}

/* 进入 Montgomery 域：a -> a*R mod m */
uint64_t mont_enter(const mont_ctx *ctx, uint64_t a)
{
    return mont_mul(ctx, a, ctx->r2);
}

/* 离开 Montgomery 域：aR -> a mod m */
uint64_t mont_leave(const mont_ctx *ctx, uint64_t aR)
{
    return mont_redc(ctx, (unsigned __int128)aR);
}

/* Montgomery 模幂 */
uint64_t mont_pow(const mont_ctx *ctx, uint64_t base, uint64_t exp)
{
    uint64_t result = mont_enter(ctx, 1);
    uint64_t b = mont_enter(ctx, base);

    while (exp > 0) {
        if (exp & 1)
            result = mont_mul(ctx, result, b);
        b = mont_mul(ctx, b, b);
        exp >>= 1;
    }

    return mont_leave(ctx, result);
}

性能分析

Montgomery 乘法将模约减中的除法替换为一次额外的乘法和一次移位。在以下场景中显著更快：

• 多精度整数：乘法可以用 Karatsuba 或 NTT 加速，而除法加速困难。
• 重复模乘：进出 Montgomery 域的开销被摊薄。模幂运算（RSA、Diffie-Hellman）是典型场景。
• 硬件实现：FPGA/ASIC 上乘法器远比除法器高效。

在 OpenSSL 和 GnuPG 中，所有 RSA 运算都使用 Montgomery 乘法。

九、Lehmer 算法：大整数的 GCD 加速

经典 GCD 的瓶颈

对于 \(n\) 位大整数，经典辗转相除法的每步需要一次 \(O(n^2)\) 的除法（或 \(O(n \log n)\) 使用 FFT 除法），总共 \(O(n)\) 步，总复杂度 \(O(n^3)\)（或 \(O(n^2 \log n)\)）。对于 RSA 密钥生成中的 2048 位整数，这可能太慢。

Lehmer 的观察

Lehmer（1938）注意到：辗转相除法的商序列主要由高位决定。具体来说，如果我们只看 \(a, b\) 的前几位，可以预测接下来若干步的商。

算法：

1. 取 \(a, b\) 的前 64 位（或一个机器字长），记为 \(\hat{a}, \hat{b}\)。
2. 在单精度上运行辗转相除法若干步，累积变换矩阵 \(\begin{pmatrix} \alpha & \beta \\ \gamma & \delta \end{pmatrix}\)。
3. 一旦高位近似不再可靠（通过检查中间结果是否可能因低位不同而改变商），停止。
4. 用多精度乘法执行累积的变换：\(a' = \alpha a + \beta b\)，\(b' = \gamma a + \delta b\)。
5. 如果 \(b'\) 不够小，用一步经典除法将其缩小，然后回到步骤 1。

渐进改进

Lehmer 算法将 GCD 的复杂度从 \(O(n^3)\) 降低到 \(O(n^2)\)。原因是：步骤 2 中的单精度运算在 \(O(1)\) 时间内消除了 \(O(\text{word\_size})\) 位，步骤 4 的多精度乘法是 \(O(n)\)，总共需要 \(O(n / \text{word\_size})\) 轮，所以总复杂度为 \(O(n^2 / \text{word\_size})\)。

更进一步，Schonhage（1971）和后续工作者发展出半 GCD 算法（Half-GCD），利用分治法将复杂度降低到 \(O(n \log^2 n \log \log n)\)，但实现极其复杂，通常只在非常大的整数（上万位）时才值得使用。

实践中的选择

GMP（GNU Multiple Precision Arithmetic Library）在不同大小的整数上自动切换算法。

十、完整实现：ext-GCD + CRT + Montgomery

以下是一个自包含的 C 实现，集成了扩展欧几里得、CRT 和 Montgomery 模乘。约 250 行，可以直接编译运行。

/* ext_gcd_demo.c
 * 编译：gcc -O2 -Wall -o ext_gcd_demo ext_gcd_demo.c
 * 运行：./ext_gcd_demo
 */

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <assert.h>

/*======================================================================
 * 第一部分：扩展欧几里得算法
 *====================================================================*/

typedef struct { int64_t g, x, y; } egcd_t;

/* 迭代版扩展 GCD */
egcd_t ext_gcd(int64_t a, int64_t b)
{
    int64_t old_r = a, r = b;
    int64_t old_s = 1, s = 0;
    int64_t old_t = 0, t = 1;

    while (r != 0) {
        int64_t q = old_r / r;
        int64_t tmp;

        tmp = r;    r = old_r - q * r;    old_r = tmp;
        tmp = s;    s = old_s - q * s;    old_s = tmp;
        tmp = t;    t = old_t - q * t;    old_t = tmp;
    }

    return (egcd_t){ old_r, old_s, old_t };
}

/* 模逆元：返回 a^{-1} mod m，不存在时返回 -1 */
int64_t mod_inv(int64_t a, int64_t m)
{
    egcd_t r = ext_gcd(a % m, m);
    if (r.g != 1) return -1;
    return ((r.x % m) + m) % m;
}

/*======================================================================
 * 第二部分：中国剩余定理
 *====================================================================*/

/* CRT：解同余方程组 x = a[i] (mod m[i])，返回 x mod M */
int64_t crt_solve(const int64_t a[], const int64_t m[], int k)
{
    int64_t M = 1;
    for (int i = 0; i < k; i++)
        M *= m[i];

    int64_t x = 0;
    for (int i = 0; i < k; i++) {
        int64_t Mi = M / m[i];
        int64_t yi = mod_inv(Mi, m[i]);
        assert(yi != -1);
        x = (x + (__int128)a[i] * ((__int128)Mi * yi % M)) % M;
    }

    return (x % M + M) % M;
}

/* Garner 算法（两模数版本，用于 RSA-CRT） */
int64_t garner2(int64_t a1, int64_t m1, int64_t a2, int64_t m2)
{
    int64_t m2_inv = mod_inv(m2, m1);
    assert(m2_inv != -1);
    int64_t h = (__int128)(a1 - a2 + m1) * m2_inv % m1;
    if (h < 0) h += m1;
    return a2 + h * m2;
}

/*======================================================================
 * 第三部分：Montgomery 模乘
 *====================================================================*/

typedef struct {
    uint64_t mod;     /* 奇数模数 */
    uint64_t inv;     /* -mod^{-1} mod 2^64 */
    uint64_t r2;      /* R^2 mod mod */
} mont_t;

static uint64_t neg_inv64(uint64_t m)
{
    uint64_t x = 1;
    for (int i = 0; i < 6; i++)
        x *= 2 - m * x;
    return -x;
}

static uint64_t compute_r2(uint64_t m)
{
    uint64_t r = ((unsigned __int128)1 << 64) % m;
    return (unsigned __int128)r * r % m;
}

void mont_setup(mont_t *ctx, uint64_t m)
{
    assert(m & 1);  /* m 必须为奇数 */
    ctx->mod = m;
    ctx->inv = neg_inv64(m);
    ctx->r2 = compute_r2(m);
}

uint64_t redc(const mont_t *ctx, unsigned __int128 T)
{
    uint64_t lo = (uint64_t)T;
    uint64_t q = lo * ctx->inv;
    uint64_t t = (uint64_t)((T + (unsigned __int128)q * ctx->mod) >> 64);
    return (t >= ctx->mod) ? (t - ctx->mod) : t;
}

uint64_t mmul(const mont_t *ctx, uint64_t a, uint64_t b)
{
    return redc(ctx, (unsigned __int128)a * b);
}

uint64_t menter(const mont_t *ctx, uint64_t a)
{
    return mmul(ctx, a, ctx->r2);
}

uint64_t mleave(const mont_t *ctx, uint64_t aR)
{
    return redc(ctx, (unsigned __int128)aR);
}

uint64_t mpow(const mont_t *ctx, uint64_t base, uint64_t exp)
{
    uint64_t r = menter(ctx, 1);
    uint64_t b = menter(ctx, base);

    for (uint64_t e = exp; e > 0; e >>= 1) {
        if (e & 1) r = mmul(ctx, r, b);
        b = mmul(ctx, b, b);
    }

    return mleave(ctx, r);
}

/*======================================================================
 * 第四部分：RSA-CRT 演示
 *====================================================================*/

typedef struct {
    int64_t p, q, n, e;
    int64_t dp, dq, qinv;
} rsa_key_t;

int64_t modpow64(int64_t base, int64_t exp, int64_t mod)
{
    int64_t result = 1;
    base %= mod;
    if (base < 0) base += mod;
    while (exp > 0) {
        if (exp & 1)
            result = (__int128)result * base % mod;
        base = (__int128)base * base % mod;
        exp >>= 1;
    }
    return result;
}

void rsa_keygen(rsa_key_t *key, int64_t p, int64_t q, int64_t e)
{
    key->p = p;
    key->q = q;
    key->n = p * q;
    key->e = e;
    key->dp = mod_inv(e, p - 1);
    key->dq = mod_inv(e, q - 1);
    key->qinv = mod_inv(q, p);
}

int64_t rsa_encrypt(const rsa_key_t *key, int64_t m)
{
    return modpow64(m, key->e, key->n);
}

int64_t rsa_decrypt_crt(const rsa_key_t *key, int64_t c)
{
    int64_t mp = modpow64(c, key->dp, key->p);
    int64_t mq = modpow64(c, key->dq, key->q);
    int64_t h = (__int128)(mp - mq + key->p) * key->qinv % key->p;
    if (h < 0) h += key->p;
    return mq + h * key->q;
}

/*======================================================================
 * 第五部分：测试
 *====================================================================*/

static void test_ext_gcd(void)
{
    printf("=== ext_gcd tests ===\n");

    egcd_t r = ext_gcd(240, 46);
    printf("ext_gcd(240, 46) = {g=%ld, x=%ld, y=%ld}\n", r.g, r.x, r.y);
    assert(r.g == 2);
    assert(240 * r.x + 46 * r.y == 2);

    r = ext_gcd(35, 15);
    printf("ext_gcd(35, 15) = {g=%ld, x=%ld, y=%ld}\n", r.g, r.x, r.y);
    assert(r.g == 5);
    assert(35 * r.x + 15 * r.y == 5);

    printf("  PASSED\n\n");
}

static void test_mod_inv(void)
{
    printf("=== mod_inv tests ===\n");

    int64_t inv = mod_inv(3, 11);
    printf("3^{-1} mod 11 = %ld\n", inv);
    assert(inv == 4);  /* 3*4 = 12 = 1 mod 11 */

    inv = mod_inv(17, 3120);
    printf("17^{-1} mod 3120 = %ld\n", inv);
    assert((__int128)17 * inv % 3120 == 1);

    assert(mod_inv(6, 9) == -1);  /* gcd(6,9)=3, 逆元不存在 */
    printf("  PASSED\n\n");
}

static void test_crt(void)
{
    printf("=== CRT tests ===\n");

    int64_t a[] = {2, 3, 2};
    int64_t m[] = {3, 5, 7};
    int64_t x = crt_solve(a, m, 3);
    printf("CRT({2,3,2}, {3,5,7}) = %ld\n", x);
    assert(x == 23);

    printf("  PASSED\n\n");
}

static void test_montgomery(void)
{
    printf("=== Montgomery tests ===\n");

    mont_t ctx;
    mont_setup(&ctx, 97);

    /* 测试模幂：2^10 mod 97 = 1024 mod 97 = 54 */
    uint64_t r = mpow(&ctx, 2, 10);
    printf("2^10 mod 97 = %lu\n", r);
    assert(r == 54);

    /* 测试模幂：3^50 mod 97 */
    r = mpow(&ctx, 3, 50);
    uint64_t expected = modpow64(3, 50, 97);
    printf("3^50 mod 97 = %lu (expected %lu)\n", r, expected);
    assert(r == expected);

    printf("  PASSED\n\n");
}

static void test_rsa_crt(void)
{
    printf("=== RSA-CRT tests ===\n");

    rsa_key_t key;
    rsa_keygen(&key, 61, 53, 17);
    printf("RSA key: n=%ld, e=%ld, dp=%ld, dq=%ld, qinv=%ld\n",
           key.n, key.e, key.dp, key.dq, key.qinv);

    int64_t plaintext = 65;
    int64_t cipher = rsa_encrypt(&key, plaintext);
    int64_t decrypted = rsa_decrypt_crt(&key, cipher);
    printf("encrypt(%ld) = %ld, decrypt_crt = %ld\n",
           plaintext, cipher, decrypted);
    assert(decrypted == plaintext);

    /* 测试多个明文 */
    for (int64_t m = 2; m < 100; m++) {
        int64_t c = rsa_encrypt(&key, m);
        int64_t d = rsa_decrypt_crt(&key, c);
        assert(d == m);
    }
    printf("  All 98 encrypt/decrypt roundtrips PASSED\n\n");
}

int main(void)
{
    test_ext_gcd();
    test_mod_inv();
    test_crt();
    test_montgomery();
    test_rsa_crt();

    printf("All tests passed.\n");
    return 0;
}

编译并运行：

gcc -O2 -Wall -o ext_gcd_demo ext_gcd_demo.c && ./ext_gcd_demo

十一、工程实践：从教科书到生产环境

RSA 私钥操作中的完整流程

一次 TLS 握手中的 RSA 签名操作大致经过以下步骤：

1. 接收：收到待签名的消息哈希 \(h\)。
2. PKCS#1 v1.5 填充或 PSS 填充：将 \(h\) 编码为与 \(n\) 等长的整数 \(m\)。
3. CRT 分解：\(m_p = m \bmod p\)，\(m_q = m \bmod q\)。
4. Montgomery 域转换：将 \(m_p, m_q\) 分别转入 Montgomery 域。
5. 模幂：用 Montgomery 乘法计算 \(s_p = m_p^{d_p} \bmod p\)，\(s_q = m_q^{d_q} \bmod q\)。
6. 离开 Montgomery 域。
7. Garner 合并：用 CRT 重组 \(s = s_q + h \cdot q\)。
8. 验证：计算 \(s^e \bmod n\) 并检查是否等于 \(m\)。如果不等，说明发生了故障攻击（Bellcore attack），丢弃结果。
9. 返回签名 \(s\)。

步骤 8 至关重要。1997 年 Boneh、DeMillo 和 Lipton 证明：如果攻击者能在 CRT 的一个分支中注入故障（例如通过电压毛刺），使 \(s_p\) 出错但 \(s_q\) 正确，那么 \(\gcd(s^e - m, n)\) 可以直接分解 \(n\)。一次故障就足以泄露私钥。

常数时间的重要性

密码学实现的一个核心要求是”常数时间”——执行时间不依赖于秘密数据。以下是常见陷阱：

在扩展欧几里得算法中，迭代次数取决于输入值，这本身就是一个时序侧信道。解决方案：

• 固定步数：总是执行最大可能步数（对于 \(n\) 位整数，约 \(2n\) 步），多余的步用空操作填充。
• Bernstein-Yang divstep：设计为天然常数时间的 GCD 算法，不需要额外填充。
• Fermat 法求逆：固定 \(\log_2 p\) 次平方-乘操作，更容易实现常数时间（但需要素数模数）。

TLS 握手中的数论

一次 TLS 1.3 握手中涉及的数论运算包括：

1. ECDHE 密钥交换：椭圆曲线点乘，底层需要模逆元（用于仿射坐标转换或 Jacobian 坐标的最终归一化）。
2. ECDSA 签名验证（如果使用）：需要模逆元计算 \(s^{-1}\)。
3. RSA-PSS 签名验证（如果使用）：需要模幂运算（Montgomery 乘法 + 平方乘法链）。
4. 证书链验证：每个证书的签名都需要上述运算。

一次典型的握手可能涉及 3-5 个签名验证，每个都在底层调用了扩展欧几里得或其变体。

十二、个人思考与总结

算法选择的决策树

经过这么多年的实践，我的个人决策流程大致如下：

1. 需要求 GCD 吗？
   • 64 位以内：经典辗转相除，一个 while 循环搞定。
   • 多精度（几百到几千位）：Lehmer GCD。
   • 超大整数（万位以上）：Half-GCD，但大概率不需要自己实现——用 GMP。
2. 需要模逆元吗？
   • 模数是素数且不需要常数时间：Fermat 小定理 + 快速幂。
   • 模数是合数或需要最快速度：扩展欧几里得。
   • 密码学场景需要常数时间：Bernstein-Yang divstep 或 Fermat + 常数时间模幂。
3. 需要大量模乘吗？
   • 单次模乘：直接计算。
   • 模幂运算（RSA、DH）：Montgomery 乘法。
   • 批量点乘（椭圆曲线）：Montgomery 乘法 + 批量求逆。

容易犯的错误

我在实际项目中见过（或亲自犯过）的错误：

负数取模。C 语言中 -7 % 3 == -1，不是 2。求模逆元后忘记处理负数是最常见的 bug。解决方案：始终用 ((x % m) + m) % m。

溢出。两个 64 位整数相乘会溢出。在 C 中必须用 __int128 或多精度库。我见过一个生产环境的 RSA 实现因为少了一个类型转换而默默产生错误结果，持续了两年才被发现。

非互素输入。调用 mod_inverse(a, m) 前必须检查 \(\gcd(a, m) = 1\)。不检查的后果是返回一个看起来合理但完全错误的值。

CRT 的模数必须两两互素。如果 \(m_1\) 和 \(m_2\) 不互素，CRT 不适用。不检查这个前提就使用 CRT，结果是未定义行为。

Montgomery 乘法只适用于奇数模数。因为需要 \(\gcd(R, m) = 1\)，当 \(R = 2^k\) 时要求 \(m\) 为奇数。RSA 的模数 \(n = pq\)（两个大素数之积）天然是奇数，所以这在 RSA 中不是问题，但在其他场景需要注意。

两千年的连续性

回到开头。欧几里得不可能预见到他的算法会用在什么地方。但数学的美妙之处在于：一个在公元前三世纪为了理解整除性而发明的算法，经过 Bezout 的扩展、Garner 的优化、Montgomery 的变换、Bernstein 和 Yang 的常数时间改造，至今仍然是全球信息安全基础设施的核心组件。

这不是巧合。GCD 和模逆元是整数运算中最基本的操作之一。任何涉及”在有限域上做算术”的应用——密码学、纠错编码、计算机代数——最终都会回到这里。

如果你只记住这篇文章的一句话，我希望是这个：扩展欧几里得算法不是数论课上的习题，它是你每天使用的互联网的一块基石。

上一篇: 素性测试 下一篇: 格基规约 相关阅读: - 素性测试 - 有限域算术