On Wed, Feb 27, 2008 at 5:10 PM, Stuart Brorson <[EMAIL PROTECTED]> 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 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion