Extremely pretty logo by J-D Harrington
 
 

Floating Point Differences Between C and F90


Return to ToyFDTD home



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





Am I the only sane dot on this page? Contacting the perpetrators: lemiller@borg.umn.edu


ToyFDTD home