x**n and pow(x,n) for integer n

[David G. Hough on validgh]

A bug report against Fortran 1.4 made me look into this again.  I had
figured it out once but it had become unfigured subsequently.  The
issue is that routines for x**n are sometimes inefficient and sometimes
inaccurate, particularly for the case of IEEE arithmetic: consider
1.0d155**-2 and 10.0**-30; these should not be computed by
reciprocating first and then powering the negated exponent, because of
gratuitous overflow in the first case and roundoff accumulation in the
second case.

Note that if you can perform the computation in higher precision and
exponent range then none of these complications arise.  Assume you have
to compute entirely in working precision.

In Fortran, use of x**n rather than x**y in source code might be a
conscious coding decision that suggests but does not prove that
computation by repeated multiplication rather than by exp and log was
intended.  In traditional C pow(x, i) suggests a coding error, while in
ANSI-C use of an integer for the exponent may be more or less
accidental.  If you want x**y to be monotonic and the general pow()
implementation is not correctly rounded, you might not want to treat
pow(x,{-1,1,2}) or x**{-1.0,1.0,2.0} as special cases.

However if you want to treat special cases in the front end, converting
	x**-1	to 1/x
	x**0	to 1
	x**1	to x
	x**2	to square(x)
	x**3	to x*square(x)
	x**4	to square(square(x))
seem worthwhile, although you lose out on x**0 if x is a signaling NaN.
What isn't OK is converting x**n to (1/x)**-n or 1/(x**-n),
or x**-n to (1/x)**n or 1/(x**n) for variable n, since you may have undone
somebody's careful plotting against overflow or underflow.

In the run-time library (pow_di in ATT f77 progeny) the code should do 
something like this.  The following avoids checking underflow or overflow
flags by doing scaling instead.  You'd use something like frexp and ldexp
instead, "something like" rather than the real thing to avoid gratuitous
entanglements with matherr() or errno.

Another consideration is whether for sufficiently large |n| you wouldn't
be better off in efficiency and performance to call x**y, but consider
monotonicity and how often this is really going to happen.  Shouldn't the 
programmer take responsibility for this decision?

pow_di(x,n):

if n>0
	s=x
	k=n
else
	if n=0
		signal invalid if x a signaling nan, e.g. by x+0,
		but look out for optimizers
		return 1
	k=-n
	s="frexp"(x,&e)

p=x
while even(k)
	p=square(p)
	k=k/2
t=p
k=k/2

while k>0
	p=square(p)
	if odd(k)
		t=t*p
	k=k/2

if n>0
	return t
else
	return "ldexp"(1/t, e*k)


Any errors are exercises for the student.