Error Propagation upon Summing Single-Precision (float) Values

Question

I'm learning single precision and would like to understand the error propagation. According to this nice website, addition is a dangerous operation.

So I wrote a small C program to test how quickly the errors add up. I'm not entirely sure if this is a valid way of testing. If it is, I'm unsure how to interpret the result, see below.

#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;
}

The result of executing

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

is

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

If the test is valid

• Why does the error change sign? If 0.1 is represented as e.g. 0.100000001, shouldn't this accumulate always to the same bias, irrespective of the number of summations?
• What's special about 168600 summations (see *)? The error becomes very small. Might be a coincidence.
• Which wall is being hit at incVal = 32.00 (see **, last two lines). I'm still well below the unsigned long limit.

Thanks in advance for your effort.

Answer 1

First, it's important to know that 0.1 can't be represented exactly, in binary it's has periodically repeating digits. The value would be 0.0001100110011.... Compare to how 1/3 and 1/7 are represented with decimal digits. It's worth repeating your test with increment 0.25, which can be represented exactly as 0.01.

I'll illustrate the errors in decimal, that's what we humans are used to. Let's work with decimal, and assume we can have 4 digits of precision. Those are the things happening here.

• Division: let's calculate 1/11:

  1/11 equals 0.090909..., which is probably rounded to 0.09091. This is, as expected, correct to 4 significant digits (in bold).

• magnitude differences: suppose we calculate 10 + 1/11.

  When adding 1/11 to 10, we have to do more rounding, since 10.09091 are 7 significant digits, and we have only four. We have to round 1/11 to two digits after the point, and the calculated sum is 10.09. That's a underestimation. Note how only one significant digit of 1/11 is retained. If you add a lot of small values together, this will limit the precision of your final result.

• Now calculate 100 + 1/11. Now we round 1/11 to 0.1 and represent the sum as 100.1. Now we have a slight overestimation instead of a slight underestimation.

  My guess is the pattern of sign changes in your test are the effect of systematic slight underestimation vs. overestimation depending on the magnitude of base.

• What about 1000 + 1/11? Now we can't have any digits after the point, as we have 4 significant digits before the point already. 1/11 is now rounded to 0, and the sum is still 1000. That's the wall you're seeing.

• Another important thing you're not seeing in your test is: what happens if the two values have a different sign. Calculate 1.234 – 1.243: both numbers have 4 significant digits. The result is -0.009. Now the result has only one correct significant digit instead of four.

An answer to a similar question here: How does floating point error propagate when doing mathematical operations in C++? . It has a few links to more information.

Answer 2

To answer your questions...

1 - IEEE float rounds to even mantissas. This was done specifically in order to prevent error accumulation from always biasing in one way or the other; if it always rounded down, or rounded up, your errors would be much larger.

2 - Nothing in particular is special about 168600 in and of itself. I haven't mathed it out but it's entirely likely that it ends up making a cleaner value in binary representation (i.e. a rational/non-repeating value). Look at the values in binary, not decimal, and see if that theory holds up.

3 - The limiting factor might be due to the float mantissa being 23 bits long. Once base gets to be a certain size, increment is so small in comparison to base that computing base + increment and then rounding the mantissa back down to 23 bits completely erases the change. That is, the difference between base and base + increment is rounding error.

Answer 3

The "wall" you are hitting has nothing to do with the increment value, if it is constant through the addition and you start at zero. It has to with the iters. 2^23 = 8 million, and you are doing 86 million additions. So once the accumulator is 2^23 bigger than the increment, you hit the wall.

Try running the code with 86323200 iterations, but an increment of 1 or 0.0000152587890625 (or any power of 2). It should have the same relative problem as an increment of 32.