Print correct value of Real number with gfortran?

[Atr cheema]

The following program is printing wrong value of RWTSED. How can I print correct value??

      program inpdat
c

      IMPLICIT NONE
      REAL       RHOMN,RWTSED,VOLSED

     
      VOLSED = 17424.0
      RHOMN = 2.42   !0000076293945
      RWTSED= VOLSED*RHOMN*1.0E6

2000 FORMAT(/,3F40.5)
      WRITE(*,2000) RWTSED,VOLSED,RHOMN

      end

The answer printed on screen is
"42166083584.00000 17424.00000 2.42000"
The left most value is is definitely wrong.
How can I print the correct value of RWTSED?

 

[Vanadium 50]

The left most value is is definitely wrong.

Why do you say that? That's what my calculator gets.

 

[Atr cheema]

Why do you say that? That's what my calculator gets.

How can you get 17424*2.42*1e6 = 42166083584.00000?? It is 42166080000 indeed!

 

[Mark44]

The following program is printing wrong value of RWTSED. How can I print correct value??

      program inpdat
c

      IMPLICIT NONE
      REAL       RHOMN,RWTSED,VOLSED

   
      VOLSED = 17424.0
      RHOMN = 2.42   !0000076293945
      RWTSED= VOLSED*RHOMN*1.0E6

2000 FORMAT(/,3F40.5)
      WRITE(*,2000) RWTSED,VOLSED,RHOMN

      end

The answer printed on screen is
"42166083584.00000 17424.00000 2.42000"
The left most value is is definitely wrong.
How can I print the correct value of RWTSED?

Your program uses the REAL data type, which is four bytes in size. Precision is limited to about 7 decimal digits. For higher precision, use DOUBLE PRECISION (or REAL*8), which is good for about 17 decimal digits.

BTW, your variable names are horrible! In fairness, your code is typical (unfortunately) of people studying physics who have had no formal computer science classes, and who have never been exposed to good coding style. Fortran no longer has an 8-character limit on variable names, and variable names and commands can be in upper or lower cases. If you ever want anyone else to be able to understand your code, use variable names that are self explanatory. I would guess that RHOMN has something to do with $\rho$ (rho), and maybe VOLSED has something to do with volume, but I can't even guess what RWTSED is supposed to represent.

 

[Atr cheema]

Your program uses the REAL data type, which is four bytes in size. Precision is limited to about 7 decimal digits. For higher precision, use DOUBLE PRECISION (or REAL*8), which is good for about 17 decimal digits.

BTW, your variable names are horrible! In fairness, your code is typical (unfortunately) of people studying physics who have had no formal computer science classes, and who have never been exposed to good coding style. Fortran no longer has an 8-character limit on variable names, and variable names and commands can be in upper or lower cases. If you ever want anyone else to be able to understand your code, use variable names that are self explanatory. I would guess that RHOMN has something to do with $\rho$ (rho), and maybe VOLSED has something to do with volume, but I can't even guess what RWTSED is supposed to represent.

Thank you for the reply. I have understood that the problem is related to precision to which real number is represented by floating point number. I changed the variable definition from REAL to DOUBLE PRECISION, and I get the answer

      program inpdat
       IMPLICIT NONE

       DOUBLE PRECISION A,B,C,C1
2001 FORMAT (F40.10)
       A = 17424.0
       B = 2.42
       C = 17424.0*2.42*1.0E6
       C1 = A*B*1.0E6

       WRITE(*, 2000) C, A, B
       WRITE(*, 2000) C1

I get following output

42166083584.000000000000000, 17424.000000000000000, 2.4200000762933945
42166081329.3457

My question is how can I get the same answer from Fortran which I get from calculator/Python i.e. 42166083584.00000.

 

[Mark44]

My question is how can I get the same answer from Fortran which I get from calculator/Python i.e. 42166083584.00000.

Possibly this will help: Write 1.0E6 as 1.0D6. This difference is that 1.0E6 is a real*4 (REAL) while 1.0D6 is a real*8 (DOUBLE PRECISION).

Also, your format statement isn't being used. The label here is 2001, but your WRITE statements refer to a nonexistent label, 2000.
You should have two FORMAT statements, one for the first write of three values, and a different one for the second write.

You might try these:

2000  FORMAT (3F40.20)
2001  FORMAT (F40.20)

In the first format statement, the 3 indicates that three numbers will be printed, each in a field of total width of 40 characters, with 20 digits to the right of the decimal point.[/code]

 

[Atr cheema]

Possibly this will help: Write 1.0E6 as 1.0D6. This difference is that 1.0E6 is a real*4 (REAL) while 1.0D6 is a real*8 (DOUBLE PRECISION).

Also, your format statement isn't being used. The label here is 2001, but your WRITE statements refer to a nonexistent label, 2000.
You should have two FORMAT statements, one for the first write of three values, and a different one for the second write.

You might try these:

2000  FORMAT (3F40.20)
2001  FORMAT (F40.20)

In the first format statement, the 3 indicates that three numbers will be printed, each in a field of total width of 40 characters, with 20 digits to the right of the decimal point.[/code]

I applied your suggestions but still unable to get the result what I get from calculator..

      program inpdat
c

      IMPLICIT NONE
      DOUBLE PRECISION      A,B,C,C1,C2
C      REAL A,B,C,C1
2000 FORMAT (' ',3F40.15)
2001 FORMAT (' ',F40.10)

      A = 17424.0
      B = 2.42   !0000076293945
      C= 17424.0*2.42*1.0D6
      C1 = A*B*1.0E6
      C2 = A*B*1.0D6


      WRITE(*,2000) C, A, B
      WRITE(*,2001) C1
      WRITE(*,2001) C2

      end

OUTPUT

42166082031.250000000000000               17424.000000000000000                       2.420000076293945
42166081329.3457031250
42166081329.3457031250

 

[FactChecker]

Your accuracy problem is coming from the 2.42 not being represented with enough accuracy. As an experiment, try setting B = 242.0 and change the 1.0D6 to 1.0D4.

Or I think you can just set B = 2.42D0 so B is set to the double precision constant.

 

[Hugh McCutchen]

Please verify what you believe to be the correct answer. 17424.0*2.42 = 42,166.08. The answer can have only two places after the decimal. Multiplying by 1.0E6 gives 42,166,080,000.00 which is the answer from PCalc. I think both your calculator and Python are giving you the wrong answer.

 

[FactChecker]

The printout of B in post #7 shows that there are garbage digits in B. I verified that those garbage digits cause the error in the answer given. The assignment of a single precision constant to the double precision B is the cause of the problem. I suggest using B = 2.42D0

 

[Mark44]

Your accuracy problem is coming from the 2.42 not being represented with enough accuracy.

I agree, and this could be the fault of the compiler, gfortran.

I wrote essentially the same code in C, and got much better results.

// Test.c

#include <stdio.h>

int main(void)

{
   double a = 17424.0, b = 2.42, c;

   c = a * b * 1.0e6;
   printf(" c = %0.18f\n a = %0.18f\n b = %0.18f\n", c, a, b);
}

Output
c = 42166080000.000000000000000000
a = 17424.000000000000000000
b = 2.419999999999999929

My value for b is off in the 17th decimal place, and is much closer to 2.42 than the value in the OP's fortran program, 2.420000076293945, which is off in the 8th decimal place. It's as if floating point literals are single precision by default, and I am not able to find a workaround. This lack of precision in the GNU fortran compiler is very disappointing to me.

Multiplying by 1.0E6 gives 42,166,080,000.00 which is the answer from PCalc. I think both your calculator and Python are giving you the wrong answer.

My C code produced the correct value. Note that @Atr cheema is using Fortran, not Python.

 

[Atr cheema]

Thank you @FactChecker and @Hugh McCutchen . The correct answer is indeed obtained by defining B=2.42D0. I want to ask now that this code snippet is part of large code and I wanted to check value of C (RWTSED in first post), while A (VOSED in post #1) and B (RHOMN in post #1) are already calculated by program. When I print C in original program it prints wrong value but I am sure internally it is considering right value of C. How can I print true value of C in case when A, and B are already being calculated by program?

 

[FactChecker]

In FORTRAN, a double precision variable, B, that is initialized in a declaration "double b=2.42" will have the correct double precision value. In the declaration section, conversion of single precision numbers to double is automatically done when needed. But in the executed part, care must be taken to represent double precision numbers as double: 2.42D0. Otherwise, there will be garbage lower-precision digits. The code "C= 17424.0*2.42*1.0D6" is wrong. It should be "C= 17424.0*2.42D0*1.0D6". And "C= 17424.0D0*2.42D0*1.0D6" would be preferred, even though it is only the fractional part of numbers that causes the problem.

 

[jbriggs444]

      VOLSED = 17424.0
      RHOMN = 2.42   !0000076293945
      RWTSED= VOLSED*RHOMN*1.0E6

If you factor that 1.0E6 and into 1.0E3 * 1.0E3 and use it to scale the values of VOLSED and RHOMN to 17424000 and 2420 then they will be integers. Their product will be an integer. If you promote VOLSED and RHOMN to double precision both of those integers and the product thereof will be exactly expressible in double precision floating point. Fortran will be able to compute and express the product exactly. No error, not even in the 17th digit.