Fourmilog: None Dare Call It Reason

« Yotta Yotta Universe | Main | Reading List: La Force de la Raison »

Tuesday, November 30, 2004

FORTRAN Double Precision Constant Gotcha

What do you think this FORTRAN program prints, when compiled with g77 and run?

        DOUBLE PRECISION X
	
        X = 0.0025

        PRINT 1000, X
1000    FORMAT (1X, F14.12)

        END

If you guessed "  0.002500000000", you've fallen into the same trap I did with the FORTRAN version of my floating point benchmark. In fact, this program prints "  0.002499999944". Why? Because the statement "X = 0.0025" assigns a REAL constant of 0.0025, which has only 32 bits of precision, to the DOUBLE PRECISION variable X, resulting in the round-off error when the value is printed to 12 decimal places. To cause the constant assigned to X to be treated as DOUBLE PRECISION, you must write "X = 0.0025D0", where the "D" exponent denotes a DOUBLE PRECISION constant. Earlier FORTRAN compilers with which the floating point benchmark was tested (for example SGI MIPSpro FORTRAN version 7.2.1) appear to promote decimal constants to DOUBLE PRECISION when they appear in a DOUBLE PRECISION expression or are used to initialise a DOUBLE PRECISION variable, but the FORTRAN standard does not specify this behaviour and the g77/gcc compiler strictly follows the standard and requires the "D" in all DOUBLE PRECISION constants.

Posted at November 30, 2004 20:48