Draft of libnumarray C-API.

1 files changed, 1554 insertions, 2564 deletions

diff --git a/numpy/numarray/_capi.c b/numpy/numarray/_capi.c
index 0e6b593c2..7c3ac1721 100644
--- a/numpy/numarray/_capi.c
+++ b/numpy/numarray/_capi.c
@@ -2,181 +2,24 @@
 #include <Python.h>
 
 #define _libnumarray_MODULE
-
-#include "Python.h"
-#include "pystate.h"
 #include "libnumarray.h"
-#include <stdio.h>
-#include <float.h>
-
-static PyObject *pNDArrayModule;
-static PyObject *pNDArrayMDict;
-static PyObject *pNDArrayClass;
-
-static PyObject *pNumArrayModule;
-static PyObject *pNumArrayMDict;
-static PyObject *pNumArrayClass;
-static PyObject *pNumArrayNewFunc;
-static PyObject *pNumArrayArrayFunc;
-
-static PyObject *pNumericTypesModule;
-static PyObject *pNumericTypesMDict;
-static PyObject *pNumericTypeClass;
-static PyObject *pNumericTypesTDict;
-
-static PyObject *pUfuncModule;
-static PyObject *pUfuncMDict;
-static PyObject *pUfuncClass;
-
-static PyObject *pCfuncClass;
-
-static PyObject *pConverterModule;
-static PyObject *pConverterMDict;
-static PyObject *pConverterClass;
-
-static PyObject *pOperatorModule;
-static PyObject *pOperatorMDict;
-static PyObject *pOperatorClass;
-
-static PyObject *pNewMemoryFunc;
-static PyObject *pHandleErrorFunc;
 
 static PyObject *pNumType[nNumarrayType];
-
+static PyObject *pCfuncClass;
 static PyTypeObject CfuncType;
+static PyObject *pHandleErrorFunc;
 
-static PyObject *pEmptyDict;
-static PyObject *pEmptyTuple;
-
-static PyObject *dealloc_list;  /* list of global references to DECREF at
-				 unload time, i.e.  when the module dict is
-				 destructed. */
-
-enum {
-	BOOL_SCALAR,
-	INT_SCALAR,
-	LONG_SCALAR,
-	FLOAT_SCALAR,
-	COMPLEX_SCALAR
-};
-
-static int initialized = 0;
-
-/* custom init function generally unuseable due to circular references */
-static int 
-libnumarray_init(void) 
-{
-	PyObject *m, *d;
-	initialized = 0;
-	if (!(dealloc_list = PyList_New(0)))
-		return -1;
-	if (!(m = PyImport_ImportModule("numarray.libnumarray")))
-		return -1;
-	d = PyModule_GetDict(m);
-	if (PyDict_SetItemString(d, "_dealloc_list", dealloc_list) < 0)
-		return -1;
-	Py_DECREF(dealloc_list);
-	Py_DECREF(m);		       				       
-	return 0; 
-}
-
-static PyObject *
-init_module(char *modulename, PyObject **pMDict)
-{
-	PyObject *pModule = PyImport_ImportModule(modulename);
-	if (!pModule) return NULL;
-	PyList_Append(dealloc_list, pModule);
-	Py_DECREF(pModule);
-	*pMDict = PyModule_GetDict(pModule);
-	PyList_Append(dealloc_list, *pMDict);
-	return pModule;
-}
-
-static PyObject *
-init_object(char *objectname, PyObject *pMDict)
-{
-	PyObject *object = PyDict_GetItemString(pMDict, objectname);
-	if (!object) return NULL;
-	PyList_Append(dealloc_list, object);
-	return object;
-}
-
-static int
-init_module_class(char *modulename, PyObject **pModule, 
-		  PyObject **pMDict, 
-		  char *classname,  PyObject **pClass)
-{
-	if ((*pModule = init_module(modulename, pMDict)))
-		*pClass = init_object(classname, *pMDict);
-	else
-		return -1;
-	return 0;
-}
-
-
-extern void *libnumarray_API[];
-
- 
 static int 
 deferred_libnumarray_init(void)
 {
 	int i;
+        static initialized=0;
 
 	if (initialized) return 0;
 
-	import_libtc();
-
-	if (init_module_class("numarray.generic", &pNDArrayModule, 
-			      &pNDArrayMDict,
-			      "NDArray", &pNDArrayClass) < 0)
-		goto _fail;
-	
-	if (init_module_class("numarray", &pNumArrayModule, 
-			      &pNumArrayMDict,
-			      "NumArray", &pNumArrayClass) < 0)
-		goto _fail;
-	
-	if (init_module_class("numarray.numerictypes", &pNumericTypesModule, 
-			      &pNumericTypesMDict, 
-			      "NumericType", &pNumericTypeClass) < 0)
-		goto _fail;
-	
-	if (init_module_class("numarray._ufunc", &pUfuncModule, 
-			      &pUfuncMDict,
-			      "_ufunc", &pUfuncClass) < 0)
-		goto _fail;
-
 	pCfuncClass = (PyObject *) &CfuncType;
 	Py_INCREF(pCfuncClass);
-	
-	if (init_module_class("numarray._operator", &pOperatorModule, 
-			      &pOperatorMDict,
-			      "_operator", &pOperatorClass) < 0)
-		goto _fail;
-	
-	if (init_module_class("numarray._converter", &pConverterModule, 
-			      &pConverterMDict,
-			      "_converter", &pConverterClass) < 0)
-		goto _fail;
-
-	if (!(pNumArrayNewFunc = PyObject_GetAttrString(
-		      pNumArrayClass, "__new__")))
-		goto _fail;
-	
-	if (!(pNumArrayArrayFunc = init_object( "array", pNumArrayMDict)))
-		goto _fail;
-
-	if (!(pNumericTypesTDict = init_object( "typeDict", pNumericTypesMDict)))
-		goto _fail;
-
-	pNewMemoryFunc = NA_initModuleGlobal("numarray.memory","new_memory");
-	if (!pNewMemoryFunc) goto _fail;
-
-	pHandleErrorFunc = 
-		NA_initModuleGlobal("numarray.ufunc", "handleError");
-	if (!pHandleErrorFunc) goto _fail;
 		
-	
 	/* Set up table of type objects */
 	for(i=0; i<ELEM(pNumType); i++) {
 		PyObject *typeobj = init_object(NA_typeNoToName(i),
@@ -209,13 +52,11 @@ deferred_libnumarray_init(void)
 		}
 	}
 
-	libnumarray_API[ 0 ]  = (void *) pNumArrayClass;
+	pHandleErrorFunc = 
+		NA_initModuleGlobal("numpy.numarray.util", "handleError");
 
-	pEmptyDict = PyDict_New();
-	if (!pEmptyDict) goto _fail;
+	if (!pHandleErrorFunc) goto _fail;
 
-	pEmptyTuple = PyTuple_New(0);
-	if (!pEmptyTuple) goto _fail;
 
 	/* _exit: */
 	initialized = 1;
@@ -227,560 +68,271 @@ deferred_libnumarray_init(void)
 }
 
 
-/* Finalize this module. */
-void 
-fini_module_class(PyObject *module, PyObject *mdict, PyObject *class)
-{
-	Py_DECREF(module);
-	Py_DECREF(mdict);
-	Py_DECREF(class);
-}
-
-static void
-NA_Done(void)
-{
-	int i;
-
-	fini_module_class(pNDArrayModule, pNDArrayMDict, pNDArrayClass);
-
-	fini_module_class(pNumArrayModule, pNumArrayMDict, pNumArrayClass);
-	Py_DECREF(pNumArrayArrayFunc);
-
-	fini_module_class(pOperatorModule, pOperatorMDict, pOperatorClass);
-
-	fini_module_class(pConverterModule, pConverterMDict, pConverterClass);
 
-	fini_module_class(pUfuncModule, pUfuncMDict, pUfuncClass);
-	Py_DECREF(pCfuncClass);
-	
-	fini_module_class(pNumericTypesModule, pNumericTypesMDict, 
-			  pNumericTypeClass);
-	Py_DECREF(pNumericTypesTDict);
+/**********************************************************************/
+/*  Buffer Utility Functions                                          */
+/**********************************************************************/
 
-	for(i=0; i<ELEM(pNumType); i++) {
-		Py_DECREF(pNumType[i]);
+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;
 	}
 }
 
-#ifdef MS_WIN32
-#pragma warning(once : 4244)
-#endif
-
-#define ELEM(x) (sizeof(x)/sizeof(x[0]))
-
-typedef struct
-{
-	char *name;
-	int typeno;
-} NumarrayTypeNameMapping;
-
-static PyArray_Descr descriptors[ ] = {
-	{ tAny,       0,                 '*'}, 
-	
-	{ tBool,      sizeof(Bool),      '?'},
-
-	{ tInt8,      sizeof(Int8),      '1'},
-	{ tUInt8,     sizeof(UInt8),     'b'},
-
-	{ tInt16,     sizeof(Int16),     's'},
-	{ tUInt16,    sizeof(UInt16),    'w'},
-
-	{ tInt32,     sizeof(Int32),     'i'},
-	{ tUInt32,    sizeof(UInt32),    'u'},
-
-	{ tInt64,     sizeof(Int64),     'N'},
-	{ tUInt64,    sizeof(UInt64),    'U'},
-
-	{ tFloat32,   sizeof(Float32),   'f'},
-	{ tFloat64,   sizeof(Float64),   'd'},
-
-	{ tComplex32, sizeof(Complex32), 'F'},
-	{ tComplex64, sizeof(Complex64), 'D'}
-};
-
-static PyArray_Descr *
-NA_DescrFromType(int type)
+/* 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) 
 {
-	if ((type >= tAny) && (type <= tComplex64)) {
-		return &descriptors[ type ];
+	PyObject *buf = getBuffer(obj);
+	int ans = 0;
+	if (buf) {
+		ans = buf->ob_type->tp_as_buffer != NULL;
+		Py_DECREF(buf);
 	} else {
-		int i;
-		for(i=0; i<ELEM(descriptors); i++)
-			if (descriptors[i].type == type)
-				return &descriptors[i];
+		PyErr_Clear();
 	}
-	PyErr_Format(
-		PyExc_TypeError, 
-		"NA_DescrFromType: unknown type: %d", type);
-	return NULL;
+	return ans;
 }
 
-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},
-};
-
-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 PyObject *
-setTypeException(int type)
+static int 
+getWriteBufferDataPtr(PyObject *buffobj, void **buff) 
 {
-	/* Check if it is a printable character */
-	if ((type >= 32) && (type <= 126)) 
-		PyErr_Format(_Error, 
-			     "Type object lookup returned"
-			     " NULL for type \'%c\'", type);
-	else
-		PyErr_Format(_Error,
-			     "Type object lookup returned"
-			     " NULL for type %d", type);
-	return NULL;
+  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 PyObject *
-getTypeObject(NumarrayType type) 
-{
-        char strcharcode[2];
-        PyObject *typeobj;
-	
-	if (deferred_libnumarray_init() < 0) return NULL;
-
-        if ((type >= tAny) && (type <= tObject)) {
-		return pNumType[type];
-	} else  {
-	      /* Test if it is a Numeric charcode */
-		strcharcode[0] = type; strcharcode[1] = 0;
-		typeobj = PyDict_GetItemString(
-			pNumericTypesTDict, strcharcode);
-		return typeobj ? typeobj : setTypeException(type);
-	}
-}
+/**********************************************************************/
 
-static PyObject *
-NA_getType( PyObject *type)
+static int 
+isBufferWriteable (PyObject *buffobj) 
 {
-	PyObject *typeobj = NULL;
-	if (deferred_libnumarray_init() < 0) goto _exit;
-	if (!type) goto _exit;
-	if (PyObject_IsInstance(type, pNumericTypeClass)) {
-		Py_INCREF(type);
-		typeobj = type;
-		goto _exit;
-	}
-	if ((typeobj = PyDict_GetItem(pNumericTypesTDict, type))) {
-		Py_INCREF(typeobj);
-	} else {
-	       PyErr_Format(PyExc_ValueError, "NA_getType: unknown type.");
-	}
-  _exit:
-	return typeobj;
+  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;
 }
 
-/* Look up the NumarrayType which corresponds to typename */
+/**********************************************************************/
 
 static int 
-NA_nameToTypeNo(char *typename)
+getReadBufferDataPtr(PyObject *buffobj, void **buff) 
 {
-	int i;
-	for(i=0; i<ELEM(NumarrayTypeNameMap); i++)
-		if (!strcmp(typename, NumarrayTypeNameMap[i].name))
-			return NumarrayTypeNameMap[i].typeno;
-	return -1;
+  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;
 }
 
-/* Convert NumarrayType 'typeno' into the string of the type's name. */
+/**********************************************************************/
 
-static char *
-NA_typeNoToName(int typeno) 
+static int 
+getBufferSize(PyObject *buffobj) 
 {
-	int i;
-	PyObject *typeObj;
-	int typeno2;
+  int segcount, size=0;
+  PyObject *buff2;
+  if ((buff2 = getBuffer(buffobj))) 
+    {
+      segcount = buff2->ob_type->tp_as_buffer->bf_getsegcount(buff2, 
+								&size);
+      Py_DECREF(buff2);
+    }
+  else
+    size = -1;
+  return size;
+}
 
-	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);
+static double numarray_zero = 0.0;
 
-	return NA_typeNoToName(typeno2);
+static double raiseDivByZero(void)
+{
+  return 1.0/numarray_zero;
 }
 
-static PyObject *
-NA_typeNoToTypeObject(int typeno)
+static double raiseNegDivByZero(void)
 {
-	PyObject *o;
-	o = getTypeObject(typeno);	
-	if (o) Py_INCREF(o);
-	return o;
+  return -1.0/numarray_zero;
 }
 
-static int 
-NA_typeObjectToTypeNo(PyObject *typeObj)
+static double num_log(double x)
 {
-	int i;
-	if (deferred_libnumarray_init() < 0) return -1;
-	for(i=0; i<ELEM(pNumType); i++)
-		if (pNumType[i] == typeObj)
-			break;
-	if (i == ELEM(pNumType)) i = -1;
-	return i;
+   if (x == 0.0)
+       return raiseNegDivByZero();
+   else
+      return log(x);
 }
 
-static NumarrayType
-_scipy_typekind_to_typeNo(char typekind, int itemsize)
+static double num_log10(double x)
 {
-	int i;
-	for(i=0; i<ELEM(scipy_descriptors); i++) {
-		scipy_typestr *ts = &scipy_descriptors[i];
-		if ((typekind == ts->suffix[0]) && 
-		    (itemsize == ts->itemsize))
-			return i;
-	}
-	PyErr_Format(PyExc_TypeError,
-		     "Unknown __array_struct__ typekind");
-	return  -1;
+   if (x == 0.0)
+     return raiseNegDivByZero();
+   else
+     return log10(x);
 }
 
-static int
-NA_scipy_typestr(NumarrayType t, int byteorder, char *typestr)
+static double num_pow(double x, double y)
 {
-	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;
+   int z = (int) y;
+   if ((x < 0.0) && (y != z))
+     return raiseDivByZero();
+   else
+     return pow(x, y);
 }
 
-static long 
-NA_isIntegerSequence(PyObject *sequence)
+/* Inverse hyperbolic trig functions from Numeric */
+static double num_acosh(double x)
 {
-	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;
+    return log(x + sqrt((x-1.0)*(x+1.0)));
 }
 
-static PyObject *
-NA_intTupleFromMaybeLongs(int len, maybelong *Longs)
+static double num_asinh(double xx)
 {
-	int i;
-	PyObject *intTuple = PyTuple_New(len);
-	if (!intTuple) goto _exit;
-	for(i=0; i<len; i++) {
-		PyObject *o = PyInt_FromLong(Longs[i]);
-		if (!o) {
-			Py_DECREF(intTuple);
-			intTuple = NULL;
-			goto _exit;
-		}
-		PyTuple_SET_ITEM(intTuple, i, o);
-	}
-  _exit:
-	return intTuple;
+    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 long
-NA_maybeLongsFromIntTuple(int len, maybelong *arr, PyObject *sequence)
+static double num_atanh(double x)
 {
-	long i, size = -1;
-	if (!PySequence_Check(sequence)) {
-		PyErr_Format(PyExc_TypeError, 
-			     "NA_maybeLongsFromIntTuple: must be a sequence of integers.");
-		goto _exit;
-	}
-	size = PySequence_Length(sequence);
-	if (size < 0) {
-		PyErr_Format(PyExc_RuntimeError, 
-			     "NA_maybeLongsFromIntTuple: error getting sequence length.");
-		size = -1;
-		goto _exit;
-	}
-	if (size > len) {
-		PyErr_Format(PyExc_ValueError, 
-			     "NA_maybeLongsFromIntTuple: sequence is too long");
-		size = -1;
-		goto _exit;
-	}
-	for(i=0; i<size; i++) {
-		PyObject *o = PySequence_GetItem(sequence, i);
-		long value;
-		if (!o || !(PyInt_Check(o) || PyLong_Check(o))) {
-			PyErr_Format(PyExc_TypeError, 
-				     "NA_maybeLongsFromIntTuple: non-integer in sequence.");
-			Py_XDECREF(o);
-			size = -1;
-			goto _exit;
-		}
-		arr[i] = value = PyInt_AsLong(o);
-		if (arr[i] != value) {
-			PyErr_Format(PyExc_ValueError, 
-				     "NA_maybeLongsFromIntTuple: integer value too large: %ld",
-				     value);
-			size = -1;
-			goto _exit;
-		}
-		if (PyErr_Occurred()) {
-			Py_DECREF(o);
-			size = -1;
-			goto _exit;
-		}
-		Py_DECREF(o);
-	}
-  _exit:
-	return size;
+    return 0.5*log((1.0+x)/(1.0-x));
 }
 
-static int
-NA_intTupleProduct(PyObject  *shape, long *prod)
+/* NUM_CROUND (in numcomplex.h) also calls num_round */
+static double num_round(double x)
 {
-	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;
+    return (x >= 0) ? floor(x+0.5) : ceil(x-0.5);
 }
 
-/* NA_updateDataPtr updates the "working" data buffer pointer from the array's
-buffer object.  Since objects which meet the buffer API can potentially
-relocate their data as a result of executing arbitrary Python code
-(i.e. array.resize), NA_updateDataPtr needs to be called each time control
-flow returns to C, prior to accessing and numarray data.
-
-_data  points to an object which must meet the buffer API
-data   points to the array data gotten from _data via the buffer API
 
+/* 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 PyArrayObject *
-NA_updateDataPtr(PyArrayObject *me)
-{
-	if (!me) return me;
-
-	if (me->_data != Py_None) {
-
-		if (getReadBufferDataPtr (me->_data, 
-					  (void **) &me->data) < 0) {
-			return (PyArrayObject *) PyErr_Format(
-				_Error, 
-				"NA_updateDataPtr: error getting read buffer data ptr");
-		}
-		if (isBufferWriteable( me->_data ))
-			me->flags |= WRITABLE;
-		else
-			me->flags &= ~WRITABLE;
-	} else {
-		me->data = NULL;
-	}
-
-	me->data += me->byteoffset;
-
-	return me;
+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;
 }
 
-/* Count the number of elements in a 1D static array. */
-#define ELEM(x) (sizeof(x)/sizeof(x[0]))
-
-static int 
-NA_ByteOrder(void)
-{
-	unsigned long byteorder_test;
-	byteorder_test = 1;
-	if (*((char *) &byteorder_test))
-		return NUM_LITTLE_ENDIAN;
+/* Likewise for Integer overflows */
+#if defined(linux)
+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 NUM_BIG_ENDIAN;
+		return 1;
 }
+#endif
 
-/* Create a new numarray specifying all attribute values and using an object which
-   meets the buffer API to store the array data.
-*/
-static void
-_stridesFromShape(PyArrayObject *self)
+static int umult64_overflow(UInt64 a, UInt64 b)
 {
-    int i;
-    if (self->nd > 0) {
-	for(i=0; i<self->nd; i++)
-	    self->strides[i] = self->bytestride;
-	for(i=self->nd-2; i>=0; i--)
-	    self->strides[i] = 
-		self->strides[i+1]*self->dimensions[i+1];
-	self->nstrides = self->nd;
-    } else 
-	self->nstrides = 0;
-}
+        UInt64 ah, al, bh, bl, w, x, y, z;
 
- 
-static PyArrayObject *
-NA_NewAllFromBuffer(int ndim, maybelong *shape, NumarrayType type,
-		    PyObject *bufferObject, maybelong byteoffset, maybelong bytestride,
-		    int byteorder, int aligned, int writeable)
-{
-	PyObject *typeObject;
-	PyArrayObject *self = NULL;
-	long i;
+	ah = (a >> 32);
+	al = (a & 0xFFFFFFFFL);
+	bh = (b >> 32);
+	bl = (b & 0xFFFFFFFFL);
 
-	if (deferred_libnumarray_init() < 0) goto _fail;
+        /* 128-bit product:  z*2**64 + (x+y)*2**32 + w  */
+	w = al*bl;
+	x = bh*al;
+	y = ah*bl;
+	z = ah*bh;
 
-	if (type == tAny)
-		type = tDefault;
-	
-	if (ndim > MAXDIM) goto _fail;
+	/* *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;
 
-	{
-		PyTypeObject *typ = (PyTypeObject *) pNumArrayClass;
-		self = (PyArrayObject *) typ->tp_new(
-			typ, pEmptyTuple, pEmptyDict);
-		if (!self) goto _fail;
-	}
+        /* Convert to non-negative quantities */
+	if (a0 < 0) { a = -a0; } else { a = a0; }
+	if (b0 < 0) { b = -b0; } else { b = b0; }
 
-	typeObject = getTypeObject(type);
-        if (!typeObject) {
-		setTypeException(type);
-		goto _fail;
-	}
-	if (!(self->descr = NA_DescrFromType(type))) {
-		goto _fail;
-	}
+	ah = (a >> 32);
+	al = (a & 0xFFFFFFFFL);
+	bh = (b >> 32);
+	bl = (b & 0xFFFFFFFFL);
 
-	self->nd = self->nstrides = ndim;
-	for(i=0; i<ndim; i++) {
-		self->dimensions[i] = shape[i];
-	}
-	if (bytestride == 0)
-		self->bytestride = self->descr->elsize;
-	else
-		self->bytestride = bytestride;
-	_stridesFromShape(self);
+	w = al*bl;
+	x = bh*al;
+	y = ah*bl;
+	z = ah*bh;
 
-	self->byteoffset = byteoffset;
-	self->byteorder = byteorder;
-	self->itemsize = self->descr->elsize;
+        /* 
+         UInt64 c = ((x + y)<<32) + w;
+	 if ((a0 < 0) ^ (b0 < 0))
+	    *c = -c;
+	 else 
+	    *c = c
+	*/
 
-	Py_XDECREF(self->_data);
-	if ((bufferObject == Py_None) || (bufferObject == NULL)) {
-		long size = self->descr->elsize;
-		for(i=0; i<self->nd; i++) {
-			size *= self->dimensions[i];
-		}
-		self->_data = PyObject_CallFunction(
-			pNewMemoryFunc, "(l)", size);
-		if (!self->_data) goto _fail;
-	} else {
-		self->_data = bufferObject;
-		Py_INCREF(self->_data);
-	}
+	return z || (x>>31) || (y>>31) ||
+               (((x & 0xFFFFFFFFL) + (y & 0xFFFFFFFFL) + (w >> 32)) >> 31);
+}
 
-	if (!NA_updateDataPtr(self))
-		goto _fail;
-	NA_updateStatus(self);
-	return self;
 
-  _fail:
-	Py_XDECREF(self);
-	return NULL;
+static void
+NA_Done(void)
+{
+        return;
 }
 
-/* Create a new numarray specifying all attribute values but with a C-array as buffer
-   which will be copied to a Python buffer. 
-*/
 static PyArrayObject *
 NA_NewAll(int ndim, maybelong *shape, NumarrayType type, 
 	  void *buffer, maybelong byteoffset, maybelong bytestride, 
@@ -821,93 +373,29 @@ NA_NewAllStrides(int ndim, maybelong *shape, maybelong *strides,
 	return result;
 }
 
-static PyArrayObject *
-NA_FromDimsStridesDescrAndData(int nd, maybelong *d, maybelong *s, PyArray_Descr *descr, char *data) 
-{
-	maybelong i, nelements, breadth, bsize, boffset, dimensions[MAXDIM], strides[MAXDIM];
-	PyArrayObject *a;
-	PyObject *buf;
-	
-	if (!descr) return NULL;
-	
-	if (nd < 0) {
-		PyErr_SetString(PyExc_ValueError, 
-				"number of dimensions must be >= 0");
-		return NULL;
-	}
-	
-	if (nd > MAXDIM) {
-		return (PyArrayObject *) PyErr_Format(PyExc_ValueError,
-				    "too many dimensions: %d", nd);
-	}
-
-	if (!s) { /* no strides specified so assume minimal contiguous array 
-		     and compute */
-		if (nd) {
-			for(i=0; i<nd; i++)
-				strides[i] = descr->elsize;
-			for(i=nd-2; i>=0; i--)
-				strides[i] = strides[i+1]*d[i+1];
-		}
-	} else {
-		for(i=0; i<nd; i++)
-			strides[i] = s[i];
-	}
-
-	bsize = descr->elsize;   /* find buffer size implied by array 
-				    dimensions and strides */
-	boffset = 0;
-	for(i=0; i<nd; i++) {
-		breadth = llabs(strides[i]) * d[i];
-		bsize = MAX(bsize, breadth);
-		if (strides[i] < 0)
-			boffset += llabs(strides[i]) * (d[i]-1);
-	}
-	
-	nelements = 1;
-	for(i=0; i<nd; i++) {
-		dimensions[i] = d[i];
-		nelements *= d[i];
-	}
-	
-	if (data) {
-		buf = PyBuffer_FromReadWriteMemory(data-boffset, bsize);
-		if (!buf) return NULL;
-	} else {
-		buf = Py_None;
-	}
-	
-	a = NA_NewAllFromBuffer( nd, dimensions, descr->type_num, buf, 
-				 boffset, descr->elsize, NA_ByteOrder(), 1, 1);
-	if (!a) return NULL;
-	
-	for(i=0; i<nd; i++)
-		a->strides[i] = strides[i];
-
-	if (!data && !s) {
-		memset(a->data, 0, bsize);
-	}
-
-	NA_updateStatus(a);
-
-	return a;
-}
 
 static PyArrayObject *
-NA_FromDimsTypeAndData(int nd, maybelong *d, int type, char *data) 
+NA_New(void *buffer, NumarrayType type, int ndim, ...)
 {
-    PyArray_Descr *descr = NA_DescrFromType(type);
-    return NA_FromDimsStridesDescrAndData(nd, d, NULL, descr, data);
+	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_FromDimsStridesTypeAndData(int nd, maybelong *shape, maybelong *strides, 
-			     int type, char *data) 
+NA_Empty(int ndim, maybelong *shape, NumarrayType type)
 {
-    PyArray_Descr *descr = NA_DescrFromType(type);
-    return NA_FromDimsStridesDescrAndData(nd, shape, strides, descr, data);
+	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.
@@ -932,548 +420,7 @@ NA_NewArray(void *buffer, NumarrayType type, int ndim, ...)
 	return NA_vNewArray(buffer, type, ndim, shape);
 }
 
-/* Original deprecated versions of new array and empty array */
-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);
-}
-
-
-/* getArray creates a new array of type 't' from the given array 'a'
-using the specified 'method', probably 'new' or 'astype'. */
-
-static PyArrayObject *
-getArray(PyArrayObject *a, NumarrayType t, char *method)
-{
-	char *name;
-
-	if (deferred_libnumarray_init() < 0) return NULL;
-
-	if (t == tAny)
-		t = a->descr->type_num;
-	name = NA_typeNoToName(t);
-	if (!name) return (PyArrayObject *) setTypeException(t);
-	return (PyArrayObject *) 
-		PyObject_CallMethod((PyObject *) a, method, "s", name);
-}
-
-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;
-}
-
-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;
-}
-
-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;
-}
-
-/* sequenceAsArray converts a python sequence (list or tuple)
-into an array of the specified type and returns it.
-*/
-static PyArrayObject*
-sequenceAsArray(PyObject *s, NumarrayType *t)
-{
-	maybelong shape[MAXDIM];
-	int dims = getShape(s, shape, 0);
-	PyArrayObject *array;
-	
-	if (dims < 0) return NULL;
-	
-	if (*t == tAny) {
-		*t = NA_NumarrayType(s);
-	}
-	
-	if (!(array = NA_vNewArray(NULL, *t, dims, shape))) 
-		return NULL;
-	
-	if (setArrayFromSequence(array, s, 0, 0) < 0) {
-		return (PyArrayObject *) PyErr_Format(
-			_Error, 
-			"sequenceAsArray: can't convert sequence to array");
-	}
-	return array;
-}
-
-/* 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;
-}
-
-/* NA_InputArray is the main input conversion routine.  NA_InputArray
- converts input array 'a' as necessary to NumarrayType 't' also guaranteeing
- that either 'a' or the converted result is contigous, aligned, and not
- byteswapped. NA_InputArray returns a pointer to a Numarray for which
- C_array is 1, and fills in 'ainfo' with the array information of the return
- value.  The return value provides a means to later deallocate any temporary
- array created by NA_InputArray, while 'ainfo' provides direct access to the
- array's "metadata" from C.  Since the reference count of the input array 'a'
- is incremented when it is directly useable by C, the return value (either 'a'
- or a temporary) should always be passed to Py_XDECREF by the calller.  Note
- that at failed sequence conversion, getNumInfo, or getArray results in the
- value NULL being returned.
-
-1. 'a' is already c-usesable.
-2. 'a' is an array, but needs conversion to be c-useable.
-3. 'a' is a numeric sequence, not an array.
-4. 'a' is a numeric scalar, not an array.
-
-The contents of the resulting array are either 'a' or 'a.astype(t)'.
-The return value should always be DECREF'ed by the caller.
-
-requires is a bitmask specifying a set of requirements on the converted array.
-*/
-
-
-static PyArrayObject *
-NA_FromArrayStruct(PyObject *obj)
-{
-	PyArrayInterface *arrayif;
-	maybelong i, shape[MAXDIM], strides[MAXDIM];
-	NumarrayType t;
-	PyObject *cobj;
-	PyArrayObject *a;
-
-	cobj = PyObject_GetAttrString(obj, "__array_struct__"); /* does Py_INCREF */
-	if (!cobj) goto _fail;
-
-	if (!PyCObject_Check(cobj)) {
-		PyErr_Format(
-			PyExc_TypeError, 
-			"__array_struct__ returned non-CObject.");
-		goto _fail;
-	}
-
-	arrayif = PyCObject_AsVoidPtr(cobj);
-	if (arrayif->nd > MAXDIM) {
-		PyErr_Format( PyExc_ValueError, 
-			      "__array_struct__ too many dimensions: %d", 
-			      arrayif->nd);
-		goto _fail;
-	}
-
-	for(i=0; i<arrayif->nd; i++) {
-		shape[i] = arrayif->shape[i];
-		strides[i] = arrayif->strides[i];
-	}
-
-	t = _scipy_typekind_to_typeNo( arrayif->typekind, arrayif->itemsize );
-	if (t < 0) goto _fail;
-
-	a = NA_FromDimsStridesTypeAndData(arrayif->nd, shape, strides, t, arrayif->data);
-	if (!a) goto _fail;
-
-	a->base = cobj;
-
-	return a;
-
-  _fail:
-	Py_XDECREF(cobj);
-	return NULL;
-}
-
-static PyArrayObject*
-NA_InputArray(PyObject *a, NumarrayType t, int requires)
-{
-	PyArrayObject *wa = NULL;
-	if (NA_isPythonScalar(a)) {
-		if (t == tAny) 
-			t = NA_NumarrayType(a);
-		if (t < 0) goto _exit;
-		wa = NA_vNewArray( NULL, t, 0, NULL);
-		if (!wa) goto _exit;
-		if (NA_setFromPythonScalar(wa, 0, a) < 0) {
-			Py_DECREF(wa);
-			wa = NULL;
-		}
-		goto _exit;
-	} else if (NA_NumArrayCheck(a)) {
-		wa = (PyArrayObject *) a;
-		Py_INCREF(a);
-	} else if (PyObject_HasAttrString(a, "__array_struct__")) {
-		wa = NA_FromArrayStruct(a);
-	} else if (PyObject_HasAttrString(a, "__array_typestr__")) {
-		wa = (PyArrayObject *) PyObject_CallFunction( pNumArrayArrayFunc, "(O)", a );
-	} else {
-		wa = sequenceAsArray(a, &t);
-	}
-	if (!wa) goto _exit;
-	if (!satisfies(wa, requires, t)) {
-		PyArrayObject *wa2 = getArray(wa, t, "astype");
-		Py_DECREF(wa);
-		wa = wa2;
-	}
-	NA_updateDataPtr(wa);
-  _exit:
-	return wa;
-}
-
-/* NA_OutputArray creates a C-usable temporary array similar to 'a' but of
-type 't' as necessary.  If 'a' is already C-useable and of type 't', then 'a'
-is returned.  In either case, 'ainfo' is filled in with the array information
-for the return value.
-
-The contents of the resulting array are undefined, assumed to be filled in by
-the caller.
-*/
-static PyArrayObject *
-NA_OutputArray(PyObject *a0, NumarrayType t, int requires)
-{
-	PyArrayObject *a = (PyArrayObject *) a0;
-
-	if (!NA_NumArrayCheck(a0)  || !PyArray_ISWRITABLE(a)) {
-		PyErr_Format(PyExc_TypeError, 
-			     "NA_OutputArray: only writable NumArrays work for output.");
-		a = NULL;
-		goto _exit;
-	}
-
-	if (satisfies(a, requires, t)) {
-		Py_INCREF(a0);
-		NA_updateDataPtr(a);
-		goto _exit;
-	} else {
-		PyArrayObject *shadow = getArray(a, t, "new");
-		if (shadow) {
-			Py_INCREF(a0);
-			shadow->_shadows = a0;
-		}
-		a = shadow;
-	}
-  _exit:
-	return a;
-}
-
-/* NA_OptionalOutputArray works like NA_ShadowOutput, 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)) {
-		PyArrayObject *rval;
-		rval = getArray(master, t, "new");
-		return rval;
-	} else {
-		return NA_OutputArray(optional, t, requires);
-	}
-}
-
-/* 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);
-
-	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 numarray must be writable NumArrays.");
-		Py_DECREF(shadow);
-		shadow = NULL;
-		goto _exit;
-	}
 
-	if ((shadow != (PyArrayObject *) a) && NA_NumArrayCheck(a)) {
-		Py_INCREF(a);
-		shadow->_shadows = a;
-	}
-  _exit:
-	return shadow;
-}
-
-/* NA_ReturnOutput handles returning a possibly unspecified output array.  If
-the array 'out' was specified on the original call to the Python wrapper
-function, then the contents of any 'shadow' array are copied into 'out' as
-required.  The function then returns Py_None.  If no output was specified in
-the original call, then the 'shadow' array *becomes* the output and is
-returned.  This results in extension functions which return Py_None when you
-specify an output array, and return an array value otherwise.  These functions
-also correctly handle data typing, alignment, byteswapping, and contiguity
-issues. 
-*/
 static PyObject*
 NA_ReturnOutput(PyObject *out, PyArrayObject *shadow)
 {
@@ -1491,596 +438,35 @@ NA_ReturnOutput(PyObject *out, PyArrayObject *shadow)
 	}
 }
 
-/* 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;
-}
-
-#define MakeChecker(name, classpointer)                    \
-static int                                                 \
-name##Exact(PyObject *op) {                                \
-        return ((PyObject *) op->ob_type) == classpointer; \
-}                                                          \
-static int                                                 \
-name(PyObject *op) {                                       \
-      int rval = -1;                                       \
-      if (deferred_libnumarray_init() < 0) goto _exit;     \
-      rval = PyObject_IsInstance(op, classpointer);        \
-  _exit:                                                   \
-      return rval;                                         \
-}
-
-MakeChecker(NA_NDArrayCheck, pNDArrayClass)
-MakeChecker(NA_NumArrayCheck, pNumArrayClass)
-MakeChecker(NA_OperatorCheck, pOperatorClass)
-MakeChecker(NA_ConverterCheck, pConverterClass)
-MakeChecker(NA_UfuncCheck, pUfuncClass)
-MakeChecker(NA_CfuncCheck, pCfuncClass)
-
-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 PyObject * 
-NA_Cast(PyArrayObject *a, int type)
-{
-	PyObject *rval = NULL;
-	if (deferred_libnumarray_init() < 0) 
-		goto _exit;
-	rval = (PyObject *) getArray(a, type, "astype");
-  _exit:
-	return rval;
-}
-
-static int
-NA_copyArray(PyArrayObject *to, const PyArrayObject *from)
-{
-	int rval = -1;
-	PyObject *result;
-	result = PyObject_CallMethod((PyObject *) to, 
-				     "_copyFrom","(O)", from);
-	if (!result) goto _exit;
-	Py_DECREF(result);
-	rval = 0;
-  _exit:
-	return rval;
-}
-
-static PyArrayObject *
-NA_copy(PyArrayObject *from)
+static long NA_getBufferPtrAndSize(PyObject *buffobj, int readonly, void **ptr)
 {
-	PyArrayObject * rval;
-	rval = (PyArrayObject *)
-		PyObject_CallMethod((PyObject *) from, "copy", NULL);
+	long rval;
+	if (readonly)
+		rval = getReadBufferDataPtr(buffobj, ptr);
+	else 
+		rval = getWriteBufferDataPtr(buffobj, ptr);
 	return rval;
 }
 
-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_getByteOffset(PyArrayObject *array, int nindices, maybelong *indices, 
-		 long *offset)
+static int NA_checkIo(char *name, 
+		      int wantIn, int wantOut, int gotIn, int gotOut)
 {
-	int i;
-
-        /* rank0 or _UBuffer */
-	if ((array->nd == 0) || (array->nstrides < 0)) {
-		*offset = array->byteoffset;
-		return 0;
-	}
-	
-	/* Check for indices/shape mismatch when not rank-0.
-	 */
-	if ((nindices > array->nd) && 
-	    !((nindices == 1) && (array->nd == 0))) {
-		PyErr_Format(PyExc_IndexError, "too many indices.");
+	if (wantIn != gotIn) {
+		PyErr_Format(_Error,
+		   "%s: wrong # of input buffers. Expected %d.  Got %d.", 
+		    name, wantIn, gotIn);
 		return -1;
 	}
-
-	*offset = array->byteoffset;
-	for(i=0; i<nindices; i++) {
-		long ix = indices[i];
-		long limit = i < array->nd ? array->dimensions[i] : 0;
-		if (ix < 0) ix += limit;
-		if (ix < 0 || ix >= limit) {
-			PyErr_Format(PyExc_IndexError, "Index out of range");
-			return -1;
-		}
-		*offset += ix*array->strides[i];
-	}
-	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");
+	if (wantOut != gotOut) {
+		PyErr_Format(_Error,
+		   "%s: wrong # of output buffers. Expected %d.  Got %d.", 
+		    name, wantOut, gotOut);
 		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;
-
-	NA_updateStatus(array);
-	
 	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;
-}
-
-
-static long
-_isaligned(PyArrayObject *self)
-{
-	long i, ptr, alignment, aligned = 1;
-
-	alignment = MAX(MIN(self->itemsize, MAX_ALIGN), 1);
-	ptr = (long) self->data;
-	aligned = (ptr % alignment) == 0;
-	for (i=0; i <self->nd; i++)
-		aligned &= ((self->strides[i] % alignment) == 0);
-	return aligned != 0;
-}
-
-static void
-NA_updateAlignment(PyArrayObject *self)
-{
-	if (_isaligned(self))
-		self->flags |= ALIGNED;
-	else
-		self->flags &= ~ALIGNED;
-
-}
-
-static long
-_is_contiguous(PyArrayObject *self, maybelong elements)
-{
-	long i, ndim, nstrides;
-
-	ndim = self->nd;
-	nstrides = self->nstrides;
-
-	/* rank-0 numarray are always contiguous */
-	if (ndim == 0) return 1;
-
-	/* zero-length arrays are also contiguous */
-	if  (elements == 0) return 1;
-
-	/* Strides must be in decreasing order. ndim >= 1 */
-	for(i=0; i<ndim-1; i++)
-		if (self->strides[i] != 
-		    self->strides[i+1]*self->dimensions[i+1])
-			return 0;
-
-	/* Broadcast numarray have 0 in some stride and are discontiguous */
-	for(i=0; i<nstrides-1; i++)
-		if (!self->strides[i])
-			return 0;
-
-	if ((self->strides[nstrides-1] == self->itemsize) &&
-	    (self->bytestride == self->itemsize))
-		return 1;
-
-	if ((self->strides[nstrides-1] == 0) && (nstrides > 1))
-		return 1;
-	
-	return 0;
-}
-
-static long
-_is_fortran_contiguous(PyArrayObject *self, maybelong elements)
-{
-	long i, sd;
-
-	/* rank-0 numarray are always fortran_contiguous */
-	if (self->nd == 0) return 1;
-
-	/* zero-length arrays are also fortran_contiguous */
-	if  (elements == 0) return 1;
-
-	sd = self->descr->elsize;
-	for (i=0;i<self->nd;++i) { /* fortran == increasing order */
-		if (self->dimensions[i] == 0) return 0;  /* broadcast array */
-		if (self->strides[i] != sd) return 0; 
-		sd *= self->dimensions[i];
-	}
-
-	return 1;
-}
-
-static void
-NA_updateContiguous(PyArrayObject *self)
-{
-	maybelong elements = NA_elements(self);
-	if (_is_contiguous(self, elements))
-		self->flags |= CONTIGUOUS;
-	else
-		self->flags &= ~CONTIGUOUS;
-	if (_is_fortran_contiguous(self, elements))
-		self->flags |= FORTRAN_CONTIGUOUS;
-	else
-		self->flags &= ~FORTRAN_CONTIGUOUS;
-}
-
-static int 
-_isbyteswapped(PyArrayObject *self)
-{
-	int syslittle = (NA_ByteOrder() == NUM_LITTLE_ENDIAN);
-	int selflittle = (self->byteorder == NUM_LITTLE_ENDIAN);
-	int byteswapped = (syslittle != selflittle);
-	return byteswapped;
-}
-
-static void
-NA_updateByteswap(PyArrayObject *self)
-{
-	if (!_isbyteswapped(self))
-		self->flags |= NOTSWAPPED;
-	else
-		self->flags &= ~NOTSWAPPED;
-}
-
-static void 
-NA_updateStatus(PyArrayObject *self)
-{
-	NA_updateAlignment(self);
-	NA_updateContiguous(self);
-	NA_updateByteswap(self);
-}
-
-static char *
-NA_getArrayData(PyArrayObject *obj)
-{
-	if (!NA_NDArrayCheck((PyObject *) obj)) {
-		PyErr_Format(PyExc_TypeError, 
-			     "expected an NDArray");
-	}
-	if (!NA_updateDataPtr(obj))
-		return NULL;
-	return obj->data;
-}
-
-/* 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(linux) || defined(darwin) || defined(__CYGWIN__)
-
-#if defined(__GLIBC__) || defined(darwin) || defined(__MINGW32__)
-#include <fenv.h>
-#elif defined(__CYGWIN__)
-#include <mingw/fenv.h>
-#endif
-
-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();
-}
-
-static int
-NA_checkAndReportFPErrors(char *name)
-{
-	int error = NA_checkFPErrors();
-	if (error) {
-		PyObject *ans;
-		char msg[128];
-		if (deferred_libnumarray_init() < 0) 
-			return -1;
-		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;
-}
-
-/**********************************************************************/
-/*  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) 
-{
-  int segcount, size=0;
-  PyObject *buff2;
-  if ((buff2 = getBuffer(buffobj))) 
-    {
-      segcount = buff2->ob_type->tp_as_buffer->bf_getsegcount(buff2, 
-								&size);
-      Py_DECREF(buff2);
-    }
-  else
-    size = -1;
-  return size;
-}
-
-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_checkOneCBuffer(char *name, long niter, 
 		     void *buffer, long bsize, size_t typesize)
 {
@@ -2101,23 +487,6 @@ static int NA_checkOneCBuffer(char *name, long niter,
 	return 0;
 }
 
-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_checkNCBuffers(char *name, int N, long niter, 
 			    void **buffers, long *bsizes,
@@ -2131,20 +500,6 @@ static int NA_checkNCBuffers(char *name, int N, long niter,
 	return 0;
 }
 
-
-#if 0
-static void
-_dump_hex(char *name, int ndata, maybelong *data)
-{
-	int i;
-	fprintf(stderr, name);
-	for(i=0; i<ndata; i++)
-		fprintf(stderr, "%08x ", data[i]);
-	fprintf(stderr, "\n");
-	fflush(stderr);
-}
-#endif
-
 static int NA_checkOneStriding(char *name, long dim, maybelong *shape,
 	       long offset, maybelong *stride, long buffersize, long itemsize, 
 	       int align)
@@ -2173,10 +528,6 @@ static int NA_checkOneStriding(char *name, long dim, maybelong *shape,
 	return -1;
       }
       if (omax + itemsize > buffersize) {
-#if 0
-	      _dump_hex("shape:", dim, shape);
-	      _dump_hex("strides:", dim, stride);
-#endif
 	PyErr_Format(_Error, 
            "%s: access beyond buffer. offset=%d buffersize=%d",
 		     name, (int) (omax+itemsize-1), (int) buffersize);
@@ -2193,153 +544,309 @@ static int NA_checkOneStriding(char *name, long dim, maybelong *shape,
   return 0;
 }
 
-Float64 NA_get_Float64(PyArrayObject *a, long offset)
+/* 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)
 {
-	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); 
+	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;
 	}
-	return 0; /* suppress warning */
 }
 
-void NA_set_Float64(PyArrayObject *a, long offset, Float64 v)
-{
-	Bool b;
+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);
+}
 
-	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);
+
+/* 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 tFloat64: 
-		NA_SETP(a, Float64, (NA_PTR(a)+offset), v);
+	case CFUNC_STRIDING:
+		return callStrideConvCFunc(self, argsTuple);
 		break;
-	case tComplex32: {
-		NA_SETP(a, Float32, (NA_PTR(a)+offset), v);
-		NA_SETP(a, Float32, (NA_PTR(a)+offset+sizeof(Float32)), 0);
+	case CFUNC_NSTRIDING:
+		return callStridingCFunc(self, argsTuple);
+	case CFUNC_FROM_PY_VALUE:
+		return NumTypeFromPyValue(self, argsTuple);
 		break;
-	}
-	case tComplex64: {
-		NA_SETP(a, Float64, (NA_PTR(a)+offset), v);
-		NA_SETP(a, Float64, (NA_PTR(a)+offset+sizeof(Float64)), 0);
+	case CFUNC_AS_PY_VALUE:
+		return NumTypeAsPyValue(self, argsTuple);
 		break;
-	}
 	default:
-		PyErr_Format( PyExc_TypeError, 
-			      "Unknown type %d in NA_set_Float64",
-			      a->descr->type_num ); 
-		PyErr_Print();
+		return PyErr_Format( _Error,
+		     "cfunc_call: Can't dispatch cfunc '%s' with type: %d.", 
+		     me->descr.name, me->descr.type);
 	}
 }
 
-static int 
-NA_overflow(PyArrayObject *a, Float64 v)
+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)
 {
-	if ((a->flags & CHECKOVERFLOW) == 0) return 0;
+    CfuncObject* cfunc;
+    
+    CfuncType.ob_type = &PyType_Type;  /* Should be done once at init.
+					  Do now since there is no init. */
 
-	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;
+    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 (type == tAny) descr = NULL;
+        else descr = PyArray_DescrFromType(type);
+        return 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;
-	#endif
-	case tFloat32: 
-		if ((v < -FLT_MAX) || (v > FLT_MAX)) goto _fail;
+	if (!PyArray_ISALIGNED(a) && (requirements & NUM_ALIGNED)) 
 		return 0;
-	case tFloat64: 
+	if (!PyArray_ISCONTIGUOUS(a) && (requirements & NUM_CONTIGUOUS))
 		return 0;
-	case tComplex32: 
-		if ((v < -FLT_MAX) || (v > FLT_MAX)) goto _fail;
+	if (!PyArray_ISWRITABLE(a) && (requirements & NUM_WRITABLE))
 		return 0;
-	case tComplex64: 
+	if (requirements & NUM_COPY)
 		return 0;
-	default:
-		PyErr_Format( PyExc_TypeError, 
-			      "Unknown type %d in NA_overflow",
-			      a->descr->type_num ); 
-		PyErr_Print();
-		return -1;
+	return type_ok;
+}
+
+
+static PyArrayObject *
+NA_OutputArray(PyObject *a, NumarrayType t, int requires)
+{
+        PyObject *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 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 |= UPDATEIFCOPY;
+        ret->base = a;
+        PyArray_FLAGS(a) &= ~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 | 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)) {
+		PyArrayObject *rval;
+		rval = PyArray_FromArray(master, descr, 0);
+		return rval;
+	} else {
+		return NA_OutputArray(optional, t, requires);
 	}
-  _fail:
-	PyErr_Format(PyExc_OverflowError, "value out of range for array");
-	return -1;
 }
 
 Complex64 NA_get_Complex64(PyArrayObject *a, long offset)
@@ -2491,6 +998,97 @@ long NA_get_offset(PyArrayObject *a, int N, ...)
 	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];
@@ -2917,7 +1515,488 @@ NA_set1D_Complex64(PyArrayObject *a, long offset, int cnt, Complex64*in)
 	return 0;
 }
 
-#if LP64
+
+/* 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 *
+setTypeException(int type)
+{
+	/* Check if it is a printable character */
+	if ((type >= 32) && (type <= 126)) 
+		PyErr_Format(_Error, 
+			     "Type object lookup returned"
+			     " NULL for type \'%c\'", type);
+	else
+		PyErr_Format(_Error,
+			     "Type object lookup returned"
+			     " NULL for type %d", type);
+	return NULL;
+}
+
+static PyObject *
+getTypeObject(NumarrayType type) 
+{
+        return 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 (SIZEOF_INTP == 8)
 #define PlatBigInt PyInt_FromLong
 #define PlatBigUInt PyLong_FromUnsignedLong
 #else
@@ -3082,15 +2161,30 @@ NA_setFromPythonScalar(PyArrayObject *a, long offset, PyObject *value)
 	}
 }
 
+
 static int
-NA_isPythonScalar(PyObject *o)
+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;
-	rval =  PyInt_Check(o) || 
-		PyLong_Check(o) || 
-		PyFloat_Check(o) || 
-		PyComplex_Check(o) ||
-		(PyString_Check(o) && (PyString_Size(o) == 1));
+	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;
 }
 
@@ -3104,26 +2198,42 @@ NA_elements(PyArrayObject  *a)
 	return n;
 }
 
-staticforward PyTypeObject CfuncType;
+static int 
+NA_typeObjectToTypeNo(PyObject *typeObj)
+{
+        PyArray_Descr *dtype;
+        int i;
+        if (PyArray_DescrConverter(typeObj, &dtype) == PY_FAIL) i=-1;
+        else i=dtype->type_num;
+	return i;
+}
 
-static void
-cfunc_dealloc(PyObject* self)
+static int
+NA_copyArray(PyArrayObject *to, const PyArrayObject *from)
 {
-	PyObject_Del(self);
+        PyArrayObject *from0 = from;
+        return PyArray_CopyInto(to, from0);
 }
 
+static PyArrayObject *
+NA_copy(PyArrayObject *from)
+{
+        return (PyArrayObject *)PyArray_NewCopy(from, 0);
+}
+
+
 static PyObject *
-cfunc_repr(PyObject *self)
+NA_getType( PyObject *type)
 {
-	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);
+	PyArray_Descr *typeobj = NULL;
+        if (!type && PyArray_DescrConverter(type, &typeobj) == PY_FAIL) {
+                PyErr_Format(PyExc_ValueError, "NA_getType: unknown type.");
+                typeobj = NULL;
+        }
+        return typeobj;
 }
 
+
 /* Call a standard "stride" function
 **
 ** Stride functions always take one input and one output array.
@@ -3239,453 +2349,402 @@ NA_callStrideConvCFuncCore(
 	}
 }
 
-static PyObject *
-callStrideConvCFunc(PyObject *self, PyObject *args) {
-    PyObject *inbuffObj, *outbuffObj, *shapeObj;
-    PyObject *inbstridesObj, *outbstridesObj;
-    CfuncObject *me = (CfuncObject *) self;
-    int  nshape, ninbstrides, noutbstrides, nargs;
-    maybelong shape[MAXDIM], inbstrides[MAXDIM], 
-	    outbstrides[MAXDIM], *outbstrides1 = outbstrides;
-    long inboffset, outboffset, nbytes=0;
-
-    nargs = PyObject_Length(args);
-    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);
-    }
+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];
+	} 
+}
 
-    nshape = NA_maybeLongsFromIntTuple(MAXDIM, shape, shapeObj);
-    if (nshape < 0) return NULL;
-    
-    ninbstrides = NA_maybeLongsFromIntTuple(MAXDIM, inbstrides, inbstridesObj);
-    if (ninbstrides < 0) return NULL;
+static int
+NA_OperatorCheck(PyObject *op) {
+        return 0;
+}
 
-    noutbstrides=  NA_maybeLongsFromIntTuple(MAXDIM, outbstrides, outbstridesObj);
-    if (noutbstrides < 0) return NULL;
+static int
+NA_ConverterCheck(PyObject *op) {
+        return 0;
+}
 
-    if (nshape && (nshape != ninbstrides)) {
-        return PyErr_Format(_Error,
-            "%s: Missmatch between input iteration and strides tuples",
-	    me->descr.name);
-    }
+static int
+NA_UfuncCheck(PyObject *op) {
+        return 0;
+}
 
-    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);
-    }
-    
-#if 0    /* reductions slow mode hack...  wrong place to do it. */
-    _dump_hex("shape: ", nshape, shape);
-    _dump_hex("instrides: ", ninbstrides, inbstrides);
-    _dump_hex("outstrides: ", noutbstrides, outbstrides);
-
-    if (ninbstrides != noutbstrides) {
-	    outbstrides1 = outbstrides + (noutbstrides - ninbstrides);
-	    noutbstrides = ninbstrides;
-    } 
-#endif
+static int 
+NA_CfuncCheck(PyObject *op) {
+        return PyObject_TypeCheck(ob, &CfuncType);
+}
 
-    return NA_callStrideConvCFuncCore(
-	    self, nshape, shape,
-	    inbuffObj,  inboffset,  ninbstrides, inbstrides,
-	    outbuffObj, outboffset, noutbstrides, outbstrides1, nbytes);
+static int
+NA_getByteOffset(PyArrayObject *array, int nindices, maybelong *indices, 
+		 long *offset)
+{
+        return 0;
 }
 
-static int 
-_NA_callStridingHelper(PyObject *aux, long dim, 
-		       long nnumarray, PyArrayObject *numarray[], char *data[],
-		       CFUNC_STRIDED_FUNC f)
+static int
+NA_swapAxes(PyArrayObject *array, int x, int y)
 {
-	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;
+	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;
 	}
-	return status;
+
+	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, UPDATE_ALL_FLAGS);	
+
+	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;
+}
 
-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;
-    }
+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;
 }
 
-/* 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.
-*/
+/* ignores bytestride */
+static PyArrayObject *
+NA_NewAllFromBuffer(int ndim, maybelong *shape, NumarrayType type,
+		    PyObject *bufferObject, maybelong byteoffset, maybelong bytestride,
+		    int byteorder, int aligned, int writeable)
+{
+	PyObject *typeObject;
+	PyArrayObject *self = NULL;
+        PyArray_Descr *dtype;
 
-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.");
+	if (type == tAny)
+		type = tDefault;
 
-    /* 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);
-    }
+        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 = PyArray_NewFromDescr(&PyArray_Type, dtype,
+                                            ndim, shape, NULL, NULL, 0, NULL);
+        }
+        else {
+		intp size = dtype->elsize;
+                int i;
+		for(i=0; i<self->nd; i++) {
+			size *= self->dimensions[i];
+		}
+                self = PyArray_FromBuffer(bufferObject, dtype, size, byteoffset);
+        }
 
-    funcptr = (CFUNCasPyValue) me->descr.fptr;
+        return self;
+}
 
-    /* Call function to build PyObject.  Bad parameters to this function
-       may render call meaningless, but "temp" guarantees that its safe.  */
-    return (*funcptr)((void *)(&temp));
+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;
 }
 
-/* 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 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 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;
+static void
+NA_updateAlignment(PyArrayObject *self)
+{
+        return PyArray_UpdateFlags(self, ALIGNED);
 }
 
-/* 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 void
+NA_updateContiguous(PyArrayObject *self)
+{
+        return PyArray_UpdateFlags(self, CONTIGUOUS | FORTRAN);
+}
 
-static PyObject *
-NA_callCUFuncCore(PyObject *self, 
-		  long niter, long ninargs, long noutargs,
-		  PyObject **BufferObj, long *offset)
+
+static void 
+NA_updateStatus(PyArrayObject *self)
 {
-	CfuncObject *me = (CfuncObject *) self;
-	char *buffers[MAXARGS];
-	long bsizes[MAXARGS];
-	long i, pnargs = ninargs + noutargs;
-	UFUNC ufuncptr;
+        return PyArray_UpdateFlags(self, UPDATE_ALL_FLAGS);
+}
 
-	if (pnargs > MAXARGS) 
-		return PyErr_Format(PyExc_RuntimeError, "NA_callCUFuncCore: too many parameters");
+static int
+NA_NumArrayCheckExact(PyObject *op) {
+        return (op->ob_type == &PyArray_Type);
+}
 
-	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 int
+NA_NDArrayCheckExact(PyObject *op) {
+        return (op->ob_type == &PyArray_Type);
 }
 
-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);
+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 NA_callCUFuncCore(self, niter, ninargs, noutargs, BufferObj, offset);
+        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;
+}
 
-/* 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 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
 */
-static PyObject *
-cfunc_call(PyObject *self, PyObject *argsTuple, PyObject *argsDict)
+
+/*  OSF/Alpha (Tru64)  ---------------------------------------------*/
+#if defined(__osf__) && defined(__alpha)
+
+static int 
+NA_checkFPErrors(void)
 {
-	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);
-	}
+	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;
 }
 
-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 */
-};
+/* MS Windows -----------------------------------------------------*/
+#elif defined(_MSC_VER)
 
+#include <float.h>
 
-/* 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)
+static int 
+NA_checkFPErrors(void)
 {
-    CfuncObject* cfunc;
-    
-    CfuncType.ob_type = &PyType_Type;  /* Should be done once at init.
-					  Do now since there is no init. */
+	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);
 
-    cfunc = PyObject_New(CfuncObject, &CfuncType);
-    
-    if (!cfunc) {
-	    return PyErr_Format(_Error,
-			       "NA_new_cfunc: failed creating '%s'",
-			       cfd->name);
-    }
 
-    cfunc->descr = *cfd;
+	return retstatus;
+}
 
-    return (PyObject*)cfunc;
+/* 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;
 }
 
-static int NA_add_cfunc(PyObject *dict, char *keystr, CfuncDescriptor *descr) 
+#elif defined(linux) || defined(darwin) || defined(__CYGWIN__)
+
+#if defined(__GLIBC__) || defined(darwin) || defined(__MINGW32__)
+#include <fenv.h>
+#elif defined(__CYGWIN__)
+#include <mingw/fenv.h>
+#endif
+
+static int 
+NA_checkFPErrors(void)
 {
-	PyObject *c = (PyObject *) NA_new_cfunc(descr);
-	if (!c) return -1;
-	return PyDict_SetItemString(dict, keystr, c);
+	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;
+
 }
 
 
@@ -3780,174 +2839,107 @@ NA_IeeeMask64( Float64 f, Int32 mask)
 	return (category & mask) != 0;	
 }
 
-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 double numarray_zero = 0.0;
-
-static double raiseDivByZero(void)
-{
-  return 1.0/numarray_zero;
-}
-
-static double raiseNegDivByZero(void)
+static int 
+_NA_callStridingHelper(PyObject *aux, long dim, 
+		       long nnumarray, PyArrayObject *numarray[], char *data[],
+		       CFUNC_STRIDED_FUNC f)
 {
-  return -1.0/numarray_zero;
+	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 double num_log(double x)
+static PyArrayObject *
+NA_FromDimsStridesDescrAndData(int nd, maybelong *d, maybelong *s, PyArray_Descr *descr, char *data) 
 {
-   if (x == 0.0)
-       return raiseNegDivByZero();
-   else
-      return log(x);
+        return PyArray_NewFromDescr(&PyArray_Type, descr, nd, d, s, data, 0, NULL);
 }
 
-static double num_log10(double x)
+static PyArrayObject *
+NA_FromDimsTypeAndData(int nd, maybelong *d, int type, char *data) 
 {
-   if (x == 0.0)
-     return raiseNegDivByZero();
-   else
-     return log10(x);
+    PyArray_Descr *descr = NA_DescrFromType(type);
+    return NA_FromDimsStridesDescrAndData(nd, d, NULL, descr, data);
 }
 
-static double num_pow(double x, double y)
+static PyArrayObject *
+NA_FromDimsStridesTypeAndData(int nd, maybelong *shape, maybelong *strides, 
+			     int type, char *data) 
 {
-   int z = (int) y;
-   if ((x < 0.0) && (y != z))
-     return raiseDivByZero();
-   else
-     return pow(x, y);
+    PyArray_Descr *descr = NA_DescrFromType(type);
+    return NA_FromDimsStridesDescrAndData(nd, shape, strides, descr, data);
 }
 
-/* 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)
+typedef struct 
 {
-    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));
-}
+	NumarrayType type_num;
+	char suffix[5];
+	int  itemsize;
+} scipy_typestr;
 
-static double num_atanh(double x)
-{
-    return 0.5*log((1.0+x)/(1.0-x));
-}
+static scipy_typestr scipy_descriptors[ ] = {
+	{ tAny,       ""}, 
+	
+	{ tBool,      "b1", 1},
 
-/* 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);
-}
+	{ tInt8,      "i1", 1},
+	{ tUInt8,     "u1", 1},
 
-/* 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                                               *
-*/
+	{ tInt16,     "i2", 2},
+	{ tUInt16,    "u2", 2},
 
-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;
-}
+	{ tInt32,     "i4", 4},
+	{ tUInt32,    "u4", 4},
 
-/* Likewise for Integer overflows */
-#if defined(linux)
-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
+	{ tInt64,     "i8", 8},
+	{ tUInt64,    "u8", 8},
 
-static int umult64_overflow(UInt64 a, UInt64 b)
-{
-        UInt64 ah, al, bh, bl, w, x, y, z;
+	{ tFloat32,   "f4", 4},
+	{ tFloat64,   "f8", 8},
 
-	ah = (a >> 32);
-	al = (a & 0xFFFFFFFFL);
-	bh = (b >> 32);
-	bl = (b & 0xFFFFFFFFL);
+	{ tComplex32, "c8", 8},
+	{ tComplex64, "c16", 16}
+};
 
-        /* 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
+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 int smult64_overflow(Int64 a0, Int64 b0)
+static PyArrayObject *
+NA_FromArrayStruct(PyObject *obj)
 {
-	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);
+        return (PyArrayObject *)PyArray_FromStructInterface(obj);
 }
 
 
-
 static PyObject *_Error;
 
 void *libnumarray_API[] = {
@@ -4090,11 +3082,6 @@ static PyMethodDef _libnumarrayMethods[] = {
 };
 #endif
 
-/* platform independent*/
-#ifdef MS_WIN32
-__declspec(dllexport)
-#endif
-
 /* boiler plate API init */
 PyMODINIT_FUNC init_capi(void)
 {
@@ -4118,6 +3105,9 @@ PyMODINIT_FUNC init_capi(void)
     }
     if (PyModule_AddObject(m, "__version__", 
                            PyString_FromString("0.9")) < 0) return;
+
+    if (_import_array() < 0) return;
+    deferred_libnumarray_init();
     return;
 }