FORTRAN floating point accuracy problems

[Dieter]

(03-28-2016 09:29 AM)HP67 Wrote:  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?

The first line (after the declarations) of your code, the one that reads SS = (DMS - IDINT(DMS)) * 100.0D0.

(03-28-2016 09:29 AM)HP67 Wrote:  Yes, the numbers are generated internally so they should be well-formed.

That's good to know. However, the third and preferred solution should handle either case.

(03-28-2016 09:29 AM)HP67 Wrote:  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?

Let's assume the working precision is 15 decimal digits (like IEEE754 double precision). If 2-3 decimal digits are kept for the integer portion of the argument (whole degrees), the remaining 12-13 hold the mmss part. So the last digit of mmss is (roughly) 1E-12. Or (base 2) about 2^–39. The idea is to get this value large enough to do the rounding but small enough to minimize its mathematical impact, as it affects the last digit(s) of the final result. If you want to implement this idea I would recommend a value for "tiny" that is based on the argument dms, e.g. tiny=dms/1E14 for 15-digit working precision.

So this is how it is supposed to work:

Code:

entered argument dms:
exact:  30,2000000000000
stored: 30,199999999999997...

fractional part of dms:
exact:   0,2000000000000
stored:  0,199999999999997...
                        ^^
now add a tiny amount here
so that the stored value rounds up
while the 13 digits of the exact value 0,2
are not affected, e.g. tiny=1E-14

exact:   0,2000000000000
stored:  0,199999999999997...
+ tiny:  0,00000000000001
rounded: 0,200000000000007

But the third approach should work even better:

(03-28-2016 09:29 AM)HP67 Wrote:  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.

Sure. Write your own rounding function. ;-)

CAVE: I DO NOT HAVE ANY EXPERIENCE WITH FORTRAN! But from what I understand there is a NINT resp. IDNINT function that rounds a value to its closest integer. For a function that rounds to n digits, simply multiply the argument by 10^n, round this to the next closest integer and divide by 10^n again. Here is the general idea for rounding a positive value X to D decimal places:

Code:

THIS IS NOT FORTRAN CODE, IT JUST SHOWS THE GENERAL IDEA!

R = IDNINT(10.0D0**D)
X = IDNINT(R*X) / R

If IDNINT is not available and IDINT is all you got:

R = IDINT(10.0D0**D + 0.5D0)
X = IDINT(R*X + 0.5D0) / R

Please note that even R is rounded since we cannot be sure that calculating 10^13 will return exactly 10,000,000,000,000. So just to be sure...

Once you know the number of digits to round (13 in the following example), rounding is easy. Simply define the factor R somewhere before (or use a constant - don't know how this is done in Fortran) and the rest is trivial:

Code:

R = 1.0D13
MMSS = 100D0 * (IDNINT(R * (DMS - DD)) / R)

or with a simple INT

R = 1.0D13
MMSS = 100D0 * (IDINT(R * (DMS - DD) + 0.5D0) / R)

A closer-to-perfect solution would round to a variable number of digits, depending on the magnitude of the dms argument. But 13 digits works well in my Excel VBA implementation, at least for angles up to 999 degrees.

Now try this and post your improved code here.

Dieter