Merge pull request #8231 from charris/integer-to-neg-integer-power

ENH: Refactor numpy ** operators for numpy scalar integer powers

2 files changed, 51 insertions, 69 deletions

diff --git a/doc/release/1.12.0-notes.rst b/doc/release/1.12.0-notes.rst
index 473511c3d..e752c1ef1 100644
--- a/doc/release/1.12.0-notes.rst
+++ b/doc/release/1.12.0-notes.rst
@@ -76,15 +76,29 @@ DeprecationWarning to error
 
 ``power`` and ``**`` raise errors for integer to negative integer powers
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-The previous behavior was
+The previous behavior depended on whether numpy scalar integers or numpy
+integer arrays were involved.
+
+For arrays
+
+* Zero to negative integer powers returned least integral value.
+* Both 1, -1 to negative integer powers returned correct values.
+* The remaining integers returned zero when raised to negative integer powers.
+
+For scalars
+
+* Zero to negative integer powers returned least integral value.
+* Both 1, -1 to negative integer powers returned correct values.
+* The remaining integers sometimes returned zero, sometimes the
+  correct float depending on the integer type combination.
+
+All of these cases now raise a ``ValueError`` except for those integer
+combinations whose common type is float, for instance uint64 and int8. It was
+felt that a simple rule was the best way to go rather than have special
+exceptions for the integer units. If you need negative powers, use an inexact
+type.
 
-* zero to negative integer powers returned least integral value.
-* 1, -1 to negative integer powers returned correct values
-* all remaining integers returned zero when raised to negative integer powers.
 
-All of these cases now raise a ``ValueError``. It was felt that a simple rule
-was the best way to go rather than have special exceptions for the units. If
-you need negative powers, use an inexact type.
 
 Relaxed stride checking is the default
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/numpy/core/src/umath/scalarmath.c.src b/numpy/core/src/umath/scalarmath.c.src
index 32a77b6e9..ed6553f69 100644
--- a/numpy/core/src/umath/scalarmath.c.src
+++ b/numpy/core/src/umath/scalarmath.c.src
@@ -220,27 +220,27 @@ static void
  */
 static void
 @name@_ctype_power(@type@ a, @type@ b, @type@ *out) {
-    @type@ temp, ix, mult;
-    /* code from Python's intobject.c, with overflow checking removed. */
-    temp = a;
-    ix = 1;
+    @type@ tmp;
+
+    if (b == 0) {
+        *out = 1;
+        return;
+    }
+    if (a == 1) {
+        *out = 1;
+        return;
+    }
+
+    tmp = b & 1 ? a : 1;
+    b >>= 1;
     while (b > 0) {
+        a *= a;
         if (b & 1) {
-            @name@_ctype_multiply(ix, temp, &mult);
-            ix = mult;
-            if (temp == 0) {
-                break;
-            }
-        }
-        b >>= 1;        /* Shift exponent down by 1 bit */
-        if (b==0) {
-            break;
+            tmp *= a;
         }
-        /* Square the value of temp */
-        @name@_ctype_multiply(temp, temp, &mult);
-        temp = mult;
+        b >>= 1;
     }
-    *out = ix;
+    *out = tmp;
 }
 /**end repeat**/
 
@@ -977,6 +977,7 @@ static PyObject *
             }
             return PyGenericArrType_Type.tp_as_number->nb_power(a,b,NULL);
         case -3:
+        default:
             /*
              * special case for longdouble and clongdouble
              * because they have a recursive getitem in their dtype
@@ -1032,11 +1033,7 @@ static PyObject *
 @name@_power(PyObject *a, PyObject *b, PyObject *NPY_UNUSED(c))
 {
     PyObject *ret;
-    @type@ arg1, arg2;
-    int retstatus;
-    int first;
-    @type@ out = @zero@;
-    @otype@ out1 = @zero@;
+    @type@ arg1, arg2, out;
 
     switch(_@name@_convert2_to_ctypes(a, &arg1, b, &arg2)) {
         case 0:
@@ -1051,6 +1048,7 @@ static PyObject *
             }
             return PyGenericArrType_Type.tp_as_number->nb_power(a,b,NULL);
         case -3:
+        default:
             /*
              * special case for longdouble and clongdouble
              * because they have a recursive getitem in their dtype
@@ -1058,55 +1056,24 @@ static PyObject *
             Py_INCREF(Py_NotImplemented);
             return Py_NotImplemented;
     }
-
     PyUFunc_clearfperr();
 
     /*
      * here we do the actual calculation with arg1 and arg2
      * as a function call.
      */
-    if (@iszero@(arg2)) {
-        out1 = out = @one@;
-    }
-    else if (arg2 < 0) {
-        @oname@_ctype_power(arg1, arg2, &out1);
-    }
-    else {
-        @name@_ctype_power(arg1, arg2, &out);
-    }
-
-    /* Check status flag.  If it is set, then look up what to do */
-    retstatus = PyUFunc_getfperr();
-    if (retstatus) {
-        int bufsize, errmask;
-        PyObject *errobj;
-
-        if (PyUFunc_GetPyValues("@name@_scalars", &bufsize, &errmask,
-                                &errobj) < 0) {
-            return NULL;
-        }
-        first = 1;
-        if (PyUFunc_handlefperr(errmask, errobj, retstatus, &first)) {
-            Py_XDECREF(errobj);
-            return NULL;
-        }
-        Py_XDECREF(errobj);
-    }
-
     if (arg2 < 0) {
-        ret = PyArrayScalar_New(@OName@);
-        if (ret == NULL) {
-            return NULL;
-        }
-        PyArrayScalar_ASSIGN(ret, @OName@, out1);
+        PyErr_SetString(PyExc_ValueError,
+                "Integers to negative integer powers are not allowed.");
+        return NULL;
     }
-    else {
-        ret = PyArrayScalar_New(@Name@);
-        if (ret == NULL) {
-            return NULL;
-        }
-        PyArrayScalar_ASSIGN(ret, @Name@, out);
+    @name@_ctype_power(arg1, arg2, &out);
+
+    ret = PyArrayScalar_New(@Name@);
+    if (ret == NULL) {
+        return NULL;
     }
+    PyArrayScalar_ASSIGN(ret, @Name@, out);
 
     return ret;
 }
@@ -1135,6 +1102,7 @@ static PyObject *
             }
             return PyGenericArrType_Type.tp_as_number->nb_power(a,b,NULL);
         case -3:
+        default:
             /*
              * special case for longdouble and clongdouble
              * because they have a recursive getitem in their dtype