快速求模2^16 + 1的乘法

Question

IDEA密码使用模乘运算，模数为2^16 + 1。是否有一种算法可以在不使用通用模运算符（仅截断模2^16）的情况下执行此操作？在IDEA中，零被解释为2^16（这意味着零不是我们乘法的参数，也不能成为结果，因此我们可以节省一个位并将值2^16存储为位模式0000000000000000）。我想知道如何高效地实现它（或者是否可能）而不使用标准模运算符。

- ciechowoj

---

我已经添加了一些澄清。 - ciechowoj

Answer 1

您可以利用这个事实：(N-1) % N == -1。

因此，(65536 * a) % 65537 == -a % 65537。 同时，当将0解释为65536时，-a % 65537 == -a + 1 (mod 65536)。

uint16_t fastmod65537(uint16_t a, uint16_t b)
{
    uint32_t c;
    uint16_t hi, lo;
    if (a == 0)
        return -b + 1;
    if (b == 0)
        return -a + 1;

    c = (uint32_t)a * (uint32_t)b;
    hi = c >> 16;
    lo = c;
    if (lo > hi)
        return lo-hi;
    return lo-hi+1;
}

这里唯一的问题是如果hi == lo，结果将为0。幸运的是，测试套件确认了它实际上不可能发生...

int main()
{
    uint64_t a, b;
    for (a = 1; a <= 65536; a++)
       for (b = 1; b <= 65536; b++)
        { 
            uint64_t c = a*b;
            uint32_t d = (c % 65537) & 65535;
            uint32_t e = m(a & 65535, b & 65535);
            if (d != e)
                printf("a * b % 65537 != m(%d, %d) real=%d m()=%d\n",
                       (uint32_t)a, (uint32_t)b, d, e);
        }
    }

输出：无

- Aki Suihkonen

---

我仍然不理解最后5行代码是如何工作的。 我猜这可能由于两个if语句导致速度较慢，或者编译器在幕后进行了一些优化。 - ciechowoj

将第一个测试组合为 (a*b == 0) 并将最后一个表达式包装为 c = lo-hi; c += c >> 16;，这样可以与 i5 上的 % 方法相媲美。当矢量化时，无分支版本当然会优于 div 方法。 - Aki Suihkonen

能否从基本原理证明 hi != lo？做 c += c >> 16 是正确的吗？我认为当 lo > hi 时，c 的上位比特都是零，当 lo < hi 时，上位比特都是一。因此，c += c >> 16 实际上意味着 c = c + 2^16 - 1，当取模 2^16 时为 c - 1，而在第一个版本中我们添加了一个而不是减去它（return lo-hi+1;）。 - ciechowoj

是的。hi == lo 要求 a * b = 0x10001 * x。0x10001 是一个质数（实际上是我们取模的相同质数65537），而且a和b都比它小。因此，a*b不能形成这样的数字。另一方面，对于M=2^32+1来说并非如此，它是一个复合数，防止将IDEA扩展到32位整数。 - Aki Suihkonen

最终的包装必须等于lo-hi+1，这实际上相当于c -= c >> 16; // mod 65536，而不是c += xxx，就像在其他答案中一样。 - Aki Suihkonen

Answer 2

首先，当a或b为零的情况。在这种情况下，它被解释为具有值2^16，因此基本模算术告诉我们：

result = -a - b + 1;

因为（在IDEA的上下文中）$2^{16}$的乘法逆元仍然是$2^{16}$，它的最低16位都是零。

一般情况比看起来容易得多，既然我们已经解决了“0”特殊情况（$2^{16}+1$是$0x10001$）：

/* This operation can overflow: */
unsigned result = (product & 0xFFFF) - (product >> 16);
/* ..so account for cases of overflow: */
result -= result >> 16;

综合起来：

/* All types must be sufficiently wide unsigned, e.g. uint32_t: */
unsigned long long product = a * b;
if (product == 0) {
    return -a - b + 1;
} else {
    result = (product & 0xFFFF) - (product >> 16);
    result -= result >> 16;
    return result & 0xFFFF;
}

- Michael Foukarakis