2007年4月 的存档

damn numerical error in Fortran

2007年4月21日 2 条评论
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…
分类: 技术相关 标签: