[Numpy-discussion] Handling of numpy.power(0, <something>)

[Timothy Hochberg]

On Wed, Feb 27, 2008 at 4:10 PM, Stuart Brorson <sdb at cloud9.net> wrote:

> I have been poking at the limits of NumPy's handling of powers of
> zero.   I find some results which are disturbing, at least to me.
> Here they are:


[SNIP]

>
> **  0^(x+y*i):  This one is tricky; please bear with me and I'll walk
> through the reason it should be NaN.
>
> In general, one can write a^(x+y*i) = (r exp(i*theta))^(x+y*i) where
> r, theta, x, and y are all reals.  Then, this expression can be
> rearranged as:
>
> (r^x) * (r^i*y) * exp(i*theta*(x+y*i))
>
> = (r^x) * (r^i*y) * exp(i*theta*x) * exp(-theta*y)
>
> Now consider what happens to each term if r = 0.
>
> -- r^x is either 0^<positive> = 1, or 0^<negative> = Inf.
>
> -- r^(i*y) = exp(i*y*ln(r)).  If y != 0 (i.e. complex power), then taking
> the ln of r = 0 is -Inf.  But what's exp(i*-Inf)?  It's probably NaN,
> since nothing else makes sense.
>
> Note that if y == 0 (real power), then this term is still NaN (y*ln(r)
> = 0*ln(0) = Nan).  However, by convention, 0^<real> is something other
> than NaN.
>
> -- exp(i*theta*x) is just a complex number.
>
> -- exp(-theta*y) is just a real number.
>
> Therefore, for 0^<complex> we have Inf * NaN * <complex> * <real>,
> which is NaN.


[SNIP]

I don't buy this argument. For instance it is claimed  0*ln(0) is NaN, but
IIRC the correct way to determine a value for cases like this is to take the
limit. (It's been a long time since I worked through any limits, so take
this with a grain of salt.)  We are interested in:

Limit (x^(a+bj))  as x -> 0.

First consider the case where a > 0.

x^(a+bj) = x^a x^bj

We guess that the limit is zero as x->0.  In other words |x^a x^bj - 0| -> 0
(I'll leave out all the epsilon and delta stuff and just wave my hands here.
Thus we are interested in the behaviour of |x^a| |x^bj| = |x^a| which does
indeed go to zero as x -> 0 for a > 0.

For a < 0, the limit diverges and it could be argued that the appropriate
value is either Infinity or undefined (NaN), depending on your preference.
Personally, I'd opt for Inf. For a = 0, the limit neither converges nor
diverges; it simply oscillates around the unit circle, so in that case the
value is indeed undefined. Unless b = 0 as well, in which case, the value is
1.

In summary:

   0^(a+bj) = 0 for a > 0
   0^(a+bj) = 1 for a == b == 0
   0^(a+bj) = NaN for a == 0, b != 0
   0^(a+bj) = Inf for a < 0

I suppose one could use NaN + j NaN, as some have proposed, but is seems
unnecessary; an undefined number is undefined and it's not usually necessary
to state that both the real and imaginary parts are undefined. Similarly I
suspect that someone will complain that the result is not necessarily real
for the (a < 0) case.  However, as I recall, in complex analysis, infinity
usually refers to a point infinitely distant from the origin with no
particular phase associated with it.

There is a certain conflict between this mathematical view and an numerical
point of view where Inf and NaN are typically viewed as pure real numbers,
so the above results could be rejected for that reason. However, I suspect
that the phase issues don't matter -- the only time you get back into the
normal numbers when dealing with Inf, is by dividing Inf by a number and in
that case one always get's zero.  Similarly, once you have a NaN, you never
get back to the normal real / imaginary numbers when one of the inputs is a
NaN, so NaN verus NaN + j NaN, are essentially equivalent and the simple NaN
is cleaner.


-- 
.  __
.   |-\
.
.  tim.hochberg at ieee.org
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20080228/aec732a3/attachment.html>