FORTRAN floating point accuracy problems

[Dieter]

(03-31-2016 02:45 PM)HP67 Wrote:  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.

OK, so ddd actually is hhh. ;-) What is the largest value you have to handle? Also, possibly two more digits could be obtained if the integer and fractional part of dms would be handled separately. But if 0.001 s is fine, the posted code is an easy solution.

(03-31-2016 02:45 PM)HP67 Wrote:  FORTRAN IV doesn't have any integer division except for integer types and doesn't have 8 byte integers.

That's exactly how the VBA code works. "Long" is a signed four-byte integer, and "\" is what Fortran does with a simple "/" division of two integers (12/5 = 2).

(03-31-2016 02:45 PM)HP67 Wrote:  In understanding your declarations, Excel VBA long means 8 byte int?

No, that's the same four bytes, i.e. integers within ±2^31. You should be able to translate this more or less 1:1 into Fortran. "Long" just means "Long integer", compared to a standard integer within ±32768.

Dieter

[Dieter]

(03-31-2016 07:30 PM)Dieter Wrote:  Also, possibly two more digits could be obtained if the integer and fractional part of dms would be handled separately.

This could be a possible solution then:

Code:

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

ddd = Int(dms)
mmss = (dms - ddd) * 1000000000
mm = mmss \ 10000000
ss = mmss Mod 10000000

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

End Function

Now nine decimals of dms are used for the conversion, i.e. the resolution is 0.00001 s. There is (almost) no limit for dms (max. ~2 billion hours).
The result should have eight valid digits after the decimal point.

Please note that VBA automatically rounds the result of "mmss = (dms - ddd) * 1000000000" to the nearest integer.
You may have to code this as MMSS = IDINT((DMS–DDD)*1.0D9 + 0.5D0).

Dieter

[HP67]

Thanks, can't respond now, will get back to you next week. Have a great weekend!

[HP67]

Still having problems with specific input values. I tried my own version before I looked at yours and I tried yours also. Unsurprisingly yours worked a little better. But we are still having problems with specific values. BTW the output is shown to 7 decimal places but the input and expected values are only to 6 values. I showed extra digits since last time it seemed important.

Aside from specific problematic values (.15, .3) this works better than required. So, in that sense much progress. Thanks.

Code:


 LEVEL 21.8 ( JUN 74 )                          OS/360  FORTRAN H                                 DATE  16.094/15.37.38
           COMPILER OPTIONS - NAME=  MAIN,OPT=02,LINECNT=50,SIZE=0000K,
                              SOURCE,EBCDIC,LIST,DECK,NOLOAD,MAP,EDIT,NOID,XREF
   ISN 0002           DOUBLE PRECISION FUNCTION DMS2DD (DMS)                            00010013
   ISN 0003           REAL*8 DMS, DDD                                                   00020013
   ISN 0004           INTEGER MMSS, MM, SS                                              00030013
                C                                                                       00040013
   ISN 0005           DDD = IDINT(DMS)                                                  00050013
   ISN 0006           MMSS = (DMS - DDD) * 1000000000                                   00060013
   ISN 0007           MM = MMSS / 10000000                                              00070013
   ISN 0008           SS = MOD(MMSS, 10000000)                                          00080013
                C                                                                       00090013
   ISN 0009           DMS2DD = DDD + MM / 60.0D0 + SS / 360000000.0D0                   00100013
   ISN 0010           RETURN                                                            00110013
   ISN 0011           END                                                               00120013

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.0083330   71.0083333    0.0000003
    5   42.2453000   42.4147220   42.4147222    0.0000002
    6   38.4225000   38.7069440   38.7069444    0.0000004
    7   29.3030000   29.5083330   29.5083333    0.0000003
    8    0.4949000    0.8302770    0.8302778    0.0000008
    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.0000000
   13   13.0724420   13.1234500   13.1234500    0.0000000
   14   21.3000000   21.5000000   21.5111111    0.0111111
   15   59.4721120   59.7892000   59.7892000    0.0000000
   16   65.1100960   65.1836000   65.1836000    0.0000000
******************************** BOTTOM OF DATA ***************************************************************************​*********
 LEVEL 21.8 ( JUN 74 )                          OS/360  FORTRAN H                                 DATE  16.094/15.38.43
           COMPILER OPTIONS - NAME=  MAIN,OPT=02,LINECNT=50,SIZE=0000K,
                              SOURCE,EBCDIC,LIST,DECK,NOLOAD,MAP,EDIT,NOID,XREF
   ISN 0002           DOUBLE PRECISION FUNCTION DMS2DD (DMS)                            00060013
   ISN 0003           REAL*8 DMS                                                        00070013
                C                                                                       00080013
   ISN 0004           INTEGER IDMS, DD, MM, SS                                          00090013
                C                                                                       00100013
   ISN 0005           IDMS = DMS * 1.0D7                                                00110013
                C                                                                       00120013
   ISN 0006           DD = IDMS / 1.0D7                                                 00130013
   ISN 0007           MM = (IDMS / 1.0D5) - (DD * 100)                                  00140013
   ISN 0008           SS = IDMS - (DD * 1.0D7) - (MM * 1.0D5)                           00150013
                C                                                                       00160013
   ISN 0009           DMS2DD = DD + MM / 60.0D0 + SS / 36.0D5                           00170013
   ISN 0010           RETURN                                                            00190013
   ISN 0011           END                                                               00200013

T1DMS2DD

 TEST          DMS       EXPECT       RESULT         DIFF

    1   89.1115000   89.1875000   89.1874997    0.0000003
    2   12.1500000   12.2500000   12.2611108    0.0111108
    3   33.3000000   33.5000000   33.5111108    0.0111108
    4   71.0030000   71.0083330   71.0083331    0.0000001
    5   42.2453000   42.4147220   42.4147219    0.0000001
    6   38.4225000   38.7069440   38.7069442    0.0000002
    7   29.3030000   29.5083330   29.5083331    0.0000001
    8    0.4949000    0.8302770    0.8302775    0.0000005
    9   75.1500000   75.2500000   75.2611108    0.0111108
   10   43.2230000   43.3750000   43.3749997    0.0000003
   11    9.3345000    9.5625000    9.5624997    0.0000003
   12   33.5752200   33.9645000   33.9644997    0.0000003
   13   13.0724420   13.1234500   13.1234497    0.0000003
   14   21.3000000   21.5000000   21.5111108    0.0111108
   15   59.4721120   59.7892000   59.7891997    0.0000003
   16   65.1100960   65.1836000   65.1835997    0.0000003
******************************** BOTTOM OF DATA ***************************************************************************​*********

This floating point hardware implementation uses a base 16 internal representation. Do you have any idea how we might change the algorithm to use powers of 16 instead of powers of 10 or if that has any chance of working better?

I want to figure out how to code a double precision (8 byte) truncation function as I mentioned earlier. Your rounding routine is really great. If my implementation wasn't hampered by the existing 4 byte truncation function maybe all this would work a lot better?

I have figured out how to convert the internal representation to decimal but there is some loss of precision on my HP calcs since they only carry 12 significant digits. I have not figured out enough about the math behind this to understand what I need to do to truncate the decimal portion while it's still in the internal base 16 format. I'm sure anybody mathematically inclined could see it right away but I'm still banging my head.

Update: I see an instruction that could help with the the 8 byte truncation function. Will update on that in the next few days hopefully.

[Dieter]

(04-03-2016 03:49 PM)HP67 Wrote:  Still having problems with specific input values. I tried my own version before I looked at yours and I tried yours also. Unsurprisingly yours worked a little better. But we are still having problems with specific values.

The problem seems to be the integer values mmss, mm and ss which are truncated so that a true result of mmss in case 3 may actually be 29.999999... instead of 30, which is then truncated to 29 minutes and 99.999... seconds. In VBA this problem does not occur as an assignment to an integer variable automatically rounds a real up or down to the next integer. Adding some explicit rounding should fix that:

You could replace line 0006
MMSS = (DMS - DDD) * 1000000000
with
MMSS = (DMS - DDD) * 1000000000 + 0.5D0
or even better
MMSS = IDINT((DMS - DDD) * 1000000000 + 0.5D0)

And once again the "expected" values are wrong. For instance

Code:

 TEST          DMS       EXPECT       RESULT         DIFF
   ...
    4   71.0030000   71.0083330   71.0083333    0.0000003
    5   42.2453000   42.4147220   42.4147222    0.0000002
    6   38.4225000   38.7069440   38.7069444    0.0000004
    7   29.3030000   29.5083330   29.5083333    0.0000003
    8    0.4949000    0.8302770    0.8302778    0.0000008
   ...

...should actually read...

Code:

 TEST          DMS       EXPECT       RESULT         DIFF
   ...
    4   71.0030000   71.0083333   71.0083333    0.0000000
    5   42.2453000   42.4147222   42.4147222    0.0000000
    6   38.4225000   38.7069444   38.7069444    0.0000000
    7   29.3030000   29.5083333   29.5083333    0.0000000
    8    0.4949000    0.8302778    0.8302778    0.0000000
   ...

All these results are dead on. Actually only case 2, 3, 9 and 14 are off, and I think this can be fixed by adding the suggested integer rounding.

The second program version has essentially the same issues:
Line 0005
IDMS = DMS * 1.0D7
should better be rounded to
IDMS = DMS * 1.0D7 + 0.5D0
or even better
IDMS = IDINT(DMS * 1.0D7 + 0.5D0)

(04-03-2016 03:49 PM)HP67 Wrote:  I want to figure out how to code a double precision (8 byte) truncation function as I mentioned earlier. Your rounding routine is really great. If my implementation wasn't hampered by the existing 4 byte truncation function maybe all this would work a lot better?

Most probably, yes.

Dieter

[HP67]

Again, I am only carrying 6 decimal places of expected results although I am displaying 7.

Thanks, I'll try your suggestions. But first I'm going to try to write an assembler subroutine for DINT. Then I can use your rounding routine which would really be great.

[Dieter]

(04-03-2016 05:17 PM)HP67 Wrote:  Again, I am only carrying 6 decimal places of expected results although I am displaying 7.

This means that "differences" in the 7th place actually are not present. ;-)
Maybe you can update the "expect" values to seven decimals.

(04-03-2016 05:17 PM)HP67 Wrote:  Thanks, I'll try your suggestions. But first I'm going to try to write an assembler subroutine for DINT. Then I can use your rounding routine which would really be great.

Hmmm... I would love to know whether the suggested rounding fixes the issue. Maybe you can add these little changes and report the results... ?-)

Dieter

[Werner]

Can I chime in?

The following formulas, courtesy of John H. Meyers on comp.sys.hp48 do not require any fiddling at all:

HR2HMS

T := HR - FP(HR)*0.4;
HMS := T - FP(T*100)*0.004;

HMS2HR

T := HMS + FP(HMS*100)/150;
HR := T + FP(T)/1.5;

Tried them on Free42S Binary. Both work like a charm.

Code:

00 { 55-Byte Prgm }
01>LBL "HMS"
02 ENTER
03 FP
04 0.4
05 ×
06 -
07 1E4
08 %
09 FP
10 4E-3
11 ×
12 -
13 RTN
14>LBL "HR"
15 1E4
16 %
17 FP
18 150
19 ÷
20 +
21 ENTER
22 FP
23 1.5
24 ÷
25 +
26 .END.

Cheers, Werner

[HP67]

(04-03-2016 05:21 PM)Dieter Wrote:  

(04-03-2016 05:17 PM)HP67 Wrote:  Again, I am only carrying 6 decimal places of expected results although I am displaying 7.

This means that "differences" in the 7th place actually are not present. ;-)

Yes I know. Therefore I was not troubled by the diff values ;-)

(04-03-2016 05:21 PM)Dieter Wrote:  Maybe you can update the "expect" values to seven decimals.

I'll have to look at it. It could be with 7 decimal places I can't fit 4 values on one line of FORTRAN code...

(04-03-2016 05:17 PM)HP67 Wrote:  Thanks, I'll try your suggestions. But first I'm going to try to write an assembler subroutine for DINT. Then I can use your rounding routine which would really be great.

(04-03-2016 05:21 PM)Dieter Wrote:  Hmmm... I would love to know whether the suggested rounding fixes the issue. Maybe you can add these little changes and report the results... ?-)

The good news is your rounding routine works perfectly up to 15 decimal places (possibly more, I will check) with my brand new shiny double precision (8 byte) DINT routine in assembler. I didn't have time to go back to our old versions of the driver code to see if it works better. I'll get back to that as soon as I can.

Thanks for the help! I am really excited about having the rounding routine in my FORTRAN IV toolbag. I have wanted one for a long time and just couldn't figure out how to do it. I do have some rounding code I use for assembler control blocks but I just couldn't extrapolate it to real numbers.

[HP67]

(04-04-2016 09:02 AM)Werner Wrote:  Can I chime in?

Hi. I know how to do the conversions by hand and with a calculator, etc. I need this for some code I am writing in circa 1974 FORTRAN IV. The math isn't the issue. The problem is dealing with loss of precision caused by the floating point hardware implementation- which as far as I know isn't a problem on HP calculators or calculators generally since they use decimal arithmetic.

[Werner]

That's why I specifically tested them on Free42 Binary, which uses the binary floating point of the computer it's running on... all your examples ran perfectly.

Another issue is that since HH.MMSSssss cannot be correctly represented in binary floating-point, you shouldn't be using it? Use HHMMSS.ssss instead..

Cheers, Werner

[HP67]

The hardware implementation I'm using is not binary but I see where you are coming from and thanks for the effort. We have tried a few things and so far have been hampered by a few missing functions. We have got that mostly resolved. I just haven't had time to test again with the old versions of code that should work now.

[Dieter]

(04-04-2016 09:02 AM)Werner Wrote:  Can I chime in?

Sure. 8-)

(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;

Ah, great. Yes, now I remember these formulas, too. Although I cannot definitely say (yet) if this solves all problems due to roundoff by non-decimal number representations, I think this should work very well.

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.

Dieter

[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

[HP67]

I have not had time to get to try the piece of code Werner posted from John Meyers and will probably not have time to get back to it for several weeks now.

I did want to post an update for Dieter on his rounding routine suggestion. It works beautifully with my new double precision truncation function.

Code:


 LEVEL 21.8 ( JUN 74 )                          OS/360  FORTRAN H                                 DATE  16.096/12.47.47

           COMPILER OPTIONS - NAME=  MAIN,OPT=02,LINECNT=50,SIZE=0000K,
                              SOURCE,EBCDIC,LIST,DECK,NOLOAD,MAP,EDIT,NOID,XREF
   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 DINT, R                                                    00080003
                C                                                                       00090003
   ISN 0006           R = DINT(10.0D0 ** P + 0.5D0)                                     00100001
   ISN 0007           DROUND = DINT(R * X + 0.5D0) / R                                  00110001
   ISN 0008           RETURN                                                            00120000
   ISN 0009           END                                                               00130000


   0005              V = 1.0D0 / 3.0D0                                                 00100000

T1ROUND

PLACES                INPUT                ROUNDED
   1      0.3333333333333333      0.3000000000000000
   2      0.3333333333333333      0.3300000000000000
   3      0.3333333333333333      0.3330000000000000
   4      0.3333333333333333      0.3333000000000000
   5      0.3333333333333333      0.3333300000000000
   6      0.3333333333333333      0.3333330000000000
   7      0.3333333333333333      0.3333333000000000
   8      0.3333333333333333      0.3333333300000000
   9      0.3333333333333333      0.3333333330000000
  10      0.3333333333333333      0.3333333333000000
  11      0.3333333333333333      0.3333333333300000
  12      0.3333333333333333      0.3333333333330000
  13      0.3333333333333333      0.3333333333333000
  14      0.3333333333333333      0.3333333333333300
  15      0.3333333333333333      0.3333333333333330
  16      0.3333333333333333      0.3333333333333333


   0005              V = 1.25D0                                                        00100000

T2ROUND

PLACES                INPUT                ROUNDED
   1      1.2500000000000000      1.3000000000000000
   2      1.2500000000000000      1.2500000000000000
   3      1.2500000000000000      1.2500000000000000
   4      1.2500000000000000      1.2500000000000000
   5      1.2500000000000000      1.2500000000000000
   6      1.2500000000000000      1.2500000000000000
   7      1.2500000000000000      1.2500000000000000
   8      1.2500000000000000      1.2500000000000000
   9      1.2500000000000000      1.2500000000000000
  10      1.2500000000000000      1.2500000000000000
  11      1.2500000000000000      1.2500000000000000
  12      1.2500000000000000      1.2500000000000000
  13      1.2500000000000000      1.2500000000000000
  14      1.2500000000000000      1.2500000000000000
  15      1.2500000000000000      1.2500000000000000
  16      1.2500000000000000      1.2500000000000000

Thanks so much for the help! This is really a great tool to have, especially for FORTRAN IV where there is a very limited selection of builtins and no obvious way to do this for double precision values.

I tried using this routine with our old code from earlier in the post but it did not help. We were still off badly on the same values we missed before. For some reason several of the results seem to have been rounded in the wrong direction. I tested the 1.25D0 value specifically because this was one of the values that seems to have gotten misadjusted upwards instead of rounded. It clearly works fine here. This doesn't make sense right now so I have to track it down later when I have time. I'm sorry I haven't been able to get back to the Meyer's code yet. I'll be very busy for another several weeks before I can get to it.

Thanks guys.

[Werner]

I was waiting till you came up with a problematic case ;-)
Actually, however, depending on how you look at it:

- this is a case of wrong input (99 for seconds)
- the result is correct (try it on Free42 Decimal, for instance)
But I agree that you can't see it as it is rounded in the display..

BTW with x = 1e13/(5e13 + 1), on Free42 Binary, INT(100*x) delivers 19 and not 20 as you stated, and there are no negatives during the computation?

Cheers, Werner

[Dieter]

(04-07-2016 01:36 PM)Werner Wrote:  - this is a case of wrong input (99 for seconds)
- the result is correct (try it on Free42 Decimal, for instance)

In decimal arithmetics everything is easy.

The point here is that already the binary representation of the argument may be off (or "wrong", as you put it). Consider the case dms=30,2 or 151/5. The closest possible binary double-precision representation of this translates to 30,199 99999 99999 993 in decimal notation. Even if you manually enter 30,2 the program will work with 30,1999... so that all the mentioned problems occur. Rounding (up) the input to 15 significant digits solves this.

(04-07-2016 01:36 PM)Werner Wrote:  BTW with x = 1e13/(5e13 + 1), on Free42 Binary, INT(100*x) delivers 19 and not 20 as you stated,

That's an easy case that is handled correctly even by Excel. Replace the "E13" part by E14 or E15, depending on Free42's working precision, so that the "+1" occurs in the least significant digit.

Dieter

[Pekis]

Hello,

Perhaps you should try these formulas:

For x>0:
DEG->DMS(x)=(90*x+INT(60*x)+100*INT(x))/250

It's mathematically exact, but special care is needed in evaluating INT(60*x) on a calculator.
=> It seems better to round 60*x to last internal digit before applying INT.

DMS->DEG(x)=(250*x-INT(100*x)-60*INT(x))/90

See you.

[Werner]

(04-07-2016 02:28 PM)Dieter Wrote:  

(04-07-2016 01:36 PM)Werner Wrote:  BTW with x = 1e13/(5e13 + 1), on Free42 Binary, INT(100*x) delivers 19 and not 20 as you stated,

That's an easy case that is handled correctly even by Excel. Replace the "E13" part by E14 or E15, depending on Free42's working precision, so that the "+1" occurs in the least significant digit.

Dieter

But it's your example a few posts up :

Quote:Consider the case dms = 1/5,000 00000 00001

of which you claimed Excel's INT function of 100*dms yielded 20:

Quote:Now INT of this yields 20 exactly (not 19!)

Well, never mind. If you do want to use d.ms format, which can't be represented exactly in binary floating point, you'll need rounding.

Cheers, Werner

[Dieter]

(04-11-2016 11:41 AM)Werner Wrote:  But it's your example a few posts up :

There was one zero missing. #-)

The example is supposed to have the 1 in the last significant (i.e. 15th) digit. So that should be 1/5,0000 00000 00001. Or 1E14/(5E14+1). That's why I said you should replace the "E13" parts by "E14". OTOH the decimal value 0,19999 99999 99999 6 was correct.

In this case Excel indeed returns 20 for INT(100*dms).

(04-11-2016 11:41 AM)Werner Wrote:  Well, never mind. If you do want to use d.ms format, which can't be represented exactly in binary floating point, you'll need rounding.

That's what I said. And that's why I definitely prefer BCD arithmetics.

BCD is the standard on pocket calculators where this d.ms notation is common, and here it works well. On the other hand most programming languages use binary arithmetics, so here a different approach should be preferred. For instance with a record, a struct or whatever it's called in your preferred language.

Dieter