2

我有一个 C 函数：

int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d)
{
    /* should return (a * b * c)/d */   
}

a 可能接近 INT64_MAX，但最终结果不会溢出，例如如果 b = 1，c = d = 40。但是，我无法弄清楚如何计算它，这样我就不会丢失数据舍入（通过先进行除法）或中间结果溢出。

如果我可以访问一个足够大的数据类型来容纳 a、b 和 c 的整个乘积，我会在该类型中进行数学运算然后截断，但是有没有什么方法可以在没有大整数的情况下做到这一点？

4

4 回答 4

5

写a = q*d + r与|r| < |d|（我假设d != 0，否则计算无论如何都是没有意义的）。然后(a*b*c)/d = q*b*c + (r*b*c)/d。如果q*b*c溢出，整个计算无论如何都会溢出，所以要么你不在乎，要么你必须检查溢出。r*b*c可能仍然会溢出，所以我们再次使用相同的方法来避免溢出，

int64_t q = a/d, r = a%d;
int64_t part1 = q*b*c;
int64_t q1 = (r*b)/d, r1 = (r*b)%d;
return part1 + q1*c + (r1*c)/d;
于 2013-08-02T17:38:40.460 回答
1

很容易看出，某些输入会产生无法由返回类型表示的输出int64_t。例如，fn(INT64_MAX, 2, 1, 1)。但是，以下方法应该允许您为实际上适合范围内的任何输入组合返回正确答案int64_t。

int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d)
{
    /* find the integer and remainder portions of a/d */
    int64_t leftI = a / d;
    int64_t leftR = a % d;

    /* multiply the integer portion of the result by b and c */
    int64_t resultI = leftI * b * c;

    /* multiply the remainder portion by b */
    int64_t resultR = leftR * b;
    resultI = resultI + (resultR / d) * c;

    /* multiply the remainder portion by c */
    resultR = (resultR % d) * c;

    return resultI + (resultR / d);
}
于 2013-08-02T17:38:24.427 回答
1

我建议找到 d 与 a、b 和 c 中的每一个的最大公约数，并在进行过程中划分公因数：

common = gcd(a,d) // You can implement GCD using Euclid's algorithm

a=a/common
d=d/common

common = gcd(b,d)
b=b/common
d=d/common

common = gcd(c,d)
c=c/common
d=d/common

然后a*b*c/d在去除所有公因数的情况下进行计算。Euclid 的 GCD 算法以对数时间运行，所以这应该是相当有效的。

于 2013-08-02T17:46:42.433 回答
0

如果您使用的是 x86_64，那么 asm 支持 128 位整数：

int64_t fn(uint64_t a, uint64_t b, uint64_t c, uint64_t d) {

    asm (
        "mulq %1\n"          // a *= b
        "movq %%rbx, %%rdx\n"// rbx = upper 64 bit of the multiplication
        "mulq %2\n"          // multiply the lower 64 bits by c
        "push %%rax\n"       // temporarily save the lowest 64 bits on the stack
        "mov %%rcx, %%rdx\n" // rcx = upper 64 bits of the multiplication
        "movq %%rax, %%rbx\n"// 
        "mulq %2\n"          // multiply the upper 64 bits by c
        "addq %%rax, %%rcx\n"// combine the middle 64 bits
        "addcq %%rdx, $0\n"  // transfer carry tp the higest 64 bits if present
        "divq %3\n"          // divide the upper 128 (of 192) bits by d
        "mov %%rbx, %%rax\n" // rbx = result
        "pop %%rax\n"
        "divq %3\n"          // divide remainder:lower 64 bits by d
        : "+a" (a)           // assigns a to rax register as in/out
        , "+b" (b)           // assigns b to rbx register
        : "g" (c)            // assigns c to random register
        , "g" (d)            // assigns d to random register
        : "edx", "rdx"       // tells the compiler that edx/rdx will be used internally, but does not need any input
    );

    // b now holds the upper 64 bit if (a * b * c / d) > UINT64_MAX
    return a;
}

请注意，所有输入整数的长度必须相同。工作长度将是输入的两倍。仅适用于未签名。

x86 上的本机指令div和mul指令完全使用双倍长度以允许溢出。可悲的是，我不知道使用它们的内在编译器。

于 2013-08-02T17:39:04.850 回答