单精度(浮点)值求和时的误差传播

Question

我正在学习单精度，并想了解错误传播。根据这个不错的网站的说法，加法是一种危险的操作。

所以我编写了一个小的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个重要数字，而我们只有4个。在点之后，我们必须把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。现在的结果只有一个正确的有效数字，而不是四个。

这里有一个类似问题的答案：在C++中做数学运算时浮点误差是如何传播的？。它有一些链接到更多的信息。

Answer 2

回答你的问题..。

IEEE浮动轮到甚至mantissas。这样做是为了防止错误积累总是以一种或另一种方式出现偏差；如果它总是四舍五入，那么您的错误就会大得多。

没有什么是特别的，它本身就是168600。我还没有弄清楚它，但很可能它最终会在二进制表示中产生一个更清晰的值(即，一个rational/非重复的值)。看看二进制值，而不是十进制，看看这个理论是否成立。

3-限制因素可能是由于浮子尾数为23位长。一旦base达到一定的大小，与base相比，increment是如此之小，以至于计算base + increment，然后将尾数舍入到23位，完全消除了变化。也就是说，base和base + increment之间的区别是舍入误差。

Answer 3

你所打的“墙”与增量值无关，如果它是通过加法常量，并且从零开始。这与iters有关。2^23 =800万，而您正在进行8600万个加法。因此，一旦累加器比增量大2^23，就会碰到墙。

尝试使用86323200次迭代来运行代码，但增量为1或0.0000152587890625 (或任何2的幂)。它应该有与32的增量相同的相对问题。