方程中的多重/除法困境答案

【问题标题】：Multiple/Division dilemma in equation方程中的多重/除法困境
【发布时间】：2012-10-20 13:56:43
【问题描述】：

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

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

使用 __int64 类型的 a,b,c,x (a,b,c,x 但是，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。

【问题讨论】：

执行 128 位算术作为四个 32 位数字？这不是那么难。甚至是汇编程序：MOV, MUL, long decrement, DIV.
他们（VC 等）还没有int128？

标签： c++ c

【解决方案1】：

您可以手动实现长算法，使用 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)

【讨论】：

我最初把它写成算法而不是实现。 [ab] 表示 a*2^32+b（以 base-2^32 串联）
@user1704182 请检查第一部分 - 我修复了一个“错误”
@Serge 我将“long”用作“两个 32 位字”，而不是 C++ 意义上的“long”；-)
@Jan any long 在结果中右移 32 位将产生 0。并且您将在计算过程中丢失进位。所以 e,f,g 必须是 long long（或者 VC 中代表 64 位整数的任何类型）
一个问题是 r 可能大到足以变成 64 位并且 [rg]=[yzg]。在这种情况下 l=[yzg]/[hi]。我们回到原来的问题

【解决方案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

【讨论】：

当我添加此代码时，编译器投诉：错误 C2415：不正确的操作数类型，在 2 行包含 RAX。我认为VS不允许RAX。
你是为 64 位目标构建这个吗？
@user1704182 即使你指定了目标CPU？
@user1704182 检查 Visual Studio 中的配置管理器以获取构建目标；参看。 stackoverflow.com/questions/5343814/…
所以你想定位 32 位机器？那我的代码就不行了。

【解决方案3】：

尝试计算 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

【讨论】：

浮点运算可能会引入非一错误。 1.0/c 在我看来就像一个浮点运算。

【解决方案4】：

你也许可以重新安排你的计算：

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

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

例如，如果a 可能是从 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（来自第一次重新排列）会比(256 * b) / c 损失8 位的精度。

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

【讨论】：