(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
Search
Member List
Calendar
Help
/
FORTRAN floating point accuracy problems
