Path: utzoo!mnetor!uunet!seismo!sundc!pitstop!sun!quintus!ok
From: ok@quintus.UUCP (Richard A. O'Keefe)
Newsgroups: comp.lang.c
Subject: Re: pow again (Complex valued or partially defined?)
Message-ID: <530@cresswell.quintus.UUCP>
Date: 14 Jan 88 02:27:19 GMT
References: <629@PT.CS.CMU.EDU> <5825@sol.ARPA> <7071@brl-smoke.ARPA>
Organization: Quintus Computer Systems, Mountain View, CA
Lines: 85
Summary: No single meaning for complex**complex

In article <7071@brl-smoke.ARPA>, gwyn@brl-smoke.ARPA (Doug Gwyn ) writes:
> In article <5825@sol.ARPA>, quiroz@cs.rochester.edu (Cesar Quiroz) writes:
> >       There is no such thing as a unique, universally accepted
> >       definition for the exponentiation.
> 
> I don't buy this.
> 
> I was under the impression that in the complex number field, which is the
> natural extension of the reals and therefore also of the integers, there
> was a single, well-defined (not necessarily well-behaved!), universally
> agreed-upon mathematical meaning for "w raised to the power z".  If this
> isn't so, please elaborate.
> 
Sorry, but it *isn't* so.  Any good book on complex analysis will explain
why.  I'll try to sketch the basic problem.

We can represent any complex number (X + i*Y) in polar co-ordinates
as (R, Theta) for some real R >= 0 and real Theta.  Any value 2*pi*k+Theta
will do, for any integer k.  Now the nth root of a complex number is most
easily taken in polar form:
	(R, Theta) ** (1/n) = (R**(1/n), Theta/n + (2*pi*k)/n)
which means that every complex number other than 0 has n distinct nth roots.

If we want a single-valued function, we have to pick which of these roots
we want.  It turns out that this cannot be done without introducing a
discontinuity (look up "branch cut" in any good text).  In the case of nth
root, what we get is a ray (R,Theta) for fixed Theta and R > 0 where the
nth root function is continuous from one side only; we can place that ray
at any Theta we like.  Different choices of branch cut are useful in
different cases.  Common Lisp spells out where the branch cuts are for all
the built-in "elementary functions", which is one of the reasons why I
mentioned Common Lisp in an earlier message.  But if you want an nth root
function on complex numbers, not only is there not a single well-defined
meaning, there are uncountably many of them!

When we restrict ourselves to the reals, a real number may have 0, 1, or
2 real nth roots.  Only when we restrict ourselves to non-negative reals
does "nth root" become a function.

(R, Theta) ** n = (R ** n, Theta*n) for integer n _is_ well-defined.
So is the exponential function exp, which is defined thus:
	exp(X + i*Y) = (e ** X)*(cos Y + i*sin Y)
As you might expect, the break-down of pow() can be traced to the
break-down of log(): if you try to find an inverse of exp() you get
	log((R, Theta)) = log(R) + i*(Theta + 2*pi*k)
which has infinitely many values.  The usual choice of branch cut for
logarithm (yielding the so-called "principal value" of the logarithm)
restricts Theta to -pi < Theta < +pi, so that the logarithm of negative
numbers is not defined.  There is no reason at all why we couldn't
put the branch cut somewhere else, say -pi/2 < Theta < 3pi/2, which
*would* give (-1 + 0i) a logarithm, but not (0 + 1i).

Thus if one tried to define pow(z1, z2) as exp("log"(z1)*z2),
"log" being the principal value, pow(-27.0, 1.0/3.0) would not be
defined, which would be rather odd, as 4.3BSD and SunOS have a
cbrt() function, and cbrt(-27.0) is perfectly well defined (-3.0).
{This is why cbrt() exists, of course.}

double rtp(x, n)		/* this is the Raise To Power */
    double x;			/* function.  The binary method */
    int n;			/* is used.  Comments welcome. */
    {
	double a;

	if (n < 0) x = 1.0/x, n = -n;
	if (n == 0) return 1.0;
	while ((n&1) == 0) x *= x, n >>= 1;
	for (a = x; --n > 0; a *= x)
	    while ((n&1) == 0) x *= x, n >>= 1;
	return a;
    }

I tried to measure the accuracy of this, but it turned out that the
errors (~4e-16) were in my test data.  You might lose about one
decimal digit in relative error.  It is possible in principle to get
1/2ulp accuracy for rtp(), but I don't know how to do that in C.

C programmers on UNIX systems which also support F77 may be able to
get at Fortran's (double precision)**(integer) function by doing

	extern double pow_di(/* double, int */);

and then linking with		-lF77 -lm
Having rtp() {perhaps called pow_di, I wouldn't mind} in the C standard
would be much cleaner.