随机排序的 IEEE 754 双精度浮点数之和的结果答案

【问题标题】：Result of the sum of random-ordered IEEE 754 double precision floats随机排序的 IEEE 754 双精度浮点数之和的结果
【发布时间】：2018-04-28 00:42:25
【问题描述】：

这是我的问题的伪代码。

我有一个 IEEE 754 双精度正数数组。

数组可以以随机顺序出现，但数字始终相同，只是在它们的位置上打乱了。此外，这些数字在double 表示的有效 IEEE 范围内可以在非常宽的范围内变化。

获得列表后，我会初始化一个变量：

double sum_result = 0.0;

我在sum_result 上累积总和，在整个数组的循环中。在我做的每一步：

sum_result += my_double_array[i]

是否保证无论double的初始数组的顺序如何，如果数字相同，打印出来的总和结果总是相同的？

【问题讨论】：

请给我们一个minimal reproducible example，而不是伪代码，其中包含示例输入和预期与当前输出。
我几乎可以保证结果几乎总是不同的。
不，不能保证。例如，如果所有大值都在前，那么小值可能会完全下降。
我不能分享代码，但对我来说不需要，因为这是一个纯粹的理论问题，而不是调试请求。
@user0042 - 我认为这是对伪代码的完全合法使用，如果您阅读整个问题，您会发现没有具体示例就可以理解。

标签： floating-point precision ieee-754

【解决方案1】：

没有。

举个简单的例子，0x1p53 加 1 得到 0x1p53。（这使用十六进制浮点表示法。“p”之前的部分是有效数，以十六进制表示，与 C 十六进制整数常量相同，除了它可能有一个“.”来标记小数的开始部分。“p”后面的数字表示有效数乘以的 2 的幂。）这是因为数学上精确的和 0x1.00000000000008p+53 不能用 IEEE-754 64 位二进制浮点数表示, 所以它被四舍五入到最接近的值，其有效数字的低位为偶数，即 0x1p53。

因此，0x1p53+1 产生 0x1p53。因此，从左到右计算的 0x1p53+1+1 也产生 0x1p53。但是1+1就是2，而2+0x1p53正好可以表示为0x1.0000000000001p+53，所以1+1+0x1p53就是0x1.0000000000001p+53。

为了以十进制显示更易于可视化的示例，假设我们只有两个十进制数字。然后 100+1 产生 100，因此 100+1+1+1+1+1+1 产生 100。但是 1+1+1+1+1+1+100 累积到 6+100，然后产生 110（由于四舍五入到两位有效数字）。

【讨论】：

所以，似乎每次出现都获得相同结果的唯一方法是先对数组进行排序。
@auserdude：double 值中的潜在位跨度从 0x1p-1074（次正规值的最低位）到 0x1p1023（最大有限值的最高位）。因此，使用包含 1023 - -1074+1 = 2098 位的数据结构，您可以构造一个能够精确执行必要算术的加法器，前提是没有溢出。
是的，我知道表示的扩展，但如果我不能创建这样的结构，似乎唯一的方法是对数组进行排序，以防需要确定性结果，对吧？
@auserdude：您可能还会看到 github.com/python/cpython/blob/… 用于正确舍入（因此是顺序不变）的求和算法

【解决方案2】：

是否保证，无论double的初始数组的顺序如何，如果数字相同，打印出来的sum结果总是相同的？

不，FP 添加不是associative。请记住，它被称为 floating 点 - 绝对精度“浮动”大约相对于 1.0。任何给定的operation 类似添加 (+) 都受round-off error 约束。

但是，如果总和完成并且 inexact 标志是明确的，那么是的，顺序不相关。**

简单的反例。

#include <math.h>
#include <float.h>
#include <fenv.h>
#include <stdio.h>

int main(void) {
  double a[3] = { DBL_MAX, -DBL_MAX, 1.0 };
  fexcept_t flag;

  feclearexcept(FE_ALL_EXCEPT);
  printf("%e\n", (a[0] + a[1]) + a[2]);
  fegetexceptflag(&flag, FE_INEXACT);
  printf("Inexact %d\n", !!(flag & FE_INEXACT));

  feclearexcept(FE_ALL_EXCEPT);
  printf("%e\n", a[0] + (a[1] + a[2]));
  fegetexceptflag(&flag, FE_INEXACT);
  printf("Inexact %d\n", !!(flag & FE_INEXACT));

  printf("%d\n", FLT_EVAL_METHOD);
  return (EXIT_SUCCESS);
}

输出

1.000000e+00  // Sum is exact
Inexact 0

0.000000e+00  // Sum is inexact
Inexact 1

0    // evaluate all operations ... just to the range and precision of the type;

根据FLT_EVAL_METHOD，FP 数学可能会使用更宽的进动和范围，但上述极端示例的总和仍会有所不同。

** 除了可能是 0.0 与 -0.0 的结果

要了解原因，请尝试使用 4 位精度的基于 10 个文本的示例。同样的原则也适用于double，它通常具有 53 位二进制精度。

a[3] = +1.000e99, -1.000e99, 1.000
sum = a[0] + a[1]   // sum now exactly 0.0 
sum += a[2]         // sum now exactly 1.0 
// vs.
sum = a[1] + a[2]   // sum now inexactly -1.000e99
sum += a[0]         // sum now inexactly 0.0

Re：“打印出的总和结果将始终相同”：除非代码以足够高的精度打印 "%a" 或 "%.*e"，否则打印的文本可能缺乏意义和两个不同的总和可能看起来一样。见Printf width specifier to maintain precision of floating-point value

【讨论】：

【解决方案3】：

让我们举个例子：为了简单起见，我正在使用以 10 为底的模型转置浮点问题，只有 2 个有效数字，运算结果四舍五入到最接近。

假设我们必须将 3 个数字相加9.9 + 8.4 + 1.4
确切的结果是19.7，但我们只有两位数，所以应该四舍五入为20.

如果我们首先对9.9 + 8.4 求和，我们得到18.3，然后四舍五入为18.
然后我们将18. + 1.4 相加，得到19.4 四舍五入为19.。

如果我们首先将最后两项相加8.4 + 1.4，我们得到9.8，现在还不需要四舍五入。
然后9.9 + 9.8 我们得到19.7 舍入到20.，得到不同的结果。

(9.9 + 8.4) + 1.4 与 9.9 + (8.4 + 1.4) 不同，求和运算不是关联的，这是由于中间舍入造成的。我们也可以用其他舍入模式展示类似的例子......

问题与 53 位有效数的基数 2 完全相同：中间舍入将导致非关联性，无论基数或有效数长度如何。

要消除此问题，您可以对数字进行排序以使顺序始终相同，或者消除中间舍入并仅保留最后一个，例如使用像这样的超级累加器https://arxiv.org/pdf/1505.05571.pdf
...或者只是接受一个近似结果（由您分析平均或更严重的错误并决定是否可以接受...）。

【讨论】：