Simple math error in gfortran?
I am starting to program in Fortran. I am currently using cygwin and the gfortran compiler running on my Windows XP PC. I'm having trouble with some simple math - the program I wrote just won't do the math. Code:
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?
source to share
You are working with a single limit of precision. The result of the subtraction is real*4
, and can be safely stored in t
. The values used in the subtraction, however, are out of range real*4
. Outside in the precision range, the resulting number is rounded up 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 attached one _8
ensures that the two numbers are double precision; the result is then converted to real before being assigned t
, but the calculation is now in double precision and .25
"stored" in subtraction.
source to share
Here's an example using the High Performance Mark proposal. Type "8" is a non-portable way of specifying double repetitions. Most compilers now support the ISO Fortran environment types used in this example. If you have an older compiler, you can instead use a ISO_C_BINDING
type c_double
that has been available for longer. As Evert did, you need 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. This is not enough for the variable on the LHS to 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
source to share
While it doesn't matter for the example you gave, if there were more inconsistencies in the values of the operands, replace
real t
from
double precision t
Otherwise, the difference will be truncated into available storage.
For your example, the precision of the constants is too low. A floating point value, the both have ( 8320671.25
and 8000000.00
) has a type real
that indicates a "single precision": a 32-bit representation of floating point values, which are equivalent to 6 to 7 of decimal digits of accuracy . Double precision
(since the 1980s) is usually a 64-bit floating point representation equivalent to 15 to 17 decimal digits . However, some implementations may implement Double precision
the same real
, although this is usually larger on IEEE-754 compliant machines.
To make a floating point constant double, write D
for exponent:
t = 8320671.25D0 - 8000000.00
Note that execution of subsequent operands Double precision
will not work due to the Fortrans typing rules:
t = 8320671.25 - 8000000.00D0 ! Wrong! Result is "real"
or
t = 8320671.25 - 8D6 ! Wrong! Result is "real"
The first example prints the output 320671.250
. The rest produce 320671.000
(if t
- type real
) or 320671.00000000000
, if t
- type Double precision
.
Fortran has been implemented on a variety of architectures since the 1940s. On some of them (like VAX), command line compilation options dictate what type of hardware is used for real
and for Double precision
. The DEC architecture had at least four floating point flavors (F-Float, D-Float, G-Float, and H-Float), all of which were developed well before IEEE-754. See this for some historical perspective on the development of IEEE-754.
source to share