CODEWITHSUNDEEP
floating-pointfortrangfortran
Dexter Ferreira
Dexter Ferreira

Reputation: 31

Simple Math Error in gfortran?

I am a beginner with regards to programming in Fortran. I am currently using cygwin and the gfortran compiler running on my Windows XP PC. I'm having difficulty with some simple math - the program I wrote simply won't do the math. The code is:

program convert
real t

t=0
t=8320671.25 - 8000000.00

write(*,*) t
end

The program should give me the answer "320671.25" but gives me 320671.00 instead! What am I doing wrong?

Upvotes: 3

Views: 427

Answers (3)

wallyk
wallyk

Reputation: 57794

Though it does not matter for the example you gave, if there were a greater disparity in the magnitudes of the operands, replace

real t

with

double precision t

Otherwise, the difference would be truncated to the available storage.

For your example, the precision of the constants is too small. A floating point value given as you have (8320671.25 and 8000000.00) is of type real which specifies "single precision": that uses a 32 bit floating point value representation which are equivalent to 6 to 7 decimal digits of precision. Double precision is (since the 1980s) usually a 64 bit floating point representation, equivalent to 15 to 17 decimal digits. However some implementations might choose to implement double precision the same as a real, though on IEEE-754 compliant machines, it usually is larger.

To make a floating point constant double, write a D for the exponent:

t = 8320671.25D0 - 8000000.00

Note that making subsequent operands double precision won't work due to Fortran's expression typing rules:

t = 8320671.25 - 8000000.00D0   ! Wrong!  Result is "real"

or

t = 8320671.25 - 8D6            ! Wrong!  Result is "real"

The first example produces the output 320671.250. The others produce 320671.000 (if t is type real) or 320671.00000000000 if t is type double precision.

Since the 1940s, Fortran has been implemented on diverse architectures. On some of them (VAX for example), command line compilation options specify what hardware type is used for real and for double precision. DEC's architecture had at least four flavors of floating point (F-Float, D-Float, G-Float, and H-Float), all of which were designed long before IEEE-754. See this for some historical perspective of the development of IEEE-754.

Upvotes: 0

M. S. B.
M. S. B.

Reputation: 29401

Here is a worked example using High Performance Mark's suggestion. The type "8" is not a portable way to specify double precision reals. Most compilers now support the ISO Fortran environment types used in the example. If you have an older compiler, you can instead use the the ISO_C_BINDING and the type c_double, which have been available longer. As Evert did, it's essential to specify the types of the constants. The calculation is done on the RHS, then the assignment is made to the variable on the LHS. It's not enough that the variable on the LHS have sufficient precision.

program convert

use, intrinsic :: ISO_FORTRAN_ENV

real (real64) :: t
t=8320671.25_real64 - 8000000.00_real64
write(*,*) t

end program convert

Upvotes: 2

user707650
user707650

Reputation:

You're running into the single precision limit. The result of the subtraction is a real*4, and can be safely stored in t. The values used in the subtraction, however, are outside the real*4 range. Outside in the range of precision, with the result that the numbers are rounded to fit into a real*4 before calculating the subtraction.

Try this for example:

program convert
real t

t=0
t=8320671.25_8 - 8000000.00_8

write(*,*) t
end

The appended _8 ensures the two numbers are double precision; the result is then converted to real before assigned to t, but the calculation is now in double precision, and the .25 is "saved" in the subtraction.

Upvotes: 3

Related Questions