Cleanup ufunc loops.

At this point loops are separated into variable kinds, so there is a fair amount
of duplication. I will probably merge loops that look the same in a later
commit. There are no changes to current behavior of loops, this will also be
changed in later work to deal with nans and such.

1 files changed, 790 insertions, 890 deletions

diff --git a/numpy/core/src/umathmodule.c.src b/numpy/core/src/umathmodule.c.src
index 77675c441..6f41cd4fa 100644
--- a/numpy/core/src/umathmodule.c.src
+++ b/numpy/core/src/umathmodule.c.src
@@ -1,6 +1,10 @@
 /* -*- c -*- */
 
 /*
+ * vim:syntax=c
+ */
+
+/*
  *****************************************************************************
  **                            INCLUDES                                     **
  *****************************************************************************
@@ -13,14 +17,40 @@
 #include "config.h"
 #include <math.h>
 
+#ifndef M_PI
+#define M_PI 3.14159265358979323846264338328
+#endif
+
 /*
  *****************************************************************************
  **                     BASIC MATH FUNCTIONS                                **
  *****************************************************************************
  */
 
-/* A whole slew of basic math functions are provided originally
-   by Konrad Hinsen. */
+float degreesf(float x) {
+    return x * (float)(180.0/M_PI);
+}
+double degrees(double x) {
+    return x * (180.0/M_PI);
+}
+longdouble degreesl(longdouble x) {
+    return x * (180.0L/M_PI);
+}
+
+float radiansf(float x) {
+    return x * (float)(M_PI/180.0);
+}
+double radians(double x) {
+    return x * (M_PI/180.0);
+}
+longdouble radiansl(longdouble x) {
+    return x * (M_PI/180.0L);
+}
+
+/* 
+ * A whole slew of basic math functions are provided originally
+ * by Konrad Hinsen.
+ */
 
 #if !defined(__STDC__) && !defined(_MSC_VER)
 extern double fmod (double, double);
@@ -28,9 +58,6 @@ extern double frexp (double, int *);
 extern double ldexp (double, int);
 extern double modf (double, double *);
 #endif
-#ifndef M_PI
-#define M_PI 3.14159265358979323846264338328
-#endif
 
 
 #if defined(DISTUTILS_USE_SDK)
@@ -443,26 +470,6 @@ trunc(double x)
 #define isfinitef(x) (!(isinff((x)) || isnanf((x))))
 #define isfinitel(x) (!(isinfl((x)) || isnanl((x))))
 
-float degreesf(float x) {
-    return x * (float)(180.0/M_PI);
-}
-double degrees(double x) {
-    return x * (180.0/M_PI);
-}
-longdouble degreesl(longdouble x) {
-    return x * (180.0L/M_PI);
-}
-
-float radiansf(float x) {
-    return x * (float)(M_PI/180.0);
-}
-double radians(double x) {
-    return x * (M_PI/180.0);
-}
-longdouble radiansl(longdouble x) {
-    return x * (M_PI/180.0L);
-}
-
 /* First, the C functions that do the real work */
 
 /* if C99 extensions not available then define dummy functions that use the
@@ -624,6 +631,61 @@ static float expm1f(float x)
 }
 #endif
 
+/*
+ *****************************************************************************
+ **                        PYTHON OBJECT FUNCTIONS                          **
+ *****************************************************************************
+ */
+
+static PyObject *
+Py_square(PyObject *o)
+{
+    return PyNumber_Multiply(o, o);
+}
+
+static PyObject *
+Py_get_one(PyObject *o)
+{
+    return PyInt_FromLong(1);
+}
+
+static PyObject *
+Py_reciprocal(PyObject *o)
+{
+    PyObject *one = PyInt_FromLong(1);
+    PyObject *result;
+
+    if (!one) {
+        return NULL;
+    }
+    result = PyNumber_Divide(one, o);
+    Py_DECREF(one);
+    return result;
+}
+
+/**begin repeat
+ * #Kind = Max, Min#
+ * #OP = >=, <=#
+ */
+static PyObject *
+_npy_Object@Kind@(PyObject *i1, PyObject *i2)
+{
+    PyObject *result;
+    int cmp;
+
+    if (PyObject_Cmp(i1, i2, &cmp) < 0) {
+        return NULL;
+    }
+    if (cmp @OP@ 0) {
+        result = i1;
+    }
+    else {
+        result = i2;
+    }
+    Py_INCREF(result);
+    return result;
+}
+/**end repeat**/
 
 /*
  *****************************************************************************
@@ -632,11 +694,12 @@ static float expm1f(float x)
  */
 
 
-/* Don't pass structures between functions (only pointers) because how
-   structures are passed is compiler dependent and could cause
-   segfaults if ufuncobject.c is compiled with a different compiler
-   than an extension that makes use of the UFUNC API
-*/
+/* 
+ * Don't pass structures between functions (only pointers) because how
+ * structures are passed is compiler dependent and could cause
+ * segfaults if ufuncobject.c is compiled with a different compiler
+ * than an extension that makes use of the UFUNC API
+ */
 
 /**begin repeat
 
@@ -650,9 +713,10 @@ static c@typ@ nc_half@c@ = {0.5, 0.};
 static c@typ@ nc_i@c@ = {0., 1.};
 static c@typ@ nc_i2@c@ = {0., 0.5};
 /*
-  static c@typ@ nc_mi@c@ = {0., -1.};
-  static c@typ@ nc_pi2@c@ = {M_PI/2., 0.};
-*/
+ *   static c@typ@ nc_mi@c@ = {0., -1.};
+ *   static c@typ@ nc_pi2@c@ = {M_PI/2., 0.};
+ */
+
 
 static void
 nc_sum@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
@@ -681,7 +745,7 @@ nc_neg@c@(c@typ@ *a, c@typ@ *r)
 static void
 nc_prod@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
 {
-    register @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+    @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
     r->real = ar*br - ai*bi;
     r->imag = ar*bi + ai*br;
     return;
@@ -691,8 +755,8 @@ static void
 nc_quot@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
 {
 
-    register @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
-    register @typ@ d = br*br + bi*bi;
+    @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;
     r->imag = (ai*br - ar*bi)/d;
     return;
@@ -801,8 +865,10 @@ nc_pow@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
             return;
         }
     }
-    /* complexobect.c uses an inline version of this formula
-       investigate whether this had better performance or accuracy */
+    /*
+     * complexobect.c uses an inline version of this formula
+     * investigate whether this had better performance or accuracy
+     */
     nc_log@c@(a, r);
     nc_prod@c@(r, b, r);
     nc_exp@c@(r, r);
@@ -823,6 +889,10 @@ nc_prodi@c@(c@typ@ *x, c@typ@ *r)
 static void
 nc_acos@c@(c@typ@ *x, c@typ@ *r)
 {
+    /*
+     * return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i,
+     * nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))));
+     */
     nc_prod@c@(x,x,r);
     nc_diff@c@(&nc_1@c@, r, r);
     nc_sqrt@c@(r, r);
@@ -832,14 +902,15 @@ nc_acos@c@(c@typ@ *x, c@typ@ *r)
     nc_prodi@c@(r, r);
     nc_neg@c@(r, r);
     return;
-    /* return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i,
-       nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))));
-    */
 }
 
 static void
 nc_acosh@c@(c@typ@ *x, c@typ@ *r)
 {
+    /*
+     * return nc_log(nc_sum(x,
+     * nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1)))));
+     */
     c@typ@ t;
 
     nc_sum@c@(x, &nc_1@c@, &t);
@@ -850,15 +921,15 @@ nc_acosh@c@(c@typ@ *x, c@typ@ *r)
     nc_sum@c@(x, r, r);
     nc_log@c@(r, r);
     return;
-    /*
-      return nc_log(nc_sum(x,
-      nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1)))));
-    */
 }
 
 static void
 nc_asin@c@(c@typ@ *x, c@typ@ *r)
 {
+    /*
+     * return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x),
+     * nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))));
+     */
     c@typ@ a, *pa=&a;
     nc_prod@c@(x, x, r);
     nc_diff@c@(&nc_1@c@, r, r);
@@ -869,30 +940,29 @@ nc_asin@c@(c@typ@ *x, c@typ@ *r)
     nc_prodi@c@(r, r);
     nc_neg@c@(r, r);
     return;
-    /*
-      return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x),
-      nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))));
-    */
 }
 
 
 static void
 nc_asinh@c@(c@typ@ *x, c@typ@ *r)
 {
+    /*
+     * return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x));
+     */
     nc_prod@c@(x, x, r);
     nc_sum@c@(&nc_1@c@, r, r);
     nc_sqrt@c@(r, r);
     nc_sum@c@(r, x, r);
     nc_log@c@(r, r);
     return;
-    /*
-      return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x));
-    */
 }
 
 static void
 nc_atan@c@(c@typ@ *x, c@typ@ *r)
 {
+    /*
+     * return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
+     */
     c@typ@ a, *pa=&a;
     nc_diff@c@(&nc_i@c@, x, pa);
     nc_sum@c@(&nc_i@c@, x, r);
@@ -900,14 +970,14 @@ nc_atan@c@(c@typ@ *x, c@typ@ *r)
     nc_log@c@(r,r);
     nc_prod@c@(&nc_i2@c@, r, r);
     return;
-    /*
-      return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
-    */
 }
 
 static void
 nc_atanh@c@(c@typ@ *x, c@typ@ *r)
 {
+    /*
+     * return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
+     */
     c@typ@ a, *pa=&a;
     nc_diff@c@(&nc_1@c@, x, r);
     nc_sum@c@(&nc_1@c@, x, pa);
@@ -915,9 +985,6 @@ nc_atanh@c@(c@typ@ *x, c@typ@ *r)
     nc_log@c@(r, r);
     nc_prod@c@(&nc_half@c@, r, r);
     return;
-    /*
-      return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
-    */
 }
 
 static void
@@ -1018,1137 +1085,973 @@ nc_tanh@c@(c@typ@ *x, c@typ@ *r)
  *****************************************************************************
  */
 
+#define OUTPUT_LOOP\
+    char *op = args[1];\
+    intp os = steps[1];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, op += os)
+
+#define UNARY_LOOP\
+    char *ip1 = args[0], *op = args[1];\
+    intp is1 = steps[0], os = steps[1];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, ip1 += is1, op += os)
+
+#define UNARY_LOOP_TWO_OUT\
+    char *ip1 = args[0], *op1 = args[1], *op2 = args[2];\
+    intp is1 = steps[0], os1 = steps[1], os2 = steps[2];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, ip1 += is1, op1 += os1, op2 += os2)
+
+#define BINARY_LOOP\
+    char *ip1 = args[0], *ip2 = args[1], *op = args[2];\
+    intp is1 = steps[0], is2 = steps[1], os = steps[2];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op += os)
+
+#define BINARY_LOOP_TWO_OUT\
+    char *ip1 = args[0], *ip2 = args[1], *op1 = args[2], *op2 = args[3];\
+    intp is1 = steps[0], is2 = steps[1], os1 = steps[2], os2 = steps[3];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1, op2 += os2)
 
-/**begin repeat
+/*
+ *****************************************************************************
+ **                             BOOLEAN LOOPS                               **
+ *****************************************************************************
+ */
 
-   #TYPE=(BOOL, BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*2#
-   #OP=||, +*13, ^, -*13#
-   #kind=add*14, subtract*14#
-   #typ=(Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*2#
-*/
+#define BOOL_invert BOOL_logical_not
+#define BOOL_negative BOOL_logical_not
+#define BOOL_add BOOL_logical_or
+#define BOOL_bitwise_and BOOL_logical_and
+#define BOOL_bitwise_or BOOL_logical_or
+#define BOOL_bitwise_xor BOOL_logical_xor
+#define BOOL_multiply BOOL_logical_and
+#define BOOL_subtract BOOL_logical_xor
+
+/**begin repeat
+ * #kind = equal, not_equal, greater, greater_equal, less, less_equal,
+ *         logical_and, logical_or#
+ * #OP =  ==, !=, >, >=, <, <=, &&, ||#
+ **/
 
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((@typ@ *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+    BINARY_LOOP {
+        Bool in1 = *((Bool *)ip1) != 0;
+        Bool in2 = *((Bool *)ip2) != 0;
+        *((Bool *)op)= in1 @OP@ in2;
     }
 }
-
 /**end repeat**/
 
-/**begin repeat
-
-   #TYPE=(CFLOAT, CDOUBLE, CLONGDOUBLE)*2#
-   #OP=+*3,-*3#
-   #kind=add*3,subtract*3#
-   #typ=(float, double, longdouble)*2#
-*/
-
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+BOOL_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        ((@typ@ *)op)[0]=((@typ@ *)i1)[0] @OP@ ((@typ@ *)i2)[0];
-        ((@typ@ *)op)[1]=((@typ@ *)i1)[1] @OP@ ((@typ@ *)i2)[1];
+    BINARY_LOOP {
+        Bool in1 = *((Bool *)ip1) != 0;
+        Bool in2 = *((Bool *)ip2) != 0;
+        *((Bool *)op)= (in1 && !in2) || (!in1 && in2);
     }
 }
 
-/**end repeat**/
-
-
+/**begin repeat
+ * #kind = maximum, minimum#
+ * #OP =  >, <#
+ **/
 static void
-BOOL_multiply(char **args, intp *dimensions, intp *steps, void *func) {
-    register intp i;
-    intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((Bool *)op) = *((Bool *)i1) && *((Bool *)i2);
+BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        Bool in1 = *((Bool *)ip1) != 0;
+        Bool in2 = *((Bool *)ip2) != 0;
+        *((Bool *)op) = (in1 @OP@ in2) ? in1 : in2;
     }
 }
+/**end repeat**/
 
 /**begin repeat
-
-   #TYP= BYTE, SHORT, INT, LONG, LONGLONG, UBYTE, USHORT, UINT, ULONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE#
-   #typ= byte, short, int, long, longlong, ubyte, ushort, uint, ulong, ulonglong, float, double, longdouble#
-*/
+ * #kind = absolute, logical_not#
+ * #OP =  !=, ==#
+ **/
 static void
-@TYP@_multiply(char **args, intp *dimensions, intp *steps, void *func)
+BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((@typ@ *)op) = (*((@typ@ *)i1)) * (*((@typ@ *)i2));
+    UNARY_LOOP {
+        Bool in1 = *(Bool *)ip1;
+        *((Bool *)op) = in1 @OP@ 0;
     }
 }
 /**end repeat**/
 
-
-/**begin repeat
-   #TYP= CFLOAT, CDOUBLE, CLONGDOUBLE#
-   #typ= float, double, longdouble#
-   #c=f,,l#
-*/
 static void
-@TYP@_multiply(char **args, intp *dimensions, intp *steps, void *func)
+BOOL_ones_like(char **args, intp *dimensions, intp *steps, void *data)
 {
-    register intp i;
-    intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        register @typ@ ar=((c@typ@ *)i1)->real, \
-            ai=((c@typ@ *)i1)->imag,        \
-            br=((c@typ@ *)i2)->real,        \
-            bi=((c@typ@ *)i2)->imag;
-        ((c@typ@ *)op)->real = ar*br - ai*bi;
-        ((c@typ@ *)op)->imag = ar*bi + ai*br;
+    OUTPUT_LOOP {
+        *((Bool *)op) = 1;
     }
 }
 
+
+/*
+ *****************************************************************************
+ **                           INTEGER LOOPS
+ *****************************************************************************
+ */
+
+/**begin repeat
+ * #type = byte, short, int, long, longlong#
+ * #TYPE = BYTE, SHORT, INT, LONG, LONGLONG#
+ * #ftype = float, float, double, double, double#
+ */
+
+/**begin repeat1
+ * both signed and unsigned integer types
+ * #s = , u#
+ * #S = , U#
+ */
+
+#define @S@@TYPE@_floor_divide @S@@TYPE@_divide
+
 static void
-@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func)
+@S@@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *data)
 {
-    register intp i;
-    intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        register @typ@ ar=((c@typ@ *)i1)->real, \
-            ai=((c@typ@ *)i1)->imag,        \
-            br=((c@typ@ *)i2)->real,        \
-            bi=((c@typ@ *)i2)->imag;
-        register @typ@ d = br*br + bi*bi;
-        ((c@typ@ *)op)->real = (ar*br + ai*bi)/d;
-        ((c@typ@ *)op)->imag = (ai*br - ar*bi)/d;
+    OUTPUT_LOOP {
+        *((@s@@type@ *)op) = 1;
     }
 }
 
 static void
-@TYP@_floor_divide(char **args, intp *dimensions, intp *steps, void *func)
+@S@@TYPE@_square(char **args, intp *dimensions, intp *steps, void *data)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        register @typ@ ar=((c@typ@ *)i1)->real,    \
-            ai=((c@typ@ *)i1)->imag,           \
-            br=((c@typ@ *)i2)->real,           \
-            bi=((c@typ@ *)i2)->imag;
-        register @typ@ d = br*br + bi*bi;
-        ((c@typ@ *)op)->real = floor@c@((ar*br + ai*bi)/d);
-        ((c@typ@ *)op)->imag = 0;
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op) = in1*in1;
     }
 }
 
-#define @TYP@_true_divide @TYP@_divide
-/**end repeat**/
-
-
-/**begin repeat
-   #TYP=UBYTE,USHORT,UINT,ULONG,ULONGLONG#
-   #typ=ubyte, ushort, uint, ulong, ulonglong#
-*/
 static void
-@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func)
+@S@@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data)
 {
-    register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        if (*((@typ@ *)i2)==0) {
-            generate_divbyzero_error();
-            *((@typ@ *)op)=0;
-        }
-        else {
-            *((@typ@ *)op)= *((@typ@ *)i1) / *((@typ@ *)i2);
-        }
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op) = 1.0/in1;
     }
 }
-/**end repeat**/
 
-
-/**begin repeat
-   #TYP=BYTE,SHORT,INT,LONG,LONGLONG#
-   #typ=char, short, int, long, longlong#
-*/
 static void
-@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func)
+@S@@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    @typ@ x, y, tmp;
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        y = *((@typ@ *)i2);
-        if (y == 0) {
-            generate_divbyzero_error();
-            *((@typ@ *)op)=0;
-        }
-        else {
-            x = *((@typ@ *)i1);
-            tmp = x / y;
-            if (((x > 0) != (y > 0)) && (x % y != 0)) tmp--;
-            *((@typ@ *)op)= tmp;
-        }
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op) = in1;
     }
 }
-/**end repeat**/
 
-
-/**begin repeat
-   #TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG#
-   #typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong#
-   #otyp=float*4, double*6#
-*/
 static void
-@TYP@_true_divide(char **args, intp *dimensions, intp *steps, void *func)
+@S@@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        if (*((@typ@ *)i2)==0) {
-            generate_divbyzero_error();
-            *((@otyp@ *)op)=0;
-        }
-        else {
-            *((@otyp@ *)op)=                                \
-                (@otyp@)((double)*((@typ@ *)i1) / (double)*((@typ@ *)i2));
-        }
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op) = (@s@@type@)(-(@type@)in1);
     }
 }
-#define @TYP@_floor_divide @TYP@_divide
-/**end repeat**/
-
 
-
-/**begin repeat
-   #TYP=FLOAT,DOUBLE,LONGDOUBLE#
-   #typ=float,double,longdouble#
-*/
 static void
-@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func)
+@S@@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((@typ@ *)op)=*((@typ@ *)i1) / *((@typ@ *)i2);
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((Bool *)op) = !in1;
     }
 }
-#define @TYP@_true_divide @TYP@_divide
-/**end repeat**/
 
-/**begin repeat
-
-   #TYP=FLOAT,DOUBLE,LONGDOUBLE#
-   #typ=float,double,longdouble#
-   #c=f,,l#
-*/
 static void
-@TYP@_floor_divide(char **args, intp *dimensions, intp *steps, void *func)
+@S@@TYPE@_invert(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((@typ@ *)op)=floor@c@(*((@typ@ *)i1) / *((@typ@ *)i2));
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op) = ~in1;
     }
 }
-/**end repeat**/
 
-/**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#
-*/
+/**begin repeat2
+ * Arithmetic
+ * #kind = add, subtract, multiply, bitwise_and, bitwise_or, bitwise_xor,
+ *          left_shift, right_shift#
+ * #OP = +, -,*, &, |, ^, <<, >>#
+ */
 static void
-@TYP@_square(char **args, intp *dimensions, intp *steps, void *data)
+@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    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) = x*x;
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((@s@@type@ *)op) = in1 @OP@ in2;
     }
 }
-/**end repeat**/
+/**end repeat2**/
 
-/**begin repeat
-   #TYP=CFLOAT,CDOUBLE,CLONGDOUBLE#
-   #typ=float, double, longdouble#
-*/
+/**begin repeat2
+ * #kind = equal, not_equal, greater, greater_equal, less, less_equal,
+ *         logical_and, logical_or#
+ * #OP =  ==, !=, >, >=, <, <=, &&, ||#
+ */
 static void
-@TYP@_square(char **args, intp *dimensions, intp *steps, void *data)
+@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    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;
-
-    for (i = 0; i < n; i++, i1 += is1, op += os) {
-        x = (c@typ@ *)i1;
-        y = (c@typ@ *)op;
-        xr = x->real;
-        xi = x->imag;
-        y->real = xr*xr - xi*xi;
-        y->imag = 2*xr*xi;
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((Bool *)op) = in1 @OP@ in2;
     }
 }
-/**end repeat**/
+/**end repeat2**/
 
-static PyObject *
-Py_square(PyObject *o)
+static void
+@S@@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
 {
-    return PyNumber_Multiply(o, o);
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((Bool *)op)= (in1 && !in2) || (!in1 && in2);
+    }
 }
 
-
-/**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#
-*/
+/**begin repeat2
+ * #kind = maximum, minimum#
+ * #OP =  >, <#
+ **/
 static void
-@TYP@_reciprocal(char **args, intp *dimensions, intp *steps, void *data)
+@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    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) = (@typ@) (1.0 / x);
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((@s@@type@ *)op) = (in1 @OP@ in2) ? in1 : in2;
     }
 }
-/**end repeat**/
+/**end repeat2**/
 
-/**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, r, denom;
-
-    for (i = 0; i < n; i++, i1 += is1, op += os) {
-        x = (c@typ@ *)i1;
-        y = (c@typ@ *)op;
-        xr = x->real;
-        xi = x->imag;
-        if (fabs(xi) <= fabs(xr)) {
-            r = xi / xr;
-            denom = xr + xi * r;
-            y->real = 1 / denom;
-            y->imag = -r / denom;
-        } else {
-            r = xr / xi;
-            denom = xr * r + xi;
-            y->real = r / denom;
-            y->imag = -1 / denom;
+@S@@TYPE@_true_divide(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@ftype@ *)op) = 0;
+        }
+        else {
+            *((@ftype@ *)op) = (@ftype@)in1 / (@ftype@)in2;
         }
     }
 }
-/**end repeat**/
-
 
-static PyObject *
-Py_reciprocal(PyObject *o)
+static void
+@S@@TYPE@_power(char **args, intp *dimensions, intp *steps, void *func)
 {
-    PyObject *one, *result;
-    one = PyInt_FromLong(1);
-    if (!one) return NULL;
-    result = PyNumber_Divide(one, o);
-    Py_DECREF(one);
-    return result;
+    BINARY_LOOP {
+        const @ftype@ in1 = (@ftype@)*(@s@@type@ *)ip1;
+        const @ftype@ in2 = (@ftype@)*(@s@@type@ *)ip2;
+        *((@s@@type@ *)op) = (@s@@type@) pow(in1, in2);
+    }
 }
 
-static PyObject *
-_npy_ObjectMax(PyObject *i1, PyObject *i2)
+static void
+@S@@TYPE@_fmod(char **args, intp *dimensions, intp *steps, void *func)
 {
-    int cmp;
-    PyObject *res;
-    if (PyObject_Cmp(i1, i2, &cmp) < 0) return NULL;
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@s@@type@ *)op) = 0;
+        }
+        else {
+            *((@s@@type@ *)op)= in1 % in2;
+        }
 
-    if (cmp >= 0) {
-        res = i1;
     }
-    else {
-        res = i2;
-    }
-    Py_INCREF(res);
-    return res;
 }
 
-static PyObject *
-_npy_ObjectMin(PyObject *i1, PyObject *i2)
-{
-    int cmp;
-    PyObject *res;
-    if (PyObject_Cmp(i1, i2, &cmp) < 0) return NULL;
+/**end repeat1**/
 
-    if (cmp <= 0) {
-        res = i1;
-    }
-    else {
-        res = i2;
+static void
+U@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+{
+    UNARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        *((u@type@ *)op) = in1;
     }
-    Py_INCREF(res);
-    return res;
 }
 
-/* ones_like is defined here because it's used for x**0 */
-
-/**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)
+@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
 {
-    intp i, os = steps[1], n = dimensions[0];
-    char *op = args[1];
-
-    for (i = 0; i < n; i++, op += os) {
-        *((@typ@ *)op) = 1;
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op) = (in1 >= 0) ? in1 : -in1;
     }
 }
-/**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)
+U@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func)
 {
-    intp i, is1 = steps[0], os = steps[1], n = dimensions[0];
-    char *i1 = args[0], *op = args[1];
-    c@typ@ *y;
-
-    for (i = 0; i < n; i++, i1 += is1, op += os) {
-        y = (c@typ@ *)op;
-        y->real = 1.0;
-        y->imag = 0.0;
+    UNARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        *((u@type@ *)op) = in1 > 0 ? 1 : 0;
     }
 }
-/**end repeat**/
 
-static PyObject *
-Py_get_one(PyObject *o)
+static void
+@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func)
 {
-    return PyInt_FromLong(1);
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0);
+    }
 }
 
-
-/**begin repeat
-   #TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG#
-   #typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong#
-   #btyp=float*4, double*6#
-*/
 static void
-@TYP@_power(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, is1=steps[0],is2=steps[1];
-    register intp os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    @btyp@ x, y;
-
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        x = (@btyp@)*((@typ@ *)i1);
-        y = (@btyp@)*((@typ@ *)i2);
-        *((@typ@ *)op) = (@typ@) pow(x,y);
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@type@ *)op) = 0;
+        }
+        else if (((in1 > 0) != (in2 > 0)) && (in1 % in2 != 0)) {
+            *((@type@ *)op) = in1/in2 - 1;
+        }
+        else {
+            *((@type@ *)op) = in1/in2;
+        }
     }
 }
-/**end repeat**/
 
-/**begin repeat
-   #TYP=UBYTE, BYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE#
-   #typ=ubyte, char, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
-*/
 static void
-@TYP@_conjugate(char **args, intp *dimensions, intp *steps, void *func)
+U@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register 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@ *)op)=*((@typ@ *)i1);
+    BINARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        const u@type@ in2 = *(u@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((u@type@ *)op) = 0;
+        }
+        else {
+            *((u@type@ *)op)= in1/in2;
+        }
     }
 }
-/**end repeat**/
 
-/**begin repeat
-
-   #TYP=CFLOAT, CDOUBLE, CLONGDOUBLE#
-   #typ=float, double, longdouble#
-*/
 static void
-@TYP@_conjugate(char **args, intp *dimensions, intp *steps, void *func) {
-    register intp i, is1=steps[0], os=steps[1], n=dimensions[0];
-    char *i1=args[0], *op=args[1];
+@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@type@ *)op) = 0;
+        }
+        else {
+            /* handle mixed case the way Python does */
+            const @type@ rem = in1 % in2;
+            if ((in1 > 0) == (in2 > 0) || rem == 0) {
+                *((@type@ *)op) = rem;
+            }
+            else {
+                *((@type@ *)op) = rem + in2;
+            }
+        }
+    }
+}
 
-    for(i=0; i<n; i++, i1+=is1, op+=os) {
-        ((@typ@ *)op)[0]=((@typ@ *)i1)[0];
-        ((@typ@ *)op)[1]=-(((@typ@ *)i1)[1]);
+static void
+U@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        const u@type@ in2 = *(u@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@type@ *)op) = 0;
+        }
+        else {
+            *((@type@ *)op) = in1 % in2;
+        }
     }
 }
+
 /**end repeat**/
 
+/*
+ *****************************************************************************
+ **                             FLOAT LOOPS                                 **
+ *****************************************************************************
+ */
+
 
 /**begin repeat
+ * Float types
+ *  #type = float, double, longdouble#
+ *  #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
+ *  #c = f, , l#
+ *  #C = F, , L#
+ */
 
-   #TYPE=BOOL,UBYTE,USHORT,UINT,ULONG,ULONGLONG#
-   #typ=Bool, ubyte, ushort, uint, ulong, ulonglong#
-*/
+/**begin repeat1
+ * Arithmetic
+ * # kind = add, subtract, multiply, divide#
+ * # OP = +, -, *, /#
+ */
 static void
-@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, n;
-    intp is1=steps[0], os=steps[1];
-    char *i1=args[0], *op=args[1];
-
-    n=dimensions[0];
-
-    for(i=0; i<n; i++, i1+=is1, op+=os) {
-        *((@typ@ *)op) = *((@typ@*)i1);
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op) = in1 @OP@ in2;
     }
 }
-/**end repeat**/
-
-/**begin repeat
+/**end repeat1**/
 
-   #TYPE=BYTE,SHORT,INT,LONG,LONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
-   #typ=byte, short, int, long, longlong, float, double, longdouble#
-*/
+/**begin repeat1
+ * #kind = equal, not_equal, less, less_equal, greater, greater_equal,
+ *        logical_and, logical_or#
+ * #OP = ==, !=, <, <=, >, >=, &&, ||#
+ */
 static void
-@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, n;
-    intp is1=steps[0], os=steps[1];
-    char *i1=args[0], *op=args[1];
-
-    n=dimensions[0];
-    for(i=0; i<n; i++, i1+=is1, op+=os) {
-        *((@typ@ *)op) = *((@typ@ *)i1) < 0 ? -*((@typ@ *)i1) : *((@typ@ *)i1);
-        *((@typ@ *)op) += 0; /* clear sign-bit */
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((Bool *)op) = in1 @OP@ in2;
     }
 }
-/**end repeat**/
+/**end repeat1**/
 
-/**begin repeat
-   #TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
-   #typ= float, double, longdouble#
-   #c= f,,l#
-*/
 static void
-@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i, n;
-    register intp is1=steps[0], os=steps[1];
-    char *i1=args[0], *op=args[1];
-    n=dimensions[0];
-    for(i=0; i<n; i++, i1+=is1, op+=os) {
-        *((@typ@ *)op) = (@typ@)sqrt@c@(((@typ@ *)i1)[0]*((@typ@ *)i1)[0] + ((@typ@ *)i1)[1]*((@typ@ *)i1)[1]);
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((Bool *)op)= (in1 && !in2) || (!in1 && in2);
     }
 }
-/**end repeat**/
-
-/**begin repeat
 
-   #kind=greater, greater_equal, less, less_equal, equal, not_equal, logical_and, logical_or, bitwise_and, bitwise_or, bitwise_xor#
-   #OP=>, >=, <, <=, ==, !=, &&, ||, &, |, ^#
-**/
 static void
-BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    Bool in1, in2;
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        in1 = (*((Bool *)i1) != 0);
-        in2 = (*((Bool *)i2) != 0);
-        *((Bool *)op)= in1 @OP@ in2;
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((Bool *)op) = !in1;
     }
 }
-/**end repeat**/
-
-/**begin repeat
-
-   #TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*4#
-   #OP= >*13, >=*13, <*13, <=*13#
-   #typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*4#
-   #kind= greater*13, greater_equal*13, less*13, less_equal*13#
-*/
 
+/**begin repeat1
+ * #kind = isnan, isinf, isfinite, signbit#
+ * #func = isnan, isinf, isfinite, signbit#
+ **/
 static void
 @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((Bool *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((Bool *)op) = @func@(in1) != 0;
     }
 }
-/**end repeat**/
-
-
-/**begin repeat
-   #TYPE=(CFLOAT,CDOUBLE,CLONGDOUBLE)*4#
-   #OP= >*3, >=*3, <*3, <=*3#
-   #typ=(cfloat, cdouble, clongdouble)*4#
-   #kind= greater*3, greater_equal*3, less*3, less_equal*3#
-*/
+/**end repeat1**/
 
+/**begin repeat1
+ * #kind = maximum, minimum#
+ * #OP =  >, <#
+ **/
 static void
 @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        if (((@typ@ *)i1)->real == ((@typ@ *)i2)->real)
-            *((Bool *)op)=((@typ@ *)i1)->imag @OP@ \
-                ((@typ@ *)i2)->imag;
-        else
-            *((Bool *)op)=((@typ@ *)i1)->real @OP@ \
-                ((@typ@ *)i2)->real;
+    /*  */
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op) = in1 @OP@ in2 ? in1 : in2;
     }
 }
-/**end repeat**/
-
+/**end repeat1**/
 
-/**begin repeat
-   #TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*4#
-   #typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*4#
-   #OP= ==*13, !=*13, &&*13, ||*13#
-   #kind=equal*13, not_equal*13, logical_and*13, logical_or*13#
-*/
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((Bool *)op) = *((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op) = floor@c@(in1/in2);
     }
 }
-/**end repeat**/
 
-
-/**begin repeat
-
-   #TYPE=(CFLOAT, CDOUBLE, CLONGDOUBLE)*4#
-   #typ=(float, double, longdouble)*4#
-   #OP= ==*3, !=*3, &&*3, ||*3#
-   #OP2= &&*3, ||*3, &&*3, ||*3#
-   #kind=equal*3, not_equal*3, logical_and*3, logical_or*3#
-*/
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((Bool *)op) = (*((@typ@ *)i1) @OP@ *((@typ@ *)i2)) @OP2@ (*((@typ@ *)i1+1) @OP@ *((@typ@ *)i2+1));
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        const @type@ res = fmod@c@(in1,in2);
+        if (res && ((in2 < 0) != (res < 0))) {
+            *((@type@ *)op) = res + in2;
+        }
+        else {
+            *((@type@ *)op) = res;
+        }
     }
 }
-/**end repeat**/
-
 
-/** OBJECT comparison for OBJECT arrays **/
-
-/**begin repeat
+static void
+@TYPE@_square(char **args, intp *dimensions, intp *steps, void *data)
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op) = in1*in1;
+    }
+}
 
-   #kind=greater, greater_equal, less, less_equal, equal, not_equal#
-   #op=GT, GE, LT, LE, EQ, NE#
-*/
 static void
-OBJECT_@kind@(char **args, intp *dimensions, intp *steps, void *func) {
-    register intp i, is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((Bool *)op)=PyObject_RichCompareBool(*((PyObject **)i1),
-                                               *((PyObject **)i2),
-                                               Py_@op@);
+@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data)
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op) = 1.0/in1;
     }
 }
-/**end repeat**/
 
-/**begin repeat
+static void
+@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *data)
+{
+    OUTPUT_LOOP {
+        *((@type@ *)op) = 1;
+    }
+}
 
-   #TYPE=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
-   #typ=byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
-   #styp=byte, byte, short, short, int, int, long, long, longlong, longlong, float, double, longdouble#
-*/
 static void
-@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp 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@ *)op) = (@typ@) (- (@styp@)*((@typ@ *)i1));
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op) = in1;
     }
 }
-/**end repeat**/
 
-#define BOOL_negative BOOL_logical_not
+static void
+@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ tmp = (in1 > 0) ? in1 : -in1;
+        /* add 0 to clear -0.0 */
+        *((@type@ *)op) = tmp + 0;
+    }
+}
 
-#define _SIGN1(x) ((x) > 0 ? 1 : ((x) < 0 ? -1 : 0))
-#define _SIGN2(x) ((x) == 0 ? 0 : 1)
-#define _SIGNC(x) (((x).real > 0) ? 1 : ((x).real < 0 ? -1 : ((x).imag > 0 ? 1 : ((x).imag < 0) ? -1 : 0)))
-/**begin repeat
-   #TYPE=BYTE,SHORT,INT,LONG,LONGLONG,FLOAT,DOUBLE,LONGDOUBLE,UBYTE,USHORT,UINT,ULONG,ULONGLONG#
-   #typ=byte,short,int,long,longlong,float,double,longdouble,ubyte,ushort,uint,ulong,ulonglong#
-   #func=_SIGN1*8,_SIGN2*5#
-*/
 static void
-@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],os=steps[1], n=dimensions[0];
-    char *i1=args[0], *op=args[1];
-    @typ@ t1;
-    for(i=0; i<n; i++, i1+=is1, op+=os) {
-        t1 = *((@typ@ *)i1);
-        *((@typ@ *)op) = (@typ@) @func@(t1);
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op) = -in1;
     }
 }
-/**end repeat**/
 
-/**begin repeat
-   #TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
-   #typ=cfloat,cdouble,clongdouble#
-   #rtyp=float,double,longdouble#
-*/
 static void
 @TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],os=steps[1], n=dimensions[0];
-    char *i1=args[0], *op=args[1];
-    @typ@ t1;
-    for(i=0; i<n; i++, i1+=is1, op+=os) {
-        t1 = *((@typ@ *)i1);
-        (*((@typ@ *)op)).real = (@rtyp@)_SIGNC(t1);
-        (*((@typ@ *)op)).imag = (@rtyp@)0;
+    /*  */
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        if (in1 > 0) {
+            *((@type@ *)op) = 1;
+        }
+        else if (in1 < 0) {
+            *((@type@ *)op) = -1;
+        }
+        else {
+            *((@type@ *)op) = 0;
+        }
     }
 }
-/**end repeat**/
-
-#undef _SIGN1
-#undef _SIGN2
-#undef _SIGNC
-
 
 static void
-OBJECT_sign(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_modf(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],os=steps[1], n=dimensions[0];
-    char *i1=args[0], *op=args[1];
-    PyObject *t1, *zero, *res;
-    zero = PyInt_FromLong(0);
-    for(i=0; i<n; i++, i1+=is1, op+=os) {
-        t1 = *((PyObject **)i1);
-        res = PyInt_FromLong((long) PyObject_Compare(t1, zero));
-        *((PyObject **)op) = res;
+    UNARY_LOOP_TWO_OUT {
+        const @type@ in1 = *(@type@ *)ip1;
+        *(@type@ *)op1 = modf@c@(in1, (@type@ *)op2);
     }
-    Py_DECREF(zero);
 }
 
-
-/**begin repeat
-   #TYPE=BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
-   #typ=Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
-*/
+#define HAVE_DOUBLE_FUNCS
+#ifdef HAVE_@TYPE@_FUNCS
 static void
-@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp 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) {
-        *((Bool *)op) = ! *((@typ@ *)i1);
+    UNARY_LOOP_TWO_OUT {
+        const @type@ in1 = *(@type@ *)ip1;
+        *(@type@ *)op1 = frexp@c@(in1, (int *)op2);
     }
 }
-/**end repeat**/
 
-/**begin repeat
-   #TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
-   #typ=cfloat, cdouble, clongdouble#
-*/
 static void
-@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func)
+@TYPE@_ldexp(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp 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) {
-        *((Bool *)op) = ! (((@typ@ *)i1)->real ||       \
-                           ((@typ@ *)i1)->imag);
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const int in2 = *(int *)ip2;
+        *(@type@ *)op = ldexp@c@(in1, in2);
     }
 }
-/**end repeat**/
+#endif
+#undef HAVE_DOUBLE_FUNCS
 
+#define @TYPE@_true_divide @TYPE@_divide
 
+/**end repeat**/
 
 
+/*
+ *****************************************************************************
+ **                           COMPLEX LOOPS                                 **
+ *****************************************************************************
+ */
+
 /**begin repeat
-   #TYPE=BYTE,SHORT,INT,LONG,LONGLONG#
-   #typ=byte, short, int, long, longlong#
-   #c=f*2,,,l*1#
-*/
+ * complex types
+ * #ctype= cfloat, cdouble, clongdouble#
+ * #CTYPE= CFLOAT, CDOUBLE, CLONGDOUBLE#
+ * #type = float, double, longdouble#
+ * #c = f, , l#
+ */
+
+/**begin repeat1
+ * arithmetic
+ * #kind = add, subtract#
+ * #OP = +, -#
+ */
 static void
-@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    register @typ@ ix,iy, tmp;
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        ix = *((@typ@ *)i1);
-        iy = *((@typ@ *)i2);
-        if (iy == 0 || ix == 0) {
-            if (iy == 0) generate_divbyzero_error();
-            *((@typ@ *)op) = 0;
-        }
-        else if ((ix > 0) == (iy > 0)) {
-            *((@typ@ *)op) = ix % iy;
-        }
-        else {  /* handle mixed case the way Python does */
-            tmp = ix % iy;
-            if (tmp) tmp += iy;
-            *((@typ@ *)op)= tmp;
-        }
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        ((@type@ *)op)[0] = in1r @OP@ in2r;
+        ((@type@ *)op)[1] = in1i @OP@ in2i;
     }
 }
-/**end repeat**/
+/**end repeat1**/
 
-/**begin repeat
-   #TYPE=UBYTE,USHORT,UINT,ULONG,ULONGLONG#
-   #typ=ubyte, ushort, uint, ulong, ulonglong#
-*/
 static void
-@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_multiply(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    register @typ@ ix,iy;
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        ix = *((@typ@ *)i1);
-        iy = *((@typ@ *)i2);
-        if (iy == 0) {
-            generate_divbyzero_error();
-            *((@typ@ *)op) = 0;
-        }
-        *((@typ@ *)op) = ix % iy;
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        ((@type@ *)op)[0] = in1r*in2r - in1i*in2i;
+        ((@type@ *)op)[1] = in1r*in2i + in1i*in2r;
     }
 }
-/**end repeat**/
 
-/**begin repeat
-   #TYPE=FLOAT,DOUBLE,LONGDOUBLE#
-   #typ=float,double,longdouble#
-   #c=f,,l#
-*/
 static void
-@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_divide(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    @typ@ x, y, res;
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        x = *((@typ@ *)i1);
-        y = *((@typ@ *)i2);
-        res = fmod@c@(x, y);
-        if (res && ((y < 0) != (res < 0))) {
-            res += y;
-        }
-        *((@typ@ *)op)= res;
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        @type@ d = in2r*in2r + in2i*in2i;
+        ((@type@ *)op)[0] = (in1r*in2r + in1i*in2i)/d;
+        ((@type@ *)op)[1] = (in1i*in2r - in1r*in2i)/d;
     }
 }
-/**end repeat**/
 
+static void
+@CTYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        @type@ d = in2r*in2r + in2i*in2i;
+        ((@type@ *)op)[0] = floor@c@((in1r*in2r + in1i*in2i)/d);
+        ((@type@ *)op)[1] = 0;
+    }
+}
 
-/**begin repeat
-
-   #TYPE=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG#
-   #typ=byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong#
+/**begin repeat1
+   #kind = equal, not_equal#
+   #OP1 = ==, !=#
+   #OP2 = &&, ||#
 */
 static void
-@TYPE@_fmod(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    @typ@ x, y;
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        x = *((@typ@ *)i1);
-        y = *((@typ@ *)i2);
-        if (y == 0) {
-            generate_divbyzero_error();
-            *((@typ@ *)op) = 0;
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        *((Bool *)op) = (in1r @OP1@ in2r) @OP2@ (in1i @OP1@ in2i);
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind= greater, greater_equal, less, less_equal#
+ * #OP = >, >=, <, <=#
+ */
+static void
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        if (in1r != in2r) {
+            *((Bool *)op) = in1r @OP@ in2r ? 1 : 0;
         }
         else {
-            *((@typ@ *)op)= x % y;
+            *((Bool *)op) = in1i @OP@ in2i ? 1 : 0;
         }
-
     }
 }
-/**end repeat**/
-
-/**begin repeat
-
-   #TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG)*5#
-   #typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong)*5#
-   #OP= &*10, |*10, ^*10, <<*10, >>*10#
-   #kind=bitwise_and*10, bitwise_or*10, bitwise_xor*10, left_shift*10, right_shift*10#
+/**end repeat1**/
 
+/**begin repeat1
+   #kind = logical_and, logical_or#
+   #OP1 = ||, ||#
+   #OP2 = &&, ||#
 */
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    register char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((@typ@ *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        *((Bool *)op) = (in1r @OP1@ in1i) @OP2@ (in2r @OP1@ in2i);
     }
 }
-/**end repeat**/
-
+/**end repeat1**/
 
-/**begin repeat
-   #TYPE=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG#
-   #typ=byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong#
-*/
 static void
-@TYPE@_invert(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp 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@ *)op) = ~ *((@typ@*)i1);
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        const Bool tmp1 = (in1r || in1i);
+        const Bool tmp2 = (in2r || in2i);
+        *((Bool *)op) = (tmp1 && !tmp2) || (!tmp1 && tmp2);
     }
 }
-/**end repeat**/
 
 static void
-BOOL_invert(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp 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) {
-        *((Bool *)op) = (*((Bool *)i1) ? FALSE : TRUE);
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        *((Bool *)op) = !(in1r || in1i);
     }
 }
 
-
-/**begin repeat
-   #TYPE=BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
-   #typ=Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
-
-*/
+/**begin repeat1
+ * #kind = isnan, isinf, isfinite#
+ * #func = isnan, isinf, isfinite#
+ * #OP = ||, ||, &&#
+ **/
 static void
-@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((Bool *)op)=(*((@typ@ *)i1) || *((@typ@ *)i2)) && !(*((@typ@ *)i1) && *((@typ@ *)i2));
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        *((Bool *)op) = @func@(in1r) @OP@ @func@(in1i);
     }
 }
-/**end repeat**/
-
+/**end repeat1**/
 
-/**begin repeat
-   #TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
-   #typ=cfloat, cdouble, clongdouble#
-*/
 static void
-@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_square(char **args, intp *dimensions, intp *steps, void *data)
 {
-    Bool p1, p2;
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        p1 = ((@typ@ *)i1)->real || ((@typ@ *)i1)->imag;
-        p2 = ((@typ@ *)i2)->real || ((@typ@ *)i2)->imag;
-        *((Bool *)op)= (p1 || p2) && !(p1 && p2);
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        ((@type@ *)op)[0] = in1r*in1r - in1i*in1i;
+        ((@type@ *)op)[1] = in1r*in1i + in1i*in1r;
     }
 }
-/**end repeat**/
-
 
-
-/**begin repeat
-
-   #TYPE=(BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*2#
-   #OP= >*14, <*14#
-   #typ=(Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*2#
-   #kind= maximum*14, minimum*14#
-*/
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        *((@typ@ *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2) ? *((@typ@ *)i1) : *((@typ@ *)i2);
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        if (fabs@c@(in1i) <= fabs@c@(in1r)) {
+            const @type@ r = in1i/in1r;
+            const @type@ d = in1r + in1i*r;
+            ((@type@ *)op)[0] = 1/d;
+            ((@type@ *)op)[1] = -r/d;
+        } else {
+            const @type@ r = in1r/in1i;
+            const @type@ d = in1r*r + in1i;
+            ((@type@ *)op)[0] = r/d;
+            ((@type@ *)op)[1] = -1/d;
+        }
     }
 }
-/**end repeat**/
-
-/**begin repeat
 
-   #TYPE=(CFLOAT,CDOUBLE,CLONGDOUBLE)*2#
-   #OP= >*3, <*3#
-   #typ=(cfloat, cdouble, clongdouble)*2#
-   #kind= maximum*3, minimum*3#
-*/
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *data)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    @typ@ *i1c, *i2c;
-    for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        i1c = (@typ@ *)i1;
-        i2c = (@typ@ *)i2;
-        if ((i1c->real @OP@ i2c->real) ||                       \
-            ((i1c->real==i2c->real) && (i1c->imag @OP@ i2c->imag)))
-            memmove(op, i1, sizeof(@typ@));
-        else
-            memmove(op, i2, sizeof(@typ@));
+    OUTPUT_LOOP {
+        ((@type@ *)op)[0] = 1;
+        ((@type@ *)op)[1] = 0;
     }
 }
-/**end repeat**/
-
 
+static void
+@CTYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func) {
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        ((@type@ *)op)[0] = in1r;
+        ((@type@ *)op)[1] = -in1i;
+    }
+}
 
-/*** isinf, isinf, isfinite, signbit ***/
-/**begin repeat
-   #kind=isnan*3, isinf*3, isfinite*3, signbit*3#
-   #TYPE=(FLOAT, DOUBLE, LONGDOUBLE)*4#
-   #typ=(float, double, longdouble)*4#
-   #c=(f,,l)*4#
-*/
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is=steps[0], os=steps[1], n=dimensions[0];
-    char *ip=args[0], *op=args[1];
-    for(i=0; i<n; i++, ip+=is, op+=os) {
-        *((Bool *)op) = (Bool) (@kind@@c@(*((@typ@ *)ip)) != 0);
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        *((@type@ *)op) = sqrt@c@(in1r*in1r + in1i*in1i);
     }
 }
-/**end repeat**/
-
 
-/**begin repeat
-   #kind=isnan*3, isinf*3, isfinite*3#
-   #TYPE=(CFLOAT, CDOUBLE, CLONGDOUBLE)*3#
-   #typ=(float, double, longdouble)*3#
-   #c=(f,,l)*3#
-   #OP=||*6,&&*3#
-*/
 static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_sign(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is=steps[0], os=steps[1], n=dimensions[0];
-    char *ip=args[0], *op=args[1];
-    for(i=0; i<n; i++, ip+=is, op+=os) {
-        *((Bool *)op) = @kind@@c@(((@typ@ *)ip)[0]) @OP@        \
-            @kind@@c@(((@typ@ *)ip)[1]);
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        if (in1r > 0) {
+            ((@type@ *)op)[0] = 1;
+        }
+        else if (in1r < 0) {
+            ((@type@ *)op)[0] = -1;
+        }
+        else {
+            if (in1i > 0) {
+                ((@type@ *)op)[0] = 1;
+            }
+            else if (in1i < 0) {
+                ((@type@ *)op)[0] = -1;
+            }
+            else {
+                ((@type@ *)op)[0] = 0;
+            }
+        }
+        ((@type@ *)op)[1] = 0;
     }
 }
-/**end repeat**/
-
-
-
-
-/****** modf ****/
 
-/**begin repeat
-   #TYPE=FLOAT, DOUBLE, LONGDOUBLE#
-   #typ=float, double, longdouble#
-   #c=f,,l#
-*/
+/**begin repeat1
+ * #kind = maximum, minimum#
+ * #OP = >, <#
+ */
 static void
-@TYPE@_modf(char **args, intp *dimensions, intp *steps, void *func)
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],os1=steps[1],os2=steps[2],n=dimensions[0];
-    char *i1=args[0], *op1=args[1], *op2=args[2];
-    @typ@ x1, y1, y2;
-    for (i=0; i<n; i++, i1+=is1, op1+=os1, op2+=os2) {
-        x1 = *((@typ@ *)i1);
-        y1 = modf@c@(x1, &y2);
-        *((@typ@ *)op1) = y1;
-        *((@typ@ *)op2) = y2;
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        if (in1r @OP@ in2r || ((in1r == in2r) && (in1i @OP@ in2i))) {
+            ((@type@ *)op)[0] = in1r;
+            ((@type@ *)op)[1] = in1i;
+        }
+        else {
+            ((@type@ *)op)[0] = in2r;
+            ((@type@ *)op)[1] = in2i;
+        }
     }
 }
+/**end repeat1**/
+
+#define @CTYPE@_true_divide @CTYPE@_divide
 /**end repeat**/
 
-#define HAVE_DOUBLE_FUNCS
+
+/*
+ *****************************************************************************
+ **                            OBJECT LOOPS                                 **
+ *****************************************************************************
+ */
+
 /**begin repeat
-   #TYPE=FLOAT, DOUBLE, LONGDOUBLE#
-   #typ=float, double, longdouble#
-   #c=f,,l#
-*/
-#ifdef HAVE_@TYPE@_FUNCS
+ * #kind = equal, not_equal, greater, greater_equal, less, less_equal#
+ * #OP = EQ, NE, GT, GE, LT, LE#
+ */
 static void
-@TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *func)
-{
-    register intp i;
-    intp is1=steps[0],os1=steps[1],os2=steps[2],n=dimensions[0];
-    char *i1=args[0], *op1=args[1], *op2=args[2];
-    @typ@ x1, y1;
-    int y2;
-    for (i=0; i<n; i++, i1+=is1, op1+=os1, op2+=os2) {
-        x1 = *((@typ@ *)i1);
-        y1 = frexp@c@(x1, &y2);
-        *((@typ@ *)op1) = y1;
-        *((int *) op2) = y2;
+OBJECT_@kind@(char **args, intp *dimensions, intp *steps, void *func) {
+    BINARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        PyObject *in2 = *(PyObject **)ip2;
+        *(Bool *)op = (Bool) PyObject_RichCompareBool(in1, in2, Py_@OP@);
     }
 }
+/**end repeat**/
 
 static void
-@TYPE@_ldexp(char **args, intp *dimensions, intp *steps, void *func)
+OBJECT_sign(char **args, intp *dimensions, intp *steps, void *func)
 {
-    register intp i;
-    intp is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
-    char *i1=args[0], *i2=args[1], *op=args[2];
-    @typ@ x1, y1;
-    int x2;
-    for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
-        x1 = *((@typ@ *)i1);
-        x2 = *((int *)i2);
-        y1 = ldexp@c@(x1, x2);
-        *((@typ@ *)op) = y1;
+    PyObject *zero = PyInt_FromLong(0);
+    UNARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        *(PyObject **)op = PyInt_FromLong(PyObject_Compare(in1, zero));
     }
+    Py_DECREF(zero);
 }
-#endif
-/**end repeat**/
-#undef HAVE_DOUBLE_FUNCS
+
+/*
+ *****************************************************************************
+ **                              END LOOPS                                  **
+ *****************************************************************************
+ */
 
 
 static PyUFuncGenericFunction frexp_functions[] = {
@@ -2194,12 +2097,8 @@ static char ldexp_signatures[] = {
 };
 
 
-
 #include "__umath_generated.c"
-
-
 #include "ufuncobject.c"
-
 #include "__ufunc_api.c"
 
 static double
@@ -2370,6 +2269,7 @@ PyMODINIT_FUNC initumath(void) {
     PyDict_SetItemString(d, "mod", s2);
 
     return;
+
  err:
     /* Check for errors */
     if (!PyErr_Occurred()) {