将两个64位整数相乘并将结果右移至64位的最快方法是什么？

Question

将两个64位整数相乘并将结果右移至64位的最快方法是什么？

4

我需要将两个带符号的64位整数a和b相乘，然后将（128位）结果向右移动以得到一个带符号的64位整数。有什么最快的方法可以做到这一点？

我的64位整数实际上表示具有fmt小数位的定点数。 fmt被选择为使a * b >> fmt不会溢出，例如，对于fmt == 56，abs（a）<64 << fmt和abs（b）<2 << fmt不会在64位中溢出，因为最终结果应该是<128 << fmt，因此适合int64。

我想这样做的原因是快速准确地评估形式为((((c5 * x + c4) * x + c3) * x + c2) * x + c1) * x + c0的五次多项式，以固定点格式输出，每个数字都是带符号的 64位固定点数字，并且有fmt小数位。我正在寻找实现此目标的最有效方法。

- Michel Rouzic

2

你提出问题的陈述表明你可能已经尝试过实现。如果是这样，你能否发布你的代码？ - ryyker

@Oliver Charlesworth 这段代码是为了可移植性而设计的，我不知道是否有广泛可用的int128实现。我想，既然无论编译器可能做什么，我都可以在没有int128类型的情况下完成，那么不需要int128类型的解决方案应该是可行的，对吗？我认为，通过进行位移来获得int64结果可能会带来一些巧妙的技巧。 - Michel Rouzic

1

了解ISA的一些信息会很有帮助。通常编写非可移植代码会更容易。 - user3528438

@user3528438 好的，一般来说是现代 PC，主要是 x86_64 架构。如果需要的话，我可以用非便携式方式完成，然后再提供一个便携式的备选方案。 - Michel Rouzic

看这里的建议，使用SSE4，https://dev59.com/gnTYa4cB1Zd3GeqPz-b3 - Jens Munk

显示剩余3条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- njuffa · Accepted Answer

正如一个评论者指出的那样，这最容易通过机器相关代码高效地实现，而不是通过可移植代码。提问者表示主要平台是x86_64，并且该平台具有执行64 ✕ 64→128位乘法的内置指令。这可以使用一小段内联汇编轻松访问。请注意，内联汇编的细节可能会因编译器而异，下面的代码是使用Intel C/C++编译器构建的。

#include <stdint.h>

/* compute mul_wide (a, b) >> s, for s in [0,63] */
int64_t mulshift (int64_t a, int64_t b, int s)
{
    int64_t res;
    __asm__ (
        "movq  %1, %%rax;\n\t"          // rax = a
        "movl  %3, %%ecx;\n\t"          // ecx = s
        "imulq %2;\n\t"                 // rdx:rax = a * b
        "shrdq %%cl, %%rdx, %%rax;\n\t" // rax = int64_t (rdx:rax >> s)
        "movq  %%rax, %0;\n\t"          // res = rax
        : "=rm" (res)
        : "rm"(a), "rm"(b), "rm"(s)
        : "%rax", "%rdx", "%ecx");
    return res;
}

下面展示了一个便携式的C99版本，与上述代码等效。我已经对其进行了广泛测试，没有发现不匹配的情况。

void umul64wide (uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo)
{
    uint64_t a_lo = (uint64_t)(uint32_t)a;
    uint64_t a_hi = a >> 32;
    uint64_t b_lo = (uint64_t)(uint32_t)b;
    uint64_t b_hi = b >> 32;

    uint64_t p0 = a_lo * b_lo;
    uint64_t p1 = a_lo * b_hi;
    uint64_t p2 = a_hi * b_lo;
    uint64_t p3 = a_hi * b_hi;

    uint32_t cy = (uint32_t)(((p0 >> 32) + (uint32_t)p1 + (uint32_t)p2) >> 32);

    *lo = p0 + (p1 << 32) + (p2 << 32);
    *hi = p3 + (p1 >> 32) + (p2 >> 32) + cy;
}

void mul64wide (int64_t a, int64_t b, int64_t *hi, int64_t *lo)
{
    umul64wide ((uint64_t)a, (uint64_t)b, (uint64_t *)hi, (uint64_t *)lo);
    if (a < 0LL) *hi -= b;
    if (b < 0LL) *hi -= a;
}

/* compute mul_wide (a, b) >> s, for s in [0,63] */
int64_t mulshift (int64_t a, int64_t b, int s)
{
    int64_t res;
    int64_t hi, lo;
    mul64wide (a, b, &hi, &lo);
    if (s) {
        res = ((uint64_t)hi << (64 - s)) | ((uint64_t)lo >> s);
    } else {
        res = lo;
    }
    return res;
}