[Matrix-SIG] Exponentiation (was: [Matrix-SIG] Warning in re 1.5.2b1)

[Charles G Waldman]

Paul Dubois writes:

 > I want to warn people that we received a complaint that a real array
 > raised to a real power gets a TypeError. We have discovered, somewhat
 > suprisingly, that this happens only with the 1.5.2b1 release of Python
 > (NumPy, release 8).  With 1.5.1 and Release 8, it works.

The culprit is some code in Objects/abstract.c which is intended to
prevent us from attempting to raise a negative number to a float
power.

There is a check in the do_pow function which attempts to ensure that
when we compute x**y and y is float, x is nonnegative.  When x is an
array, this check cannot be performed, and an error is raised.

There are some interesting differences in behavior in Python between
the 1.5, 1.5.1 and 1.5.2 version with respect to exponentiation, as
these examples show:


Version 1.5 refuses to compute x**y if x<0 and y is real, even if y is 
actually integral:

  Python 1.5
  >>> (-2)**2.0
  ValueError: negative number to float power
  >>> (-2)**2.1
  ValueError: negative number to float power
  >>> (-2.0)**2.0
  ValueError: negative number to float power
  >>> (-2.0)**2.1
  ValueError: negative number to float power



Version 1.5.1 lets us compute -2**2.0 but not (-2.0)**2.0:

  Python 1.5.1 
  >>> (-2)**2.0
  4.0
  >>> (-2.0)**2.0
  ValueError: negative number to float power
  >>> (-2)**2.1
  ValueError: negative number to a float power


Version 1.5.2b2 goes back to the same behavior as 1.5:

  Python 1.5.2b2 
  >>> (-2)**2.0
  ValueError: negative number to float power
  >>> (-2)**2.1
  ValueError: negative number to float power


Anyhow, this may not seem related to the problem of raising arrays to
real powers, but it is.  The difference between 1.5 and 1.5.1 is due
to the following code change:

diff -u Python-1.5/Objects/abstract.c Python-1.5.1/Objects/
--- Python-1.5/Objects/abstract.c       Tue May 20 18:09:08 1997
+++ Python-1.5.1/Objects/abstract.c     Fri Apr  3 18:38:59 1998
@@ -463,7 +463,8 @@
                                "pow() requires numeric arguments");
                return NULL;
        }
-       if (PyFloat_Check(w) && PyFloat_AsDouble(v) < 0.0) {
+       if (PyFloat_Check(v) && PyFloat_Check(w) &&
+           PyFloat_AsDouble(v) < 0.0) {
                if (!PyErr_Occurred())
                    PyErr_SetString(PyExc_ValueError,
                                    "negative number to float power");

Here "v" is the base and "w" is the exponent.  The PyFloat_Check(v) in
the 1.5.1 code short-circuits the v<0 test if v is not a float; that's
why in 1.5.1 we can compute (-2)**2.0 but not (-2.0)**2.0, and also
why we can raise arrays to real powers.  The fact that we can do this
in 1.5.1 is really a case of "slipping between the cracks" in this
test.


In 1.5.2 this code has been changed to:

        if (
#ifndef WITHOUT_COMPLEX
            !PyComplex_Check(v) && 
#endif
            PyFloat_Check(w) && PyFloat_AsDouble(v) < 0.0) {
                if (!PyErr_Occurred())
                    PyErr_SetString(PyExc_ValueError,
                                    "negative number to float power");
                return NULL

so we can no longer get away with raising arrays to real powers.


An easy fix for this is simply to remove this chunk of code from
abstract.c entirely.  There is already code in
floatobject.c:float_pow() to raise an exception in the bad cases.

With the checks removed from abstract.c:do_pow (the "if" block
beginning at line 589) I get what I consider to be more reasonable
behavior:

Python 1.5.2b2 (patched)

>>> (-2)**2.0
4.0
>>> (-2.0)**2.0
4.0
>>> (-2.0)**2.1
ValueError: negative number to a float power
>>> from Numeric import *
>>> a = array((1,2,3))
>>> a**2.0
array([ 1.,  4.,  9.])
>>> a = array((1,-2,3))
>>> a**2.0
array([ 1.,  4.,  9.])
>>> a**2.1
ValueError: math domain error



Here's a patch for abstract.c
===================================================================
RCS file: /projects/cvsroot/python/dist/src/Objects/abstract.c,v
retrieving revision 2.26
diff -u -u -r2.26 abstract.c
--- abstract.c  1998/12/10 16:54:48     2.26
+++ abstract.c  1999/01/10 07:51:16
@@ -586,16 +586,6 @@
                                "pow(x, y) requires numeric arguments");
                return NULL;
        }
-       if (
-#ifndef WITHOUT_COMPLEX
-            !PyComplex_Check(v) && 
-#endif
-            PyFloat_Check(w) && PyFloat_AsDouble(v) < 0.0) {
-               if (!PyErr_Occurred())
-                   PyErr_SetString(PyExc_ValueError,
-                                   "negative number to float power");
-               return NULL;
-       }
        if (PyNumber_Coerce(&v, &w) != 0)
                return NULL;
        if ((f = v->ob_type->tp_as_number->nb_power) != NULL)