8

我需要用 C/C++ 计算这个方程： 

x=(a*b-1)/c;

使用 __int64 类型的 a,b,c,x (a,b,c,x<10^13)。在所有情况下，选择 a,b,c 以使 x 适合 __int64。
但是，a*b 非常大，会导致溢出，x 是错误的。
我尝试通过类型转换来分隔 a*b： 

x=(__int64)(((double)a/c)*(double)b - 1.0/c);

这样，首先计算 a/c 并且不会发生溢出错误。
但是，问题是 ((double)a/c)*(double)b​​ 有时值很大（大约数十亿）并且精度降低，因此 1.0/c（非常小）不起作用并导致错误+-1。

例如：(__int64)(((double)a/c)*(double)b​​=123456789.01更有可能变成123456789.0和1.0/c=0.02。这种情况下有+1的误差。

有没有办法在没有外部库（如 Boost 或 Bignum）的情况下计算 x 精确？即使有错误 +-1 也会搞砸我的代码。
提前致谢。
顺便说一句，我使用 Visual Studio 10。

4

4 回答 4

3

您可以手动实现长算法，使用 32 位块：

长乘法：

  (ax + b) * (cx + d)
= ac x^2 + (ad+bc) x + bd

  [ab]*[cd]=[efg]:

   //long long e,f,g
   g = (long long)b*d;
   f = (long long)a*d+b*c;
   e = (long long)a*c;
   //perform carry
   f += g>>32; g &= 0xFFFFFFFF
   e += f>>32; f &= 0xFFFFFFFF

学分，假设无符号算术：

   [efg]/[hi]=[jkl]:

    [jk] = [ef]/[hi];
    r = [ef]-j*[hi];
    l = [rg]/[hi];
    if j > 0, result doesn't fit
    x = [kl];

如果a和b有符号，则首先修复符号并使用绝对值进行计算，如@Serge 所建议的那样：如果 a 或 b 为零，x=(-1)/c 否则，sign(x)=sign(a)*sign(b)*sign(c)

于 2012-10-20T04:01:28.717 回答
2

如果您的代码可能依赖于 CPU，最简单的方法可能是使用汇编程序来保留高位 8 个字节。x64，假设结果适合 8 个字节：

__asm{
  MOV RAX, a
  MUL b
  SUB RAX, 1
  SBB RDX, 0
  DIV c
  MOV x, RAX
}

[1] http://en.wikibooks.org/wiki/X86_Assembly/Arithmetic

于 2012-10-20T03:03:23.947 回答
0

尝试找出 a 的最大乘数，它不会溢出。例如，如果 a*4 会溢出，则部分执行：

(a*5) = (a*3) + (a*2)

因此，如果您发现中等值，例如

b1, b2, b3, where b1+b2+b3 == b

然后

x = a*b1/c + a*b2/c + a*b3/c - 1.0/c

您将从 b1 = floor(MAX_INT / a) 中找到最大的可用值，然后 b2 将是 b 的其余部分，仅当 b-b1 < b1 如果不是，则必须重复。

于 2012-10-20T03:04:40.433 回答
0

您也许可以重新排列您的计算：

x = (a*b - 1) / c
  = a*b/c - 1 / c
  = (a/c)*b - 1 / c       ;If a is likely to be huge
  = a*(b/c) - 1 / c       ;If b is likely to be huge

问题将是精度损失。这可以通过跟踪值的范围和引入缩放因子（选择使中间值尽可能大而不会溢出）来最小化。

例如，ifa可以是0到1024，b可以是0到16777215，c也可以是4096到8192；然后：

x  = (a * ((256*b) /c) - 256 / c) / 256

在这种情况下(256*b) <= 0xFFFFFF00（尽可能大以避免溢出 32 位无符号整数）( (256*b) / c) <= 0x000FFFFF、 和a * ((256*b) /c) <= 0x3FFFFC00.

同样对于这种情况，a*b（根据原始公式）将溢出一个 32 位无符号整数；和b/c（从第一次重新排列开始）会比现在损失 8 位的精度(256 * b) / c。

当然，针对您的具体情况的最佳公式（在没有溢出的情况下给出最小精度损失的公式）取决于您具体情况下可能的变量范围。

于 2012-10-20T05:56:43.783 回答