FORTRAN floating point accuracy problems

[HP67]

Hi Dieter, thanks so much for your time and explanations. I'm not sure I got it all. I understand floating point representation is inexact. But I don't know how to fix it.

I did see in some cases adding a small correction would help but when I added it in all cases it made the results much worse overall. Below when you said to add the correction in the first line, did you mean the first line of my code (I think there it is already too late) or in the first line of your example code?

(03-27-2016 05:15 PM)Dieter Wrote:  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.

(03-27-2016 05:15 PM)Dieter Wrote:  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".

Yes, the numbers are generated internally so they should be well-formed. I think I understand your second suggestion but I'll have to look at your code and my code further.

(03-27-2016 05:15 PM)Dieter Wrote:  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.

I will have to check. The computer I'm using does not use IEEE binary, it uses a proprietary hexadecimal format and I don't remember how many bits of precision are carried in the mantissa. I'll get back to you on that. How do you choose a correction value close to one ULP? You said for 52 bit binary we should use 2^-40 or 1E-12. What is the connection between these values and the bit precision?

(03-27-2016 05:15 PM)Dieter Wrote:  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:


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

Maybe this is the preferred solution.

I don't have a rounding function and I don't understand how to write one. Is this solution still possible? Is it possible to write a general purpose rounding function? All I have is trig functions, integer (which does not simply take the integer portion and store in the 8 byte floating point value but rather it demotes the 8 byte floating point integer value to a 4 byte integer type.

(03-27-2016 05:15 PM)Dieter Wrote:  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.

Thanks, that will teach me to rely on test data from the web. I should have checked this with an HP calculator...

Thank you very much.