2686阅读
5回复

[讨论]Floating Point Differences Between C and F90

离线febi

UID ：2537

注册：2007-05-11
登录：2016-05-05
发帖：123
等级：退休版主

0楼发表于: 2008-03-25 07:09:45

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