Dispatch to reciprocal, ones_like, copy, sqrt, square inside array_power and array_inplace_power when power is a scalar in [-1, 0, 1, 0.5, 1, 2]. Also, added the ufuncs reciprocal and ones_like.

2 files changed, 200 insertions, 18 deletions

diff --git a/numpy/core/src/arrayobject.c b/numpy/core/src/arrayobject.c
index a90fb7cca..43c148548 100644
--- a/numpy/core/src/arrayobject.c
+++ b/numpy/core/src/arrayobject.c
@@ -2311,7 +2311,11 @@ typedef struct {
                 *multiply,
                 *divide,
                 *remainder,
-                *power,
+                *power,

+                *square,

+                *reciprocal,

+                *ones_like,

+                *copy,
 		*sqrt,
                 *negative,
                 *absolute,
@@ -2341,7 +2345,8 @@ typedef struct {
 static NumericOps n_ops = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 
                            NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
                            NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
-                           NULL, NULL, NULL, NULL, NULL};
+                           NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 

+                           NULL};
 
 /* Dictionary can contain any of the numeric operations, by name. 
   Those not present will not be changed
@@ -2367,7 +2372,11 @@ PyArray_SetNumericOps(PyObject *dict)
         SET(multiply);
         SET(divide);
         SET(remainder);
-        SET(power);
+        SET(power);

+        SET(square);

+        SET(reciprocal);

+        SET(ones_like);

+        SET(copy);
 	SET(sqrt);
         SET(negative);
         SET(absolute);
@@ -2412,7 +2421,11 @@ PyArray_GetNumericOps(void)
         GET(multiply);
         GET(divide);
         GET(remainder);
-        GET(power);
+        GET(power);

+        GET(square);

+        GET(reciprocal);

+        GET(ones_like);

+        GET(copy);
 	GET(sqrt);
         GET(negative);
         GET(absolute);
@@ -2527,6 +2540,16 @@ PyArray_GenericInplaceBinaryFunction(PyArrayObject *m1,
         }
         return PyObject_CallFunction(op, "OOO", m1, m2, m1);
 }
+

+static PyObject *
+PyArray_GenericInplaceUnaryFunction(PyArrayObject *m1, PyObject *op) 
+{
+        if (op == NULL) {
+                Py_INCREF(Py_NotImplemented);
+                return Py_NotImplemented; 
+        }
+        return PyObject_CallFunction(op, "OO", m1, m1);
+}

 
 static PyObject *
 array_add(PyArrayObject *m1, PyObject *m2) 
@@ -2557,12 +2580,61 @@ array_remainder(PyArrayObject *m1, PyObject *m2)
 {
         return PyArray_GenericBinaryFunction(m1, m2, n_ops.remainder);
 }
-
-static PyObject *
-array_power(PyArrayObject *m1, PyObject *m2) 
-{
-        return PyArray_GenericBinaryFunction(m1, m2, n_ops.power);
-}
+

+

+static PyObject *array_float(PyArrayObject *v);

+

+

+static int
+array_power_is_scalar(PyObject *o2, longdouble* exp) 
+{

+    PyObject *temp;

+    const int optimize_fpexps = 1;

+    if (PyInt_Check(o2)) {

+        *exp = (longdouble)PyInt_AsLong(o2);

+        return 1;

+    }

+    if (optimize_fpexps && PyFloat_Check(o2)) {

+        *exp = PyFloat_AsDouble(o2);

+        return 1;

+    }

+    if (PyArray_CheckScalar(o2)) {

+        if (PyArray_ISINTEGER(o2) || (optimize_fpexps && PyArray_ISFLOAT(o2))) {

+            temp = array_float(o2);

+            if (temp != NULL) {

+                *exp = PyFloat_AsDouble(o2);

+                Py_DECREF(temp);

+                return 1;

+            }

+        }

+    }

+    return 0;

+}

+

+static PyObject *

+fast_scalar_power_op(PyArrayObject *a1, PyObject *o2) {

+    double exp;

+    if (PyArray_Check(a1) && (PyArray_ISFLOAT(a1) || PyArray_ISCOMPLEX(a1))) {

+        if (array_power_is_scalar(o2, &exp)) {

+            if (exp == -1.0) return n_ops.reciprocal;

+            if (exp ==  0.0) return n_ops.ones_like;

+            if (exp ==  0.5) return n_ops.sqrt;

+            if (exp ==  1.0) return n_ops.copy;

+            if (exp ==  2.0) return n_ops.square;

+        }

+    }

+    return NULL;

+}

+
+static PyObject *
+array_power(PyArrayObject *a1, PyObject *o2) 
+{

+    PyObject *fastop;

+    fastop = fast_scalar_power_op(a1, o2);

+    if (fastop) return PyArray_GenericUnaryFunction(a1, fastop);

+    return PyArray_GenericBinaryFunction(a1, o2, n_ops.power);

+}

+    
 
 static PyObject *
 array_negative(PyArrayObject *m1) 
@@ -2640,12 +2712,15 @@ static PyObject *
 array_inplace_remainder(PyArrayObject *m1, PyObject *m2) 
 {
         return PyArray_GenericInplaceBinaryFunction(m1, m2, n_ops.remainder);
-}
-
+}

+

 static PyObject *
-array_inplace_power(PyArrayObject *m1, PyObject *m2) 
-{
-        return PyArray_GenericInplaceBinaryFunction(m1, m2, n_ops.power);
+array_inplace_power(PyArrayObject *a1, PyObject *o2) 
+{

+    PyObject *fastop;

+    fastop = fast_scalar_power_op(a1, o2);

+    if (fastop) return PyArray_GenericInplaceUnaryFunction(a1, fastop);

+    return PyArray_GenericInplaceBinaryFunction(a1, o2, n_ops.power);

 }
 
 static PyObject *
@@ -2808,8 +2883,8 @@ array_float(PyArrayObject *v)
         pv = v->descr->f->getitem(v->data, v);
         if (pv == NULL) return NULL;
         if (pv->ob_type->tp_as_number == 0) {
-                PyErr_SetString(PyExc_TypeError, "cannot convert to an "\
-				"int; scalar object is not a number");
+                PyErr_SetString(PyExc_TypeError, "cannot convert to a "\
+				"float; scalar object is not a number");
                 Py_DECREF(pv);
                 return NULL;
         }
diff --git a/numpy/core/src/umathmodule.c.src b/numpy/core/src/umathmodule.c.src
index be7d6ee2f..14e7e4d82 100644
--- a/numpy/core/src/umathmodule.c.src
+++ b/numpy/core/src/umathmodule.c.src
@@ -563,7 +563,7 @@ nc_prod@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
 static void
 nc_quot@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
 {
-
+ 
 	@typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
 	@typ@ d = br*br + bi*bi;
 	r->real = (ar*br + ai*bi)/d;
@@ -1261,6 +1261,113 @@ Py_square(PyObject *o)
 {
     return PyNumber_Multiply(o, o);
 }
+

+

+/**begin repeat
+#TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
+#typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
+*/
+static void
+@TYP@_reciprocal(char **args, intp *dimensions, intp *steps, void *data)
+{
+    intp i, is1 = steps[0], os = steps[1], n = dimensions[0];
+    char *i1 = args[0], *op = args[1];
+
+    for (i = 0; i < n; i++, i1 += is1, op += os) {
+        @typ@ x = *((@typ@ *)i1);
+        *((@typ@ *)op) = 1.0 / x;
+    }
+}
+/**end repeat**/
+
+/**begin repeat
+#TYP=CFLOAT,CDOUBLE,CLONGDOUBLE#
+#typ=float, double, longdouble#
+*/
+static void
+@TYP@_reciprocal(char **args, intp *dimensions, intp *steps, void *data)
+{
+    intp i, is1 = steps[0], os = steps[1], n = dimensions[0];
+    char *i1 = args[0], *op = args[1];
+    c@typ@ *x, *y;
+    @typ@ xr, xi, xmag2;
+
+    for (i = 0; i < n; i++, i1 += is1, op += os) {
+        x = (c@typ@ *)i1;
+        y = (c@typ@ *)op;
+        xr = x->real;
+        xi = x->imag;

+        xmag2 = xr*xr + xi*xi;
+        y->real = xr*xr - xi*xi;
+        y->imag = 2*xr*xi;
+    }
+}
+/**end repeat**/

+

+

+

+static PyObject *
+Py_reciprocal(PyObject *o)
+{

+    PyObject *one, *result;

+    one = PyInt_FromLong(1);

+    if (!one) return NULL;
+    result = PyNumber_Divide(one, o);

+    Py_DECREF(one);

+    return result;
+}

+

+

+/**begin repeat
+#TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
+#typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
+*/
+static void
+@TYP@_ones_like(char **args, intp *dimensions, intp *steps, void *data)
+{
+    intp i, is1 = steps[0], os = steps[1], n = dimensions[0];
+    char *i1 = args[0], *op = args[1];
+
+    for (i = 0; i < n; i++, i1 += is1, op += os) {
+        @typ@ x = *((@typ@ *)i1);
+        *((@typ@ *)op) = 1.0;
+    }
+}
+/**end repeat**/
+
+/**begin repeat
+#TYP=CFLOAT,CDOUBLE,CLONGDOUBLE#
+#typ=float, double, longdouble#
+*/
+static void
+@TYP@_ones_like(char **args, intp *dimensions, intp *steps, void *data)
+{
+    intp i, is1 = steps[0], os = steps[1], n = dimensions[0];
+    char *i1 = args[0], *op = args[1];
+    c@typ@ *x, *y;
+    @typ@ xr, xi, xmag2;
+
+    for (i = 0; i < n; i++, i1 += is1, op += os) {
+        x = (c@typ@ *)i1;
+        y = (c@typ@ *)op;
+        xr = x->real;
+        xi = x->imag;

+        xmag2 = xr*xr + xi*xi;
+        y->real = 1.0;
+        y->imag = 0.0;
+    }
+}
+/**end repeat**/

+

+

+static PyObject *
+Py_get_one(PyObject *o)
+{

+    return PyInt_FromLong(1.0);

+}

+

+

+

 
 /**begin repeat