The basic process at work here is the following

1. You've got two numbers of widely differing scale:

   printf "%.17e\n", 1.00000000000000040e+000;;
   1.00000000000000040e+000
   
   printf "%.17e\n", 5.00000000000100030e-016;;
   5.00000000000100030e-016
   [download]

2. If you just add those together, you obviously loose some precision from the smaller number due to hardware limitation:

   printf "%.17e\n", $S = 1.00000000000000040e+000 + 5.00000000000100030e
   +-016;;
   1.00000000000000090e+000
   [download]

3. So, doubledouble uses a second double to store that lost precision. But to do so,it needs to 'recover' the loss.

   By subtracting the two numbers in turn, from their sum (numerically ought to be zero); but due to the machine limitations, a delta falls out and the lost precision is recovered:

   printf "%.17e\n", $S = 1.00000000000000040e+000 + 5.00000000000100030e
   +-016;;
   1.00000000000000090e+000
   
   printf "%.17e\n", $e = ( $S = 1.00000000000000040e+000 + 5.00000000000
   +100030e-016 ) - 1.00000000000000040e+000;;
   4.44089209850062620e-016
   
   printf "%.17e\n", ( $e = ( $S = 1.00000000000000040e+000 + 5.000000000
   +00100030e-016 ) - 1.00000000000000040e+000 ) - 5.00000000000100030e-0
   +16;;
   -5.59107901500374110e-017
   [download]

   That last value is the lost precision that needs to be stored in the second double.

Of course, it doesn't end there. There might already be values in those other (low) doubles; and they need to be added together along with the spillage from the hi order doubles above.

But that calculation itself can result in what might be termed 'overflow' or 'carry-over'; and that needs to be added back into the high order part of the result. But that ...

Hopefully, you get the picture.

The sequences of additions and subtractions in that sub are meant to ensure that borrows from the high order doubles and carries from the low order doubles are sorted out and merged; with the result that you end up with 105/6 bits of precision.

It -- the C++ -- appears to work; plenty of people have used it -- but I want to port (parts of) it to C; and that's caused me to look closely at it. And there is weirdness afoot.

According to this and (for me) this, that should either not happen; or at least be controllable.

But then, that may not have been the case for which compiler the code was originally (or subsequently) targeted at; so it's a good thought. Thank you.

---

With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'

Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.

"Science is about questioning the status quo. Questioning authority".

In the absence of evidence, opinion is indistinguishable from prejudice.

I'm with torvalds on this Agile (and TDD) debunked I told'em LLVM was the way to go. But did they listen!

I took a look through some of the source code you linked, and the type cast seems to be applied in particular when there's a subtraction that's obviously close to zero. I don't have enough C++ experience to know how the compilers are expected to behave on such things (even more so older versions), but I walked through a bunch of the arithmetic and there's a lot of adding and subtracting things that should either add up to zero or add up neatly to just one of the terms that was in a previous step, so I'm actually a little surprised that there aren't more oddball things in there (or at least a note about setting a compiler directive not to simplify) to keep the epsilons from getting optimized away.

There are enough things that just cancel away in odd ways that I was partly looking to verify you hadn't made any typos!

EDIT: in your first link it actually points out that the Intel C++ compiler does allow associativity transformations by default, so it may be a bit of legacy code to prevent that.