CODEWITHSUNDEEP
fortrangfortranfortran-iso-c-binding
pangeorg
pangeorg

Reputation: 3

Behaviour of intrinsic exp function in fortran with iso_c_binding

In order to check some properties of the iso_c_binding in fortran, I wrote a little program

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

Compiling with gfortran (gfortran -o prog -Wall -I. -lm build/prog.o) and running gives the following result:

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

It seems that even the value 'array(1)' has some error in it as it already deviates from e https://en.wikipedia.org/wiki/E_(mathematical_constant)

Where does this behaviour come from?

Upvotes: 0

Views: 161

Answers (1)

amaurea
amaurea

Reputation: 5067

The problem is that your calculation for array,

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

is done in single precision. The array itself is double precision, but the expression real(i) is single precision because the default real type is single precision, and hence, you can only expect 7 digits of accuracy from exp when called on it. The other arrays are computed in full double precision.

Try replacing exp(real(i)) with exp(real(i,rk)), and the discrepancy should go away.

Upvotes: 6

Related Questions