Exponentiation in Fortran/gfortran to high precision

Question

How does gfortran handle exponentiation with a integer vs a real? I always assumed it was the same, but consider the example:

program main

implicit none

integer, parameter :: dp = selected_real_kind(33,4931)

real(kind=dp) :: x = 82.4754500815524510_dp

print *, x
print *, x**4
print *, x**4.0_dp

end program main

Compiling with gfortran gives

82.4754500815524510000000000000000003      
46269923.0191143410452125643548442147      
46269923.0191143410452125643548442211 

Now clearly these numbers almost agree - but if gfortran handles integers and reals for exponentiation in the same way I would expect them to be identical. What gives?

Answer 1

Expanding your program slightly shows what is going on:

ijb@ianbushdesktop ~/work/stack $ cat exp.f90
program main

implicit none

integer, parameter :: dp = selected_real_kind(33,4931)

real(kind=dp) :: x = 82.4754500815524510_dp

print *, x
print *, x**4
print *,(x*x)*(x*x)
print *,Nearest((x*x)*(x*x),+1.0_dp)
print *, x**4.0_dp

end program main

Compiling and running gives:

ijb@ianbushdesktop ~/work/stack $ gfortran --version
GNU Fortran (GCC) 7.4.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ijb@ianbushdesktop ~/work/stack $ gfortran exp.f90 
ijb@ianbushdesktop ~/work/stack $ ./a.out
   82.4754500815524510000000000000000003      
   46269923.0191143410452125643548442147      
   46269923.0191143410452125643548442147      
   46269923.0191143410452125643548442211      
   46269923.0191143410452125643548442211      
ijb@ianbushdesktop ~/work/stack $ 

Thus

1. It looks like the compiler is clever enough to work out that integer powers can be done by multiplies, which is much quicker than a general exponentiation function

2. Using the general exponentiation function is only 1 bit different from the multiplication answer. Given we can't say per se which is more accurate we must accepts both as equally accurate

Thus in conclusion the compiler uses simple multiplies where possible instead of a full blown exponentiation routine, but even if it has to use the more expensive route it gives the same answer, for carefully considered meaning of the word "same"