正如该问题的评论者指出的那样,这最容易通过机器相关代码而不是可移植代码有效地完成。提问者说主平台是 x86_64,并且有一个用于执行 64 ✕ 64 → 128 位乘法的内置指令。这很容易使用一小块内联汇编来访问。请注意,内联汇编的细节可能与编译器有所不同,下面的代码是使用英特尔 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;
}