/* -*- c -*- */ /* The purpose of this module is to add faster math for array scalars that does not go through the ufunc machinery but still supports error-modes. */ #define _UMATHMODULE #define NPY_NO_DEPRECATED_API NPY_API_VERSION #include "Python.h" #include "npy_config.h" #define PY_ARRAY_UNIQUE_SYMBOL _npy_umathmodule_ARRAY_API #define NO_IMPORT_ARRAY #include "numpy/arrayobject.h" #include "numpy/ufuncobject.h" #include "numpy/arrayscalars.h" #include "npy_pycompat.h" #include "numpy/halffloat.h" #include "templ_common.h" #include "binop_override.h" #include "npy_longdouble.h" /* Basic operations: * * BINARY: * * add, subtract, multiply, divide, remainder, divmod, power, * floor_divide, true_divide * * lshift, rshift, and, or, xor (integers only) * * UNARY: * * negative, positive, absolute, nonzero, invert, int, long, float, oct, hex * */ /**begin repeat * #name = byte, short, int, long, longlong# * #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong# */ static void @name@_ctype_add(@type@ a, @type@ b, @type@ *out) { *out = a + b; if ((*out^a) >= 0 || (*out^b) >= 0) { return; } npy_set_floatstatus_overflow(); return; } static void @name@_ctype_subtract(@type@ a, @type@ b, @type@ *out) { *out = a - b; if ((*out^a) >= 0 || (*out^~b) >= 0) { return; } npy_set_floatstatus_overflow(); return; } /**end repeat**/ /**begin repeat * #name = ubyte, ushort, uint, ulong, ulonglong# * #type = npy_ubyte, npy_ushort, npy_uint, npy_ulong, npy_ulonglong# */ static void @name@_ctype_add(@type@ a, @type@ b, @type@ *out) { *out = a + b; if (*out >= a && *out >= b) { return; } npy_set_floatstatus_overflow(); return; } static void @name@_ctype_subtract(@type@ a, @type@ b, @type@ *out) { *out = a - b; if (a >= b) { return; } npy_set_floatstatus_overflow(); return; } /**end repeat**/ #ifndef NPY_SIZEOF_BYTE #define NPY_SIZEOF_BYTE 1 #endif /**begin repeat * * #name = byte, ubyte, short, ushort, * int, uint, long, ulong# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, * npy_int, npy_uint, npy_long, npy_ulong# * #big = npy_int, npy_uint, npy_int, npy_uint, * npy_longlong, npy_ulonglong, npy_longlong, npy_ulonglong# * #NAME = BYTE, UBYTE, SHORT, USHORT, * INT, UINT, LONG, ULONG# * #SIZENAME = BYTE*2, SHORT*2, INT*2, LONG*2# * #SIZE = INT*4,LONGLONG*4# * #neg = (1,0)*4# */ #if NPY_SIZEOF_@SIZE@ > NPY_SIZEOF_@SIZENAME@ static void @name@_ctype_multiply(@type@ a, @type@ b, @type@ *out) { @big@ temp; temp = ((@big@) a) * ((@big@) b); *out = (@type@) temp; #if @neg@ if (temp > NPY_MAX_@NAME@ || temp < NPY_MIN_@NAME@) #else if (temp > NPY_MAX_@NAME@) #endif npy_set_floatstatus_overflow(); return; } #endif /**end repeat**/ /**begin repeat * * #name = int, uint, long, ulong, * longlong, ulonglong# * #type = npy_int, npy_uint, npy_long, npy_ulong, * npy_longlong, npy_ulonglong# * #SIZE = INT*2, LONG*2, LONGLONG*2# */ #if NPY_SIZEOF_LONGLONG == NPY_SIZEOF_@SIZE@ static void @name@_ctype_multiply(@type@ a, @type@ b, @type@ *out) { if (npy_mul_with_overflow_@name@(out, a, b)) { npy_set_floatstatus_overflow(); } return; } #endif /**end repeat**/ /**begin repeat * * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong# * #neg = (1,0)*5# */ static void @name@_ctype_divide(@type@ a, @type@ b, @type@ *out) { if (b == 0) { npy_set_floatstatus_divbyzero(); *out = 0; } #if @neg@ else if (b == -1 && a < 0 && a == -a) { npy_set_floatstatus_overflow(); *out = a / b; } #endif else { #if @neg@ @type@ tmp; tmp = a / b; if (((a > 0) != (b > 0)) && (a % b != 0)) { tmp--; } *out = tmp; #else *out = a / b; #endif } } #define @name@_ctype_floor_divide @name@_ctype_divide static void @name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) { if (a == 0 || b == 0) { if (b == 0) npy_set_floatstatus_divbyzero(); *out = 0; return; } #if @neg@ else if ((a > 0) == (b > 0)) { *out = a % b; } else { /* handled like Python does */ *out = a % b; if (*out) *out += b; } #else *out = a % b; #endif } /**end repeat**/ /**begin repeat * * #name = byte, ubyte, short, ushort, int, uint, long, * ulong, longlong, ulonglong# * #otyp = npy_float*4, npy_double*6# */ #define @name@_ctype_true_divide(a, b, out) \ *(out) = ((@otyp@) (a)) / ((@otyp@) (b)); /**end repeat**/ /* b will always be positive in this call */ /**begin repeat * * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong# * #upc = BYTE, UBYTE, SHORT, USHORT, INT, UINT, * LONG, ULONG, LONGLONG, ULONGLONG# */ static void @name@_ctype_power(@type@ a, @type@ b, @type@ *out) { @type@ tmp; if (b == 0) { *out = 1; return; } if (a == 1) { *out = 1; return; } tmp = b & 1 ? a : 1; b >>= 1; while (b > 0) { a *= a; if (b & 1) { tmp *= a; } b >>= 1; } *out = tmp; } /**end repeat**/ /* QUESTION: Should we check for overflow / underflow in (l,r)shift? */ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong# */ /**begin repeat1 * #oper = and, xor, or, lshift, rshift# * #op = &, ^, |, <<, >># */ #define @name@_ctype_@oper@(arg1, arg2, out) *(out) = (arg1) @op@ (arg2) /**end repeat1**/ /**end repeat**/ /**begin repeat * #name = float, double, longdouble# * #type = npy_float, npy_double, npy_longdouble# * #c = f, , l# */ #define @name@_ctype_add(a, b, outp) *(outp) = (a) + (b) #define @name@_ctype_subtract(a, b, outp) *(outp) = (a) - (b) #define @name@_ctype_multiply(a, b, outp) *(outp) = (a) * (b) #define @name@_ctype_divide(a, b, outp) *(outp) = (a) / (b) #define @name@_ctype_true_divide @name@_ctype_divide static void @name@_ctype_floor_divide(@type@ a, @type@ b, @type@ *out) { @type@ mod; *out = npy_divmod@c@(a, b, &mod); } static void @name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) { npy_divmod@c@(a, b, out); } static void @name@_ctype_divmod(@type@ a, @type@ b, @type@ *out1, @type@ *out2) { *out1 = npy_divmod@c@(a, b, out2); } /**end repeat**/ #define half_ctype_add(a, b, outp) *(outp) = \ npy_float_to_half(npy_half_to_float(a) + npy_half_to_float(b)) #define half_ctype_subtract(a, b, outp) *(outp) = \ npy_float_to_half(npy_half_to_float(a) - npy_half_to_float(b)) #define half_ctype_multiply(a, b, outp) *(outp) = \ npy_float_to_half(npy_half_to_float(a) * npy_half_to_float(b)) #define half_ctype_divide(a, b, outp) *(outp) = \ npy_float_to_half(npy_half_to_float(a) / npy_half_to_float(b)) #define half_ctype_true_divide half_ctype_divide static void half_ctype_floor_divide(npy_half a, npy_half b, npy_half *out) { npy_half mod; *out = npy_half_divmod(a, b, &mod); } static void half_ctype_remainder(npy_half a, npy_half b, npy_half *out) { npy_half_divmod(a, b, out); } static void half_ctype_divmod(npy_half a, npy_half b, npy_half *out1, npy_half *out2) { *out1 = npy_half_divmod(a, b, out2); } /**begin repeat * #name = cfloat, cdouble, clongdouble# * #rname = float, double, longdouble# * #rtype = npy_float, npy_double, npy_longdouble# * #c = f,,l# */ #define @name@_ctype_add(a, b, outp) do{ \ (outp)->real = (a).real + (b).real; \ (outp)->imag = (a).imag + (b).imag; \ } while(0) #define @name@_ctype_subtract(a, b, outp) do{ \ (outp)->real = (a).real - (b).real; \ (outp)->imag = (a).imag - (b).imag; \ } while(0) #define @name@_ctype_multiply(a, b, outp) do{ \ (outp)->real = (a).real * (b).real - (a).imag * (b).imag; \ (outp)->imag = (a).real * (b).imag + (a).imag * (b).real; \ } while(0) /* Algorithm identical to that in loops.c.src, for consistency */ #define @name@_ctype_divide(a, b, outp) do{ \ @rtype@ in1r = (a).real; \ @rtype@ in1i = (a).imag; \ @rtype@ in2r = (b).real; \ @rtype@ in2i = (b).imag; \ @rtype@ in2r_abs = npy_fabs@c@(in2r); \ @rtype@ in2i_abs = npy_fabs@c@(in2i); \ if (in2r_abs >= in2i_abs) { \ if (in2r_abs == 0 && in2i_abs == 0) { \ /* divide by zero should yield a complex inf or nan */ \ (outp)->real = in1r/in2r_abs; \ (outp)->imag = in1i/in2i_abs; \ } \ else { \ @rtype@ rat = in2i/in2r; \ @rtype@ scl = 1.0@c@/(in2r + in2i*rat); \ (outp)->real = (in1r + in1i*rat)*scl; \ (outp)->imag = (in1i - in1r*rat)*scl; \ } \ } \ else { \ @rtype@ rat = in2r/in2i; \ @rtype@ scl = 1.0@c@/(in2i + in2r*rat); \ (outp)->real = (in1r*rat + in1i)*scl; \ (outp)->imag = (in1i*rat - in1r)*scl; \ } \ } while(0) #define @name@_ctype_true_divide @name@_ctype_divide #define @name@_ctype_floor_divide(a, b, outp) do { \ @rname@_ctype_floor_divide( \ ((a).real*(b).real + (a).imag*(b).imag), \ ((b).real*(b).real + (b).imag*(b).imag), \ &((outp)->real)); \ (outp)->imag = 0; \ } while(0) /**end repeat**/ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, long, ulong, * longlong, ulonglong, cfloat, cdouble, clongdouble# */ #define @name@_ctype_divmod(a, b, out, out2) { \ @name@_ctype_floor_divide(a, b, out); \ @name@_ctype_remainder(a, b, out2); \ } /**end repeat**/ /**begin repeat * #name = float, double, longdouble# * #type = npy_float, npy_double, npy_longdouble# */ static npy_@name@ (*_basic_@name@_pow)(@type@ a, @type@ b); static void @name@_ctype_power(@type@ a, @type@ b, @type@ *out) { *out = _basic_@name@_pow(a, b); } /**end repeat**/ static void half_ctype_power(npy_half a, npy_half b, npy_half *out) { const npy_float af = npy_half_to_float(a); const npy_float bf = npy_half_to_float(b); const npy_float outf = _basic_float_pow(af,bf); *out = npy_float_to_half(outf); } /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * float, double, longdouble# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_float, npy_double, npy_longdouble# * #uns = (0,1)*5,0*3# */ static void @name@_ctype_negative(@type@ a, @type@ *out) { #if @uns@ npy_set_floatstatus_overflow(); #endif *out = -a; } /**end repeat**/ static void half_ctype_negative(npy_half a, npy_half *out) { *out = a^0x8000u; } /**begin repeat * #name = cfloat, cdouble, clongdouble# * #type = npy_cfloat, npy_cdouble, npy_clongdouble# */ static void @name@_ctype_negative(@type@ a, @type@ *out) { out->real = -a.real; out->imag = -a.imag; } /**end repeat**/ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_double, npy_longdouble# */ static void @name@_ctype_positive(@type@ a, @type@ *out) { *out = a; } /**end repeat**/ /* * Get the nc_powf, nc_pow, and nc_powl functions from * the data area of the power ufunc in umathmodule. */ /**begin repeat * #name = cfloat, cdouble, clongdouble# * #type = npy_cfloat, npy_cdouble, npy_clongdouble# */ static void @name@_ctype_positive(@type@ a, @type@ *out) { out->real = a.real; out->imag = a.imag; } static void (*_basic_@name@_pow)(@type@ *, @type@ *, @type@ *); static void @name@_ctype_power(@type@ a, @type@ b, @type@ *out) { _basic_@name@_pow(&a, &b, out); } /**end repeat**/ /**begin repeat * #name = ubyte, ushort, uint, ulong, ulonglong# */ #define @name@_ctype_absolute @name@_ctype_positive /**end repeat**/ /**begin repeat * #name = byte, short, int, long, longlong# * #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong# */ static void @name@_ctype_absolute(@type@ a, @type@ *out) { *out = (a < 0 ? -a : a); } /**end repeat**/ /**begin repeat * #name = float, double, longdouble# * #type = npy_float, npy_double, npy_longdouble# * #c = f,,l# */ static void @name@_ctype_absolute(@type@ a, @type@ *out) { *out = npy_fabs@c@(a); } /**end repeat**/ static void half_ctype_absolute(npy_half a, npy_half *out) { *out = a&0x7fffu; } /**begin repeat * #name = cfloat, cdouble, clongdouble# * #type = npy_cfloat, npy_cdouble, npy_clongdouble# * #rtype = npy_float, npy_double, npy_longdouble# * #c = f,,l# */ static void @name@_ctype_absolute(@type@ a, @rtype@ *out) { *out = npy_cabs@c@(a); } /**end repeat**/ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, long, * ulong, longlong, ulonglong# */ #define @name@_ctype_invert(a, out) *(out) = ~a; /**end repeat**/ /*** END OF BASIC CODE **/ /* The general strategy for commutative binary operators is to * * 1) Convert the types to the common type if both are scalars (0 return) * 2) If both are not scalars use ufunc machinery (-2 return) * 3) If both are scalars but cannot be cast to the right type * return NotImplmented (-1 return) * * 4) Perform the function on the C-type. * 5) If an error condition occurred, check to see * what the current error-handling is and handle the error. * * 6) Construct and return the output scalar. */ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, longdouble, * cfloat, cdouble, clongdouble# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble# * #Name = Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong, * Half, Float, LongDouble, * CFloat, CDouble, CLongDouble# * #TYPE = NPY_BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UINT, * NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, * NPY_HALF, NPY_FLOAT, NPY_LONGDOUBLE, * NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE# */ static int _@name@_convert_to_ctype(PyObject *a, @type@ *arg1) { PyObject *temp; if (PyArray_IsScalar(a, @Name@)) { *arg1 = PyArrayScalar_VAL(a, @Name@); return 0; } else if (PyArray_IsScalar(a, Generic)) { PyArray_Descr *descr1; if (!PyArray_IsScalar(a, Number)) { return -1; } descr1 = PyArray_DescrFromTypeObject((PyObject *)Py_TYPE(a)); if (PyArray_CanCastSafely(descr1->type_num, @TYPE@)) { PyArray_CastScalarDirect(a, descr1, arg1, @TYPE@); Py_DECREF(descr1); return 0; } else { Py_DECREF(descr1); return -1; } } else if (PyArray_GetPriority(a, NPY_PRIORITY) > NPY_PRIORITY) { return -2; } else if ((temp = PyArray_ScalarFromObject(a)) != NULL) { int retval = _@name@_convert_to_ctype(temp, arg1); Py_DECREF(temp); return retval; } return -2; } /**end repeat**/ /* Same as above but added exact checks against known python types for speed */ /**begin repeat * #name = double# * #type = npy_double# * #Name = Double# * #TYPE = NPY_DOUBLE# * #PYCHECKEXACT = PyFloat_CheckExact# * #PYEXTRACTCTYPE = PyFloat_AS_DOUBLE# */ static int _@name@_convert_to_ctype(PyObject *a, @type@ *arg1) { PyObject *temp; if (@PYCHECKEXACT@(a)){ *arg1 = @PYEXTRACTCTYPE@(a); return 0; } if (PyArray_IsScalar(a, @Name@)) { *arg1 = PyArrayScalar_VAL(a, @Name@); return 0; } else if (PyArray_IsScalar(a, Generic)) { PyArray_Descr *descr1; if (!PyArray_IsScalar(a, Number)) { return -1; } descr1 = PyArray_DescrFromTypeObject((PyObject *)Py_TYPE(a)); if (PyArray_CanCastSafely(descr1->type_num, @TYPE@)) { PyArray_CastScalarDirect(a, descr1, arg1, @TYPE@); Py_DECREF(descr1); return 0; } else { Py_DECREF(descr1); return -1; } } else if (PyArray_GetPriority(a, NPY_PRIORITY) > NPY_PRIORITY) { return -2; } else if ((temp = PyArray_ScalarFromObject(a)) != NULL) { int retval = _@name@_convert_to_ctype(temp, arg1); Py_DECREF(temp); return retval; } return -2; } /**end repeat**/ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, cfloat, cdouble# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_double, npy_cfloat, npy_cdouble# */ static int _@name@_convert2_to_ctypes(PyObject *a, @type@ *arg1, PyObject *b, @type@ *arg2) { int ret; ret = _@name@_convert_to_ctype(a, arg1); if (ret < 0) { return ret; } ret = _@name@_convert_to_ctype(b, arg2); if (ret < 0) { return ret; } return 0; } /**end repeat**/ /**begin repeat * #name = longdouble, clongdouble# * #type = npy_longdouble, npy_clongdouble# */ static int _@name@_convert2_to_ctypes(PyObject *a, @type@ *arg1, PyObject *b, @type@ *arg2) { int ret; ret = _@name@_convert_to_ctype(a, arg1); if (ret < 0) { return ret; } ret = _@name@_convert_to_ctype(b, arg2); if (ret == -2) { ret = -3; } if (ret < 0) { return ret; } return 0; } /**end repeat**/ #if defined(NPY_PY3K) #define CODEGEN_SKIP_divide_FLAG #endif /**begin repeat * * #name = (byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong)*13, * (half, float, double, longdouble, * cfloat, cdouble, clongdouble)*6, * (half, float, double, longdouble)*2# * #Name = (Byte, UByte, Short, UShort, Int, UInt, * Long, ULong,LongLong,ULongLong)*13, * (Half, Float, Double, LongDouble, * CFloat, CDouble, CLongDouble)*6, * (Half, Float, Double, LongDouble)*2# * #type = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong)*13, * (npy_half, npy_float, npy_double, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble)*6, * (npy_half, npy_float, npy_double, npy_longdouble)*2# * * #oper = add*10, subtract*10, multiply*10, divide*10, remainder*10, * divmod*10, floor_divide*10, lshift*10, rshift*10, and*10, * or*10, xor*10, true_divide*10, * add*7, subtract*7, multiply*7, divide*7, floor_divide*7, true_divide*7, * divmod*4, remainder*4# * * #fperr = 1*70,0*50,1*10, * 1*42, * 1*8# * #twoout = 0*50,1*10,0*70, * 0*42, * 1*4,0*4# * #otype = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong)*12, * npy_float*4, npy_double*6, * (npy_half, npy_float, npy_double, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble)*6, * (npy_half, npy_float, npy_double, npy_longdouble)*2# * #OName = (Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong)*12, * Float*4, Double*6, * (Half, Float, Double, LongDouble, * CFloat, CDouble, CLongDouble)*6, * (Half, Float, Double, LongDouble)*2# */ #if !defined(CODEGEN_SKIP_@oper@_FLAG) static PyObject * @name@_@oper@(PyObject *a, PyObject *b) { PyObject *ret; @type@ arg1, arg2; /* * NOTE: In gcc >= 4.1, the compiler will reorder floating point * operations and floating point error state checks. In * particular, the arithmetic operations were being reordered * so that the errors weren't caught. Declaring this output * variable volatile was the minimal fix for the issue. * (Ticket #1671) */ volatile @otype@ out; #if @twoout@ @otype@ out2; PyObject *obj; #endif #if @fperr@ int retstatus; int first; #endif BINOP_GIVE_UP_IF_NEEDED(a, b, nb_@oper@, @name@_@oper@); switch(_@name@_convert2_to_ctypes(a, &arg1, b, &arg2)) { case 0: break; case -1: /* one of them can't be cast safely must be mixed-types*/ return PyArray_Type.tp_as_number->nb_@oper@(a,b); case -2: /* use default handling */ if (PyErr_Occurred()) { return NULL; } return PyGenericArrType_Type.tp_as_number->nb_@oper@(a,b); case -3: /* * special case for longdouble and clongdouble * because they have a recursive getitem in their dtype */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } #if @fperr@ PyUFunc_clearfperr(); #endif /* * here we do the actual calculation with arg1 and arg2 * as a function call. */ #if @twoout@ @name@_ctype_@oper@(arg1, arg2, (@otype@ *)&out, &out2); #else @name@_ctype_@oper@(arg1, arg2, (@otype@ *)&out); #endif #if @fperr@ /* Check status flag. If it is set, then look up what to do */ retstatus = PyUFunc_getfperr(); if (retstatus) { int bufsize, errmask; PyObject *errobj; if (PyUFunc_GetPyValues("@name@_scalars", &bufsize, &errmask, &errobj) < 0) { return NULL; } first = 1; if (PyUFunc_handlefperr(errmask, errobj, retstatus, &first)) { Py_XDECREF(errobj); return NULL; } Py_XDECREF(errobj); } #endif #if @twoout@ ret = PyTuple_New(2); if (ret == NULL) { return NULL; } obj = PyArrayScalar_New(@OName@); if (obj == NULL) { Py_DECREF(ret); return NULL; } PyArrayScalar_ASSIGN(obj, @OName@, out); PyTuple_SET_ITEM(ret, 0, obj); obj = PyArrayScalar_New(@OName@); if (obj == NULL) { Py_DECREF(ret); return NULL; } PyArrayScalar_ASSIGN(obj, @OName@, out2); PyTuple_SET_ITEM(ret, 1, obj); #else ret = PyArrayScalar_New(@OName@); if (ret == NULL) { return NULL; } PyArrayScalar_ASSIGN(ret, @OName@, out); #endif return ret; } #endif /**end repeat**/ #undef CODEGEN_SKIP_divide_FLAG #define _IS_ZERO(x) (x == 0) /**begin repeat * * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble# * * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_double, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble# * * #Name = Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong, * Half, Float, Double, LongDouble, * CFloat, CDouble, CLongDouble# * * #oname = float*4, double*6, half, float, double, longdouble, * cfloat, cdouble, clongdouble# * * #otype = npy_float*4, npy_double*6, npy_half, npy_float, * npy_double, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble# * * #OName = Float*4, Double*6, Half, Float, * Double, LongDouble, * CFloat, CDouble, CLongDouble# * * #isint = (1,0)*5,0*7# * #cmplx = 0*14,1*3# * #iszero = _IS_ZERO*10, npy_half_iszero, _IS_ZERO*6# * #zero = 0*10, NPY_HALF_ZERO, 0*6# * #one = 1*10, NPY_HALF_ONE, 1*6# */ #if @cmplx@ static PyObject * @name@_power(PyObject *a, PyObject *b, PyObject *modulo) { PyObject *ret; @type@ arg1, arg2; int retstatus; int first; @type@ out = {@zero@, @zero@}; BINOP_GIVE_UP_IF_NEEDED(a, b, nb_power, @name@_power); switch(_@name@_convert2_to_ctypes(a, &arg1, b, &arg2)) { case 0: break; case -1: /* can't cast both safely mixed-types? */ return PyArray_Type.tp_as_number->nb_power(a,b,modulo); case -2: /* use default handling */ if (PyErr_Occurred()) { return NULL; } return PyGenericArrType_Type.tp_as_number->nb_power(a,b,modulo); case -3: default: /* * special case for longdouble and clongdouble * because they have a recursive getitem in their dtype */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } if (modulo != Py_None) { /* modular exponentiation is not implemented (gh-8804) */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } PyUFunc_clearfperr(); /* * here we do the actual calculation with arg1 and arg2 * as a function call. */ if (@iszero@(arg2.real) && @iszero@(arg2.imag)) { out.real = @one@; out.imag = @zero@; } else { @name@_ctype_power(arg1, arg2, &out); } /* Check status flag. If it is set, then look up what to do */ retstatus = PyUFunc_getfperr(); if (retstatus) { int bufsize, errmask; PyObject *errobj; if (PyUFunc_GetPyValues("@name@_scalars", &bufsize, &errmask, &errobj) < 0) { return NULL; } first = 1; if (PyUFunc_handlefperr(errmask, errobj, retstatus, &first)) { Py_XDECREF(errobj); return NULL; } Py_XDECREF(errobj); } ret = PyArrayScalar_New(@Name@); if (ret == NULL) { return NULL; } PyArrayScalar_ASSIGN(ret, @Name@, out); return ret; } #elif @isint@ static PyObject * @name@_power(PyObject *a, PyObject *b, PyObject *modulo) { PyObject *ret; @type@ arg1, arg2, out; BINOP_GIVE_UP_IF_NEEDED(a, b, nb_power, @name@_power); switch(_@name@_convert2_to_ctypes(a, &arg1, b, &arg2)) { case 0: break; case -1: /* can't cast both safely mixed-types? */ return PyArray_Type.tp_as_number->nb_power(a,b,modulo); case -2: /* use default handling */ if (PyErr_Occurred()) { return NULL; } return PyGenericArrType_Type.tp_as_number->nb_power(a,b,modulo); case -3: default: /* * special case for longdouble and clongdouble * because they have a recursive getitem in their dtype */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } if (modulo != Py_None) { /* modular exponentiation is not implemented (gh-8804) */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } PyUFunc_clearfperr(); /* * here we do the actual calculation with arg1 and arg2 * as a function call. */ if (arg2 < 0) { PyErr_SetString(PyExc_ValueError, "Integers to negative integer powers are not allowed."); return NULL; } @name@_ctype_power(arg1, arg2, &out); ret = PyArrayScalar_New(@Name@); if (ret == NULL) { return NULL; } PyArrayScalar_ASSIGN(ret, @Name@, out); return ret; } #else static PyObject * @name@_power(PyObject *a, PyObject *b, PyObject *modulo) { PyObject *ret; @type@ arg1, arg2; int retstatus; int first; @type@ out = @zero@; BINOP_GIVE_UP_IF_NEEDED(a, b, nb_power, @name@_power); switch(_@name@_convert2_to_ctypes(a, &arg1, b, &arg2)) { case 0: break; case -1: /* can't cast both safely mixed-types? */ return PyArray_Type.tp_as_number->nb_power(a,b,modulo); case -2: /* use default handling */ if (PyErr_Occurred()) { return NULL; } return PyGenericArrType_Type.tp_as_number->nb_power(a,b,modulo); case -3: default: /* * special case for longdouble and clongdouble * because they have a recursive getitem in their dtype */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } if (modulo != Py_None) { /* modular exponentiation is not implemented (gh-8804) */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } PyUFunc_clearfperr(); /* * here we do the actual calculation with arg1 and arg2 * as a function call. */ if (@iszero@(arg2)) { out = @one@; } else { @name@_ctype_power(arg1, arg2, &out); } /* Check status flag. If it is set, then look up what to do */ retstatus = PyUFunc_getfperr(); if (retstatus) { int bufsize, errmask; PyObject *errobj; if (PyUFunc_GetPyValues("@name@_scalars", &bufsize, &errmask, &errobj) < 0) { return NULL; } first = 1; if (PyUFunc_handlefperr(errmask, errobj, retstatus, &first)) { Py_XDECREF(errobj); return NULL; } Py_XDECREF(errobj); } ret = PyArrayScalar_New(@Name@); if (ret == NULL) { return NULL; } PyArrayScalar_ASSIGN(ret, @Name@, out); return ret; } #endif /**end repeat**/ #undef _IS_ZERO /**begin repeat * * #name = cfloat, cdouble, clongdouble# * */ /**begin repeat1 * * #oper = divmod, remainder# * */ #define @name@_@oper@ NULL /**end repeat1**/ /**end repeat**/ /**begin repeat * * #name = half, float, double, longdouble, cfloat, cdouble, clongdouble# * */ /**begin repeat1 * * #oper = lshift, rshift, and, or, xor# * */ #define @name@_@oper@ NULL /**end repeat1**/ /**end repeat**/ /**begin repeat * #name = (byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble)*3, * * byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong# * * #type = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_double, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble)*3, * * npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong# * * #otype = (npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_double, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble)*2, * npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_double, npy_longdouble, * npy_float, npy_double, npy_longdouble, * * npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, npy_uint, * npy_long, npy_ulong, npy_longlong, npy_ulonglong# * * #OName = (Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong, * Half, Float, Double, LongDouble, * CFloat, CDouble, CLongDouble)*2, * Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong, * Half, Float, Double, LongDouble, * Float, Double, LongDouble, * * Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong# * * #oper = negative*17, positive*17, absolute*17, invert*10# */ static PyObject * @name@_@oper@(PyObject *a) { @type@ arg1; @otype@ out; PyObject *ret; switch(_@name@_convert_to_ctype(a, &arg1)) { case 0: break; case -1: /* can't cast both safely use different add function */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; case -2: /* use default handling */ if (PyErr_Occurred()) { return NULL; } return PyGenericArrType_Type.tp_as_number->nb_@oper@(a); } /* * here we do the actual calculation with arg1 and arg2 * make it a function call. */ @name@_ctype_@oper@(arg1, &out); ret = PyArrayScalar_New(@OName@); PyArrayScalar_ASSIGN(ret, @OName@, out); return ret; } /**end repeat**/ /**begin repeat * * #name = half, float, double, longdouble, cfloat, cdouble, clongdouble# */ #define @name@_invert NULL /**end repeat**/ #if defined(NPY_PY3K) #define NONZERO_NAME(prefix) prefix##bool #else #define NONZERO_NAME(prefix) prefix##nonzero #endif #define _IS_NONZERO(x) (x != 0) /**begin repeat * * #name = byte, ubyte, short, ushort, int, * uint, long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble# * #type = npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int, * npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong, * npy_half, npy_float, npy_double, npy_longdouble, * npy_cfloat, npy_cdouble, npy_clongdouble# * #simp = 1*14, 0*3# * #nonzero = _IS_NONZERO*10, !npy_half_iszero, _IS_NONZERO*6# */ static int NONZERO_NAME(@name@_)(PyObject *a) { int ret; @type@ arg1; if (_@name@_convert_to_ctype(a, &arg1) < 0) { if (PyErr_Occurred()) { return -1; } return PyGenericArrType_Type.tp_as_number->NONZERO_NAME(nb_)(a); } /* * here we do the actual calculation with arg1 and arg2 * make it a function call. */ #if @simp@ ret = @nonzero@(arg1); #else ret = (@nonzero@(arg1.real) || @nonzero@(arg1.imag)); #endif return ret; } /**end repeat**/ #undef _IS_NONZERO static int emit_complexwarning(void) { static PyObject *cls = NULL; if (cls == NULL) { PyObject *mod; mod = PyImport_ImportModule("numpy.core"); assert(mod != NULL); cls = PyObject_GetAttrString(mod, "ComplexWarning"); assert(cls != NULL); Py_DECREF(mod); } return PyErr_WarnEx(cls, "Casting complex values to real discards the imaginary part", 1); } /**begin repeat * * #name = byte, ubyte, short, ushort, int, * uint, long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble# * * #Name = Byte, UByte, Short, UShort, Int, * UInt, Long, ULong, LongLong, ULongLong, * Half, Float, Double, LongDouble, * CFloat, CDouble, CLongDouble# * * #cmplx = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1# * #sign = (signed, unsigned)*5, , , , , , , # * #ctype = long*8, PY_LONG_LONG*2, * double*3, npy_longdouble, double*2, npy_longdouble# * #to_ctype = , , , , , , , , , , npy_half_to_double, , , , , , # * #func = (PyLong_FromLong, PyLong_FromUnsignedLong)*4, * PyLong_FromLongLong, PyLong_FromUnsignedLongLong, * PyLong_FromDouble*3, npy_longdouble_to_PyLong, * PyLong_FromDouble*2, npy_longdouble_to_PyLong# */ static PyObject * @name@_int(PyObject *obj) { PyObject *long_result; #if @cmplx@ @sign@ @ctype@ x = @to_ctype@(PyArrayScalar_VAL(obj, @Name@).real); #else @sign@ @ctype@ x = @to_ctype@(PyArrayScalar_VAL(obj, @Name@)); #endif #if @cmplx@ if (emit_complexwarning() < 0) { return NULL; } #endif long_result = @func@(x); if (long_result == NULL){ return NULL; } #ifndef NPY_PY3K /* Invoke long.__int__ to try to downcast */ { PyObject *before_downcast = long_result; long_result = Py_TYPE(long_result)->tp_as_number->nb_int(long_result); Py_DECREF(before_downcast); } #endif return long_result; } /**end repeat**/ /**begin repeat * * #name = (byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble)*2# * #Name = (Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong, * Half, Float, Double, LongDouble, * CFloat, CDouble, CLongDouble)*2# * #cmplx = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1)*2# * #to_ctype = (, , , , , , , , , , npy_half_to_double, , , , , , )*2# * #which = long*17, float*17# * #func = (PyLong_FromLongLong, PyLong_FromUnsignedLongLong)*5, * PyLong_FromDouble*3, npy_longdouble_to_PyLong, * PyLong_FromDouble*2, npy_longdouble_to_PyLong, * PyFloat_FromDouble*17# */ static NPY_INLINE PyObject * @name@_@which@(PyObject *obj) { #if @cmplx@ if (emit_complexwarning() < 0) { return NULL; } return @func@(@to_ctype@(PyArrayScalar_VAL(obj, @Name@).real)); #else return @func@(@to_ctype@(PyArrayScalar_VAL(obj, @Name@))); #endif } /**end repeat**/ #if !defined(NPY_PY3K) /**begin repeat * * #name = (byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble)*2# * #oper = oct*17, hex*17# * #kind = (int*5, long*5, int*2, long*2, int, long*2)*2# * #cap = (Int*5, Long*5, Int*2, Long*2, Int, Long*2)*2# */ static PyObject * @name@_@oper@(PyObject *obj) { PyObject *pyint; pyint = @name@_@kind@(obj); if (pyint == NULL) { return NULL; } return Py@cap@_Type.tp_as_number->nb_@oper@(pyint); } /**end repeat**/ #endif /**begin repeat * #oper = le, ge, lt, gt, eq, ne# * #op = <=, >=, <, >, ==, !=# * #halfop = npy_half_le, npy_half_ge, npy_half_lt, * npy_half_gt, npy_half_eq, npy_half_ne# */ #define def_cmp_@oper@(arg1, arg2) (arg1 @op@ arg2) #define cmplx_cmp_@oper@(arg1, arg2) ((arg1.real == arg2.real) ? \ arg1.imag @op@ arg2.imag : \ arg1.real @op@ arg2.real) #define def_half_cmp_@oper@(arg1, arg2) @halfop@(arg1, arg2) /**end repeat**/ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble# * #simp = def*10, def_half, def*3, cmplx*3# */ static PyObject* @name@_richcompare(PyObject *self, PyObject *other, int cmp_op) { npy_@name@ arg1, arg2; int out=0; RICHCMP_GIVE_UP_IF_NEEDED(self, other); switch(_@name@_convert2_to_ctypes(self, &arg1, other, &arg2)) { case 0: break; case -1: /* can't cast both safely use different add function */ case -2: /* use ufunc */ if (PyErr_Occurred()) { return NULL; } return PyGenericArrType_Type.tp_richcompare(self, other, cmp_op); case -3: /* * special case for longdouble and clongdouble * because they have a recursive getitem in their dtype */ Py_INCREF(Py_NotImplemented); return Py_NotImplemented; } /* here we do the actual calculation with arg1 and arg2 */ switch (cmp_op) { case Py_EQ: out = @simp@_cmp_eq(arg1, arg2); break; case Py_NE: out = @simp@_cmp_ne(arg1, arg2); break; case Py_LE: out = @simp@_cmp_le(arg1, arg2); break; case Py_GE: out = @simp@_cmp_ge(arg1, arg2); break; case Py_LT: out = @simp@_cmp_lt(arg1, arg2); break; case Py_GT: out = @simp@_cmp_gt(arg1, arg2); break; } if (out) { PyArrayScalar_RETURN_TRUE; } else { PyArrayScalar_RETURN_FALSE; } } /**end repeat**/ /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble# **/ static PyNumberMethods @name@_as_number = { (binaryfunc)@name@_add, /*nb_add*/ (binaryfunc)@name@_subtract, /*nb_subtract*/ (binaryfunc)@name@_multiply, /*nb_multiply*/ #if defined(NPY_PY3K) #else (binaryfunc)@name@_divide, /*nb_divide*/ #endif (binaryfunc)@name@_remainder, /*nb_remainder*/ (binaryfunc)@name@_divmod, /*nb_divmod*/ (ternaryfunc)@name@_power, /*nb_power*/ (unaryfunc)@name@_negative, (unaryfunc)@name@_positive, /*nb_pos*/ (unaryfunc)@name@_absolute, /*nb_abs*/ #if defined(NPY_PY3K) (inquiry)@name@_bool, /*nb_bool*/ #else (inquiry)@name@_nonzero, /*nb_nonzero*/ #endif (unaryfunc)@name@_invert, /*nb_invert*/ (binaryfunc)@name@_lshift, /*nb_lshift*/ (binaryfunc)@name@_rshift, /*nb_rshift*/ (binaryfunc)@name@_and, /*nb_and*/ (binaryfunc)@name@_xor, /*nb_xor*/ (binaryfunc)@name@_or, /*nb_or*/ #if defined(NPY_PY3K) #else 0, /*nb_coerce*/ #endif (unaryfunc)@name@_int, /*nb_int*/ #if defined(NPY_PY3K) (unaryfunc)0, /*nb_reserved*/ #else (unaryfunc)@name@_long, /*nb_long*/ #endif (unaryfunc)@name@_float, /*nb_float*/ #if defined(NPY_PY3K) #else (unaryfunc)@name@_oct, /*nb_oct*/ (unaryfunc)@name@_hex, /*nb_hex*/ #endif 0, /*inplace_add*/ 0, /*inplace_subtract*/ 0, /*inplace_multiply*/ #if defined(NPY_PY3K) #else 0, /*inplace_divide*/ #endif 0, /*inplace_remainder*/ 0, /*inplace_power*/ 0, /*inplace_lshift*/ 0, /*inplace_rshift*/ 0, /*inplace_and*/ 0, /*inplace_xor*/ 0, /*inplace_or*/ (binaryfunc)@name@_floor_divide, /*nb_floor_divide*/ (binaryfunc)@name@_true_divide, /*nb_true_divide*/ 0, /*nb_inplace_floor_divide*/ 0, /*nb_inplace_true_divide*/ (unaryfunc)NULL, /*nb_index*/ }; /**end repeat**/ NPY_NO_EXPORT void add_scalarmath(void) { /**begin repeat * #name = byte, ubyte, short, ushort, int, uint, * long, ulong, longlong, ulonglong, * half, float, double, longdouble, * cfloat, cdouble, clongdouble# * #NAME = Byte, UByte, Short, UShort, Int, UInt, * Long, ULong, LongLong, ULongLong, * Half, Float, Double, LongDouble, * CFloat, CDouble, CLongDouble# **/ @name@_as_number.nb_index = Py@NAME@ArrType_Type.tp_as_number->nb_index; Py@NAME@ArrType_Type.tp_as_number = &(@name@_as_number); Py@NAME@ArrType_Type.tp_richcompare = @name@_richcompare; /**end repeat**/ } static int get_functions(PyObject * mm) { PyObject *obj; void **funcdata; char *signatures; int i, j; int ret = -1; /* Get the nc_pow functions */ /* Get the pow functions */ obj = PyObject_GetAttrString(mm, "power"); if (obj == NULL) { goto fail; } funcdata = ((PyUFuncObject *)obj)->data; signatures = ((PyUFuncObject *)obj)->types; i = 0; j = 0; while (signatures[i] != NPY_FLOAT) { i += 3; j++; } _basic_float_pow = funcdata[j]; _basic_double_pow = funcdata[j + 1]; _basic_longdouble_pow = funcdata[j + 2]; _basic_cfloat_pow = funcdata[j + 3]; _basic_cdouble_pow = funcdata[j + 4]; _basic_clongdouble_pow = funcdata[j + 5]; Py_DECREF(obj); return ret = 0; fail: Py_DECREF(mm); return ret; } NPY_NO_EXPORT int initscalarmath(PyObject * m) { if (get_functions(m) < 0) { return -1; } add_scalarmath(); return 0; }