具有iso_c_binding的fortran内部exp函数的行为

Question

为了检查fortran中iso_c_binding的一些属性,我写了一个小程序

program double
    use iso_c_binding, only : C_DOUBLE, C_INT32_T, C_FLOAT
    implicit none
    integer, parameter :: rk  = C_DOUBLE
    integer, parameter :: rk4 = C_FLOAT
    integer, parameter :: ik  = C_INT32_T
    integer(kind=ik) :: i, j
    integer(kind=ik), parameter :: n = 10
    real(kind=rk), dimension (n) :: array, narray, oarray
    real(kind=rk) :: barray(n) = (/(i, i=1,n,1)/)

    do i = 1, n
      oarray(i) = 1.d0 * real(i)
      array(i) = exp(real(i))
    end do
    narray(:) = exp(oarray(:))
    barray(:) = exp(barray(:))
    do i = 1, n
      write(*,*), array(i), narray(i), barray(i)
    end do

  end program double

使用gfortran(gfortran -o prog -Wall -I.-lm build/prog.o)进行编译并运行会得到以下结果:

2.7182817459106445  2.7182818284590451  2.7182818284590451    
7.3890562057495117  7.3890560989306504  7.3890560989306504     
20.085536956787109  20.085536923187668  20.085536923187668     
54.598148345947266  54.598150033144236  54.598150033144236     
148.41316223144531  148.41315910257660  148.41315910257660     
403.42880249023438  403.42879349273511  403.42879349273511     
1096.6331787109375  1096.6331584284585  1096.6331584284585     
2980.9580078125000  2980.9579870417283  2980.9579870417283     
8103.0839843750000  8103.0839275753842  8103.0839275753842     
22026.464843750000  22026.465794806718  22026.465794806718

似乎即使值'array(1)'也有一些错误,因为它已经偏离了e https://en.wikipedia.org/wiki/E_(mathematical_constant)

这种行为来自哪里？

Answer 1

问题是你的计算array,

array(i) = exp(real(i))

以单精度完成.数组本身是双精度的,但表达式real(i)是单精度,因为默认的实数类型是单精度,因此,exp在调用它时,您只能期望7位数的精度.其他数组以完全双精度计算.

尝试更换exp(real(i))与exp(real(i,rk))和差异就会消失.