Consistent floating point arithmetic in fortran with different compilers on different computers?

Question

How can I get the same results with floating-point arithmetic using different compilers and possibly also different computers?

This is p.f90

program fl
implicit none
real(kind=8) :: a

a=9d0/10d0
write(*,*) a

end program

With gfortran -o p p.f90 I get 0.90000000000000002

With ifort -g -o p p.f90 I get 0.900000000000000

Is there a way to achieve consistency? The small things propagate and become larger.

Best regards Alessandro

Answer 1

We'll turn to floating-point arithmetic below, but let me deal first with a misconception you have. Your use of list-directed formatting in the write statement means that the format chosen by the compiler to write out a variable is, well, the compiler's choice; it's not mandated by the language standard. Your use of the second * in write(*,*) tells the compiler to write out the value of a variable as it wishes. So what you have is not evidence that there is something wrong with the arithmetic, but that there are differences between gfortran and ifort. If I modify your write statement to

write(*,'(f21.18)') a

then my Intel Fortran program writes

0.900000000000000022

to the console.

SO is littered with questions arising from unfamiliarity with the details of floating-point arithmetic so I'm not going to write a treatise, just a few observations relevant to your question.

IEEE-754 64-bit floating-point numbers (which is probably what you get with the declaration real(kind=8)) only provide approximately 16 decimal digits of useful information. Actually, since they're binary and there is no easy correspondence between the number of digits in one base and in another, it's actually 15.95 decimal digits and many users round this down and never look at anything after the 15th significant figure in a floating-point number's decimal representation. So both ifort and gfortran are misleading you with their trailing 2s.

IEEE-754 defines not only the formats for f-p numbers, but also some rules for rounding and for arithmetic operations. A carefully-written program which uses only those arithmetic operations (I think square root is also specified) and which takes care about rounding modes and rounding operations should produce the same results on two different processors. Of course, not many useful numerical programs confine themselves to the basic arithmetic operations.

Since the 2003 standard Fortran has included an intrinsic module called ieee_arithmetic which gives the programmer direct access to the underlying IEEE-754 capabilities of the hardware on which it is running -- but note that it does not require that the hardware have any such capabilities. Using ieee_arithmetic and another intrinsic module called ieee_exceptions you should, if your hardware provides the necessary support, be able to write programs which compile under both gfortran and ifort and which produce, when executed, the same results to the last bit of every f-p number.

You'll also need to familiarise yourself with the optimisation options and numerical arithmetic options of the compiler(s) you're using. Most compilers have an option whose meaning is sacrifice IEEE compliance for speed (for ifort I think it's fp-model). In general adherence to IEEE-754 for operations will slow your programs.

Answer 2

You cannot.

See under Systems Aspects in What Every Computer Scientist Should Know About Floating-Point Arithmetic.