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

Question

考虑以下函数：

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：当然结果将取决于处理器。假设在最后一台可用的至强超级计算机上使用。

Answer 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。

Answer 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)（似乎有效），但我不敢相信它。

Answer 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 与模加法逆相乘来避免这种情况。