diff options
Diffstat (limited to 'numpy/numarray/_capi.c')
-rw-r--r-- | numpy/numarray/_capi.c | 3337 |
1 files changed, 3337 insertions, 0 deletions
diff --git a/numpy/numarray/_capi.c b/numpy/numarray/_capi.c new file mode 100644 index 000000000..0ba1be221 --- /dev/null +++ b/numpy/numarray/_capi.c @@ -0,0 +1,3337 @@ +#include <Python.h> + +#define _libnumarray_MODULE +#include "numpy/libnumarray.h" +#include <float.h> + +#if defined(__GLIBC__) || defined(__APPLE__) || defined(__MINGW32__) +#include <fenv.h> +#elif defined(__CYGWIN__) +#include "numpy/fenv/fenv.h" +#include "numpy/fenv/fenv.c" +#endif + +static PyObject *pCfuncClass; +static PyTypeObject CfuncType; +static PyObject *pHandleErrorFunc; + +static int +deferred_libnumarray_init(void) +{ + static int initialized=0; + + if (initialized) return 0; + + pCfuncClass = (PyObject *) &CfuncType; + Py_INCREF(pCfuncClass); + + pHandleErrorFunc = + NA_initModuleGlobal("numpy.numarray.util", "handleError"); + + if (!pHandleErrorFunc) goto _fail; + + + /* _exit: */ + initialized = 1; + return 0; + + _fail: + initialized = 0; + return -1; +} + + + +/**********************************************************************/ +/* Buffer Utility Functions */ +/**********************************************************************/ + +static PyObject * +getBuffer( PyObject *obj) +{ + if (!obj) return PyErr_Format(PyExc_RuntimeError, + "NULL object passed to getBuffer()"); + if (obj->ob_type->tp_as_buffer == NULL) { + return PyObject_CallMethod(obj, "__buffer__", NULL); + } else { + Py_INCREF(obj); /* Since CallMethod returns a new object when it + succeeds, We'll need to DECREF later to free it. + INCREF ordinary buffers here so we don't have to + remember where the buffer came from at DECREF time. + */ + return obj; + } +} + +/* Either it defines the buffer API, or it is an instance which returns +a buffer when obj.__buffer__() is called */ +static int +isBuffer (PyObject *obj) +{ + PyObject *buf = getBuffer(obj); + int ans = 0; + if (buf) { + ans = buf->ob_type->tp_as_buffer != NULL; + Py_DECREF(buf); + } else { + PyErr_Clear(); + } + return ans; +} + +/**********************************************************************/ + +static int +getWriteBufferDataPtr(PyObject *buffobj, void **buff) +{ + int rval = -1; + PyObject *buff2; + if ((buff2 = getBuffer(buffobj))) + { + if (buff2->ob_type->tp_as_buffer->bf_getwritebuffer) + rval = buff2->ob_type->tp_as_buffer->bf_getwritebuffer(buff2, + 0, buff); + Py_DECREF(buff2); + } + return rval; +} + +/**********************************************************************/ + +static int +isBufferWriteable (PyObject *buffobj) +{ + void *ptr; + int rval = -1; + rval = getWriteBufferDataPtr(buffobj, &ptr); + if (rval == -1) + PyErr_Clear(); /* Since we're just "testing", it's not really an error */ + return rval != -1; +} + +/**********************************************************************/ + +static int +getReadBufferDataPtr(PyObject *buffobj, void **buff) +{ + int rval = -1; + PyObject *buff2; + if ((buff2 = getBuffer(buffobj))) { + if (buff2->ob_type->tp_as_buffer->bf_getreadbuffer) + rval = buff2->ob_type->tp_as_buffer->bf_getreadbuffer(buff2, + 0, buff); + Py_DECREF(buff2); + } + return rval; +} + +/**********************************************************************/ + +static int +getBufferSize(PyObject *buffobj) +{ + Py_ssize_t size=0; + PyObject *buff2; + if ((buff2 = getBuffer(buffobj))) + { + (void) buff2->ob_type->tp_as_buffer->bf_getsegcount(buff2, &size); + Py_DECREF(buff2); + } + else + size = -1; + return size; +} + + +static double numarray_zero = 0.0; + +static double raiseDivByZero(void) +{ + return 1.0/numarray_zero; +} + +static double raiseNegDivByZero(void) +{ + return -1.0/numarray_zero; +} + +static double num_log(double x) +{ + if (x == 0.0) + return raiseNegDivByZero(); + else + return log(x); +} + +static double num_log10(double x) +{ + if (x == 0.0) + return raiseNegDivByZero(); + else + return log10(x); +} + +static double num_pow(double x, double y) +{ + int z = (int) y; + if ((x < 0.0) && (y != z)) + return raiseDivByZero(); + else + return pow(x, y); +} + +/* Inverse hyperbolic trig functions from Numeric */ +static double num_acosh(double x) +{ + return log(x + sqrt((x-1.0)*(x+1.0))); +} + +static double num_asinh(double xx) +{ + double x; + int sign; + if (xx < 0.0) { + sign = -1; + x = -xx; + } + else { + sign = 1; + x = xx; + } + return sign*log(x + sqrt(x*x+1.0)); +} + +static double num_atanh(double x) +{ + return 0.5*log((1.0+x)/(1.0-x)); +} + +/* NUM_CROUND (in numcomplex.h) also calls num_round */ +static double num_round(double x) +{ + return (x >= 0) ? floor(x+0.5) : ceil(x-0.5); +} + + +/* The following routine is used in the event of a detected integer * +** divide by zero so that a floating divide by zero is generated. * +** This is done since numarray uses the floating point exception * +** sticky bits to detect errors. The last bit is an attempt to * +** prevent optimization of the divide by zero away, the input value * +** should always be 0 * +*/ + +static int int_dividebyzero_error(long value, long unused) { + double dummy; + dummy = 1./numarray_zero; + if (dummy) /* to prevent optimizer from eliminating expression */ + return 0; + else + return 1; +} + +/* Likewise for Integer overflows */ +#if defined(__GLIBC__) || defined(__APPLE__) || defined(__CYGWIN__) || defined(__MINGW32__) +static int int_overflow_error(Float64 value) { /* For x86_64 */ + feraiseexcept(FE_OVERFLOW); + return (int) value; +} +#else +static int int_overflow_error(Float64 value) { + double dummy; + dummy = pow(1.e10, fabs(value/2)); + if (dummy) /* to prevent optimizer from eliminating expression */ + return (int) value; + else + return 1; +} +#endif + +static int umult64_overflow(UInt64 a, UInt64 b) +{ + UInt64 ah, al, bh, bl, w, x, y, z; + + ah = (a >> 32); + al = (a & 0xFFFFFFFFL); + bh = (b >> 32); + bl = (b & 0xFFFFFFFFL); + + /* 128-bit product: z*2**64 + (x+y)*2**32 + w */ + w = al*bl; + x = bh*al; + y = ah*bl; + z = ah*bh; + + /* *c = ((x + y)<<32) + w; */ + return z || (x>>32) || (y>>32) || + (((x & 0xFFFFFFFFL) + (y & 0xFFFFFFFFL) + (w >> 32)) >> 32); +} + +static int smult64_overflow(Int64 a0, Int64 b0) +{ + UInt64 a, b; + UInt64 ah, al, bh, bl, w, x, y, z; + + /* Convert to non-negative quantities */ + if (a0 < 0) { a = -a0; } else { a = a0; } + if (b0 < 0) { b = -b0; } else { b = b0; } + + ah = (a >> 32); + al = (a & 0xFFFFFFFFL); + bh = (b >> 32); + bl = (b & 0xFFFFFFFFL); + + w = al*bl; + x = bh*al; + y = ah*bl; + z = ah*bh; + + /* + UInt64 c = ((x + y)<<32) + w; + if ((a0 < 0) ^ (b0 < 0)) + *c = -c; + else + *c = c + */ + + return z || (x>>31) || (y>>31) || + (((x & 0xFFFFFFFFL) + (y & 0xFFFFFFFFL) + (w >> 32)) >> 31); +} + + +static void +NA_Done(void) +{ + return; +} + +static PyArrayObject * +NA_NewAll(int ndim, maybelong *shape, NumarrayType type, + void *buffer, maybelong byteoffset, maybelong bytestride, + int byteorder, int aligned, int writeable) +{ + PyArrayObject *result = NA_NewAllFromBuffer( + ndim, shape, type, Py_None, byteoffset, bytestride, + byteorder, aligned, writeable); + + if (result) { + if (!NA_NumArrayCheck((PyObject *) result)) { + PyErr_Format( PyExc_TypeError, + "NA_NewAll: non-NumArray result"); + result = NULL; + } else { + if (buffer) { + memcpy(result->data, buffer, NA_NBYTES(result)); + } else { + memset(result->data, 0, NA_NBYTES(result)); + } + } + } + return result; +} + +static PyArrayObject * +NA_NewAllStrides(int ndim, maybelong *shape, maybelong *strides, + NumarrayType type, void *buffer, maybelong byteoffset, + int byteorder, int aligned, int writeable) +{ + int i; + PyArrayObject *result = NA_NewAll(ndim, shape, type, buffer, + byteoffset, 0, + byteorder, aligned, writeable); + for(i=0; i<ndim; i++) + result->strides[i] = strides[i]; + return result; +} + + +static PyArrayObject * +NA_New(void *buffer, NumarrayType type, int ndim, ...) +{ + int i; + maybelong shape[MAXDIM]; + va_list ap; + va_start(ap, ndim); + for(i=0; i<ndim; i++) + shape[i] = va_arg(ap, int); + va_end(ap); + return NA_NewAll(ndim, shape, type, buffer, 0, 0, + NA_ByteOrder(), 1, 1); +} + +static PyArrayObject * +NA_Empty(int ndim, maybelong *shape, NumarrayType type) +{ + return NA_NewAll(ndim, shape, type, NULL, 0, 0, + NA_ByteOrder(), 1, 1); +} + + +/* Create a new numarray which is initially a C_array, or which +references a C_array: aligned, !byteswapped, contiguous, ... +Call with buffer==NULL to allocate storage. +*/ +static PyArrayObject * +NA_vNewArray(void *buffer, NumarrayType type, int ndim, maybelong *shape) +{ + return (PyArrayObject *) NA_NewAll(ndim, shape, type, buffer, 0, 0, + NA_ByteOrder(), 1, 1); +} + +static PyArrayObject * +NA_NewArray(void *buffer, NumarrayType type, int ndim, ...) +{ + int i; + maybelong shape[MAXDIM]; + va_list ap; + va_start(ap, ndim); + for(i=0; i<ndim; i++) + shape[i] = va_arg(ap, int); /* literals will still be ints */ + va_end(ap); + return NA_vNewArray(buffer, type, ndim, shape); +} + + +static PyObject* +NA_ReturnOutput(PyObject *out, PyArrayObject *shadow) +{ + if ((out == Py_None) || (out == NULL)) { + /* default behavior: return shadow array as the result */ + return (PyObject *) shadow; + } else { + PyObject *rval; + /* specified output behavior: return None */ + /* del(shadow) --> out.copyFrom(shadow) */ + Py_DECREF(shadow); + Py_INCREF(Py_None); + rval = Py_None; + return rval; + } +} + +static long NA_getBufferPtrAndSize(PyObject *buffobj, int readonly, void **ptr) +{ + long rval; + if (readonly) + rval = getReadBufferDataPtr(buffobj, ptr); + else + rval = getWriteBufferDataPtr(buffobj, ptr); + return rval; +} + + +static int NA_checkIo(char *name, + int wantIn, int wantOut, int gotIn, int gotOut) +{ + if (wantIn != gotIn) { + PyErr_Format(_Error, + "%s: wrong # of input buffers. Expected %d. Got %d.", + name, wantIn, gotIn); + return -1; + } + if (wantOut != gotOut) { + PyErr_Format(_Error, + "%s: wrong # of output buffers. Expected %d. Got %d.", + name, wantOut, gotOut); + return -1; + } + return 0; +} + +static int NA_checkOneCBuffer(char *name, long niter, + void *buffer, long bsize, size_t typesize) +{ + Int64 lniter = niter, ltypesize = typesize; + + if (lniter*ltypesize > bsize) { + PyErr_Format(_Error, + "%s: access out of buffer. niter=%d typesize=%d bsize=%d", + name, (int) niter, (int) typesize, (int) bsize); + return -1; + } + if ((typesize <= sizeof(Float64)) && (((long) buffer) % typesize)) { + PyErr_Format(_Error, + "%s: buffer not aligned on %d byte boundary.", + name, (int) typesize); + return -1; + } + return 0; +} + + +static int NA_checkNCBuffers(char *name, int N, long niter, + void **buffers, long *bsizes, + Int8 *typesizes, Int8 *iters) +{ + int i; + for (i=0; i<N; i++) + if (NA_checkOneCBuffer(name, iters[i] ? iters[i] : niter, + buffers[i], bsizes[i], typesizes[i])) + return -1; + return 0; +} + +static int NA_checkOneStriding(char *name, long dim, maybelong *shape, + long offset, maybelong *stride, long buffersize, long itemsize, + int align) +{ + int i; + long omin=offset, omax=offset; + long alignsize = (itemsize <= sizeof(Float64) ? itemsize : sizeof(Float64)); + + if (align && (offset % alignsize)) { + PyErr_Format(_Error, + "%s: buffer not aligned on %d byte boundary.", + name, (int) alignsize); + return -1; + } + for(i=0; i<dim; i++) { + long strideN = stride[i] * (shape[i]-1); + long tmax = omax + strideN; + long tmin = omin + strideN; + if (shape[i]-1 >= 0) { /* Skip dimension == 0. */ + omax = MAX(omax, tmax); + omin = MIN(omin, tmin); + if (align && (ABS(stride[i]) % alignsize)) { + PyErr_Format(_Error, + "%s: stride %d not aligned on %d byte boundary.", + name, (int) stride[i], (int) alignsize); + return -1; + } + if (omax + itemsize > buffersize) { + PyErr_Format(_Error, + "%s: access beyond buffer. offset=%d buffersize=%d", + name, (int) (omax+itemsize-1), (int) buffersize); + return -1; + } + if (omin < 0) { + PyErr_Format(_Error, + "%s: access before buffer. offset=%d buffersize=%d", + name, (int) omin, (int) buffersize); + return -1; + } + } + } + return 0; +} + +/* Function to call standard C Ufuncs +** +** The C Ufuncs expect contiguous 1-d data numarray, input and output numarray +** iterate with standard increments of one data element over all numarray. +** (There are some exceptions like arrayrangexxx which use one or more of +** the data numarray as parameter or other sources of information and do not +** iterate over every buffer). +** +** Arguments: +** +** Number of iterations (simple integer value). +** Number of input numarray. +** Number of output numarray. +** Tuple of tuples, one tuple per input/output array. Each of these +** tuples consists of a buffer object and a byte offset to start. +** +** Returns None +*/ + + +static PyObject * +NA_callCUFuncCore(PyObject *self, + long niter, long ninargs, long noutargs, + PyObject **BufferObj, long *offset) +{ + CfuncObject *me = (CfuncObject *) self; + char *buffers[MAXARGS]; + long bsizes[MAXARGS]; + long i, pnargs = ninargs + noutargs; + UFUNC ufuncptr; + + if (pnargs > MAXARGS) + return PyErr_Format(PyExc_RuntimeError, "NA_callCUFuncCore: too many parameters"); + + if (!PyObject_IsInstance(self, (PyObject *) &CfuncType) + || me->descr.type != CFUNC_UFUNC) + return PyErr_Format(PyExc_TypeError, + "NA_callCUFuncCore: problem with cfunc."); + + for (i=0; i<pnargs; i++) { + int readonly = (i < ninargs); + if (offset[i] < 0) + return PyErr_Format(_Error, + "%s: invalid negative offset:%d for buffer[%d]", + me->descr.name, (int) offset[i], (int) i); + if ((bsizes[i] = NA_getBufferPtrAndSize(BufferObj[i], readonly, + (void *) &buffers[i])) < 0) + return PyErr_Format(_Error, + "%s: Problem with %s buffer[%d].", + me->descr.name, + readonly ? "read" : "write", (int) i); + buffers[i] += offset[i]; + bsizes[i] -= offset[i]; /* "shorten" buffer size by offset. */ + } + + ufuncptr = (UFUNC) me->descr.fptr; + + /* If it's not a self-checking ufunc, check arg count match, + buffer size, and alignment for all buffers */ + if (!me->descr.chkself && + (NA_checkIo(me->descr.name, + me->descr.wantIn, me->descr.wantOut, ninargs, noutargs) || + NA_checkNCBuffers(me->descr.name, pnargs, + niter, (void **) buffers, bsizes, + me->descr.sizes, me->descr.iters))) + return NULL; + + /* Since the parameters are valid, call the C Ufunc */ + if (!(*ufuncptr)(niter, ninargs, noutargs, (void **)buffers, bsizes)) { + Py_INCREF(Py_None); + return Py_None; + } else { + return NULL; + } +} + +static PyObject * +callCUFunc(PyObject *self, PyObject *args) { + PyObject *DataArgs, *ArgTuple; + long pnargs, ninargs, noutargs, niter, i; + CfuncObject *me = (CfuncObject *) self; + PyObject *BufferObj[MAXARGS]; + long offset[MAXARGS]; + + if (!PyArg_ParseTuple(args, "lllO", + &niter, &ninargs, &noutargs, &DataArgs)) + return PyErr_Format(_Error, + "%s: Problem with argument list", me->descr.name); + + /* check consistency of stated inputs/outputs and supplied buffers */ + pnargs = PyObject_Length(DataArgs); + if ((pnargs != (ninargs+noutargs)) || (pnargs > MAXARGS)) + return PyErr_Format(_Error, + "%s: wrong buffer count for function", me->descr.name); + + /* Unpack buffers and offsets, get data pointers */ + for (i=0; i<pnargs; i++) { + ArgTuple = PySequence_GetItem(DataArgs, i); + Py_DECREF(ArgTuple); + if (!PyArg_ParseTuple(ArgTuple, "Ol", &BufferObj[i], &offset[i])) + return PyErr_Format(_Error, + "%s: Problem with buffer/offset tuple", me->descr.name); + } + return NA_callCUFuncCore(self, niter, ninargs, noutargs, BufferObj, offset); +} + +static PyObject * +callStrideConvCFunc(PyObject *self, PyObject *args) { + PyObject *inbuffObj, *outbuffObj, *shapeObj; + PyObject *inbstridesObj, *outbstridesObj; + CfuncObject *me = (CfuncObject *) self; + int nshape, ninbstrides, noutbstrides; + maybelong shape[MAXDIM], inbstrides[MAXDIM], + outbstrides[MAXDIM], *outbstrides1 = outbstrides; + long inboffset, outboffset, nbytes=0; + + if (!PyArg_ParseTuple(args, "OOlOOlO|l", + &shapeObj, &inbuffObj, &inboffset, &inbstridesObj, + &outbuffObj, &outboffset, &outbstridesObj, + &nbytes)) { + return PyErr_Format(_Error, + "%s: Problem with argument list", + me->descr.name); + } + + nshape = NA_maybeLongsFromIntTuple(MAXDIM, shape, shapeObj); + if (nshape < 0) return NULL; + + ninbstrides = NA_maybeLongsFromIntTuple(MAXDIM, inbstrides, inbstridesObj); + if (ninbstrides < 0) return NULL; + + noutbstrides= NA_maybeLongsFromIntTuple(MAXDIM, outbstrides, outbstridesObj); + if (noutbstrides < 0) return NULL; + + if (nshape && (nshape != ninbstrides)) { + return PyErr_Format(_Error, + "%s: Missmatch between input iteration and strides tuples", + me->descr.name); + } + + if (nshape && (nshape != noutbstrides)) { + if (noutbstrides < 1 || + outbstrides[ noutbstrides - 1 ])/* allow 0 for reductions. */ + return PyErr_Format(_Error, + "%s: Missmatch between output " + "iteration and strides tuples", + me->descr.name); + } + + return NA_callStrideConvCFuncCore( + self, nshape, shape, + inbuffObj, inboffset, ninbstrides, inbstrides, + outbuffObj, outboffset, noutbstrides, outbstrides1, nbytes); +} + +static int +_NA_callStridingHelper(PyObject *aux, long dim, + long nnumarray, PyArrayObject *numarray[], char *data[], + CFUNC_STRIDED_FUNC f) +{ + int i, j, status=0; + dim -= 1; + for(i=0; i<numarray[0]->dimensions[dim]; i++) { + for (j=0; j<nnumarray; j++) + data[j] += numarray[j]->strides[dim]*i; + if (dim == 0) + status |= f(aux, nnumarray, numarray, data); + else + status |= _NA_callStridingHelper( + aux, dim, nnumarray, numarray, data, f); + for (j=0; j<nnumarray; j++) + data[j] -= numarray[j]->strides[dim]*i; + } + return status; +} + + +static PyObject * +callStridingCFunc(PyObject *self, PyObject *args) { + CfuncObject *me = (CfuncObject *) self; + PyObject *aux; + PyArrayObject *numarray[MAXARRAYS]; + char *data[MAXARRAYS]; + CFUNC_STRIDED_FUNC f; + int i; + + int nnumarray = PySequence_Length(args)-1; + if ((nnumarray < 1) || (nnumarray > MAXARRAYS)) + return PyErr_Format(_Error, "%s, too many or too few numarray.", + me->descr.name); + + aux = PySequence_GetItem(args, 0); + if (!aux) + return NULL; + + for(i=0; i<nnumarray; i++) { + PyObject *otemp = PySequence_GetItem(args, i+1); + if (!otemp) + return PyErr_Format(_Error, "%s couldn't get array[%d]", + me->descr.name, i); + if (!NA_NDArrayCheck(otemp)) + return PyErr_Format(PyExc_TypeError, + "%s arg[%d] is not an array.", + me->descr.name, i); + numarray[i] = (PyArrayObject *) otemp; + data[i] = numarray[i]->data; + Py_DECREF(otemp); + if (!NA_updateDataPtr(numarray[i])) + return NULL; + } + + /* Cast function pointer and perform stride operation */ + f = (CFUNC_STRIDED_FUNC) me->descr.fptr; + + if (_NA_callStridingHelper(aux, numarray[0]->nd, + nnumarray, numarray, data, f)) { + return NULL; + } else { + Py_INCREF(Py_None); + return Py_None; + } +} + +/* Convert a standard C numeric value to a Python numeric value. +** +** Handles both nonaligned and/or byteswapped C data. +** +** Input arguments are: +** +** Buffer object that contains the C numeric value. +** Offset (in bytes) into the buffer that the data is located at. +** The size of the C numeric data item in bytes. +** Flag indicating if the C data is byteswapped from the processor's +** natural representation. +** +** Returns a Python numeric value. +*/ + +static PyObject * +NumTypeAsPyValue(PyObject *self, PyObject *args) { + PyObject *bufferObj; + long offset, itemsize, byteswap, i, buffersize; + Py_complex temp; /* to hold copies of largest possible type */ + void *buffer; + char *tempptr; + CFUNCasPyValue funcptr; + CfuncObject *me = (CfuncObject *) self; + + if (!PyArg_ParseTuple(args, "Olll", + &bufferObj, &offset, &itemsize, &byteswap)) + return PyErr_Format(_Error, + "NumTypeAsPyValue: Problem with argument list"); + + if ((buffersize = NA_getBufferPtrAndSize(bufferObj, 1, &buffer)) < 0) + return PyErr_Format(_Error, + "NumTypeAsPyValue: Problem with array buffer"); + + if (offset < 0) + return PyErr_Format(_Error, + "NumTypeAsPyValue: invalid negative offset: %d", (int) offset); + + /* Guarantee valid buffer pointer */ + if (offset+itemsize > buffersize) + return PyErr_Format(_Error, + "NumTypeAsPyValue: buffer too small for offset and itemsize."); + + /* Do byteswapping. Guarantee double alignment by using temp. */ + tempptr = (char *) &temp; + if (!byteswap) { + for (i=0; i<itemsize; i++) + *(tempptr++) = *(((char *) buffer)+offset+i); + } else { + tempptr += itemsize-1; + for (i=0; i<itemsize; i++) + *(tempptr--) = *(((char *) buffer)+offset+i); + } + + funcptr = (CFUNCasPyValue) me->descr.fptr; + + /* Call function to build PyObject. Bad parameters to this function + may render call meaningless, but "temp" guarantees that its safe. */ + return (*funcptr)((void *)(&temp)); +} + +/* Convert a Python numeric value to a standard C numeric value. +** +** Handles both nonaligned and/or byteswapped C data. +** +** Input arguments are: +** +** The Python numeric value to be converted. +** Buffer object to contain the C numeric value. +** Offset (in bytes) into the buffer that the data is to be copied to. +** The size of the C numeric data item in bytes. +** Flag indicating if the C data is byteswapped from the processor's +** natural representation. +** +** Returns None +*/ + +static PyObject * +NumTypeFromPyValue(PyObject *self, PyObject *args) { + PyObject *bufferObj, *valueObj; + long offset, itemsize, byteswap, i, buffersize; + Py_complex temp; /* to hold copies of largest possible type */ + void *buffer; + char *tempptr; + CFUNCfromPyValue funcptr; + CfuncObject *me = (CfuncObject *) self; + + if (!PyArg_ParseTuple(args, "OOlll", + &valueObj, &bufferObj, &offset, &itemsize, &byteswap)) + return PyErr_Format(_Error, + "%s: Problem with argument list", me->descr.name); + + if ((buffersize = NA_getBufferPtrAndSize(bufferObj, 0, &buffer)) < 0) + return PyErr_Format(_Error, + "%s: Problem with array buffer (read only?)", me->descr.name); + + funcptr = (CFUNCfromPyValue) me->descr.fptr; + + /* Convert python object into "temp". Always safe. */ + if (!((*funcptr)(valueObj, (void *)( &temp)))) + return PyErr_Format(_Error, + "%s: Problem converting value", me->descr.name); + + /* Check buffer offset. */ + if (offset < 0) + return PyErr_Format(_Error, + "%s: invalid negative offset: %d", me->descr.name, (int) offset); + + if (offset+itemsize > buffersize) + return PyErr_Format(_Error, + "%s: buffer too small(%d) for offset(%d) and itemsize(%d)", + me->descr.name, (int) buffersize, (int) offset, (int) itemsize); + + /* Copy "temp" to array buffer. */ + tempptr = (char *) &temp; + if (!byteswap) { + for (i=0; i<itemsize; i++) + *(((char *) buffer)+offset+i) = *(tempptr++); + } else { + tempptr += itemsize-1; + for (i=0; i<itemsize; i++) + *(((char *) buffer)+offset+i) = *(tempptr--); + } + Py_INCREF(Py_None); + return Py_None; +} + +/* Handle "calling" the cfunc object at the python level. + Dispatch the call to the appropriate python-c wrapper based + on the cfunc type. Do this dispatch to avoid having to + check that python code didn't somehow create a mismatch between + cfunc and wrapper. +*/ +static PyObject * +cfunc_call(PyObject *self, PyObject *argsTuple, PyObject *argsDict) +{ + CfuncObject *me = (CfuncObject *) self; + switch(me->descr.type) { + case CFUNC_UFUNC: + return callCUFunc(self, argsTuple); + break; + case CFUNC_STRIDING: + return callStrideConvCFunc(self, argsTuple); + break; + case CFUNC_NSTRIDING: + return callStridingCFunc(self, argsTuple); + case CFUNC_FROM_PY_VALUE: + return NumTypeFromPyValue(self, argsTuple); + break; + case CFUNC_AS_PY_VALUE: + return NumTypeAsPyValue(self, argsTuple); + break; + default: + return PyErr_Format( _Error, + "cfunc_call: Can't dispatch cfunc '%s' with type: %d.", + me->descr.name, me->descr.type); + } +} + +staticforward PyTypeObject CfuncType; + +static void +cfunc_dealloc(PyObject* self) +{ + PyObject_Del(self); +} + +static PyObject * +cfunc_repr(PyObject *self) +{ + char buf[256]; + CfuncObject *me = (CfuncObject *) self; + sprintf(buf, "<cfunc '%s' at %08lx check-self:%d align:%d io:(%d, %d)>", + me->descr.name, (unsigned long ) me->descr.fptr, + me->descr.chkself, me->descr.align, + me->descr.wantIn, me->descr.wantOut); + return PyString_FromString(buf); +} + +static PyTypeObject CfuncType = { + PyObject_HEAD_INIT(NULL) + 0, + "Cfunc", + sizeof(CfuncObject), + 0, + cfunc_dealloc, /*tp_dealloc*/ + 0, /*tp_print*/ + 0, /*tp_getattr*/ + 0, /*tp_setattr*/ + 0, /*tp_compare*/ + cfunc_repr, /*tp_repr*/ + 0, /*tp_as_number*/ + 0, /*tp_as_sequence*/ + 0, /*tp_as_mapping*/ + 0, /*tp_hash */ + cfunc_call, /* tp_call */ +}; + +/* CfuncObjects are created at the c-level only. They ensure that each +cfunc is called via the correct python-c-wrapper as defined by its +CfuncDescriptor. The wrapper, in turn, does conversions and buffer size +and alignment checking. Allowing these to be created at the python level +would enable them to be created *wrong* at the python level, and thereby +enable python code to *crash* python. +*/ +static PyObject* +NA_new_cfunc(CfuncDescriptor *cfd) +{ + CfuncObject* cfunc; + + CfuncType.ob_type = &PyType_Type; /* Should be done once at init. + Do now since there is no init. */ + + cfunc = PyObject_New(CfuncObject, &CfuncType); + + if (!cfunc) { + return PyErr_Format(_Error, + "NA_new_cfunc: failed creating '%s'", + cfd->name); + } + + cfunc->descr = *cfd; + + return (PyObject*)cfunc; +} + +static int NA_add_cfunc(PyObject *dict, char *keystr, CfuncDescriptor *descr) +{ + PyObject *c = (PyObject *) NA_new_cfunc(descr); + if (!c) return -1; + return PyDict_SetItemString(dict, keystr, c); +} + +static PyArrayObject* +NA_InputArray(PyObject *a, NumarrayType t, int requires) +{ + PyArray_Descr *descr; + if (t == tAny) descr = NULL; + else descr = PyArray_DescrFromType(t); + return (PyArrayObject *) \ + PyArray_CheckFromAny(a, descr, 0, 0, requires, NULL); +} + +/* satisfies ensures that 'a' meets a set of requirements and matches +the specified type. +*/ +static int +satisfies(PyArrayObject *a, int requirements, NumarrayType t) +{ + int type_ok = (a->descr->type_num == t) || (t == tAny); + + if (PyArray_ISCARRAY(a)) + return type_ok; + if (PyArray_ISBYTESWAPPED(a) && (requirements & NUM_NOTSWAPPED)) + return 0; + if (!PyArray_ISALIGNED(a) && (requirements & NUM_ALIGNED)) + return 0; + if (!PyArray_ISCONTIGUOUS(a) && (requirements & NUM_CONTIGUOUS)) + return 0; + if (!PyArray_ISWRITABLE(a) && (requirements & NUM_WRITABLE)) + return 0; + if (requirements & NUM_COPY) + return 0; + return type_ok; +} + + +static PyArrayObject * +NA_OutputArray(PyObject *a, NumarrayType t, int requires) +{ + PyArray_Descr *dtype; + PyArrayObject *ret; + + if (!PyArray_Check(a) || !PyArray_ISWRITEABLE(a)) { + PyErr_Format(PyExc_TypeError, + "NA_OutputArray: only writeable arrays work for output."); + return NULL; + } + + if (satisfies((PyArrayObject *)a, requires, t)) { + Py_INCREF(a); + return (PyArrayObject *)a; + } + if (t == tAny) { + dtype = PyArray_DESCR(a); + Py_INCREF(dtype); + } + else { + dtype = PyArray_DescrFromType(t); + } + ret = (PyArrayObject *)PyArray_Empty(PyArray_NDIM(a), PyArray_DIMS(a), + dtype, 0); + ret->flags |= NPY_UPDATEIFCOPY; + ret->base = a; + PyArray_FLAGS(a) &= ~NPY_WRITEABLE; + Py_INCREF(a); + return ret; +} + + +/* NA_IoArray is a combination of NA_InputArray and NA_OutputArray. + +Unlike NA_OutputArray, if a temporary is required it is initialized to a copy +of the input array. + +Unlike NA_InputArray, deallocating any resulting temporary array results in a +copy from the temporary back to the original. +*/ +static PyArrayObject * +NA_IoArray(PyObject *a, NumarrayType t, int requires) +{ + PyArrayObject *shadow = NA_InputArray(a, t, requires | NPY_UPDATEIFCOPY ); + + if (!shadow) return NULL; + + /* Guard against non-writable, but otherwise satisfying requires. + In this case, shadow == a. + */ + if (!PyArray_ISWRITABLE(shadow)) { + PyErr_Format(PyExc_TypeError, + "NA_IoArray: I/O array must be writable array"); + PyArray_XDECREF_ERR(shadow); + return NULL; + } + + return shadow; +} + +/* NA_OptionalOutputArray works like NA_OutputArray, but handles the case +where the output array 'optional' is omitted entirely at the python level, +resulting in 'optional'==Py_None. When 'optional' is Py_None, the return +value is cloned (but with NumarrayType 't') from 'master', typically an input +array with the same shape as the output array. +*/ +static PyArrayObject * +NA_OptionalOutputArray(PyObject *optional, NumarrayType t, int requires, + PyArrayObject *master) +{ + if ((optional == Py_None) || (optional == NULL)) { + PyObject *rval; + PyArray_Descr *descr; + if (t == tAny) descr=NULL; + else descr = PyArray_DescrFromType(t); + rval = PyArray_FromArray( + master, descr, NUM_C_ARRAY | NUM_COPY | NUM_WRITABLE); + return (PyArrayObject *)rval; + } else { + return NA_OutputArray(optional, t, requires); + } +} + +Complex64 NA_get_Complex64(PyArrayObject *a, long offset) +{ + Complex32 v0; + Complex64 v; + + switch(a->descr->type_num) { + case tComplex32: + v0 = NA_GETP(a, Complex32, (NA_PTR(a)+offset)); + v.r = v0.r; + v.i = v0.i; + break; + case tComplex64: + v = NA_GETP(a, Complex64, (NA_PTR(a)+offset)); + break; + default: + v.r = NA_get_Float64(a, offset); + v.i = 0; + break; + } + return v; +} + +void NA_set_Complex64(PyArrayObject *a, long offset, Complex64 v) +{ + Complex32 v0; + + switch(a->descr->type_num) { + case tComplex32: + v0.r = v.r; + v0.i = v.i; + NA_SETP(a, Complex32, (NA_PTR(a)+offset), v0); + break; + case tComplex64: + NA_SETP(a, Complex64, (NA_PTR(a)+offset), v); + break; + default: + NA_set_Float64(a, offset, v.r); + break; + } +} + +Int64 NA_get_Int64(PyArrayObject *a, long offset) +{ + switch(a->descr->type_num) { + case tBool: + return NA_GETP(a, Bool, (NA_PTR(a)+offset)) != 0; + case tInt8: + return NA_GETP(a, Int8, (NA_PTR(a)+offset)); + case tUInt8: + return NA_GETP(a, UInt8, (NA_PTR(a)+offset)); + case tInt16: + return NA_GETP(a, Int16, (NA_PTR(a)+offset)); + case tUInt16: + return NA_GETP(a, UInt16, (NA_PTR(a)+offset)); + case tInt32: + return NA_GETP(a, Int32, (NA_PTR(a)+offset)); + case tUInt32: + return NA_GETP(a, UInt32, (NA_PTR(a)+offset)); + case tInt64: + return NA_GETP(a, Int64, (NA_PTR(a)+offset)); + case tUInt64: + return NA_GETP(a, UInt64, (NA_PTR(a)+offset)); + case tFloat32: + return NA_GETP(a, Float32, (NA_PTR(a)+offset)); + case tFloat64: + return NA_GETP(a, Float64, (NA_PTR(a)+offset)); + case tComplex32: + return NA_GETP(a, Float32, (NA_PTR(a)+offset)); + case tComplex64: + return NA_GETP(a, Float64, (NA_PTR(a)+offset)); + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_get_Int64", + a->descr->type_num); + PyErr_Print(); + } + return 0; /* suppress warning */ +} + +void NA_set_Int64(PyArrayObject *a, long offset, Int64 v) +{ + Bool b; + + switch(a->descr->type_num) { + case tBool: + b = (v != 0); + NA_SETP(a, Bool, (NA_PTR(a)+offset), b); + break; + case tInt8: NA_SETP(a, Int8, (NA_PTR(a)+offset), v); + break; + case tUInt8: NA_SETP(a, UInt8, (NA_PTR(a)+offset), v); + break; + case tInt16: NA_SETP(a, Int16, (NA_PTR(a)+offset), v); + break; + case tUInt16: NA_SETP(a, UInt16, (NA_PTR(a)+offset), v); + break; + case tInt32: NA_SETP(a, Int32, (NA_PTR(a)+offset), v); + break; + case tUInt32: NA_SETP(a, UInt32, (NA_PTR(a)+offset), v); + break; + case tInt64: NA_SETP(a, Int64, (NA_PTR(a)+offset), v); + break; + case tUInt64: NA_SETP(a, UInt64, (NA_PTR(a)+offset), v); + break; + case tFloat32: + NA_SETP(a, Float32, (NA_PTR(a)+offset), v); + break; + case tFloat64: + NA_SETP(a, Float64, (NA_PTR(a)+offset), v); + break; + case tComplex32: + NA_SETP(a, Float32, (NA_PTR(a)+offset), v); + NA_SETP(a, Float32, (NA_PTR(a)+offset+sizeof(Float32)), 0); + break; + case tComplex64: + NA_SETP(a, Float64, (NA_PTR(a)+offset), v); + NA_SETP(a, Float64, (NA_PTR(a)+offset+sizeof(Float64)), 0); + break; + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_set_Int64", + a->descr->type_num); + PyErr_Print(); + } +} + +/* NA_get_offset computes the offset specified by the set of indices. +If N > 0, the indices are taken from the outer dimensions of the array. +If N < 0, the indices are taken from the inner dimensions of the array. +If N == 0, the offset is 0. +*/ +long NA_get_offset(PyArrayObject *a, int N, ...) +{ + int i; + long offset = 0; + va_list ap; + va_start(ap, N); + if (N > 0) { /* compute offset of "outer" indices. */ + for(i=0; i<N; i++) + offset += va_arg(ap, long) * a->strides[i]; + } else { /* compute offset of "inner" indices. */ + N = -N; + for(i=0; i<N; i++) + offset += va_arg(ap, long) * a->strides[a->nd-N+i]; + } + va_end(ap); + return offset; +} + +Float64 NA_get_Float64(PyArrayObject *a, long offset) +{ + switch(a->descr->type_num) { + case tBool: + return NA_GETP(a, Bool, (NA_PTR(a)+offset)) != 0; + case tInt8: + return NA_GETP(a, Int8, (NA_PTR(a)+offset)); + case tUInt8: + return NA_GETP(a, UInt8, (NA_PTR(a)+offset)); + case tInt16: + return NA_GETP(a, Int16, (NA_PTR(a)+offset)); + case tUInt16: + return NA_GETP(a, UInt16, (NA_PTR(a)+offset)); + case tInt32: + return NA_GETP(a, Int32, (NA_PTR(a)+offset)); + case tUInt32: + return NA_GETP(a, UInt32, (NA_PTR(a)+offset)); + case tInt64: + return NA_GETP(a, Int64, (NA_PTR(a)+offset)); + #if HAS_UINT64 + case tUInt64: + return NA_GETP(a, UInt64, (NA_PTR(a)+offset)); + #endif + case tFloat32: + return NA_GETP(a, Float32, (NA_PTR(a)+offset)); + case tFloat64: + return NA_GETP(a, Float64, (NA_PTR(a)+offset)); + case tComplex32: /* Since real value is first */ + return NA_GETP(a, Float32, (NA_PTR(a)+offset)); + case tComplex64: /* Since real value is first */ + return NA_GETP(a, Float64, (NA_PTR(a)+offset)); + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_get_Float64", + a->descr->type_num); + } + return 0; /* suppress warning */ +} + +void NA_set_Float64(PyArrayObject *a, long offset, Float64 v) +{ + Bool b; + + switch(a->descr->type_num) { + case tBool: + b = (v != 0); + NA_SETP(a, Bool, (NA_PTR(a)+offset), b); + break; + case tInt8: NA_SETP(a, Int8, (NA_PTR(a)+offset), v); + break; + case tUInt8: NA_SETP(a, UInt8, (NA_PTR(a)+offset), v); + break; + case tInt16: NA_SETP(a, Int16, (NA_PTR(a)+offset), v); + break; + case tUInt16: NA_SETP(a, UInt16, (NA_PTR(a)+offset), v); + break; + case tInt32: NA_SETP(a, Int32, (NA_PTR(a)+offset), v); + break; + case tUInt32: NA_SETP(a, UInt32, (NA_PTR(a)+offset), v); + break; + case tInt64: NA_SETP(a, Int64, (NA_PTR(a)+offset), v); + break; + #if HAS_UINT64 + case tUInt64: NA_SETP(a, UInt64, (NA_PTR(a)+offset), v); + break; + #endif + case tFloat32: + NA_SETP(a, Float32, (NA_PTR(a)+offset), v); + break; + case tFloat64: + NA_SETP(a, Float64, (NA_PTR(a)+offset), v); + break; + case tComplex32: { + NA_SETP(a, Float32, (NA_PTR(a)+offset), v); + NA_SETP(a, Float32, (NA_PTR(a)+offset+sizeof(Float32)), 0); + break; + } + case tComplex64: { + NA_SETP(a, Float64, (NA_PTR(a)+offset), v); + NA_SETP(a, Float64, (NA_PTR(a)+offset+sizeof(Float64)), 0); + break; + } + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_set_Float64", + a->descr->type_num ); + PyErr_Print(); + } +} + + +Float64 NA_get1_Float64(PyArrayObject *a, long i) +{ + long offset = i * a->strides[0]; + return NA_get_Float64(a, offset); +} + +Float64 NA_get2_Float64(PyArrayObject *a, long i, long j) +{ + long offset = i * a->strides[0] + + j * a->strides[1]; + return NA_get_Float64(a, offset); +} + +Float64 NA_get3_Float64(PyArrayObject *a, long i, long j, long k) +{ + long offset = i * a->strides[0] + + j * a->strides[1] + + k * a->strides[2]; + return NA_get_Float64(a, offset); +} + +void NA_set1_Float64(PyArrayObject *a, long i, Float64 v) +{ + long offset = i * a->strides[0]; + NA_set_Float64(a, offset, v); +} + +void NA_set2_Float64(PyArrayObject *a, long i, long j, Float64 v) +{ + long offset = i * a->strides[0] + + j * a->strides[1]; + NA_set_Float64(a, offset, v); +} + +void NA_set3_Float64(PyArrayObject *a, long i, long j, long k, Float64 v) +{ + long offset = i * a->strides[0] + + j * a->strides[1] + + k * a->strides[2]; + NA_set_Float64(a, offset, v); +} + +Complex64 NA_get1_Complex64(PyArrayObject *a, long i) +{ + long offset = i * a->strides[0]; + return NA_get_Complex64(a, offset); +} + +Complex64 NA_get2_Complex64(PyArrayObject *a, long i, long j) +{ + long offset = i * a->strides[0] + + j * a->strides[1]; + return NA_get_Complex64(a, offset); +} + +Complex64 NA_get3_Complex64(PyArrayObject *a, long i, long j, long k) +{ + long offset = i * a->strides[0] + + j * a->strides[1] + + k * a->strides[2]; + return NA_get_Complex64(a, offset); +} + +void NA_set1_Complex64(PyArrayObject *a, long i, Complex64 v) +{ + long offset = i * a->strides[0]; + NA_set_Complex64(a, offset, v); +} + +void NA_set2_Complex64(PyArrayObject *a, long i, long j, Complex64 v) +{ + long offset = i * a->strides[0] + + j * a->strides[1]; + NA_set_Complex64(a, offset, v); +} + +void NA_set3_Complex64(PyArrayObject *a, long i, long j, long k, Complex64 v) +{ + long offset = i * a->strides[0] + + j * a->strides[1] + + k * a->strides[2]; + NA_set_Complex64(a, offset, v); +} + +Int64 NA_get1_Int64(PyArrayObject *a, long i) +{ + long offset = i * a->strides[0]; + return NA_get_Int64(a, offset); +} + +Int64 NA_get2_Int64(PyArrayObject *a, long i, long j) +{ + long offset = i * a->strides[0] + + j * a->strides[1]; + return NA_get_Int64(a, offset); +} + +Int64 NA_get3_Int64(PyArrayObject *a, long i, long j, long k) +{ + long offset = i * a->strides[0] + + j * a->strides[1] + + k * a->strides[2]; + return NA_get_Int64(a, offset); +} + +void NA_set1_Int64(PyArrayObject *a, long i, Int64 v) +{ + long offset = i * a->strides[0]; + NA_set_Int64(a, offset, v); +} + +void NA_set2_Int64(PyArrayObject *a, long i, long j, Int64 v) +{ + long offset = i * a->strides[0] + + j * a->strides[1]; + NA_set_Int64(a, offset, v); +} + +void NA_set3_Int64(PyArrayObject *a, long i, long j, long k, Int64 v) +{ + long offset = i * a->strides[0] + + j * a->strides[1] + + k * a->strides[2]; + NA_set_Int64(a, offset, v); +} + +/* SET_CMPLX could be made faster by factoring it into 3 seperate loops. +*/ +#define NA_SET_CMPLX(a, type, base, cnt, in) \ +{ \ + int i; \ + int stride = a->strides[ a->nd - 1]; \ + NA_SET1D(a, type, base, cnt, in); \ + base = NA_PTR(a) + offset + sizeof(type); \ + for(i=0; i<cnt; i++) { \ + NA_SETP(a, Float32, base, 0); \ + base += stride; \ + } \ +} + +static int +NA_get1D_Float64(PyArrayObject *a, long offset, int cnt, Float64*out) +{ + char *base = NA_PTR(a) + offset; + + switch(a->descr->type_num) { + case tBool: + NA_GET1D(a, Bool, base, cnt, out); + break; + case tInt8: + NA_GET1D(a, Int8, base, cnt, out); + break; + case tUInt8: + NA_GET1D(a, UInt8, base, cnt, out); + break; + case tInt16: + NA_GET1D(a, Int16, base, cnt, out); + break; + case tUInt16: + NA_GET1D(a, UInt16, base, cnt, out); + break; + case tInt32: + NA_GET1D(a, Int32, base, cnt, out); + break; + case tUInt32: + NA_GET1D(a, UInt32, base, cnt, out); + break; + case tInt64: + NA_GET1D(a, Int64, base, cnt, out); + break; + #if HAS_UINT64 + case tUInt64: + NA_GET1D(a, UInt64, base, cnt, out); + break; + #endif + case tFloat32: + NA_GET1D(a, Float32, base, cnt, out); + break; + case tFloat64: + NA_GET1D(a, Float64, base, cnt, out); + break; + case tComplex32: + NA_GET1D(a, Float32, base, cnt, out); + break; + case tComplex64: + NA_GET1D(a, Float64, base, cnt, out); + break; + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_get1D_Float64", + a->descr->type_num); + PyErr_Print(); + return -1; + } + return 0; +} + +static Float64 * +NA_alloc1D_Float64(PyArrayObject *a, long offset, int cnt) +{ + Float64 *result = PyMem_New(Float64, cnt); + if (!result) return NULL; + if (NA_get1D_Float64(a, offset, cnt, result) < 0) { + PyMem_Free(result); + return NULL; + } + return result; +} + +static int +NA_set1D_Float64(PyArrayObject *a, long offset, int cnt, Float64*in) +{ + char *base = NA_PTR(a) + offset; + + switch(a->descr->type_num) { + case tBool: + NA_SET1D(a, Bool, base, cnt, in); + break; + case tInt8: + NA_SET1D(a, Int8, base, cnt, in); + break; + case tUInt8: + NA_SET1D(a, UInt8, base, cnt, in); + break; + case tInt16: + NA_SET1D(a, Int16, base, cnt, in); + break; + case tUInt16: + NA_SET1D(a, UInt16, base, cnt, in); + break; + case tInt32: + NA_SET1D(a, Int32, base, cnt, in); + break; + case tUInt32: + NA_SET1D(a, UInt32, base, cnt, in); + break; + case tInt64: + NA_SET1D(a, Int64, base, cnt, in); + break; + #if HAS_UINT64 + case tUInt64: + NA_SET1D(a, UInt64, base, cnt, in); + break; + #endif + case tFloat32: + NA_SET1D(a, Float32, base, cnt, in); + break; + case tFloat64: + NA_SET1D(a, Float64, base, cnt, in); + break; + case tComplex32: + NA_SET_CMPLX(a, Float32, base, cnt, in); + break; + case tComplex64: + NA_SET_CMPLX(a, Float64, base, cnt, in); + break; + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_set1D_Float64", + a->descr->type_num); + PyErr_Print(); + return -1; + } + return 0; +} + +static int +NA_get1D_Int64(PyArrayObject *a, long offset, int cnt, Int64*out) +{ + char *base = NA_PTR(a) + offset; + + switch(a->descr->type_num) { + case tBool: + NA_GET1D(a, Bool, base, cnt, out); + break; + case tInt8: + NA_GET1D(a, Int8, base, cnt, out); + break; + case tUInt8: + NA_GET1D(a, UInt8, base, cnt, out); + break; + case tInt16: + NA_GET1D(a, Int16, base, cnt, out); + break; + case tUInt16: + NA_GET1D(a, UInt16, base, cnt, out); + break; + case tInt32: + NA_GET1D(a, Int32, base, cnt, out); + break; + case tUInt32: + NA_GET1D(a, UInt32, base, cnt, out); + break; + case tInt64: + NA_GET1D(a, Int64, base, cnt, out); + break; + case tUInt64: + NA_GET1D(a, UInt64, base, cnt, out); + break; + case tFloat32: + NA_GET1D(a, Float32, base, cnt, out); + break; + case tFloat64: + NA_GET1D(a, Float64, base, cnt, out); + break; + case tComplex32: + NA_GET1D(a, Float32, base, cnt, out); + break; + case tComplex64: + NA_GET1D(a, Float64, base, cnt, out); + break; + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_get1D_Int64", + a->descr->type_num); + PyErr_Print(); + return -1; + } + return 0; +} + +static Int64 * +NA_alloc1D_Int64(PyArrayObject *a, long offset, int cnt) +{ + Int64 *result = PyMem_New(Int64, cnt); + if (!result) return NULL; + if (NA_get1D_Int64(a, offset, cnt, result) < 0) { + PyMem_Free(result); + return NULL; + } + return result; +} + +static int +NA_set1D_Int64(PyArrayObject *a, long offset, int cnt, Int64*in) +{ + char *base = NA_PTR(a) + offset; + + switch(a->descr->type_num) { + case tBool: + NA_SET1D(a, Bool, base, cnt, in); + break; + case tInt8: + NA_SET1D(a, Int8, base, cnt, in); + break; + case tUInt8: + NA_SET1D(a, UInt8, base, cnt, in); + break; + case tInt16: + NA_SET1D(a, Int16, base, cnt, in); + break; + case tUInt16: + NA_SET1D(a, UInt16, base, cnt, in); + break; + case tInt32: + NA_SET1D(a, Int32, base, cnt, in); + break; + case tUInt32: + NA_SET1D(a, UInt32, base, cnt, in); + break; + case tInt64: + NA_SET1D(a, Int64, base, cnt, in); + break; + case tUInt64: + NA_SET1D(a, UInt64, base, cnt, in); + break; + case tFloat32: + NA_SET1D(a, Float32, base, cnt, in); + break; + case tFloat64: + NA_SET1D(a, Float64, base, cnt, in); + break; + case tComplex32: + NA_SET_CMPLX(a, Float32, base, cnt, in); + break; + case tComplex64: + NA_SET_CMPLX(a, Float64, base, cnt, in); + break; + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_set1D_Int64", + a->descr->type_num); + PyErr_Print(); + return -1; + } + return 0; +} + +static int +NA_get1D_Complex64(PyArrayObject *a, long offset, int cnt, Complex64*out) +{ + char *base = NA_PTR(a) + offset; + + switch(a->descr->type_num) { + case tComplex64: + NA_GET1D(a, Complex64, base, cnt, out); + break; + default: + PyErr_Format( PyExc_TypeError, + "Unsupported type %d in NA_get1D_Complex64", + a->descr->type_num); + PyErr_Print(); + return -1; + } + return 0; +} + +static int +NA_set1D_Complex64(PyArrayObject *a, long offset, int cnt, Complex64*in) +{ + char *base = NA_PTR(a) + offset; + + switch(a->descr->type_num) { + case tComplex64: + NA_SET1D(a, Complex64, base, cnt, in); + break; + default: + PyErr_Format( PyExc_TypeError, + "Unsupported type %d in NA_set1D_Complex64", + a->descr->type_num); + PyErr_Print(); + return -1; + } + return 0; +} + + +/* NA_ShapeEqual returns 1 if 'a' and 'b' have the same shape, 0 otherwise. +*/ +static int +NA_ShapeEqual(PyArrayObject *a, PyArrayObject *b) +{ + int i; + + if (!NA_NDArrayCheck((PyObject *) a) || + !NA_NDArrayCheck((PyObject*) b)) { + PyErr_Format( + PyExc_TypeError, + "NA_ShapeEqual: non-array as parameter."); + return -1; + } + if (a->nd != b->nd) + return 0; + for(i=0; i<a->nd; i++) + if (a->dimensions[i] != b->dimensions[i]) + return 0; + return 1; +} + +/* NA_ShapeLessThan returns 1 if a.shape[i] < b.shape[i] for all i, else 0. +If they have a different number of dimensions, it compares the innermost +overlapping dimensions of each. +*/ +static int +NA_ShapeLessThan(PyArrayObject *a, PyArrayObject *b) +{ + int i; + int mindim, aoff, boff; + if (!NA_NDArrayCheck((PyObject *) a) || + !NA_NDArrayCheck((PyObject *) b)) { + PyErr_Format(PyExc_TypeError, + "NA_ShapeLessThan: non-array as parameter."); + return -1; + } + mindim = MIN(a->nd, b->nd); + aoff = a->nd - mindim; + boff = b->nd - mindim; + for(i=0; i<mindim; i++) + if (a->dimensions[i+aoff] >= b->dimensions[i+boff]) + return 0; + return 1; +} + +static int +NA_ByteOrder(void) +{ + unsigned long byteorder_test; + byteorder_test = 1; + if (*((char *) &byteorder_test)) + return NUM_LITTLE_ENDIAN; + else + return NUM_BIG_ENDIAN; +} + +static Bool +NA_IeeeSpecial32( Float32 *f, Int32 *mask) +{ + return NA_IeeeMask32(*f, *mask); +} + +static Bool +NA_IeeeSpecial64( Float64 *f, Int32 *mask) +{ + return NA_IeeeMask64(*f, *mask); +} + +static PyArrayObject * +NA_updateDataPtr(PyArrayObject *me) +{ + return me; +} + + +#define ELEM(x) (sizeof(x)/sizeof(x[0])) + +typedef struct +{ + char *name; + int typeno; +} NumarrayTypeNameMapping; + +static NumarrayTypeNameMapping NumarrayTypeNameMap[] = { + {"Any", tAny}, + {"Bool", tBool}, + {"Int8", tInt8}, + {"UInt8", tUInt8}, + {"Int16", tInt16}, + {"UInt16", tUInt16}, + {"Int32", tInt32}, + {"UInt32", tUInt32}, + {"Int64", tInt64}, + {"UInt64", tUInt64}, + {"Float32", tFloat32}, + {"Float64", tFloat64}, + {"Complex32", tComplex32}, + {"Complex64", tComplex64}, + {"Object", tObject}, + {"Long", tLong}, +}; + + +/* Convert NumarrayType 'typeno' into the string of the type's name. */ +static char * +NA_typeNoToName(int typeno) +{ + int i; + PyObject *typeObj; + int typeno2; + + for(i=0; i<ELEM(NumarrayTypeNameMap); i++) + if (typeno == NumarrayTypeNameMap[i].typeno) + return NumarrayTypeNameMap[i].name; + + /* Handle Numeric typecodes */ + typeObj = NA_typeNoToTypeObject(typeno); + if (!typeObj) return 0; + typeno2 = NA_typeObjectToTypeNo(typeObj); + Py_DECREF(typeObj); + + return NA_typeNoToName(typeno2); +} + +/* Look up the NumarrayType which corresponds to typename */ + +static int +NA_nameToTypeNo(char *typename) +{ + int i; + for(i=0; i<ELEM(NumarrayTypeNameMap); i++) + if (!strcmp(typename, NumarrayTypeNameMap[i].name)) + return NumarrayTypeNameMap[i].typeno; + return -1; +} + +static PyObject * +getTypeObject(NumarrayType type) +{ + return (PyObject *)PyArray_DescrFromType(type); +} + + +static PyObject * +NA_typeNoToTypeObject(int typeno) +{ + PyObject *o; + o = getTypeObject(typeno); + if (o) Py_INCREF(o); + return o; +} + + +static PyObject * +NA_intTupleFromMaybeLongs(int len, maybelong *Longs) +{ + return PyArray_IntTupleFromIntp(len, Longs); +} + +static long +NA_maybeLongsFromIntTuple(int len, maybelong *arr, PyObject *sequence) +{ + return PyArray_IntpFromSequence(sequence, arr, len); +} + + +static int +NA_intTupleProduct(PyObject *shape, long *prod) +{ + int i, nshape, rval = -1; + + if (!PySequence_Check(shape)) { + PyErr_Format(PyExc_TypeError, + "NA_intSequenceProduct: object is not a sequence."); + goto _exit; + } + nshape = PySequence_Size(shape); + + for(i=0, *prod=1; i<nshape; i++) { + PyObject *obj = PySequence_GetItem(shape, i); + if (!obj || !(PyInt_Check(obj) || PyLong_Check(obj))) { + PyErr_Format(PyExc_TypeError, + "NA_intTupleProduct: non-integer in shape."); + Py_XDECREF(obj); + goto _exit; + } + *prod *= PyInt_AsLong(obj); + Py_DECREF(obj); + if (PyErr_Occurred()) + goto _exit; + } + rval = 0; + _exit: + return rval; +} + +static long +NA_isIntegerSequence(PyObject *sequence) +{ + PyObject *o; + long i, size, isInt = 1; + if (!sequence) { + isInt = -1; + goto _exit; + } + if (!PySequence_Check(sequence)) { + isInt = 0; + goto _exit; + } + if ((size = PySequence_Length(sequence)) < 0) { + isInt = -1; + goto _exit; + } + for(i=0; i<size; i++) { + o = PySequence_GetItem(sequence, i); + if (!PyInt_Check(o) && !PyLong_Check(o)) { + isInt = 0; + Py_XDECREF(o); + goto _exit; + } + Py_XDECREF(o); + } + _exit: + return isInt; +} + +static int +getShape(PyObject *a, maybelong *shape, int dims) +{ + long slen; + + if (PyString_Check(a)) { + PyErr_Format(PyExc_TypeError, + "getShape: numerical sequences can't contain strings."); + return -1; + } + + if (!PySequence_Check(a) || + (NA_NDArrayCheck(a) && (PyArray(a)->nd == 0))) + return dims; + slen = PySequence_Length(a); + if (slen < 0) { + PyErr_Format(_Error, + "getShape: couldn't get sequence length."); + return -1; + } + if (!slen) { + *shape = 0; + return dims+1; + } else if (dims < MAXDIM) { + PyObject *item0 = PySequence_GetItem(a, 0); + if (item0) { + *shape = PySequence_Length(a); + dims = getShape(item0, ++shape, dims+1); + Py_DECREF(item0); + } else { + PyErr_Format(_Error, + "getShape: couldn't get sequence item."); + return -1; + } + } else { + PyErr_Format(_Error, + "getShape: sequence object nested more than MAXDIM deep."); + return -1; + } + return dims; +} + + + +typedef enum { + NOTHING, + NUMBER, + SEQUENCE +} SequenceConstraint; + +static int +setArrayFromSequence(PyArrayObject *a, PyObject *s, int dim, long offset) +{ + SequenceConstraint mustbe = NOTHING; + int i, seqlen=-1, slen = PySequence_Length(s); + + if (dim > a->nd) { + PyErr_Format(PyExc_ValueError, + "setArrayFromSequence: sequence/array dimensions mismatch."); + return -1; + } + + if (slen != a->dimensions[dim]) { + PyErr_Format(PyExc_ValueError, + "setArrayFromSequence: sequence/array shape mismatch."); + return -1; + } + + for(i=0; i<slen; i++) { + PyObject *o = PySequence_GetItem(s, i); + if (!o) { + PyErr_SetString(_Error, + "setArrayFromSequence: Can't get a sequence item"); + return -1; + } else if ((NA_isPythonScalar(o) || + (NA_NumArrayCheck(o) && PyArray(o)->nd == 0)) && + ((mustbe == NOTHING) || (mustbe == NUMBER))) { + if (NA_setFromPythonScalar(a, offset, o) < 0) + return -2; + mustbe = NUMBER; + } else if (PyString_Check(o)) { + PyErr_SetString( PyExc_ValueError, + "setArrayFromSequence: strings can't define numeric numarray."); + return -3; + } else if (PySequence_Check(o)) { + if ((mustbe == NOTHING) || (mustbe == SEQUENCE)) { + if (mustbe == NOTHING) { + mustbe = SEQUENCE; + seqlen = PySequence_Length(o); + } else if (PySequence_Length(o) != seqlen) { + PyErr_SetString( + PyExc_ValueError, + "Nested sequences with different lengths."); + return -5; + } + setArrayFromSequence(a, o, dim+1, offset); + } else { + PyErr_SetString(PyExc_ValueError, + "Nested sequences with different lengths."); + return -4; + } + } else { + PyErr_SetString(PyExc_ValueError, "Invalid sequence."); + return -6; + } + Py_DECREF(o); + offset += a->strides[dim]; + } + return 0; +} + +static PyObject * +NA_setArrayFromSequence(PyArrayObject *a, PyObject *s) +{ + maybelong shape[MAXDIM]; + + if (!PySequence_Check(s)) + return PyErr_Format( PyExc_TypeError, + "NA_setArrayFromSequence: (array, seq) expected."); + + if (getShape(s, shape, 0) < 0) + return NULL; + + if (!NA_updateDataPtr(a)) + return NULL; + + if (setArrayFromSequence(a, s, 0, 0) < 0) + return NULL; + + Py_INCREF(Py_None); + return Py_None; +} + +enum { + BOOL_SCALAR, + INT_SCALAR, + LONG_SCALAR, + FLOAT_SCALAR, + COMPLEX_SCALAR +}; + + +static int +_NA_maxType(PyObject *seq, int limit) +{ + if (limit > MAXDIM) { + PyErr_Format( PyExc_ValueError, + "NA_maxType: sequence nested too deep." ); + return -1; + } + if (NA_NumArrayCheck(seq)) { + switch(PyArray(seq)->descr->type_num) { + case tBool: + return BOOL_SCALAR; + case tInt8: + case tUInt8: + case tInt16: + case tUInt16: + case tInt32: + case tUInt32: + return INT_SCALAR; + case tInt64: + case tUInt64: + return LONG_SCALAR; + case tFloat32: + case tFloat64: + return FLOAT_SCALAR; + case tComplex32: + case tComplex64: + return COMPLEX_SCALAR; + default: + PyErr_Format(PyExc_TypeError, + "Expecting a python numeric type, got something else."); + return -1; + } + } else if (PySequence_Check(seq) && !PyString_Check(seq)) { + long i, maxtype=BOOL_SCALAR, slen; + + slen = PySequence_Length(seq); + if (slen < 0) return -1; + + if (slen == 0) return INT_SCALAR; + + for(i=0; i<slen; i++) { + long newmax; + PyObject *o = PySequence_GetItem(seq, i); + if (!o) return -1; + newmax = _NA_maxType(o, limit+1); + if (newmax < 0) + return -1; + else if (newmax > maxtype) { + maxtype = newmax; + } + Py_DECREF(o); + } + return maxtype; + } else { +#if PY_VERSION_HEX >= 0x02030000 + if (PyBool_Check(seq)) + return BOOL_SCALAR; + else +#endif + if (PyInt_Check(seq)) + return INT_SCALAR; + else if (PyLong_Check(seq)) + return LONG_SCALAR; + else if (PyFloat_Check(seq)) + return FLOAT_SCALAR; + else if (PyComplex_Check(seq)) + return COMPLEX_SCALAR; + else { + PyErr_Format(PyExc_TypeError, + "Expecting a python numeric type, got something else."); + return -1; + } + } +} + +static int +NA_maxType(PyObject *seq) +{ + int rval; + rval = _NA_maxType(seq, 0); + return rval; +} + +static int +NA_isPythonScalar(PyObject *o) +{ + int rval; + rval = PyInt_Check(o) || + PyLong_Check(o) || + PyFloat_Check(o) || + PyComplex_Check(o) || + (PyString_Check(o) && (PyString_Size(o) == 1)); + return rval; +} + +#if (NPY_SIZEOF_INTP == 8) +#define PlatBigInt PyInt_FromLong +#define PlatBigUInt PyLong_FromUnsignedLong +#else +#define PlatBigInt PyLong_FromLongLong +#define PlatBigUInt PyLong_FromUnsignedLongLong +#endif + + +static PyObject * +NA_getPythonScalar(PyArrayObject *a, long offset) +{ + int type = a->descr->type_num; + PyObject *rval = NULL; + + switch(type) { + case tBool: + case tInt8: + case tUInt8: + case tInt16: + case tUInt16: + case tInt32: { + Int64 v = NA_get_Int64(a, offset); + rval = PyInt_FromLong(v); + break; + } + case tUInt32: { + Int64 v = NA_get_Int64(a, offset); + rval = PlatBigUInt(v); + break; + } + case tInt64: { + Int64 v = NA_get_Int64(a, offset); + rval = PlatBigInt( v); + break; + } + case tUInt64: { + Int64 v = NA_get_Int64(a, offset); + rval = PlatBigUInt( v); + break; + } + case tFloat32: + case tFloat64: { + Float64 v = NA_get_Float64(a, offset); + rval = PyFloat_FromDouble( v ); + break; + } + case tComplex32: + case tComplex64: + { + Complex64 v = NA_get_Complex64(a, offset); + rval = PyComplex_FromDoubles(v.r, v.i); + break; + } + default: + rval = PyErr_Format(PyExc_TypeError, + "NA_getPythonScalar: bad type %d\n", + type); + } + return rval; +} + +static int +NA_overflow(PyArrayObject *a, Float64 v) +{ + if ((a->flags & CHECKOVERFLOW) == 0) return 0; + + switch(a->descr->type_num) { + case tBool: + return 0; + case tInt8: + if ((v < -128) || (v > 127)) goto _fail; + return 0; + case tUInt8: + if ((v < 0) || (v > 255)) goto _fail; + return 0; + case tInt16: + if ((v < -32768) || (v > 32767)) goto _fail; + return 0; + case tUInt16: + if ((v < 0) || (v > 65535)) goto _fail; + return 0; + case tInt32: + if ((v < -2147483648.) || + (v > 2147483647.)) goto _fail; + return 0; + case tUInt32: + if ((v < 0) || (v > 4294967295.)) goto _fail; + return 0; + case tInt64: + if ((v < -9223372036854775808.) || + (v > 9223372036854775807.)) goto _fail; + return 0; + #if HAS_UINT64 + case tUInt64: + if ((v < 0) || + (v > 18446744073709551615.)) goto _fail; + return 0; + #endif + case tFloat32: + if ((v < -FLT_MAX) || (v > FLT_MAX)) goto _fail; + return 0; + case tFloat64: + return 0; + case tComplex32: + if ((v < -FLT_MAX) || (v > FLT_MAX)) goto _fail; + return 0; + case tComplex64: + return 0; + default: + PyErr_Format( PyExc_TypeError, + "Unknown type %d in NA_overflow", + a->descr->type_num ); + PyErr_Print(); + return -1; + } + _fail: + PyErr_Format(PyExc_OverflowError, "value out of range for array"); + return -1; +} + +static int +_setFromPythonScalarCore(PyArrayObject *a, long offset, PyObject*value, int entries) +{ + Int64 v; + if (entries >= 100) { + PyErr_Format(PyExc_RuntimeError, + "NA_setFromPythonScalar: __tonumtype__ conversion chain too long"); + return -1; + } else if (PyInt_Check(value)) { + v = PyInt_AsLong(value); + if (NA_overflow(a, v) < 0) + return -1; + NA_set_Int64(a, offset, v); + } else if (PyLong_Check(value)) { + if (a->descr->type_num == tInt64) { + v = (Int64) PyLong_AsLongLong( value ); + } else if (a->descr->type_num == tUInt64) { + v = (UInt64) PyLong_AsUnsignedLongLong( value ); + } else if (a->descr->type_num == tUInt32) { + v = PyLong_AsUnsignedLong(value); + } else { + v = PyLong_AsLongLong(value); + } + if (PyErr_Occurred()) + return -1; + if (NA_overflow(a, v) < 0) + return -1; + NA_set_Int64(a, offset, v); + } else if (PyFloat_Check(value)) { + Float64 v = PyFloat_AsDouble(value); + if (NA_overflow(a, v) < 0) + return -1; + NA_set_Float64(a, offset, v); + } else if (PyComplex_Check(value)) { + Complex64 vc; + vc.r = PyComplex_RealAsDouble(value); + vc.i = PyComplex_ImagAsDouble(value); + if (NA_overflow(a, vc.r) < 0) + return -1; + if (NA_overflow(a, vc.i) < 0) + return -1; + NA_set_Complex64(a, offset, vc); + } else if (PyObject_HasAttrString(value, "__tonumtype__")) { + int rval; + PyObject *type = NA_typeNoToTypeObject(a->descr->type_num); + if (!type) return -1; + value = PyObject_CallMethod( + value, "__tonumtype__", "(N)", type); + if (!value) return -1; + rval = _setFromPythonScalarCore(a, offset, value, entries+1); + Py_DECREF(value); + return rval; + } else if (PyString_Check(value)) { + long size = PyString_Size(value); + if ((size <= 0) || (size > 1)) { + PyErr_Format( PyExc_ValueError, + "NA_setFromPythonScalar: len(string) must be 1."); + return -1; + } + NA_set_Int64(a, offset, *PyString_AsString(value)); + } else { + PyErr_Format(PyExc_TypeError, + "NA_setFromPythonScalar: bad value type."); + return -1; + } + return 0; +} + +static int +NA_setFromPythonScalar(PyArrayObject *a, long offset, PyObject *value) +{ + if (a->flags & WRITABLE) + return _setFromPythonScalarCore(a, offset, value, 0); + else { + PyErr_Format( + PyExc_ValueError, "NA_setFromPythonScalar: assigment to readonly array buffer"); + return -1; + } +} + + +static int +NA_NDArrayCheck(PyObject *obj) { + return PyArray_Check(obj); +} + +static int +NA_NumArrayCheck(PyObject *obj) { + return PyArray_Check(obj); +} + +static int +NA_ComplexArrayCheck(PyObject *a) +{ + int rval = NA_NumArrayCheck(a); + if (rval > 0) { + PyArrayObject *arr = (PyArrayObject *) a; + switch(arr->descr->type_num) { + case tComplex64: case tComplex32: + return 1; + default: + return 0; + } + } + return rval; +} + +static unsigned long +NA_elements(PyArrayObject *a) +{ + int i; + unsigned long n = 1; + for(i = 0; i<a->nd; i++) + n *= a->dimensions[i]; + return n; +} + +static int +NA_typeObjectToTypeNo(PyObject *typeObj) +{ + PyArray_Descr *dtype; + int i; + if (PyArray_DescrConverter(typeObj, &dtype) == NPY_FAIL) i=-1; + else i=dtype->type_num; + return i; +} + +static int +NA_copyArray(PyArrayObject *to, const PyArrayObject *from) +{ + return PyArray_CopyInto(to, (PyArrayObject *)from); +} + +static PyArrayObject * +NA_copy(PyArrayObject *from) +{ + return (PyArrayObject *)PyArray_NewCopy(from, 0); +} + + +static PyObject * +NA_getType( PyObject *type) +{ + PyArray_Descr *typeobj = NULL; + if (!type && PyArray_DescrConverter(type, &typeobj) == NPY_FAIL) { + PyErr_Format(PyExc_ValueError, "NA_getType: unknown type."); + typeobj = NULL; + } + return (PyObject *)typeobj; +} + + +/* Call a standard "stride" function +** +** Stride functions always take one input and one output array. +** They can handle n-dimensional data with arbitrary strides (of +** either sign) for both the input and output numarray. Typically +** these functions are used to copy data, byteswap, or align data. +** +** +** It expects the following arguments: +** +** Number of iterations for each dimension as a tuple +** Input Buffer Object +** Offset in bytes for input buffer +** Input strides (in bytes) for each dimension as a tuple +** Output Buffer Object +** Offset in bytes for output buffer +** Output strides (in bytes) for each dimension as a tuple +** An integer (Optional), typically the number of bytes to copy per +* element. +** +** Returns None +** +** The arguments expected by the standard stride functions that this +** function calls are: +** +** Number of dimensions to iterate over +** Long int value (from the optional last argument to +** callStrideConvCFunc) +** often unused by the C Function +** An array of long ints. Each is the number of iterations for each +** dimension. NOTE: the previous argument as well as the stride +** arguments are reversed in order with respect to how they are +** used in Python. Fastest changing dimension is the first element +** in the numarray! +** A void pointer to the input data buffer. +** The starting offset for the input data buffer in bytes (long int). +** An array of long int input strides (in bytes) [reversed as with +** the iteration array] +** A void pointer to the output data buffer. +** The starting offset for the output data buffer in bytes (long int). +** An array of long int output strides (in bytes) [also reversed] +*/ + + +static PyObject * +NA_callStrideConvCFuncCore( + PyObject *self, int nshape, maybelong *shape, + PyObject *inbuffObj, long inboffset, + int ninbstrides, maybelong *inbstrides, + PyObject *outbuffObj, long outboffset, + int noutbstrides, maybelong *outbstrides, + long nbytes) +{ + CfuncObject *me = (CfuncObject *) self; + CFUNC_STRIDE_CONV_FUNC funcptr; + void *inbuffer, *outbuffer; + long inbsize, outbsize; + maybelong i, lshape[MAXDIM], in_strides[MAXDIM], out_strides[MAXDIM]; + maybelong shape_0, inbstr_0, outbstr_0; + + if (nshape == 0) { /* handle rank-0 numarray. */ + nshape = 1; + shape = &shape_0; + inbstrides = &inbstr_0; + outbstrides = &outbstr_0; + shape[0] = 1; + inbstrides[0] = outbstrides[0] = 0; + } + + for(i=0; i<nshape; i++) + lshape[i] = shape[nshape-1-i]; + for(i=0; i<nshape; i++) + in_strides[i] = inbstrides[nshape-1-i]; + for(i=0; i<nshape; i++) + out_strides[i] = outbstrides[nshape-1-i]; + + if (!PyObject_IsInstance(self , (PyObject *) &CfuncType) + || me->descr.type != CFUNC_STRIDING) + return PyErr_Format(PyExc_TypeError, + "NA_callStrideConvCFuncCore: problem with cfunc"); + + if ((inbsize = NA_getBufferPtrAndSize(inbuffObj, 1, &inbuffer)) < 0) + return PyErr_Format(_Error, + "%s: Problem with input buffer", me->descr.name); + + if ((outbsize = NA_getBufferPtrAndSize(outbuffObj, 0, &outbuffer)) < 0) + return PyErr_Format(_Error, + "%s: Problem with output buffer (read only?)", + me->descr.name); + + /* Check buffer alignment and bounds */ + if (NA_checkOneStriding(me->descr.name, nshape, lshape, + inboffset, in_strides, inbsize, + (me->descr.sizes[0] == -1) ? + nbytes : me->descr.sizes[0], + me->descr.align) || + NA_checkOneStriding(me->descr.name, nshape, lshape, + outboffset, out_strides, outbsize, + (me->descr.sizes[1] == -1) ? + nbytes : me->descr.sizes[1], + me->descr.align)) + return NULL; + + /* Cast function pointer and perform stride operation */ + funcptr = (CFUNC_STRIDE_CONV_FUNC) me->descr.fptr; + if ((*funcptr)(nshape-1, nbytes, lshape, + inbuffer, inboffset, in_strides, + outbuffer, outboffset, out_strides) == 0) { + Py_INCREF(Py_None); + return Py_None; + } else { + return NULL; + } +} + +static void +NA_stridesFromShape(int nshape, maybelong *shape, maybelong bytestride, + maybelong *strides) +{ + int i; + if (nshape > 0) { + for(i=0; i<nshape; i++) + strides[i] = bytestride; + for(i=nshape-2; i>=0; i--) + strides[i] = strides[i+1]*shape[i+1]; + } +} + +static int +NA_OperatorCheck(PyObject *op) { + return 0; +} + +static int +NA_ConverterCheck(PyObject *op) { + return 0; +} + +static int +NA_UfuncCheck(PyObject *op) { + return 0; +} + +static int +NA_CfuncCheck(PyObject *op) { + return PyObject_TypeCheck(op, &CfuncType); +} + +static int +NA_getByteOffset(PyArrayObject *array, int nindices, maybelong *indices, + long *offset) +{ + return 0; +} + +static int +NA_swapAxes(PyArrayObject *array, int x, int y) +{ + long temp; + + if (((PyObject *) array) == Py_None) return 0; + + if (array->nd < 2) return 0; + + if (x < 0) x += array->nd; + if (y < 0) y += array->nd; + + if ((x < 0) || (x >= array->nd) || + (y < 0) || (y >= array->nd)) { + PyErr_Format(PyExc_ValueError, + "Specified dimension does not exist"); + return -1; + } + + temp = array->dimensions[x]; + array->dimensions[x] = array->dimensions[y]; + array->dimensions[y] = temp; + + temp = array->strides[x]; + array->strides[x] = array->strides[y]; + array->strides[y] = temp; + + PyArray_UpdateFlags(array, NPY_UPDATE_ALL); + + return 0; +} + +static PyObject * +NA_initModuleGlobal(char *modulename, char *globalname) +{ + PyObject *module, *dict, *global = NULL; + module = PyImport_ImportModule(modulename); + if (!module) { + PyErr_Format(PyExc_RuntimeError, + "Can't import '%s' module", + modulename); + goto _exit; + } + dict = PyModule_GetDict(module); + global = PyDict_GetItemString(dict, globalname); + if (!global) { + PyErr_Format(PyExc_RuntimeError, + "Can't find '%s' global in '%s' module.", + globalname, modulename); + goto _exit; + } + Py_DECREF(module); + Py_INCREF(global); + _exit: + return global; +} + +NumarrayType +NA_NumarrayType(PyObject *seq) +{ + int maxtype = NA_maxType(seq); + int rval; + switch(maxtype) { + case BOOL_SCALAR: + rval = tBool; + goto _exit; + case INT_SCALAR: + case LONG_SCALAR: + rval = tLong; /* tLong corresponds to C long int, + not Python long int */ + goto _exit; + case FLOAT_SCALAR: + rval = tFloat64; + goto _exit; + case COMPLEX_SCALAR: + rval = tComplex64; + goto _exit; + default: + PyErr_Format(PyExc_TypeError, + "expecting Python numeric scalar value; got something else."); + rval = -1; + } + _exit: + return rval; +} + +/* ignores bytestride */ +static PyArrayObject * +NA_NewAllFromBuffer(int ndim, maybelong *shape, NumarrayType type, + PyObject *bufferObject, maybelong byteoffset, maybelong bytestride, + int byteorder, int aligned, int writeable) +{ + PyArrayObject *self = NULL; + PyArray_Descr *dtype; + + if (type == tAny) + type = tDefault; + + dtype = PyArray_DescrFromType(type); + if (dtype == NULL) return NULL; + + if (byteorder != NA_ByteOrder()) { + PyArray_Descr *temp; + temp = PyArray_DescrNewByteorder(dtype, PyArray_SWAP); + Py_DECREF(dtype); + if (temp == NULL) return NULL; + dtype = temp; + } + + if (bufferObject == Py_None || bufferObject == NULL) { + self = (PyArrayObject *) \ + PyArray_NewFromDescr(&PyArray_Type, dtype, + ndim, shape, NULL, NULL, + 0, NULL); + } + else { + npy_intp size = 1; + int i; + PyArrayObject *newself; + PyArray_Dims newdims; + for(i=0; i<ndim; i++) { + size *= shape[i]; + } + self = (PyArrayObject *)\ + PyArray_FromBuffer(bufferObject, dtype, + size, byteoffset); + + if (self == NULL) return self; + newdims.len = ndim; + newdims.ptr = shape; + newself = (PyArrayObject *)\ + PyArray_Newshape(self, &newdims, PyArray_CORDER); + Py_DECREF(self); + self = newself; + } + + return self; +} + +static void +NA_updateAlignment(PyArrayObject *self) +{ + PyArray_UpdateFlags(self, NPY_ALIGNED); +} + +static void +NA_updateContiguous(PyArrayObject *self) +{ + PyArray_UpdateFlags(self, NPY_CONTIGUOUS | NPY_FORTRAN); +} + + +static void +NA_updateStatus(PyArrayObject *self) +{ + PyArray_UpdateFlags(self, NPY_UPDATE_ALL); +} + +static int +NA_NumArrayCheckExact(PyObject *op) { + return (op->ob_type == &PyArray_Type); +} + +static int +NA_NDArrayCheckExact(PyObject *op) { + return (op->ob_type == &PyArray_Type); +} + +static int +NA_OperatorCheckExact(PyObject *op) { + return 0; +} + +static int +NA_ConverterCheckExact(PyObject *op) { + return 0; +} + +static int +NA_UfuncCheckExact(PyObject *op) { + return 0; +} + + +static int +NA_CfuncCheckExact(PyObject *op) { + return op->ob_type == &CfuncType; +} + +static char * +NA_getArrayData(PyArrayObject *obj) +{ + if (!NA_NDArrayCheck((PyObject *) obj)) { + PyErr_Format(PyExc_TypeError, + "expected an NDArray"); + } + return obj->data; +} + +/* Byteswap is not a flag of the array --- it is implicit in the data-type */ +static void +NA_updateByteswap(PyArrayObject *self) +{ + return; +} + +static PyArray_Descr * +NA_DescrFromType(int type) +{ + if (type == tAny) + type = tDefault; + return PyArray_DescrFromType(type); +} + +static PyObject * +NA_Cast(PyArrayObject *a, int type) +{ + return PyArray_Cast(a, type); +} + + +/* The following function has much platform dependent code since +** there is no platform-independent way of checking Floating Point +** status bits +*/ + +/* OSF/Alpha (Tru64) ---------------------------------------------*/ +#if defined(__osf__) && defined(__alpha) + +static int +NA_checkFPErrors(void) +{ + unsigned long fpstatus; + int retstatus; + +#include <machine/fpu.h> /* Should migrate to global scope */ + + fpstatus = ieee_get_fp_control(); + /* clear status bits as well as disable exception mode if on */ + ieee_set_fp_control( 0 ); + retstatus = + pyFPE_DIVIDE_BY_ZERO* (int)((IEEE_STATUS_DZE & fpstatus) != 0) + + pyFPE_OVERFLOW * (int)((IEEE_STATUS_OVF & fpstatus) != 0) + + pyFPE_UNDERFLOW * (int)((IEEE_STATUS_UNF & fpstatus) != 0) + + pyFPE_INVALID * (int)((IEEE_STATUS_INV & fpstatus) != 0); + + return retstatus; +} + +/* MS Windows -----------------------------------------------------*/ +#elif defined(_MSC_VER) + +#include <float.h> + +static int +NA_checkFPErrors(void) +{ + int fpstatus = (int) _clear87(); + int retstatus = + pyFPE_DIVIDE_BY_ZERO * ((SW_ZERODIVIDE & fpstatus) != 0) + + pyFPE_OVERFLOW * ((SW_OVERFLOW & fpstatus) != 0) + + pyFPE_UNDERFLOW * ((SW_UNDERFLOW & fpstatus) != 0) + + pyFPE_INVALID * ((SW_INVALID & fpstatus) != 0); + + + return retstatus; +} + +/* Solaris --------------------------------------------------------*/ +/* --------ignoring SunOS ieee_flags approach, someone else can +** deal with that! */ +#elif defined(sun) +#include <ieeefp.h> + +static int +NA_checkFPErrors(void) +{ + int fpstatus; + int retstatus; + + fpstatus = (int) fpgetsticky(); + retstatus = pyFPE_DIVIDE_BY_ZERO * ((FP_X_DZ & fpstatus) != 0) + + pyFPE_OVERFLOW * ((FP_X_OFL & fpstatus) != 0) + + pyFPE_UNDERFLOW * ((FP_X_UFL & fpstatus) != 0) + + pyFPE_INVALID * ((FP_X_INV & fpstatus) != 0); + (void) fpsetsticky(0); + + return retstatus; +} + +#elif defined(__GLIBC__) || defined(__APPLE__) || defined(__CYGWIN__) || defined(__MINGW32__) + +static int +NA_checkFPErrors(void) +{ + int fpstatus = (int) fetestexcept( + FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW | FE_INVALID); + int retstatus = + pyFPE_DIVIDE_BY_ZERO * ((FE_DIVBYZERO & fpstatus) != 0) + + pyFPE_OVERFLOW * ((FE_OVERFLOW & fpstatus) != 0) + + pyFPE_UNDERFLOW * ((FE_UNDERFLOW & fpstatus) != 0) + + pyFPE_INVALID * ((FE_INVALID & fpstatus) != 0); + (void) feclearexcept(FE_DIVBYZERO | FE_OVERFLOW | + FE_UNDERFLOW | FE_INVALID); + return retstatus; +} + +#else + +static int +NA_checkFPErrors(void) +{ + return 0; +} + +#endif + +static void +NA_clearFPErrors() +{ + NA_checkFPErrors(); +} + +/* Not supported yet */ +static int +NA_checkAndReportFPErrors(char *name) +{ + int error = NA_checkFPErrors(); + if (error) { + PyObject *ans; + char msg[128]; + strcpy(msg, " in "); + strncat(msg, name, 100); + ans = PyObject_CallFunction(pHandleErrorFunc, "(is)", error, msg); + if (!ans) return -1; + Py_DECREF(ans); /* Py_None */ + } + return 0; + +} + + +#define WITHIN32(v, f) (((v) >= f##_MIN32) && ((v) <= f##_MAX32)) +#define WITHIN64(v, f) (((v) >= f##_MIN64) && ((v) <= f##_MAX64)) + +static Bool +NA_IeeeMask32( Float32 f, Int32 mask) +{ + Int32 category; + UInt32 v = *(UInt32 *) &f; + + if (v & BIT(31)) { + if (WITHIN32(v, NEG_NORMALIZED)) { + category = MSK_NEG_NOR; + } else if (WITHIN32(v, NEG_DENORMALIZED)) { + category = MSK_NEG_DEN; + } else if (WITHIN32(v, NEG_SIGNAL_NAN)) { + category = MSK_NEG_SNAN; + } else if (WITHIN32(v, NEG_QUIET_NAN)) { + category = MSK_NEG_QNAN; + } else if (v == NEG_INFINITY_MIN32) { + category = MSK_NEG_INF; + } else if (v == NEG_ZERO_MIN32) { + category = MSK_NEG_ZERO; + } else if (v == INDETERMINATE_MIN32) { + category = MSK_INDETERM; + } else { + category = MSK_BUG; + } + } else { + if (WITHIN32(v, POS_NORMALIZED)) { + category = MSK_POS_NOR; + } else if (WITHIN32(v, POS_DENORMALIZED)) { + category = MSK_POS_DEN; + } else if (WITHIN32(v, POS_SIGNAL_NAN)) { + category = MSK_POS_SNAN; + } else if (WITHIN32(v, POS_QUIET_NAN)) { + category = MSK_POS_QNAN; + } else if (v == POS_INFINITY_MIN32) { + category = MSK_POS_INF; + } else if (v == POS_ZERO_MIN32) { + category = MSK_POS_ZERO; + } else { + category = MSK_BUG; + } + } + return (category & mask) != 0; +} + +static Bool +NA_IeeeMask64( Float64 f, Int32 mask) +{ + Int32 category; + UInt64 v = *(UInt64 *) &f; + + if (v & BIT(63)) { + if (WITHIN64(v, NEG_NORMALIZED)) { + category = MSK_NEG_NOR; + } else if (WITHIN64(v, NEG_DENORMALIZED)) { + category = MSK_NEG_DEN; + } else if (WITHIN64(v, NEG_SIGNAL_NAN)) { + category = MSK_NEG_SNAN; + } else if (WITHIN64(v, NEG_QUIET_NAN)) { + category = MSK_NEG_QNAN; + } else if (v == NEG_INFINITY_MIN64) { + category = MSK_NEG_INF; + } else if (v == NEG_ZERO_MIN64) { + category = MSK_NEG_ZERO; + } else if (v == INDETERMINATE_MIN64) { + category = MSK_INDETERM; + } else { + category = MSK_BUG; + } + } else { + if (WITHIN64(v, POS_NORMALIZED)) { + category = MSK_POS_NOR; + } else if (WITHIN64(v, POS_DENORMALIZED)) { + category = MSK_POS_DEN; + } else if (WITHIN64(v, POS_SIGNAL_NAN)) { + category = MSK_POS_SNAN; + } else if (WITHIN64(v, POS_QUIET_NAN)) { + category = MSK_POS_QNAN; + } else if (v == POS_INFINITY_MIN64) { + category = MSK_POS_INF; + } else if (v == POS_ZERO_MIN64) { + category = MSK_POS_ZERO; + } else { + category = MSK_BUG; + } + } + return (category & mask) != 0; +} + +static PyArrayObject * +NA_FromDimsStridesDescrAndData(int nd, maybelong *d, maybelong *s, PyArray_Descr *descr, char *data) +{ + return (PyArrayObject *)\ + PyArray_NewFromDescr(&PyArray_Type, descr, nd, d, + s, data, 0, NULL); +} + +static PyArrayObject * +NA_FromDimsTypeAndData(int nd, maybelong *d, int type, char *data) +{ + PyArray_Descr *descr = NA_DescrFromType(type); + return NA_FromDimsStridesDescrAndData(nd, d, NULL, descr, data); +} + +static PyArrayObject * +NA_FromDimsStridesTypeAndData(int nd, maybelong *shape, maybelong *strides, + int type, char *data) +{ + PyArray_Descr *descr = NA_DescrFromType(type); + return NA_FromDimsStridesDescrAndData(nd, shape, strides, descr, data); +} + + +typedef struct +{ + NumarrayType type_num; + char suffix[5]; + int itemsize; +} scipy_typestr; + +static scipy_typestr scipy_descriptors[ ] = { + { tAny, ""}, + + { tBool, "b1", 1}, + + { tInt8, "i1", 1}, + { tUInt8, "u1", 1}, + + { tInt16, "i2", 2}, + { tUInt16, "u2", 2}, + + { tInt32, "i4", 4}, + { tUInt32, "u4", 4}, + + { tInt64, "i8", 8}, + { tUInt64, "u8", 8}, + + { tFloat32, "f4", 4}, + { tFloat64, "f8", 8}, + + { tComplex32, "c8", 8}, + { tComplex64, "c16", 16} +}; + + +static int +NA_scipy_typestr(NumarrayType t, int byteorder, char *typestr) +{ + int i; + if (byteorder) + strcpy(typestr, ">"); + else + strcpy(typestr, "<"); + for(i=0; i<ELEM(scipy_descriptors); i++) { + scipy_typestr *ts = &scipy_descriptors[i]; + if (ts->type_num == t) { + strncat(typestr, ts->suffix, 4); + return 0; + } + } + return -1; +} + +static PyArrayObject * +NA_FromArrayStruct(PyObject *obj) +{ + return (PyArrayObject *)PyArray_FromStructInterface(obj); +} + + +static PyObject *_Error; + +void *libnumarray_API[] = { + (void*) getBuffer, + (void*) isBuffer, + (void*) getWriteBufferDataPtr, + (void*) isBufferWriteable, + (void*) getReadBufferDataPtr, + (void*) getBufferSize, + (void*) num_log, + (void*) num_log10, + (void*) num_pow, + (void*) num_acosh, + (void*) num_asinh, + (void*) num_atanh, + (void*) num_round, + (void*) int_dividebyzero_error, + (void*) int_overflow_error, + (void*) umult64_overflow, + (void*) smult64_overflow, + (void*) NA_Done, + (void*) NA_NewAll, + (void*) NA_NewAllStrides, + (void*) NA_New, + (void*) NA_Empty, + (void*) NA_NewArray, + (void*) NA_vNewArray, + (void*) NA_ReturnOutput, + (void*) NA_getBufferPtrAndSize, + (void*) NA_checkIo, + (void*) NA_checkOneCBuffer, + (void*) NA_checkNCBuffers, + (void*) NA_checkOneStriding, + (void*) NA_new_cfunc, + (void*) NA_add_cfunc, + (void*) NA_InputArray, + (void*) NA_OutputArray, + (void*) NA_IoArray, + (void*) NA_OptionalOutputArray, + (void*) NA_get_offset, + (void*) NA_get_Float64, + (void*) NA_set_Float64, + (void*) NA_get_Complex64, + (void*) NA_set_Complex64, + (void*) NA_get_Int64, + (void*) NA_set_Int64, + (void*) NA_get1_Float64, + (void*) NA_get2_Float64, + (void*) NA_get3_Float64, + (void*) NA_set1_Float64, + (void*) NA_set2_Float64, + (void*) NA_set3_Float64, + (void*) NA_get1_Complex64, + (void*) NA_get2_Complex64, + (void*) NA_get3_Complex64, + (void*) NA_set1_Complex64, + (void*) NA_set2_Complex64, + (void*) NA_set3_Complex64, + (void*) NA_get1_Int64, + (void*) NA_get2_Int64, + (void*) NA_get3_Int64, + (void*) NA_set1_Int64, + (void*) NA_set2_Int64, + (void*) NA_set3_Int64, + (void*) NA_get1D_Float64, + (void*) NA_set1D_Float64, + (void*) NA_get1D_Int64, + (void*) NA_set1D_Int64, + (void*) NA_get1D_Complex64, + (void*) NA_set1D_Complex64, + (void*) NA_ShapeEqual, + (void*) NA_ShapeLessThan, + (void*) NA_ByteOrder, + (void*) NA_IeeeSpecial32, + (void*) NA_IeeeSpecial64, + (void*) NA_updateDataPtr, + (void*) NA_typeNoToName, + (void*) NA_nameToTypeNo, + (void*) NA_typeNoToTypeObject, + (void*) NA_intTupleFromMaybeLongs, + (void*) NA_maybeLongsFromIntTuple, + (void*) NA_intTupleProduct, + (void*) NA_isIntegerSequence, + (void*) NA_setArrayFromSequence, + (void*) NA_maxType, + (void*) NA_isPythonScalar, + (void*) NA_getPythonScalar, + (void*) NA_setFromPythonScalar, + (void*) NA_NDArrayCheck, + (void*) NA_NumArrayCheck, + (void*) NA_ComplexArrayCheck, + (void*) NA_elements, + (void*) NA_typeObjectToTypeNo, + (void*) NA_copyArray, + (void*) NA_copy, + (void*) NA_getType, + (void*) NA_callCUFuncCore, + (void*) NA_callStrideConvCFuncCore, + (void*) NA_stridesFromShape, + (void*) NA_OperatorCheck, + (void*) NA_ConverterCheck, + (void*) NA_UfuncCheck, + (void*) NA_CfuncCheck, + (void*) NA_getByteOffset, + (void*) NA_swapAxes, + (void*) NA_initModuleGlobal, + (void*) NA_NumarrayType, + (void*) NA_NewAllFromBuffer, + (void*) NA_alloc1D_Float64, + (void*) NA_alloc1D_Int64, + (void*) NA_updateAlignment, + (void*) NA_updateContiguous, + (void*) NA_updateStatus, + (void*) NA_NumArrayCheckExact, + (void*) NA_NDArrayCheckExact, + (void*) NA_OperatorCheckExact, + (void*) NA_ConverterCheckExact, + (void*) NA_UfuncCheckExact, + (void*) NA_CfuncCheckExact, + (void*) NA_getArrayData, + (void*) NA_updateByteswap, + (void*) NA_DescrFromType, + (void*) NA_Cast, + (void*) NA_checkFPErrors, + (void*) NA_clearFPErrors, + (void*) NA_checkAndReportFPErrors, + (void*) NA_IeeeMask32, + (void*) NA_IeeeMask64, + (void*) _NA_callStridingHelper, + (void*) NA_FromDimsStridesDescrAndData, + (void*) NA_FromDimsTypeAndData, + (void*) NA_FromDimsStridesTypeAndData, + (void*) NA_scipy_typestr, + (void*) NA_FromArrayStruct +}; + +#if (!defined(METHOD_TABLE_EXISTS)) +static PyMethodDef _libnumarrayMethods[] = { + {NULL, NULL} /* Sentinel */ +}; +#endif + +/* boiler plate API init */ +PyMODINIT_FUNC init_capi(void) +{ + PyObject *m = Py_InitModule("_capi", _libnumarrayMethods); + PyObject *c_api_object; + + _Error = PyErr_NewException("numpy.numarray._capi.error", NULL, NULL); + + /* Create a CObject containing the API pointer array's address */ + c_api_object = PyCObject_FromVoidPtr((void *)libnumarray_API, NULL); + + if (c_api_object != NULL) { + /* Create a name for this object in the module's namespace */ + PyObject *d = PyModule_GetDict(m); + + PyDict_SetItemString(d, "_C_API", c_api_object); + PyDict_SetItemString(d, "error", _Error); + Py_DECREF(c_api_object); + } else { + return; + } + if (PyModule_AddObject(m, "__version__", + PyString_FromString("0.9")) < 0) return; + + if (_import_array() < 0) return; + deferred_libnumarray_init(); + return; +} + + |