/* -*- c -*- */ /* * vim:syntax=c */ /* ***************************************************************************** ** INCLUDES ** ***************************************************************************** */ #include "Python.h" #include "numpy/noprefix.h" #define _UMATHMODULE #include "numpy/ufuncobject.h" #include "abstract.h" #include "config.h" #include #ifndef M_PI #define M_PI 3.14159265358979323846264338328 #endif #include "math_c99.inc" float degreesf(float x) { return x * (float)(180.0/M_PI); } double degrees(double x) { return x * (180.0/M_PI); } longdouble degreesl(longdouble x) { return x * (180.0L/M_PI); } float radiansf(float x) { return x * (float)(M_PI/180.0); } double radians(double x) { return x * (M_PI/180.0); } longdouble radiansl(longdouble x) { return x * (M_PI/180.0L); } /* ***************************************************************************** ** PYTHON OBJECT FUNCTIONS ** ***************************************************************************** */ static PyObject * Py_square(PyObject *o) { return PyNumber_Multiply(o, o); } static PyObject * Py_get_one(PyObject *o) { return PyInt_FromLong(1); } static PyObject * Py_reciprocal(PyObject *o) { PyObject *one = PyInt_FromLong(1); PyObject *result; if (!one) { return NULL; } result = PyNumber_Divide(one, o); Py_DECREF(one); return result; } /**begin repeat * #Kind = Max, Min# * #OP = >=, <=# */ static PyObject * _npy_Object@Kind@(PyObject *i1, PyObject *i2) { PyObject *result; int cmp; if (PyObject_Cmp(i1, i2, &cmp) < 0) { return NULL; } if (cmp @OP@ 0) { result = i1; } else { result = i2; } Py_INCREF(result); return result; } /**end repeat**/ /* ***************************************************************************** ** COMPLEX FUNCTIONS ** ***************************************************************************** */ /* * Don't pass structures between functions (only pointers) because how * structures are passed is compiler dependent and could cause * segfaults if ufuncobject.c is compiled with a different compiler * than an extension that makes use of the UFUNC API */ /**begin repeat #typ=float, double, longdouble# #c=f,,l# */ /* constants */ static c@typ@ nc_1@c@ = {1., 0.}; static c@typ@ nc_half@c@ = {0.5, 0.}; static c@typ@ nc_i@c@ = {0., 1.}; static c@typ@ nc_i2@c@ = {0., 0.5}; /* * static c@typ@ nc_mi@c@ = {0., -1.}; * static c@typ@ nc_pi2@c@ = {M_PI/2., 0.}; */ static void nc_sum@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) { r->real = a->real + b->real; r->imag = a->imag + b->imag; return; } static void nc_diff@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) { r->real = a->real - b->real; r->imag = a->imag - b->imag; return; } static void nc_neg@c@(c@typ@ *a, c@typ@ *r) { r->real = -a->real; r->imag = -a->imag; return; } static void nc_prod@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) { @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; r->real = ar*br - ai*bi; r->imag = ar*bi + ai*br; return; } static void nc_quot@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) { @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; @typ@ d = br*br + bi*bi; r->real = (ar*br + ai*bi)/d; r->imag = (ai*br - ar*bi)/d; return; } static void nc_sqrt@c@(c@typ@ *x, c@typ@ *r) { @typ@ s,d; if (x->real == 0. && x->imag == 0.) *r = *x; else { s = sqrt@c@((fabs@c@(x->real) + hypot@c@(x->real,x->imag))/2); d = x->imag/(2*s); if (x->real > 0) { r->real = s; r->imag = d; } else if (x->imag >= 0) { r->real = d; r->imag = s; } else { r->real = -d; r->imag = -s; } } return; } static void nc_rint@c@(c@typ@ *x, c@typ@ *r) { r->real = rint@c@(x->real); r->imag = rint@c@(x->imag); } static void nc_log@c@(c@typ@ *x, c@typ@ *r) { @typ@ l = hypot@c@(x->real,x->imag); r->imag = atan2@c@(x->imag, x->real); r->real = log@c@(l); return; } static void nc_log1p@c@(c@typ@ *x, c@typ@ *r) { @typ@ l = hypot@c@(x->real + 1,x->imag); r->imag = atan2@c@(x->imag, x->real + 1); r->real = log@c@(l); return; } static void nc_exp@c@(c@typ@ *x, c@typ@ *r) { @typ@ a = exp@c@(x->real); r->real = a*cos@c@(x->imag); r->imag = a*sin@c@(x->imag); return; } static void nc_expm1@c@(c@typ@ *x, c@typ@ *r) { @typ@ a = exp@c@(x->real); r->real = a*cos@c@(x->imag) - 1; r->imag = a*sin@c@(x->imag); return; } static void nc_pow@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) { intp n; @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; if (br == 0. && bi == 0.) { r->real = 1.; r->imag = 0.; return; } if (ar == 0. && ai == 0.) { r->real = 0.; r->imag = 0.; return; } if (bi == 0 && (n=(intp)br) == br) { if (n > -100 && n < 100) { c@typ@ p, aa; intp mask = 1; if (n < 0) n = -n; aa = nc_1@c@; p.real = ar; p.imag = ai; while (1) { if (n & mask) nc_prod@c@(&aa,&p,&aa); mask <<= 1; if (n < mask || mask <= 0) break; nc_prod@c@(&p,&p,&p); } r->real = aa.real; r->imag = aa.imag; if (br < 0) nc_quot@c@(&nc_1@c@, r, r); return; } } /* * complexobect.c uses an inline version of this formula * investigate whether this had better performance or accuracy */ nc_log@c@(a, r); nc_prod@c@(r, b, r); nc_exp@c@(r, r); return; } static void nc_prodi@c@(c@typ@ *x, c@typ@ *r) { @typ@ xr = x->real; r->real = -x->imag; r->imag = xr; return; } static void nc_acos@c@(c@typ@ *x, c@typ@ *r) { /* * return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i, * nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))))); */ nc_prod@c@(x,x,r); nc_diff@c@(&nc_1@c@, r, r); nc_sqrt@c@(r, r); nc_prodi@c@(r, r); nc_sum@c@(x, r, r); nc_log@c@(r, r); nc_prodi@c@(r, r); nc_neg@c@(r, r); return; } static void nc_acosh@c@(c@typ@ *x, c@typ@ *r) { /* * return nc_log(nc_sum(x, * nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1))))); */ c@typ@ t; nc_sum@c@(x, &nc_1@c@, &t); nc_sqrt@c@(&t, &t); nc_diff@c@(x, &nc_1@c@, r); nc_sqrt@c@(r, r); nc_prod@c@(&t, r, r); nc_sum@c@(x, r, r); nc_log@c@(r, r); return; } static void nc_asin@c@(c@typ@ *x, c@typ@ *r) { /* * return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x), * nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))); */ c@typ@ a, *pa=&a; nc_prod@c@(x, x, r); nc_diff@c@(&nc_1@c@, r, r); nc_sqrt@c@(r, r); nc_prodi@c@(x, pa); nc_sum@c@(pa, r, r); nc_log@c@(r, r); nc_prodi@c@(r, r); nc_neg@c@(r, r); return; } static void nc_asinh@c@(c@typ@ *x, c@typ@ *r) { /* * return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x)); */ nc_prod@c@(x, x, r); nc_sum@c@(&nc_1@c@, r, r); nc_sqrt@c@(r, r); nc_sum@c@(r, x, r); nc_log@c@(r, r); return; } static void nc_atan@c@(c@typ@ *x, c@typ@ *r) { /* * return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x)))); */ c@typ@ a, *pa=&a; nc_diff@c@(&nc_i@c@, x, pa); nc_sum@c@(&nc_i@c@, x, r); nc_quot@c@(r, pa, r); nc_log@c@(r,r); nc_prod@c@(&nc_i2@c@, r, r); return; } static void nc_atanh@c@(c@typ@ *x, c@typ@ *r) { /* * return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x)))); */ c@typ@ a, *pa=&a; nc_diff@c@(&nc_1@c@, x, r); nc_sum@c@(&nc_1@c@, x, pa); nc_quot@c@(pa, r, r); nc_log@c@(r, r); nc_prod@c@(&nc_half@c@, r, r); return; } static void nc_cos@c@(c@typ@ *x, c@typ@ *r) { @typ@ xr=x->real, xi=x->imag; r->real = cos@c@(xr)*cosh@c@(xi); r->imag = -sin@c@(xr)*sinh@c@(xi); return; } static void nc_cosh@c@(c@typ@ *x, c@typ@ *r) { @typ@ xr=x->real, xi=x->imag; r->real = cos@c@(xi)*cosh@c@(xr); r->imag = sin@c@(xi)*sinh@c@(xr); return; } #define M_LOG10_E 0.434294481903251827651128918916605082294397 static void nc_log10@c@(c@typ@ *x, c@typ@ *r) { nc_log@c@(x, r); r->real *= (@typ@) M_LOG10_E; r->imag *= (@typ@) M_LOG10_E; return; } static void nc_sin@c@(c@typ@ *x, c@typ@ *r) { @typ@ xr=x->real, xi=x->imag; r->real = sin@c@(xr)*cosh@c@(xi); r->imag = cos@c@(xr)*sinh@c@(xi); return; } static void nc_sinh@c@(c@typ@ *x, c@typ@ *r) { @typ@ xr=x->real, xi=x->imag; r->real = cos@c@(xi)*sinh@c@(xr); r->imag = sin@c@(xi)*cosh@c@(xr); return; } static void nc_tan@c@(c@typ@ *x, c@typ@ *r) { @typ@ sr,cr,shi,chi; @typ@ rs,is,rc,ic; @typ@ d; @typ@ xr=x->real, xi=x->imag; sr = sin@c@(xr); cr = cos@c@(xr); shi = sinh@c@(xi); chi = cosh@c@(xi); rs = sr*chi; is = cr*shi; rc = cr*chi; ic = -sr*shi; d = rc*rc + ic*ic; r->real = (rs*rc+is*ic)/d; r->imag = (is*rc-rs*ic)/d; return; } static void nc_tanh@c@(c@typ@ *x, c@typ@ *r) { @typ@ si,ci,shr,chr; @typ@ rs,is,rc,ic; @typ@ d; @typ@ xr=x->real, xi=x->imag; si = sin@c@(xi); ci = cos@c@(xi); shr = sinh@c@(xr); chr = cosh@c@(xr); rs = ci*shr; is = si*chr; rc = ci*chr; ic = si*shr; d = rc*rc + ic*ic; r->real = (rs*rc+is*ic)/d; r->imag = (is*rc-rs*ic)/d; return; } /**end repeat**/ /* ***************************************************************************** ** UFUNC LOOPS ** ***************************************************************************** */ #define OUTPUT_LOOP\ char *op = args[1];\ intp os = steps[1];\ intp n = dimensions[0];\ intp i;\ for(i = 0; i < n; i++, op += os) #define UNARY_LOOP\ char *ip1 = args[0], *op = args[1];\ intp is1 = steps[0], os = steps[1];\ intp n = dimensions[0];\ intp i;\ for(i = 0; i < n; i++, ip1 += is1, op += os) #define UNARY_LOOP_TWO_OUT\ char *ip1 = args[0], *op1 = args[1], *op2 = args[2];\ intp is1 = steps[0], os1 = steps[1], os2 = steps[2];\ intp n = dimensions[0];\ intp i;\ for(i = 0; i < n; i++, ip1 += is1, op1 += os1, op2 += os2) #define BINARY_LOOP\ char *ip1 = args[0], *ip2 = args[1], *op = args[2];\ intp is1 = steps[0], is2 = steps[1], os = steps[2];\ intp n = dimensions[0];\ intp i;\ for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op += os) #define BINARY_LOOP_TWO_OUT\ char *ip1 = args[0], *ip2 = args[1], *op1 = args[2], *op2 = args[3];\ intp is1 = steps[0], is2 = steps[1], os1 = steps[2], os2 = steps[3];\ intp n = dimensions[0];\ intp i;\ for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1, op2 += os2) /* ***************************************************************************** ** 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 /**begin repeat * #kind = equal, not_equal, greater, greater_equal, less, less_equal, * logical_and, logical_or# * #OP = ==, !=, >, >=, <, <=, &&, ||# **/ static void BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { Bool in1 = *((Bool *)ip1) != 0; Bool in2 = *((Bool *)ip2) != 0; *((Bool *)op)= in1 @OP@ in2; } } /**end repeat**/ static void BOOL_logical_xor(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { Bool in1 = *((Bool *)ip1) != 0; Bool in2 = *((Bool *)ip2) != 0; *((Bool *)op)= (in1 && !in2) || (!in1 && in2); } } /**begin repeat * #kind = maximum, minimum# * #OP = >, <# **/ static void BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { Bool in1 = *((Bool *)ip1) != 0; Bool in2 = *((Bool *)ip2) != 0; *((Bool *)op) = (in1 @OP@ in2) ? in1 : in2; } } /**end repeat**/ /**begin repeat * #kind = absolute, logical_not# * #OP = !=, ==# **/ static void BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { Bool in1 = *(Bool *)ip1; *((Bool *)op) = in1 @OP@ 0; } } /**end repeat**/ static void BOOL_ones_like(char **args, intp *dimensions, intp *steps, void *data) { OUTPUT_LOOP { *((Bool *)op) = 1; } } /* ***************************************************************************** ** INTEGER LOOPS ***************************************************************************** */ /**begin repeat * #type = byte, short, int, long, longlong# * #TYPE = BYTE, SHORT, INT, LONG, LONGLONG# * #ftype = float, float, double, double, double# */ /**begin repeat1 * both signed and unsigned integer types * #s = , u# * #S = , U# */ #define @S@@TYPE@_floor_divide @S@@TYPE@_divide static void @S@@TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *data) { OUTPUT_LOOP { *((@s@@type@ *)op) = 1; } } static void @S@@TYPE@_square(char **args, intp *dimensions, intp *steps, void *data) { UNARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; *((@s@@type@ *)op) = in1*in1; } } static void @S@@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data) { UNARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; *((@s@@type@ *)op) = 1.0/in1; } } static void @S@@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; *((@s@@type@ *)op) = in1; } } static void @S@@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; *((@s@@type@ *)op) = (@s@@type@)(-(@type@)in1); } } static void @S@@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; *((Bool *)op) = !in1; } } static void @S@@TYPE@_invert(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; *((@s@@type@ *)op) = ~in1; } } /**begin repeat2 * Arithmetic * #kind = add, subtract, multiply, bitwise_and, bitwise_or, bitwise_xor, * left_shift, right_shift# * #OP = +, -,*, &, |, ^, <<, >># */ static void @S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; const @s@@type@ in2 = *(@s@@type@ *)ip2; *((@s@@type@ *)op) = in1 @OP@ in2; } } /**end repeat2**/ /**begin repeat2 * #kind = equal, not_equal, greater, greater_equal, less, less_equal, * logical_and, logical_or# * #OP = ==, !=, >, >=, <, <=, &&, ||# */ static void @S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; const @s@@type@ in2 = *(@s@@type@ *)ip2; *((Bool *)op) = in1 @OP@ in2; } } /**end repeat2**/ static void @S@@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; const @s@@type@ in2 = *(@s@@type@ *)ip2; *((Bool *)op)= (in1 && !in2) || (!in1 && in2); } } /**begin repeat2 * #kind = maximum, minimum# * #OP = >, <# **/ static void @S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; const @s@@type@ in2 = *(@s@@type@ *)ip2; *((@s@@type@ *)op) = (in1 @OP@ in2) ? in1 : in2; } } /**end repeat2**/ static void @S@@TYPE@_true_divide(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; const @s@@type@ in2 = *(@s@@type@ *)ip2; if (in2 == 0) { generate_divbyzero_error(); *((@ftype@ *)op) = 0; } else { *((@ftype@ *)op) = (@ftype@)in1 / (@ftype@)in2; } } } static void @S@@TYPE@_power(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @ftype@ in1 = (@ftype@)*(@s@@type@ *)ip1; const @ftype@ in2 = (@ftype@)*(@s@@type@ *)ip2; *((@s@@type@ *)op) = (@s@@type@) pow(in1, in2); } } static void @S@@TYPE@_fmod(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @s@@type@ in1 = *(@s@@type@ *)ip1; const @s@@type@ in2 = *(@s@@type@ *)ip2; if (in2 == 0) { generate_divbyzero_error(); *((@s@@type@ *)op) = 0; } else { *((@s@@type@ *)op)= in1 % in2; } } } /**end repeat1**/ static void U@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const u@type@ in1 = *(u@type@ *)ip1; *((u@type@ *)op) = in1; } } static void @TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((@type@ *)op) = (in1 >= 0) ? in1 : -in1; } } static void U@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const u@type@ in1 = *(u@type@ *)ip1; *((u@type@ *)op) = in1 > 0 ? 1 : 0; } } static void @TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((@type@ *)op) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0); } } static void @TYPE@_divide(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; if (in2 == 0) { generate_divbyzero_error(); *((@type@ *)op) = 0; } else if (((in1 > 0) != (in2 > 0)) && (in1 % in2 != 0)) { *((@type@ *)op) = in1/in2 - 1; } else { *((@type@ *)op) = in1/in2; } } } static void U@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const u@type@ in1 = *(u@type@ *)ip1; const u@type@ in2 = *(u@type@ *)ip2; if (in2 == 0) { generate_divbyzero_error(); *((u@type@ *)op) = 0; } else { *((u@type@ *)op)= in1/in2; } } } static void @TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; if (in2 == 0) { generate_divbyzero_error(); *((@type@ *)op) = 0; } else { /* handle mixed case the way Python does */ const @type@ rem = in1 % in2; if ((in1 > 0) == (in2 > 0) || rem == 0) { *((@type@ *)op) = rem; } else { *((@type@ *)op) = rem + in2; } } } } static void U@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const u@type@ in1 = *(u@type@ *)ip1; const u@type@ in2 = *(u@type@ *)ip2; if (in2 == 0) { generate_divbyzero_error(); *((@type@ *)op) = 0; } else { *((@type@ *)op) = in1 % in2; } } } /**end repeat**/ /* ***************************************************************************** ** FLOAT LOOPS ** ***************************************************************************** */ /**begin repeat * Float types * #type = float, double, longdouble# * #TYPE = FLOAT, DOUBLE, LONGDOUBLE# * #c = f, , l# * #C = F, , L# */ /**begin repeat1 * Arithmetic * # kind = add, subtract, multiply, divide# * # OP = +, -, *, /# */ static void @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; *((@type@ *)op) = in1 @OP@ in2; } } /**end repeat1**/ /**begin repeat1 * #kind = equal, not_equal, less, less_equal, greater, greater_equal, * logical_and, logical_or# * #OP = ==, !=, <, <=, >, >=, &&, ||# */ static void @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; *((Bool *)op) = in1 @OP@ in2; } } /**end repeat1**/ static void @TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; *((Bool *)op)= (in1 && !in2) || (!in1 && in2); } } static void @TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((Bool *)op) = !in1; } } /**begin repeat1 * #kind = isnan, isinf, isfinite, signbit# * #func = isnan, isinf, isfinite, signbit# **/ static void @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((Bool *)op) = @func@(in1) != 0; } } /**end repeat1**/ /**begin repeat1 * #kind = maximum, minimum# * #OP = >, <# **/ static void @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { /* */ BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; *((@type@ *)op) = in1 @OP@ in2 ? in1 : in2; } } /**end repeat1**/ static void @TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; *((@type@ *)op) = floor@c@(in1/in2); } } static void @TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ in2 = *(@type@ *)ip2; const @type@ res = fmod@c@(in1,in2); if (res && ((in2 < 0) != (res < 0))) { *((@type@ *)op) = res + in2; } else { *((@type@ *)op) = res; } } } static void @TYPE@_square(char **args, intp *dimensions, intp *steps, void *data) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((@type@ *)op) = in1*in1; } } static void @TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((@type@ *)op) = 1.0/in1; } } static void @TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *data) { OUTPUT_LOOP { *((@type@ *)op) = 1; } } static void @TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((@type@ *)op) = in1; } } static void @TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const @type@ tmp = (in1 > 0) ? in1 : -in1; /* add 0 to clear -0.0 */ *((@type@ *)op) = tmp + 0; } } static void @TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; *((@type@ *)op) = -in1; } } static void @TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func) { /* */ UNARY_LOOP { const @type@ in1 = *(@type@ *)ip1; if (in1 > 0) { *((@type@ *)op) = 1; } else if (in1 < 0) { *((@type@ *)op) = -1; } else { *((@type@ *)op) = 0; } } } static void @TYPE@_modf(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP_TWO_OUT { const @type@ in1 = *(@type@ *)ip1; *(@type@ *)op1 = modf@c@(in1, (@type@ *)op2); } } #ifdef HAVE_FREXP@C@ static void @TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP_TWO_OUT { const @type@ in1 = *(@type@ *)ip1; *(@type@ *)op1 = frexp@c@(in1, (int *)op2); } } #endif #ifdef HAVE_LDEXP@C@ static void @TYPE@_ldexp(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1 = *(@type@ *)ip1; const int in2 = *(int *)ip2; *(@type@ *)op = ldexp@c@(in1, in2); } } #endif #define @TYPE@_true_divide @TYPE@_divide /**end repeat**/ /* ***************************************************************************** ** COMPLEX LOOPS ** ***************************************************************************** */ /**begin repeat * complex types * #ctype= cfloat, cdouble, clongdouble# * #CTYPE= CFLOAT, CDOUBLE, CLONGDOUBLE# * #type = float, double, longdouble# * #c = f, , l# */ /**begin repeat1 * arithmetic * #kind = add, subtract# * #OP = +, -# */ static void @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; ((@type@ *)op)[0] = in1r @OP@ in2r; ((@type@ *)op)[1] = in1i @OP@ in2i; } } /**end repeat1**/ static void @CTYPE@_multiply(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; ((@type@ *)op)[0] = in1r*in2r - in1i*in2i; ((@type@ *)op)[1] = in1r*in2i + in1i*in2r; } } static void @CTYPE@_divide(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; @type@ d = in2r*in2r + in2i*in2i; ((@type@ *)op)[0] = (in1r*in2r + in1i*in2i)/d; ((@type@ *)op)[1] = (in1i*in2r - in1r*in2i)/d; } } static void @CTYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; @type@ d = in2r*in2r + in2i*in2i; ((@type@ *)op)[0] = floor@c@((in1r*in2r + in1i*in2i)/d); ((@type@ *)op)[1] = 0; } } /**begin repeat1 #kind = equal, not_equal# #OP1 = ==, !=# #OP2 = &&, ||# */ static void @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; *((Bool *)op) = (in1r @OP1@ in2r) @OP2@ (in1i @OP1@ in2i); } } /**end repeat1**/ /**begin repeat1 * #kind= greater, greater_equal, less, less_equal# * #OP = >, >=, <, <=# */ static void @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; if (in1r != in2r) { *((Bool *)op) = in1r @OP@ in2r ? 1 : 0; } else { *((Bool *)op) = in1i @OP@ in2i ? 1 : 0; } } } /**end repeat1**/ /**begin repeat1 #kind = logical_and, logical_or# #OP1 = ||, ||# #OP2 = &&, ||# */ static void @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; *((Bool *)op) = (in1r @OP1@ in1i) @OP2@ (in2r @OP1@ in2i); } } /**end repeat1**/ static void @CTYPE@_logical_xor(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]; const Bool tmp1 = (in1r || in1i); const Bool tmp2 = (in2r || in2i); *((Bool *)op) = (tmp1 && !tmp2) || (!tmp1 && tmp2); } } static void @CTYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; *((Bool *)op) = !(in1r || in1i); } } /**begin repeat1 * #kind = isnan, isinf, isfinite# * #func = isnan, isinf, isfinite# * #OP = ||, ||, &&# **/ static void @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; *((Bool *)op) = @func@(in1r) @OP@ @func@(in1i); } } /**end repeat1**/ static void @CTYPE@_square(char **args, intp *dimensions, intp *steps, void *data) { UNARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; ((@type@ *)op)[0] = in1r*in1r - in1i*in1i; ((@type@ *)op)[1] = in1r*in1i + in1i*in1r; } } static void @CTYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data) { UNARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; if (fabs@c@(in1i) <= fabs@c@(in1r)) { const @type@ r = in1i/in1r; const @type@ d = in1r + in1i*r; ((@type@ *)op)[0] = 1/d; ((@type@ *)op)[1] = -r/d; } else { const @type@ r = in1r/in1i; const @type@ d = in1r*r + in1i; ((@type@ *)op)[0] = r/d; ((@type@ *)op)[1] = -1/d; } } } static void @CTYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *data) { OUTPUT_LOOP { ((@type@ *)op)[0] = 1; ((@type@ *)op)[1] = 0; } } static void @CTYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; ((@type@ *)op)[0] = in1r; ((@type@ *)op)[1] = -in1i; } } static void @CTYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; *((@type@ *)op) = sqrt@c@(in1r*in1r + in1i*in1i); } } static void @CTYPE@_sign(char **args, intp *dimensions, intp *steps, void *func) { UNARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; if (in1r > 0) { ((@type@ *)op)[0] = 1; } else if (in1r < 0) { ((@type@ *)op)[0] = -1; } else { if (in1i > 0) { ((@type@ *)op)[0] = 1; } else if (in1i < 0) { ((@type@ *)op)[0] = -1; } else { ((@type@ *)op)[0] = 0; } } ((@type@ *)op)[1] = 0; } } /**begin repeat1 * #kind = maximum, minimum# * #OP = >, <# */ static void @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { const @type@ in1r = ((@type@ *)ip1)[0]; const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; if (in1r @OP@ in2r || ((in1r == in2r) && (in1i @OP@ in2i))) { ((@type@ *)op)[0] = in1r; ((@type@ *)op)[1] = in1i; } else { ((@type@ *)op)[0] = in2r; ((@type@ *)op)[1] = in2i; } } } /**end repeat1**/ #define @CTYPE@_true_divide @CTYPE@_divide /**end repeat**/ /* ***************************************************************************** ** OBJECT LOOPS ** ***************************************************************************** */ /**begin repeat * #kind = equal, not_equal, greater, greater_equal, less, less_equal# * #OP = EQ, NE, GT, GE, LT, LE# */ static void OBJECT_@kind@(char **args, intp *dimensions, intp *steps, void *func) { BINARY_LOOP { PyObject *in1 = *(PyObject **)ip1; PyObject *in2 = *(PyObject **)ip2; *(Bool *)op = (Bool) PyObject_RichCompareBool(in1, in2, Py_@OP@); } } /**end repeat**/ OBJECT_sign(char **args, intp *dimensions, intp *steps, void *func) { PyObject *zero = PyInt_FromLong(0); UNARY_LOOP { PyObject *in1 = *(PyObject **)ip1; *(PyObject **)op = PyInt_FromLong(PyObject_Compare(in1, zero)); } Py_DECREF(zero); } /* ***************************************************************************** ** END LOOPS ** ***************************************************************************** */ static PyUFuncGenericFunction frexp_functions[] = { #ifdef HAVE_FREXPF FLOAT_frexp, #endif DOUBLE_frexp #ifdef HAVE_FREXPL ,LONGDOUBLE_frexp #endif }; static void * blank3_data[] = { (void *)NULL, (void *)NULL, (void *)NULL}; static char frexp_signatures[] = { #ifdef HAVE_FREXPF PyArray_FLOAT, PyArray_FLOAT, PyArray_INT, #endif PyArray_DOUBLE, PyArray_DOUBLE, PyArray_INT #ifdef HAVE_FREXPL ,PyArray_LONGDOUBLE, PyArray_LONGDOUBLE, PyArray_INT #endif }; static PyUFuncGenericFunction ldexp_functions[] = { #ifdef HAVE_LDEXPF FLOAT_ldexp, #endif DOUBLE_ldexp #ifdef HAVE_LDEXPL ,LONGDOUBLE_ldexp #endif }; static char ldexp_signatures[] = { #ifdef HAVE_LDEXPF PyArray_FLOAT, PyArray_INT, PyArray_FLOAT, #endif PyArray_DOUBLE, PyArray_INT, PyArray_DOUBLE #ifdef HAVE_LDEXPL ,PyArray_LONGDOUBLE, PyArray_INT, PyArray_LONGDOUBLE #endif }; #include "__umath_generated.c" #include "ufuncobject.c" #include "__ufunc_api.c" 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; int num=1; #ifdef HAVE_FREXPL num += 1; #endif #ifdef HAVE_FREXPF num += 1; #endif f = PyUFunc_FromFuncAndData(frexp_functions, blank3_data, frexp_signatures, num, 1, 2, PyUFunc_None, "frexp", "Split the number, x, into a normalized"\ " fraction (y1) and exponent (y2)",0); PyDict_SetItemString(dictionary, "frexp", f); Py_DECREF(f); num = 1; #ifdef HAVE_LDEXPL num += 1; #endif #ifdef HAVE_LDEXPF num += 1; #endif f = PyUFunc_FromFuncAndData(ldexp_functions, blank3_data, ldexp_signatures, num, 2, 1, PyUFunc_None, "ldexp", "Compute y = x1 * 2**x2.",0); PyDict_SetItemString(dictionary, "ldexp", f); Py_DECREF(f); return; } static struct PyMethodDef methods[] = { {"frompyfunc", (PyCFunction) ufunc_frompyfunc, METH_VARARGS | METH_KEYWORDS, doc_frompyfunc}, {"seterrobj", (PyCFunction) ufunc_seterr, METH_VARARGS, NULL}, {"geterrobj", (PyCFunction) ufunc_geterr, METH_VARARGS, NULL}, {NULL, NULL, 0} /* sentinel */ }; PyMODINIT_FUNC initumath(void) { PyObject *m, *d, *s, *s2, *c_api; double pinf, pzero, mynan; int UFUNC_FLOATING_POINT_SUPPORT = 1; #ifdef NO_UFUNC_FLOATING_POINT_SUPPORT UFUNC_FLOATING_POINT_SUPPORT = 0; #endif /* Create the module and add the functions */ m = Py_InitModule("umath", methods); /* Import the array */ if (_import_array() < 0) { if (!PyErr_Occurred()) { PyErr_SetString(PyExc_ImportError, "umath failed: Could not import array core."); } return; } /* Initialize the types */ if (PyType_Ready(&PyUFunc_Type) < 0) return; /* Add some symbolic constants to the module */ d = PyModule_GetDict(m); c_api = PyCObject_FromVoidPtr((void *)PyUFunc_API, NULL); if (PyErr_Occurred()) goto err; PyDict_SetItemString(d, "_UFUNC_API", c_api); Py_DECREF(c_api); if (PyErr_Occurred()) goto err; s = PyString_FromString("0.4.0"); PyDict_SetItemString(d, "__version__", s); Py_DECREF(s); /* Load the ufunc operators into the array module's namespace */ InitOperators(d); InitOtherOperators(d); PyDict_SetItemString(d, "pi", s = PyFloat_FromDouble(M_PI)); Py_DECREF(s); PyDict_SetItemString(d, "e", s = PyFloat_FromDouble(exp(1.0))); Py_DECREF(s); #define ADDCONST(str) PyModule_AddIntConstant(m, #str, UFUNC_##str) #define ADDSCONST(str) PyModule_AddStringConstant(m, "UFUNC_" #str, UFUNC_##str) ADDCONST(ERR_IGNORE); ADDCONST(ERR_WARN); ADDCONST(ERR_CALL); ADDCONST(ERR_RAISE); ADDCONST(ERR_PRINT); ADDCONST(ERR_LOG); ADDCONST(ERR_DEFAULT); ADDCONST(ERR_DEFAULT2); ADDCONST(SHIFT_DIVIDEBYZERO); ADDCONST(SHIFT_OVERFLOW); ADDCONST(SHIFT_UNDERFLOW); ADDCONST(SHIFT_INVALID); ADDCONST(FPE_DIVIDEBYZERO); ADDCONST(FPE_OVERFLOW); ADDCONST(FPE_UNDERFLOW); ADDCONST(FPE_INVALID); ADDCONST(FLOATING_POINT_SUPPORT); ADDSCONST(PYVALS_NAME); #undef ADDCONST #undef ADDSCONST PyModule_AddIntConstant(m, "UFUNC_BUFSIZE_DEFAULT", (long)PyArray_BUFSIZE); pinf = pinf_init(); pzero = pzero_init(); mynan = pinf / pinf; PyModule_AddObject(m, "PINF", PyFloat_FromDouble(pinf)); PyModule_AddObject(m, "NINF", PyFloat_FromDouble(-pinf)); PyModule_AddObject(m, "PZERO", PyFloat_FromDouble(pzero)); PyModule_AddObject(m, "NZERO", PyFloat_FromDouble(-pzero)); PyModule_AddObject(m, "NAN", PyFloat_FromDouble(mynan)); s = PyDict_GetItemString(d, "conjugate"); s2 = PyDict_GetItemString(d, "remainder"); /* Setup the array object's numerical structures with appropriate ufuncs in d*/ PyArray_SetNumericOps(d); PyDict_SetItemString(d, "conj", s); PyDict_SetItemString(d, "mod", s2); return; err: /* Check for errors */ if (!PyErr_Occurred()) { PyErr_SetString(PyExc_RuntimeError, "cannot load umath module."); } return; }