用 C 计算 64x64 int 乘积的高 64 位

Question

我希望我的 C 函数能够有效地计算两个 64 位有符号整数乘积的高 64 位。我知道如何在 x86-64 程序集中使用 imulq 并将结果从 %rdx 中提取出来。但是我完全不知道如何用 C 来编写它，更不用说哄编译器高效地完成它了。

有人对用 C 编写这个有什么建议吗？这对性能很敏感，因此“手动方法”（如俄罗斯农民或 bignum 库）已被淘汰。

我写的这个笨拙的内联汇编函数可以工作，大致就是我所追求的代码生成：

static long mull_hi(long inp1, long inp2) {
    long output = -1;
    __asm__("movq %[inp1], %%rax;"
            "imulq %[inp2];"
            "movq %%rdx, %[output];"
            : [output] "=r" (output)
            : [inp1] "r" (inp1), [inp2] "r" (inp2)
            :"%rax", "%rdx");
    return output;
}

Answer 1

如果您在 x86_64 上使用相对较新的 GCC：

int64_t mulHi(int64_t x, int64_t y) {
    return (int64_t)((__int128_t)x*y >> 64);
}

在 -O1 及更高版本，这将编译为您想要的：

_mulHi:
0000000000000000    movq    %rsi,%rax
0000000000000003    imulq   %rdi
0000000000000006    movq    %rdx,%rax
0000000000000009    ret

我相信 clang 和 VC++ 也支持 __int128_t 类型，因此这也应该适用于这些平台，但通常需要注意的是自己尝试一下。

Answer 2

一般的答案是x * y 可以分解为(a + b) * (c + d)，其中a 和c 是高阶部分。

首先，扩展为ac + ad + bc + bd

现在，您将这些项乘以 32 位数字，存储为 long long（或者更好的是 uint64_t），您只需要记住，当您乘以更高的阶数时，您需要按 32 位进行缩放。然后你做加法，记住检测进位。跟踪标志。当然，您需要分段添加。

有关实现上述内容的代码，请参阅my other answer。

Answer 3

关于您的组装解决方案，不要硬编码mov 指令！让编译器为你做这件事。这是您的代码的修改版本：

static long mull_hi(long inp1, long inp2) {
    long output;
    __asm__("imulq %2"
            : "=d" (output)
            : "a" (inp1), "r" (inp2));
    return output;
}

有用的参考：Machine Constraints

Answer 4

由于您在解决自己的机器代码问题方面做得非常好，我认为您应该在可移植版本方面得到一些帮助。如果在 x86 上的 gnu 中，我会留下一个 ifdef 在你只使用程序集的地方。

无论如何，这是基于my general answer 的实现。我很确定这是正确的，但不能保证，我昨晚刚敲了这个。您可能应该摆脱静态 positive_result[] 和 result_negative - 这些只是我的单元测试的人工制品。

#include <stdlib.h>
#include <stdio.h>

// stdarg.h doesn't help much here because we need to call llabs()

typedef unsigned long long uint64_t;
typedef   signed long long  int64_t;

#define B32 0xffffffffUL

static uint64_t positive_result[2]; // used for testing
static int result_negative;         // used for testing

static void mixed(uint64_t *result, uint64_t innerTerm)
{
  // the high part of innerTerm is actually the easy part

    result[1] += innerTerm >> 32;

  // the low order a*d might carry out of the low order result

    uint64_t was = result[0];

    result[0] += (innerTerm & B32) << 32;

    if (result[0] < was) // carry!
      ++result[1];
}


static uint64_t negate(uint64_t *result)
{
  uint64_t t = result[0] = ~result[0];
  result[1] = ~result[1];
  if (++result[0] < t)
    ++result[1];
  return result[1];
}

uint64_t higherMul(int64_t sx, int64_t sy)
{
    uint64_t x, y, result[2] = { 0 }, a, b, c, d;

    x = (uint64_t)llabs(sx);
    y = (uint64_t)llabs(sy);

    a = x >> 32;
    b = x & B32;
    c = y >> 32;
    d = y & B32;

  // the highest and lowest order terms are easy

    result[1] = a * c;
    result[0] = b * d;

  // now have the mixed terms ad + bc to worry about

    mixed(result, a * d);
    mixed(result, b * c);

  // now deal with the sign

    positive_result[0] = result[0];
    positive_result[1] = result[1];
    result_negative = sx < 0 ^ sy < 0;
    return result_negative ? negate(result) : result[1];
}

Answer 5

等等，您已经有了一个完美的优化装配解决方案 为此工作，您想将其退出并尝试将其写入 不支持 128 位数学的环境？我没有关注。

您显然知道，此操作是 x86-64。显然，您所做的任何事情都不会使它工作得更好。 如果你真的想要便携式 C，你需要做类似的事情 上面的 DigitalRoss 的代码，希望你的优化器弄清楚什么 你在做。

如果您需要架构可移植性但愿意限制自己 对于 gcc 平台，有 __int128_t（和 __uint128_t）类型 编译器内部函数可以满足您的需求。