我在米勒-拉宾素数测试的代码中遇到了这个模乘法函数。这应该消除计算( 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;
}
此代码目前出现在模块化算术维基百科页面上，仅适用于最多 63 位的参数 - 见底部。
概述
计算普通乘法a * b的一种方法是添加b的左移副本 -a中每个 1 位一个副本。这类似于我们大多数人在学校做长乘法的方式，但进行了简化：由于我们只需要将b的每个副本"乘以"1或0，我们需要做的就是添加b的移位副本(当相应的a位为1时)或不执行任何操作(当它为0时)。
此代码执行类似操作。但是，为了避免溢出(主要是;见下文)，它不是移动每个b副本，然后将其添加到总数中，而是将未移位的b副本添加到总数中，并依靠以后对总数执行的左移将其移动到正确的位置。您可以认为这些变化"作用于"到目前为止添加到总数中的所有汇总。例如，第一次循环迭代检查a的最高位(即位 63)是否为 1(这就是a & 0x8000000000000000ULL所做的)，如果是，则将b的未移位副本添加到总数中;到循环完成时，前一行代码将把总d向左移动 63 次。
这样做的主要优点是，我们总是将我们已经知道小于m的两个数字(即b和d)相加，因此处理模包很便宜：我们知道b + d < 2 * m，因此为了确保到目前为止我们的总数保持在m， 检查是否b + d < m就足够了，如果没有，则减去m。如果我们使用移位后加法，我们将需要每比特的%模运算，这与除法一样昂贵 - 通常比减法贵得多。
模算术的一个特性是，每当我们想执行一系列算术运算模化一些数字m时，以通常的算术执行它们并在最后取余数模m总是产生与对每个中间结果取余数模m相同的结果(前提是没有发生溢出)。
法典
在循环体的第一行之前，我们有不变量d < m和b < m。
该行
d = (d > mp2) ? (d << 1) - m : d << 1;
是一种谨慎的方法，将总d向左移动 1 位，同时将其保持在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来解决此问题将打破m奇数和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，但代码将返回0x7FFFFFFFFFFFFFFDULL(查看正确答案的一种简单方法是使用a的值重新运行，并交换b)。具体来说，每当行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;
该代码将适用于所有输入，包括设置了顶部位的输入。神秘的第一行只是一种无条件(因此通常更快)的写作方式
unsigned long long x = (a & 0x8000000000000000ULL) ? b : 0;
测试d >= m - x在修改之前对d进行操作 - 它类似于旧的d >= m测试，但b(当a的顶部位打开时)或0(否则)已从两侧减去。这将测试d添加到x后是m还是更大。我们知道 RHSm - x永远不会下溢，因为最大的x可能是b，我们已经在函数的顶部建立了该b < m。