FORTRAN floating point accuracy problems

[Dieter]

(03-27-2016 03:02 PM)HP67 Wrote:  I'm getting some unacceptable results from rounding in old FORTRAN.
(...)
You can see in test case #3 above there is a big problem extracting the minutes. The code winds up with 29 instead of 30 because of the truncation function.

Sure. This is caused by the inexact binary representation of decimal numbers. For instance, the value "30.3" here obviously has a binary value that translates to 30.29999999999... in decimal notation. This again is interpreted as 30 degrees, 29 minutes and 99.99999... seconds, to that the final result is 30.5111... degrees.

I translated your Fortran code to Excel VBA, and indeed the same problems occured. There are several ways to overcome the problem. For a quick and dirty solution, simply add a tiny amount to the ss value in the first line of code, say, 1E-12 if 15 digits are used and the angle is below 1000° . This will cause a slight error (1E–10 seconds or 3 E–15 degrees) in the final result.

A better solution may check if the mm.ss value is very close to the next higher integer, i.e. mm.999.... Since the seconds cannot exceed 59.9999, even a test whether the ss portion of mm.ss exceeds 0.6 may be sufficient. At least if the user does not enter non-standard values like 30°40'85".

Here are two solutions in VBA. The first one checks if the ss part of mm.ss exceeds 0.99999... In this case the minutes are rounded up to the next higher integer, which also means that the seconds are zero.

Code:

Function dms2dd(dms)

dd = Int(dms)
mmss = (dms - dd) * 100
mm = Int(mmss)

If mmss - mm > 0.9999999999 Then
   mm = mm + 1
   ss = 0
Else
   ss = (mmss - mm) * 100
End If

dms2dd = dd + mm / 60 + ss / 3600

End Function

The second solution simply increases mm by a tiny amount which is close to one ULP of the input. For 52-bit binary arithmetics you may try 2^–40 or 1E–12.

Code:

Function dms2dd(dms)

const tiny = 2^(-40)

dd = Int(dms)
mmss = (dms - dd) * 100 + tiny
mm = Int(mmss)
ss = (mmss - mm) * 100

dms2dd = dd + mm / 60 + ss / 3600

End Function

A third solution would simply round the argument (i.e. dms) to the max. possible number of decimals. For instance, if dms=30.3000 and you're working with 15 digit precision, the 0.3000 part cannot have more than 13 digits following the decimal point (15 minus two for the integer part). So the fractional part cannot have 15 digits like 0.29999 99999 99997 (which causes the roundoff problem), but it has to be 0.30000 00000 000. In the final code this rounding procedure would replace the line with the "tiny" variable in the second solution:

Code:

Function dms2dd(dms)

dd = Int(dms)
mmss = 100 * Round(dms - dd, 13)
mm = Int(mmss)
ss = (mmss - mm) * 100

dms2dd = dd + mm / 60 + ss / 3600

End Function

Maybe this is the preferred solution.

BTW, your result table seems to be wrong. For instance, for dms(71.0030) the expected value is 71 + 30/3600 = 71.00833333... and not 71.00830. This is also true for several other results.

Dieter