[Numpy-svn] r6090 - in trunk/numpy/core: . code_generators src
numpy-svn at scipy.org
numpy-svn at scipy.org
Fri Nov 21 23:25:32 EST 2008
Author: charris
Date: 2008-11-21 22:25:21 -0600 (Fri, 21 Nov 2008)
New Revision: 6090
Added:
trunk/numpy/core/src/umath_funcs.inc.src
trunk/numpy/core/src/umath_loops.inc.src
Modified:
trunk/numpy/core/SConscript
trunk/numpy/core/code_generators/genapi.py
trunk/numpy/core/setup.py
trunk/numpy/core/src/umathmodule.c.src
Log:
Merge branch 'ufunc'
Modified: trunk/numpy/core/SConscript
===================================================================
--- trunk/numpy/core/SConscript 2008-11-22 01:28:52 UTC (rev 6089)
+++ trunk/numpy/core/SConscript 2008-11-22 04:25:21 UTC (rev 6090)
@@ -137,9 +137,9 @@
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 @@
# 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 @@
# 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'))
Modified: trunk/numpy/core/code_generators/genapi.py
===================================================================
--- trunk/numpy/core/code_generators/genapi.py 2008-11-22 01:28:52 UTC (rev 6089)
+++ trunk/numpy/core/code_generators/genapi.py 2008-11-22 04:25:21 UTC (rev 6090)
@@ -18,6 +18,7 @@
'multiarraymodule.c',
'scalartypes.inc.src',
'umath_ufunc_object.inc',
+ 'umath_funcs.inc.src',
'umathmodule.c.src'
]
THIS_DIR = os.path.dirname(__file__)
Modified: trunk/numpy/core/setup.py
===================================================================
--- trunk/numpy/core/setup.py 2008-11-22 01:28:52 UTC (rev 6089)
+++ trunk/numpy/core/setup.py 2008-11-22 04:25:21 UTC (rev 6090)
@@ -358,6 +358,8 @@
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,
Added: trunk/numpy/core/src/umath_funcs.inc.src
===================================================================
--- trunk/numpy/core/src/umath_funcs.inc.src 2008-11-22 01:28:52 UTC (rev 6089)
+++ trunk/numpy/core/src/umath_funcs.inc.src 2008-11-22 04:25:21 UTC (rev 6090)
@@ -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 at c@ 2.7182818284590452353602874713526625 at C@ /* e */
+#define NPY_LOG2E at c@ 1.4426950408889634073599246810018921 at C@ /* log_2 e */
+#define NPY_LOG10E at c@ 0.4342944819032518276511289189166051 at C@ /* log_10 e */
+#define NPY_LOGE2 at c@ 0.6931471805599453094172321214581766 at C@ /* log_e 2 */
+#define NPY_LOGE10 at c@ 2.3025850929940456840179914546843642 at C@ /* log_e 10 */
+#define NPY_PI at c@ 3.1415926535897932384626433832795029 at C@ /* pi */
+#define NPY_PI_2 at c@ 1.5707963267948966192313216916397514 at C@ /* pi/2 */
+#define NPY_PI_4 at c@ 0.7853981633974483096156608458198757 at C@ /* pi/4 */
+#define NPY_1_PI at c@ 0.3183098861837906715377675267450287 at C@ /* 1/pi */
+#define NPY_2_PI at c@ 0.6366197723675813430755350534900574 at C@ /* 2/pi */
+/**end repeat**/
+
+/*
+ ******************************************************************************
+ ** FLOAT FUNCTIONS **
+ ******************************************************************************
+ */
+
+/**begin repeat
+ * #type = float, double, longdouble#
+ * #c = f, ,l#
+ * #C = F, ,L#
+ */
+
+#define LOGE2 NPY_LOGE2 at c@
+#define LOG2E NPY_LOG2E at c@
+#define RAD2DEG (180.0 at c@/NPY_PI at c@)
+#define DEG2RAD (NPY_PI at c@/180.0 at c@)
+
+static @type@
+rad2deg at c@(@type@ x) {
+ return x*RAD2DEG;
+}
+
+static @type@
+deg2rad at c@(@type@ x) {
+ return x*DEG2RAD;
+}
+
+static @type@
+log2_1p at c@(@type@ x)
+{
+ @type@ u = 1 + x;
+ if (u == 1) {
+ return LOG2E*x;
+ } else {
+ return log2 at c@(u) * x / (u - 1);
+ }
+}
+
+static @type@
+exp2_1m at c@(@type@ x)
+{
+ @type@ u = exp at c@(x);
+ if (u == 1.0) {
+ return LOGE2*x;
+ } else if (u - 1 == -1) {
+ return -LOGE2;
+ } else {
+ return (u - 1) * x/log2 at c@(u);
+ }
+}
+
+static @type@
+logaddexp at c@(@type@ x, @type@ y)
+{
+ const @type@ tmp = x - y;
+ if (tmp > 0) {
+ return x + log1p at c@(exp at c@(-tmp));
+ }
+ else {
+ return y + log1p at c@(exp at c@(tmp));
+ }
+}
+
+static @type@
+logaddexp2 at c@(@type@ x, @type@ y)
+{
+ const @type@ tmp = x - y;
+ if (tmp > 0) {
+ return x + log2_1p at c@(exp2 at c@(-tmp));
+ }
+ else {
+ return y + log2_1p at c@(exp2 at c@(tmp));
+ }
+}
+
+#define degrees at c@ rad2deg at c@
+#define radians at c@ deg2rad at 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 at 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 at typ@ nc_1 at c@ = {1., 0.};
+static c at typ@ nc_half at c@ = {0.5, 0.};
+static c at typ@ nc_i at c@ = {0., 1.};
+static c at typ@ nc_i2 at c@ = {0., 0.5};
+/*
+ * static c at typ@ nc_mi at c@ = {0., -1.};
+ * static c at typ@ nc_pi2 at c@ = {M_PI/2., 0.};
+ */
+
+
+static void
+nc_sum at c@(c at typ@ *a, c at typ@ *b, c at typ@ *r)
+{
+ r->real = a->real + b->real;
+ r->imag = a->imag + b->imag;
+ return;
+}
+
+static void
+nc_diff at c@(c at typ@ *a, c at typ@ *b, c at typ@ *r)
+{
+ r->real = a->real - b->real;
+ r->imag = a->imag - b->imag;
+ return;
+}
+
+static void
+nc_neg at c@(c at typ@ *a, c at typ@ *r)
+{
+ r->real = -a->real;
+ r->imag = -a->imag;
+ return;
+}
+
+static void
+nc_prod at c@(c at typ@ *a, c at typ@ *b, c at 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 at c@(c at typ@ *a, c at typ@ *b, c at 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 at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ s,d;
+ if (x->real == 0. && x->imag == 0.)
+ *r = *x;
+ else {
+ s = sqrt at c@((fabs at c@(x->real) + hypot at 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 at c@(c at typ@ *x, c at typ@ *r)
+{
+ r->real = rint at c@(x->real);
+ r->imag = rint at c@(x->imag);
+}
+
+static void
+nc_log at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ l = hypot at c@(x->real,x->imag);
+ r->imag = atan2 at c@(x->imag, x->real);
+ r->real = log at c@(l);
+ return;
+}
+
+static void
+nc_log1p at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ l = hypot at c@(x->real + 1,x->imag);
+ r->imag = atan2 at c@(x->imag, x->real + 1);
+ r->real = log at c@(l);
+ return;
+}
+
+static void
+nc_exp at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ a = exp at c@(x->real);
+ r->real = a*cos at c@(x->imag);
+ r->imag = a*sin at c@(x->imag);
+ return;
+}
+
+static void
+nc_expm1 at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ a = exp at c@(x->real);
+ r->real = a*cos at c@(x->imag) - 1;
+ r->imag = a*sin at c@(x->imag);
+ return;
+}
+
+static void
+nc_pow at c@(c at typ@ *a, c at typ@ *b, c at 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 at typ@ p, aa;
+ intp mask = 1;
+ if (n < 0) n = -n;
+ aa = nc_1 at c@;
+ p.real = ar; p.imag = ai;
+ while (1) {
+ if (n & mask)
+ nc_prod at c@(&aa,&p,&aa);
+ mask <<= 1;
+ if (n < mask || mask <= 0) break;
+ nc_prod at c@(&p,&p,&p);
+ }
+ r->real = aa.real; r->imag = aa.imag;
+ if (br < 0) nc_quot at c@(&nc_1 at c@, r, r);
+ return;
+ }
+ }
+ /*
+ * complexobect.c uses an inline version of this formula
+ * investigate whether this had better performance or accuracy
+ */
+ nc_log at c@(a, r);
+ nc_prod at c@(r, b, r);
+ nc_exp at c@(r, r);
+ return;
+}
+
+
+static void
+nc_prodi at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ xr = x->real;
+ r->real = -x->imag;
+ r->imag = xr;
+ return;
+}
+
+
+static void
+nc_acos at c@(c at typ@ *x, c at 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 at c@(x,x,r);
+ nc_diff at c@(&nc_1 at c@, r, r);
+ nc_sqrt at c@(r, r);
+ nc_prodi at c@(r, r);
+ nc_sum at c@(x, r, r);
+ nc_log at c@(r, r);
+ nc_prodi at c@(r, r);
+ nc_neg at c@(r, r);
+ return;
+}
+
+static void
+nc_acosh at c@(c at typ@ *x, c at 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 at typ@ t;
+
+ nc_sum at c@(x, &nc_1 at c@, &t);
+ nc_sqrt at c@(&t, &t);
+ nc_diff at c@(x, &nc_1 at c@, r);
+ nc_sqrt at c@(r, r);
+ nc_prod at c@(&t, r, r);
+ nc_sum at c@(x, r, r);
+ nc_log at c@(r, r);
+ return;
+}
+
+static void
+nc_asin at c@(c at typ@ *x, c at 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 at typ@ a, *pa=&a;
+ nc_prod at c@(x, x, r);
+ nc_diff at c@(&nc_1 at c@, r, r);
+ nc_sqrt at c@(r, r);
+ nc_prodi at c@(x, pa);
+ nc_sum at c@(pa, r, r);
+ nc_log at c@(r, r);
+ nc_prodi at c@(r, r);
+ nc_neg at c@(r, r);
+ return;
+}
+
+
+static void
+nc_asinh at c@(c at typ@ *x, c at typ@ *r)
+{
+ /*
+ * return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x));
+ */
+ nc_prod at c@(x, x, r);
+ nc_sum at c@(&nc_1 at c@, r, r);
+ nc_sqrt at c@(r, r);
+ nc_sum at c@(r, x, r);
+ nc_log at c@(r, r);
+ return;
+}
+
+static void
+nc_atan at c@(c at typ@ *x, c at typ@ *r)
+{
+ /*
+ * return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
+ */
+ c at typ@ a, *pa=&a;
+ nc_diff at c@(&nc_i at c@, x, pa);
+ nc_sum at c@(&nc_i at c@, x, r);
+ nc_quot at c@(r, pa, r);
+ nc_log at c@(r,r);
+ nc_prod at c@(&nc_i2 at c@, r, r);
+ return;
+}
+
+static void
+nc_atanh at c@(c at typ@ *x, c at typ@ *r)
+{
+ /*
+ * return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
+ */
+ c at typ@ a, *pa=&a;
+ nc_diff at c@(&nc_1 at c@, x, r);
+ nc_sum at c@(&nc_1 at c@, x, pa);
+ nc_quot at c@(pa, r, r);
+ nc_log at c@(r, r);
+ nc_prod at c@(&nc_half at c@, r, r);
+ return;
+}
+
+static void
+nc_cos at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ xr=x->real, xi=x->imag;
+ r->real = cos at c@(xr)*cosh at c@(xi);
+ r->imag = -sin at c@(xr)*sinh at c@(xi);
+ return;
+}
+
+static void
+nc_cosh at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ xr=x->real, xi=x->imag;
+ r->real = cos at c@(xi)*cosh at c@(xr);
+ r->imag = sin at c@(xi)*sinh at c@(xr);
+ return;
+}
+
+static void
+nc_log10 at c@(c at typ@ *x, c at typ@ *r)
+{
+ nc_log at c@(x, r);
+ r->real *= (@typ@) M_LOG10_E;
+ r->imag *= (@typ@) M_LOG10_E;
+ return;
+}
+
+static void
+nc_sin at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ xr=x->real, xi=x->imag;
+ r->real = sin at c@(xr)*cosh at c@(xi);
+ r->imag = cos at c@(xr)*sinh at c@(xi);
+ return;
+}
+
+static void
+nc_sinh at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ xr=x->real, xi=x->imag;
+ r->real = cos at c@(xi)*sinh at c@(xr);
+ r->imag = sin at c@(xi)*cosh at c@(xr);
+ return;
+}
+
+static void
+nc_tan at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ sr,cr,shi,chi;
+ @typ@ rs,is,rc,ic;
+ @typ@ d;
+ @typ@ xr=x->real, xi=x->imag;
+ sr = sin at c@(xr);
+ cr = cos at c@(xr);
+ shi = sinh at c@(xi);
+ chi = cosh at 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 at c@(c at typ@ *x, c at typ@ *r)
+{
+ @typ@ si,ci,shr,chr;
+ @typ@ rs,is,rc,ic;
+ @typ@ d;
+ @typ@ xr=x->real, xi=x->imag;
+ si = sin at c@(xi);
+ ci = cos at c@(xi);
+ shr = sinh at c@(xr);
+ chr = cosh at 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**/
+
Copied: trunk/numpy/core/src/umath_loops.inc.src (from rev 6089, trunk/numpy/core/src/umathmodule.c.src)
===================================================================
--- trunk/numpy/core/src/umathmodule.c.src 2008-11-22 01:28:52 UTC (rev 6089)
+++ trunk/numpy/core/src/umath_loops.inc.src 2008-11-22 04:25:21 UTC (rev 6090)
@@ -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_ at 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_ at 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_ at 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 at _floor_divide @S@@TYPE at _divide
+#define @S@@TYPE at _fmax @S@@TYPE at _maximum
+#define @S@@TYPE at _fmin @S@@TYPE at _minimum
+
+static void
+ at S@@TYPE at _ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+ OUTPUT_LOOP {
+ *((@s@@type@ *)op1) = 1;
+ }
+}
+
+static void
+ at S@@TYPE at _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
+ at S@@TYPE at _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
+ at S@@TYPE at _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
+ at S@@TYPE at _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
+ at S@@TYPE at _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
+ at S@@TYPE at _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
+ at S@@TYPE at _@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
+ at S@@TYPE at _@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
+ at S@@TYPE at _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
+ at S@@TYPE at _@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
+ at S@@TYPE at _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
+ at S@@TYPE at _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
+ at S@@TYPE at _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 at TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ UNARY_LOOP {
+ const u at type@ in1 = *(u at type@ *)ip1;
+ *((u at type@ *)op1) = in1;
+ }
+}
+
+static void
+ at 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 at TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ UNARY_LOOP {
+ const u at type@ in1 = *(u at type@ *)ip1;
+ *((u at type@ *)op1) = in1 > 0 ? 1 : 0;
+ }
+}
+
+static void
+ at 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
+ at 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 at TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ BINARY_LOOP {
+ const u at type@ in1 = *(u at type@ *)ip1;
+ const u at type@ in2 = *(u at type@ *)ip2;
+ if (in2 == 0) {
+ generate_divbyzero_error();
+ *((u at type@ *)op1) = 0;
+ }
+ else {
+ *((u at type@ *)op1)= in1/in2;
+ }
+ }
+}
+
+static void
+ at 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 at TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ BINARY_LOOP {
+ const u at type@ in1 = *(u at type@ *)ip1;
+ const u at type@ in2 = *(u at 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
+ at TYPE@_ at 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
+ at TYPE@_ at 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
+ at 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
+ at 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
+ at TYPE@_ at 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
+ at TYPE@_ at 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
+ at TYPE@_ at 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
+ at 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 at c@(in1/in2);
+ }
+}
+
+static void
+ at 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 at c@(in1,in2);
+ if (res && ((in2 < 0) != (res < 0))) {
+ *((@type@ *)op1) = res + in2;
+ }
+ else {
+ *((@type@ *)op1) = res;
+ }
+ }
+}
+
+static void
+ at 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
+ at 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
+ at TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
+{
+ OUTPUT_LOOP {
+ *((@type@ *)op1) = 1;
+ }
+}
+
+static void
+ at TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ UNARY_LOOP {
+ const @type@ in1 = *(@type@ *)ip1;
+ *((@type@ *)op1) = in1;
+ }
+}
+
+static void
+ at 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
+ at TYPE@_negative(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ UNARY_LOOP {
+ const @type@ in1 = *(@type@ *)ip1;
+ *((@type@ *)op1) = -in1;
+ }
+}
+
+static void
+ at 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
+ at TYPE@_modf(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ UNARY_LOOP_TWO_OUT {
+ const @type@ in1 = *(@type@ *)ip1;
+ *((@type@ *)op1) = modf at c@(in1, (@type@ *)op2);
+ }
+}
+
+#ifdef HAVE_FREXP at C@
+static void
+ at TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
+{
+ UNARY_LOOP_TWO_OUT {
+ const @type@ in1 = *(@type@ *)ip1;
+ *((@type@ *)op1) = frexp at c@(in1, (int *)op2);
+ }
+}
+#endif
+
+#ifdef HAVE_LDEXP at C@
+static void
+ at 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 at c@(in1, in2);
+ }
+}
+#endif
+
+#define @TYPE at _true_divide @TYPE at _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 at TYPE@_ at 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 at 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 at 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 at 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 at 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 at TYPE@_ at 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 at TYPE@_ at 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 at 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 at 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 at TYPE@_ at 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 at 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 at 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 at c@(in1i) <= fabs at 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 at 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 at 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 at 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 at c@(in1r*in1r + in1i*in1i);
+ }
+}
+
+static void
+C at 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 at TYPE@_ at 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 at TYPE@_ at 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 at TYPE@_true_divide C at 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_ at 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_ at 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 **
+ *****************************************************************************
+ */
+
+
Modified: trunk/numpy/core/src/umathmodule.c.src
===================================================================
--- trunk/numpy/core/src/umathmodule.c.src 2008-11-22 01:28:52 UTC (rev 6089)
+++ trunk/numpy/core/src/umathmodule.c.src 2008-11-22 04:25:21 UTC (rev 6090)
@@ -9,1938 +9,34 @@
** 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 **
+ ** 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"
-/**begin repeat
- * #type = float, double, longdouble#
- * #c = f, ,l#
- * #C = F, ,L#
- */
-/* fixme: need more precision for LOG2 and INVLOG2 */
-
-#define PI 3.14159265358979323846264338328 at c@
-#define LOG2 0.69314718055994530943 at c@
-#define INVLOG2 1.4426950408889634074 at c@
-#define degrees at c@ rad2deg at c@
-#define radians at c@ deg2rad at c@
-
-static @type@
-rad2deg at c@(@type@ x) {
- return x*(180.0 at c@/PI);
-}
-
-static @type@
-deg2rad at c@(@type@ x) {
- return x*(PI/180.0 at c@);
-}
-
-static @type@
-log2_1p at c@(@type@ x)
-{
- @type@ u = 1 + x;
- if (u == 1) {
- return INVLOG2*x;
- } else {
- return log2 at c@(u) * x / (u - 1);
- }
-}
-
-static @type@
-exp2_1m at c@(@type@ x)
-{
- @type@ u = exp at c@(x);
- if (u == 1.0) {
- return LOG2*x;
- } else if (u - 1 == -1) {
- return -LOG2;
- } else {
- return (u - 1) * x/log2 at c@(u);
- }
-}
-
-static @type@
-logaddexp at c@(@type@ x, @type@ y)
-{
- const @type@ tmp = x - y;
- if (tmp > 0) {
- return x + log1p at c@(exp at c@(-tmp));
- }
- else {
- return y + log1p at c@(exp at c@(tmp));
- }
-}
-
-static @type@
-logaddexp2 at c@(@type@ x, @type@ y)
-{
- const @type@ tmp = x - y;
- if (tmp > 0) {
- return x + log2_1p at c@(exp2 at c@(-tmp));
- }
- else {
- return y + log2_1p at c@(exp2 at 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 at 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 at typ@ nc_1 at c@ = {1., 0.};
-static c at typ@ nc_half at c@ = {0.5, 0.};
-static c at typ@ nc_i at c@ = {0., 1.};
-static c at typ@ nc_i2 at c@ = {0., 0.5};
-/*
- * static c at typ@ nc_mi at c@ = {0., -1.};
- * static c at typ@ nc_pi2 at c@ = {M_PI/2., 0.};
- */
-
-
-static void
-nc_sum at c@(c at typ@ *a, c at typ@ *b, c at typ@ *r)
-{
- r->real = a->real + b->real;
- r->imag = a->imag + b->imag;
- return;
-}
-
-static void
-nc_diff at c@(c at typ@ *a, c at typ@ *b, c at typ@ *r)
-{
- r->real = a->real - b->real;
- r->imag = a->imag - b->imag;
- return;
-}
-
-static void
-nc_neg at c@(c at typ@ *a, c at typ@ *r)
-{
- r->real = -a->real;
- r->imag = -a->imag;
- return;
-}
-
-static void
-nc_prod at c@(c at typ@ *a, c at typ@ *b, c at 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 at c@(c at typ@ *a, c at typ@ *b, c at 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 at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ s,d;
- if (x->real == 0. && x->imag == 0.)
- *r = *x;
- else {
- s = sqrt at c@((fabs at c@(x->real) + hypot at 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 at c@(c at typ@ *x, c at typ@ *r)
-{
- r->real = rint at c@(x->real);
- r->imag = rint at c@(x->imag);
-}
-
-static void
-nc_log at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ l = hypot at c@(x->real,x->imag);
- r->imag = atan2 at c@(x->imag, x->real);
- r->real = log at c@(l);
- return;
-}
-
-static void
-nc_log1p at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ l = hypot at c@(x->real + 1,x->imag);
- r->imag = atan2 at c@(x->imag, x->real + 1);
- r->real = log at c@(l);
- return;
-}
-
-static void
-nc_exp at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ a = exp at c@(x->real);
- r->real = a*cos at c@(x->imag);
- r->imag = a*sin at c@(x->imag);
- return;
-}
-
-static void
-nc_expm1 at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ a = exp at c@(x->real);
- r->real = a*cos at c@(x->imag) - 1;
- r->imag = a*sin at c@(x->imag);
- return;
-}
-
-static void
-nc_pow at c@(c at typ@ *a, c at typ@ *b, c at 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 at typ@ p, aa;
- intp mask = 1;
- if (n < 0) n = -n;
- aa = nc_1 at c@;
- p.real = ar; p.imag = ai;
- while (1) {
- if (n & mask)
- nc_prod at c@(&aa,&p,&aa);
- mask <<= 1;
- if (n < mask || mask <= 0) break;
- nc_prod at c@(&p,&p,&p);
- }
- r->real = aa.real; r->imag = aa.imag;
- if (br < 0) nc_quot at c@(&nc_1 at c@, r, r);
- return;
- }
- }
- /*
- * complexobect.c uses an inline version of this formula
- * investigate whether this had better performance or accuracy
- */
- nc_log at c@(a, r);
- nc_prod at c@(r, b, r);
- nc_exp at c@(r, r);
- return;
-}
-
-
-static void
-nc_prodi at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ xr = x->real;
- r->real = -x->imag;
- r->imag = xr;
- return;
-}
-
-
-static void
-nc_acos at c@(c at typ@ *x, c at 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 at c@(x,x,r);
- nc_diff at c@(&nc_1 at c@, r, r);
- nc_sqrt at c@(r, r);
- nc_prodi at c@(r, r);
- nc_sum at c@(x, r, r);
- nc_log at c@(r, r);
- nc_prodi at c@(r, r);
- nc_neg at c@(r, r);
- return;
-}
-
-static void
-nc_acosh at c@(c at typ@ *x, c at 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 at typ@ t;
-
- nc_sum at c@(x, &nc_1 at c@, &t);
- nc_sqrt at c@(&t, &t);
- nc_diff at c@(x, &nc_1 at c@, r);
- nc_sqrt at c@(r, r);
- nc_prod at c@(&t, r, r);
- nc_sum at c@(x, r, r);
- nc_log at c@(r, r);
- return;
-}
-
-static void
-nc_asin at c@(c at typ@ *x, c at 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 at typ@ a, *pa=&a;
- nc_prod at c@(x, x, r);
- nc_diff at c@(&nc_1 at c@, r, r);
- nc_sqrt at c@(r, r);
- nc_prodi at c@(x, pa);
- nc_sum at c@(pa, r, r);
- nc_log at c@(r, r);
- nc_prodi at c@(r, r);
- nc_neg at c@(r, r);
- return;
-}
-
-
-static void
-nc_asinh at c@(c at typ@ *x, c at typ@ *r)
-{
- /*
- * return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x));
- */
- nc_prod at c@(x, x, r);
- nc_sum at c@(&nc_1 at c@, r, r);
- nc_sqrt at c@(r, r);
- nc_sum at c@(r, x, r);
- nc_log at c@(r, r);
- return;
-}
-
-static void
-nc_atan at c@(c at typ@ *x, c at typ@ *r)
-{
- /*
- * return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
- */
- c at typ@ a, *pa=&a;
- nc_diff at c@(&nc_i at c@, x, pa);
- nc_sum at c@(&nc_i at c@, x, r);
- nc_quot at c@(r, pa, r);
- nc_log at c@(r,r);
- nc_prod at c@(&nc_i2 at c@, r, r);
- return;
-}
-
-static void
-nc_atanh at c@(c at typ@ *x, c at typ@ *r)
-{
- /*
- * return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
- */
- c at typ@ a, *pa=&a;
- nc_diff at c@(&nc_1 at c@, x, r);
- nc_sum at c@(&nc_1 at c@, x, pa);
- nc_quot at c@(pa, r, r);
- nc_log at c@(r, r);
- nc_prod at c@(&nc_half at c@, r, r);
- return;
-}
-
-static void
-nc_cos at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ xr=x->real, xi=x->imag;
- r->real = cos at c@(xr)*cosh at c@(xi);
- r->imag = -sin at c@(xr)*sinh at c@(xi);
- return;
-}
-
-static void
-nc_cosh at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ xr=x->real, xi=x->imag;
- r->real = cos at c@(xi)*cosh at c@(xr);
- r->imag = sin at c@(xi)*sinh at c@(xr);
- return;
-}
-
-
-#define M_LOG10_E 0.434294481903251827651128918916605082294397
-
-static void
-nc_log10 at c@(c at typ@ *x, c at typ@ *r)
-{
- nc_log at c@(x, r);
- r->real *= (@typ@) M_LOG10_E;
- r->imag *= (@typ@) M_LOG10_E;
- return;
-}
-
-static void
-nc_sin at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ xr=x->real, xi=x->imag;
- r->real = sin at c@(xr)*cosh at c@(xi);
- r->imag = cos at c@(xr)*sinh at c@(xi);
- return;
-}
-
-static void
-nc_sinh at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ xr=x->real, xi=x->imag;
- r->real = cos at c@(xi)*sinh at c@(xr);
- r->imag = sin at c@(xi)*cosh at c@(xr);
- return;
-}
-
-static void
-nc_tan at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ sr,cr,shi,chi;
- @typ@ rs,is,rc,ic;
- @typ@ d;
- @typ@ xr=x->real, xi=x->imag;
- sr = sin at c@(xr);
- cr = cos at c@(xr);
- shi = sinh at c@(xi);
- chi = cosh at 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 at c@(c at typ@ *x, c at typ@ *r)
-{
- @typ@ si,ci,shr,chr;
- @typ@ rs,is,rc,ic;
- @typ@ d;
- @typ@ xr=x->real, xi=x->imag;
- si = sin at c@(xi);
- ci = cos at c@(xi);
- shr = sinh at c@(xr);
- chr = cosh at 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_ at 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_ at 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_ at 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 at _floor_divide @S@@TYPE at _divide
-#define @S@@TYPE at _fmax @S@@TYPE at _maximum
-#define @S@@TYPE at _fmin @S@@TYPE at _minimum
-
-static void
- at S@@TYPE at _ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
- OUTPUT_LOOP {
- *((@s@@type@ *)op1) = 1;
- }
-}
-
-static void
- at S@@TYPE at _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
- at S@@TYPE at _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
- at S@@TYPE at _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
- at S@@TYPE at _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
- at S@@TYPE at _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
- at S@@TYPE at _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
- at S@@TYPE at _@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
- at S@@TYPE at _@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
- at S@@TYPE at _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
- at S@@TYPE at _@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
- at S@@TYPE at _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
- at S@@TYPE at _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
- at S@@TYPE at _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 at TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- UNARY_LOOP {
- const u at type@ in1 = *(u at type@ *)ip1;
- *((u at type@ *)op1) = in1;
- }
-}
-
-static void
- at 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 at TYPE@_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- UNARY_LOOP {
- const u at type@ in1 = *(u at type@ *)ip1;
- *((u at type@ *)op1) = in1 > 0 ? 1 : 0;
- }
-}
-
-static void
- at 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
- at 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 at TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- BINARY_LOOP {
- const u at type@ in1 = *(u at type@ *)ip1;
- const u at type@ in2 = *(u at type@ *)ip2;
- if (in2 == 0) {
- generate_divbyzero_error();
- *((u at type@ *)op1) = 0;
- }
- else {
- *((u at type@ *)op1)= in1/in2;
- }
- }
-}
-
-static void
- at 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 at TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- BINARY_LOOP {
- const u at type@ in1 = *(u at type@ *)ip1;
- const u at type@ in2 = *(u at 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
- at TYPE@_ at 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
- at TYPE@_ at 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
- at 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
- at 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
- at TYPE@_ at 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
- at TYPE@_ at 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
- at TYPE@_ at 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
- at 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 at c@(in1/in2);
- }
-}
-
-static void
- at 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 at c@(in1,in2);
- if (res && ((in2 < 0) != (res < 0))) {
- *((@type@ *)op1) = res + in2;
- }
- else {
- *((@type@ *)op1) = res;
- }
- }
-}
-
-static void
- at 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
- at 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
- at TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
-{
- OUTPUT_LOOP {
- *((@type@ *)op1) = 1;
- }
-}
-
-static void
- at TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- UNARY_LOOP {
- const @type@ in1 = *(@type@ *)ip1;
- *((@type@ *)op1) = in1;
- }
-}
-
-static void
- at 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
- at TYPE@_negative(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- UNARY_LOOP {
- const @type@ in1 = *(@type@ *)ip1;
- *((@type@ *)op1) = -in1;
- }
-}
-
-static void
- at 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
- at TYPE@_modf(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- UNARY_LOOP_TWO_OUT {
- const @type@ in1 = *(@type@ *)ip1;
- *((@type@ *)op1) = modf at c@(in1, (@type@ *)op2);
- }
-}
-
-#ifdef HAVE_FREXP at C@
-static void
- at TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
-{
- UNARY_LOOP_TWO_OUT {
- const @type@ in1 = *(@type@ *)ip1;
- *((@type@ *)op1) = frexp at c@(in1, (int *)op2);
- }
-}
-#endif
-
-#ifdef HAVE_LDEXP at C@
-static void
- at 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 at c@(in1, in2);
- }
-}
-#endif
-
-#define @TYPE at _true_divide @TYPE at _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 at TYPE@_ at 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 at 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 at 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 at 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 at 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 at TYPE@_ at 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 at TYPE@_ at 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 at 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 at 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 at TYPE@_ at 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 at 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 at 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 at c@(in1i) <= fabs at 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 at 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 at 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 at 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 at c@(in1r*in1r + in1i*in1i);
- }
-}
-
-static void
-C at 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 at TYPE@_ at 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 at TYPE@_ at 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 at TYPE@_true_divide C at 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_ at 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_ at 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 **
- *****************************************************************************
- */
-
-
-/*
- *****************************************************************************
** SETUP UFUNCS **
*****************************************************************************
*/
-#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 @@
#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 @@
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},
More information about the Numpy-svn
mailing list