(n*2-1)%p：n和p为32位时避免使用64位

Posted 2023-02-22

【中文标题】(n*2-1)%p：n和p为32位时避免使用64位

【英文标题】：(n*2-1)%p: avoiding the use of 64 bits when n and p are 32 bits

【发布时间】：2015-04-20 21:36:13

【问题描述】：

考虑以下函数：

inline unsigned int f(unsigned int n, unsigned int p) 

    return (n*2-1)%p;

现在假设n（和p）大于std::numeric_limits<int>::max()。

例如f(4294967295U, 4294967291U)。

数学结果是7，但函数会返回2，因为n*2会溢出。

那么解决方案很简单：我们只需要使用 64 位整数来代替。假设函数的声明必须保持不变：

inline unsigned int f(unsigned int n, unsigned int p) 

    return (static_cast<unsigned long long int>(n)*2-1)%p;

一切都很好。至少在原则上。问题是这个函数将在我的代码中被调用数百万次（我的意思是溢出版本），并且 64 位模数比 32 位版本慢得多（例如，参见 here）。

问题如下：是否有任何技巧（数学或算法）来避免执行 64 位版本的模运算。使用这个技巧的f 的新版本会是什么？ （保持相同的声明）。

注1：n > 0 注2：p > 2 注3：n可以低于p：n=4294967289U，p=4294967291U 注4：使用的模运算次数越少越好（3 32位模太大，2很有趣，1肯定会胜过） 注意 5：当然结果将取决于处理器。假设在最后一台可用的至强超级计算机上使用。

【问题讨论】：

你应该使用一个更好的例子，否则你会得到便宜的答案，利用 n 没有预先减少

为什么不妥协？仅当答案可能溢出时（即如果 n > 0x80000000），才执行 64 位模数。如果你不经常用大n 调用它，那么这会很好用。

@nneonneo：我在帖子里加了个注释，函数经常会被大n调用。

顺便问一下，您愿意接受 x86 汇编中的答案吗？你可以回避这个问题的大部分，因为 32 位模确实需要 64 位除法（受到一些限制，但如果 p > 1 和 n != 0 在这里可以解决）

请注意，这将非常依赖处理器...我使用以下一些方法进行了快速基准测试（即使它们有缺陷，只是为了了解额外的复杂性是否更多成本高于 64 位操作）。在我的旧 Phenom II 955 上，仅以 64 位进行计算甚至比 32 位略好；如果 n*2 会溢出，则切换到 64 位的成本要高出约 20%，执行 (n%p + n%p) 的成本是两倍。 OTOH，在最近的i7上，如果32位需要1，64位需要2，检查高位成本2.2，做（n％p + n％p）成本1.4。

【参考方案1】：

我们知道p 小于max，那么n % p 小于最大值。它们都是无符号的，这意味着n % p 是正数，并且小于p。无符号溢出是明确定义的，所以如果n % p * 2超过p，我们可以计算为n % p - p + n % p，不会溢出，所以加起来是这样的：

unsigned m = n % p;
unsigned r;
if (p - m < m) // m * 2 > p
    r = m - p + m;
else // m * 2 <= p
    r = m * 2;

// subtract 1, account for the fact that r can be 0
if (r == 0) r = p - 1;
else r = r - 1;
return r % p;

注意你可以避免最后一个模数，因为我们知道r不超过p * 2（最多是m * 2，而m不超过p），所以最后一行可以改写为

return r >= p ? r - p : r

这使模数运算的数量为 1。

【讨论】：

"在 64 位平台上，64 位整数不应比 32 位整数慢。"在谈论除法时，这是不正确的。

同意上述两个 cmets，编辑了帖子，因为它与我的回答无关

我想知道如何评估r = m - p + m。 m - p 低于零，所以这不会导致溢出，因为r 是无符号的吗？我对c++不熟悉。

是的，m - p 确实会导致有意的下溢，但无符号整数的下溢是明确定义的。添加m 会溢出回来，因为m + m 在该分支中至少与p 一样大。

不错。我不知道它会溢出来。好主意。

【参考方案2】：

尽管我不喜欢处理 AT&T 语法和 GCC 的“扩展 asm 约束”，但我认为这很有效（它在我的、公认的有限测试中有效）

uint32_t f(uint32_t n, uint32_t p)

    uint32_t res;
    asm (
      "xorl %%edx, %%edx\n\t"
      "addl %%eax, %%eax\n\t"
      "adcl %%edx, %%edx\n\t"
      "subl $1, %%eax\n\t"
      "sbbl $0, %%edx\n\t"
      "divl %1"
      : "=d"(res)
      : "S"(p), "a"(n)
      : 
      );
  return res;

我不知道，约束可能是不必要的严格或错误。它似乎奏效了。

这里的想法是进行常规的 32 位除法，实际上需要 64 位除法。它仅在商适合 32 位时才有效（否则会发出溢出信号），在这种情况下总是如此（p 至少为 2，n 不为零）。除法之前的东西处理时间2（溢出到edx，“高半部分”），然后是带有潜在借位的“减1”。 "=d" 输出事物使其将余数作为结果。 "a"(n) 将 n 放入 eax （让它选择其他寄存器没有帮助，除法将在 edx:eax 中输入）。 "S"(p) 可能是 "r"(p)（似乎有效），但我不敢相信它。

【讨论】：

有没有办法将结果存储在 n 而不是 res 以避免不必要的分配？ + 你能解释一下“S”（d）是什么意思吗？

@Vincent 也许可以，但这并不能真正避免任何事情，没有内存写入或类似的东西，只是一个 movl %edx, %eax 来返回它。 "S"(d) 将 d 放入 esi（它已经存在）编辑：哦等等，它是 p。好吧，无论如何，它会将 p 放入 esi 中。

【参考方案3】：

FWIW，这个版本似乎避免了任何溢出：

std::uint32_t f(std::uint32_t n, std::uint32_t p) 

    auto m = n%p;
    if (m <= p/2) 
        return (m==0)*p+2*m-1;
    
    return p-2*(p-m)-1;

Demo。这个想法是，如果在2*m-1 中发生溢出，我们可以使用p-2*(p-m)-1，通过将2 与模加法逆相乘来避免这种情况。

【讨论】：

如果 p 大于该类型可以表示的最大值的一半，则可能会失败，因为 (n%p)*2 仍然可以溢出。

仅在该示例中。翻转n 和p... 这样n%p == n。

@Peter 现在看来代码是正确的，有没有溢出？

@Columbo 这个想法是正确的。但是有几个小细节。 f(0, ...) 返回 p - 1，而原始函数返回 ((ulong)-1) % p。另外，如果m == p/2，我们得到(uint)-1。

@AlexD 第一部分是有意。 -1 等价于 p-1 模 p。但是，是的，我知道后者有问题，我修复了它。现在怎么样？