模乘函数：将两个整数乘以一个模

Posted 2023-02-22

【中文标题】模乘函数：将两个整数乘以一个模

【英文标题】：Modulo Multiplication Function: Multiplying two integers under a modulus

【发布时间】：2021-07-06 03:08:29

【问题描述】：

我在 miller-rabin 素数测试的代码中遇到了这个模乘函数。这应该是为了消除计算( a * b ) % m时发生的整数溢出。

我需要一些帮助来了解这里发生了什么。为什么这行得通？那个数字文字0x8000000000000000ULL的意义是什么？

unsigned long long mul_mod(unsigned long long a, unsigned long long b, unsigned long long m) 
    unsigned long long d = 0, mp2 = m >> 1;
    if (a >= m) a %= m;
    if (b >= m) b %= m;
    for (int i = 0; i < 64; i++)
    
        d = (d > mp2) ? (d << 1) - m : d << 1;
        if (a & 0x8000000000000000ULL)
            d += b;
        if (d >= m) d -= m;
        a <<= 1;
    
    return d;

【问题讨论】：

次要注意：不需要0x8000000000000000ULL 中的LL。

我们将位向左移动，直到高位变为 1 - 然后我们知道我们不能再移动，否则我们将移走我们需要的配合。

经典的例子，作者的几个 cmets 会有很大帮助。

虽然这很有趣且可移植，但如果您的工具链提供bigger integer type，您可以使用may want。

【参考方案1】：

此代码目前出现在 the modular arithmetic Wikipedia page 上，仅适用于最多 63 位的参数 -- 见底部。

概述

计算普通乘法a * b 的一种方法是添加b 的左移副本——a 中的每个 1 位一个。这类似于我们大多数人在学校做长乘法的方式，但简化了：因为我们只需要将b 的每个副本“乘”1 或0，我们需要做的就是添加@ 的移位副本987654328@（a 对应位为 1 时）或什么都不做（为 0 时）。

这段代码做了类似的事情。但是，为了避免溢出（主要是；见下文），它不是移动 b 的每个副本然后将其添加到总数中，而是将 b 的未移动副本添加到总数中，并依赖于稍后执行的左移总将其移到正确的位置。您可以认为这些转变“作用于”迄今为止添加到总数中的所有总和。例如，第一个循环迭代检查a 的最高位，即位63，是否为1（这就是a & 0x8000000000000000ULL 所做的），如果是，则将b 的未移位副本添加到总数中；到循环完成时，前一行代码会将d 的总数向左移动 1 位 63 次以上。

这样做的主要优点是我们总是添加两个我们已经知道小于m的数字（即b和d，因此处理模环绕很便宜：我们知道那b + d < 2 * m，所以为了确保我们到目前为止的总数仍然小于m，检查是否b + d < m就足够了，如果不是，减去m。如果我们改用 shift-then-add 方法，我们将需要一个 % 每位模运算，这与除法一样昂贵——而且通常比减法昂贵得多。

模算术的一个特性是，每当我们想要执行一系列算术运算以某个数 m 为模时，以通常的算术执行它们并在最后取余数模 m 总是产生与对每个中间结果取余数模 m 相同的结果（前提是没有发生溢出）。

代码

在循环体的第一行之前，我们有不变量d < m 和b < m。

线

d = (d > mp2) ? (d << 1) - m : d << 1;

是一种将d 的总数左移一位的谨慎方法，同时将其保持在0 .. m 的范围内并避免溢出。我们不是先移动它然后测试结果是否为m 或更大，而是测试它当前是否严格高于RoundDown(m/2)——因为如果是这样，加倍后，肯定会严格高于2 * RoundDown(m/2) >= m - 1，因此需要减去m 才能回到范围内。请注意，即使(d << 1) - m 中的(d << 1) 可能会溢出并丢失d 的最高位，但这并没有什么坏处，因为它不会影响减法结果的最低64 位，这是我们唯一感兴趣的（另请注意，如果 d == m/2 准确无误，我们会以 d == m 结束，这有点超出范围 - 但是将测试从 d > mp2 更改为 d >= mp2 以解决此问题将打破 @ 987654363@是奇数，d == RoundDown(m/2)，所以我们不得不忍受这个。没关系，因为它会在下面修复。）

为什么不直接写d <<= 1; if (d >= m) d -= m;呢？假设，在无限精度算术中，d << 1 >= m，所以我们应该执行减法——但d 的高位打开，d << 1 的其余部分小于m：在这种情况下，初始移位将丢失高位，if 将无法执行。

限制为 63 位或更少的输入

上述边缘情况只能在d 的高位打开时发生，这种情况只能在m 的高位也打开时发生（因为我们保持不变的d < m）。因此，即使m 的值非常高，看起来代码也很难正常工作。不幸的是，它仍然可以在其他地方溢出，从而导致某些设置最高位的输入的答案不正确。例如，当a = 3、b = 0x7FFFFFFFFFFFFFFFULL 和m = 0xFFFFFFFFFFFFFFFFULL 时，正确答案应该是0x7FFFFFFFFFFFFFFEULL，但the code will return 0x7FFFFFFFFFFFFFFDULL（查看正确答案的简单方法是使用a 和@ 的值重新运行987654381@ 已交换）。具体来说，每当d += b 行溢出并使截断的d 小于m 时，就会发生此行为，从而导致错误地跳过减法。

只要记录了这种行为（就像在***页面上一样），这只是一个限制，而不是错误。

解除限制

如果我们替换这些行

    if (a & 0x8000000000000000ULL)
        d += b;
    if (d >= m) d -= m;

与

    unsigned long long x = -(a >> 63) & b;
    if (d >= m - x) d -= m;
    d += x;

the code will work for all inputs，包括设置了最高位的那些。神秘的第一行只是一种无条件（因此通常更快）的写作方式

    unsigned long long x = (a & 0x8000000000000000ULL) ? b : 0;

d >= m - x 测试在 d 上运行 之前它已被修改 -- 这就像旧的 d >= m 测试，但 b（当 a 的最高位是on) 或 0 (否则) 已从两边减去。这将测试 d 是否为 m 或更大的 x 被添加到它。我们知道 RHS m - x 永远不会下溢，因为最大的 x 可以是 b，我们已经在函数顶部确定了 b < m。

【讨论】：

***指出该算法“适用于不大于 63 位的无符号整数”，在这种情况下，不会设置高位并且不会发生错误。

@interjay：感谢您注意到这一点！我已经更新了我的答案。