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

[Robert Kern]

On Wed, Feb 27, 2008 at 5: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:
>
>  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.

I believe the failure is occurring because of the coercion to complex.
With plain floats:

In [14]: zeros(2) ** array([-1.0, 0.0])
Out[14]: array([ Inf,   1.])

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

Python gives 1.0:

In [12]: 0.0 ** 0.0
Out[12]: 1.0

I'm not sure about the reasons for this, but I'm willing to assume
that they're acceptable.

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

You could probably stop the analysis here. If a=0, then theta is
already undefined. I believe that NaN+NaN*j is the correct answer.

The relevant function is nc_pow() in numpy/core/src/umathmodule.c. The
problem is that a=(0+0j) is special-cased incorrectly:

    if (ar == 0. && ai == 0.) {
        r->real = 0.;
        r->imag = 0.;
        return;
    }

The preceding if clause (br == 0. && bi == 0.) takes care of the
(0+0j)**(0+0j) case. It's worth noting that the general case at the
bottom returns the expected (NaN+NaN*j). However, we can't just remove
this if-clause; it makes (0+0j)**(-1+0j) return (NaN+NaN*j). It also
makes (0+0j)**(1.5+0j) give (NaN+NaN*j), too.

>  -- 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?

In this case, I think Matlab looks about right.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco