Floating Point Differences Between C and F90

When Max Smirnoff and Matt Rundquist were translating ToyFDTD1 into F90, they found that the floating point outputs were diverging from those of the C version. Max traced this effect to a difference in the assembly; in the end, they found this to be due to a difference in the language standards. Max wrote this page to illustrate exactly what happens.

To see the differences between C and F90 floating point, consider a sample program, and its C equivalent.  These programs simply assign 0.7 to a number and than print it on screen.

F90:

program foo
real*8 :: a = 0.7
print*,a
end program foo

C:

#include <stdio.h>

void main(void)
{
double a = 0.7;
printf(" %.10f \n",a);
}

C output :  0.70000000000
F90 output : 0.69999998807907104

Upon examination of the assembly machine code produced by F90, it appears to ignore that a is declared as an 8 byte real, and instead stores only the first 4 bytes of 0.7 in IEEE format.

C assembly:

.dword 0x3fe6666666666666 # double 0.700000

F90 assembly:

.dword 0x3fe6666660000000 # double 0.700000

The solution to that problem is to store every number by using a "d" instead of an "e" to denote its exponent. For example:

store 0.7 as 0.7d0, or 3.455e-5 as 3.455d-5.

This improves accuracy, sometimes significantly. The example below stores light speed squared in the old and new notation. Notice the improvement in precision:

old notation : a =   89875517873681764.0    prints as     89875514973487104
new notation : a =  89875517873681764.0d0   prints as  89875517873681760

The minor discrepancy might not seem significant, and for a few calculations, it probably isn't, but notice what happens when we do many calculations, where the errors can accumulate. The program below simply starts with 0.0, adds 0.01 to the variable a 100 times, and prints out the value of the variable a at each iteration. The output from both C and F90 code is shown for comparison.

F90:

program foo
real*8 :: a = 0
integer :: n
do n=0,100,1
print*,a
a = a + 0.01
end do
end program foo

C:

#include <stdio.h>
void main(void)
{
double a = 0;
int n;
for(n=0; n<=100; n++)
{
printf(" %.10f \n",a);
a += 0.01;
}
}

 F90 output 0.E+0  9.99999977648258209E-3  1.99999995529651642E-2  2.99999993294477463E-2  3.99999991059303284E-2  4.99999988824129105E-2  5.99999986588954926E-2  6.99999984353780746E-2  7.99999982118606567E-2  8.99999979883432388E-2  9.99999977648258209E-2  0.1099999975413084  0.11999999731779099  0.12999999709427357  0.13999999687075615  0.14999999664723873  0.15999999642372131  0.1699999962002039  0.17999999597668648  0.18999999575316906  0.19999999552965164  0.20999999530613422  0.21999999508261681           .....    .....    .....   0.80999998189508915  0.81999998167157173  0.82999998144805431  0.8399999812245369  0.84999998100101948  0.85999998077750206  0.86999998055398464  0.87999998033046722  0.88999998010694981  0.89999997988343239  0.90999997965991497  0.91999997943639755  0.92999997921288013  0.93999997898936272  0.9499999787658453  0.95999997854232788  0.96999997831881046  0.97999997809529305  0.98999997787177563  0.99999997764825821 C output  0.0000000000   0.0100000000   0.0200000000   0.0300000000   0.0400000000   0.0500000000   0.0600000000   0.0700000000   0.0800000000   0.0900000000   0.1000000000   0.1100000000   0.1200000000   0.1300000000   0.1400000000   0.1500000000   0.1600000000   0.1700000000   0.1800000000   0.1900000000   0.2000000000   0.2100000000   0.2200000000     ....        ....    ....  0.8100000000   0.8200000000   0.8300000000   0.8400000000   0.8500000000   0.8600000000   0.8700000000   0.8800000000   0.8900000000   0.9000000000   0.9100000000   0.9200000000   0.9300000000   0.9400000000   0.9500000000   0.9600000000   0.9700000000   0.9800000000   0.9900000000   1.0000000000

From this partial listing of the outputs, you can see that errors do accumulate, and over many iterations, such as in our simulations, the results can diverge. If we modify the F90 code so that 0.01d0 is added to a instead of 0.01, the precision will get better, but 1.0 will still be printed as 0.9999... by F90. This problem might not exist on other compiler versions, but you should be aware of it.

-- Max Smirnoff

 Contacting the perpetrators: lemiller@borg.umn.edu