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

[Robert Kern]

On Thu, Feb 28, 2008 at 3:37 PM, Christian Heimes <lists at cheimes.de> wrote:
> Stuart Brorson wrote:

>  > Also, what do these specs say about 0^<complex>?
>
>  See for yourself
>  http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf The
>  interesting information are in Annex F.9 and Annex G.6.

Thank you! This is quite helpful. G.6.3 tells us most of the special
cases for cexp() and clog(). In G.6.4, cpow(z,c) is more-or-less
defined as cexp(c * clog(z)).

One thing that I was reminded of is that decomposing 0 into the form
r*cexp(theta*i) is *not* indeterminate in theta in IEEE floating point
arithmetic as I claimed earlier. IEEE floating point models an "affine
extension" of the reals. I'm not particularly clear on the more
abstract mathematical consequences of that, but the practical upshot
is that we have separate +Inf and -Inf as well as +0 and -0. For
complex numbers, we also have (+0+0j), (-0+0j), (+0-0j), and (-0-0j)
as well as a similar permutation set for the Infs. The practical use
of the signed zeros is that clog(+0+0j) can be defined as the limit of
clog(z) coming to the origin from a particular direction and thus have
a somewhat-arbitrary, but well-defined limit. Namely, the standard
defines clog(+0+0j) as (-Inf+pi*j).

I'll get around to working through all of the cases at some point, and
make a table. Suffice it to say, for the original problem in the
thread:

  power(+0.0, (1+1j)) == cexp((1+1j) * clog(+0.0+0.0j))
        == cexp((1+1j) * (-inf+pi*1j))
        == cexp((-inf-inf*1j))
        == (+-0 + (+-0j))

where the signs are unspecified by the standard.

So numpy's behavior happens to be correct for power(0.0, (1+1j)). It
is incorrect for power(0.0, (-1.0+0j)).

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