Montgomery 模乘

背景

在密码学中，经常会涉及到模乘操作：\(a * b \bmod N\)。朴素的实现方法是，先求出 \(a * b\)，再对 N 进行除法，那么余数就是模乘的结果。

但由于此时的 \(a\) \(b\) \(N\) 三个数都很大，在计算机上需要用大整数来表示，而大整数的乘法和除法都是需要耗比较多的时间的。如果用 Schönhage–Strassen 算法，计算两个 \(n\) 位大整数的乘法需要的时间是 \(O(n \log(n) \log(\log(n)))\)。

定义

Montgomery 模乘是一种提高模乘的性能的方法。具体地，Montgomery 模乘需要一个参数 \(R\) 满足 \(R\) 和 \(N\) 互质，且 \(R > N\)，那么 Montgomery 模乘实现的是如下计算：

\[ M(a, b) = a * b * R^{-1} \bmod N \]

用途

这看起来很奇怪，为什么要在原来的 \(a*b\) 上再多乘一项：这样可以把原来对 \(N\) 的除法，转化为对 \(R\) 的除法，而如果把 \(R\) 选择为二的幂次，硬件上实现 \(R\) 的除法就很简单，只需要位运算即可。这样就免去了原来很慢的除法操作。

但是，这种计算方法会引入额外的 \(R^{-1}\) 系数，导致计算结果不是想要的 \(a * b\)。解决方法是，预先计算好 \(R^2 \bmod N\)，然后先用 Montgomery 模乘把要计算的数都乘以 \(R\)：

\[\begin{align} M(a, R^2) \equiv a * R^2 * R^{-1} \equiv aR \pmod N \\ M(b, R^2) \equiv b * R^2 * R^{-1} \equiv bR \pmod N \\ \end{align}\]

此时所有数都自带一个 \(R\) 的系数，此时再进行 Montgomery 乘，会发现：

\[ M(aR, bR) \equiv (aR) * (bR) * R^{-1} \equiv abR \pmod N \]

得到的结果是 \(ab\) 乘以 \(R\) 的形式，那么这个结果可以继续和其他带有 \(R\) 系数的值进行运算。当最后运算完成以后，要把结果恢复到原来的值的时候，再和 \(1\) 进行 Montgomery 乘，即可得到最终结果

\[ M(abR, 1) \equiv abR * 1 * R^{-1} \equiv ab \bmod N \]

因此，如果要用 Montgomery 模乘来进行加速，需要经过三个步骤：

1. 第一步，把所有数都进行 Montgomery 模乘，乘以 \(R^2 \bmod N\)，把所有数都添加上 \(R\) 的系数
2. 第二步，按照正常的流程计算，所有值都带有 \(R\) 系数，同时所有的模乘都被替换为了更加高效的 Montgomery 模乘
3. 第三步，把结果和 \(1\) 进行 Montgomery 模乘，还原为真实值

因此，只要中间计算的过程是大头，初始化和最后的处理时间就可以忽略，可以享受到 Montgomery 模乘带来的性能提升。

算法

接下来介绍 Montgomery 的具体算法，看看它如何提高模乘性能。

REDC

首先介绍 Montogomery 的 REDC 算法，它的步骤是：

1. 预先计算 \(N'\)，满足 \(NN' \equiv -1 \bmod R\)
2. 计算 \(T=a*b\)
3. 计算 \(m=((T \bmod R)N') \bmod R\)
4. 计算 \(t=(T + mN) / R\)
5. 如果 \(t \ge N\)，那么 \(a * b * R^{-1} \bmod N = t - N\)
6. 否则 \(a * b * R^{-1} \bmod N = t\)

可以看到这个过程中只涉及到关于 \(R\) 的除法，不涉及到关于 \(N\) 的除法。下面推导一下为什么上面的公式是正确的：

首先计算过程中出现了 \(t = (T+mN) / R\)，为了保证 \(t\) 是整数，需要证明 \(T+mN\) 可以整除 \(R\)，也就是要求 \(T + mN \equiv 0 \pmod R\)：

\[\begin{align} T + mN &\equiv T + (((T \bmod R)N') \bmod R)N \pmod R \\ &\equiv T + TN'N \pmod R \\ &\equiv T + T(-1) \pmod R \\ &\equiv 0 \pmod R \end{align}\]

说明 \(t\) 是整数。其次，需要证明 \(a * b * R^{-1} \equiv t \pmod N\)：

\[\begin{align} t &= (T + mN) / R \\ &\equiv (T + mN) R^{-1} \pmod N \\ &\equiv TR^{-1} + mNR^{-1} \pmod N \\ &\equiv TR^{-1} + 0 \pmod N \\ &\equiv a * b * R^{-1} \pmod N \end{align}\]

说明 \(t\) 和答案 \(a * b * R^{-1}\) 在模 \(N\) 意义下相等。

根据 \(m\) 的定义，\(m\) 的范围是 \([0, R-1]\)，同时 \(a, b \in [0, N-1]\)，且已知 \(N < R\)，计算 \(t\) 的范围：

\[\begin{align} t &= (T+mN) / R \\ &= (ab + mN) / R \\ &\le ((N-1)(N-1) + (R-1)N) / R \\ &< (RN + RN) / R \\ &= N + N \\ &= 2N \end{align}\]

因此 \(t \in [0, 2N-1]\)。前面已经证明，\(t\) 和答案 \(a * b * R^{-1}\) 在模 \(N\) 意义下相等。在 \([0, 2N-1]\) 范围内，模 \(N\) 意义下相等只有两种可能：相等或者差一个 \(N\)。所以 REDC 算法的最后一步就是：如果 \(t \ge N\)，只可能 \(t\) 和答案差一个 \(N\)，所以答案 \(a * b * R^{-1} \bmod N = t - N\)；否则 \(t < N\)，此时答案 \(a * b * R^{-1} \bmod N = t\)。

大整数 Montgomery 模乘

实际在计算机上运行 Montgomery 算法的时候，由于这些数都很大，因此为了表示大整数，需要用固定位数的整数数组来表示，例如用多个 64 位整数来表示一个大整数。此时，把大整数运算拆成多个 64 位整数的运算，然后把大整数的运算和 Montgomery 模乘结合在一起，得到更高性能的 Montgomery 模乘。

论文 Analyzing and Comparing Montgomery Multiplication Algorithms 分析了几种混合了大整数运算和 Montgomery 模乘的算法。下面讲解论文中提到的部分算法。

在下面的讨论中，假设机器整数的宽度是 \(w\) 位，例如 \(w=64\) 表示用 64 位整数进行运算，此时 \(R=2^{sw}\)，也就是说，\(R\) 等于 \(s\) 个 \(w\) 位整数可以表示的最大值加一，那么除以 \(R\) 相当于舍弃最低的 \(s\) 个 \(w\) 位整数。同时也意味着，\(a\) \(b\) \(N\) 都可以用 \(s\) 个 \(w\) 位整数表示。

Separated Operand Scanning

第一种方法是 Separated Operand Scanning 方法（同时也是 Wikipedia 中提到的 MultiPrecisionREDC 算法），它的步骤是：

第一步：按照传统方式进行大整数乘法，计算出 \(T = a*b\)：

# t = a * b
for i=0 to s-1
    C := 0
    for j=0 to s-1
        (C, S) := t[i+j] + a[j]*b[i] + C
        t[i+j] := S
    t[i+s] := C

得到的结果放在 \(t\) 数组中。

第二步，求 \(t = (T + mN) / R\)，此时 \(T\) 已经计算出来，接下来首先要计算出 \(m=((T \bmod R)N') \bmod R\)，在这里 \(\bmod R\) 就是取大数的最低 \(s\) 个 \(w\) 位整数，因此可以简化大整数乘法为：

# m = ((T % R) * N') % R
for i=0 to s-1
    C := 0
    for j=0 to s-i-1
        (C, S) := m[i+j] + t[i]*n'[j] + C
        m[i+j] := S

接下来求大整数 \(m\) 乘以大整数 \(N\) 的积，求积的同时把结果累加到 \(T\) 上。伪代码：

# t += m * N
for i=0 to s-1
    C := 0
    for j=0 to s-1
        (C, S) := t[i+j] + m[i]*n[j] + C
        t[i+j] := S
    ADD (t[i+s], C)

这里的 ADD 函数指的是大整数加法运算里面，求和后不断进位直到不再进位为止的函数，这里就不展开了。

在这里有一个重要的优化：实际上，不需要把 \(m\) 整个大整数计算出来，而是可以直接求 \(T + mN\)：回忆一下，最初计算 \(T + mN\) 的目的是让结果整除 \(R\)，现在把 \(T + mN\) 的计算拆成 \(s\) 个小步：第 \(i\) 步让结果整除 \(2^{(i+1)w}\)：

第一步：\(T + m_1N \equiv 0 \pmod{2^w}\)，此时 \(m_1 = T * N' \bmod 2^w\)，对应在代码上，就是 m_1 = t[0] * n'[0]，舍去溢出的部分。

第二步：\(T + m_1N + m_22^wN \equiv 0 \pmod{2^{2w}}\)，此时 \(m_2 * 2^w = (T + m_1N) * N' \bmod 2^{2w}\)，此时会惊喜地发现，由于 \(T + m_1N\) 是 \(2^w\) 的倍数，因此计算 \((T + m_1N) * N'\) 的时候，\(T + m_1N\) 的低 \(w\) 位全是 0，也意味着实际上 \((T + m_1N) * N'\) 就是拿 \((T + m_1N) \bmod 2^{2w}\) 的高 \(w\) 位乘以 \(N'\) 的低 \(w\) 位，再左移 \(w\) 位，结果等于 \(m_2 * 2 ^ w\)，所以 \(m_2\) 就等于 \(((T+m_1N) \bmod 2^{2w}) / 2^w * (N' \bmod 2^w)\)，对应在代码上，就是 m_2 = t[1] * n'[0]。

这个过程可以一直继续下去，每一步的 m_i 都可以用 m_i = t[i] * n'[0] 计算。因此不再需要先求 \(m\)，再求 \(T + mN\)，而是可以同时计算：

# t += m * N
for i=0 to s-1
    C := 0
    # W = 2^w
    m := t[i] * n'[0] mod W
    for j=0 to s-1
        (C, S) := t[i+j] + m*n[j] + C
        t[i+j] := S
    ADD (t[i+s], C)

计算出 \(T + mN\) 以后，最后就是除以 \(R\) 了，实际上也非常简单，直接去掉 t 数组的低 \(s\) 项即可：

# u = t / R
for i=0 to s
    u[j] := t[j+s]

最后再用大整数减法和比较，使得结果 \(u\) 落在 \([0, N-1]\) 的范围内，这里就不单独列代码了。把上面的代码合在一起，就得到最终完整 Separated Operand Scanning 算法的伪代码：

# t = a * b
for i=0 to s-1
    C := 0
    for j=0 to s-1
        (C, S) := t[i+j] + a[j]*b[i] + C
        t[i+j] := S
    t[i+s] := C

# t += m * N
for i=0 to s-1
    C := 0
    # W = 2^w
    m := t[i] * n'[0] mod W
    for j=0 to s-1
        (C, S) := t[i+j] + m*n[j] + C
        t[i+j] := S
    ADD (t[i+s], C)

# u = t / R
for i=0 to s
    u[j] := t[j+s]

# return u or u - N
B := 0
for i=0 to s-1
    (B,D) := u[i] - n[i] - B
    t[i] := D
(B,D) := u[s] - B
t[s] := D
if B=0 then
    return t[0], t[1], ... , t[s-1]
else
    return u[0], u[1], ... , u[s-1]

Coarse Integrated Operand Scanning

第二种算法 Coarse Integrated Operand Scanning 是在 Separated Operand Scanning 的基础上，把 \(a*b\) 和后面的计算过程交错进行，放在同一个大循环中，因为后面使用到 t 数组的时候，只会依赖已经计算出来的部分。同时，每次循环结束的时候就把整个 t 数组右移一次，因此原来的 t[i] 就会变成 t[0]，t[i+j] 变成 t[j]。

for i=0 to s-1

    # t = a * b
    C := 0
    for j=0 to s-1
        (C, S) := t[j] + a[j]*b[i] + C
        t[j] := S
    (C, S) := t[s] + C
    t[s] := S
    t[s+1] := C

    # t += m * N
    C := 0
    # W = 2^w
    m := t[0] * n'[0] mod W
    for j=0 to s-1
        (C, S) := t[j] + m*n[j] + C
        t[j] := S
    (C, S) := t[s] + C
    t[s] := S
    t[s+1] := t[s+1] + C

    # t /= W
    for j=0 to s
        t[j] := t[j+1]

# return t or t - N
# save as above, omitted

每次循环都右移一次，那么循环 \(s\) 次就是除以 \(R\)，因此原来的 \(u = t / R\) 一步就不需要了。同时，\(t\) 数组需要的存储空间也缩小了，因为不需要保存完整的 \(a*b\) 的结果。更进一步，还可以把 t += m * N 和移位两步合并在一起进行：

for i=0 to s-1

    # t = a * b
    C := 0
    for j=0 to s-1
        (C, S) := t[j] + a[j]*b[i] + C
        t[j] := S
    (C, S) := t[s] + C
    t[s] := S
    t[s+1] := C

    # t = (t + m * N) / W
    # W = 2^w
    m := t[0] * n'[0] mod W
    (C, S) := t[0] + m*n[0]
    for j=1 to s-1
        (C, S) := t[j] + m*n[j] + C
        t[j-1] := S
    (C, S) := t[s] + C
    t[s-1] := S
    t[s] := t[s+1] + C

# return t or t - N
# save as above, omitted

这样就得到了最终的 Coarse Integrated Operand Scanning 算法，下面是一段用 Rust 语言编写的实现：

// https://github.com/jiegec/rust-monty-comparison/blob/6c941d5c95d37dd9ee8f12aa57df577e0f2b623b/src/lib.rs#L89-L139
let mut res = [0u32; WORDS + 2];
// for i=0 to s-1
for i in 0..WORDS {
    // C := 0
    let mut c = 0;
    // for j = 0 to s-1
    for j in 0..WORDS {
        // (C, S) := t[j] + a[j] * b[i] + C
        let mut cs = res[j] as u64;
        cs += self.num[j] as u64 * other.num[i] as u64;
        cs += c as u64;
        c = (cs >> 32) as u32;
        // t[j] := S
        res[j] = cs as u32;
    }
    // (C, S) := t[s] + C
    let cs = res[WORDS] as u64 + c as u64;
    // t[s] := S
    res[WORDS] = cs as u32;
    // t[s+1] := C
    res[WORDS + 1] = (cs >> 32) as u32;

    // m := t[0]*n'[0] mod W
    let m: u32 = (res[0] as u64 * N_INV as u64) as u32;
    // (C, S) := t[0] + m*n[0]
    let mut cs = res[0] as u64 + m as u64 * N[0] as u64;
    c = (cs >> 32) as u32;
    // for j=1 to s-1
    for j in 1..WORDS {
        // (C, S) := t[j] + m*n[j] + C
        cs = res[j] as u64;
        cs += m as u64 * N[j] as u64;
        cs += c as u64;
        c = (cs >> 32) as u32;
        // t[j-1] := S
        res[j - 1] = cs as u32;
    }
    // (C, S) := t[s] + C
    cs = res[WORDS] as u64 + c as u64;
    // t[s-1] := S
    res[WORDS - 1] = cs as u32;
    // t[s] := t[s+1] + C
    res[WORDS] = res[WORDS + 1] + (cs >> 32) as u32;
}

let res_scalar = MontyBigNum::from_u32_slice(&res[0..WORDS]);
let mut res_scalar_sub = res_scalar;
let borrow = bignum_sub(&mut res_scalar_sub, &MODULO);
if res[WORDS] != 0 || borrow == 0 {
    res_scalar_sub
} else {
    res_scalar
}

OpenSSL 也在函数 bn_mul_mont 中实现了这个算法。

常数时间

在 Montgomery 模乘的最后一步，需要把计算结果和 \(N\) 比较，然后进行减法，这一步会出现一个条件分支，可能会导致运行时间和数据相关，成为一个潜在的测信道攻击的点。因此为了解决这个问题，可以有如下的解决方法（参考论文 Montgomery Arithmetic from a Software Perspective）：

既然条件分支是为了和 \(N\) 比较，那如果去掉这个限制，也就是说让结果在 \([0, 2N-1]\) 的范围而不是 \([0, N-1]\)，看看能否继续把结果传给下一次的 Montgomery 模乘。首先要求 \(R > 2N\)，因为输入参数的范围是模 \(2N\)，而不是模 \(N\)；其次，重新考虑 \(t = (T+mN) / R\) 的放缩：

\[\begin{align} t &= (T+mN) / R \\ &= (ab + mN) / R \\ &\le ((2N-1)(2N-1) + (R-1)N) / R \end{align}\]

此时为了让 \(t < 2N\) 成立，需要额外添加条件 \(R > 4N\)，此时：

\[\begin{align} t &= (T+mN) / R \\ &= (ab + mN) / R \\ &\le ((2N-1)(2N-1) + (R-1)N) / R \\ &< ((\frac{1}{2}R)(2N) + RN) / R \\ &= N + N \\ &= 2N \end{align}\]

也就是说，如果规定 \(R > 4N\)，那么最后一步可以不和 \(N\) 比较，把结果保留在 \([0, 2N-1]\) 的范围，继续进行后续的 Montgomery 模乘。但是这样有一个问题，就是在 RSA 场景下，通常 \(N\) 的位数就是二的幂次，例如 \(N\) 是一个 2048 位的大整数，为了满足 \(R > 4N\)，不得不多存一个整数，相应地 Montgomery 模乘计算中的循环次数也要增加。

另一种思路是，既然和 \(N\) 比较比较费时间，而且需要条件分支，那就改成和 \(R\) 比较，让结果保持在 \([0, R-1]\) 的区间内，而不是原来的 \([0, N-1]\) 区间。此时重新考虑 \(t = (T+mN) / R\) 的放缩：

\[\begin{align} t &= (T+mN) / R \\ &= (ab + mN) / R \\ &\le ((R-1)(R-1) + (R-1)N) / R \\ &< (R^2 + RN) / R \\ &= R + N \end{align}\]

也就是说，结果会在 \([0, R+N-1]\) 的范围内，此时最后一步变成和 \(R\) 比较大小，如果比 \(R\) 大，就减去 \(N\)。和 \(R\) 比较大小很简单，直接看最高位是否为 \(0\) 即可。这种方法并没有消除条件分支，但是把条件分支的开销降到了单次整数比较。如果要保证常数时间的话，可以在结果小于 \(R\) 时，减去零。

逆推

下面尝试逆推出 Montgomery 是如何设计 Montgomery 模乘算法的。

首先目标是求 \(t = a * b \bmod N\)，也就是要找到一个 \(m\) 满足 \(t = a * b + mN \in [0, N-1]\)，但是求 \(m\) 的这一步需要计算除法。既然 \(N\) 不方便计算除法，那就尝试换一个除数 \(R\)，把问题转化为求 \(m\) 使得 \(t = (a * b + mN) / R \in [0, N-1]\)，此时 \(t = a * b * R^{-1} \bmod N\)，但是不要紧，只要能高效地计算 \(a * b * R^{-1} \bmod N\)，就可以高效地计算 \(a * b \bmod N\)。

于是下面的问题就来到了怎么找到 \(m\) 使得 \(t = (a * b + mN) / R \in [0, N-1]\)。把 \(R\) 挪到等式左边，得到

\[ t * R = a * b + m * N \]

由于 \(t\) 未知，所以等式左边的值未知，为了消除未知数，等式两边同时对 \(R\) 取余：

\[ 0 \equiv a * b + m * N \pmod R \]

于是 \(m \equiv - a * b * N^{-1} \pmod R\)，其中 \(-N^{-1} \bmod R\) 的部分可以预先计算好。但是这里依然是在模 \(R\) 的意义下相等，\(m\) 有多种可能。

很幸运的是，前面已经推导过，当 \(m = ((a * b) \bmod R * -N^{-1}) \bmod R\) 时，恰好可以保证最终结果 \(t \in [0, 2N-1]\)，此时只需要简单判断一下 \(t\) 和 \(N\) 的大小关系，就可以求出答案。这样就把 Montgomery 模乘的流程逆推出来了。

所以两点观察很重要：一是换一个容易做除法的除数 \(R\)，想到要设 \(t = (a * b + mN) / R\)；二是即使让模乘加上一个系数 \(R^{-1}\)，也不妨碍它的使用。

参考资料

• Wikipedia
• Analyzing and Comparing Montgomery Multiplication Algorithms

Comments