adding numbers and floating point errors

[smeenz]

Hi,

Today I discovered that floating point errors seem to be much more of a problem than I had anticipated.

Specifically, it seems that if I set up a simple loop to add 0.05 to a variable and then display the result, it quickly goes wrong.

For example, the following code:

  my $c = 0;
  for ($x = 0; $x <= 1; $x += 0.05){
    $c++;
    print "[$c] x= $x (".sprintf("%20.40f",$x).") \n";

    }

produces this:

[1] x= 0 (0.0000000000000000000000000000000000000000)
[2] x= 0.05 (0.0500000000000000027755575615628913510591)
[3] x= 0.1 (0.1000000000000000055511151231257827021182)
[4] x= 0.15 (0.1500000000000000222044604925031308084726)
[5] x= 0.2 (0.2000000000000000111022302462515654042363)
[6] x= 0.25 (0.2500000000000000000000000000000000000000)
[7] x= 0.3 (0.2999999999999999888977697537484345957637)
[8] x= 0.35 (0.3499999999999999777955395074968691915274)
[9] x= 0.4 (0.3999999999999999666933092612453037872910)
[10] x= 0.45 (0.4499999999999999555910790149937383830547)
[11] x= 0.5 (0.4999999999999999444888487687421729788184)
[12] x= 0.55 (0.5499999999999999333866185224906075745821)
[13] x= 0.6 (0.5999999999999999777955395074968691915274)
[14] x= 0.65 (0.6500000000000000222044604925031308084726)
[15] x= 0.7 (0.7000000000000000666133814775093924254179)
[16] x= 0.75 (0.7500000000000001110223024625156540423632)
[17] x= 0.8 (0.8000000000000001554312234475219156593084)
[18] x= 0.85 (0.8500000000000001998401444325281772762537)
[19] x= 0.9 (0.9000000000000002442490654175344388931990)
[20] x= 0.95 (0.9500000000000002886579864025407005101442)

which is... okayish.. because the errors only show up when I force sprintf to display to 40dp... the problem is that even though I'm only adding 0.05 (and not 0.00000000000000000000000005), the error still creeps into the results pretty quickly.

If I change the loop to run from 0 to 4, which is 100 interations, perl tells me at iteration 74 that 3.65 + 0.05 = 3.69999999999999.

I'm not sure what to think of this.. I wouldn't have expected that adding 0.05 to a number 74 times would start to produce errors as soon as the 1st decimal place.

[imp]

BrowserUK provided an excellent explanation for why this occurs a few years ago.

To (mostly) fix the issue for your script you can use bignum.

use bignum;

my $c = 0;
for ($x = 0; $x <= 4; $x += 0.05){
    $c++;
    print "[$c] x= $x (".sprintf("%20.40f",$x).") \n";
}

Read more... See all output (5 kB)

[78] x= 3.85 (3.8500000000000000888178419700125232338905)
[79] x= 3.9 (3.8999999999999999111821580299874767661095)
[80] x= 3.95 (3.9500000000000001776356839400250464677811)
[81] x= 4 (4.0000000000000000000000000000000000000000)

[smeenz]

Thanks. That explanation does indeed explain the issue.

[blazar]

perldoc -q decimals

[smeenz]

Thanks.. I had tried -q float and -q precision.. hadn't thought of -q decimal.

[philcrow]

Unless your interval is a fraction of a power of two or can be represented by a small number of them added together, your number can't be represented exactly in binary. Think about what you are asking your poor binary computer to do. .05 is 1/20. The first number that isn't bigger than your number is 1/32. Subtract that from 1/20. Now look for 1/2**n which is smaller than what's left. Continue until it comes out even or you run out of precision.

Decimal fractions are rarely representable exactly in computing memory.

To prevent errors from building up, round more often.

As for

3.65 + 0.05 = 3.69999999999999

I wouldn't say the error is in the first place. It's actually in the last place you showed. Any prior rounding would work. If the difference was significant in the first digit, you could round in the second place without getting the correct answer.

Decimals tie people in knots even when computers aren't involved. I remember a heated discussion in a math class I took, where the teacher was making no headway in convincing some student that .3 repeating is 1/3.

Phil

[smeenz]

Thanks to everyone.. I see what the problem is now... I had thought that these sorts of numbers would be represented as whole numbers and then the decimal place added back in for the final result, so 0.0005 + 0.0005 would be the same as adding 5+5 (or 00000101 and 00000101), and then shifting the decimal point appropriately at the output stage, but now that you've explained that the computer is trying to approximate the actual fraction using powers of 2, it makes a lot more sense as to what is happening.

[ikegami]

If you really have a loop and you want to avoid the accumulation of error, you can avoid use an integer as the loop counter, and compute the current value of the real from the integer.

# Loop over [0.00..1.00] by 0.05
for (my $i = 0; $i <= 100; $i += 5) {
   my $x = $i / 100;
   printf("[%d] x= %g (%20.40f)\n", ++$c, $x, $x);
}

# Loop over [0.00..1.00] by 0.05
for (my $i = 0; $i <= 20; $i++) {
   my $x = $i / 20;
   printf("[%d] x= %g (%20.40f)\n", $i+1, $x, $x);
}

# Loop over [0.00..1.00] by 0.05
for (0..20) {
   my $x = $_ / 20;
   printf("[%d] x= %g (%20.40f)\n", $_+1, $x, $x);
}

All of the above output the following:

[1] x= 0 (0.0000000000000000000000000000000000000000)
[2] x= 0.05 (0.0500000000000000030000000000000000000000)
[3] x= 0.1 (0.1000000000000000100000000000000000000000)
[4] x= 0.15 (0.1499999999999999900000000000000000000000)
[5] x= 0.2 (0.2000000000000000100000000000000000000000)
[6] x= 0.25 (0.2500000000000000000000000000000000000000)
[7] x= 0.3 (0.2999999999999999900000000000000000000000)
[8] x= 0.35 (0.3499999999999999800000000000000000000000)
[9] x= 0.4 (0.4000000000000000200000000000000000000000)
[10] x= 0.45 (0.4500000000000000100000000000000000000000)
[11] x= 0.5 (0.5000000000000000000000000000000000000000)
[12] x= 0.55 (0.5500000000000000400000000000000000000000)
[13] x= 0.6 (0.5999999999999999800000000000000000000000)
[14] x= 0.65 (0.6500000000000000200000000000000000000000)
[15] x= 0.7 (0.6999999999999999600000000000000000000000)
[16] x= 0.75 (0.7500000000000000000000000000000000000000)
[17] x= 0.8 (0.8000000000000000400000000000000000000000)
[18] x= 0.85 (0.8499999999999999800000000000000000000000)
[19] x= 0.9 (0.9000000000000000200000000000000000000000)
[20] x= 0.95 (0.9499999999999999600000000000000000000000)
[21] x= 1 (1.0000000000000000000000000000000000000000)

[spiritway]

Computers don't handle fractional floating point values very well. Anything other than powers of 1/2 (or their multiples) results in a form that repeats endlessly. This results in truncation or rounding errors, and these will tend to accumulate when you iterate.

You might want to consider keeping your data in rational form (e.g., bigrat) while you are performing your calculations. Doing this can save you lots of problems with rounding or truncation errors. Once you've finished your calculations you can convert your rationals into floating point format, if you wish, for display.

[ikegami]

Anything other than powers of 1/2 (or their multiples) results in a form that repeats endlessly.

Not quite. Anything other than sums of powers of 2 result in a form that repeats endlessly.

Sums of powers of 2 do not repeat endlessly. For example, 0.3125 is not a power of 2, yet it does not repeat because it is a sum of powers of 2 (2^-2 + 2^-4).

Note that "powers of 1/2" and "powers of 2" are synonymous.

Update: Nevermind. What you said and what I said are mathematically equivalent. 0.3125 is a multiple of a power of 2 (5 * 2^-4).

[traveler]

You might try Math::FixedPrecision if this is for a real application.

[mattr]

Well, I just saw the exact same question posted somewhere else a couple days ago so it is interesting you came up with the same exact question. Anyway I found some more info about it (see end of my post).

If you just print it without using sprintf, there is no problem. If you do use sprintf you are intentionally trying to see how the internal math library is working. IIRC 13 significant digits should be fine though your example looks like 15 digits. Maybe that's a 16 bit precision C float. It happens due to the math routines and IIRC that floating point numbers are generally represented in scientific notation, which means it only guarantees a fixed number of digits precision. Also shortcuts are made for speed sometimes. If you needed as many digits as you are requesting you should be using one of the modules that does that for you as other people have noted for arbitrary precision or more precision (I remember BigNum myself), PDL or others might help too.

If you had a problem showing up without using sprintf I'd be a little more worried. Conceivably a 64bit computer might be using higher precision routines in their base libraries, but I could not say myself.

Bottom line is, this is a well known artifact common to computers and is not a bug in Perl.

HOWEVER, You want to be really careful about numbers the computer spits out when you are doing anything important like financial, scientific or statistical calculations. A rounding error or something like this can get ugly. Some interesting info along here.. also I am not sure about whether there is a perl facility to (looks like Math::BigFloat would) provide access to 80 bit extended precision values.

Wikipedia on floating point addition
General Decimal Arithmetic website at IBM and the FAQ , especially the part on how many digits precision are needed for decimal arithmetic. (It says business requirements are typically 20-30 digits but interest rate calculations could require over 2000 digits.)
What Every Computer Scientist Should Know About Floating-Point Arithmetic (html) and pdf (neat.)
perlnumber (the pod on perl's numerics)

UPDATED:Boneheaded mistake made unusable links! Sorry. Now fixed.

[Anonymous Monk]

I htink this looks horrible, but is not :-)
You noticed the first time the error gets less than zero! If you calculate the error, you'll see it's around 10E-15 at this point (at last on my pc).

#!/usr/bin/perl
my $x=0.0;
my $y=0;
my $e_min=1.0;
my $e_max=0.0;
my $n=0;
my $e;
for ( 1 .. 999 ) {
    $x += 0.05;
    $y += 5;
    $e = abs($x - $y / 100);
    $n++ if $e == 0.0;
    $e_min = $e if $e > 0.0 && $e < $e_min;
    $e_max = $e if $e > $e_max;

    printf "%4d x = %.30f e = %E\n", $_, $x, $e;
}

printf "min %E max %E exact %d\n", $e_min, $e_max, $n;

result

 998 x = 49.899999999999302247033483581617 e = 6.963319E-13
 999 x = 49.949999999999299404862540541217 e = 7.034373E-13
min 2.775558E-17 max 7.034373E-13 exact 16

You'll find the lines arround 190 interesting:

 189 x = 9.449999999999999289457264239900 e = 0.000000E+00
 190 x = 9.500000000000000000000000000000 e = 0.000000E+00
 191 x = 9.550000000000000710542735760100 e = 0.000000E+00