获得64位整数乘法的高分

Question

在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);

Answer 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位数而不会溢出.但是,上述计算的中间结果不会.

以下代码将实现mulhi(a,b),当必须以2 ^ 64模数执行数学化验时:

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关闭,你可以省略进位的计算.

Answer 2

这是我今晚提出的经过单元测试的版本，提供完整的 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

Answer 3

用于 64 位 ISA 的带有 GCC 的 TL:DR：可以(a * (unsigned __int128)b) >> 64很好地编译为单个全乘法或高半乘法指令。 无需搞乱内联汇编。

不幸的是，当前的编译器没有优化@craigster0 的便携版本，所以如果你想利用 64 位 CPU，你不能使用它，除非作为你没有 for 的目标的后备#ifdef。（我没有看到优化它的通用方法；您需要 128 位类型或内在类型。）

GNU C（gcc、clang 或 ICC）在大多数 64 位平台上都有unsigned __int128。（或在旧版本中，__uint128_t）。不过，GCC 不会在 32 位平台上实现这种类型。

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

MSVC 也有一个__umulh64 位高半乘法的内在函数，但它同样只在 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指令，并编译成一对movs 来处理调用约定（应该在此内联之后优化掉）。来自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（使用 gccunsigned __int128或 MSVC __umulh）：

test_var:
    umulh   x0, x0, x1
    ret

使用 2 倍乘法器的编译时常数幂，我们通常会得到预期的右移以获取一些高位。但是 gcc 使用shld得很有趣（参见 Godbolt 链接）。

不幸的是，当前的编译器没有优化@craigster0 的便携版本。你得到 8x shr r64,32、 4ximul r64,r64和一堆x86-64的add/mov指令。即它编译为很多 32x32 => 64 位乘法和结果的解包。因此，如果您想要利用 64 位 CPU 的功能，则需要一些#ifdefs。

一条完整的乘法mul 64指令在 Intel CPU 上是 2 uop，但仍然只有 3 个周期的延迟，与imul r64,r64仅产生 64 位结果相同。因此，__int128根据基于http://agner.org/optimize/的快速眼球猜测，/ 内在版本在现代 x86-64 上的延迟和吞吐量（对周围代码的影响）比便携式版本便宜 5 到 10 倍。

在上述链接的 Godbolt 编译器资源管理器中查看。

不过，gcc 在乘以 16 时确实完全优化了这个函数：你得到一个右移，比unsigned __int128乘法更有效。

Answer 4

长乘法性能应该不错。

分裂a*b成(hia+loa)*(hib+lob). 这给出了 4 个 32 位乘法加上一些移位。以 64 位进行，并手动进行进位，您将得到高位部分。

请注意，高部分的近似可以通过更少的乘法来完成 - 1 次乘法的精确度在 2^33 左右，而 3 次乘法的精确度在 1 以内。

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