获取 64 位整数乘法的高位部分
Getting the high part of 64 bit integer multiplication
在 C++ 中，表示：
uint64_t i;
uint64_t j;
然后 i * j 将产生一个 uint64_t，其值为 i 和 j 之间的乘积的较低部分，即 (i * j) mod 2^64。
现在，如果我想要乘法的较高部分怎么办？我知道在使用 32 位整数时存在类似的汇编指令，但我对汇编一点都不熟悉，所以我希望得到帮助。
制作如下内容的最有效方法是什么：
uint64_t k = mulhi(i, j);
长乘法性能应该没问题。
将 a*b 拆分为 (hia+loa)*(hib+lob)。这给出了 4 个 32 位乘法加上一些移位。用64位做，手动做进位，你会得到高的部分。
请注意，高部分的近似值可以用更少的乘法来完成——1 次乘法精确到 2^33 左右，3 次乘法精确到 1 以内。
我认为没有可移植的替代品。
如果您使用的是 gcc，并且您的版本支持 128 位数字（尝试使用 __uint128_t），那么执行 128 乘法并提取高 64 位可能是最有效的获取方式结果。
如果您的编译器不支持 128 位数字，那么 Yakk 的回答是正确的。但是，对于一般消费而言，它可能过于简短。特别是，实际的实现必须小心溢出 64 位整数。
他提出的简单且可移植的解决方案是将 a 和 b 中的每一个分解为 2 个 32 位数字，然后使用 64 位乘法运算将这些 32 位数字相乘。如果我们写：
uint64_t a_lo = (uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint32_t)b;
uint64_t b_hi = b >> 32;
那么很明显：
a = (a_hi << 32) + a_lo;
b = (b_hi << 32) + b_lo;
和：
a * b = ((a_hi << 32) + a_lo) * ((b_hi << 32) + b_lo)
= ((a_hi * b_hi) << 64) +
((a_hi * b_lo) << 32) +
((b_hi * a_lo) << 32) +
a_lo * b_lo
前提是使用 128 位（或更高）算法执行计算。
但是这道题需要我们用64位算法进行所有的计算，所以我们不得不担心溢出。
由于 a_hi、a_lo、b_hi 和 b_lo 都是无符号的 32 位数字，它们的乘积将适合无符号的 64 位数字而不会溢出。但是上面计算的中间结果不会。
当数学必须以 2^64 为模执行时，以下代码将实现 mulhi(a, b)：
uint64_t a_lo = (uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint32_t)b;
uint64_t b_hi = b >> 32;
uint64_t a_x_b_hi = a_hi * b_hi;
uint64_t a_x_b_mid = a_hi * b_lo;
uint64_t b_x_a_mid = b_hi * a_lo;
uint64_t a_x_b_lo = a_lo * b_lo;
uint64_t carry_bit = ((uint64_t)(uint32_t)a_x_b_mid +
(uint64_t)(uint32_t)b_x_a_mid +
(a_x_b_lo >> 32) ) >> 32;
uint64_t multhi = a_x_b_hi +
(a_x_b_mid >> 32) + (b_x_a_mid >> 32) +
carry_bit;
return multhi;
正如 Yakk 指出的那样，如果您不介意在高 64 位中被 +1 偏移，则可以省略进位位的计算。
TL:DR 与 GCC 用于 64 位 ISA：(a * (unsigned __int128)b) >> 64 编译得很好，可以编译为单个全乘或高半乘指令。 不需要乱用内联汇编。
不幸的是当前的编译器不优化@craigster0 的便携版本，所以如果你想利用 64 位 CPU ，你不能使用它，除非作为你没有 #ifdef 的目标的后备。 （我没有看到优化它的通用方法；您需要 128 位类型或内部类型。）
大多数 64 位平台上的 GNU C（gcc、clang 或 ICC）has unsigned __int128。 （或者在旧版本中，__uint128_t）。不过，GCC 并未在 32 位平台上实现此类型。
这是让编译器发出 64 位全乘指令并保留高半部分的简单而有效的方法。 （GCC 知道 uint64_t 转换为 128 位整数的上半部分仍然全为零，因此您不会使用三个 64 位乘法得到 128 位乘法。）
MSVC also has a __umulh intrinsic 用于 64 位高半乘法，但同样它仅适用于 64 位平台（特别是 x86-64 和 AArch64。文档还提到 IPF (IA-64) 具有 _umul128 可用，但我没有可用的 Itanium MSVC。（可能无论如何都不相关。）
#define HAVE_FAST_mul64 1
#ifdef __SIZEOF_INT128__ // GNU C
static inline
uint64_t mulhi64(uint64_t a, uint64_t b) {
unsigned __int128 prod = a * (unsigned __int128)b;
return prod >> 64;
}
#elif defined(_M_X64) || defined(_M_ARM64) // MSVC
// MSVC for x86-64 or AArch64
// possibly also || defined(_M_IA64) || defined(_WIN64)
// but the docs only guarantee x86-64! Don't use *just* _WIN64; it doesn't include AArch64 Android / Linux
// https://docs.microsoft.com/en-gb/cpp/intrinsics/umulh
#include <intrin.h>
#define mulhi64 __umulh
#elif defined(_M_IA64) // || defined(_M_ARM) // MSVC again
// https://docs.microsoft.com/en-gb/cpp/intrinsics/umul128
// incorrectly say that _umul128 is available for ARM
// which would be weird because there's no single insn on AArch32
#include <intrin.h>
static inline
uint64_t mulhi64(uint64_t a, uint64_t b) {
unsigned __int64 HighProduct;
(void)_umul128(a, b, &HighProduct);
return HighProduct;
}
#else
# undef HAVE_FAST_mul64
uint64_t mulhi64(uint64_t a, uint64_t b); // non-inline prototype
// or you might want to define @craigster0's version here so it can inline.
#endif
对于 x86-64、AArch64 和 PowerPC64（以及其他），这会编译成一个 mul 指令 ，以及一对 mov 到处理调用约定（在内联之后应该优化掉）。
来自 the Godbolt compiler explorer（使用 x86-64、PowerPC64 和 AArch64 的源代码 + asm）：
# x86-64 gcc7.3. clang and ICC are the same. (x86-64 System V calling convention)
# MSVC makes basically the same function, but with different regs for x64 __fastcall
mov rax, rsi
mul rdi # RDX:RAX = RAX * RDI
mov rax, rdx
ret
（或使用 clang -march=haswell 启用 BMI2：mov rdx, rsi / mulx rax, rcx, rdi 将高半部分直接放入 RAX。gcc 很笨，仍然使用额外的 mov.)
对于 AArch64（使用 gcc unsigned __int128 或使用 __umulh 的 MSVC）：
test_var:
umulh x0, x0, x1
ret
使用编译时常量 2 的乘数，我们通常会得到预期的右移以获取几个高位。但是 gcc 有趣地使用 shld（参见 Godbolt link）。
不幸的是，当前的编译器不优化@craigster0 的便携版本。你得到 8x shr r64,32、4x imul r64,r64 和一堆针对 x86-64 的 add/mov 指令。即它编译成很多 32x32 => 64 位乘法和解包结果。所以如果你想要一些利用 64 位 CPU 的东西，你需要一些 #ifdefs.
一个全乘 mul 64 指令在 Intel CPU 上是 2 微指令，但仍然只有 3 个周期延迟，与 imul r64,r64 相同，它只产生 64 位结果。因此，__int128 / intrinsic 版本在现代 x86-64 上的延迟和吞吐量（对周围代码的影响）比便携式版本便宜 5 到 10 倍，这是基于 http://agner.org/optimize/ 的快速眼球猜测。
在上面 link.
上的 Godbolt 编译器资源管理器中查看
gcc 在乘以 16 时确实完全优化了这个函数，但是：你得到一个右移，比 unsigned __int128 乘法更有效。
这是我今晚想出的单元测试版本，提供完整的 128 位产品。经过检查，它似乎比大多数其他在线解决方案（例如 Botan 库和此处的其他答案）更简单，因为它利用了代码注释中解释的中间部分不会溢出的方式。
为了上下文，我为这个 github 项目编写了它：https://github.com/catid/fp61
//------------------------------------------------------------------------------
// Portability Macros
// Compiler-specific force inline keyword
#ifdef _MSC_VER
# define FP61_FORCE_INLINE inline __forceinline
#else
# define FP61_FORCE_INLINE inline __attribute__((always_inline))
#endif
//------------------------------------------------------------------------------
// Portable 64x64->128 Multiply
// CAT_MUL128: r{hi,lo} = x * y
// Returns low part of product, and high part is set in r_hi
FP61_FORCE_INLINE uint64_t Emulate64x64to128(
uint64_t& r_hi,
const uint64_t x,
const uint64_t y)
{
const uint64_t x0 = (uint32_t)x, x1 = x >> 32;
const uint64_t y0 = (uint32_t)y, y1 = y >> 32;
const uint64_t p11 = x1 * y1, p01 = x0 * y1;
const uint64_t p10 = x1 * y0, p00 = x0 * y0;
/*
This is implementing schoolbook multiplication:
x1 x0
X y1 y0
-------------
00 LOW PART
-------------
00
10 10 MIDDLE PART
+ 01
-------------
01
+ 11 11 HIGH PART
-------------
*/
// 64-bit product + two 32-bit values
const uint64_t middle = p10 + (p00 >> 32) + (uint32_t)p01;
/*
Proof that 64-bit products can accumulate two more 32-bit values
without overflowing:
Max 32-bit value is 2^32 - 1.
PSum = (2^32-1) * (2^32-1) + (2^32-1) + (2^32-1)
= 2^64 - 2^32 - 2^32 + 1 + 2^32 - 1 + 2^32 - 1
= 2^64 - 1
Therefore it cannot overflow regardless of input.
*/
// 64-bit product + two 32-bit values
r_hi = p11 + (middle >> 32) + (p01 >> 32);
// Add LOW PART and lower half of MIDDLE PART
return (middle << 32) | (uint32_t)p00;
}
#if defined(_MSC_VER) && defined(_WIN64)
// Visual Studio 64-bit
# include <intrin.h>
# pragma intrinsic(_umul128)
# define CAT_MUL128(r_hi, r_lo, x, y) \
r_lo = _umul128(x, y, &(r_hi));
#elif defined(__SIZEOF_INT128__)
// Compiler supporting 128-bit values (GCC/Clang)
# define CAT_MUL128(r_hi, r_lo, x, y) \
{ \
unsigned __int128 w = (unsigned __int128)x * y; \
r_lo = (uint64_t)w; \
r_hi = (uint64_t)(w >> 64); \
}
#else
// Emulate 64x64->128-bit multiply with 64x64->64 operations
# define CAT_MUL128(r_hi, r_lo, x, y) \
r_lo = Emulate64x64to128(r_hi, x, y);
#endif // End CAT_MUL128
这是 ARMv8 或 Aarch64 版本的 asm：
// High (p1) and low (p0) product
uint64_t p0, p1;
// multiplicand and multiplier
uint64_t a = ..., b = ...;
p0 = a*b; asm ("umulh %0,%1,%2" : "=r"(p1) : "r"(a), "r"(b));
这里是旧 DEC 编译器的 asm：
p0 = a*b; p1 = asm("umulh %a0, %a1, %v0", a, b);
如果你有 x86 的 BMI2 并且想使用 mulxq:
asm ("mulxq %3, %0, %1" : "=r"(p0), "=r"(p1) : "d"(a), "r"(b));
通用 x86 乘以 mulq:
asm ("mulq %3" : "=a"(p0), "=d"(p1) : "a"(a), "g"(b) : "cc");
在 C++ 中，表示：
uint64_t i;
uint64_t j;
然后 i * j 将产生一个 uint64_t，其值为 i 和 j 之间的乘积的较低部分，即 (i * j) mod 2^64。
现在，如果我想要乘法的较高部分怎么办？我知道在使用 32 位整数时存在类似的汇编指令，但我对汇编一点都不熟悉，所以我希望得到帮助。
制作如下内容的最有效方法是什么：
uint64_t k = mulhi(i, j);
长乘法性能应该没问题。
将 a*b 拆分为 (hia+loa)*(hib+lob)。这给出了 4 个 32 位乘法加上一些移位。用64位做，手动做进位，你会得到高的部分。
请注意，高部分的近似值可以用更少的乘法来完成——1 次乘法精确到 2^33 左右，3 次乘法精确到 1 以内。
我认为没有可移植的替代品。
如果您使用的是 gcc，并且您的版本支持 128 位数字（尝试使用 __uint128_t），那么执行 128 乘法并提取高 64 位可能是最有效的获取方式结果。
如果您的编译器不支持 128 位数字，那么 Yakk 的回答是正确的。但是，对于一般消费而言，它可能过于简短。特别是，实际的实现必须小心溢出 64 位整数。
他提出的简单且可移植的解决方案是将 a 和 b 中的每一个分解为 2 个 32 位数字，然后使用 64 位乘法运算将这些 32 位数字相乘。如果我们写：
uint64_t a_lo = (uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint32_t)b;
uint64_t b_hi = b >> 32;
那么很明显：
a = (a_hi << 32) + a_lo;
b = (b_hi << 32) + b_lo;
和：
a * b = ((a_hi << 32) + a_lo) * ((b_hi << 32) + b_lo)
= ((a_hi * b_hi) << 64) +
((a_hi * b_lo) << 32) +
((b_hi * a_lo) << 32) +
a_lo * b_lo
前提是使用 128 位（或更高）算法执行计算。
但是这道题需要我们用64位算法进行所有的计算，所以我们不得不担心溢出。
由于 a_hi、a_lo、b_hi 和 b_lo 都是无符号的 32 位数字，它们的乘积将适合无符号的 64 位数字而不会溢出。但是上面计算的中间结果不会。
当数学必须以 2^64 为模执行时，以下代码将实现 mulhi(a, b)：
uint64_t a_lo = (uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint32_t)b;
uint64_t b_hi = b >> 32;
uint64_t a_x_b_hi = a_hi * b_hi;
uint64_t a_x_b_mid = a_hi * b_lo;
uint64_t b_x_a_mid = b_hi * a_lo;
uint64_t a_x_b_lo = a_lo * b_lo;
uint64_t carry_bit = ((uint64_t)(uint32_t)a_x_b_mid +
(uint64_t)(uint32_t)b_x_a_mid +
(a_x_b_lo >> 32) ) >> 32;
uint64_t multhi = a_x_b_hi +
(a_x_b_mid >> 32) + (b_x_a_mid >> 32) +
carry_bit;
return multhi;
正如 Yakk 指出的那样，如果您不介意在高 64 位中被 +1 偏移，则可以省略进位位的计算。
TL:DR 与 GCC 用于 64 位 ISA：(a * (unsigned __int128)b) >> 64 编译得很好，可以编译为单个全乘或高半乘指令。 不需要乱用内联汇编。
不幸的是当前的编译器不优化@craigster0 的便携版本，所以如果你想利用 64 位 CPU ，你不能使用它，除非作为你没有 #ifdef 的目标的后备。 （我没有看到优化它的通用方法；您需要 128 位类型或内部类型。）
大多数 64 位平台上的 GNU C（gcc、clang 或 ICC）has unsigned __int128。 （或者在旧版本中，__uint128_t）。不过，GCC 并未在 32 位平台上实现此类型。
这是让编译器发出 64 位全乘指令并保留高半部分的简单而有效的方法。 （GCC 知道 uint64_t 转换为 128 位整数的上半部分仍然全为零，因此您不会使用三个 64 位乘法得到 128 位乘法。）
MSVC also has a __umulh intrinsic 用于 64 位高半乘法，但同样它仅适用于 64 位平台（特别是 x86-64 和 AArch64。文档还提到 IPF (IA-64) 具有 _umul128 可用，但我没有可用的 Itanium MSVC。（可能无论如何都不相关。）
#define HAVE_FAST_mul64 1
#ifdef __SIZEOF_INT128__ // GNU C
static inline
uint64_t mulhi64(uint64_t a, uint64_t b) {
unsigned __int128 prod = a * (unsigned __int128)b;
return prod >> 64;
}
#elif defined(_M_X64) || defined(_M_ARM64) // MSVC
// MSVC for x86-64 or AArch64
// possibly also || defined(_M_IA64) || defined(_WIN64)
// but the docs only guarantee x86-64! Don't use *just* _WIN64; it doesn't include AArch64 Android / Linux
// https://docs.microsoft.com/en-gb/cpp/intrinsics/umulh
#include <intrin.h>
#define mulhi64 __umulh
#elif defined(_M_IA64) // || defined(_M_ARM) // MSVC again
// https://docs.microsoft.com/en-gb/cpp/intrinsics/umul128
// incorrectly say that _umul128 is available for ARM
// which would be weird because there's no single insn on AArch32
#include <intrin.h>
static inline
uint64_t mulhi64(uint64_t a, uint64_t b) {
unsigned __int64 HighProduct;
(void)_umul128(a, b, &HighProduct);
return HighProduct;
}
#else
# undef HAVE_FAST_mul64
uint64_t mulhi64(uint64_t a, uint64_t b); // non-inline prototype
// or you might want to define @craigster0's version here so it can inline.
#endif
对于 x86-64、AArch64 和 PowerPC64（以及其他），这会编译成一个 mul 指令 ，以及一对 mov 到处理调用约定（在内联之后应该优化掉）。
来自 the Godbolt compiler explorer（使用 x86-64、PowerPC64 和 AArch64 的源代码 + asm）：
# x86-64 gcc7.3. clang and ICC are the same. (x86-64 System V calling convention)
# MSVC makes basically the same function, but with different regs for x64 __fastcall
mov rax, rsi
mul rdi # RDX:RAX = RAX * RDI
mov rax, rdx
ret
（或使用 clang -march=haswell 启用 BMI2：mov rdx, rsi / mulx rax, rcx, rdi 将高半部分直接放入 RAX。gcc 很笨，仍然使用额外的 mov.)
对于 AArch64（使用 gcc unsigned __int128 或使用 __umulh 的 MSVC）：
test_var:
umulh x0, x0, x1
ret
使用编译时常量 2 的乘数，我们通常会得到预期的右移以获取几个高位。但是 gcc 有趣地使用 shld（参见 Godbolt link）。
不幸的是，当前的编译器不优化@craigster0 的便携版本。你得到 8x shr r64,32、4x imul r64,r64 和一堆针对 x86-64 的 add/mov 指令。即它编译成很多 32x32 => 64 位乘法和解包结果。所以如果你想要一些利用 64 位 CPU 的东西，你需要一些 #ifdefs.
一个全乘 mul 64 指令在 Intel CPU 上是 2 微指令，但仍然只有 3 个周期延迟，与 imul r64,r64 相同，它只产生 64 位结果。因此，__int128 / intrinsic 版本在现代 x86-64 上的延迟和吞吐量（对周围代码的影响）比便携式版本便宜 5 到 10 倍，这是基于 http://agner.org/optimize/ 的快速眼球猜测。
在上面 link.
上的 Godbolt 编译器资源管理器中查看gcc 在乘以 16 时确实完全优化了这个函数，但是：你得到一个右移，比 unsigned __int128 乘法更有效。
这是我今晚想出的单元测试版本，提供完整的 128 位产品。经过检查，它似乎比大多数其他在线解决方案（例如 Botan 库和此处的其他答案）更简单，因为它利用了代码注释中解释的中间部分不会溢出的方式。
为了上下文，我为这个 github 项目编写了它：https://github.com/catid/fp61
//------------------------------------------------------------------------------
// Portability Macros
// Compiler-specific force inline keyword
#ifdef _MSC_VER
# define FP61_FORCE_INLINE inline __forceinline
#else
# define FP61_FORCE_INLINE inline __attribute__((always_inline))
#endif
//------------------------------------------------------------------------------
// Portable 64x64->128 Multiply
// CAT_MUL128: r{hi,lo} = x * y
// Returns low part of product, and high part is set in r_hi
FP61_FORCE_INLINE uint64_t Emulate64x64to128(
uint64_t& r_hi,
const uint64_t x,
const uint64_t y)
{
const uint64_t x0 = (uint32_t)x, x1 = x >> 32;
const uint64_t y0 = (uint32_t)y, y1 = y >> 32;
const uint64_t p11 = x1 * y1, p01 = x0 * y1;
const uint64_t p10 = x1 * y0, p00 = x0 * y0;
/*
This is implementing schoolbook multiplication:
x1 x0
X y1 y0
-------------
00 LOW PART
-------------
00
10 10 MIDDLE PART
+ 01
-------------
01
+ 11 11 HIGH PART
-------------
*/
// 64-bit product + two 32-bit values
const uint64_t middle = p10 + (p00 >> 32) + (uint32_t)p01;
/*
Proof that 64-bit products can accumulate two more 32-bit values
without overflowing:
Max 32-bit value is 2^32 - 1.
PSum = (2^32-1) * (2^32-1) + (2^32-1) + (2^32-1)
= 2^64 - 2^32 - 2^32 + 1 + 2^32 - 1 + 2^32 - 1
= 2^64 - 1
Therefore it cannot overflow regardless of input.
*/
// 64-bit product + two 32-bit values
r_hi = p11 + (middle >> 32) + (p01 >> 32);
// Add LOW PART and lower half of MIDDLE PART
return (middle << 32) | (uint32_t)p00;
}
#if defined(_MSC_VER) && defined(_WIN64)
// Visual Studio 64-bit
# include <intrin.h>
# pragma intrinsic(_umul128)
# define CAT_MUL128(r_hi, r_lo, x, y) \
r_lo = _umul128(x, y, &(r_hi));
#elif defined(__SIZEOF_INT128__)
// Compiler supporting 128-bit values (GCC/Clang)
# define CAT_MUL128(r_hi, r_lo, x, y) \
{ \
unsigned __int128 w = (unsigned __int128)x * y; \
r_lo = (uint64_t)w; \
r_hi = (uint64_t)(w >> 64); \
}
#else
// Emulate 64x64->128-bit multiply with 64x64->64 operations
# define CAT_MUL128(r_hi, r_lo, x, y) \
r_lo = Emulate64x64to128(r_hi, x, y);
#endif // End CAT_MUL128
这是 ARMv8 或 Aarch64 版本的 asm：
// High (p1) and low (p0) product
uint64_t p0, p1;
// multiplicand and multiplier
uint64_t a = ..., b = ...;
p0 = a*b; asm ("umulh %0,%1,%2" : "=r"(p1) : "r"(a), "r"(b));
这里是旧 DEC 编译器的 asm：
p0 = a*b; p1 = asm("umulh %a0, %a1, %v0", a, b);
如果你有 x86 的 BMI2 并且想使用 mulxq:
asm ("mulxq %3, %0, %1" : "=r"(p0), "=r"(p1) : "d"(a), "r"(b));
通用 x86 乘以 mulq:
asm ("mulq %3" : "=a"(p0), "=d"(p1) : "a"(a), "g"(b) : "cc");