获取 64 位整数乘法的高部分

Getting the high part of 64 bit integer multiplication

本文关键字:高部 整数 获取      更新时间:2023-10-16

在C++中,说:

uint64_t i;
uint64_t j;

然后i * j将产生一个uint64_t,其值为 ij 之间乘法的下部,即 (i * j) mod 2^64 。现在,如果我想要乘法的较高部分怎么办?我知道在使用 32 位整数时存在一个汇编指令可以做类似的事情,但我根本不熟悉汇编,所以我希望得到帮助。

制作以下内容的最有效方法是什么:

uint64_t k = mulhi(i, j);

如果您使用的是 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很好地编译为单个全乘法或高半乘法指令。 无需弄乱内联 asm。


不幸的是,当前的编译器并没有优化@craigster0漂亮的便携式版本,所以如果你想利用 64 位 CPU,除了作为你没有#ifdef的目标的后备之外,你不能使用它。 (我没有看到优化它的通用方法;你需要一个 128 位类型或内部函数。


GNU C(gcc,clang或ICC(在大多数64位平台上都有unsigned __int128。 (或在旧版本中,__uint128_t (。 但是,GCC 不会在 32 位平台上实现此类型。

这是让编译器发出 64 位全乘法指令并保持高半部分的一种简单有效的方法。 (GCC 知道转换为 128 位整数的uint64_t仍然具有上半部分全部为零,因此您不会使用三个 64 位乘法获得 128 位乘法。

MSVC 还具有 64 位高半乘法的__umulh内在函数,但它仅在 64 位平台上可用(特别是 x86-64 和 AArch64(。 文档还提到IPF(IA-64(有可用的_umul128,但我没有用于安腾的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://learn.microsoft.com/en-gb/cpp/intrinsics/umulh
  #include <intrin.h>
  #define mulhi64 __umulh
#elif defined(_M_IA64) // || defined(_M_ARM)       // MSVC again
  // https://learn.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来处理调用约定(在此内联之后应该会优化(。 来自 Godbolt 编译器资源管理器(使用 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 或 MSVC 带__umulh(:

test_var:
    umulh   x0, x0, x1
    ret

在编译时常功率为 2 乘法器的情况下,我们通常会得到预期的右移来抓取几个高位。 但是gcc很有趣地使用了shld(参见Godbolt链接(。


不幸的是,当前的编译器并没有优化@craigster0漂亮的便携式版本。 你会得到 8x shr r64,32 , 4x imul r64,r64 ,以及一堆 x86-64 的 add/mov 指令。 即它编译为很多 32x32 => 64 位乘法和解压缩结果。 因此,如果您想要利用 64 位 CPU 的东西,您需要一些#ifdef

在英特尔 CPU 上,全乘法mul 64指令为 2 uops,但仍然只有 3 个周期延迟,与仅产生 64 位结果的 imul r64,r64 相同。 因此,根据基于 http://agner.org/optimize/的快速眼球猜测,现代 x86-64 的 __int128/内部版本在延迟和吞吐量(对周围代码的影响(方面的延迟和吞吐量(对周围代码的影响(比便携式版本便宜 1 到 64 倍。

在上面的链接上的 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

长乘法应该没问题。

a*b拆分为(hia+loa)*(hib+lob)。 这给出了 4 个 32 位乘法加上一些移位。 以 64 位执行它们,并手动执行进位,您将获得高位。

请注意,高部分的近似可以用更少的乘法来完成 - 1 次乘法在 2^33 左右,3 次乘法在 1 以内准确。

我认为没有便携式替代品。

以下是 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");