Merge branch 'ufunc'

6 files changed, 2001 insertions, 1959 deletions

diff --git a/numpy/core/SConscript b/numpy/core/SConscript
index 517bbfdcd..0ffbb444a 100644
--- a/numpy/core/SConscript
+++ b/numpy/core/SConscript
@@ -137,9 +137,9 @@ mfuncs = ('expl', 'expf', 'log1p', 'expm1', 'asinh', 'atanhf', 'atanhl',
 mfuncs_defined = dict([(f, 0) for f in mfuncs])
 
 # Check for mandatory funcs: we barf if a single one of those is not there
-mandatory_funcs = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs",
-"floor", "ceil", "sqrt", "log10", "log", "exp", "asin", "acos", "atan", "fmod",
-'modf', 'frexp', 'ldexp']
+    mandatory_funcs = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs",
+		"floor", "ceil", "sqrt", "log10", "log", "exp", "asin",
+		"acos", "atan", "fmod", 'modf', 'frexp', 'ldexp']
 
 if not config.CheckFuncsAtOnce(mandatory_funcs):
     raise SystemError("One of the required function to build numpy is not"
@@ -159,16 +159,17 @@ def check_funcs(funcs):
 
 # XXX: we do not test for hypot because python checks for it (HAVE_HYPOT in
 # python.h... I wish they would clean their public headers someday)
-optional_stdfuncs = ["expm1", "log1p", "acosh", "asinh", "atanh",
-                     "rint", "trunc"]
+    optional_stdfuncs = ["expm1", "log1p", "acosh", "asinh", "atanh",
+                         "rint", "trunc", "exp2", "log2"]
 
 check_funcs(optional_stdfuncs)
 
 # C99 functions: float and long double versions
-c99_funcs = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs", "floor",
-             "ceil", "rint", "trunc", "sqrt", "log10", "log", "exp",
-             "expm1", "asin", "acos", "atan", "asinh", "acosh", "atanh",
-             "hypot", "atan2", "pow", "fmod", "modf", 'frexp', 'ldexp']
+    c99_funcs = ["sin", "cos", "tan", "sinh", "cosh", "tanh", "fabs", "floor",
+                 "ceil", "rint", "trunc", "sqrt", "log10", "log", "log1p", "exp",
+                 "expm1", "asin", "acos", "atan", "asinh", "acosh", "atanh",
+                 "hypot", "atan2", "pow", "fmod", "modf", 'frexp', 'ldexp',
+                 "exp2", "log2"]
 
 for prec in ['l', 'f']:
     fns = [f + prec for f in c99_funcs]
@@ -251,7 +252,9 @@ env.Append(BUILDERS = {'GenerateMultiarrayApi' : array_api_gen_bld,
 # Generate generated code
 #------------------------
 scalartypes_src = env.GenerateFromTemplate(pjoin('src', 'scalartypes.inc.src'))
-math_c99_src = env.GenerateFromTemplate(pjoin('src', 'umath_funcs_c99.inc.src'))
+umath_funcs_c99_src = env.GenerateFromTemplate(pjoin('src', 'umath_funcs_c99.inc.src'))
+umath_funcs_src = env.GenerateFromTemplate(pjoin('src', 'umath_funcs.inc.src'))
+umath_loops_src = env.GenerateFromTemplate(pjoin('src', 'umath_loops.inc.src'))
 arraytypes_src = env.GenerateFromTemplate(pjoin('src', 'arraytypes.inc.src'))
 sortmodule_src = env.GenerateFromTemplate(pjoin('src', '_sortmodule.c.src'))
 umathmodule_src = env.GenerateFromTemplate(pjoin('src', 'umathmodule.c.src'))
diff --git a/numpy/core/code_generators/genapi.py b/numpy/core/code_generators/genapi.py
index 59789b3bb..11c88b2cb 100644
--- a/numpy/core/code_generators/genapi.py
+++ b/numpy/core/code_generators/genapi.py
@@ -18,6 +18,7 @@ API_FILES = ['arraymethods.c',
              'multiarraymodule.c',
              'scalartypes.inc.src',
              'umath_ufunc_object.inc',
+             'umath_funcs.inc.src',
              'umathmodule.c.src'
             ]
 THIS_DIR = os.path.dirname(__file__)
diff --git a/numpy/core/setup.py b/numpy/core/setup.py
index 463d1dc29..28e1eded3 100644
--- a/numpy/core/setup.py
+++ b/numpy/core/setup.py
@@ -358,6 +358,8 @@ def configuration(parent_package='',top_path=None):
                                     join('src','scalartypes.inc.src'),
                                     join('src','arraytypes.inc.src'),
                                     join('src','umath_funcs_c99.inc.src'),
+                                    join('src','umath_funcs.inc.src'),
+                                    join('src','umath_loops.inc.src'),
                                     ],
                          depends = [join('src','umath_ufunc_object.inc'),
                                     generate_umath_py,
diff --git a/numpy/core/src/umath_funcs.inc.src b/numpy/core/src/umath_funcs.inc.src
new file mode 100644
index 000000000..1e5337afe
--- /dev/null
+++ b/numpy/core/src/umath_funcs.inc.src
@@ -0,0 +1,570 @@
+/*
+ * This file is for the definitions of the non-c99 functions used in ufuncs.
+ * All the complex ufuncs are defined here along with a smattering of real and
+ * object functions.
+ */
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846264338328
+#endif
+
+#define M_LOG10_E       0.434294481903251827651128918916605082294397
+
+/* Useful constants in three precisions.*/
+
+/**begin repeat
+ * #c = f, ,l#
+ * #C = F, ,L#
+ */
+#define NPY_E@c@        2.7182818284590452353602874713526625@C@ /* e */
+#define NPY_LOG2E@c@    1.4426950408889634073599246810018921@C@ /* log_2 e */
+#define NPY_LOG10E@c@   0.4342944819032518276511289189166051@C@ /* log_10 e */
+#define NPY_LOGE2@c@    0.6931471805599453094172321214581766@C@ /* log_e 2 */
+#define NPY_LOGE10@c@   2.3025850929940456840179914546843642@C@ /* log_e 10 */
+#define NPY_PI@c@       3.1415926535897932384626433832795029@C@ /* pi */
+#define NPY_PI_2@c@     1.5707963267948966192313216916397514@C@ /* pi/2 */
+#define NPY_PI_4@c@     0.7853981633974483096156608458198757@C@ /* pi/4 */
+#define NPY_1_PI@c@     0.3183098861837906715377675267450287@C@ /* 1/pi */
+#define NPY_2_PI@c@     0.6366197723675813430755350534900574@C@ /* 2/pi */
+/**end repeat**/
+
+/*
+ ******************************************************************************
+ **                            FLOAT FUNCTIONS                               **
+ ******************************************************************************
+ */
+
+/**begin repeat
+ * #type = float, double, longdouble#
+ * #c = f, ,l#
+ * #C = F, ,L#
+ */
+
+#define LOGE2    NPY_LOGE2@c@
+#define LOG2E    NPY_LOG2E@c@
+#define RAD2DEG  (180.0@c@/NPY_PI@c@)
+#define DEG2RAD  (NPY_PI@c@/180.0@c@)
+
+static @type@
+rad2deg@c@(@type@ x) {
+    return x*RAD2DEG;
+}
+
+static @type@
+deg2rad@c@(@type@ x) {
+    return x*DEG2RAD;
+}
+
+static @type@
+log2_1p@c@(@type@ x)
+{
+    @type@ u = 1 + x;
+    if (u == 1) {
+        return LOG2E*x;
+    } else {
+        return log2@c@(u) * x / (u - 1);
+    }
+}
+
+static @type@
+exp2_1m@c@(@type@ x)
+{
+    @type@ u = exp@c@(x);
+    if (u == 1.0) {
+        return LOGE2*x;
+    } else if (u - 1 == -1) {
+        return -LOGE2;
+    } else {
+        return (u - 1) * x/log2@c@(u);
+    }
+}
+
+static @type@
+logaddexp@c@(@type@ x, @type@ y)
+{
+    const @type@ tmp = x - y;
+    if (tmp > 0) {
+        return x + log1p@c@(exp@c@(-tmp));
+    }
+    else {
+        return y + log1p@c@(exp@c@(tmp));
+    }
+}
+
+static @type@
+logaddexp2@c@(@type@ x, @type@ y)
+{
+    const @type@ tmp = x - y;
+    if (tmp > 0) {
+        return x + log2_1p@c@(exp2@c@(-tmp));
+    }
+    else {
+        return y + log2_1p@c@(exp2@c@(tmp));
+    }
+}
+
+#define degrees@c@ rad2deg@c@
+#define radians@c@ deg2rad@c@
+
+#undef LOGE2
+#undef LOG2E
+#undef RAD2DEG
+#undef DEG2RAD
+
+/**end repeat**/
+
+/*
+ *****************************************************************************
+ **                        PYTHON OBJECT FUNCTIONS                          **
+ *****************************************************************************
+ */
+
+static PyObject *
+Py_square(PyObject *o)
+{
+    return PyNumber_Multiply(o, o);
+}
+
+static PyObject *
+Py_get_one(PyObject *NPY_UNUSED(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;
+}
+
+/*
+ * Define numpy version of PyNumber_Power as binary function.
+ */
+static PyObject *
+npy_ObjectPower(PyObject *x, PyObject *y)
+{
+    return PyNumber_Power(x, y, Py_None);
+}
+
+/**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**/
+
+
+/*
+ *****************************************************************************
+ **                           COMPLEX FUNCTIONS                             **
+ *****************************************************************************
+ */
+
+
+/*
+ * Don't pass structures between functions (only pointers) because how
+ * structures are passed is compiler dependent and could cause segfaults if
+ * umath_ufunc_object.inc is compiled with a different compiler than an
+ * extension that makes use of the UFUNC API
+ */
+
+/**begin repeat
+
+   #typ=float, double, longdouble#
+   #c=f,,l#
+*/
+
+/* constants */
+static c@typ@ nc_1@c@ = {1., 0.};
+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 void
+nc_sum@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+    r->real = a->real + b->real;
+    r->imag = a->imag + b->imag;
+    return;
+}
+
+static void
+nc_diff@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+    r->real = a->real - b->real;
+    r->imag = a->imag - b->imag;
+    return;
+}
+
+static void
+nc_neg@c@(c@typ@ *a, c@typ@ *r)
+{
+    r->real = -a->real;
+    r->imag = -a->imag;
+    return;
+}
+
+static void
+nc_prod@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+    @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;
+}
+
+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;
+    r->imag = (ai*br - ar*bi)/d;
+    return;
+}
+
+static void
+nc_sqrt@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ s,d;
+    if (x->real == 0. && x->imag == 0.)
+        *r = *x;
+    else {
+        s = sqrt@c@((fabs@c@(x->real) + hypot@c@(x->real,x->imag))/2);
+        d = x->imag/(2*s);
+        if (x->real > 0) {
+            r->real = s;
+            r->imag = d;
+        }
+        else if (x->imag >= 0) {
+            r->real = d;
+            r->imag = s;
+        }
+        else {
+            r->real = -d;
+            r->imag = -s;
+        }
+    }
+    return;
+}
+
+static void
+nc_rint@c@(c@typ@ *x, c@typ@ *r)
+{
+    r->real = rint@c@(x->real);
+    r->imag = rint@c@(x->imag);
+}
+
+static void
+nc_log@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ l = hypot@c@(x->real,x->imag);
+    r->imag = atan2@c@(x->imag, x->real);
+    r->real = log@c@(l);
+    return;
+}
+
+static void
+nc_log1p@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ l = hypot@c@(x->real + 1,x->imag);
+    r->imag = atan2@c@(x->imag, x->real + 1);
+    r->real = log@c@(l);
+    return;
+}
+
+static void
+nc_exp@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ a = exp@c@(x->real);
+    r->real = a*cos@c@(x->imag);
+    r->imag = a*sin@c@(x->imag);
+    return;
+}
+
+static void
+nc_expm1@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ a = exp@c@(x->real);
+    r->real = a*cos@c@(x->imag) - 1;
+    r->imag = a*sin@c@(x->imag);
+    return;
+}
+
+static void
+nc_pow@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+    intp n;
+    @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+
+    if (br == 0. && bi == 0.) {
+        r->real = 1.;
+        r->imag = 0.;
+        return;
+    }
+    if (ar == 0. && ai == 0.) {
+        r->real = 0.;
+        r->imag = 0.;
+        return;
+    }
+    if (bi == 0 && (n=(intp)br) == br) {
+        if (n > -100 && n < 100) {
+            c@typ@ p, aa;
+            intp mask = 1;
+            if (n < 0) n = -n;
+            aa = nc_1@c@;
+            p.real = ar; p.imag = ai;
+            while (1) {
+                if (n & mask)
+                    nc_prod@c@(&aa,&p,&aa);
+                mask <<= 1;
+                if (n < mask || mask <= 0) break;
+                nc_prod@c@(&p,&p,&p);
+            }
+            r->real = aa.real; r->imag = aa.imag;
+            if (br < 0) nc_quot@c@(&nc_1@c@, r, r);
+            return;
+        }
+    }
+    /*
+     * 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);
+    return;
+}
+
+
+static void
+nc_prodi@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ xr = x->real;
+    r->real = -x->imag;
+    r->imag = xr;
+    return;
+}
+
+
+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);
+    nc_prodi@c@(r, r);
+    nc_sum@c@(x, r, r);
+    nc_log@c@(r, r);
+    nc_prodi@c@(r, r);
+    nc_neg@c@(r, r);
+    return;
+}
+
+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);
+    nc_sqrt@c@(&t, &t);
+    nc_diff@c@(x, &nc_1@c@, r);
+    nc_sqrt@c@(r, r);
+    nc_prod@c@(&t, r, r);
+    nc_sum@c@(x, r, r);
+    nc_log@c@(r, r);
+    return;
+}
+
+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);
+    nc_sqrt@c@(r, r);
+    nc_prodi@c@(x, pa);
+    nc_sum@c@(pa, r, r);
+    nc_log@c@(r, r);
+    nc_prodi@c@(r, r);
+    nc_neg@c@(r, r);
+    return;
+}
+
+
+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;
+}
+
+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);
+    nc_quot@c@(r, pa, r);
+    nc_log@c@(r,r);
+    nc_prod@c@(&nc_i2@c@, r, r);
+    return;
+}
+
+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);
+    nc_quot@c@(pa, r, r);
+    nc_log@c@(r, r);
+    nc_prod@c@(&nc_half@c@, r, r);
+    return;
+}
+
+static void
+nc_cos@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ xr=x->real, xi=x->imag;
+    r->real = cos@c@(xr)*cosh@c@(xi);
+    r->imag = -sin@c@(xr)*sinh@c@(xi);
+    return;
+}
+
+static void
+nc_cosh@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ xr=x->real, xi=x->imag;
+    r->real = cos@c@(xi)*cosh@c@(xr);
+    r->imag = sin@c@(xi)*sinh@c@(xr);
+    return;
+}
+
+static void
+nc_log10@c@(c@typ@ *x, c@typ@ *r)
+{
+    nc_log@c@(x, r);
+    r->real *= (@typ@) M_LOG10_E;
+    r->imag *= (@typ@) M_LOG10_E;
+    return;
+}
+
+static void
+nc_sin@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ xr=x->real, xi=x->imag;
+    r->real = sin@c@(xr)*cosh@c@(xi);
+    r->imag = cos@c@(xr)*sinh@c@(xi);
+    return;
+}
+
+static void
+nc_sinh@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ xr=x->real, xi=x->imag;
+    r->real = cos@c@(xi)*sinh@c@(xr);
+    r->imag = sin@c@(xi)*cosh@c@(xr);
+    return;
+}
+
+static void
+nc_tan@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ sr,cr,shi,chi;
+    @typ@ rs,is,rc,ic;
+    @typ@ d;
+    @typ@ xr=x->real, xi=x->imag;
+    sr = sin@c@(xr);
+    cr = cos@c@(xr);
+    shi = sinh@c@(xi);
+    chi = cosh@c@(xi);
+    rs = sr*chi;
+    is = cr*shi;
+    rc = cr*chi;
+    ic = -sr*shi;
+    d = rc*rc + ic*ic;
+    r->real = (rs*rc+is*ic)/d;
+    r->imag = (is*rc-rs*ic)/d;
+    return;
+}
+
+static void
+nc_tanh@c@(c@typ@ *x, c@typ@ *r)
+{
+    @typ@ si,ci,shr,chr;
+    @typ@ rs,is,rc,ic;
+    @typ@ d;
+    @typ@ xr=x->real, xi=x->imag;
+    si = sin@c@(xi);
+    ci = cos@c@(xi);
+    shr = sinh@c@(xr);
+    chr = cosh@c@(xr);
+    rs = ci*shr;
+    is = si*chr;
+    rc = ci*chr;
+    ic = si*shr;
+    d = rc*rc + ic*ic;
+    r->real = (rs*rc+is*ic)/d;
+    r->imag = (is*rc-rs*ic)/d;
+    return;
+}
+
+/**end repeat**/
+
diff --git a/numpy/core/src/umath_loops.inc.src b/numpy/core/src/umath_loops.inc.src
new file mode 100644
index 000000000..c754cef1f
--- /dev/null
+++ b/numpy/core/src/umath_loops.inc.src
@@ -0,0 +1,1369 @@
+/* -*- c -*- */
+
+/*
+ *****************************************************************************
+ **                             UFUNC LOOPS                                 **
+ *****************************************************************************
+ */
+
+#define OUTPUT_LOOP\
+    char *op1 = args[1];\
+    intp os1 = steps[1];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, op1 += os1)
+
+#define UNARY_LOOP\
+    char *ip1 = args[0], *op1 = args[1];\
+    intp is1 = steps[0], os1 = steps[1];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, ip1 += is1, op1 += os1)
+
+#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], *op1 = args[2];\
+    intp is1 = steps[0], is2 = steps[1], os1 = steps[2];\
+    intp n = dimensions[0];\
+    intp i;\
+    for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1)
+
+#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)
+
+/******************************************************************************
+ **                          GENERIC FLOAT LOOPS                             **
+ *****************************************************************************/
+
+
+typedef float floatUnaryFunc(float x);
+typedef double doubleUnaryFunc(double x);
+typedef longdouble longdoubleUnaryFunc(longdouble x);
+typedef float floatBinaryFunc(float x, float y);
+typedef double doubleBinaryFunc(double x, double y);
+typedef longdouble longdoubleBinaryFunc(longdouble x, longdouble y);
+
+
+/*UFUNC_API*/
+static void
+PyUFunc_f_f(char **args, intp *dimensions, intp *steps, void *func)
+{
+    floatUnaryFunc *f = (floatUnaryFunc *)func;
+    UNARY_LOOP {
+        const float in1 = *(float *)ip1;
+        *(float *)op1 = f(in1);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_f_f_As_d_d(char **args, intp *dimensions, intp *steps, void *func)
+{
+    doubleUnaryFunc *f = (doubleUnaryFunc *)func;
+    UNARY_LOOP {
+        const float in1 = *(float *)ip1;
+        *(float *)op1 = (float)f((double)in1);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_ff_f(char **args, intp *dimensions, intp *steps, void *func)
+{
+    floatBinaryFunc *f = (floatBinaryFunc *)func;
+    BINARY_LOOP {
+        float in1 = *(float *)ip1;
+        float in2 = *(float *)ip2;
+        *(float *)op1 = f(in1, in2);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_ff_f_As_dd_d(char **args, intp *dimensions, intp *steps, void *func)
+{
+    doubleBinaryFunc *f = (doubleBinaryFunc *)func;
+    BINARY_LOOP {
+        float in1 = *(float *)ip1;
+        float in2 = *(float *)ip2;
+        *(float *)op1 = (double)f((double)in1, (double)in2);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_d_d(char **args, intp *dimensions, intp *steps, void *func)
+{
+    doubleUnaryFunc *f = (doubleUnaryFunc *)func;
+    UNARY_LOOP {
+        double in1 = *(double *)ip1;
+        *(double *)op1 = f(in1);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_dd_d(char **args, intp *dimensions, intp *steps, void *func)
+{
+    doubleBinaryFunc *f = (doubleBinaryFunc *)func;
+    BINARY_LOOP {
+        double in1 = *(double *)ip1;
+        double in2 = *(double *)ip2;
+        *(double *)op1 = f(in1, in2);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_g_g(char **args, intp *dimensions, intp *steps, void *func)
+{
+    longdoubleUnaryFunc *f = (longdoubleUnaryFunc *)func;
+    UNARY_LOOP {
+        longdouble in1 = *(longdouble *)ip1;
+        *(longdouble *)op1 = f(in1);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_gg_g(char **args, intp *dimensions, intp *steps, void *func)
+{
+    longdoubleBinaryFunc *f = (longdoubleBinaryFunc *)func;
+    BINARY_LOOP {
+        longdouble in1 = *(longdouble *)ip1;
+        longdouble in2 = *(longdouble *)ip2;
+        *(longdouble *)op1 = f(in1, in2);
+    }
+}
+
+
+
+/******************************************************************************
+ **                          GENERIC COMPLEX LOOPS                           **
+ *****************************************************************************/
+
+
+typedef void cdoubleUnaryFunc(cdouble *x, cdouble *r);
+typedef void cfloatUnaryFunc(cfloat *x, cfloat *r);
+typedef void clongdoubleUnaryFunc(clongdouble *x, clongdouble *r);
+typedef void cdoubleBinaryFunc(cdouble *x, cdouble *y, cdouble *r);
+typedef void cfloatBinaryFunc(cfloat *x, cfloat *y, cfloat *r);
+typedef void clongdoubleBinaryFunc(clongdouble *x, clongdouble *y,
+                                   clongdouble *r);
+
+/*UFUNC_API*/
+static void
+PyUFunc_F_F(char **args, intp *dimensions, intp *steps, void *func)
+{
+    cfloatUnaryFunc *f = (cfloatUnaryFunc *)func;
+    UNARY_LOOP {
+        cfloat in1 = *(cfloat *)ip1;
+        cfloat *out = (cfloat *)op1;
+        f(&in1, out);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_F_F_As_D_D(char **args, intp *dimensions, intp *steps, void *func)
+{
+    cdoubleUnaryFunc *f = (cdoubleUnaryFunc *)func;
+    UNARY_LOOP {
+        const float *in1 = (float *)ip1;
+        cdouble tmp = {(double)(in1[0]),(double)in1[1]};
+        cdouble out;
+        f(&tmp, &out);
+        ((float *)op1)[0] = (float)out.real;
+        ((float *)op1)[1] = (float)out.imag;
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_FF_F(char **args, intp *dimensions, intp *steps, void *func)
+{
+    cfloatBinaryFunc *f = (cfloatBinaryFunc *)func;
+    BINARY_LOOP {
+        cfloat in1 = *(cfloat *)ip1;
+        cfloat in2 = *(cfloat *)ip2;
+        cfloat *out = (cfloat *)op1;
+        f(&in1, &in2, out);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_FF_F_As_DD_D(char **args, intp *dimensions, intp *steps, void *func)
+{
+    cdoubleBinaryFunc *f = (cdoubleBinaryFunc *)func;
+    BINARY_LOOP {
+        const float *in1 = (float *)ip1;
+        const float *in2 = (float *)ip2;
+        cdouble tmp1 = {(double)(in1[0]),(double)in1[1]};
+        cdouble tmp2 = {(double)(in2[0]),(double)in2[1]};
+        cdouble out;
+        f(&tmp1, &tmp2, &out);
+        ((float *)op1)[0] = (float)out.real;
+        ((float *)op1)[1] = (float)out.imag;
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_D_D(char **args, intp *dimensions, intp *steps, void *func)
+{
+    cdoubleUnaryFunc *f = (cdoubleUnaryFunc *)func;
+    UNARY_LOOP {
+        cdouble in1 = *(cdouble *)ip1;
+        cdouble *out = (cdouble *)op1;
+        f(&in1, out);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_DD_D(char **args, intp *dimensions, intp *steps, void *func)
+{
+    cdoubleBinaryFunc *f = (cdoubleBinaryFunc *)func;
+    BINARY_LOOP {
+        cdouble in1 = *(cdouble *)ip1;
+        cdouble in2 = *(cdouble *)ip2;
+        cdouble *out = (cdouble *)op1;
+        f(&in1, &in2, out);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_G_G(char **args, intp *dimensions, intp *steps, void *func)
+{
+    clongdoubleUnaryFunc *f = (clongdoubleUnaryFunc *)func;
+    UNARY_LOOP {
+        clongdouble in1 = *(clongdouble *)ip1;
+        clongdouble *out = (clongdouble *)op1;
+        f(&in1, out);
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_GG_G(char **args, intp *dimensions, intp *steps, void *func)
+{
+    clongdoubleBinaryFunc *f = (clongdoubleBinaryFunc *)func;
+    BINARY_LOOP {
+        clongdouble in1 = *(clongdouble *)ip1;
+        clongdouble in2 = *(clongdouble *)ip2;
+        clongdouble *out = (clongdouble *)op1;
+        f(&in1, &in2, out);
+    }
+}
+
+
+/******************************************************************************
+ **                         GENERIC OBJECT lOOPS                             **
+ *****************************************************************************/
+
+/*UFUNC_API*/
+static void
+PyUFunc_O_O(char **args, intp *dimensions, intp *steps, void *func)
+{
+    unaryfunc f = (unaryfunc)func;
+    UNARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        PyObject **out = (PyObject **)op1;
+        PyObject *ret = f(in1);
+        if ((ret == NULL) || PyErr_Occurred()) {
+            return;
+        }
+        Py_XDECREF(*out);
+        *out = ret;
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_O_O_method(char **args, intp *dimensions, intp *steps, void *func)
+{
+    char *meth = (char *)func;
+    UNARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        PyObject **out = (PyObject **)op1;
+        PyObject *ret = PyObject_CallMethod(in1, meth, NULL);
+        if (ret == NULL) {
+            return;
+        }
+        Py_XDECREF(*out);
+        *out = ret;
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_OO_O(char **args, intp *dimensions, intp *steps, void *func)
+{
+    binaryfunc f = (binaryfunc)func;
+    BINARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        PyObject *in2 = *(PyObject **)ip2;
+        PyObject **out = (PyObject **)op1;
+        PyObject *ret = f(in1, in2);
+        if (PyErr_Occurred()) {
+            return;
+        }
+        Py_XDECREF(*out);
+        *out = ret;
+    }
+}
+
+/*UFUNC_API*/
+static void
+PyUFunc_OO_O_method(char **args, intp *dimensions, intp *steps, void *func)
+{
+    char *meth = (char *)func;
+    BINARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        PyObject *in2 = *(PyObject **)ip2;
+        PyObject **out = (PyObject **)op1;
+        PyObject *ret = PyObject_CallMethod(in1, meth, "(O)", in2);
+        if (ret == NULL) {
+            return;
+        }
+        Py_XDECREF(*out);
+        *out = ret;
+    }
+}
+
+/*
+ * A general-purpose ufunc that deals with general-purpose Python callable.
+ * func is a structure with nin, nout, and a Python callable function
+ */
+
+/*UFUNC_API*/
+static void
+PyUFunc_On_Om(char **args, intp *dimensions, intp *steps, void *func)
+{
+    intp n =  dimensions[0];
+    PyUFunc_PyFuncData *data = (PyUFunc_PyFuncData *)func;
+    int nin = data->nin;
+    int nout = data->nout;
+    PyObject *tocall = data->callable;
+    char *ptrs[NPY_MAXARGS];
+    PyObject *arglist, *result;
+    PyObject *in, **op;
+    intp i, j, ntot;
+
+    ntot = nin+nout;
+
+    for(j = 0; j < ntot; j++) {
+        ptrs[j] = args[j];
+    }
+    for(i = 0; i < n; i++) {
+        arglist = PyTuple_New(nin);
+        if (arglist == NULL) {
+            return;
+        }
+        for(j = 0; j < nin; j++) {
+            in = *((PyObject **)ptrs[j]);
+            if (in == NULL) {
+                Py_DECREF(arglist);
+                return;
+            }
+            PyTuple_SET_ITEM(arglist, j, in);
+            Py_INCREF(in);
+        }
+        result = PyEval_CallObject(tocall, arglist);
+        Py_DECREF(arglist);
+        if (result == NULL) {
+            return;
+        }
+        if PyTuple_Check(result) {
+            if (nout != PyTuple_Size(result)) {
+                Py_DECREF(result);
+                return;
+            }
+            for(j = 0; j < nout; j++) {
+                op = (PyObject **)ptrs[j+nin];
+                Py_XDECREF(*op);
+                *op = PyTuple_GET_ITEM(result, j);
+                Py_INCREF(*op);
+            }
+            Py_DECREF(result);
+        }
+        else {
+            op = (PyObject **)ptrs[nin];
+            Py_XDECREF(*op);
+            *op = result;
+        }
+        for(j = 0; j < ntot; j++) {
+            ptrs[j] += steps[j];
+        }
+    }
+}
+
+/*
+ *****************************************************************************
+ **                             BOOLEAN LOOPS                               **
+ *****************************************************************************
+ */
+
+#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
+#define BOOL_fmax BOOL_maximum
+#define BOOL_fmin BOOL_minimum
+
+/**begin repeat
+ * #kind = equal, not_equal, greater, greater_equal, less, less_equal,
+ *         logical_and, logical_or#
+ * #OP =  ==, !=, >, >=, <, <=, &&, ||#
+ **/
+
+static void
+BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        Bool in1 = *((Bool *)ip1) != 0;
+        Bool in2 = *((Bool *)ip2) != 0;
+        *((Bool *)op1)= in1 @OP@ in2;
+    }
+}
+/**end repeat**/
+
+static void
+BOOL_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        Bool in1 = *((Bool *)ip1) != 0;
+        Bool in2 = *((Bool *)ip2) != 0;
+        *((Bool *)op1)= (in1 && !in2) || (!in1 && in2);
+    }
+}
+
+/**begin repeat
+ * #kind = maximum, minimum#
+ * #OP =  >, <#
+ **/
+static void
+BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        Bool in1 = *((Bool *)ip1) != 0;
+        Bool in2 = *((Bool *)ip2) != 0;
+        *((Bool *)op1) = (in1 @OP@ in2) ? in1 : in2;
+    }
+}
+/**end repeat**/
+
+/**begin repeat
+ * #kind = absolute, logical_not#
+ * #OP =  !=, ==#
+ **/
+static void
+BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        Bool in1 = *(Bool *)ip1;
+        *((Bool *)op1) = in1 @OP@ 0;
+    }
+}
+/**end repeat**/
+
+static void
+BOOL_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    OUTPUT_LOOP {
+        *((Bool *)op1) = 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
+#define @S@@TYPE@_fmax @S@@TYPE@_maximum
+#define @S@@TYPE@_fmin @S@@TYPE@_minimum
+
+static void
+@S@@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    OUTPUT_LOOP {
+        *((@s@@type@ *)op1) = 1;
+    }
+}
+
+static void
+@S@@TYPE@_square(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op1) = in1*in1;
+    }
+}
+
+static void
+@S@@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op1) = (@s@@type@)(1.0/in1);
+    }
+}
+
+static void
+@S@@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op1) = in1;
+    }
+}
+
+static void
+@S@@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op1) = (@s@@type@)(-(@type@)in1);
+    }
+}
+
+static void
+@S@@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((Bool *)op1) = !in1;
+    }
+}
+
+static void
+@S@@TYPE@_invert(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        *((@s@@type@ *)op1) = ~in1;
+    }
+}
+
+/**begin repeat2
+ * Arithmetic
+ * #kind = add, subtract, multiply, bitwise_and, bitwise_or, bitwise_xor,
+ *          left_shift, right_shift#
+ * #OP = +, -,*, &, |, ^, <<, >>#
+ */
+static void
+@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((@s@@type@ *)op1) = in1 @OP@ in2;
+    }
+}
+/**end repeat2**/
+
+/**begin repeat2
+ * #kind = equal, not_equal, greater, greater_equal, less, less_equal,
+ *         logical_and, logical_or#
+ * #OP =  ==, !=, >, >=, <, <=, &&, ||#
+ */
+static void
+@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((Bool *)op1) = in1 @OP@ in2;
+    }
+}
+/**end repeat2**/
+
+static void
+@S@@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((Bool *)op1)= (in1 && !in2) || (!in1 && in2);
+    }
+}
+
+/**begin repeat2
+ * #kind = maximum, minimum#
+ * #OP =  >, <#
+ **/
+static void
+@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        *((@s@@type@ *)op1) = (in1 @OP@ in2) ? in1 : in2;
+    }
+}
+/**end repeat2**/
+
+static void
+@S@@TYPE@_true_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@ftype@ *)op1) = 0;
+        }
+        else {
+            *((@ftype@ *)op1) = (@ftype@)in1 / (@ftype@)in2;
+        }
+    }
+}
+
+static void
+@S@@TYPE@_power(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @ftype@ in1 = (@ftype@)*(@s@@type@ *)ip1;
+        const @ftype@ in2 = (@ftype@)*(@s@@type@ *)ip2;
+        *((@s@@type@ *)op1) = (@s@@type@) pow(in1, in2);
+    }
+}
+
+static void
+@S@@TYPE@_fmod(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @s@@type@ in1 = *(@s@@type@ *)ip1;
+        const @s@@type@ in2 = *(@s@@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@s@@type@ *)op1) = 0;
+        }
+        else {
+            *((@s@@type@ *)op1)= in1 % in2;
+        }
+
+    }
+}
+
+/**end repeat1**/
+
+static void
+U@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        *((u@type@ *)op1) = in1;
+    }
+}
+
+static void
+@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = (in1 >= 0) ? in1 : -in1;
+    }
+}
+
+static void
+U@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        *((u@type@ *)op1) = in1 > 0 ? 1 : 0;
+    }
+}
+
+static void
+@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0);
+    }
+}
+
+static void
+@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@type@ *)op1) = 0;
+        }
+        else if (((in1 > 0) != (in2 > 0)) && (in1 % in2 != 0)) {
+            *((@type@ *)op1) = in1/in2 - 1;
+        }
+        else {
+            *((@type@ *)op1) = in1/in2;
+        }
+    }
+}
+
+static void
+U@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        const u@type@ in2 = *(u@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((u@type@ *)op1) = 0;
+        }
+        else {
+            *((u@type@ *)op1)= in1/in2;
+        }
+    }
+}
+
+static void
+@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@type@ *)op1) = 0;
+        }
+        else {
+            /* handle mixed case the way Python does */
+            const @type@ rem = in1 % in2;
+            if ((in1 > 0) == (in2 > 0) || rem == 0) {
+                *((@type@ *)op1) = rem;
+            }
+            else {
+                *((@type@ *)op1) = rem + in2;
+            }
+        }
+    }
+}
+
+static void
+U@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const u@type@ in1 = *(u@type@ *)ip1;
+        const u@type@ in2 = *(u@type@ *)ip2;
+        if (in2 == 0) {
+            generate_divbyzero_error();
+            *((@type@ *)op1) = 0;
+        }
+        else {
+            *((@type@ *)op1) = in1 % in2;
+        }
+    }
+}
+
+/**end repeat**/
+
+/*
+ *****************************************************************************
+ **                             FLOAT LOOPS                                 **
+ *****************************************************************************
+ */
+
+
+/**begin repeat
+ * Float types
+ *  #type = float, double, longdouble#
+ *  #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
+ *  #c = f, , l#
+ *  #C = F, , L#
+ */
+
+
+/**begin repeat1
+ * Arithmetic
+ * # kind = add, subtract, multiply, divide#
+ * # OP = +, -, *, /#
+ */
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op1) = in1 @OP@ in2;
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = equal, not_equal, less, less_equal, greater, greater_equal,
+ *        logical_and, logical_or#
+ * #OP = ==, !=, <, <=, >, >=, &&, ||#
+ */
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((Bool *)op1) = in1 @OP@ in2;
+    }
+}
+/**end repeat1**/
+
+static void
+@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((Bool *)op1)= (in1 && !in2) || (!in1 && in2);
+    }
+}
+
+static void
+@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((Bool *)op1) = !in1;
+    }
+}
+
+/**begin repeat1
+ * #kind = isnan, isinf, isfinite, signbit#
+ * #func = isnan, isinf, isfinite, signbit#
+ **/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((Bool *)op1) = @func@(in1) != 0;
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = maximum, minimum#
+ * #OP =  >=, <=#
+ **/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    /*  */
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op1) = (in1 @OP@ in2 || isnan(in1)) ? in1 : in2;
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = fmax, fmin#
+ * #OP =  >=, <=#
+ **/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+    /*  */
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op1) = (in1 @OP@ in2 || isnan(in2)) ? in1 : in2;
+    }
+}
+/**end repeat1**/
+
+static void
+@TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op1) = floor@c@(in1/in2);
+    }
+}
+
+static void
+@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    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@ *)op1) = res + in2;
+        }
+        else {
+            *((@type@ *)op1) = res;
+        }
+    }
+}
+
+static void
+@TYPE@_square(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = in1*in1;
+    }
+}
+
+static void
+@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = 1/in1;
+    }
+}
+
+static void
+@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    OUTPUT_LOOP {
+        *((@type@ *)op1) = 1;
+    }
+}
+
+static void
+@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = in1;
+    }
+}
+
+static void
+@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ tmp = in1 > 0 ? in1 : -in1;
+        /* add 0 to clear -0.0 */
+        *((@type@ *)op1) = tmp + 0;
+    }
+}
+
+static void
+@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = -in1;
+    }
+}
+
+static void
+@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    /* Sign of nan is currently 0 */
+    UNARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0);
+    }
+}
+
+static void
+@TYPE@_modf(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP_TWO_OUT {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = modf@c@(in1, (@type@ *)op2);
+    }
+}
+
+#ifdef HAVE_FREXP@C@
+static void
+@TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP_TWO_OUT {
+        const @type@ in1 = *(@type@ *)ip1;
+        *((@type@ *)op1) = frexp@c@(in1, (int *)op2);
+    }
+}
+#endif
+
+#ifdef HAVE_LDEXP@C@
+static void
+@TYPE@_ldexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const int in2 = *(int *)ip2;
+        *((@type@ *)op1) = ldexp@c@(in1, in2);
+    }
+}
+#endif
+
+#define @TYPE@_true_divide @TYPE@_divide
+
+/**end repeat**/
+
+
+/*
+ *****************************************************************************
+ **                           COMPLEX LOOPS                                 **
+ *****************************************************************************
+ */
+
+#define CGE(xr,xi,yr,yi) (xr > yr || (xr == yr && xi >= yi))
+#define CLE(xr,xi,yr,yi) (xr < yr || (xr == yr && xi <= yi))
+#define CGT(xr,xi,yr,yi) (xr > yr || (xr == yr && xi > yi))
+#define CLT(xr,xi,yr,yi) (xr < yr || (xr == yr && xi < yi))
+#define CEQ(xr,xi,yr,yi) (xr == yr && xi == yi)
+#define CNE(xr,xi,yr,yi) (xr != yr || xi != yi)
+
+/**begin repeat
+ * complex types
+ * #type = float, double, longdouble#
+ * #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
+ * #c = f, , l#
+ */
+
+/**begin repeat1
+ * arithmetic
+ * #kind = add, subtract#
+ * #OP = +, -#
+ */
+static void
+C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = in1r @OP@ in2r;
+        ((@type@ *)op1)[1] = in1i @OP@ in2i;
+    }
+}
+/**end repeat1**/
+
+static void
+C@TYPE@_multiply(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = in1r*in2r - in1i*in2i;
+        ((@type@ *)op1)[1] = in1r*in2i + in1i*in2r;
+    }
+}
+
+static void
+C@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = (in1r*in2r + in1i*in2i)/d;
+        ((@type@ *)op1)[1] = (in1i*in2r - in1r*in2i)/d;
+    }
+}
+
+static void
+C@TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = floor@c@((in1r*in2r + in1i*in2i)/d);
+        ((@type@ *)op1)[1] = 0;
+    }
+}
+
+/**begin repeat1
+ * #kind= greater, greater_equal, less, less_equal, equal, not_equal#
+ * #OP = CGT, CGE, CLT, CLE, CEQ, CNE#
+ */
+static void
+C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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];
+        *((Bool *)op1) = @OP@(in1r,in1i,in2r,in2i);
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+   #kind = logical_and, logical_or#
+   #OP1 = ||, ||#
+   #OP2 = &&, ||#
+*/
+static void
+C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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];
+        *((Bool *)op1) = (in1r @OP1@ in1i) @OP2@ (in2r @OP1@ in2i);
+    }
+}
+/**end repeat1**/
+
+static void
+C@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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];
+        const Bool tmp1 = (in1r || in1i);
+        const Bool tmp2 = (in2r || in2i);
+        *((Bool *)op1) = (tmp1 && !tmp2) || (!tmp1 && tmp2);
+    }
+}
+
+static void
+C@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        *((Bool *)op1) = !(in1r || in1i);
+    }
+}
+
+/**begin repeat1
+ * #kind = isnan, isinf, isfinite#
+ * #func = isnan, isinf, isfinite#
+ * #OP = ||, ||, &&#
+ **/
+static void
+C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        *((Bool *)op1) = @func@(in1r) @OP@ @func@(in1i);
+    }
+}
+/**end repeat1**/
+
+static void
+C@TYPE@_square(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        ((@type@ *)op1)[0] = in1r*in1r - in1i*in1i;
+        ((@type@ *)op1)[1] = in1r*in1i + in1i*in1r;
+    }
+}
+
+static void
+C@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    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@ *)op1)[0] = 1/d;
+            ((@type@ *)op1)[1] = -r/d;
+        } else {
+            const @type@ r = in1r/in1i;
+            const @type@ d = in1r*r + in1i;
+            ((@type@ *)op1)[0] = r/d;
+            ((@type@ *)op1)[1] = -1/d;
+        }
+    }
+}
+
+static void
+C@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+    OUTPUT_LOOP {
+        ((@type@ *)op1)[0] = 1;
+        ((@type@ *)op1)[1] = 0;
+    }
+}
+
+static void
+C@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) {
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        ((@type@ *)op1)[0] = in1r;
+        ((@type@ *)op1)[1] = -in1i;
+    }
+}
+
+static void
+C@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        *((@type@ *)op1) = sqrt@c@(in1r*in1r + in1i*in1i);
+    }
+}
+
+static void
+C@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    /* fixme: sign of nan is currently 0 */
+    UNARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        ((@type@ *)op1)[0] = CGT(in1r, in1i, 0, 0) ? 1 :
+                            (CLT(in1r, in1i, 0, 0) ? -1 : 0);
+        ((@type@ *)op1)[1] = 0;
+    }
+}
+
+/**begin repeat1
+ * #kind = maximum, minimum#
+ * #OP = CGE, CLE#
+ */
+static void
+C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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 (@OP@(in1r, in1i, in2r, in2i) || isnan(in1r) || isnan(in1i)) {
+            ((@type@ *)op1)[0] = in1r;
+            ((@type@ *)op1)[1] = in1i;
+        }
+        else {
+            ((@type@ *)op1)[0] = in2r;
+            ((@type@ *)op1)[1] = in2i;
+        }
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = fmax, fmin#
+ * #OP = CGE, CLE#
+ */
+static void
+C@TYPE@_@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 (@OP@(in1r, in1i, in2r, in2i) || isnan(in2r) || isnan(in2i)) {
+            ((@type@ *)op1)[0] = in1r;
+            ((@type@ *)op1)[1] = in1i;
+        }
+        else {
+            ((@type@ *)op1)[0] = in2r;
+            ((@type@ *)op1)[1] = in2i;
+        }
+    }
+}
+/**end repeat1**/
+
+#define C@TYPE@_true_divide C@TYPE@_divide
+
+/**end repeat**/
+
+#undef CGE
+#undef CLE
+#undef CGT
+#undef CLT
+#undef CEQ
+#undef CNE
+
+/*
+ *****************************************************************************
+ **                            OBJECT LOOPS                                 **
+ *****************************************************************************
+ */
+
+/**begin repeat
+ * #kind = equal, not_equal, greater, greater_equal, less, less_equal#
+ * #OP = EQ, NE, GT, GE, LT, LE#
+ */
+static void
+OBJECT_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) {
+    BINARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        PyObject *in2 = *(PyObject **)ip2;
+        int ret = PyObject_RichCompareBool(in1, in2, Py_@OP@);
+        if (ret == -1) {
+            return;
+        }
+        *((Bool *)op1) = (Bool)ret;
+    }
+}
+/**end repeat**/
+
+static void
+OBJECT_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+    PyObject *zero = PyInt_FromLong(0);
+    UNARY_LOOP {
+        PyObject *in1 = *(PyObject **)ip1;
+        PyObject **out = (PyObject **)op1;
+        PyObject *ret = PyInt_FromLong(PyObject_Compare(in1, zero));
+        if (PyErr_Occurred()) {
+            return;
+        }
+        Py_XDECREF(*out);
+        *out = ret;
+    }
+    Py_DECREF(zero);
+}
+
+/*
+ *****************************************************************************
+ **                              END LOOPS                                  **
+ *****************************************************************************
+ */
+
+
diff --git a/numpy/core/src/umathmodule.c.src b/numpy/core/src/umathmodule.c.src
index d07bffa2a..8ac518d4a 100644
--- a/numpy/core/src/umathmodule.c.src
+++ b/numpy/core/src/umathmodule.c.src
@@ -9,1927 +9,25 @@
  **                            INCLUDES                                     **
  *****************************************************************************
  */
+#define _UMATHMODULE /* needed in one of the numpy include files */
+#include <math.h>
 #include "Python.h"
 #include "numpy/noprefix.h"
-#define _UMATHMODULE
 #include "numpy/ufuncobject.h"
 #include "abstract.h"
 #include "config.h"
-#include <math.h>
-
-#ifndef M_PI
-#define M_PI 3.14159265358979323846264338328
-#endif
-
-#include "umath_funcs_c99.inc"
-
-/*
- *****************************************************************************
- **                            FLOAT FUNCTIONS                              **
- *****************************************************************************
- */
-
-/**begin repeat
- * #type = float, double, longdouble#
- * #c = f, ,l#
- * #C = F, ,L#
- */
-
-/* fixme: need more precision for LOG2 and INVLOG2 */
-
-#define PI 3.14159265358979323846264338328@c@
-#define LOG2 0.69314718055994530943@c@
-#define INVLOG2 1.4426950408889634074@c@
-#define degrees@c@ rad2deg@c@
-#define radians@c@ deg2rad@c@
-
-static @type@
-rad2deg@c@(@type@ x) {
-    return x*(180.0@c@/PI);
-}
-
-static @type@
-deg2rad@c@(@type@ x) {
-    return x*(PI/180.0@c@);
-}
-
-static @type@
-log2_1p@c@(@type@ x)
-{
-    @type@ u = 1 + x;
-    if (u == 1) {
-        return INVLOG2*x;
-    } else {
-        return log2@c@(u) * x / (u - 1);
-    }
-}
-
-static @type@
-exp2_1m@c@(@type@ x)
-{
-    @type@ u = exp@c@(x);
-    if (u == 1.0) {
-        return LOG2*x;
-    } else if (u - 1 == -1) {
-        return -LOG2;
-    } else {
-        return (u - 1) * x/log2@c@(u);
-    }
-}
-
-static @type@
-logaddexp@c@(@type@ x, @type@ y)
-{
-    const @type@ tmp = x - y;
-    if (tmp > 0) {
-        return x + log1p@c@(exp@c@(-tmp));
-    }
-    else {
-        return y + log1p@c@(exp@c@(tmp));
-    }
-}
-
-static @type@
-logaddexp2@c@(@type@ x, @type@ y)
-{
-    const @type@ tmp = x - y;
-    if (tmp > 0) {
-        return x + log2_1p@c@(exp2@c@(-tmp));
-    }
-    else {
-        return y + log2_1p@c@(exp2@c@(tmp));
-    }
-}
-
-#undef PI
-#undef LOG2
-#undef INVLOG2
-
-/**end repeat**/
-
-/*
- *****************************************************************************
- **                        PYTHON OBJECT FUNCTIONS                          **
- *****************************************************************************
- */
-
-static PyObject *
-Py_square(PyObject *o)
-{
-    return PyNumber_Multiply(o, o);
-}
-
-static PyObject *
-Py_get_one(PyObject *NPY_UNUSED(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;
-}
-
-/*
- * Define numpy version of PyNumber_Power as binary function.
- */
-static PyObject *
-npy_ObjectPower(PyObject *x, PyObject *y)
-{
-    return PyNumber_Power(x, y, Py_None);
-}
-
-/**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**/
-
-
-/*
- *****************************************************************************
- **                           COMPLEX FUNCTIONS                             **
- *****************************************************************************
- */
-
-
-/*
- * Don't pass structures between functions (only pointers) because how
- * structures are passed is compiler dependent and could cause segfaults if
- * umath_ufunc_object.inc is compiled with a different compiler than an
- * extension that makes use of the UFUNC API
- */
-
-/**begin repeat
-
-   #typ=float, double, longdouble#
-   #c=f,,l#
-*/
-
-/* constants */
-static c@typ@ nc_1@c@ = {1., 0.};
-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 void
-nc_sum@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
-{
-    r->real = a->real + b->real;
-    r->imag = a->imag + b->imag;
-    return;
-}
-
-static void
-nc_diff@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
-{
-    r->real = a->real - b->real;
-    r->imag = a->imag - b->imag;
-    return;
-}
-
-static void
-nc_neg@c@(c@typ@ *a, c@typ@ *r)
-{
-    r->real = -a->real;
-    r->imag = -a->imag;
-    return;
-}
-
-static void
-nc_prod@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
-{
-    @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;
-}
-
-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;
-    r->imag = (ai*br - ar*bi)/d;
-    return;
-}
-
-static void
-nc_sqrt@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ s,d;
-    if (x->real == 0. && x->imag == 0.)
-        *r = *x;
-    else {
-        s = sqrt@c@((fabs@c@(x->real) + hypot@c@(x->real,x->imag))/2);
-        d = x->imag/(2*s);
-        if (x->real > 0) {
-            r->real = s;
-            r->imag = d;
-        }
-        else if (x->imag >= 0) {
-            r->real = d;
-            r->imag = s;
-        }
-        else {
-            r->real = -d;
-            r->imag = -s;
-        }
-    }
-    return;
-}
-
-static void
-nc_rint@c@(c@typ@ *x, c@typ@ *r)
-{
-    r->real = rint@c@(x->real);
-    r->imag = rint@c@(x->imag);
-}
-
-static void
-nc_log@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ l = hypot@c@(x->real,x->imag);
-    r->imag = atan2@c@(x->imag, x->real);
-    r->real = log@c@(l);
-    return;
-}
-
-static void
-nc_log1p@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ l = hypot@c@(x->real + 1,x->imag);
-    r->imag = atan2@c@(x->imag, x->real + 1);
-    r->real = log@c@(l);
-    return;
-}
-
-static void
-nc_exp@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ a = exp@c@(x->real);
-    r->real = a*cos@c@(x->imag);
-    r->imag = a*sin@c@(x->imag);
-    return;
-}
-
-static void
-nc_expm1@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ a = exp@c@(x->real);
-    r->real = a*cos@c@(x->imag) - 1;
-    r->imag = a*sin@c@(x->imag);
-    return;
-}
-
-static void
-nc_pow@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
-{
-    intp n;
-    @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
-
-    if (br == 0. && bi == 0.) {
-        r->real = 1.;
-        r->imag = 0.;
-        return;
-    }
-    if (ar == 0. && ai == 0.) {
-        r->real = 0.;
-        r->imag = 0.;
-        return;
-    }
-    if (bi == 0 && (n=(intp)br) == br) {
-        if (n > -100 && n < 100) {
-            c@typ@ p, aa;
-            intp mask = 1;
-            if (n < 0) n = -n;
-            aa = nc_1@c@;
-            p.real = ar; p.imag = ai;
-            while (1) {
-                if (n & mask)
-                    nc_prod@c@(&aa,&p,&aa);
-                mask <<= 1;
-                if (n < mask || mask <= 0) break;
-                nc_prod@c@(&p,&p,&p);
-            }
-            r->real = aa.real; r->imag = aa.imag;
-            if (br < 0) nc_quot@c@(&nc_1@c@, r, r);
-            return;
-        }
-    }
-    /*
-     * 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);
-    return;
-}
-
-
-static void
-nc_prodi@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ xr = x->real;
-    r->real = -x->imag;
-    r->imag = xr;
-    return;
-}
-
-
-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);
-    nc_prodi@c@(r, r);
-    nc_sum@c@(x, r, r);
-    nc_log@c@(r, r);
-    nc_prodi@c@(r, r);
-    nc_neg@c@(r, r);
-    return;
-}
-
-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);
-    nc_sqrt@c@(&t, &t);
-    nc_diff@c@(x, &nc_1@c@, r);
-    nc_sqrt@c@(r, r);
-    nc_prod@c@(&t, r, r);
-    nc_sum@c@(x, r, r);
-    nc_log@c@(r, r);
-    return;
-}
-
-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);
-    nc_sqrt@c@(r, r);
-    nc_prodi@c@(x, pa);
-    nc_sum@c@(pa, r, r);
-    nc_log@c@(r, r);
-    nc_prodi@c@(r, r);
-    nc_neg@c@(r, r);
-    return;
-}
-
-
-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;
-}
-
-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);
-    nc_quot@c@(r, pa, r);
-    nc_log@c@(r,r);
-    nc_prod@c@(&nc_i2@c@, r, r);
-    return;
-}
-
-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);
-    nc_quot@c@(pa, r, r);
-    nc_log@c@(r, r);
-    nc_prod@c@(&nc_half@c@, r, r);
-    return;
-}
-
-static void
-nc_cos@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ xr=x->real, xi=x->imag;
-    r->real = cos@c@(xr)*cosh@c@(xi);
-    r->imag = -sin@c@(xr)*sinh@c@(xi);
-    return;
-}
-
-static void
-nc_cosh@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ xr=x->real, xi=x->imag;
-    r->real = cos@c@(xi)*cosh@c@(xr);
-    r->imag = sin@c@(xi)*sinh@c@(xr);
-    return;
-}
-
-
-#define M_LOG10_E 0.434294481903251827651128918916605082294397
-
-static void
-nc_log10@c@(c@typ@ *x, c@typ@ *r)
-{
-    nc_log@c@(x, r);
-    r->real *= (@typ@) M_LOG10_E;
-    r->imag *= (@typ@) M_LOG10_E;
-    return;
-}
-
-static void
-nc_sin@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ xr=x->real, xi=x->imag;
-    r->real = sin@c@(xr)*cosh@c@(xi);
-    r->imag = cos@c@(xr)*sinh@c@(xi);
-    return;
-}
-
-static void
-nc_sinh@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ xr=x->real, xi=x->imag;
-    r->real = cos@c@(xi)*sinh@c@(xr);
-    r->imag = sin@c@(xi)*cosh@c@(xr);
-    return;
-}
-
-static void
-nc_tan@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ sr,cr,shi,chi;
-    @typ@ rs,is,rc,ic;
-    @typ@ d;
-    @typ@ xr=x->real, xi=x->imag;
-    sr = sin@c@(xr);
-    cr = cos@c@(xr);
-    shi = sinh@c@(xi);
-    chi = cosh@c@(xi);
-    rs = sr*chi;
-    is = cr*shi;
-    rc = cr*chi;
-    ic = -sr*shi;
-    d = rc*rc + ic*ic;
-    r->real = (rs*rc+is*ic)/d;
-    r->imag = (is*rc-rs*ic)/d;
-    return;
-}
-
-static void
-nc_tanh@c@(c@typ@ *x, c@typ@ *r)
-{
-    @typ@ si,ci,shr,chr;
-    @typ@ rs,is,rc,ic;
-    @typ@ d;
-    @typ@ xr=x->real, xi=x->imag;
-    si = sin@c@(xi);
-    ci = cos@c@(xi);
-    shr = sinh@c@(xr);
-    chr = cosh@c@(xr);
-    rs = ci*shr;
-    is = si*chr;
-    rc = ci*chr;
-    ic = si*shr;
-    d = rc*rc + ic*ic;
-    r->real = (rs*rc+is*ic)/d;
-    r->imag = (is*rc-rs*ic)/d;
-    return;
-}
-
-/**end repeat**/
-
-/*
- *****************************************************************************
- **                             UFUNC LOOPS                                 **
- *****************************************************************************
- */
-
-#define OUTPUT_LOOP\
-    char *op1 = args[1];\
-    intp os1 = steps[1];\
-    intp n = dimensions[0];\
-    intp i;\
-    for(i = 0; i < n; i++, op1 += os1)
-
-#define UNARY_LOOP\
-    char *ip1 = args[0], *op1 = args[1];\
-    intp is1 = steps[0], os1 = steps[1];\
-    intp n = dimensions[0];\
-    intp i;\
-    for(i = 0; i < n; i++, ip1 += is1, op1 += os1)
-
-#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], *op1 = args[2];\
-    intp is1 = steps[0], is2 = steps[1], os1 = steps[2];\
-    intp n = dimensions[0];\
-    intp i;\
-    for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1)
-
-#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)
-
-/******************************************************************************
- **                          GENERIC FLOAT LOOPS                             **
- *****************************************************************************/
-
-
-typedef float floatUnaryFunc(float x);
-typedef double doubleUnaryFunc(double x);
-typedef longdouble longdoubleUnaryFunc(longdouble x);
-typedef float floatBinaryFunc(float x, float y);
-typedef double doubleBinaryFunc(double x, double y);
-typedef longdouble longdoubleBinaryFunc(longdouble x, longdouble y);
-
-
-/*UFUNC_API*/
-static void
-PyUFunc_f_f(char **args, intp *dimensions, intp *steps, void *func)
-{
-    floatUnaryFunc *f = (floatUnaryFunc *)func;
-    UNARY_LOOP {
-        const float in1 = *(float *)ip1;
-        *(float *)op1 = f(in1);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_f_f_As_d_d(char **args, intp *dimensions, intp *steps, void *func)
-{
-    doubleUnaryFunc *f = (doubleUnaryFunc *)func;
-    UNARY_LOOP {
-        const float in1 = *(float *)ip1;
-        *(float *)op1 = (float)f((double)in1);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_ff_f(char **args, intp *dimensions, intp *steps, void *func)
-{
-    floatBinaryFunc *f = (floatBinaryFunc *)func;
-    BINARY_LOOP {
-        float in1 = *(float *)ip1;
-        float in2 = *(float *)ip2;
-        *(float *)op1 = f(in1, in2);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_ff_f_As_dd_d(char **args, intp *dimensions, intp *steps, void *func)
-{
-    doubleBinaryFunc *f = (doubleBinaryFunc *)func;
-    BINARY_LOOP {
-        float in1 = *(float *)ip1;
-        float in2 = *(float *)ip2;
-        *(float *)op1 = (double)f((double)in1, (double)in2);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_d_d(char **args, intp *dimensions, intp *steps, void *func)
-{
-    doubleUnaryFunc *f = (doubleUnaryFunc *)func;
-    UNARY_LOOP {
-        double in1 = *(double *)ip1;
-        *(double *)op1 = f(in1);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_dd_d(char **args, intp *dimensions, intp *steps, void *func)
-{
-    doubleBinaryFunc *f = (doubleBinaryFunc *)func;
-    BINARY_LOOP {
-        double in1 = *(double *)ip1;
-        double in2 = *(double *)ip2;
-        *(double *)op1 = f(in1, in2);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_g_g(char **args, intp *dimensions, intp *steps, void *func)
-{
-    longdoubleUnaryFunc *f = (longdoubleUnaryFunc *)func;
-    UNARY_LOOP {
-        longdouble in1 = *(longdouble *)ip1;
-        *(longdouble *)op1 = f(in1);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_gg_g(char **args, intp *dimensions, intp *steps, void *func)
-{
-    longdoubleBinaryFunc *f = (longdoubleBinaryFunc *)func;
-    BINARY_LOOP {
-        longdouble in1 = *(longdouble *)ip1;
-        longdouble in2 = *(longdouble *)ip2;
-        *(longdouble *)op1 = f(in1, in2);
-    }
-}
-
-
-
-/******************************************************************************
- **                          GENERIC COMPLEX LOOPS                           **
- *****************************************************************************/
-
-
-typedef void cdoubleUnaryFunc(cdouble *x, cdouble *r);
-typedef void cfloatUnaryFunc(cfloat *x, cfloat *r);
-typedef void clongdoubleUnaryFunc(clongdouble *x, clongdouble *r);
-typedef void cdoubleBinaryFunc(cdouble *x, cdouble *y, cdouble *r);
-typedef void cfloatBinaryFunc(cfloat *x, cfloat *y, cfloat *r);
-typedef void clongdoubleBinaryFunc(clongdouble *x, clongdouble *y,
-                                   clongdouble *r);
-
-/*UFUNC_API*/
-static void
-PyUFunc_F_F(char **args, intp *dimensions, intp *steps, void *func)
-{
-    cfloatUnaryFunc *f = (cfloatUnaryFunc *)func;
-    UNARY_LOOP {
-        cfloat in1 = *(cfloat *)ip1;
-        cfloat *out = (cfloat *)op1;
-        f(&in1, out);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_F_F_As_D_D(char **args, intp *dimensions, intp *steps, void *func)
-{
-    cdoubleUnaryFunc *f = (cdoubleUnaryFunc *)func;
-    UNARY_LOOP {
-        const float *in1 = (float *)ip1;
-        cdouble tmp = {(double)(in1[0]),(double)in1[1]};
-        cdouble out;
-        f(&tmp, &out);
-        ((float *)op1)[0] = (float)out.real;
-        ((float *)op1)[1] = (float)out.imag;
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_FF_F(char **args, intp *dimensions, intp *steps, void *func)
-{
-    cfloatBinaryFunc *f = (cfloatBinaryFunc *)func;
-    BINARY_LOOP {
-        cfloat in1 = *(cfloat *)ip1;
-        cfloat in2 = *(cfloat *)ip2;
-        cfloat *out = (cfloat *)op1;
-        f(&in1, &in2, out);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_FF_F_As_DD_D(char **args, intp *dimensions, intp *steps, void *func)
-{
-    cdoubleBinaryFunc *f = (cdoubleBinaryFunc *)func;
-    BINARY_LOOP {
-        const float *in1 = (float *)ip1;
-        const float *in2 = (float *)ip2;
-        cdouble tmp1 = {(double)(in1[0]),(double)in1[1]};
-        cdouble tmp2 = {(double)(in2[0]),(double)in2[1]};
-        cdouble out;
-        f(&tmp1, &tmp2, &out);
-        ((float *)op1)[0] = (float)out.real;
-        ((float *)op1)[1] = (float)out.imag;
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_D_D(char **args, intp *dimensions, intp *steps, void *func)
-{
-    cdoubleUnaryFunc *f = (cdoubleUnaryFunc *)func;
-    UNARY_LOOP {
-        cdouble in1 = *(cdouble *)ip1;
-        cdouble *out = (cdouble *)op1;
-        f(&in1, out);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_DD_D(char **args, intp *dimensions, intp *steps, void *func)
-{
-    cdoubleBinaryFunc *f = (cdoubleBinaryFunc *)func;
-    BINARY_LOOP {
-        cdouble in1 = *(cdouble *)ip1;
-        cdouble in2 = *(cdouble *)ip2;
-        cdouble *out = (cdouble *)op1;
-        f(&in1, &in2, out);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_G_G(char **args, intp *dimensions, intp *steps, void *func)
-{
-    clongdoubleUnaryFunc *f = (clongdoubleUnaryFunc *)func;
-    UNARY_LOOP {
-        clongdouble in1 = *(clongdouble *)ip1;
-        clongdouble *out = (clongdouble *)op1;
-        f(&in1, out);
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_GG_G(char **args, intp *dimensions, intp *steps, void *func)
-{
-    clongdoubleBinaryFunc *f = (clongdoubleBinaryFunc *)func;
-    BINARY_LOOP {
-        clongdouble in1 = *(clongdouble *)ip1;
-        clongdouble in2 = *(clongdouble *)ip2;
-        clongdouble *out = (clongdouble *)op1;
-        f(&in1, &in2, out);
-    }
-}
-
-
-/******************************************************************************
- **                         GENERIC OBJECT lOOPS                             **
- *****************************************************************************/
-
-/*UFUNC_API*/
-static void
-PyUFunc_O_O(char **args, intp *dimensions, intp *steps, void *func)
-{
-    unaryfunc f = (unaryfunc)func;
-    UNARY_LOOP {
-        PyObject *in1 = *(PyObject **)ip1;
-        PyObject **out = (PyObject **)op1;
-        PyObject *ret = f(in1);
-        if ((ret == NULL) || PyErr_Occurred()) {
-            return;
-        }
-        Py_XDECREF(*out);
-        *out = ret;
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_O_O_method(char **args, intp *dimensions, intp *steps, void *func)
-{
-    char *meth = (char *)func;
-    UNARY_LOOP {
-        PyObject *in1 = *(PyObject **)ip1;
-        PyObject **out = (PyObject **)op1;
-        PyObject *ret = PyObject_CallMethod(in1, meth, NULL);
-        if (ret == NULL) {
-            return;
-        }
-        Py_XDECREF(*out);
-        *out = ret;
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_OO_O(char **args, intp *dimensions, intp *steps, void *func)
-{
-    binaryfunc f = (binaryfunc)func;
-    BINARY_LOOP {
-        PyObject *in1 = *(PyObject **)ip1;
-        PyObject *in2 = *(PyObject **)ip2;
-        PyObject **out = (PyObject **)op1;
-        PyObject *ret = f(in1, in2);
-        if (PyErr_Occurred()) {
-            return;
-        }
-        Py_XDECREF(*out);
-        *out = ret;
-    }
-}
-
-/*UFUNC_API*/
-static void
-PyUFunc_OO_O_method(char **args, intp *dimensions, intp *steps, void *func)
-{
-    char *meth = (char *)func;
-    BINARY_LOOP {
-        PyObject *in1 = *(PyObject **)ip1;
-        PyObject *in2 = *(PyObject **)ip2;
-        PyObject **out = (PyObject **)op1;
-        PyObject *ret = PyObject_CallMethod(in1, meth, "(O)", in2);
-        if (ret == NULL) {
-            return;
-        }
-        Py_XDECREF(*out);
-        *out = ret;
-    }
-}
-
-/*
- * A general-purpose ufunc that deals with general-purpose Python callable.
- * func is a structure with nin, nout, and a Python callable function
- */
-
-/*UFUNC_API*/
-static void
-PyUFunc_On_Om(char **args, intp *dimensions, intp *steps, void *func)
-{
-    intp n =  dimensions[0];
-    PyUFunc_PyFuncData *data = (PyUFunc_PyFuncData *)func;
-    int nin = data->nin;
-    int nout = data->nout;
-    PyObject *tocall = data->callable;
-    char *ptrs[NPY_MAXARGS];
-    PyObject *arglist, *result;
-    PyObject *in, **op;
-    intp i, j, ntot;
-
-    ntot = nin+nout;
-
-    for(j = 0; j < ntot; j++) {
-        ptrs[j] = args[j];
-    }
-    for(i = 0; i < n; i++) {
-        arglist = PyTuple_New(nin);
-        if (arglist == NULL) {
-            return;
-        }
-        for(j = 0; j < nin; j++) {
-            in = *((PyObject **)ptrs[j]);
-            if (in == NULL) {
-                Py_DECREF(arglist);
-                return;
-            }
-            PyTuple_SET_ITEM(arglist, j, in);
-            Py_INCREF(in);
-        }
-        result = PyEval_CallObject(tocall, arglist);
-        Py_DECREF(arglist);
-        if (result == NULL) {
-            return;
-        }
-        if PyTuple_Check(result) {
-            if (nout != PyTuple_Size(result)) {
-                Py_DECREF(result);
-                return;
-            }
-            for(j = 0; j < nout; j++) {
-                op = (PyObject **)ptrs[j+nin];
-                Py_XDECREF(*op);
-                *op = PyTuple_GET_ITEM(result, j);
-                Py_INCREF(*op);
-            }
-            Py_DECREF(result);
-        }
-        else {
-            op = (PyObject **)ptrs[nin];
-            Py_XDECREF(*op);
-            *op = result;
-        }
-        for(j = 0; j < ntot; j++) {
-            ptrs[j] += steps[j];
-        }
-    }
-}
-
-/*
- *****************************************************************************
- **                             BOOLEAN LOOPS                               **
- *****************************************************************************
- */
-
-#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
-#define BOOL_fmax BOOL_maximum
-#define BOOL_fmin BOOL_minimum
-
-/**begin repeat
- * #kind = equal, not_equal, greater, greater_equal, less, less_equal,
- *         logical_and, logical_or#
- * #OP =  ==, !=, >, >=, <, <=, &&, ||#
- **/
-
-static void
-BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        Bool in1 = *((Bool *)ip1) != 0;
-        Bool in2 = *((Bool *)ip2) != 0;
-        *((Bool *)op1)= in1 @OP@ in2;
-    }
-}
-/**end repeat**/
-
-static void
-BOOL_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        Bool in1 = *((Bool *)ip1) != 0;
-        Bool in2 = *((Bool *)ip2) != 0;
-        *((Bool *)op1)= (in1 && !in2) || (!in1 && in2);
-    }
-}
-
-/**begin repeat
- * #kind = maximum, minimum#
- * #OP =  >, <#
- **/
-static void
-BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        Bool in1 = *((Bool *)ip1) != 0;
-        Bool in2 = *((Bool *)ip2) != 0;
-        *((Bool *)op1) = (in1 @OP@ in2) ? in1 : in2;
-    }
-}
-/**end repeat**/
-
-/**begin repeat
- * #kind = absolute, logical_not#
- * #OP =  !=, ==#
- **/
-static void
-BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        Bool in1 = *(Bool *)ip1;
-        *((Bool *)op1) = in1 @OP@ 0;
-    }
-}
-/**end repeat**/
-
-static void
-BOOL_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    OUTPUT_LOOP {
-        *((Bool *)op1) = 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
-#define @S@@TYPE@_fmax @S@@TYPE@_maximum
-#define @S@@TYPE@_fmin @S@@TYPE@_minimum
-
-static void
-@S@@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    OUTPUT_LOOP {
-        *((@s@@type@ *)op1) = 1;
-    }
-}
-
-static void
-@S@@TYPE@_square(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    UNARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        *((@s@@type@ *)op1) = in1*in1;
-    }
-}
-
-static void
-@S@@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    UNARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        *((@s@@type@ *)op1) = (@s@@type@)(1.0/in1);
-    }
-}
-
-static void
-@S@@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        *((@s@@type@ *)op1) = in1;
-    }
-}
-
-static void
-@S@@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        *((@s@@type@ *)op1) = (@s@@type@)(-(@type@)in1);
-    }
-}
-
-static void
-@S@@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        *((Bool *)op1) = !in1;
-    }
-}
-
-static void
-@S@@TYPE@_invert(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        *((@s@@type@ *)op1) = ~in1;
-    }
-}
-
-/**begin repeat2
- * Arithmetic
- * #kind = add, subtract, multiply, bitwise_and, bitwise_or, bitwise_xor,
- *          left_shift, right_shift#
- * #OP = +, -,*, &, |, ^, <<, >>#
- */
-static void
-@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        const @s@@type@ in2 = *(@s@@type@ *)ip2;
-        *((@s@@type@ *)op1) = in1 @OP@ in2;
-    }
-}
-/**end repeat2**/
-
-/**begin repeat2
- * #kind = equal, not_equal, greater, greater_equal, less, less_equal,
- *         logical_and, logical_or#
- * #OP =  ==, !=, >, >=, <, <=, &&, ||#
- */
-static void
-@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        const @s@@type@ in2 = *(@s@@type@ *)ip2;
-        *((Bool *)op1) = in1 @OP@ in2;
-    }
-}
-/**end repeat2**/
-
-static void
-@S@@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        const @s@@type@ in2 = *(@s@@type@ *)ip2;
-        *((Bool *)op1)= (in1 && !in2) || (!in1 && in2);
-    }
-}
-
-/**begin repeat2
- * #kind = maximum, minimum#
- * #OP =  >, <#
- **/
-static void
-@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        const @s@@type@ in2 = *(@s@@type@ *)ip2;
-        *((@s@@type@ *)op1) = (in1 @OP@ in2) ? in1 : in2;
-    }
-}
-/**end repeat2**/
-
-static void
-@S@@TYPE@_true_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        const @s@@type@ in2 = *(@s@@type@ *)ip2;
-        if (in2 == 0) {
-            generate_divbyzero_error();
-            *((@ftype@ *)op1) = 0;
-        }
-        else {
-            *((@ftype@ *)op1) = (@ftype@)in1 / (@ftype@)in2;
-        }
-    }
-}
-
-static void
-@S@@TYPE@_power(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @ftype@ in1 = (@ftype@)*(@s@@type@ *)ip1;
-        const @ftype@ in2 = (@ftype@)*(@s@@type@ *)ip2;
-        *((@s@@type@ *)op1) = (@s@@type@) pow(in1, in2);
-    }
-}
-
-static void
-@S@@TYPE@_fmod(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        const @s@@type@ in2 = *(@s@@type@ *)ip2;
-        if (in2 == 0) {
-            generate_divbyzero_error();
-            *((@s@@type@ *)op1) = 0;
-        }
-        else {
-            *((@s@@type@ *)op1)= in1 % in2;
-        }
-
-    }
-}
-
-/**end repeat1**/
-
-static void
-U@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const u@type@ in1 = *(u@type@ *)ip1;
-        *((u@type@ *)op1) = in1;
-    }
-}
-
-static void
-@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = (in1 >= 0) ? in1 : -in1;
-    }
-}
-
-static void
-U@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const u@type@ in1 = *(u@type@ *)ip1;
-        *((u@type@ *)op1) = in1 > 0 ? 1 : 0;
-    }
-}
-
-static void
-@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0);
-    }
-}
-
-static void
-@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        if (in2 == 0) {
-            generate_divbyzero_error();
-            *((@type@ *)op1) = 0;
-        }
-        else if (((in1 > 0) != (in2 > 0)) && (in1 % in2 != 0)) {
-            *((@type@ *)op1) = in1/in2 - 1;
-        }
-        else {
-            *((@type@ *)op1) = in1/in2;
-        }
-    }
-}
-
-static void
-U@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const u@type@ in1 = *(u@type@ *)ip1;
-        const u@type@ in2 = *(u@type@ *)ip2;
-        if (in2 == 0) {
-            generate_divbyzero_error();
-            *((u@type@ *)op1) = 0;
-        }
-        else {
-            *((u@type@ *)op1)= in1/in2;
-        }
-    }
-}
-
-static void
-@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        if (in2 == 0) {
-            generate_divbyzero_error();
-            *((@type@ *)op1) = 0;
-        }
-        else {
-            /* handle mixed case the way Python does */
-            const @type@ rem = in1 % in2;
-            if ((in1 > 0) == (in2 > 0) || rem == 0) {
-                *((@type@ *)op1) = rem;
-            }
-            else {
-                *((@type@ *)op1) = rem + in2;
-            }
-        }
-    }
-}
-
-static void
-U@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const u@type@ in1 = *(u@type@ *)ip1;
-        const u@type@ in2 = *(u@type@ *)ip2;
-        if (in2 == 0) {
-            generate_divbyzero_error();
-            *((@type@ *)op1) = 0;
-        }
-        else {
-            *((@type@ *)op1) = in1 % in2;
-        }
-    }
-}
-
-/**end repeat**/
-
-/*
- *****************************************************************************
- **                             FLOAT LOOPS                                 **
- *****************************************************************************
- */
-
-
-/**begin repeat
- * Float types
- *  #type = float, double, longdouble#
- *  #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
- *  #c = f, , l#
- *  #C = F, , L#
- */
-
-
-/**begin repeat1
- * Arithmetic
- * # kind = add, subtract, multiply, divide#
- * # OP = +, -, *, /#
- */
-static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        *((@type@ *)op1) = in1 @OP@ in2;
-    }
-}
-/**end repeat1**/
-
-/**begin repeat1
- * #kind = equal, not_equal, less, less_equal, greater, greater_equal,
- *        logical_and, logical_or#
- * #OP = ==, !=, <, <=, >, >=, &&, ||#
- */
-static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        *((Bool *)op1) = in1 @OP@ in2;
-    }
-}
-/**end repeat1**/
-
-static void
-@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        *((Bool *)op1)= (in1 && !in2) || (!in1 && in2);
-    }
-}
-
-static void
-@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((Bool *)op1) = !in1;
-    }
-}
-
-/**begin repeat1
- * #kind = isnan, isinf, isfinite, signbit#
- * #func = isnan, isinf, isfinite, signbit#
- **/
-static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((Bool *)op1) = @func@(in1) != 0;
-    }
-}
-/**end repeat1**/
-
-/**begin repeat1
- * #kind = maximum, minimum#
- * #OP =  >=, <=#
- **/
-static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    /*  */
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        *((@type@ *)op1) = (in1 @OP@ in2 || isnan(in1)) ? in1 : in2;
-    }
-}
-/**end repeat1**/
-
-/**begin repeat1
- * #kind = fmax, fmin#
- * #OP =  >=, <=#
- **/
-static void
-@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
-{
-    /*  */
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        *((@type@ *)op1) = (in1 @OP@ in2 || isnan(in2)) ? in1 : in2;
-    }
-}
-/**end repeat1**/
-
-static void
-@TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ in2 = *(@type@ *)ip2;
-        *((@type@ *)op1) = floor@c@(in1/in2);
-    }
-}
-
-static void
-@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    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@ *)op1) = res + in2;
-        }
-        else {
-            *((@type@ *)op1) = res;
-        }
-    }
-}
-
-static void
-@TYPE@_square(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = in1*in1;
-    }
-}
-
-static void
-@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = 1/in1;
-    }
-}
-
-static void
-@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    OUTPUT_LOOP {
-        *((@type@ *)op1) = 1;
-    }
-}
-
-static void
-@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = in1;
-    }
-}
-
-static void
-@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const @type@ tmp = in1 > 0 ? in1 : -in1;
-        /* add 0 to clear -0.0 */
-        *((@type@ *)op1) = tmp + 0;
-    }
-}
-
-static void
-@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = -in1;
-    }
-}
-
-static void
-@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    /* Sign of nan is currently 0 */
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0);
-    }
-}
-
-static void
-@TYPE@_modf(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP_TWO_OUT {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = modf@c@(in1, (@type@ *)op2);
-    }
-}
-
-#ifdef HAVE_FREXP@C@
-static void
-@TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP_TWO_OUT {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = frexp@c@(in1, (int *)op2);
-    }
-}
-#endif
-
-#ifdef HAVE_LDEXP@C@
-static void
-@TYPE@_ldexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    BINARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        const int in2 = *(int *)ip2;
-        *((@type@ *)op1) = ldexp@c@(in1, in2);
-    }
-}
-#endif
-
-#define @TYPE@_true_divide @TYPE@_divide
-
-/**end repeat**/
-
-
-/*
- *****************************************************************************
- **                           COMPLEX LOOPS                                 **
- *****************************************************************************
- */
-
-#define CGE(xr,xi,yr,yi) (xr > yr || (xr == yr && xi >= yi))
-#define CLE(xr,xi,yr,yi) (xr < yr || (xr == yr && xi <= yi))
-#define CGT(xr,xi,yr,yi) (xr > yr || (xr == yr && xi > yi))
-#define CLT(xr,xi,yr,yi) (xr < yr || (xr == yr && xi < yi))
-#define CEQ(xr,xi,yr,yi) (xr == yr && xi == yi)
-#define CNE(xr,xi,yr,yi) (xr != yr || xi != yi)
-
-/**begin repeat
- * complex types
- * #type = float, double, longdouble#
- * #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
- * #c = f, , l#
- */
-
-/**begin repeat1
- * arithmetic
- * #kind = add, subtract#
- * #OP = +, -#
- */
-static void
-C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = in1r @OP@ in2r;
-        ((@type@ *)op1)[1] = in1i @OP@ in2i;
-    }
-}
-/**end repeat1**/
-
-static void
-C@TYPE@_multiply(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = in1r*in2r - in1i*in2i;
-        ((@type@ *)op1)[1] = in1r*in2i + in1i*in2r;
-    }
-}
-
-static void
-C@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = (in1r*in2r + in1i*in2i)/d;
-        ((@type@ *)op1)[1] = (in1i*in2r - in1r*in2i)/d;
-    }
-}
-
-static void
-C@TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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@ *)op1)[0] = floor@c@((in1r*in2r + in1i*in2i)/d);
-        ((@type@ *)op1)[1] = 0;
-    }
-}
-
-/**begin repeat1
- * #kind= greater, greater_equal, less, less_equal, equal, not_equal#
- * #OP = CGT, CGE, CLT, CLE, CEQ, CNE#
- */
-static void
-C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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];
-        *((Bool *)op1) = @OP@(in1r,in1i,in2r,in2i);
-    }
-}
-/**end repeat1**/
-
-/**begin repeat1
-   #kind = logical_and, logical_or#
-   #OP1 = ||, ||#
-   #OP2 = &&, ||#
-*/
-static void
-C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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];
-        *((Bool *)op1) = (in1r @OP1@ in1i) @OP2@ (in2r @OP1@ in2i);
-    }
-}
-/**end repeat1**/
-
-static void
-C@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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];
-        const Bool tmp1 = (in1r || in1i);
-        const Bool tmp2 = (in2r || in2i);
-        *((Bool *)op1) = (tmp1 && !tmp2) || (!tmp1 && tmp2);
-    }
-}
-
-static void
-C@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1r = ((@type@ *)ip1)[0];
-        const @type@ in1i = ((@type@ *)ip1)[1];
-        *((Bool *)op1) = !(in1r || in1i);
-    }
-}
-
-/**begin repeat1
- * #kind = isnan, isinf, isfinite#
- * #func = isnan, isinf, isfinite#
- * #OP = ||, ||, &&#
- **/
-static void
-C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1r = ((@type@ *)ip1)[0];
-        const @type@ in1i = ((@type@ *)ip1)[1];
-        *((Bool *)op1) = @func@(in1r) @OP@ @func@(in1i);
-    }
-}
-/**end repeat1**/
-
-static void
-C@TYPE@_square(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    UNARY_LOOP {
-        const @type@ in1r = ((@type@ *)ip1)[0];
-        const @type@ in1i = ((@type@ *)ip1)[1];
-        ((@type@ *)op1)[0] = in1r*in1r - in1i*in1i;
-        ((@type@ *)op1)[1] = in1r*in1i + in1i*in1r;
-    }
-}
-
-static void
-C@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    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@ *)op1)[0] = 1/d;
-            ((@type@ *)op1)[1] = -r/d;
-        } else {
-            const @type@ r = in1r/in1i;
-            const @type@ d = in1r*r + in1i;
-            ((@type@ *)op1)[0] = r/d;
-            ((@type@ *)op1)[1] = -1/d;
-        }
-    }
-}
-
-static void
-C@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
-    OUTPUT_LOOP {
-        ((@type@ *)op1)[0] = 1;
-        ((@type@ *)op1)[1] = 0;
-    }
-}
-
-static void
-C@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) {
-    UNARY_LOOP {
-        const @type@ in1r = ((@type@ *)ip1)[0];
-        const @type@ in1i = ((@type@ *)ip1)[1];
-        ((@type@ *)op1)[0] = in1r;
-        ((@type@ *)op1)[1] = -in1i;
-    }
-}
-
-static void
-C@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    UNARY_LOOP {
-        const @type@ in1r = ((@type@ *)ip1)[0];
-        const @type@ in1i = ((@type@ *)ip1)[1];
-        *((@type@ *)op1) = sqrt@c@(in1r*in1r + in1i*in1i);
-    }
-}
-
-static void
-C@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    /* fixme: sign of nan is currently 0 */
-    UNARY_LOOP {
-        const @type@ in1r = ((@type@ *)ip1)[0];
-        const @type@ in1i = ((@type@ *)ip1)[1];
-        ((@type@ *)op1)[0] = CGT(in1r, in1i, 0, 0) ? 1 :
-                            (CLT(in1r, in1i, 0, 0) ? -1 : 0);
-        ((@type@ *)op1)[1] = 0;
-    }
-}
-
-/**begin repeat1
- * #kind = maximum, minimum#
- * #OP = CGE, CLE#
- */
-static void
-C@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(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 (@OP@(in1r, in1i, in2r, in2i) || isnan(in1r) || isnan(in1i)) {
-            ((@type@ *)op1)[0] = in1r;
-            ((@type@ *)op1)[1] = in1i;
-        }
-        else {
-            ((@type@ *)op1)[0] = in2r;
-            ((@type@ *)op1)[1] = in2i;
-        }
-    }
-}
-/**end repeat1**/
-
-/**begin repeat1
- * #kind = fmax, fmin#
- * #OP = CGE, CLE#
- */
-static void
-C@TYPE@_@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 (@OP@(in1r, in1i, in2r, in2i) || isnan(in2r) || isnan(in2i)) {
-            ((@type@ *)op1)[0] = in1r;
-            ((@type@ *)op1)[1] = in1i;
-        }
-        else {
-            ((@type@ *)op1)[0] = in2r;
-            ((@type@ *)op1)[1] = in2i;
-        }
-    }
-}
-/**end repeat1**/
-
-#define C@TYPE@_true_divide C@TYPE@_divide
-
-/**end repeat**/
-
-#undef CGE
-#undef CLE
-#undef CGT
-#undef CLT
-#undef CEQ
-#undef CNE
-
-/*
- *****************************************************************************
- **                            OBJECT LOOPS                                 **
- *****************************************************************************
- */
-
-/**begin repeat
- * #kind = equal, not_equal, greater, greater_equal, less, less_equal#
- * #OP = EQ, NE, GT, GE, LT, LE#
- */
-static void
-OBJECT_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) {
-    BINARY_LOOP {
-        PyObject *in1 = *(PyObject **)ip1;
-        PyObject *in2 = *(PyObject **)ip2;
-        int ret = PyObject_RichCompareBool(in1, in2, Py_@OP@);
-        if (ret == -1) {
-            return;
-        }
-        *((Bool *)op1) = (Bool)ret;
-    }
-}
-/**end repeat**/
-
-static void
-OBJECT_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
-    PyObject *zero = PyInt_FromLong(0);
-    UNARY_LOOP {
-        PyObject *in1 = *(PyObject **)ip1;
-        PyObject **out = (PyObject **)op1;
-        PyObject *ret = PyInt_FromLong(PyObject_Compare(in1, zero));
-        if (PyErr_Occurred()) {
-            return;
-        }
-        Py_XDECREF(*out);
-        *out = ret;
-    }
-    Py_DECREF(zero);
-}
 
 /*
  *****************************************************************************
- **                              END LOOPS                                  **
+ **                    INCLUDE GENERATED CODE                               **
  *****************************************************************************
  */
+#include "umath_funcs_c99.inc"
+#include "umath_funcs.inc"
+#include "umath_loops.inc"
+#include "umath_ufunc_object.inc"
+#include "__umath_generated.c"
+#include "__ufunc_api.c"
 
 
 /*
@@ -1938,9 +36,7 @@ OBJECT_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
  *****************************************************************************
  */
 
-#include "__umath_generated.c"
-#include "umath_ufunc_object.inc"
-#include "__ufunc_api.c"
+/* Less automated additions to the ufuncs */
 
 static PyUFuncGenericFunction frexp_functions[] = {
 #ifdef HAVE_FREXPF
@@ -1983,41 +79,6 @@ static char ldexp_signatures[] = {
 #endif
 };
 
-
-static double
-pinf_init(void)
-{
-    double mul = 1e10;
-    double tmp = 0.0;
-    double pinf;
-
-    pinf = mul;
-    for (;;) {
-        pinf *= mul;
-        if (pinf == tmp) break;
-        tmp = pinf;
-    }
-    return pinf;
-}
-
-static double
-pzero_init(void)
-{
-    double div = 1e10;
-    double tmp = 0.0;
-    double pinf;
-
-    pinf = div;
-    for (;;) {
-        pinf /= div;
-        if (pinf == tmp) break;
-        tmp = pinf;
-    }
-    return pinf;
-}
-
-/* Less automated additions to the ufuncs */
-
 static void
 InitOtherOperators(PyObject *dictionary) {
     PyObject *f;
@@ -2052,6 +113,42 @@ InitOtherOperators(PyObject *dictionary) {
     return;
 }
 
+/* Setup +inf and +0 */ 
+
+static double
+pinf_init(void)
+{
+    double mul = 1e10;
+    double tmp = 0.0;
+    double pinf;
+
+    pinf = mul;
+    for (;;) {
+        pinf *= mul;
+        if (pinf == tmp) break;
+        tmp = pinf;
+    }
+    return pinf;
+}
+
+static double
+pzero_init(void)
+{
+    double div = 1e10;
+    double tmp = 0.0;
+    double pinf;
+
+    pinf = div;
+    for (;;) {
+        pinf /= div;
+        if (pinf == tmp) break;
+        tmp = pinf;
+    }
+    return pinf;
+}
+
+/* Setup the umath module */
+
 static struct PyMethodDef methods[] = {
     {"frompyfunc", (PyCFunction) ufunc_frompyfunc,
      METH_VARARGS | METH_KEYWORDS, doc_frompyfunc},