FORTRAN floating point accuracy problems

[Dieter]

(04-04-2016 06:42 PM)Dieter Wrote:  

(04-04-2016 09:02 AM)Werner Wrote:  The following formulas, courtesy of John H. Meyers on comp.sys.hp48 do not require any fiddling at all:
...
T := HMS + FP(HMS*100)/150;
HR := T + FP(T)/1.5;

(...)
So a Fortran version of this sounds like a good idea. I would give it a try. Maybe (or even "probably") this will solve the whole issue very elegantly.

I now took a closer look at this, at here is what I found:

1. In most cases the code works surprisingly well. The key seems to be the division by 150 that shifts roundoff errors in the number representation by two digits to the right. Here they usually will not seriously affect the final result.

2. However, there are cases where errors may occur. I did some tests in Excel VBA. Here is an example:

Consider the case dms = 1/5,000 00000 00001 = 0,19999 99999 99999 6000... which is displayed as 0,2 since Excel does not show more than 15 significant digits. Try 10*this – 1 and you see the correct value 0,99999 99999 99996. In the above code, the term 100*dms rounds to 20 even with "all" 15 digits displayed. Now INT of this yields 20 exactly (not 19!) and FP(100*dms) = 100*dms – INT(dms) returns –0,00000 00000 00039, i.e. a negative value! This messes up the whole calculation and the final result is 0,34444... hours instead of 0,33333...

The key problem in the above example is the INT function that obviously returns the integer portion of what you see with 15 digits. This might avoid confusions for many users ("why is int(20) = 19?"), but numerically everything is going the wrong way from here.

Well... "wrong" is a hard word. What is the conversion code supposed to do? Should it take 0,19999 99999 99999 6 as 19 minutes and 99,9999... seconds? Then the "wrong" result is correct. Or should it assume that this is supposed to be 0,2 exactly since the 16th digit can be off anytime? Then the mentioned problem occurs.

So it seems the only reliable way is rounding the input value to an integer with as many digits as the working precision reliably allows. For instance 15 digits in IEEE754 double precision. Or at least round the fractional part of the input value to the maximum number of decimals that is reliably supported. So an input value like 30,2 would be rounded to 13 decimals (15 minus 2 for the integer part). Code for this approach has already been shown in this thread.

Dieter