Fortran x**n vs. x**y

[David G. Hough at validgh]

Like many other, but not all, Fortran compilers, Sun's evaluates x**n by
O(logn) multiplications, and x**y using exp and log.   Of course the
results vary between x**n and x**dble(n).

Sun got a bug report on this, to which I responded as follows.   I'm interested
in comments from this group.

> Not a bug.
> 
> The underlying problem is expressed in Fortran (not in C since x**n has no
> analog there) is that
> 
>         double precision x, y
>         integer n
> 
>         x= 1.0000003205183601d0
>         y= -2147483648.0d0
>         n=y
>         print *,x,y,' x**y ',x**y
>         print *,x,n,' x**n ',x**n
>         print *,x-x,y-n,' diff ',x**y-x**n
>         end
> 
> produces
> 
>      1.0000003205184   -2147483648.0000 x**y     1.1795298980653-299
>      1.0000003205184  -2147483648 x**n     1.1795296972117-299
>    0.  0. diff     2.0085364200187-306
> 
> 
> This is pretty much an inevitable difference in the way x**y and x**n are
> customarily evaluated, by log+exp and by repeated multiplication.   x**n is a poor choice for large |n|, as it's usually
> evaluated by repeated multiplication; x**y should be used instead, in the
> program source code.   The following expanded example demonstrates the
> deteriorating accuracy due to the rounding accompanying each multiplication:
> 
>         double precision x, y
>         integer n
> 
>         x= 1.0000003205183601d0
>         n=-1
>         do 1 i = 1, 32
>         y=n
>         print *,x,y,' x**y ',x**y
>         print *,x,n,' x**n ',x**n
>         print *,x-x,y-n,' rel diff ',(x**y-x**n)/x**y
>         n=2*n
>  1      continue
>         end
> 
> 
> This should be clearly documented in the Sun Numerical Computation Guide,
> and the reasons explained.
> 
> The only alternative is to implement x**n in a more complicated way, switching
> to x**y if appropriate based on a-priori error bounds, or when
> more efficient, but that
> defeats the purpose of x**n.
> 
> I suspect the original problem involves interest rates or other exponential
> growth computations.   There the appropriate solution is to use
> compound(3m) or annuity(3m) or their Fortran analogs.
> 
> Note that if there is any uncertainty at all in x close to 1, when raised to
> a large |n|, 
> the relative uncertainty in x**n or x**y will be much greater than
> the relative uncertainty in x.   So even the x**n algorithm provides a 
> satisfactory backward error analysis.    
> But as an extra bonus, we provide x**y which uses
> implicit higher precision to provide a satisfactory forward error analysis,
> which is as good as we can hope for, short of correct rounding.   It's slower,
> of course.
> 
> The foregoing addresses the differences between x**n and x**y as expressed in
> source code.   Just in case the issue is the compiler front end converting
> x**y to a fixed series of multiplications and divisions
> when y is known at compile time to be a fixed integer K, 
> we have a more serious problem.    That optimization is safe
> enough when y is known at compile time to be -1, 0, 1, or 2, 
> and perhaps tolerable for 3 and 4, but should probably be foregone otherwise.