logo标签列表

计算 (a*b) mod c，其中 c=2^N +-1，快速实现。

mathintegermoduloprngknuth
10
在32位整数运算中，基本的加法和乘法操作会隐式地对2^32求模，这意味着你的结果将是加法或乘法的最低位。
如果想要使用不同的模数计算结果，可以使用各种语言中的BigInt类。对于a、b、c < 2^32的值，可以使用64位long int计算中间值，并使用内置的%运算符减少到正确的答案。
但是，我听说有一些特殊技巧可以高效地计算a*b mod C，其中C的形式为(2^N)-1或(2^N)+1，这些技巧不使用64位数学或BigInt库，效率比任意模数求值方法更高，并且可以正确计算那些包括中间乘法的32位int会溢出的情况。
不幸的是，尽管听说这样的特殊情况有快速的求值方法，但我并没有找到该方法的描述。 "Knuth的书上有吗？" "维基百科上有吗？" 这是我听到的喃喃自语。
显然，在进行a*b mod 2147483647（因为2147483647是等于2^31 -1的质数）的乘法时，这是一种常见的技术方法。
所以，我向专家们请教。这个聪明的特殊情况乘以模方法是什么，我找不到任何讨论呢？
6个回答
13

我认为诀窍在于以下内容（我将使用十进制进行演示，因为更容易理解，但原则应该相同）

假设您要计算 a*b mod 10000-1，并且

a = 1234 = 12 * 100 + 34
b = 5432 = 54 * 100 + 32

现在 a*b = 12 * 54 * 10000 + 34 * 54 * 100 + 12 * 32 * 100 + 34 * 32

12 * 54 * 10000 =  648 * 10000
34 * 54 * 100   = 1836 * 100
12 * 32 * 100   =  384 * 100
34 * 32         = 1088

由于 x * 10000 ≡ x (mod 10000-1) [1]，因此第一项和最后一项变成了648+1088。第二项和第三项是“诀窍”的关键。请注意：

1836 = 18 * 100 + 36
1836 * 100 ≡ 18 * 10000 + 3600 ≡ 3618 (mod 10000-1).

本质上这是一个循环移位。计算648 + 3618 + 8403 + 1088的结果。同时请注意，在所有情况下，乘法的数字都小于10000（因为a < 100且b < 100），因此如果只能将2位数相乘并相加，则可以计算出结果。

在二进制中，计算方式类似。

从a和b开始，它们都是32位的。假设你想要对它们进行模2^31-1的乘法，但你只有一个16位的乘数（得到32位）。算法大致如下：

 a = 0x12345678
 b = 0xfedbca98
 accumulator = 0
 for (x = 0; x < 32; x += 16)
     for (y = 0; y < 32; y += 16)
         // do the multiplication, 16-bit * 16-bit = 32-bit
         temp = ((a >> x) & 0xFFFF) * ((b >> y) & 0xFFFF)

         // add the bits to the accumulator, shifting over the right amount
         total_bits_shifted = x + y
         for (bits = 0; bits < total_bits_shifted + 32; bits += 31)
             accumulator += (temp >> (bits - total_bits_shifted)) & 0x7FFFFFFF

         // do modulus if it overflows
         if (accumulator > 0x7FFFFFFFF)
             accumulator = (accumulator >> 31) + (accumulator & 0x7FFFFFFF);

现在已经很晚了，因此累加器部分可能无法正常工作。但我认为原则上是正确的。有人可以进行编辑以使其正确。

展开后，这也相当快，这应该是PRNG使用的方式，我猜测。

[1]: x*10000 ≡ x*(9999+1) ≡ 9999*x + x ≡ x (mod 9999)
1
而且我仍然不理解这背后的数学，这就是为什么我在大学放弃了数学专业的原因... - Jonathan Rupp
1
好的，这有点像在除以9（10-1）时得到余数。你只需要把数字加起来。现在，在这种情况下，你不是在使用十进制或二进制，而是“基于”2^N。 - FryGuy
4
假设你可以将a*b计算为

p*2^N+q

。这可能需要64位计算，或者您可以将a和b分成16位部分并在32位上进行计算。
然后，

a*b mod 2^N-1 = p+q mod 2^N-1

，因为

2^N mod 2^N-1 = 1

。
而且

a*b mod 2^N+1 = -p+q mod 2^N+1

，因为

2^N mod 2^N+1 = -1

。
在两种情况下，都没有除以

2^N-1

或

2^N+1

。
2
这篇论文使用浮点算术而不是 N 的属性来使计算更有效。我自己在浮点计算周围会有点紧张，但除此之外我没有深入检查过...它可能足够好用。 - Pontus Gagge
有趣的论文，值得一读！它提供了一种更通用的方法来处理任意模数值。不幸的是，在计算过程中，它将这些值转换为64位双精度浮点数。虽然这在一般情况下可能是非常高效的计算方式，但对于特殊的c=2^N +-1情况，还有更快的方法。无论如何，因为这是一个很棒的链接，我还是点了个赞！ - SPWorley
1

与其在每一步进行模数减法，你可以使用蒙哥马利约简（还有其他描述）来降低模数乘法计算的成本。然而，这仍然没有利用N是正负二的幂的属性。

1

你要找的身份是 x mod N = (x mod 2^q)- c*floor(x/2^q)，其中 N = 2^q + c，c 是任意整数（但通常为 ±1）。

你可能想阅读 Richard Crandall 和 Carl Pomerance 的《Prime Numbers: A Computational Perspective》中的 第9.2.3节：“特殊形式的模数”。除了理论外，它还包含实现上述关系的算法伪代码。

你确定这是正确的吗？例如，wolfram alpha 原模数计算示例 与 你的公式示例 对于 2^16+1 和 68519263247*1564832。 - JohannesB
我找到了你提到的那本书，我认为你在结尾处忘记了一个重要部分：你的计算的第二部分应该添加mod N，如下所示：x mod N = (x mod 2^q)- c*floor(x/2^q) mod N（[例子](https://www.wolframalpha.com/input/?i=%2868519263247*1564832+mod+%282%5E16%29%29+-+1+*+%28%28floor%28%2868519263247*1564832%29+%2F+2%5E16%29%29+mod+%282%5E16%2B1%29%29）。由于这仍然包含`mod N`，我不确定使用此公式的优点是什么（？） - JohannesB
你可以递归地应用这个公式，直到 c*floor(x/2^q) < N。 - JohannesB
0

我在这个非常相关的主题上找到了一个相当广泛的页面，不仅讨论了算法，甚至还涉及了问题和解决方案的具体历史以及人们使用解决方案的方式。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接