如何准确地乘以和除以64位整数？

我有一个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的整个产品，我只会在该类型中进行数学运算然后进行截断，但是有一些方法可以在没有大整数的情况下执行此操作吗？

Answer 1

用|r| < |d|写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; 

Answer 2

我建议用a，b和c中的每一个找到d的最大公约数，除以常用因子：

 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算法以对数时间运行，因此这应该是相当有效的。

Answer 3

很容易看出一些输入会产生无法用返回类型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); } 

Answer 4

如果您正在使用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指令完全以双倍长度工作以允许溢出。 可悲的是，我不知道使用它们的内在编译器。