计算可能溢出的整数运算的最安全,最有效的方法

Question

假设我们有2个常数A&B和一个变量i,所有64位整数.我们想要计算一个简单的通用算术运算,例如:

i * A / B    (1)

为了简化问题,我们假设变量i总是在范围内[INT64_MIN*B/A, INT64_MAX*B/A],因此算术运算(1)的最终结果不会溢出(即:适合该范围[INT64_MIN, INT64_MAX]).

另外,i假设更有可能在友好范围Range1 = [INT64_MIN/A, INT64_MAX/A](即:接近0),但是i可能(不太可能)在该范围之外.在第一种情况下,一个平凡的整数计算i * A不会溢出(这就是我们称之为范围友好的原因); 并且在后一种情况下,i * A将会溢出的平凡整数计算,导致(1)的计算中的错误结果.

什么是"最安全"和"最有效"的计算操作方法(1)(其中"最安全"意味着:保持准确性或至少相当精确,"最有效"意味着:最低平均计算时间),提供i更有可能在友谊范围Range1.

目前,代码中当前实现的解决方案如下:

(int64_t)((double)A / B * i)

哪个解决方案非常安全(没有溢出)虽然不准确(由于双重有效位和53位限制导致的精度损失)并且非常快,因为(double)A / B在编译时预分配了双重除法,只允许在运行时计算双乘法.

Answer 1

如果你不能得到所涉及的范围更好界限,那么你最好关闭以下iammilind的建议来使用__int128.

原因在于,否则你必须实现双字乘法和双字逐字的完整逻辑.英特尔和AMD处理器手册包含有用的信息和现成的代码,但它非常复杂,使用C/C++而不是汇编程序会使事情变得更加复杂.

所有优秀的编译器都将有用的原语作为内在函数.微软的列表似乎不包含类似muldiv的原语,但__mul128内在函数将128位产品的两半作为两个64位整数.基于此,您可以执行两位数的长除以一位数,其中一个"数字"将是64位整数(通常称为"肢体",因为大于数字但仍然只是整数的一部分).仍然相当复杂,但比使用纯C/C++要好很多.但是,便携性方面并不比__int128直接使用更好.至少在这种情况下,编译器实现者已经为您完成了所有艰苦的工作.

如果您的应用程序域可以为您提供有用的边界,(u % d) * v那么不会溢出,那么您可以使用该标识

(u * v) / d = (u / d) * v + ((u % d) * v) / d

其中,/表示整数除法,只要u为非负数且d为正数(否则您可能会与运算符的语义允许的余地相冲突%).

在任何情况下,您可能必须分离操作数的符号并使用无符号运算,以便找到可以利用的更有用的机制 - 或者避免编译器的破坏,例如您提到的饱和乘法.有符号整数运算的溢出会调用未定义的行为,编译器可以随心所欲地执行任何操作.相比之下,无符号类型的溢出是明确定义的.

此外,对于无符号类型,您可以依赖于这样的规则s = a (+) b(其中(+)可能是溢出的无符号加法),您将拥有s == a + b或者s < a && s < b,这可以让您在事后使用廉价操作检测溢出.

但是,你不太可能在这条道路上走得更远,因为所需的努力很快接近 - 甚至超过 - 我实施前面提到的双肢操作的努力.只有对应用程序域进行全面分析才能提供规划/部署此类快捷方式所需的信息.在一般情况下,你给你的界限几乎没有运气.

Answer 2

为了提供问题的量化答案，我对不同解决方案进行了基准测试，作为本文提出的解决方案的一部分（感谢评论和答案）。

基准测试测量不同实现的计算时间，何时在友好范围Range1 =i内，何时在友好范围之外（但在安全范围Range2 =内）。[INT64_MIN/A, INT64_MAX/A]i[INT64_MIN*B/A, INT64_MAX*B/A]

每个实现都执行操作的“安全”（即没有溢出）计算：（i * A / B第一个实现除外，作为参考计算时间给出）。然而，某些实现可能会返回不常见的不准确的计算结果（通知哪些行为）。

提出的一些解决方案尚未经过测试或未在下文列出；这些是：使用解决方案__int128（ms vc编译器不支持），但int128_t已使用boost；使用扩展 80 位的解决方案long double（ms vc 编译器不支持）；使用解决方案InfInt（工作和测试过，但速度太慢，无法成为像样的竞争对手）。

时间测量以 ps/op（每次操作皮秒）为单位指定。基准平台是 Windows 7 x64 下的 Intel Q6600@3GHz，使用 MS vc14、x64/Release 目标编译的可执行文件。以下引用的变量、常量和函数定义为：

int64_t       i;
const int64_t A     = 1234567891;
const int64_t B     = 4321987;
inline bool   in_safe_range(int64_t i) { return (INT64_MIN/A <= i) && (i <= INT64_MAX/A); }

1. (i * A / B) [参考] Range1内的
   i：1469 ps/op，Range1外的i：不相关（溢出）
   
2. ((int64_t)((double)i * A / B))
   Range1内的i：10613 ps/op，Range1外的i：10606 ps/op 注意：在整个Range2范围内，不常见的不准确结果（最大误差 = 1 位）
   
   
3. ((int64_t)((double)A / B * i))
   Range1内的i：1073 ps/op，Range1外的i：1071 ps/op 注意：整个范围Range2中很少出现不准确的结果（最大错误 = 1 位） 注意：编译器可能会进行预计算，从而导致与以前的解决方案相比所观察到的性能提升。
   
   (double)A / B
   
4. (!in_safe_range(i) ? (int64_t)((double)A / B * i) : (i * A / B))
   i在Range1中：2009 ps/op，i在Range1之外：1606 ps/op
   注意：在Range1之外很少有不准确的结果（最大误差 = 1 位）
   
5. ((int64_t)((int128_t)i * A / B)) [boost int128_t]
   i在Range1中：89924 ps/op，i在Range1之外：89289 ps/op
   注意：boostint128_t在工作台平台上的表现非常糟糕（不知道为什么）
   
6. ((i / B) * A + ((i % B) * A) / B)
   i在Range1内：5876 ps/op，i在Range1外：5879 ps/op
   
7. (!in_safe_range(i) ? ((i / B) * A + ((i % B) * A) / B) : (i * A / B))
   i在Range1内：1999 ps/op，i在Range1外：6135 ps/op
   

结论
a) 如果在整个Range2范围内轻微的计算误差是可以接受的，那么解决方案（3）是最快的，甚至比作为参考给出的直接整数计算更快。
b) 如果计算误差在友好范围Range1内是不可接受的，但在该范围之外是可接受的，则解决方案（4）是最快的。
c) 如果计算误差在整个范围Range2中是不可接受的，则解决方案（7）在友好范围Range1中的表现与解决方案（4）一样，并且在该范围之外保持相当快的速度。