c++ - 计算（64 位无符号整数）*（64 位无符号整数）除以 2^64 的商

Question

我想计算“两个 64 位无符号整数相乘”除以 2^64 的商。

有人说可以使用 128 位整数，但在某些编程语言或某些平台（如 Visual Studio C++）中不支持内置 128 位整数。

但是我们不想用除法，因为除法太费时间了。我认为它可以通过加法/减法、乘法和按位运算（如位移）来完成。

Answer 1

将您的数字分成两部分（使用位移位和位掩码）并应用一些代数。

• 第一个数字：A*2^32 + C，其中A和C均小于 2^32。
• 第二个数字：B*2^32 + D，其中B和D均小于 2^32。
• (A*2^32 + C) * (B*2^32 + D) = (A*B)*2^64 + (A*D)*2^32 + (B*C)*2^32 + (C*D)
• 除以 2^64：(A*B) + (A*D)/2^32 + (B*C)/2^32 + (C*D)/2^64

所以答案几乎是(A*B) + (A*D)>>32 + (B*C)>>32，但这可能会导致舍入误差。错误是什么？从实数（浮点数）商中减去这个几乎是答案：

• (A*D)&0xFFFFFFFF/2^32 + (B*C)&0xFFFFFFFF/2^32 + (C*D)/2^64（请将除法视为“实数”或浮点数）。
• = [(A*D)&0xFFFFFFFF + (B*C)&0xFFFFFFFF + (C*D)/2^32] / 2^32（再次，真正的划分）
• = [(A*D)&0xFFFFFFFF + (B*C)&0xFFFFFFFF + (C*D)>>32] >> 32加上小于 1 的东西。

所以你可以得到想要的数字 (A*B) + (A*D)>>32 + (B*C)>>32 + [(A*D)&0xFFFFFFFF + (B*C)&0xFFFFFFFF + (C*D)>>32] >> 32

Answer 2

将您的数字 a 和 b 分解为 32 位部分：

a = a1 * 2**32 + a0
b = b1 * 2**32 + b0
a * b = (a1 * b1) * 2**64 + (a1 * b0 + a0 * b1) * 2**32 + a0 * b0

要得到结果的高64位，a1 * b1部分是显而易见的，其他部分应该分解：首先将a1 * b0的32个高位与a0 * b1的32个高位相加；其次将 (a1 * b0 + a0 * b1) 的 32 个低位添加到 a0 * b0 的 32 个高位，并保留此中间结果的 32 个高位（实际上是 1 个有效位），以考虑从低位溢出位之前丢弃它们。

下面是对结果进行基本检查的代码。

#include <cstdint>
#include <iostream>

using namespace std;

inline uint64_t high32(uint64_t x) {
    return x >> 32;
}

inline uint64_t low32(uint64_t x) {
    return static_cast<uint32_t>(x);
}

uint64_t mul64(uint64_t a, uint64_t b)
{
    uint64_t a1 = high32(a);
    uint64_t a0 = low32(a);
    uint64_t b1 = high32(b);
    uint64_t b0 = low32(b);

    uint64_t a1_b0 = a1 * b0;
    uint64_t a0_b1 = a0 * b1;

    return a1 * b1 + high32(a1_b0) + high32(a0_b1)
         + high32(low32(a1_b0 + a0_b1) + high32(a0 * b0));
}

int main()
{
    cout << mul64(UINT64_MAX, UINT64_MAX) << endl;
    cout << UINT64_MAX - 1 << endl;

    cout << endl;
    cout << mul64(UINT64_MAX - 1, UINT64_MAX) << endl;
    cout << UINT64_MAX - 2 << endl;

    cout << endl;
    cout << mul64(UINT64_MAX - 2, UINT64_MAX) << endl;
    cout << UINT64_MAX - 3 << endl;
}