FORTRAN floating point accuracy problems

[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