FORTRAN floating point accuracy problems

[HP67]

I'm getting some unacceptable results from rounding in old FORTRAN. I know there are some very talented mathematicians on the board here so I wonder if you guys would be willing to help me. To say I'm bad at math would be an understatement! I have no idea where things are going wrong.

I'm trying to accomplish something very simple in a language which has very little in the way of helpful functions to do what I want. Given a floating point input value for degrees, minutes, seconds in the form dd.mmsst... I want to produce the value in decimal degrees. This is no problem of math even for me. The problem is extracting the digits with only a truncation function and nothing else helpful in the repetoire as far as I can see.

Here is the piece of code:

Code:


      DOUBLE PRECISION FUNCTION DMS2DD (DMS)
      REAL*8 DMS, DD, MM, SS
      SS = (DMS - IDINT(DMS)) * 100.0D0
      SS = ((SS - IDINT(SS))  * 100.0D0) / 60.0D0
      MM = (IDINT((DMS - IDINT(DMS)) * 100.0D0) + SS) / 60.0D0
      DMS2DD = IDINT(DMS) + MM
      RETURN
      END

and here are the results of some tests:

Code:


  TEST    DMS         EXPECTED    RESULT       DIFFERENCE

     1    89.11150    89.18750    89.18750     0.00000
     2    12.15000    12.25000    12.25000     0.00000
     3    33.30000    33.50000    33.51111     0.01111
     4    71.00300    71.00830    71.00833     0.00003
     5    42.24530    42.41470    42.41472     0.00002
     6    38.42250    38.70690    38.70694     0.00004
     7    29.30300    29.50830    29.50833     0.00003
     8     0.49490     0.83030     0.83028     0.00002
     9    75.15000    75.25000    75.25000     0.00000
    10    43.22300    43.37500    43.37500     0.00000
    11     9.33450     9.56250     9.56250     0.00000
    12    33.57522    33.96450    33.96450     0.00000
    13    13.07244    13.12345    13.12345     0.00000
    14    21.30000    21.50000    21.50000     0.00000
    15    59.47211    59.78920    59.78920     0.00000
    16    65.11010    65.18360    65.18360     0.00000

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.

I find this kind of problem is coming up a lot in various calculations I'm doing and I don't know how to code around it.

Can anybody please explain how I can rearrange the computations or do something differently to avoid this problem?

Thanks.

[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

[Paul Dale]

For this kind of problem, I find first dealing with the sign and then getting into an integer domain as fast as possible and then applying corrections there is usually the best approach.

You've got integer degrees already, multiply by 60 and get integer minutes, then again by sixty to get seconds (which will have to stay real). Then sanity check everything to deal with the corner cases.

You'll end up with a pile of code along the lines of:

Code:

if seconds .ge. 60 then
  seconds = seconds - 60
  minutes = minutes + 1
endif

- Pauli

[Dieter]

(03-27-2016 10:47 PM)Paul Dale Wrote:  For this kind of problem, I find first dealing with the sign and then getting into an integer domain as fast as possible and then applying corrections there is usually the best approach.

OK, but how do you prevent the roundoff problems due to binary arithmetics? Consider the example with 33,3000. Isolating the minutes leads to 29,99999... instead of 30, which then is interpreted as 29 minutes and 99,99999... seconds, i.e. 30 minutes and 39,9999... seconds. So some kind of rounding has to be applied.

Dieter

[Paul Dale]

The most recent time I did this was using decimal arithmetic and it took several goes.

Yes, rounding will be required. At least it is easy enough to detect the errant case -- when the seconds value is greater than sixty there is possibly a problem. Figuring out where and how to round is another matter entirely -- this is getting close to fudging the answer.


- Pauli

[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.

[HP67]

Hi and thanks for helping with this.

(03-27-2016 10:47 PM)Paul Dale Wrote:  For this kind of problem, I find first dealing with the sign and then getting into an integer domain as fast as possible and then applying corrections there is usually the best approach.

That would be ideal but here the truncation and loss of precision occurs during extracting the minutes and seconds. And once in integer values the results of division also produce truncated values. So, I can't see how to do this.

(03-27-2016 10:47 PM)Paul Dale Wrote:  You've got integer degrees already, multiply by 60 and get integer minutes, then again by sixty to get seconds (which will have to stay real). Then sanity check everything to deal with the corner cases.

If I multiply by 60 to move the minutes to the left of the decimal there is already some loss of precision for certain values like 30 minutes for example. I don't know how to extract the minutes or seconds. This is exactly the problem.

[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

[HP67]

Thanks again, Dieter. Modern Fortran has many more functions including ceil, floor, and nearest integer but I'm using FORTRAN IV which has only truncation. I want to look over the hardware doc before I decide what to do. I don't understand the floating point instructions sine I never used them. I do remember there is some kind of hardware rounding function. It might be better to write whatever primitives I need in assembler and call them from FORTRAN. I'll post again after I have a better idea what hardware support there is.

What I should have mentioned: Modern Fortran has IEEE support and depending on the compiler that should include all IEEE rounding modes, etc. But, a long time ago (F90, possibly F77) Fortran added a feature called internal reads where you can "read" and "write" into variables without any I/O. Since FORTRAN always rounds on input/output this was a way to get a standard rounding mode for whatever precision you wanted. However it's not available as far back as I am working, so...

[HP67]

Dieter, something is wrong here.

The rounding function works to some extent when I specify a smallish value for decimal digits. If I round to 4 or 5 places it seems to work better than 13. I don't know why. In any case, the results when rounded are worse overall than without the rounding calculation. In many cases, it appears to be rounding down instead of up.

I think the problem is caused by the lack of a nearest integer function but I am not sure. The IDINT function is a double precision (8 byte) that simply returns the integer part like IP in RPL.

The MMSS value is for the following line (out of order because I added a WRITE statement to the subroutine to debug it. For example in the first test the input value is 89.11150 but MMSS is extracted and rounded to 13 digits according to your example above and comes out to 11.1499999649823)

Code:


  TEST    DMS       EXPECT       RESULT         DIFF

 MMSS   11.1499999649823     
     1    89.11150    89.18750    89.18750     0.00000
 MMSS   14.9999999441206     
     2    12.15000    12.25000    12.26111     0.01111
 MMSS   29.9999999348074     
     3    33.30000    33.50000    33.51111     0.01111
 MMSS  0.299999956041574     
     4    71.00300    71.00830    71.00833     0.00003
 MMSS   24.5299999602139     
     5    42.24530    42.41470    42.41472     0.00002
 MMSS   42.2499999403954     
     6    38.42250    38.70690    38.70694     0.00004
 MMSS   30.2999999374151     
     7    29.30300    29.50830    29.50833     0.00003
 MMSS   49.4899999350309     
     8     0.49490     0.83020     0.83028     0.00008
 MMSS   14.9999999441206     
     9    75.15000    75.25000    75.26111     0.01111
 MMSS   22.2999999765307     
    10    43.22300    43.37500    43.37500     0.00000
 MMSS   33.4499999415129     
    11     9.33450     9.56250     9.56250     0.00000
 MMSS   57.5219999533147     
    12    33.57522    33.96450    33.96450     0.00000
 MMSS   7.24419993348420     
    13    13.07244    13.12345    13.12345     0.00000
 MMSS   29.9999999348074     
    14    21.30000    21.50000    21.51111     0.01111
 MMSS   47.2111999522895     
    15    59.47211    59.78920    59.78920     0.00000
 MMSS   11.0095999669284     
    16    65.11010    65.18600    65.18360     0.00240

If I specify rounding to 4 places the results improve a lot but are still worse than when unrounded:

Code:


  TEST    DMS       EXPECT       RESULT         DIFF

 MMSS   11.1500000000000     
     1    89.11150    89.18750    89.18750     0.00000
 MMSS   15.0000000000000     
     2    12.15000    12.25000    12.25000     0.00000
 MMSS   30.0000000000000     
     3    33.30000    33.50000    33.50000     0.00000
 MMSS  0.300000000000000     
     4    71.00300    71.00830    71.00833     0.00003
 MMSS   24.5300000000000     
     5    42.24530    42.41470    42.41472     0.00002
 MMSS   42.2500000000000     
     6    38.42250    38.70690    38.70694     0.00004
 MMSS   30.3000000000000     
     7    29.30300    29.50830    29.50833     0.00003
 MMSS   49.4900000000000     
     8     0.49490     0.83020     0.83028     0.00008
 MMSS   15.0000000000000     
     9    75.15000    75.25000    75.25000     0.00000
 MMSS   22.3000000000000     
    10    43.22300    43.37500    43.37500     0.00000
 MMSS   33.4500000000000     
    11     9.33450     9.56250     9.56250     0.00000
 MMSS   57.5200000000000     
    12    33.57522    33.96450    33.96444     0.00006
 MMSS   7.24000000000000     
    13    13.07244    13.12345    13.12333     0.00012
 MMSS   30.0000000000000     
    14    21.30000    21.50000    21.50000     0.00000
 MMSS   47.2100000000000     
    15    59.47211    59.78920    59.78917     0.00003
 MMSS   11.0100000000000     
    16    65.11010    65.18600    65.18361     0.00239

Do you have any idea what's going on? Why should rounding to more places cause it to round down instead of up?

Here is the code at this point:

Code:


      DOUBLE PRECISION FUNCTION DMS2DD (DMS)
      REAL*8 DMS
C
      REAL*8 DROUND
      REAL*8 DD, MM, SS
      REAL*8 MMSS
C
      DD = IDINT(DMS)
C
      MMSS = 100.0D0 * DROUND(DMS - DD, 4)
      WRITE (*,*) 'MMSS', MMSS
      MM = IDINT(MMSS)
      SS = (MMSS - MM) * 100.0D0
      DMS2DD = DD + MM / 60.0D0 + SS / 3600.0D0
      RETURN
      END

Code:


      DOUBLE PRECISION FUNCTION DROUND (X, P)
      REAL*8 X
      INTEGER P
C
      REAL*8 R
C
      R = IDINT(10.0D0 ** P + 0.5D0)
      DROUND = IDINT(R * X + 0.5D0) / R
      RETURN
      END

[Dieter]

(03-30-2016 01:10 PM)HP67 Wrote:  Do you have any idea what's going on? Why should rounding to more places cause it to round down instead of up?

Yes, I think I know what's going on here: Your Fortran seems to support a precision of about ten decimal digits. That's why rounding to 13 digits does not help. I think your DROUND function is fine, just replace the 4 digit to 8 digit rounding. Or, if you are sure that nothing beyond 1/100 second occurs, simply set it to 6 digit rounding.

You may test this easily: simply write a short program that prints a periodic value like 1/3 or 1/9 with 15 digits. For instance, simply let the program print 1.0D0/3.0D0 with 15 decimal places and see what you get. If it's something like 0,33333 33334 8855... you know that 31 bits or 9-10 decimal places is all you got. ;-)

(03-30-2016 01:10 PM)HP67 Wrote:  If I specify rounding to 4 places the results improve a lot but are still worse than when unrounded:

Sure. Rounding to 4 places will work fine if 0.mmss does not have more than 4 digits, i.e. no fractions of a second. And indeed for these cases your results are fine, only the "EXPECT" values are wrong again: dms2dd(71,0030) = 71,0083333.... and dms2dd(0,4949) = 0,83027777...., or dms2dd(65,11010) = 65,18361111... which is exactly what your program returns. The "expected" values 71,0080, 0,83020 and 65.18600 are WRONG!

Actually all cases with not more than 4 decimals are perfectly correct. Even those with 5 places (i.e. fractions of a second) are exact if the argument is rounded to full seconds (4 digits). For instance, dms2dd(13,07244) actually evaluates dms2dd(13,0724), so mm.ss = 7,2400000 and the result for this is 13,12333. Adjust the roundoff function to 5 places (or 6, or 8) and these will return exact results as well.

Now please try rounding to 6 or 8 digits and report what you got.
And, more important, check your results with a reliable reference, e.g. an HP. ;-)

Dieter

[HP67]

(03-30-2016 03:41 PM)Dieter Wrote:  Yes, I think I know what's going on here: Your Fortran seems to support a precision of about ten decimal digits. That's why rounding to 13 digits does not help. I think your DROUND function is fine, just replace the 4 digit to 8 digit rounding. Or, if you are sure that nothing beyond 1/100 second occurs, simply set it to 6 digit rounding.

According to the doc it's 14 hexadecimal digits which is supposed to be approximately equivalent to 16.8 decimal digits. I don't understand how it's encoded yet so I can't verify this. But it's IBM and the doc is never wrong. It's just not very helpful in this manual. I have another couple places I want to check. So, the 13 digits you mentioned should have been perfect.

(03-30-2016 03:41 PM)Dieter Wrote:  You may test this easily: simply write a short program that prints a periodic value like 1/3 or 1/9 with 15 digits. For instance, simply let the program print 1.0D0/3.0D0 with 15 decimal places and see what you get. If it's something like 0,33333 33334 8855... you know that 31 bits or 9-10 decimal places is all you got. ;-)

Yes, and more to the point I'll run through a loop printing the output of the rounding function for 1/3 for 1..16 digits. [face palm] why didn't I think of this..

(03-30-2016 01:10 PM)HP67 Wrote:  If I specify rounding to 4 places the results improve a lot but are still worse than when unrounded:

(03-30-2016 03:41 PM)Dieter Wrote:  Sure. Rounding to 4 places will work fine if 0.mmss does not have more than 4 digits, i.e. no fractions of a second. And indeed for these cases your results are fine, only the "EXPECT" values are wrong again: dms2dd(71,0030) = 71,0083333.... and dms2dd(0,4949) = 0,83027777...., or dms2dd(65,11010) = 65,18361111... which is exactly what your program returns. The "expected" values 71,0080, 0,83020 and 65.18600 are WRONG!

They should be correct in the latest run. I checked all the calculations on an HP 48 and made a few changes before I started working on this today. I think the problem is lack of precision in the printout and the expected values are really correct but I'll read your post again and check it against my test. I'm only displaying 5 decimal places and the output routines round. I'll extend the size of the input values in the next printout. I was only using 5 decimal places because that will be enough for my end use of this routine. For debugging we should probably show more.

(03-30-2016 03:41 PM)Dieter Wrote:  Even those with 5 places (i.e. fractions of a second) are exact if the argument is rounded to full seconds (4 digits). For instance, dms2dd(13,07244) actually evaluates dms2dd(13,0724), so mm.ss = 7,2400000 and the result for this is 13,12333. Adjust the roundoff function to 5 places (or 6, or 8) and these will return exact results as well.

I don't think this is working. Although the input, expected, and actual results are shown they are rounded by the output routines. The difference value is valid at least to 5 decimal places rounded. In many cases the difference is non-zero within that precision so I think there is a still a problem.

(03-30-2016 03:41 PM)Dieter Wrote:  Now please try rounding to 6 or 8 digits and report what you got.

Yes, I'll run the test and post again tomorrow hopefully.

(03-30-2016 03:41 PM)Dieter Wrote:  And, more important, check your results with a reliable reference, e.g. an HP. ;-)

If I had a fullsize keyboard on my HP calcs I would probably throw my PC out the window!

Thank you very much for helping me with this. It's so frustrating not to get it when it comes to math.

[Dieter]

(03-30-2016 04:21 PM)HP67 Wrote:  According to the doc it's 14 hexadecimal digits which is supposed to be approximately equivalent to 16.8 decimal digits. I don't understand how it's encoded yet so I can't verify this. But it's IBM and the doc is never wrong.

I am quite confident that there's not more than 31 or 32 bit. Take a look at the slightly inaccurate MMSS values you listed. These can be accurately represented in binary notation with 31 bits or less:

Code:

30.299999999348074 = 65,068,754,533 / 2^31

89.111499999649823 = 95,682,744,549 / 2^30

12.149999999441206 = 13,045,963,161 / 2^30

38.422499999403954 =  2,578,490,327 / 2^26

29.302999999374151 =  7,865,964,167 / 2^28

32 bits mean 9,3 decimal digits, which matches the results. Since 2 digits are used for the integer portion, rounding the mmss part to 7 decimals shoud be OK.

(03-30-2016 04:21 PM)HP67 Wrote:  Yes, and more to the point I'll run through a loop printing the output of the rounding function for 1/3 for 1..16 digits. [face palm] why didn't I think of this..

Sounds like a good idea.

Ideally the number of rounded digits should depend on the integer part of the dms argument. 2 digits => round to 7 decimals. 3 digits => round to 6 decimals. 0 digits => round to 9 decimals. More generally: round to 8 – floor(log10(dms)) decimals. In other words: round to 9 significant (!) decimal digits. The dms argument obviously is not represented more accurately anyway.

(03-30-2016 04:21 PM)HP67 Wrote:  They should be correct in the latest run. I checked all the calculations on an HP 48 and made a few changes before I started working on this today. I think the problem is lack of precision in the printout and the expected values are really correct but I'll read your post again and check it against my test.

Well, 71°00'30" is 71 + 30/3600 = 71,008333333... degrees and not 71,0080°. There's not much room for different interpretations. ;-)

(03-30-2016 04:21 PM)HP67 Wrote:  

(03-30-2016 03:41 PM)Dieter Wrote:  Even those with 5 places (i.e. fractions of a second) are exact if the argument is rounded to full seconds (4 digits). For instance, dms2dd(13,07244) actually evaluates dms2dd(13,0724), so mm.ss = 7,2400000 and the result for this is 13,12333. Adjust the roundoff function to 5 places (or 6, or 8) and these will return exact results as well.

I don't think this is working.

I do. :-)

Dieter

[rprosperi]

(03-30-2016 06:38 PM)Dieter Wrote:  32 bits mean 9,3 decimal digits, which matches the results. Since 2 digits are used for the integer portion, rounding the mmss part to 7 decimals shoud be OK.

Dieter - what is the meaning of "9,3 decimal digits" ? I get that due to internal representation inaccuracies you can only get somewhere between 9 and 10 reliable decimal digits, but what is the meaning or interpretation of the ",3" ?

These threads are interesting to follow, I'm always learning something new, but I've not seen this usage before. Or maybe have simply forgotten it in the last 30+ years since doing this stuff on a daily basis.

[Dieter]

(03-30-2016 08:32 PM)rprosperi Wrote:  Dieter - what is the meaning of "9,3 decimal digits" ? I get that due to internal representation inaccuracies you can only get somewhere between 9 and 10 reliable decimal digits, but what is the meaning or interpretation of the ",3" ?

In mathematics, results here and there may not exactly match real life. ;-)

The meaning is exactly as you suggested: "somewhere between 9 and 10 reliable decimal digits". Say, 9 digits with an error less than 1 ULP. I think it should be about 0,5 ULP in this case, as 2^–31 is approx. 4,6 E–10. Or 2^–32 = 2,3 E–10 for 32 bits. Which BTW should translate to 9,6 digits (sorry).

(03-30-2016 08:32 PM)rprosperi Wrote:  These threads are interesting to follow, I'm always learning something new, but I've not seen this usage before. Or maybe have simply forgotten it in the last 30+ years since doing this stuff on a daily basis.

Specifying a number of valid digits with fractional values is quite common. For instance in discussions on the accuracy of approximation methods. I remember such notations ("10,7 digits") in a thread on various methods for the Gamma function. In that case it was simply the base-10 log of the relative error. ;-)

Dieter

[HP67]

(03-30-2016 06:38 PM)Dieter Wrote:  I am quite confident that there's not more than 31 or 32 bit. Take a look at the slightly inaccurate MMSS values you listed. These can be accurately represented in binary notation with 31 bits or less:

I certainly cannot and am not arguing that your observation is not correct. Still, there is 56 bits of precision here for an 8 byte real. 56 bits for the fraction, 1 for the sign, and 7 for the exponent. I'm testing on two different compilers on two different hardware architectures (although the former being emulated on the latter) and both get similar results with what is supposed to be 8 byte reals but totally different data formats. I suspect an error in my code like an incorrect default declaration being used (resulting in a 4 byte real) but I can't find it. This should not be happening. It has been bothering me all along that the results should be this far off. Something is definitely wrong here.

(03-30-2016 04:21 PM)HP67 Wrote:  Yes, and more to the point I'll run through a loop printing the output of the rounding function for 1/3 for 1..16 digits. [face palm] why didn't I think of this..

Results shown below.

(03-30-2016 04:21 PM)HP67 Wrote:  They should be correct in the latest run. I checked all the calculations on an HP 48 and made a few changes before I started working on this today. I think the problem is lack of precision in the printout and the expected values are really correct but I'll read your post again and check it against my test.

(03-30-2016 06:38 PM)Dieter Wrote:  Well, 71°00'30" is 71 + 30/3600 = 71,008333333... degrees and not 71,0080°. There's not much room for different interpretations. ;-)

What I was saying is that the output shown was rounded by the FORTRAN output routines. I increased the length of the output format and you can see this is not a real problem. There was one error remaining in the last test (dyslexia on my part) but that is now fixed.

(03-30-2016 04:21 PM)HP67 Wrote:  I don't think this is working.

(03-30-2016 06:38 PM)Dieter Wrote:  I do. :-)

Please see below. With rounding to 5 places some of the values are significantly off. I also tried using 9 places and 10 and it didn't help.

Code:


                C
                C TEST DMSDD
                C
   ISN 0002           REAL*8 DMS2DD                                                     00080000
                C                                                                       00090000
   ISN 0003           REAL*8 DMS(16) /89.1115D0, 12.1500D0, 33.3000D0, 71.0030D0,       00100000
                     X                42.2453D0, 38.4225D0, 29.3030D0, 00.4949D0,       00110000
                     X                75.1500D0, 43.2230D0, 09.3345D0, 33.57522D0,      00120000
                     X                13.072442D0, 21.3000D0, 59.472112D0, 65.110096D0/ 00130000
                C                                                                       00140005
   ISN 0004           REAL*8 EXPECT(16) /89.1875D0, 12.2500D0, 33.5000D0, 71.0083D0,    00150000
                     X                   42.4147D0, 38.7069D0, 29.5083D0, 00.8302D0,    00160000
                     X                   75.25D0, 43.375D0, 09.5625D0, 33.9645D0,       00170000
                     X                   13.12345D0, 21.5D0, 59.7892D0, 65.1836D0/      00180005
                C                                                                       00190000
   ISN 0005           REAL*8  DIFF, RESULT                                              00200000
   ISN 0006           INTEGER I                                                         00210000
                C                                                                       00220000
   ISN 0007           WRITE (6,5)                                                       00230000
   ISN 0008           WRITE (6,10)                                                      00240000
   ISN 0009         5 FORMAT (/' ', 'T1DMS2DD'/)                                        00250000
   ISN 0010        10 FORMAT(' ',1X,'TEST',10X,'DMS',7X,'EXPECT',7X,'RESULT',9X,'DIFF'/)00260004
                C                                                                       00270000
   ISN 0011           DO 100 I = 1, 16                                                  00280000
   ISN 0012           RESULT =  DMS2DD(DMS(I))                                          00290000
   ISN 0013           DIFF = DABS(RESULT - EXPECT(I))                                   00300000
   ISN 0014           WRITE (6,50) I, DMS(I), EXPECT(I), RESULT, DIFF                   00310000
   ISN 0015        50 FORMAT (' ',1X,I4,4(3X,F10.7))                                    00320000
   ISN 0016       100 CONTINUE                                                          00330000
   ISN 0017           STOP                                                              00340000
   ISN 0018           END                                                               00350000

   ISN 0002           DOUBLE PRECISION FUNCTION DMS2DD (DMS)                            00060000
   ISN 0003           REAL*8 DMS                                                        00070000
                C                                                                       00080000
   ISN 0004           REAL*8 DROUND                                                     00090000
   ISN 0005           REAL*8 DD, MM, SS                                                 00100000
   ISN 0006           REAL*8 MMSS                                                       00110000
                C                                                                       00120000
   ISN 0007           DD = IDINT(DMS)                                                   00130000
                C                                                                       00140000
   ISN 0008           MMSS = 100.0D0 * DROUND(DMS - DD, 5)                              00150005
   ISN 0009           MM = IDINT(MMSS)                                                  00160000
   ISN 0010           SS = (MMSS - MM) * 100.0D0                                        00170000
   ISN 0011           DMS2DD = DD + MM / 60.0D0 + SS / 3600.0D0                         00180000
   ISN 0012           RETURN                                                            00190000
   ISN 0013           END                                                               00200000

   ISN 0002           DOUBLE PRECISION FUNCTION DROUND (X, P)                           00040000
   ISN 0003           REAL*8 X                                                          00050000
   ISN 0004           INTEGER P                                                         00060000
                C                                                                       00070000
   ISN 0005           REAL*8 R                                                          00080000
                C                                                                       00090000
   ISN 0006           R = IDINT(10.0D0 ** P + 0.5D0)                                    00100000
   ISN 0007           DROUND = IDINT(R * X + 0.5D0) / R                                 00110000
   ISN 0008           RETURN                                                            00120000
   ISN 0009           END                                                               00130000

T1DMS2DD

 TEST          DMS       EXPECT       RESULT         DIFF

    1   89.1115000   89.1875000   89.1875000    0.0000000
    2   12.1500000   12.2500000   12.2611111    0.0111111
    3   33.3000000   33.5000000   33.5111111    0.0111111
    4   71.0030000   71.0083000   71.0083333    0.0000333
    5   42.2453000   42.4147000   42.4147222    0.0000222
    6   38.4225000   38.7069000   38.7069444    0.0000444
    7   29.3030000   29.5083000   29.5083333    0.0000333
    8    0.4949000    0.8302000    0.8302778    0.0000778
    9   75.1500000   75.2500000   75.2611111    0.0111111
   10   43.2230000   43.3750000   43.3750000    0.0000000
   11    9.3345000    9.5625000    9.5625000    0.0000000
   12   33.5752200   33.9645000   33.9645000    0.0
   13   13.0724420   13.1234500   13.1234444    0.0000056
   14   21.3000000   21.5000000   21.5111111    0.0111111
   15   59.4721120   59.7892000   59.7891944    0.0000056
   16   65.1100960   65.1836000   65.1836111    0.0000111


------------------------------------------------------------------------------------------------
   ISN 0002           REAL*8 DROUND                                                     00040000
                C                                                                       00050000
   ISN 0003           REAL*8 V                                                          00060000
   ISN 0004           REAL*8 A                                                          00070000
   ISN 0005           INTEGER P                                                         00080000
                C                                                                       00090000
   ISN 0006           V = 1.0D0 / 3.0D0                                                 00100000
                C                                                                       00110000
   ISN 0007           WRITE (6,5)                                                       00120000
   ISN 0008         5 FORMAT (' ', 'T1ROUND'/)                                          00130000
   ISN 0009           WRITE (6,10)                                                      00140000
   ISN 0010        10 FORMAT (' ','PLACES                INPUT                ROUNDED') 00150001
   ISN 0011           DO 100 P = 1, 16                                                  00160000
   ISN 0012           A = DROUND (V, P)                                                 00170000
   ISN 0013           WRITE (6,20) P, V, A                                              00180000
   ISN 0014        20 FORMAT (' ',I4,2(5X,F18.15))                                      00190000
   ISN 0015       100 CONTINUE                                                          00200000
   ISN 0016           STOP                                                              00210000
   ISN 0017           END                                                               00220000


T1ROUND

PLACES                INPUT                ROUNDED
   1      0.333333333333333      0.300000000000000
   2      0.333333333333333      0.330000000000000
   3      0.333333333333333      0.333000000000000
   4      0.333333333333333      0.333300000000000
   5      0.333333333333333      0.333330000000000
   6      0.333333333333333      0.333333000000000
   7      0.333333333333333      0.333333300000000
   8      0.333333333333333      0.333333330000000
   9      0.333333333333333      0.333333333000000
  10      0.333333333333333      0.333333333569729
  11      0.333333333333333      0.333333333607512
  12      0.333333333333333      0.333333331500270
  13      0.333333333333333      0.333333333080067
  14      0.333333333333333      0.333333332127558
  15      0.333333333333333      0.333333332679950
  16      0.333333333333333      0.333333333155548

[Dieter]

(03-31-2016 09:00 AM)HP67 Wrote:  [... a lot ]

Thanks for all the input. Here and there I see clearer now.

• Your Fortran compiler supports extended precision, at least 15 digits. This can be seen in the input column of the 1/3 test. The variable V that holds 1,0D0/3,0D0 has the correct value 0,33333 33333 33333.
• The rounded output however does not hold more that 9 valid digits. So some precision was lost inside the DROUND function. So the limited accuracy is lost by something Fortran-specific.
  

As I said before, I do not have any experience with Fortran. But from what I read, mixing several data types in a calculation can cause weird results. In this case the DROUND function has two major critical points:

• R = IDINT(10.0D0 ** P + 0.5D0)
  Here R is declared as a double precision real. P is an (two-byte?) integer and IDINT is said to return a 4-byte integer (i.e. something with at most 10 digits). So what is the resulting precision for R? Maybe you can do the 1/3 test again and have DROUND return R to that it is listed in the table where now the rounded 1/3 values appear.
• DROUND = IDINT(R * X + 0.5D0) / R
  Again, at this moment we do not know the precision of R. And again IDINT should return a 4-byte integer with at most 10 digits. This could be the reason for the reduced accuracy.
  

This is the point where I fear I cannot be of much help any longer, as further debugging requires knowledge in Fortran. But maybe there are some experts here.

Dieter

[HP67]

Hi and thank you very much. I was out for a walk and I realized the problem was staring me in the face. I was hoping I could update the thread before you saw it. Yes, it's the IDINT function and it happens in DROUND as you suspected.

FORTRAN IV's IDINT returns a 4 byte integer. There is no 8 byte INT() function in FORTRAN IV.

I see three possible approaches:

1) If you can define ROUND in terms of MOD (there is an 8 byte mod function available) or some other way that doesn't require taking the integer portion with the IDINT function then we can fix the rounding routine. The magnitude of the dd.mmss or hh.mmss data itself never gets close to the 2^31 signed integer limit.

2) If you know an algorithm we could substitute for the existing IDINT function. If we can go through some process to find the integer portion we can do it in 8 byte arithmetic.

3) If I can figure out the instruction and data formats for this I should be able to write an 8 byte truncation function and I wouldn't need IDINT.

Some further updates on specific points:

(03-31-2016 12:39 PM)Dieter Wrote:  But from what I read, mixing several data types in a calculation can cause weird results.

This is generally safe if you pay attention and specify the correct constants and don't use probably the one builtin that downgrades precision, IDINT! Certainly the value of the RHS is promoted to the type of the LHS. That doesn't always help, obviously.

(03-31-2016 12:39 PM)Dieter Wrote:  R = IDINT(10.0D0 ** P + 0.5D0)
Here R is declared as a double precision real. P is an (two-byte?) integer and IDINT is said to return a 4-byte integer (i.e. something with at most 10 digits). So what is the resulting precision for R?

The exponent is a 4 byte integer- I checked earlier that the result of an 8 byte real raised to any legal variable is still an 8 byte real. That part is fine since we're using integral powers less than 2^31 here.

What is the precision of R? R is a double precision (8 byte) floating point value that was promoted from a single precision 4 byte integer :-(

[Dieter]

(03-31-2016 12:51 PM)HP67 Wrote:  FORTRAN IV's IDINT returns a 4 byte integer. There is no 8 byte INT() function in FORTRAN IV.

That's the problem.

(03-31-2016 12:51 PM)HP67 Wrote:  I see three possible approaches:

I see a fourth one: simply do it all in integer arithmetics.

Right at the beginning the dms input is converted to an integer. If the largest 4-bit integer is 2147483647, the largest dms input is 214,7483647. So we just multiply the input by 10.000.000 and use that as an integer.

I tried this in Excel VBA and it works well:

Code:

Function dms2dd(dms as double)
Dim x, ddd, mmss, mm, ss As Long

x = dms * 10000000
ddd = x \ 10000000
mmss = x Mod 10000000
mm = mmss \ 100000
ss = mmss Mod 100000

dms2dd = ddd + mm / 60# + ss / 3600000#

End Function

Note: the "\" operator means integer division. The "#" suffix marks a variable or value as double precision. I think you would use 60.0D0 and 3.6D6.

The major limitation of this method is the domain of the dms argument. The largest angle is 214 degrees and the precision of dms is 0.001 seconds. In other words, everything beyond the 7th decimal is neglected. This way the result has 5-6 valid decimals (0,001/3600 =1,7E–6).

You may check if Fortran truncates or rounds x when it is converted in the first line (x=dms*10000000). Maybe you will have to code this as X = IDINT(DMS*1.0D7+0.5D0).

What do you think?

Dieter

[HP67]

Thanks very much.

Mostly I'm using this to convert decimal time back and forth. .001 seconds should be plenty of resolution. I'll be happy with a reliable hh.mmsst for now.

I think a rounding function and a double precision truncation function would be more generally useful but I'll study what you came up with. I can anticipate a few problems but fortunately I am almost always wrong about math so maybe this time it will work in my favor.

FORTRAN IV doesn't have any integer division except for integer types and doesn't have 8 byte integers. When we get back to dividing again won't we go through the same or similar conversions/promotions as before? So I am expecting problems. Again, I'm probably wrong.

I'll play around with your suggestion. Eventually I would like to come up with an 8 byte truncation function for reals. The rounding function is very handy also.

In understanding your declarations, Excel VBA long means 8 byte int? No comparable type in FORTRAN IV. Will a 4 byte int work? I think you explained that it will, so I don't understand why you chose long in your code.