Tracing a rounding error 3

[WebDrake]

I'm wondering if anyone can give any advice on tracing a rounding error in a program.

I have two slightly different implementations of a process, which basically come down to the following.

In Version 1, various increments are added to a [tt]double[/tt] variable, and at the end I check to see if it is greater or equal than a threshold variable.

In Version 2, the increments are summed separately and the total is then added to the variable. The variable is then compared to the threshold as normal.

In some cases the values produced are very slightly different (on the order of 2.220446e-16). I've been trying to trace exactly why but have not been able to come up with a clear-cut example of what it is that makes the difference.

Any suggestions/comments?

Thanks in advance!

 

[cpjust]

It's strange that there's no FAQ on this, since the question comes up so often... But anyways, floating point numbers are never going to be 100% accurate because of the way they're stored:

http://en.wikipedia.org/wiki/Floating_point

 

[WebDrake]

Sure. I know the principle, I was just wondering if there was any advice on how to trace exactly where such an error comes in.

My impression was that the error was occurring in the event of something like this happening:

    x += a;
    x -= a;

Where [tt]a[/tt] is of much greater order than [tt]x[/tt] (e.g. [tt]a[/tt] ~ 0.1, [tt]x[/tt] ~ 10^-20).

Of course, if you sum those increments first you get [tt]x += 0[/tt] so [tt]x[/tt] is unaffected, but doing the above will wipe out the fine detail.

I tried to put some statements in the code for it to alert me when that happened but had no luck...

 

[DaveSinkula]

http://en.wikipedia.org/wiki/Floating_point#Peculiar_properties_of_floating_point_arithmetic

 

[ArkM]

It's not an error. It's a normal floating point calculations behaviour.
That's why math programming is not so easy.
To catch the moment:

for (;;)
{
  // 1st add thread
  // 2nd add thread
  // Compare temp results with ==
}

Don't compare floats with == or != (except == 0) in production codes.
Remember: 53 bits of double mantissa, right shift of added item (where is its right bits?). Try to simulate this situation on 2-3 digit precision arithmetics (take a paper sheet and a pencil.

 

[WebDrake]

Belated thanks for the useful replies. And yes, I know it's normal floating point behaviour. I'm just interested in tracing when numerical methods diverge.

 

[lionelhill]

it sort of depends on your numbers, but in a sense it's fairly obvious where the difference lies.

Version 2 is usually better. Think about a 2-place decimal number with floating point. Say it holds the number 80.

Now make your increment; say the increment is 0.1

80 + 0.1 = 80.1, but oops! That needs 3 digits to store. So the result will actually be 80 in our 2-digit register... no change.

Do this ten times for ten increments, and the answer is still 80.

Now try Version 2; start with zero, and add 0.1, the first increment to sum. The answer is 0.1. Now add the second, and so on until you've done all ten. The answer is 1.0; we haven't had any problems with space in the register.
Now add this to the 8.0 we needed to start with, and the answer is 81, which fits in our 2-place register and therefore works fine.

What's happening in version 1 is that each time you add an increment you're having to slide the increment number to the right in order to line up its digits with those of the main sum number, and digits are dropping off the right hand side. In version 2 there's a lot less sliding to sum the increments, especially if all the increments are fairly similar in size. Then the result is slid less far, and only once, to line up with the double variable. There is therefore only one error, and it's smaller.

Have a look in any decent textbook on numerical analysis for further details.

And if anyone is reading this and thinking "I'm safe, because I always use doubles", stop right there! The best defence is to use good, robust algorithms, and you're on thin ice if you think you can save yourself from numerical catastrophe by racking up the precision and hoping...