Archive for April, 2007

damn numerical error in Fortran

April 21st, 2007 2 comments
if we run the following simple code:
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;
}
}

we may find different outputs from different compilers:
Intel Fortran g95 Gcc
0.909999979659915
0.919999979436398
0.929999979212880
0.939999978989363
0.949999978765845
0.959999978542328
0.969999978318810
0.979999978095293
0.989999977871776
0.999999977648258
0.909999979659915
0.9199999794363976
0.9299999792128801
0.9399999789893627
0.9499999787658453
0.9599999785423279
0.9699999783188105
0.979999978095293
0.9899999778717756
0.9999999776482582
0.9100000000
0.9200000000
0.9300000000
0.9400000000
0.9500000000
0.9600000000
0.9700000000
0.9800000000
0.9900000000
1.0000000000
if we change that +0.01 to +0.01d0, it will improve the accuracy but still unsatisfactory
Ifort real*16 Ifort real*8 G95 real*8
0.910000000000000018943180357666733
0.920000000000000019151347174783950 0.930000000000000019359513991901167 0.940000000000000019567680809018384 0.950000000000000019775847626135601 0.960000000000000019984014443252818 0.970000000000000020192181260370035 0.980000000000000020400348077487251 0.990000000000000020608514894604468 1.00000000000000002081668171172169
0.910000000000001
0.920000000000001
0.930000000000001
0.940000000000001
0.950000000000001
0.960000000000001
0.970000000000001
0.980000000000001
0.990000000000001
1.00000000000000
0.9100000000000006
0.9200000000000006
0.9300000000000006
0.9400000000000006
0.9500000000000006
0.9600000000000006
0.9700000000000006
0.9800000000000006
0.9900000000000007
1.0000000000000007
so.. god bless fortran…
Categories: Technical Tags: