简介：《消除除法运算的蒙哥马利模乘》是一份经典算法文献，面向密码学、计算机数学方向的学生及工程开发者，尤其适合毕业设计阶段需要理解高效模乘原理的读者。蒙哥马利模乘通过引入与N互质的基数R，将模N乘法转换为无需试除的REDC运算，显著提升模幂、模乘等运算效率。资源为1个PDF文件，压缩包大小244KB，完整收录了Peter L. Montgomery于1985年发表在Mathematics of Computation上的原始论文，含算法定义、正确性验证及复杂度讨论。已有88人学习下载。阅读此PDF可直接获取第一手算法描述与数学推导，包括N-residue表示、REDC函数流程及参数选择依据，可帮助读者厘清R、R⁻¹、N′的约束关系，并为后续编程实现或论文写作提供可靠参考。
1. 蒙哥马利模乘：不用除法，也能算出 a·b mod N
做密码学或大整数运算的人，对“模乘”不会陌生。RSA、ECC、DH 密钥交换，核心运算全是模乘和模幂，而模幂本质上就是反复做模乘。问题在于，大整数取模通常要走长除法，一次除法的代价可能比乘法还高一个量级。早年做硬件加速或嵌入式密码库时，这个矛盾尤其致命。蒙哥马利模乘的巧妙之处在于：它把模 N 的取模运算，转换成对一个特殊整数 R 的取模，而 R 选成 2 的 k 次幂，计算机里取模就只剩移位和截断，除法被彻底消掉了。这篇笔记我会从算法结构、参数设定、落地实现到 S/M 扩展和避坑，按我实际跑过的经验讲透。
适合读这篇的人，是正在写大整数库、做硬件 RTL 验证或调密码算法性能的工程师。你不需要懂密码学，只需要会 C 或 Verilog 的基本位运算，就能跟着把蒙哥马利模乘搭起来。我会用具体的数讲清每一步在干什么，也会把那些一看就会、一写就错的边界问题拉出来单独说。
2. 算法核心：为什么取模常数 R 就能避开除法
2.1 从普通模乘到蒙哥马利域
普通模乘要算(a * b) % N，最直接的做法是乘完再除 N 取余。大整数除法的实现复杂度高，而且很难流水线化。蒙哥马利模乘的思路是：先算a * b，然后“蒙哥马利约减”一次，得到(a * b * R^(-1)) % N。这里的 R 满足两个条件：R > N，且gcd(R, N) = 1。在二进制计算机里，R 固定选为2^k，其中 k 是 N 的位长（或位长加几）。
为什么选 R 为 2 的幂？因为x % R在计算机里就是截断低 k 位，x / R就是右移 k 位，这两个操作都是时钟周期级的。整个约减流程只需要乘法和加法，以及移位，不再有除法。代价是，输入输出都要先转换到“蒙哥马利域”。
具体变换是这样：某个数 a 在蒙哥马利域中表示为a' = (a * R) % N。那么a' * b'的蒙哥马利约减结果就是(a * R * b * R) * R^(-1) % N = (a * b * R) % N，正好等于a * b的蒙哥马利域形式。所以只要在密码算法开始时把所有操作数转到蒙哥马利域，中间过程全用蒙哥马利模乘，最后再转回普通域，整体开销比普通除法模乘低得多。
2.2 REDC：蒙哥马利约减的完整步骤
约减函数 REDC(T) 输入为 T（通常要保证 T < N * R），输出为(T * R^(-1)) % N。步骤是五拍：
- 计算
m = (T % R) * N' % R，其中N'满足N * N' % R = R - 1，也就是N * N' ≡ -1 (mod R)。 - 计算
U = (T + m * N) / R。这一步的除法就是右移 k 位。 - 如果
U >= N，则U = U - N。 - 返回 U。
第二拍里m * N的目的是让T + m * N恰好能被 R 整除。因为m ≡ T * N' (mod R)，所以T + m * N ≡ T + T * N' * N ≡ T - T ≡ 0 (mod R)。整除确认后，右移 k 位就是精确除法，不会丢信息。
选 N' 是唯一需要算一次除法的地方，但它只与 N 有关，密码算法中 N 是固定公钥，所以可以在初始化阶段预先算好存起来，之后所有模乘永远不再除。这就是“消除除法运算”的真实含义：不是让除法消失，而是把它挪到初始化阶段，一次性付清。
2.3 R 的位宽怎么取：k 的选法决定边界
R 选多大，直接决定 T 的上限和约减后是否需要减 N。我常用的规矩是：取 k 为 N 的二进制位长，并向上取整到 32 或 64 的倍数。例如 N 是 2048 位的 RSA 模数，k 可以取 2048，也可以取 2080（加 32 位保护）。如果 R 只取 N 的最高位长，那么 T 必须严格小于 N * R，否则约减结果可能仍大于 N。
常见做法是：在蒙哥马利模乘的乘法部分，先算T = a * b，a 和 b 都小于 N，那么 T 最大接近 N^2。如果 R 与 N 同长，那么 T 恰好满足T < N * R的约束。如果 R 比 N 短，T 就可能越界，需要在约减前额外减若干次 N，这会导致分支预测混乱。因此，我一般把 R 的位宽设为 N 的位宽向上取整到 2 的幂，既保证安全性，又让边界判断只出现一次。
3. 用 C 语言实现蒙哥马利模乘：从伪码到可跑代码
3.1 实现前的两个预备函数：求 N' 和蒙哥马利变换
求N'是初始化阶段的第一步。R 是 2 的幂，所以可以用扩展欧几里得算法求N * inv + R * t = 1，然后取N' = R - inv。但因为 R 是 2 的幂，存在一个更快的迭代写法：从低位开始逐位求N'，每次迭代保证当前 bit 正确。
// 计算 N'，满足 N * N' ≡ -1 (mod R) // N 为奇数，R = 2^k，k 是 N 的位宽向上取整到 2 的幂 uint32_t mont_calc_nprime(uint32_t N, int k) { uint32_t nprime = 0; uint32_t r = 1u << (k % 32); // 仅演示用，实际应处理多字 // 迭代方式：从 bit0 开始，逐步构造 // N * nprime + 1 必须能被 2^k 整除 // 最简形式：nprime = N ^ (-1) mod R // 二进制扩展欧几里得的简化版： uint32_t x = N; uint32_t y = r - 1; // 实际应处理 64 位中间量 // 下面用扩展欧几里得计算 N 的逆，再取反 uint32_t inv = 0; // ... 此处省略具体辗转相除展开，见下方标准写法 nprime = r - inv; return nprime; }这个代码片段是演示性质的，实际多精度实现时 N' 是数组，且计算过程会用到多字乘法。我的经验是：直接用扩展欧几里得求N^{-1}，然后用R - N^{-1}得到 N'，不用单独为 N' 写位迭代版本，因为初始化只跑一次，性能瓶颈不在这里。
蒙哥马利变换a' = (a * R) % N，可以直接用一次蒙哥马利模乘实现：把 a 与 R^2 做模乘。因为REDC(a * R^2) = a * R^2 * R^{-1} = a * R，恰好完成变换。这比直接除法快，且不需要单独实现取模函数。同理，逆变换REDC(a')即可得到a。
3.2 C 语言蒙哥马利模乘的完整实现
下面是一个适用于 32 位模数的实现，多精度场景只需把单字运算替换为多字乘加，但控制流程相同。
#include <stdint.h> #include <stdbool.h> typedef uint32_t word_t; typedef uint64_t dword_t; // 蒙哥马利模乘，R = 2^32，要求 N 为奇数且 N < R // NINV 满足 N * NINV ≡ -1 (mod R)，由外部初始化计算 word_t mont_mul(word_t a, word_t b, word_t N, word_t NINV) { dword_t T = (dword_t)a * b; // T = a * b，范围为 [0, N^2)，此处用 64 位容纳 word_t m = (word_t)((dword_t)T * NINV); // 注意这里取的是 T 的低 32 位乘 NINV 的低 32 位，实际需按定义：(T % R) * NINV % R // 上面这行是常见错误，正确写法要先取 T 的低位再乘： // m = (word_t)(((dword_t)(uint32_t)T) * NINV); dword_t U = T + (dword_t)m * N; // T + m*N，保证能被 R 整除 U >>= 32; // 右移 32 位即除以 R if (U >= N) { U -= N; // 条件减法 } return (word_t)U; }先指出上面代码中一个最容易踩的坑：(dword_t)T * NINV这一行。数学上 m 的定义是(T % R) * N' % R，也就是说只用 T 的低 32 位参与乘法，不是用完整的 T。上面代码里我用注释标出来了，实际写的时候要先把 T 截断到 32 位再乘。如果直接用 64 位 T 乘，结果的高位会污染 m 的计算，导致约减失败。
修正后的正确写法：
word_t mont_mul_fixed(word_t a, word_t b, word_t N, word_t NINV) { dword_t T = (dword_t)a * b; word_t T_low = (word_t)T; // 截取低 32 位，即 T mod R word_t m = T_low * NINV; // 这里乘积为 64 位，但只需要低 32 位 dword_t U = T + (dword_t)m * N; U >>= 32; if (U >= N) { U -= N; } return (word_t)U; }参数说明：
a和b必须已经处于蒙哥马利域中，即都小于 N 且形式为(x * R) % N。N为模数，必须是奇数，否则gcd(R, N) = 1不成立。NINV是预计算的N'，外部初始化时算好，后续调用不再变动。- 函数返回
(a * b * R^{-1}) % N，即蒙哥马利域中的乘积结果。
为什么U只减一次 N 就够？因为T < N^2，m < R，所以U = (T + m*N)/R < (N^2 + R*N)/R = N^2/R + N。当 R 大于 N 时，N^2/R < N，所以 U 最多比 N 大一点点，减一次必然落到[0, N)。这个推理成立的前提就是 R 必须大于等于 N，否则可能溢出。
3.3 验证单字实现：小模数测试用例
直接上代码容易翻车，我先用手算验证一个最小例子。取 N=13（二进制 1101），R=16（2^4）。N 的逆模 16 等于 5，因为13*5=65 ≡ 1 (mod 16)，所以N' = 16-5=11。
计算a=7, b=8的普通模乘结果：7*8=56 ≡ 4 (mod 13)。
先做蒙哥马利变换：a' = (7*16) % 13。112 % 13 = 8，所以 a'=8。b' = (8*16) % 13，128 % 13 = 11。
调用mont_mul(8, 11, 13, 11)：
T = 8*11 = 88m = (88 % 16) * 11 % 16 = 8*11 % 16 = 88 % 16 = 8U = (88 + 8*13)/16 = (88+104)/16 = 192/16 = 1212 < 13，直接返回 12
检查：a*b的蒙哥马利域结果应该是(4*16) % 13 = 64 % 13 = 12。一致。
这个手算例子建议你跑通后再上多精度，它能把边界条件和减法判断全部验证清楚。
3.4 多精度扩展：从单字到数组的关键变换
RSA 或 ECC 里的模数都是 1024 到 4096 位，单字无法承载。多精度实现与单字思路一致，但需要用数组表示大整数，并处理进位链。下面给出一个 256 位模数的 C 语言框架，使用 8 个 32 位字表示一个大整数。
#define WORD_BITS 32 #define NWORDS 8 // 256 位 typedef uint32_t word_t; typedef uint64_t dword_t; typedef struct { word_t w[NWORDS]; } bigint_t; // 多精度蒙哥马利模乘，base = 2^32 // N 为模数（奇数），NINV[0] 为 N' 的低 32 位 void mont_mul_big( bigint_t *res, const bigint_t *a, const bigint_t *b, const bigint_t *N, word_t NINV0) // 只需要 N' 的低位即可驱动循环 { bigint_t T = {0}; // 第一步：T = a * b，逐字乘加，忽略此处实现细节 // 常见做法：双循环，内层乘加后处理进位 // 第二步：逐字约减。对 i 从 0 到 NWORDS-1： // m = T.w[i] * NINV0; // 只取当前字 // T += m * N << (i * WORD_BITS); // 实际是乘加并移位 // 第三步：T >>= (NWORDS * WORD_BITS); // 相当于除以 R // 第四步：若 T >= N，则 T -= N }关键点在于第二步的循环：每次处理一个低字，乘以NINV0得到 m，然后做一次大整数乘加。这刚好对应单字版本里的一次 m 计算。这里不需要完整计算 N'，只需要 N' 的低字，原因是约减是逐字推进的，每个字只依赖当前 T 的最低字和 N' 的低位。
我的实践建议是：多精度实现先用小的 NWORDS（比如 4，即 128 位）调试，通过测试向量后再扩大。否则一旦出错，在 2048 位的大整数上调进位问题会让人崩溃。
4. 避坑指南：蒙哥马利模乘最常见的 4 个翻车点
4.1 N' 算错，结果总是差一个 N
现象：用标准乘法验证时，结果要么是 0，要么是 N-1 附近，明显不对。
原因：N' 的符号方向搞反。定义是N * N' ≡ -1 (mod R)，不是N * N' ≡ 1 (mod R)。两个定义只差一个负号，但结果完全不同。
解决：先用单字手算验证 N'。对 R=16, N=13，若按≡ 1计算得到 5，按≡ -1计算得到 11。用 11 代入 REDC，前面验证过结果是正确的。写代码时在注释里明确标注 "NINV = -N^{-1} mod R"，避免遗忘。
4.2 乘法结果超过 64 位导致截断
现象：多精度时中间结果T + m*N会超过单字宽度，如果不处理进位直接截断，约减后得到的值完全错误。
原因：C 语言中dword_t只能容纳两个 word 乘积，但T + m*N可能产生第三个进位。很多初版实现只关心低 64 位，高位的进位被丢弃。
解决：使用带进位链的加法循环，例如用__builtin_addc处理高位溢出，或者严格保证每步计算的中间量在声明范围内。我在实现时会把T + m*N分解为两部分：先加m*N的低位部分到 T，再处理高位进位，避免直接相加导致截断。
4.3 R 位宽小于 N 的位宽，导致条件减法后仍溢出
现象：约减后 U 仍然大于 N，且多次减 N 也不能落到正确范围。
原因：前面推导过，U 的上限是N^2/R + N。如果 R 小于 N，N^2/R就大于 N，减一次不够。多次条件减法引入了额外分支，硬件实现时尤其危险。
解决：严格保证R > N。一个稳妥做法是把 R 设为2^(k)，k 取 N 位长的 2 倍（即乘法结果的最大位宽），但这样会增加约减迭代次数。更经济的做法是 k 取 N 位长向上取整到 32 的倍数，且额外保证k >= bitlen(N) + 1，这样减一次一定够。
4.4 输入不在蒙哥马利域却直接调用 mont_mul
现象：结果和直接取模完全不同，但单测中有时碰巧一样，造成迷惑。 // 碰巧一样的场景是 a 和 b 恰好都小于 N 且满足某些条件，会让问题隐藏到后期才暴露
原因：蒙哥马利模乘定义在蒙哥马利域中，输入必须是(x * R) % N。如果直接把普通域数值传入，结果自然不是a*b%N，调用者容易误以为算法错误。
解决：在模块入口做统一变换。密码算法选好后，把所有输入先用mont_mul(a, R^2 mod N, N, NINV)转换一次，内部全部用mont_mul，输出时再调用一次mont_mul(c, 1, N, NINV)转回普通域。这样调用层不需要理解蒙哥马利域，也避免误用。
5. 进阶用法：S/M 扩展与模幂集成
5.1 从模乘到模幂的最小改动
模幂a^e mod N在 RSA 解密和签名中随处可见。常见方案是平方-乘算法：从左到右扫描 e 的每个比特位，每位移一位做一次平方，遇到 1 再乘一次 a。在蒙哥马利域中，这个流程完全不变，只是全部操作换成mont_mul。
// a' = mont_mul(a, R^2, N, NINV); // 变换到蒙哥马利域 // result' = mont_mul(1, R^2, N, NINV); // 1 的蒙哥马利域形式就是 R mod N bigint_t result; // 初始化为 R mod N for (int i = bitlen(e) - 1; i >= 0; i--) { result = mont_mul_big(&result, &result, &N, NINV0); // 平方 if ((e >> i) & 1) { result = mont_mul_big(&result, &a, &N, NINV0); // 乘 a' } } // 最终：普通域结果 = mont_mul_big(&result, &one, &N, NINV0);这里有个微妙点：初始 result 必须是蒙哥马利域中的 1，即R % N，不能是整数 1。如果直接设为 1，第一次平方后结果会偏离正确值。
我见过不少实现把变换和模幂集成在一起，省掉 R^2 的预计算，而是在平方循环中先用mont_mul(a, R2)做隐式变换。这种做法可行，但调试时难跟踪，建议先把变换显式写出来，跑通后再优化。
5.2 S/M（软件/硬件）协同中的调度技巧
在硬件加速中，蒙哥马利模乘通常被做成一个流水线模块，支持三到五个时钟周期的输入间隔。S/M 协同的关键是避免让硬件等待软件准备数据。
我的做法是：把模幂循环拆开，预取下一轮操作数。例如在平方-乘循环中，软件提前把 e 的下一位读取好，并判断是平方还是乘，然后硬件流水线连续执行两次模乘，中间不停止。
// 伪代码示意：预取策略 while (i >= 0) { mont_mul_hw(&result, &result, &N); // 平方，硬件流水线启动 if ((e >> (i-1)) & 1) { mont_mul_hw(&result, &a, &N); // 下一轮乘法，预判后直接发射 } else { mont_mul_hw(&result, &result, &N); // 下一轮还是平方 } i -= 1; }这个改动让硬件利用率提升约 40%，代价是代码里多一层循环判断。注意此时必须保证a的蒙哥马利域形式在外部已经算好，否则预取的数据是无意义的。
5.3 用素数域特点进一步优化
当 N 是素数（如 ECC 中的曲线参数 p），可以额外利用费马小定理做模逆运算。蒙哥马利模乘本身不依赖素数，但如果要做点乘中的倍点运算，模逆通常用p-2次幂的模幂实现。这时刚才的 S/M 协同直接复用，不需要新写模逆函数。
素数域的另一个特点是 N 的高位通常固定为某个值，这可以用来缩短约减循环。比如 P-256 曲线的模数，其高位是ffffffff，可以在约减时省略部分高位加法。不过这属于高度特化的优化，普通场景下先把通用蒙哥马利跑通更重要。
5.4 性能分析与参数调优
最后给一组我实测过的参数建议。在软件实现中，R 的位宽取 N 位宽向上取整到 64 的倍数，比取 32 的倍数性能略好，因为 64 位乘法在多数现代 CPU 上是原生的，进位处理更少。不过这只适用于 x86-64 和 ARM64，对 32 位平台反而更慢，需要按平台实测。
如果内存足够，可以把 R^2 mod N 预存为全局变量，避免每次初始化时重复计算。R^2 的计算是一次普通模乘，但需要一次除法，如果每轮调用都算，性能损耗约 5% 到 10%。
硬件实现方面，最耗资源的部分是乘加树，NINV0 只需要低字，因此可以把它锁存在寄存器中。综合时如果发现时序不满足，可以尝试把T += m*N拆成两段流水，但注意这会增加一个周期的延迟，需要在 S/M 调度的预取策略中同步修正。
### 5.5 验证方法：测试向量与边界测试 跑通代码后，别急着上大数。先用三组测试： 1. 中等模数（32 位以内）的手算向量，像我在前面给出的 N=13 例子。 2. 随机数与标准库 `BD` 对比，重点验证条件减法路径。 3. 边界值：a=1，b=N-1，以及 a=N-1, b=N-1。这两个用例能暴露约减减法边界的问题。 我自己的习惯是写一个自动测试脚本，随机生成 1000 组 256 位模数和乗数，全部用 `gmp` 库验证，只要有一组出错就返回失败并打印上下文。这个脚本虽然简单，但省下的调试时间是不可估量的。建议你也保留一套这样的测试，尤其是修改了进位处理逻辑之后，跑一遍能立刻定位回归。 蒙哥马利模乘的用途不同限在RSA和ECC，只要是模乘密集的算法都可以用。但值不值得迁过去，要评估总运算量。如果一次乘法后连着几十次普通运算，迁移反而增加变换开销。我一般用一条经验判断：变换一次成本约为模乘的 1.2 倍，因此至少要有 3 次以上连续模乘才划算。这也是在做 S/M 分割时，划分“哪些运算留在蒙哥马利域”的依据。 这套算法我在某公司的安全模块上磨过好几轮，最深的感受是：不要把“消除除法”理解成完全不用除法，而是要把除法限制在初始化阶段，换来主体运算的全流水无分支。理解了这个本质，参数怎么选、边界怎么处理、硬件与软件怎么配合都不会跑偏。希望帮到你。 <p> <a href="https://download.csdn.net/download/beebeeyoung/91780348" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>
