Rounding error when converting Double to Real in fortran

Question

i'm trying to convert a double into a real for some operations but i noticed an error. i have a macro for reading purposes, short is 4 bytes, and long is 8 bytes.

this is the execution code :

print*, x, real(x, kind=short) 
print*, x*2.0, real(x*2.0, kind=short) 

print*, real( real(x, kind=short), kind=long) 
print*, real( real(x*2.0, kind=short), kind=long)

this prints out :

173.43304556187957 173.433044 
346.86609112375913 346.866089

173.43304443359375 
346.86608886718750 

so first of all why do i even have any values after converting to real then back again? Shouldn't it be something like 173.43304400000000? And why are my values changing when i convert them? How can i change that ?

For info, fortran 2003, compiled with gfortran but i think we had the same problem with ifort.

Thank you!

Answer 1

No, it shouldn't. The values are stored as binary (base 2) floating point, and the extra fraction bits for double precision are binary zeroes, not decimal zeroes. Most decimal fractions don't convert exactly to binary (and vice-versa).

You can better see what is going on when you display the values in hexadecimal format. I used this program:

integer, parameter :: short = kind(0.0)
integer, parameter :: long = kind(0.0D0)
real(short) :: r4
real(long) :: r8

r8 = 173.43304556187957_long
print '(G0.16,1X,Z16.16)', r8,r8
r4 = r8
print '(G0.16,1X,Z8.8)',r4, r4
r8 = r4
print '(G0.16,1X,Z16.16)', r8,r8
end

and when I run it, I get:

173.4330455618796 4065ADDB825DBE6C
173.4330444335938 432D6EDC
173.4330444335938 4065ADDB80000000

Note that the reconversion to "long" just adds binary zeroes, not changing the value. By the way, using list-directed formatting can obscure details such as this.