2

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?

francescalus
  • 30,576
  • 16
  • 61
  • 96
user1887919
  • 829
  • 2
  • 9
  • 24
  • You can see what gfortran is doing with the -fdump-tree-original option. For the integer exponent, a call to _gfortran_pow_r16_i4() is generated. This is a part of gfortran runtime library. Its algorithm comes from DE Knuth, The Art of Computer Programming, Vol. 2. For the floating point exponent, a call to powq() is generated. powq() is a part of libquadmath, so you can go read how an integral valued exponent is handled. – evets Jul 10 '19 at 15:00

2 Answers2

5

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"

L. Scott Johnson
  • 4,213
  • 2
  • 17
  • 28
Ian Bush
  • 6,996
  • 1
  • 21
  • 27
  • It may be worth noting that exponentiation by an integer power is an exception for "numeric operands are converted to the type and kind of the result before evaluating the expression". That is, it's allowed to see a difference between `x**4` and `x**4.` because we don't need to convert the former to the latter. – francescalus Jul 10 '19 at 11:07
  • Hmmmmm interesting! So in general the best practice for exponentiation is to use integers since it is quicker and has the same accuracy? Also, it is somewhat unsatisfying that we have to accept both answers as equally accurate - clearly one of them will be closer to the true value! – user1887919 Jul 10 '19 at 12:10
  • @user1887919 For any given calculation, one of them will give a result closer to the IEEE 754 round-to-nearest of the exact answer. There may not be a general rule that says which will be more accurate in general. – Patricia Shanahan Jul 10 '19 at 13:04
1

For the sake of completeness, exponentiating the (exact) fraction leads to

46269923.019114341045212564354844220930226304938209797955723262974801
46269923.0191143410452125643548442211 is the nearest
46269923.0191143410452125643548442147

Exponentiating the nearest quadruple precision floating point (https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format) leads to

46269923.01911434104521256435484422112355946434320203876355837024902939447201788425445556640625
46269923.0191143410452125643548442211 does fit

So it looks like the exponentiation function did round to the nearest float. I don't know if it is luck or if the underlying math library has guaranteed correct rounding.

Integer exponentation:

x2 = x*x
x4 = x2*x2

cumulates two round off error, so it might be 1 ulp off, it's not unsurprising.

aka.nice
  • 9,100
  • 1
  • 28
  • 40