[Numpy-discussion] Handling of numpy.power(0, <something>)
Robert Kern
robert.kern@gmail....
Wed Feb 27 18:45:20 CST 2008
On Wed, Feb 27, 2008 at 5:10 PM, Stuart Brorson <sdb@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
More information about the Numpy-discussion
mailing list