c - 对单精度(浮点)值求和时的错误传播

Question

我正在学习单精度并想了解错误传播。根据this nice website ，加法是一个危险的操作。

所以我编写了一个小的 C 程序来测试错误累积的速度。我不完全确定这是否是一种有效的测试方法。如果是，我不确定如何解释结果，请参见下文。

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

#define TYPE float
#define NUM_IT 168600

void increment (TYPE base, const TYPE increment, const unsigned long num_iters) {

  TYPE err;
  unsigned long i;
  const TYPE ref = base + increment * num_iters;

  for (i=0; i < num_iters; i++ ) {
    base += increment; 
  }
  err = (base - ref)/ref;
  printf("%lu\t%9f\t%9f\t%+1.9f\n", i, base, ref, err);

}

int
main()
{
  int j;
  printf("iters\tincVal\trefVal\trelErr\n");

  for (j = 1; j < 20; j++ ) {
    increment(1e-1, 1e-6, (unsigned long) (pow(2, (j-10))* NUM_IT));
  }

  return 0;
}

执行结果

gcc -pedantic -Wall -Wextra -Werror -lm errorPropagation.c && ./a.out  | tee float.dat  | column -t

是

iters     incVal     refVal     relErr
329       0.100328   0.100329   -0.000005347
658       0.100657   0.100658   -0.000010585
1317      0.101315   0.101317   -0.000021105
2634      0.102630   0.102634   -0.000041596
5268      0.105259   0.105268   -0.000081182
10537     0.110520   0.110537   -0.000154624
21075     0.121041   0.121075   -0.000282393
42150     0.142082   0.142150   -0.000480946
84300     0.184163   0.184300   -0.000741986
168600    0.268600   0.268600   +0.000000222    <-- *
337200    0.439439   0.437200   +0.005120996
674400    0.781117   0.774400   +0.008673230
1348800   1.437150   1.448800   -0.008041115
2697600   2.723466   2.797600   -0.026499098
5395200   5.296098   5.495200   -0.036231972
10790400  10.441361  10.890400  -0.041232508
21580800  25.463778  21.680799  +0.174485177
43161600  32.000000  43.261597  -0.260313928    <-- **
86323200  32.000000  86.423195  -0.629729033

如果测试有效

• 为什么错误会改变符号？如果 0.1 表示为例如0.100000001，无论求和次数如何，这不应该总是累积到相同的偏差吗？
• 168600 求和有什么特别之处(参见 *)？误差变得很小。可能是巧合。
• 在 incVal = 32.00 时撞到了哪堵墙(参见 **，最后两行)。我仍然远低于 unsigned long 限制。

预先感谢您的努力。

Answer 1

首先，重要的是要知道 0.1 不能准确表示，在二进制中它有周期性重复的数字。该值将为 0.0001100110011...。比较 1/3 和 1/7 是如何用十进制数字表示的。值得用增量 0.25 重复您的测试，它可以精确表示为 0.01。

我将用十进制说明错误，这是我们人类习惯的做法。让我们使用小数，并假设我们可以有 4 位精度。这些就是这里发生的事情。

• 除法:让我们计算 1/11:

  1/11 等于 0.090909...，可能四舍五入为 0.09091。正如预期的那样，这正确到 4 位有效数字(粗体)。

• 幅度差异:假设我们计算 10 + 1/11。

  将 1/11 加到 10 时，我们必须进行更多舍入，因为 10.09091 是 7 位有效数字，而我们只有四位。我们要将 1/11 四舍五入到小数点后的两位数，计算出的总和为 10.09。这是低估了。请注意如何仅保留 1/11 的一位有效数字。如果将很多小值加在一起，这将限制最终结果的精度。

• 现在计算 100 + 1/11。现在我们将 1/11 舍入为 0.1，并将总和表示为 100.1。现在我们有轻微的高估而不是轻微的低估。

  我的猜测是，您测试中符号变化的模式是系统性轻微低估与高估的影响，具体取决于 base 的大小。

• 1000 + 1/11 呢？现在我们不能在点之后有任何数字，因为我们在点之前已经有 4 个有效数字。 1/11 现在四舍五入为 0，总和仍然是 1000。这就是您看到的墙。

• 您在测试中没有看到的另一件重要事情是:如果两个值的符号不同会发生什么。计算 1.234 – 1.243:两个数字都有 4 位有效数字。结果是 -0.009。现在结果只有一位正确的有效数字，而不是四位。

此处类似问题的答案:How does floating point error propagate when doing mathematical operations in C++? .它有一些指向更多信息的链接。