FORTRAN floating point accuracy problems

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