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

[Stuart Brorson]

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:

In [67]: A = numpy.array([0, 0, 0])

In [68]: B = numpy.array([-1, 0, 1+1j])

In [69]: numpy.power(A, B)
Out[69]: array([ 0.+0.j,  1.+0.j,  0.+0.j])

IMO, the answers should be [Inf, NaN, and NaN].  The reasons:

**  0^-1 is 1/0, which is infinity.  Not much argument here, I would
think.

**  0^0:  This is problematic.  People smarter than I have argued for
both NaN and for 1, although I understand that 1 is the preferred
value nowadays.  If the NumPy gurus also think so, then I buy it.

**  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.

Another observation to chew on.  I know NumPy != Matlab, but FWIW,
here's what Matlab says about these values:

>> A = [0, 0, 0]

A =

      0     0     0

>> B = [-1, 0, 1+1*i]

B =

   -1.0000           0      1.0000 + 1.0000i

>> A .^ B

ans =

       Inf      1.0000         NaN +    NaNi



Any reactions to this?  Does NumPy just make library calls when
computing power, or does it do any trapping of corner cases?  And
should the returns from power conform to the above suggestions?

Regards,

Stuart Brorson
Interactive Supercomputing, inc.
135 Beaver Street | Waltham | MA | 02452 | USA
http://www.interactivesupercomputing.com/