From dd96fbda20855997eb6d234a345d45c462c4e74b Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Wed, 1 Oct 2008 18:08:41 +0000 Subject: Cleanup ufunc loops. At this point loops are separated into variable kinds, so there is a fair amount of duplication. I will probably merge loops that look the same in a later commit. There are no changes to current behavior of loops, this will also be changed in later work to deal with nans and such. --- numpy/core/src/umathmodule.c.src | 1680 ++++++++++++++++++-------------------- 1 file changed, 790 insertions(+), 890 deletions(-) (limited to 'numpy/core/src/umathmodule.c.src') diff --git a/numpy/core/src/umathmodule.c.src b/numpy/core/src/umathmodule.c.src index 77675c441..6f41cd4fa 100644 --- a/numpy/core/src/umathmodule.c.src +++ b/numpy/core/src/umathmodule.c.src @@ -1,5 +1,9 @@ /* -*- c -*- */ +/* + * vim:syntax=c + */ + /* ***************************************************************************** ** INCLUDES ** @@ -13,14 +17,40 @@ #include "config.h" #include +#ifndef M_PI +#define M_PI 3.14159265358979323846264338328 +#endif + /* ***************************************************************************** ** BASIC MATH FUNCTIONS ** ***************************************************************************** */ -/* A whole slew of basic math functions are provided originally - by Konrad Hinsen. */ +float degreesf(float x) { + return x * (float)(180.0/M_PI); +} +double degrees(double x) { + return x * (180.0/M_PI); +} +longdouble degreesl(longdouble x) { + return x * (180.0L/M_PI); +} + +float radiansf(float x) { + return x * (float)(M_PI/180.0); +} +double radians(double x) { + return x * (M_PI/180.0); +} +longdouble radiansl(longdouble x) { + return x * (M_PI/180.0L); +} + +/* + * A whole slew of basic math functions are provided originally + * by Konrad Hinsen. + */ #if !defined(__STDC__) && !defined(_MSC_VER) extern double fmod (double, double); @@ -28,9 +58,6 @@ extern double frexp (double, int *); extern double ldexp (double, int); extern double modf (double, double *); #endif -#ifndef M_PI -#define M_PI 3.14159265358979323846264338328 -#endif #if defined(DISTUTILS_USE_SDK) @@ -443,26 +470,6 @@ trunc(double x) #define isfinitef(x) (!(isinff((x)) || isnanf((x)))) #define isfinitel(x) (!(isinfl((x)) || isnanl((x)))) -float degreesf(float x) { - return x * (float)(180.0/M_PI); -} -double degrees(double x) { - return x * (180.0/M_PI); -} -longdouble degreesl(longdouble x) { - return x * (180.0L/M_PI); -} - -float radiansf(float x) { - return x * (float)(M_PI/180.0); -} -double radians(double x) { - return x * (M_PI/180.0); -} -longdouble radiansl(longdouble x) { - return x * (M_PI/180.0L); -} - /* First, the C functions that do the real work */ /* if C99 extensions not available then define dummy functions that use the @@ -624,6 +631,61 @@ static float expm1f(float x) } #endif +/* + ***************************************************************************** + ** PYTHON OBJECT FUNCTIONS ** + ***************************************************************************** + */ + +static PyObject * +Py_square(PyObject *o) +{ + return PyNumber_Multiply(o, o); +} + +static PyObject * +Py_get_one(PyObject *o) +{ + return PyInt_FromLong(1); +} + +static PyObject * +Py_reciprocal(PyObject *o) +{ + PyObject *one = PyInt_FromLong(1); + PyObject *result; + + if (!one) { + return NULL; + } + result = PyNumber_Divide(one, o); + Py_DECREF(one); + return result; +} + +/**begin repeat + * #Kind = Max, Min# + * #OP = >=, <=# + */ +static PyObject * +_npy_Object@Kind@(PyObject *i1, PyObject *i2) +{ + PyObject *result; + int cmp; + + if (PyObject_Cmp(i1, i2, &cmp) < 0) { + return NULL; + } + if (cmp @OP@ 0) { + result = i1; + } + else { + result = i2; + } + Py_INCREF(result); + return result; +} +/**end repeat**/ /* ***************************************************************************** @@ -632,11 +694,12 @@ static float expm1f(float x) */ -/* Don't pass structures between functions (only pointers) because how - structures are passed is compiler dependent and could cause - segfaults if ufuncobject.c is compiled with a different compiler - than an extension that makes use of the UFUNC API -*/ +/* + * Don't pass structures between functions (only pointers) because how + * structures are passed is compiler dependent and could cause + * segfaults if ufuncobject.c is compiled with a different compiler + * than an extension that makes use of the UFUNC API + */ /**begin repeat @@ -650,9 +713,10 @@ static c@typ@ nc_half@c@ = {0.5, 0.}; static c@typ@ nc_i@c@ = {0., 1.}; static c@typ@ nc_i2@c@ = {0., 0.5}; /* - static c@typ@ nc_mi@c@ = {0., -1.}; - static c@typ@ nc_pi2@c@ = {M_PI/2., 0.}; -*/ + * static c@typ@ nc_mi@c@ = {0., -1.}; + * static c@typ@ nc_pi2@c@ = {M_PI/2., 0.}; + */ + static void nc_sum@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) @@ -681,7 +745,7 @@ nc_neg@c@(c@typ@ *a, c@typ@ *r) static void nc_prod@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) { - register @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; + @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; r->real = ar*br - ai*bi; r->imag = ar*bi + ai*br; return; @@ -691,8 +755,8 @@ static void nc_quot@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) { - register @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; - register @typ@ d = br*br + bi*bi; + @typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag; + @typ@ d = br*br + bi*bi; r->real = (ar*br + ai*bi)/d; r->imag = (ai*br - ar*bi)/d; return; @@ -801,8 +865,10 @@ nc_pow@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r) return; } } - /* complexobect.c uses an inline version of this formula - investigate whether this had better performance or accuracy */ + /* + * complexobect.c uses an inline version of this formula + * investigate whether this had better performance or accuracy + */ nc_log@c@(a, r); nc_prod@c@(r, b, r); nc_exp@c@(r, r); @@ -823,6 +889,10 @@ nc_prodi@c@(c@typ@ *x, c@typ@ *r) static void nc_acos@c@(c@typ@ *x, c@typ@ *r) { + /* + * return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i, + * nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))))); + */ nc_prod@c@(x,x,r); nc_diff@c@(&nc_1@c@, r, r); nc_sqrt@c@(r, r); @@ -832,14 +902,15 @@ nc_acos@c@(c@typ@ *x, c@typ@ *r) nc_prodi@c@(r, r); nc_neg@c@(r, r); return; - /* return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i, - nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))))); - */ } static void nc_acosh@c@(c@typ@ *x, c@typ@ *r) { + /* + * return nc_log(nc_sum(x, + * nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1))))); + */ c@typ@ t; nc_sum@c@(x, &nc_1@c@, &t); @@ -850,15 +921,15 @@ nc_acosh@c@(c@typ@ *x, c@typ@ *r) nc_sum@c@(x, r, r); nc_log@c@(r, r); return; - /* - return nc_log(nc_sum(x, - nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1))))); - */ } static void nc_asin@c@(c@typ@ *x, c@typ@ *r) { + /* + * return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x), + * nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))); + */ c@typ@ a, *pa=&a; nc_prod@c@(x, x, r); nc_diff@c@(&nc_1@c@, r, r); @@ -869,30 +940,29 @@ nc_asin@c@(c@typ@ *x, c@typ@ *r) nc_prodi@c@(r, r); nc_neg@c@(r, r); return; - /* - return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x), - nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))); - */ } static void nc_asinh@c@(c@typ@ *x, c@typ@ *r) { + /* + * return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x)); + */ nc_prod@c@(x, x, r); nc_sum@c@(&nc_1@c@, r, r); nc_sqrt@c@(r, r); nc_sum@c@(r, x, r); nc_log@c@(r, r); return; - /* - return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x)); - */ } static void nc_atan@c@(c@typ@ *x, c@typ@ *r) { + /* + * return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x)))); + */ c@typ@ a, *pa=&a; nc_diff@c@(&nc_i@c@, x, pa); nc_sum@c@(&nc_i@c@, x, r); @@ -900,14 +970,14 @@ nc_atan@c@(c@typ@ *x, c@typ@ *r) nc_log@c@(r,r); nc_prod@c@(&nc_i2@c@, r, r); return; - /* - return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x)))); - */ } static void nc_atanh@c@(c@typ@ *x, c@typ@ *r) { + /* + * return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x)))); + */ c@typ@ a, *pa=&a; nc_diff@c@(&nc_1@c@, x, r); nc_sum@c@(&nc_1@c@, x, pa); @@ -915,9 +985,6 @@ nc_atanh@c@(c@typ@ *x, c@typ@ *r) nc_log@c@(r, r); nc_prod@c@(&nc_half@c@, r, r); return; - /* - return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x)))); - */ } static void @@ -1018,1137 +1085,973 @@ nc_tanh@c@(c@typ@ *x, c@typ@ *r) ***************************************************************************** */ +#define OUTPUT_LOOP\ + char *op = args[1];\ + intp os = steps[1];\ + intp n = dimensions[0];\ + intp i;\ + for(i = 0; i < n; i++, op += os) + +#define UNARY_LOOP\ + char *ip1 = args[0], *op = args[1];\ + intp is1 = steps[0], os = steps[1];\ + intp n = dimensions[0];\ + intp i;\ + for(i = 0; i < n; i++, ip1 += is1, op += os) + +#define UNARY_LOOP_TWO_OUT\ + char *ip1 = args[0], *op1 = args[1], *op2 = args[2];\ + intp is1 = steps[0], os1 = steps[1], os2 = steps[2];\ + intp n = dimensions[0];\ + intp i;\ + for(i = 0; i < n; i++, ip1 += is1, op1 += os1, op2 += os2) + +#define BINARY_LOOP\ + char *ip1 = args[0], *ip2 = args[1], *op = args[2];\ + intp is1 = steps[0], is2 = steps[1], os = steps[2];\ + intp n = dimensions[0];\ + intp i;\ + for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op += os) + +#define BINARY_LOOP_TWO_OUT\ + char *ip1 = args[0], *ip2 = args[1], *op1 = args[2], *op2 = args[3];\ + intp is1 = steps[0], is2 = steps[1], os1 = steps[2], os2 = steps[3];\ + intp n = dimensions[0];\ + intp i;\ + for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1, op2 += os2) + +/* + ***************************************************************************** + ** 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 = ==, !=, >, >=, <, <=, &&, ||# + **/ - #TYPE=(BOOL, BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*2# - #OP=||, +*13, ^, -*13# - #kind=add*14, subtract*14# - #typ=(Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*2# -*/ +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 -@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +BOOL_logical_xor(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i, <# + **/ +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 - - #TYPE=(CFLOAT, CDOUBLE, CLONGDOUBLE)*2# - #OP=+*3,-*3# - #kind=add*3,subtract*3# - #typ=(float, double, longdouble)*2# -*/ - + * #kind = absolute, logical_not# + * #OP = !=, ==# + **/ static void -@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; ireal, \ - ai=((c@typ@ *)i1)->imag, \ - br=((c@typ@ *)i2)->real, \ - bi=((c@typ@ *)i2)->imag; - ((c@typ@ *)op)->real = ar*br - ai*bi; - ((c@typ@ *)op)->imag = ar*bi + ai*br; + UNARY_LOOP { + const @s@@type@ in1 = *(@s@@type@ *)ip1; + *((@s@@type@ *)op) = in1*in1; } } static void -@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func) +@S@@TYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data) { - register intp i; - intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for (i=0; ireal, \ - ai=((c@typ@ *)i1)->imag, \ - br=((c@typ@ *)i2)->real, \ - bi=((c@typ@ *)i2)->imag; - register @typ@ d = br*br + bi*bi; - ((c@typ@ *)op)->real = (ar*br + ai*bi)/d; - ((c@typ@ *)op)->imag = (ai*br - ar*bi)/d; + UNARY_LOOP { + const @s@@type@ in1 = *(@s@@type@ *)ip1; + *((@s@@type@ *)op) = 1.0/in1; } } static void -@TYP@_floor_divide(char **args, intp *dimensions, intp *steps, void *func) +@S@@TYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; ireal, \ - ai=((c@typ@ *)i1)->imag, \ - br=((c@typ@ *)i2)->real, \ - bi=((c@typ@ *)i2)->imag; - register @typ@ d = br*br + bi*bi; - ((c@typ@ *)op)->real = floor@c@((ar*br + ai*bi)/d); - ((c@typ@ *)op)->imag = 0; + UNARY_LOOP { + const @s@@type@ in1 = *(@s@@type@ *)ip1; + *((@s@@type@ *)op) = in1; } } -#define @TYP@_true_divide @TYP@_divide -/**end repeat**/ - - -/**begin repeat - #TYP=UBYTE,USHORT,UINT,ULONG,ULONGLONG# - #typ=ubyte, ushort, uint, ulong, ulonglong# -*/ static void -@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func) +@S@@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func) { - register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i 0) != (y > 0)) && (x % y != 0)) tmp--; - *((@typ@ *)op)= tmp; - } + UNARY_LOOP { + const @s@@type@ in1 = *(@s@@type@ *)ip1; + *((Bool *)op) = !in1; } } -/**end repeat**/ - -/**begin repeat - #TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG# - #typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong# - #otyp=float*4, double*6# -*/ static void -@TYP@_true_divide(char **args, intp *dimensions, intp *steps, void *func) +@S@@TYPE@_invert(char **args, intp *dimensions, intp *steps, void *func) { - register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i># + */ static void -@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func) +@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i, >=, <, <=, &&, ||# + */ static void -@TYP@_floor_divide(char **args, intp *dimensions, intp *steps, void *func) +@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i, <# + **/ static void -@TYP@_square(char **args, intp *dimensions, intp *steps, void *data) +@S@@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - intp i, is1 = steps[0], os = steps[1], n = dimensions[0]; - char *i1 = args[0], *op = args[1]; - c@typ@ *x, *y; - @typ@ xr, xi; - - for (i = 0; i < n; i++, i1 += is1, op += os) { - x = (c@typ@ *)i1; - y = (c@typ@ *)op; - xr = x->real; - xi = x->imag; - y->real = xr*xr - xi*xi; - y->imag = 2*xr*xi; + BINARY_LOOP { + const @s@@type@ in1 = *(@s@@type@ *)ip1; + const @s@@type@ in2 = *(@s@@type@ *)ip2; + *((@s@@type@ *)op) = (in1 @OP@ in2) ? in1 : in2; } } -/**end repeat**/ +/**end repeat2**/ -static PyObject * -Py_square(PyObject *o) +static void +@S@@TYPE@_true_divide(char **args, intp *dimensions, intp *steps, void *func) { - return PyNumber_Multiply(o, o); + BINARY_LOOP { + const @s@@type@ in1 = *(@s@@type@ *)ip1; + const @s@@type@ in2 = *(@s@@type@ *)ip2; + if (in2 == 0) { + generate_divbyzero_error(); + *((@ftype@ *)op) = 0; + } + else { + *((@ftype@ *)op) = (@ftype@)in1 / (@ftype@)in2; + } + } } - -/**begin repeat - #TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE# - #typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble# -*/ static void -@TYP@_reciprocal(char **args, intp *dimensions, intp *steps, void *data) +@S@@TYPE@_power(char **args, intp *dimensions, intp *steps, void *func) { - intp i, is1 = steps[0], os = steps[1], n = dimensions[0]; - char *i1 = args[0], *op = args[1]; - - for (i = 0; i < n; i++, i1 += is1, op += os) { - @typ@ x = *((@typ@ *)i1); - *((@typ@ *)op) = (@typ@) (1.0 / x); + BINARY_LOOP { + const @ftype@ in1 = (@ftype@)*(@s@@type@ *)ip1; + const @ftype@ in2 = (@ftype@)*(@s@@type@ *)ip2; + *((@s@@type@ *)op) = (@s@@type@) pow(in1, in2); } } -/**end repeat**/ -/**begin repeat - #TYP=CFLOAT,CDOUBLE,CLONGDOUBLE# - #typ=float, double, longdouble# -*/ static void -@TYP@_reciprocal(char **args, intp *dimensions, intp *steps, void *data) -{ - intp i, is1 = steps[0], os = steps[1], n = dimensions[0]; - char *i1 = args[0], *op = args[1]; - c@typ@ *x, *y; - @typ@ xr, xi, r, denom; - - for (i = 0; i < n; i++, i1 += is1, op += os) { - x = (c@typ@ *)i1; - y = (c@typ@ *)op; - xr = x->real; - xi = x->imag; - if (fabs(xi) <= fabs(xr)) { - r = xi / xr; - denom = xr + xi * r; - y->real = 1 / denom; - y->imag = -r / denom; - } else { - r = xr / xi; - denom = xr * r + xi; - y->real = r / denom; - y->imag = -1 / denom; +@S@@TYPE@_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 repeat**/ +/**end repeat1**/ -static PyObject * -Py_reciprocal(PyObject *o) +static void +U@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) { - PyObject *one, *result; - one = PyInt_FromLong(1); - if (!one) return NULL; - result = PyNumber_Divide(one, o); - Py_DECREF(one); - return result; + UNARY_LOOP { + const u@type@ in1 = *(u@type@ *)ip1; + *((u@type@ *)op) = in1; + } } -static PyObject * -_npy_ObjectMax(PyObject *i1, PyObject *i2) +static void +@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) { - int cmp; - PyObject *res; - if (PyObject_Cmp(i1, i2, &cmp) < 0) return NULL; - - if (cmp >= 0) { - res = i1; - } - else { - res = i2; + UNARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + *((@type@ *)op) = (in1 >= 0) ? in1 : -in1; } - Py_INCREF(res); - return res; } -static PyObject * -_npy_ObjectMin(PyObject *i1, PyObject *i2) +static void +U@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func) { - int cmp; - PyObject *res; - if (PyObject_Cmp(i1, i2, &cmp) < 0) return NULL; - - if (cmp <= 0) { - res = i1; - } - else { - res = i2; + UNARY_LOOP { + const u@type@ in1 = *(u@type@ *)ip1; + *((u@type@ *)op) = in1 > 0 ? 1 : 0; } - Py_INCREF(res); - return res; } -/* ones_like is defined here because it's used for x**0 */ - -/**begin repeat - #TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE# - #typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble# -*/ static void -@TYP@_ones_like(char **args, intp *dimensions, intp *steps, void *data) +@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func) { - intp i, os = steps[1], n = dimensions[0]; - char *op = args[1]; - - for (i = 0; i < n; i++, op += os) { - *((@typ@ *)op) = 1; + UNARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + *((@type@ *)op) = in1 > 0 ? 1 : (in1 < 0 ? -1 : 0); } } -/**end repeat**/ -/**begin repeat - #TYP=CFLOAT,CDOUBLE,CLONGDOUBLE# - #typ=float, double, longdouble# -*/ static void -@TYP@_ones_like(char **args, intp *dimensions, intp *steps, void *data) +@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *func) { - intp i, is1 = steps[0], os = steps[1], n = dimensions[0]; - char *i1 = args[0], *op = args[1]; - c@typ@ *y; - - for (i = 0; i < n; i++, i1 += is1, op += os) { - y = (c@typ@ *)op; - y->real = 1.0; - y->imag = 0.0; + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const @type@ in2 = *(@type@ *)ip2; + if (in2 == 0) { + generate_divbyzero_error(); + *((@type@ *)op) = 0; + } + else if (((in1 > 0) != (in2 > 0)) && (in1 % in2 != 0)) { + *((@type@ *)op) = in1/in2 - 1; + } + else { + *((@type@ *)op) = in1/in2; + } } } -/**end repeat**/ -static PyObject * -Py_get_one(PyObject *o) -{ - return PyInt_FromLong(1); -} - - -/**begin repeat - #TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG# - #typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong# - #btyp=float*4, double*6# -*/ static void -@TYP@_power(char **args, intp *dimensions, intp *steps, void *func) +U@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *func) { - register intp i, is1=steps[0],is2=steps[1]; - register intp os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - @btyp@ x, y; - - for(i=0; i 0) == (in2 > 0) || rem == 0) { + *((@type@ *)op) = rem; + } + else { + *((@type@ *)op) = rem + in2; + } + } } } -/**end repeat**/ - -/**begin repeat - #TYP=CFLOAT, CDOUBLE, CLONGDOUBLE# - #typ=float, double, longdouble# -*/ static void -@TYP@_conjugate(char **args, intp *dimensions, intp *steps, void *func) { - register intp i, is1=steps[0], os=steps[1], n=dimensions[0]; - char *i1=args[0], *op=args[1]; - - for(i=0; i, >=, &&, ||# + */ static void -@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) +@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i, n; - intp is1=steps[0], os=steps[1]; - char *i1=args[0], *op=args[1]; - - n=dimensions[0]; - for(i=0; i, >=, <, <=, ==, !=, &&, ||, &, |, ^# -**/ static void -BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func) +@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - Bool in1, in2; - for(i=0; i*13, >=*13, <*13, <=*13# - #typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*4# - #kind= greater*13, greater_equal*13, less*13, less_equal*13# -*/ +/**begin repeat1 + * #kind = isnan, isinf, isfinite, signbit# + * #func = isnan, isinf, isfinite, signbit# + **/ static void @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i*3, >=*3, <*3, <=*3# - #typ=(cfloat, cdouble, clongdouble)*4# - #kind= greater*3, greater_equal*3, less*3, less_equal*3# -*/ +/**end repeat1**/ +/**begin repeat1 + * #kind = maximum, minimum# + * #OP = >, <# + **/ static void @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; ireal == ((@typ@ *)i2)->real) - *((Bool *)op)=((@typ@ *)i1)->imag @OP@ \ - ((@typ@ *)i2)->imag; - else - *((Bool *)op)=((@typ@ *)i1)->real @OP@ \ - ((@typ@ *)i2)->real; + /* */ + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const @type@ in2 = *(@type@ *)ip2; + *((@type@ *)op) = in1 @OP@ in2 ? in1 : in2; } } -/**end repeat**/ +/**end repeat1**/ - -/**begin repeat - #TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*4# - #typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*4# - #OP= ==*13, !=*13, &&*13, ||*13# - #kind=equal*13, not_equal*13, logical_and*13, logical_or*13# -*/ static void -@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +@TYPE@_floor_divide(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i 0) ? in1 : -in1; + /* add 0 to clear -0.0 */ + *((@type@ *)op) = tmp + 0; + } +} -#define _SIGN1(x) ((x) > 0 ? 1 : ((x) < 0 ? -1 : 0)) -#define _SIGN2(x) ((x) == 0 ? 0 : 1) -#define _SIGNC(x) (((x).real > 0) ? 1 : ((x).real < 0 ? -1 : ((x).imag > 0 ? 1 : ((x).imag < 0) ? -1 : 0))) -/**begin repeat - #TYPE=BYTE,SHORT,INT,LONG,LONGLONG,FLOAT,DOUBLE,LONGDOUBLE,UBYTE,USHORT,UINT,ULONG,ULONGLONG# - #typ=byte,short,int,long,longlong,float,double,longdouble,ubyte,ushort,uint,ulong,ulonglong# - #func=_SIGN1*8,_SIGN2*5# -*/ static void -@TYPE@_sign(char **args, intp *dimensions, intp *steps, void *func) +@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],os=steps[1], n=dimensions[0]; - char *i1=args[0], *op=args[1]; - @typ@ t1; - for(i=0; i 0) { + *((@type@ *)op) = 1; + } + else if (in1 < 0) { + *((@type@ *)op) = -1; + } + else { + *((@type@ *)op) = 0; + } } } -/**end repeat**/ - -#undef _SIGN1 -#undef _SIGN2 -#undef _SIGNC - static void -OBJECT_sign(char **args, intp *dimensions, intp *steps, void *func) +@TYPE@_modf(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],os=steps[1], n=dimensions[0]; - char *i1=args[0], *op=args[1]; - PyObject *t1, *zero, *res; - zero = PyInt_FromLong(0); - for(i=0; ireal || \ - ((@typ@ *)i1)->imag); + BINARY_LOOP { + const @type@ in1 = *(@type@ *)ip1; + const int in2 = *(int *)ip2; + *(@type@ *)op = ldexp@c@(in1, in2); } } -/**end repeat**/ +#endif +#undef HAVE_DOUBLE_FUNCS +#define @TYPE@_true_divide @TYPE@_divide +/**end repeat**/ +/* + ***************************************************************************** + ** COMPLEX LOOPS ** + ***************************************************************************** + */ + /**begin repeat - #TYPE=BYTE,SHORT,INT,LONG,LONGLONG# - #typ=byte, short, int, long, longlong# - #c=f*2,,,l*1# -*/ + * complex types + * #ctype= cfloat, cdouble, clongdouble# + * #CTYPE= CFLOAT, CDOUBLE, CLONGDOUBLE# + * #type = float, double, longdouble# + * #c = f, , l# + */ + +/**begin repeat1 + * arithmetic + * #kind = add, subtract# + * #OP = +, -# + */ static void -@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func) +@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - register @typ@ ix,iy, tmp; - for(i=0; i 0) == (iy > 0)) { - *((@typ@ *)op) = ix % iy; - } - else { /* handle mixed case the way Python does */ - tmp = ix % iy; - if (tmp) tmp += iy; - *((@typ@ *)op)= tmp; - } + BINARY_LOOP { + const @type@ in1r = ((@type@ *)ip1)[0]; + const @type@ in1i = ((@type@ *)ip1)[1]; + const @type@ in2r = ((@type@ *)ip2)[0]; + const @type@ in2i = ((@type@ *)ip2)[1]; + ((@type@ *)op)[0] = in1r @OP@ in2r; + ((@type@ *)op)[1] = in1i @OP@ in2i; } } -/**end repeat**/ +/**end repeat1**/ -/**begin repeat - #TYPE=UBYTE,USHORT,UINT,ULONG,ULONGLONG# - #typ=ubyte, ushort, uint, ulong, ulonglong# -*/ static void -@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func) +@CTYPE@_multiply(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - register @typ@ ix,iy; - for(i=0; i, >=, <, <=# + */ +static void +@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +{ + BINARY_LOOP { + const @type@ in1r = ((@type@ *)ip1)[0]; + const @type@ in1i = ((@type@ *)ip1)[1]; + const @type@ in2r = ((@type@ *)ip2)[0]; + const @type@ in2i = ((@type@ *)ip2)[1]; + if (in1r != in2r) { + *((Bool *)op) = in1r @OP@ in2r ? 1 : 0; } else { - *((@typ@ *)op)= x % y; + *((Bool *)op) = in1i @OP@ in2i ? 1 : 0; } - } } -/**end repeat**/ - -/**begin repeat - - #TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG)*5# - #typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong)*5# - #OP= &*10, |*10, ^*10, <<*10, >>*10# - #kind=bitwise_and*10, bitwise_or*10, bitwise_xor*10, left_shift*10, right_shift*10# +/**end repeat1**/ +/**begin repeat1 + #kind = logical_and, logical_or# + #OP1 = ||, ||# + #OP2 = &&, ||# */ static void -@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - register char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; ireal || ((@typ@ *)i1)->imag; - p2 = ((@typ@ *)i2)->real || ((@typ@ *)i2)->imag; - *((Bool *)op)= (p1 || p2) && !(p1 && p2); + UNARY_LOOP { + const @type@ in1r = ((@type@ *)ip1)[0]; + const @type@ in1i = ((@type@ *)ip1)[1]; + ((@type@ *)op)[0] = in1r*in1r - in1i*in1i; + ((@type@ *)op)[1] = in1r*in1i + in1i*in1r; } } -/**end repeat**/ - - -/**begin repeat - - #TYPE=(BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*2# - #OP= >*14, <*14# - #typ=(Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*2# - #kind= maximum*14, minimum*14# -*/ static void -@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +@CTYPE@_reciprocal(char **args, intp *dimensions, intp *steps, void *data) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - for(i=0; i*3, <*3# - #typ=(cfloat, cdouble, clongdouble)*2# - #kind= maximum*3, minimum*3# -*/ static void -@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +@CTYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *data) { - register intp i; - intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0]; - char *i1=args[0], *i2=args[1], *op=args[2]; - @typ@ *i1c, *i2c; - for(i=0; ireal @OP@ i2c->real) || \ - ((i1c->real==i2c->real) && (i1c->imag @OP@ i2c->imag))) - memmove(op, i1, sizeof(@typ@)); - else - memmove(op, i2, sizeof(@typ@)); + OUTPUT_LOOP { + ((@type@ *)op)[0] = 1; + ((@type@ *)op)[1] = 0; } } -/**end repeat**/ - +static void +@CTYPE@_conjugate(char **args, intp *dimensions, intp *steps, void *func) { + UNARY_LOOP { + const @type@ in1r = ((@type@ *)ip1)[0]; + const @type@ in1i = ((@type@ *)ip1)[1]; + ((@type@ *)op)[0] = in1r; + ((@type@ *)op)[1] = -in1i; + } +} -/*** isinf, isinf, isfinite, signbit ***/ -/**begin repeat - #kind=isnan*3, isinf*3, isfinite*3, signbit*3# - #TYPE=(FLOAT, DOUBLE, LONGDOUBLE)*4# - #typ=(float, double, longdouble)*4# - #c=(f,,l)*4# -*/ static void -@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) +@CTYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is=steps[0], os=steps[1], n=dimensions[0]; - char *ip=args[0], *op=args[1]; - for(i=0; i 0) { + ((@type@ *)op)[0] = 1; + } + else if (in1r < 0) { + ((@type@ *)op)[0] = -1; + } + else { + if (in1i > 0) { + ((@type@ *)op)[0] = 1; + } + else if (in1i < 0) { + ((@type@ *)op)[0] = -1; + } + else { + ((@type@ *)op)[0] = 0; + } + } + ((@type@ *)op)[1] = 0; } } -/**end repeat**/ - - - -/****** modf ****/ - -/**begin repeat - #TYPE=FLOAT, DOUBLE, LONGDOUBLE# - #typ=float, double, longdouble# - #c=f,,l# -*/ +/**begin repeat1 + * #kind = maximum, minimum# + * #OP = >, <# + */ static void -@TYPE@_modf(char **args, intp *dimensions, intp *steps, void *func) +@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func) { - register intp i; - intp is1=steps[0],os1=steps[1],os2=steps[2],n=dimensions[0]; - char *i1=args[0], *op1=args[1], *op2=args[2]; - @typ@ x1, y1, y2; - for (i=0; i