/*
   Python Array Object -- Provide multidimensional arrays as a basic
   object type in python.
    
  Copyright (c) 1995, 1996, 1997 Jim Hugunin, hugunin@mit.edu
   See file COPYING for details.


	These arrays are primarily designed for supporting multidimensional,
	homogeneous arrays of basic C numeric types.  They also can support
	arrays of arbitrary Python Objects, if you are willing to sacrifice
	performance for heterogeneity.
*/

/* $Id: arrayobject.c,v 1.19 2000/12/06 19:54:57 teoliphant Exp $ */

#include "Python.h"
/*Silly trick to make dll library work right under Win32*/
#ifdef MS_WIN32
#undef DL_IMPORT
#define DL_IMPORT(RTYPE) __declspec(dllexport) RTYPE
#endif
#define _ARRAY_MODULE
#include "arrayobject.h"
#define _UFUNC_MODULE
#include "ufuncobject.h"

#define OBJECT(O) ((PyObject*)(O))

/* There are several places in the code where an array of dimensions is */
/* allocated statically.  This is the size of that static allocation.  I */
/* can't come up with any reasonable excuse for a larger array than this. */

#define MAX_DIMS 40

/* Helper Functions */
extern int _PyArray_multiply_list(int *l1, int n) {
	int s=1, i=0;
	while (i < n) s *= l1[i++];
	return s;
}
extern int _PyArray_compare_lists(int *l1, int *l2, int n) {
	int i;
	for(i=0;i<n;i++) {
		if (l1[i] != l2[i]) return 0;
	} 
	return 1;
}

/* These can be cleaned up */
#define SIZE(mp) (_PyArray_multiply_list((mp)->dimensions, (mp)->nd))
#define NBYTES(mp) ((mp)->descr->elsize * SIZE(mp))
/* Obviously this needs some work. */
#define ISCONTIGUOUS(m) ((m)->flags & CONTIGUOUS)

#define PyArray_CONTIGUOUS(m) (ISCONTIGUOUS(m) ? Py_INCREF(m), m : \
(PyArrayObject *)(PyArray_ContiguousFromObject((PyObject *)(m), \
(m)->descr->type_num, 0,0)))

int do_sliced_copy(char *dest, int *dest_strides, int *dest_dimensions,
		   int dest_nd, char *src, int *src_strides, 
		   int *src_dimensions, int src_nd, int elsize, 
		   int copies) {
	int i, j;
	
        if (src_nd == 0 && dest_nd == 0) {
	    if (elsize == sizeof(char)) {
		memset(dest, *src, copies);
            } else if (elsize == sizeof(short)) {
                for(j=copies; j; --j, dest += sizeof(short))
                    *(short*)dest = *(short*)src;
	    } else if (elsize == sizeof(int)) {
                for(j=copies; j; --j, dest += sizeof(int))
                    *(int*)dest = *(int*)src;
            } else if (elsize == sizeof(long)) {
                for(j=copies; j; --j, dest += sizeof(long))
                    *(long*)dest = *(long*)src;
            } else if (elsize == sizeof(float)) {
                for(j=copies; j; --j, dest += sizeof(float))
                    *(float*)dest = *(float*)src;
            } else if (elsize == sizeof(double)) {
                for(j=copies; j; --j, dest += sizeof(double))
                    *(double*)dest = *(double*)src;
            } else {
                for(j=copies; j; --j, dest += elsize)
                    memcpy(dest, src, elsize);
            }
            return 0;
   	}


	if (dest_nd > src_nd) {
		for(i=0; i<*dest_dimensions; i++, dest += *dest_strides) {
			if (do_sliced_copy(dest, dest_strides+1, 
					   dest_dimensions+1, dest_nd-1,
					   src, src_strides, 
					   src_dimensions, src_nd, 
					   elsize, copies) == -1) 
			  return -1;
		}
		return 0;
	}
	
	if (dest_nd == 1) {
		if (*dest_dimensions != *src_dimensions) {
			PyErr_SetString(PyExc_ValueError, 
			  "matrices are not aligned for copy");
			return -1;
		}
		for(i=0; i<*dest_dimensions; i++, src += *src_strides) {
			for(j=0; j<copies; j++) {
				memcpy(dest, src, elsize);
				dest += *dest_strides;
			}
		}
	} else {
		for(i=0; i<*dest_dimensions; i++, dest += *dest_strides, 
		      src += *src_strides) {
			if (do_sliced_copy(dest, dest_strides+1, 
					   dest_dimensions+1, dest_nd-1, 
					   src, src_strides+1, 
					   src_dimensions+1, src_nd-1, 
					   elsize, copies) == -1) 
			  return -1;
		}
	}
	return 0;
}

int optimize_slices(int **dest_strides, int **dest_dimensions, 
		    int *dest_nd, int **src_strides, 
		    int **src_dimensions, int *src_nd,
		    int *elsize, int *copies) {
          while (*src_nd > 0) {
	          if (((*dest_strides)[*dest_nd-1] == *elsize) && 
		      ((*src_strides)[*src_nd-1] == *elsize)) {
			if ((*dest_dimensions)[*dest_nd-1] != 
			    (*src_dimensions)[*src_nd-1]) {
				PyErr_SetString(PyExc_ValueError, 
				  "matrices are not aligned for copy");
				return -1;
			}
			*elsize *= (*dest_dimensions)[*dest_nd-1];
			*dest_nd-=1; *src_nd-=1;
		} else {
			break;
		}
	}
	if (*src_nd == 0) {
		while (*dest_nd > 0) {
			if (((*dest_strides)[*dest_nd-1] == *elsize)) {
				*copies *= (*dest_dimensions)[*dest_nd-1];
				*dest_nd-=1;
			} else {
				break;
			}
		}
	}
	return 0;
}

static char *contiguous_data(PyArrayObject *src) {
	int dest_strides[MAX_DIMS], *dest_strides_ptr;
	int *dest_dimensions=src->dimensions;
	int dest_nd=src->nd;
	int *src_strides = src->strides;
	int *src_dimensions=src->dimensions;
	int src_nd=src->nd;
	int elsize=src->descr->elsize;
	int copies=1;
	int ret, i;
	int stride=elsize;
	char *new_data;
	
	for(i=dest_nd-1; i>=0; i--) {
		dest_strides[i] = stride;
		stride *= dest_dimensions[i];
	}
	
	dest_strides_ptr = dest_strides;
	
	if (optimize_slices(&dest_strides_ptr, &dest_dimensions, &dest_nd,
			    &src_strides, &src_dimensions, &src_nd,
			    &elsize, &copies) == -1) 
		return NULL;
	
	new_data = (char *)malloc(stride);
	
	ret = do_sliced_copy(new_data, dest_strides_ptr, dest_dimensions, 
			     dest_nd, src->data, src_strides, 
			     src_dimensions, src_nd, elsize, copies);
	
	if (ret != -1) { return new_data; }
	else { free(new_data); return NULL; }
}


/* Used for arrays of python objects to increment the reference count of */
/* every python object in the array. */
extern int PyArray_INCREF(PyArrayObject *mp) {
	int i, n;
	PyObject **data, **data2;
	
	if (mp->descr->type_num != PyArray_OBJECT) return 0;
	
	if (ISCONTIGUOUS(mp)) {
		data = (PyObject **)mp->data;
	} else {
		if ((data = (PyObject **)contiguous_data(mp)) == NULL) 
			return -1;
	}
	
	n = SIZE(mp);
	data2 = data;
	for(i=0; i<n; i++, data++) Py_XINCREF(*data);
	
	if (!ISCONTIGUOUS(mp)) free(data2);
	
	return 0;
}

extern int PyArray_XDECREF(PyArrayObject *mp) {
	int i, n;
	PyObject **data, **data2;
	
	if (mp->descr->type_num != PyArray_OBJECT) return 0;
	
	if (ISCONTIGUOUS(mp)) {
		data = (PyObject **)mp->data;
	} else {
		if ((data = (PyObject **)contiguous_data(mp)) == NULL) 
			return -1;
	}
	
	n = SIZE(mp);
	data2 = data;    
	for(i=0; i<n; i++, data++) Py_XDECREF(*data);
	
	if (!ISCONTIGUOUS(mp)) free(data2);
	
	return 0;
}

/* Including a C file is admittedly ugly.  Look at the C file if you want to */
/* see something even uglier. This is computer generated code. */
#include "arraytypes.c"

char *index2ptr(PyArrayObject *mp, int i) {
	if (i==0 && (mp->nd == 0 || mp->dimensions[0] > 0)) 
		return mp->data;
	
	if (mp->nd>0 &&  i>0 && i < mp->dimensions[0]) {
		return mp->data+i*mp->strides[0];
	}
	PyErr_SetString(PyExc_IndexError,"index out of bounds");      
	return NULL;
}

extern int PyArray_Size(PyObject *op) {
	if (PyArray_Check(op)) {
		return SIZE((PyArrayObject *)op);
	} else {
		return 0;
	}
}

int PyArray_CopyArray(PyArrayObject *dest, PyArrayObject *src) {
	int *dest_strides=dest->strides;
	int *dest_dimensions=dest->dimensions;
	int dest_nd=dest->nd;
	int *src_strides = src->strides;
	int *src_dimensions=src->dimensions;
	int src_nd=src->nd;
	int elsize=src->descr->elsize;
	int copies=1;
	int ret;
	
	if (src->nd > dest->nd) {
		PyErr_SetString(PyExc_ValueError, 
				"array too large for destination");
		return -1;
	}
	if (dest->descr->type_num != src->descr->type_num) {
		PyErr_SetString(PyExc_ValueError, 
		      "can only copy from a array of the same type.");
		return -1;
	}
	
	if (optimize_slices(&dest_strides, &dest_dimensions, &dest_nd, 
			    &src_strides, &src_dimensions, &src_nd,
			    &elsize, &copies) == -1) 
		return -1;
	
	ret = do_sliced_copy(dest->data, dest_strides, dest_dimensions, 
			     dest_nd, src->data, src_strides, 
			     src_dimensions, src_nd, elsize, copies);

	if (ret != -1) { ret = PyArray_INCREF(dest); }
	return ret;
}

int PyArray_CopyObject(PyArrayObject *dest, PyObject *src_object) {
	PyArrayObject *src;
	PyObject *tmp;
	int ret, n_new, n_old;
	char *new_string;
	
	/* Special function added here to try and make arrays of strings
	   work out. */
	if ((dest->descr->type_num == PyArray_CHAR) && dest->nd > 0 
	    && PyString_Check(src_object)) {
		n_new = dest->dimensions[dest->nd-1];
		n_old = PyString_Size(src_object); 
		if (n_new > n_old) {
			new_string = (char *)malloc(n_new*sizeof(char));
			memcpy(new_string, 
			  PyString_AS_STRING((PyStringObject *)src_object),
			                    n_old);
			memset(new_string+n_old, ' ', n_new-n_old);
			tmp = PyString_FromStringAndSize(new_string, 
							 n_new);
			free(new_string);
			src_object = tmp;
		}
	}
	src = (PyArrayObject *)PyArray_FromObject(src_object,
		dest->descr->type_num, 0,
		dest->nd);
	if (src == NULL) return -1;
	
	ret = PyArray_CopyArray(dest, src);
	Py_DECREF(src);
	return ret;
}

/* This is the basic array allocation function. */
PyObject *PyArray_FromDimsAndDataAndDescr(int nd, int *d, 
					  PyArray_Descr *descr,
					  char *data) {
	PyArrayObject *self;
	int i,sd;
	int *dimensions, *strides;
	int flags=CONTIGUOUS | OWN_DIMENSIONS | OWN_STRIDES;
	
	dimensions = strides = NULL;
	
	if (nd < 0) {
		PyErr_SetString(PyExc_ValueError, 
				"number of dimensions must be >= 0");
		return NULL;
	}
	
	if (nd > 0) {
		if ((dimensions = (int *)malloc(nd*sizeof(int))) == NULL) {
			PyErr_SetString(PyExc_MemoryError, "can't allocate memory for array");
			goto fail;
		}
		if ((strides = (int *)malloc(nd*sizeof(int))) == NULL) {
			PyErr_SetString(PyExc_MemoryError, "can't allocate memory for array");
			goto fail;
		}
		memcpy(dimensions, d, sizeof(int)*nd);
	}
	
	sd = descr->elsize;
	for(i=nd-1;i>=0;i--) {
		if (flags & OWN_STRIDES) strides[i] = sd;
		if (dimensions[i] < 0) {
			PyErr_SetString(PyExc_ValueError, "negative dimensions are not allowed");
			goto fail;
		}
		/* 
		   This may waste some space, but it seems to be
		   (unsuprisingly) unhealthy to allow strides that are
		   longer than sd.
		   */
		sd *= dimensions[i] ? dimensions[i] : 1;
/*		sd *= dimensions[i]; */
	}
	
	/* Make sure we're alligned on ints. */
	sd += sizeof(int) - sd%sizeof(int); 
	
	if (data == NULL) {
		if ((data = (char *)malloc(sd)) == NULL) {
			PyErr_SetString(PyExc_MemoryError, "can't allocate memory for array");
			goto fail;
		}
		flags |= OWN_DATA;
	}

	if((self = PyObject_NEW(PyArrayObject, &PyArray_Type)) == NULL) goto fail;
	if (flags & OWN_DATA) memset(data, 0, sd);
	
	self->data=data;
	self->dimensions = dimensions;
	self->strides = strides;
	self->nd=nd;
	self->descr=descr;
	self->base = (PyObject *)NULL;
	self->flags = flags;
	
	return (PyObject*)self;
	
fail:
	
	if (flags & OWN_DATA) free(data);
	if (dimensions != NULL) free(dimensions);
	if (strides != NULL) free(strides);
	return NULL;    
}

PyObject *PyArray_FromDimsAndData(int nd, int *d, int type, char *data) {
	PyArray_Descr *descr;
	PyObject *op;
	char real_type;
		
	real_type = type;
	real_type = real_type & ~SAVESPACEBIT;
	if ((descr = PyArray_DescrFromType((int)real_type)) == NULL) return NULL;
	op = PyArray_FromDimsAndDataAndDescr(nd, d, descr, data);
	if (type & SAVESPACEBIT)
	        ((PyArrayObject *)op)->flags |= SAVESPACE;
	return op;
}

PyObject *PyArray_FromDims(int nd, int *d, int type) {
	PyArray_Descr *descr;
	PyObject *op;
        char real_type;

	real_type = type;
	real_type = real_type & ~SAVESPACEBIT; 
	if ((descr = PyArray_DescrFromType((int)real_type)) == NULL) return NULL;
	op = PyArray_FromDimsAndDataAndDescr(nd, d, descr, NULL);
	if (type & SAVESPACEBIT) 
	        ((PyArrayObject *)op)->flags |= SAVESPACE;
	return op;
}


extern PyObject *PyArray_Copy(PyArrayObject *m1) {
	PyArrayObject *ret = 
		(PyArrayObject *)PyArray_FromDims(m1->nd, m1->dimensions, m1->descr->type_num);
	
	if (PyArray_CopyArray(ret, m1) == -1) return NULL;
	
	return (PyObject *)ret;
}

static void array_dealloc(PyArrayObject *self) {
	if(self->base) Py_DECREF(self->base);
	
	if (self->flags & OWN_DATA) {
		PyArray_XDECREF(self);
		free(self->data); 
	}
	
	if (self->flags & OWN_DIMENSIONS && self->dimensions != NULL) free(self->dimensions); 
	if (self->flags & OWN_STRIDES && self->strides != NULL) free(self->strides); 
	
	PyMem_DEL(self);
}

/* Code to handle accessing Array objects as sequence objects */
static int array_length(PyArrayObject *self) {
	if (self->nd != 0) {
		return self->dimensions[0];
	} else {
		return 1; /* Because a[0] works on 0d arrays. */
	}
}

static PyObject *array_item(PyArrayObject *self, int i) {
	char *item;

	if ((item = index2ptr(self, i)) == NULL) return NULL;
	
	if(self->nd > 1) {
		PyArrayObject *r;
        if (( r= PyObject_NEW(PyArrayObject, &PyArray_Type)) == NULL) return NULL;

		r->nd = self->nd-1;
		r->dimensions = self->dimensions+1;
		r->strides = self->strides+1;
		r->descr = self->descr;
		r->data = item;
		r->base = (PyObject *)self;
		r->flags = self->flags & (CONTIGUOUS | SAVESPACE);
		Py_INCREF(self);
		return (PyObject*)r;
	} else {

	  /* I would like to do this, but it requires a fix to several places of code.
             fprintf(stderr,"Getting a Python scalar by indexing a rank-0 array is obsolete: use a.toscalar().\n");
	  */
                return self->descr->getitem(item);
	} 
}

extern PyObject * PyArray_Item(PyObject *op, int i) {
	if (PyArray_Check(op)) {
		return array_item((PyArrayObject *)op, i);
	} else {
		PyErr_SetString(PyExc_ValueError, "Not an array object");
		return NULL;
	}
}

extern PyObject *PyArray_Return(PyArrayObject *mp) {
	PyObject *op;
	
	if (PyErr_Occurred()) {
		if (mp != NULL) Py_DECREF(mp);
		return NULL;
	}
	if (mp->nd == 0 && (mp->descr->type_num == PyArray_LONG || mp->descr->type_num == PyArray_DOUBLE || mp->descr->type_num == PyArray_CDOUBLE || mp->descr->type_num == PyArray_OBJECT )) {
       	        op = mp->descr->getitem(mp->data);
		Py_DECREF(mp);
		return op;
	} else {
		return (PyObject *)mp;
	}
}

static PyObject * array_slice(PyArrayObject *self, int ilow, int ihigh) {
	PyArrayObject *r;
	int l;
	char *data;
	
	if (self->nd == 0) {
		PyErr_SetString(PyExc_ValueError, "can't slice a scalar");
		return NULL;
	}
	
	l=self->dimensions[0];
	if (ihigh < 0) ihigh += l;
	if (ilow  < 0) ilow  += l;
	if (ilow < 0) ilow = 0;
	else if (ilow > l) ilow = l;
	if (ihigh < 0) ihigh = 0;
	else if (ihigh > l) ihigh = l;
	if (ihigh < ilow) ihigh = ilow;

	if (ihigh != ilow) {
		data = index2ptr(self, ilow);
		if (data == NULL) return NULL;
	} else {
 		data = self->data;
	}

	self->dimensions[0] = ihigh-ilow;
	r = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(self->nd, self->dimensions, self->descr, data);
	self->dimensions[0] = l;
	
	if (!ISCONTIGUOUS(self)) r->flags &= ~CONTIGUOUS;
	if (self->flags & SAVESPACE) r->flags |= SAVESPACE;
	memcpy(r->strides, self->strides, sizeof(int)*self->nd);
	r->base = (PyObject *)self;
	Py_INCREF(self);
	
	return (PyObject *)r;
}

/* These will need some serious work when I feel like it. */
static int array_ass_item(PyArrayObject *self, int i, PyObject *v) {
	PyObject *c=NULL;
	PyArrayObject *tmp;
	char *item;
	int ret;
	
    if (v == NULL) {
        PyErr_SetString(PyExc_ValueError, "Can't delete array elements.");
        return -1;
    }

	if (i < 0) i = i+self->dimensions[0];

	if (self->nd > 1) {
		if((tmp = (PyArrayObject *)array_item(self, i)) == NULL) return -1;
		ret = PyArray_CopyObject(tmp, v);
		Py_DECREF(tmp);
		return ret;
	}
	
	if ((item = index2ptr(self, i)) == NULL) return -1;
	
	if(self->descr->type_num != PyArray_OBJECT && PyString_Check(v) && PyObject_Length(v) == 1) {
		char *s;
		if ((s=PyString_AsString(v)) == NULL) return -1;
		if(self->descr->type == 'c') {
			((char*)self->data)[i]=*s;
			return 0;
		}
		if(s) c=PyInt_FromLong((long)*s);
		if(c) v=c;
	}
	
	self->descr->setitem(v, item);
	
	if(c) Py_DECREF(c);
	if(PyErr_Occurred()) return -1;
	return 0;
}


static int array_ass_slice(PyArrayObject *self, int ilow, int ihigh, PyObject *v) {
	int ret;
	PyArrayObject *tmp;
	
        if (v == NULL) {
                PyErr_SetString(PyExc_ValueError, "Can't delete array elements.");
                return -1;
        }
	if ((tmp = (PyArrayObject *)array_slice(self, ilow, ihigh)) == NULL) return -1;  
	ret = PyArray_CopyObject(tmp, v);
	Py_DECREF(tmp);
	
	return ret;
}

/* -------------------------------------------------------------- */
static int slice_GetIndices(PySliceObject *r, int length, 
					   int *start, int *stop, int *step)
{
	if (r->step == Py_None) {
		*step = 1;
	} else {
		if (!PyInt_Check(r->step)) return -1;
		*step = PyInt_AsLong(r->step);
	}
	if (r->start == Py_None) {
		*start = *step < 0 ? length-1 : 0;
	} else {
		if (!PyInt_Check(r->start)) return -1;
		*start = PyInt_AsLong(r->start);
		if (*start < 0) *start += length;
	}
	if (r->stop == Py_None) {
		*stop = *step < 0 ? -1 : length;
	} else {
		if (!PyInt_Check(r->stop)) return -1;
		*stop = PyInt_AsLong(r->stop);
		if (*stop < 0) *stop += length;
	}

	if (*start > (length-1)) *start = length;
	if (*start < 0) *start = 0;
	if (*stop < -1) *stop = -1;
	else if (*stop > length) *stop = length;
	return 0;
}


static int get_slice(PyObject *op, int max, int *np, int *sp) {
	int start, stop, step;
	
	if (PySlice_Check(op)) {
		if (slice_GetIndices((PySliceObject *)op, max, 
								&start, &stop, &step) == -1) return -1;
		
		if (step != 0) {
			if (step < 0) *np = (stop-start+1+step)/step;
			else *np = (stop-start-1+step)/step;
		} else {
			if (stop == start) {
				*np = 0; step = 1;
			}
			else return -1;
		}
		if (*np < 0) *np = 0;
		*sp = step;
		return start;
	}  
	return -1;
}

#define PseudoIndex -1
#define RubberIndex -2
#define SingleIndex -3

static int parse_subindex(PyObject *op, int *step_size, int *n_steps, int max) {
	int i, tmp;
	
	if (op == Py_None) {
		*n_steps = PseudoIndex;
		return 0;
	}
	
	if (op == Py_Ellipsis) {
		*n_steps = RubberIndex;
		return 0;
	}
	
	if (PySlice_Check(op)) {
		if ((i = get_slice(op, max, n_steps, step_size)) >= 0) {
			return i;
		} else {
			PyErr_SetString(PyExc_IndexError, "invalid slice");
			return -1;
		}
	}
	
	if (PyInt_Check(op)) {
		*n_steps=SingleIndex;
		*step_size=0;
		tmp = PyInt_AsLong(op);
		if (tmp < 0) tmp += max;
		if (tmp >= max || tmp < 0) {
			PyErr_SetString(PyExc_IndexError, "invalid index");
			return -1;
		}
		return tmp;
	} 
	PyErr_SetString(PyExc_IndexError, 
		"each subindex must be either a slice, an integer, Ellipsis, or NewAxis");
	return -1;
}

static int parse_index(PyArrayObject *self, PyObject *op, 
		       int *dimensions, int *strides, int *offset_ptr) {
	int i, j, n;
	int nd_old, nd_new, start, offset, n_add, n_pseudo;
	int step_size, n_steps;
	PyObject *op1;
	int is_slice;
	
	if (PySlice_Check(op) || op == Py_Ellipsis) {
		n = 1;
		op1 = op;
		Py_INCREF(op);  
                /* this relies on the fact that n==1 for loop below */
		is_slice = 1;
	} else {
		if (!PySequence_Check(op)) {
			PyErr_SetString(PyExc_IndexError, 
                           "index must be either an int or a sequence");
			return -1;
		}
		n = PySequence_Length(op);
		is_slice = 0;
	}
	
	nd_old = nd_new = 0;
	
	offset = 0;
	for(i=0; i<n; i++) {
		if (!is_slice) {
		  if (!(op1=PySequence_GetItem(op, i))) {
		    PyErr_SetString(PyExc_IndexError, "invalid index");
		    return -1;
		  }
		}
		
		start = parse_subindex(op1, &step_size, &n_steps, 
			nd_old < self->nd ? self->dimensions[nd_old] : 0);
		Py_DECREF(op1);
		if (start == -1) break;
		
		if (n_steps == PseudoIndex) {
			dimensions[nd_new] = 1; strides[nd_new] = 0; nd_new++;
		} else {
			if (n_steps == RubberIndex) {
				for(j=i+1, n_pseudo=0; j<n; j++) {
					op1 = PySequence_GetItem(op, j);
					if (op1 == Py_None) n_pseudo++;
					Py_DECREF(op1);
				}
				n_add = self->nd-(n-i-n_pseudo-1+nd_old);
				if (n_add < 0) {
					PyErr_SetString(PyExc_IndexError, "too many indices");
					return -1;
				}
				for(j=0; j<n_add; j++) {
					dimensions[nd_new] = self->dimensions[nd_old];
					strides[nd_new] = self->strides[nd_old];
					nd_new++; nd_old++;
				}
			} else {
				if (nd_old >= self->nd) {
					PyErr_SetString(PyExc_IndexError, "too many indices");
					return -1;
				}
				offset += self->strides[nd_old]*start;
				nd_old++;
				if (n_steps != SingleIndex) {
					dimensions[nd_new] = n_steps;
					strides[nd_new] = step_size*self->strides[nd_old-1];
					nd_new++;
				}
			}
		}
	}
	if (i < n) return -1;
	n_add = self->nd-nd_old;
	for(j=0; j<n_add; j++) {
		dimensions[nd_new] = self->dimensions[nd_old];
		strides[nd_new] = self->strides[nd_old];
		nd_new++; nd_old++;
	}     
	*offset_ptr = offset;
	return nd_new;
}

/* Called when treating array object like a mapping -- called first from 
   Python when using a[object] */
static PyObject *array_subscript(PyArrayObject *self, PyObject *op) {
	int dimensions[MAX_DIMS], strides[MAX_DIMS];
	int nd, offset, i, elsize;
	PyArrayObject *other;
	
	if (PyInt_Check(op)) {
	        i = PyInt_AsLong(op);
		if (i < 0 && self->nd > 0) i = i+self->dimensions[0]; 
		return array_item(self, i);
	}
	
	if ((nd = parse_index(self, op, dimensions, strides, &offset)) 
	    == -1) {
		return NULL;
	}
	
	if ((other = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(nd, 
					dimensions,
					self->descr,
		                        self->data+offset)) == NULL) {
		return NULL;
	}
	memcpy(other->strides, strides, sizeof(int)*other->nd);
	other->base = (PyObject *)self;
	Py_INCREF(self);
	
	elsize=other->descr->elsize;
	/* Check to see if other is CONTIGUOUS:  see if strides match 
           dimensions */
	for (i=other->nd-1; i>=0; i--) {
		if (other->strides[i] == elsize) {
			elsize *= other->dimensions[i];
		} else {
			break;
		}
	}
	if (i >= 0) other->flags &= ~CONTIGUOUS; 

	/* Maintain SAVESPACE flag on selection */
	if (self->flags & SAVESPACE) other->flags |= SAVESPACE;
	
	return (PyObject *)other;
}

static PyObject *array_subscript_nice(PyArrayObject *self, PyObject *op) {
  	PyObject *ret;
        
	if ((ret = array_subscript(self, op)) == NULL) return NULL;
	if (PyArray_Check(ret)) return PyArray_Return((PyArrayObject *)ret);
	else return ret;
}


/* Another assignment hacked by using CopyObject.  */

static int array_ass_sub(PyArrayObject *self, PyObject *index, PyObject *op) {
	int ret;
	PyArrayObject *tmp;
	
	if (op == NULL) {
		PyErr_SetString(PyExc_ValueError, 
				"Can't delete array elements.");
		return -1;
	}
	
	if (PyInt_Check(index)) 
		return array_ass_item(self, PyInt_AsLong(index), op);


	
	if ((tmp = (PyArrayObject *)array_subscript(self, index)) == NULL)
		return -1; 
	ret = PyArray_CopyObject(tmp, op);
	Py_DECREF(tmp);
	
	return ret;
}

static PyMappingMethods array_as_mapping = {
	(inquiry)array_length,		/*mp_length*/
	(binaryfunc)array_subscript_nice,   /*mp_subscript*/
	(objobjargproc)array_ass_sub,	/*mp_ass_subscript*/
};

/*  These should be added
static PyBufferProcs array_as_buffer = {
        (getreadbufferproc)array_getreadbuf,
        (getwritebufferproc)array_getwritebuf,
        (getsegcountproc)array_getsegcount,
        (getcharbufferproc)array_getcharbuf,
};
*/



/* -------------------------------------------------------------- */

typedef struct {
	PyUFuncObject *add, 
		*subtract, 
		*multiply, 
		*divide, 
		*remainder, 
		*power, 
		*negative, 
		*absolute;
	PyUFuncObject *invert, 
		*left_shift, 
		*right_shift, 
		*bitwise_and, 
		*bitwise_xor,
		*bitwise_or;
} NumericOps;

static NumericOps n_ops;

#define GET(op) n_ops.op = (PyUFuncObject *)PyDict_GetItemString(dict, #op)

int PyArray_SetNumericOps(PyObject *dict) {
	GET(add);
	GET(subtract);
	GET(multiply);
	GET(divide);
	GET(remainder);
	GET(power);
	GET(negative);
	GET(absolute);
	GET(invert);
	GET(left_shift);
	GET(right_shift);
	GET(bitwise_and);
	GET(bitwise_or);
	GET(bitwise_xor);
	return 0;
}

static int array_coerce(PyArrayObject **pm, PyObject **pw) {
	PyObject *new_op;
	if ((new_op = PyArray_FromObject(*pw, PyArray_NOTYPE, 0, 0)) 
	    == NULL) 
		return -1;
	Py_INCREF(*pm);
	*pw = new_op;
	return 0;
	
/**Eliminate coercion at this stage.  Handled by umath now
	PyObject *new_op;
	char t1, t2;
	
	t1 = (*pm)->descr->type_num;
	t2 = PyArray_ObjectType(*pw, 0);
	
	if (PyArray_CanCastSafely(t2, t1)) { 
		if ((new_op = PyArray_FromObject(*pw, t1, 0, 0)) == NULL) return -1;
		Py_INCREF(*pm);
		*pw = new_op;
	} else {
		if (PyArray_CanCastSafely(t1, t2)) {
			*pm = (PyArrayObject *)PyArray_Cast(*pm, t2);
			if ((new_op = PyArray_FromObject(*pw, t2, 0, 0)) == NULL) return -1;
			*pw = new_op;
		} else {
			PyErr_SetString(PyExc_TypeError, "cannot perform this operation on these types");
			return -1;
		}
	}
	return 0;
***/
}

static PyObject *PyUFunc_BinaryFunction(PyUFuncObject *s, PyArrayObject *mp1, PyObject *mp2) {
	PyObject *arglist;
	PyArrayObject *mps[3];
	
	arglist = Py_BuildValue("(OO)", mp1, mp2);
	
	mps[0] = mps[1] = mps[2] = NULL;
	if (PyUFunc_GenericFunction(s, arglist, mps) == -1) {
		Py_XDECREF(mps[0]); Py_XDECREF(mps[1]); Py_XDECREF(mps[2]);
		return NULL;
	}
	
	Py_DECREF(mps[0]); Py_DECREF(mps[1]);
	Py_DECREF(arglist);
	return PyArray_Return(mps[2]);
}

static PyObject *PyUFunc_UnaryFunction(PyUFuncObject *s, PyArrayObject *mp1) {
	PyObject *arglist;
	PyArrayObject *mps[3];
	
	arglist = Py_BuildValue("(O)", mp1);
	
	mps[0] = mps[1] = NULL;
	if (PyUFunc_GenericFunction(s, arglist, mps) == -1) {
		Py_XDECREF(mps[0]); Py_XDECREF(mps[1]);
		return NULL;
	}
	
	Py_DECREF(mps[0]);
	Py_DECREF(arglist);
	return PyArray_Return(mps[1]);
}

static PyObject *array_add(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.add, m1, m2);
}
static PyObject *array_subtract(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.subtract, m1, m2);
}
static PyObject *array_multiply(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.multiply, m1, m2);
}
static PyObject *array_divide(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.divide, m1, m2);
}
static PyObject *array_remainder(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.remainder, m1, m2);
}
static PyObject *array_power(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.power, m1, m2);
}
static PyObject *array_negative(PyArrayObject *m1) { 
    return PyUFunc_UnaryFunction(n_ops.negative, m1);
}
static PyObject *array_absolute(PyArrayObject *m1) { 
    return PyUFunc_UnaryFunction(n_ops.absolute, m1);
}
static PyObject *array_invert(PyArrayObject *m1) { 
    return PyUFunc_UnaryFunction(n_ops.invert, m1);
}
static PyObject *array_left_shift(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.left_shift, m1, m2);
}
static PyObject *array_right_shift(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.right_shift, m1, m2);
}
static PyObject *array_bitwise_and(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.bitwise_and, m1, m2);
}
static PyObject *array_bitwise_or(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.bitwise_or, m1, m2);
}
static PyObject *array_bitwise_xor(PyArrayObject *m1, PyObject *m2) {
    return PyUFunc_BinaryFunction(n_ops.bitwise_xor, m1, m2);
}


static int array_nonzero(PyArrayObject *mp) {
	char *zero;
	PyArrayObject *self;
	char *data;
	int i, s, elsize;
	
	self = PyArray_CONTIGUOUS(mp);
	zero = self->descr->zero;

	s = SIZE(self);
	elsize = self->descr->elsize;
	data = self->data;
	for(i=0; i<s; i++, data+=elsize) {
		if (memcmp(zero, data, elsize) != 0) break;
	}
	
	Py_DECREF(self);
	
	return i!=s;
}

static PyObject *array_divmod(PyArrayObject *op1, PyObject *op2) {
	PyObject *divp, *modp;

	divp = array_divide(op1, op2);
	if (divp == NULL) return NULL;
	modp = array_remainder(op1, op2);
	if (modp == NULL) {
		Py_DECREF(divp);
		return NULL;
	}
	return Py_BuildValue("OO", divp, modp);
}


/* methods added by Travis Oliphant, Jan 2000 */
static PyObject *array_int(PyArrayObject *v) {        
        PyObject *pv, *pv2;
        if (v->nd != 0) {
	        PyErr_SetString(PyExc_TypeError, "Only rank-0 arrays can be converted to Python scalars.");
		return NULL;
	}
	pv = v->descr->getitem(v->data);
	if (pv == NULL) return NULL;
	if (pv->ob_type->tp_as_number == 0) {
	        PyErr_SetString(PyExc_TypeError, "cannot convert to an int, scalar object is not a number.");
		Py_DECREF(pv);
                return NULL;
	}
	if (pv->ob_type->tp_as_number->nb_int == 0) {
	        PyErr_SetString(PyExc_TypeError, "don't know how to convert scalar number to int");
		Py_DECREF(pv);
                return NULL;
	}

	pv2 = pv->ob_type->tp_as_number->nb_int(pv);
	Py_DECREF(pv);
	return pv2;        
}

static PyObject *array_float(PyArrayObject *v) {        
        PyObject *pv, *pv2;
        if (v->nd != 0) {
	        PyErr_SetString(PyExc_TypeError, "only rank-0 arrays can be converted to Python scalars.");
		return NULL;
	}
	pv = v->descr->getitem(v->data);
	if (pv == NULL) return NULL;
	if (pv->ob_type->tp_as_number == 0) {
	        PyErr_SetString(PyExc_TypeError, "cannot convert to an int, scalar object is not a number.");
		Py_DECREF(pv);
                return NULL;
	}
	if (pv->ob_type->tp_as_number->nb_float == 0) {
	        PyErr_SetString(PyExc_TypeError, "don't know how to convert scalar number to float");
		Py_DECREF(pv);
                return NULL;
	}
	pv2 = pv->ob_type->tp_as_number->nb_float(pv);
	Py_DECREF(pv);
	return pv2;        
}

static PyObject *array_long(PyArrayObject *v) {        
        PyObject *pv, *pv2;
        if (v->nd != 0) {
	        PyErr_SetString(PyExc_TypeError, "only rank-0 arrays can be converted to Python scalars.");
		return NULL;
	}
	pv = v->descr->getitem(v->data);
	if (pv->ob_type->tp_as_number == 0) {
	        PyErr_SetString(PyExc_TypeError, "cannot convert to an int, scalar object is not a number.");
                return NULL;
	}
	if (pv->ob_type->tp_as_number->nb_long == 0) {
	        PyErr_SetString(PyExc_TypeError, "don't know how to convert scalar number to long");
                return NULL;
	}
	pv2 = pv->ob_type->tp_as_number->nb_long(pv);
	Py_DECREF(pv);
	return pv2;        
}

static PyObject *array_oct(PyArrayObject *v) {        
        PyObject *pv, *pv2;
        if (v->nd != 0) {
	        PyErr_SetString(PyExc_TypeError, "only rank-0 arrays can be converted to Python scalars.");
		return NULL;
	}
	pv = v->descr->getitem(v->data);
	if (pv->ob_type->tp_as_number == 0) {
	        PyErr_SetString(PyExc_TypeError, "cannot convert to an int, scalar object is not a number.");
                return NULL;
	}
	if (pv->ob_type->tp_as_number->nb_oct == 0) {
	        PyErr_SetString(PyExc_TypeError, "don't know how to convert scalar number to oct");
                return NULL;
	}
	pv2 = pv->ob_type->tp_as_number->nb_oct(pv);
	Py_DECREF(pv);
	return pv2;        
}

static PyObject *array_hex(PyArrayObject *v) {        
        PyObject *pv, *pv2;
        if (v->nd != 0) {
	        PyErr_SetString(PyExc_TypeError, "only rank-0 arrays can be converted to Python scalars.");
		return NULL;
	}
	pv = v->descr->getitem(v->data);
	if (pv->ob_type->tp_as_number == 0) {
	        PyErr_SetString(PyExc_TypeError, "cannot convert to an int, scalar object is not a number.");
                return NULL;
	}
	if (pv->ob_type->tp_as_number->nb_hex == 0) {
	        PyErr_SetString(PyExc_TypeError, "don't know how to convert scalar number to hex");
                return NULL;
	}
        pv2 = pv->ob_type->tp_as_number->nb_hex(pv);
	Py_DECREF(pv);
	return pv2;        
}

/* end of methods added by Travis Oliphant */


static PyNumberMethods array_as_number = {
	(binaryfunc)array_add, /*nb_add*/
	(binaryfunc)array_subtract, /*nb_subtract*/
	(binaryfunc)array_multiply, /*nb_multiply*/
	(binaryfunc)array_divide, /*nb_divide*/
	(binaryfunc)array_remainder, /*nb_remainder*/
	(binaryfunc)array_divmod, /*nb_divmod*/
	(ternaryfunc)array_power, /*nb_power*/
	(unaryfunc)array_negative, 
	(unaryfunc)PyArray_Copy, /*nb_pos*/ 
	(unaryfunc)array_absolute, /* (unaryfunc)array_abs,  */
	(inquiry)array_nonzero, /*nb_nonzero*/
	(unaryfunc)array_invert,		/*nb_invert*/
	(binaryfunc)array_left_shift,		/*nb_lshift*/
	(binaryfunc)array_right_shift,		/*nb_rshift*/
	(binaryfunc)array_bitwise_and,		/*nb_and*/
	(binaryfunc)array_bitwise_xor,		/*nb_xor*/
	(binaryfunc)array_bitwise_or,		/*nb_or*/
	(coercion)array_coerce, /*nb_coerce*/
	(unaryfunc)array_int, /*nb_int*/
	(unaryfunc)array_long, /*nb_long*/
	(unaryfunc)array_float, /*nb_float*/
	(unaryfunc)array_oct,	/*nb_oct*/
	(unaryfunc)array_hex,	/*nb_hex*/
};

static PySequenceMethods array_as_sequence = {
	(inquiry)array_length,		/*sq_length*/
	(binaryfunc)array_add, /*nb_add, concat is numeric add*/
	(intargfunc)NULL, /*nb_multiply, repeat is numeric multiply*/
	(intargfunc)array_item,		/*sq_item*/
	(intintargfunc)array_slice,		/*sq_slice*/
	(intobjargproc)array_ass_item,	/*sq_ass_item*/
	(intintobjargproc)array_ass_slice,	/*sq_ass_slice*/
};

/* -------------------------------------------------------------- */

/*This ugly function quickly prints out an array.  It's likely to disappear */
/*in the near future. */
static int dump_data(char **string, int *n, int *max_n, 
					 char *data, int nd, int *dimensions, 
					 int *strides, PyArray_Descr *descr) {
	PyObject *op, *sp;
	char *ostring;
	int i, N;
	
#define CHECK_MEMORY if (*n >= *max_n-16) { *max_n *= 2; *string = (char *)realloc(*string, *max_n); }
	
	if (nd == 0) {
		
		if ((op = descr->getitem(data)) == NULL) return -1;
		sp = PyObject_Repr(op);
		if (sp == NULL) {Py_DECREF(op); return -1;}
		ostring = PyString_AsString(sp);
		N = PyString_Size(sp)*sizeof(char);
		*n += N;
		CHECK_MEMORY
			memcpy(*string+(*n-N), ostring, N);
		Py_DECREF(sp);
		Py_DECREF(op);
		return 0;
	} else {
		if (nd == 1 && descr->type_num == PyArray_CHAR) {
			N = (dimensions[0]+2)*sizeof(char);
			*n += N;
			CHECK_MEMORY
				(*string)[*n-N] = '"';
			memcpy(*string+(*n-N+1), data, N-2);
			(*string)[*n-1] = '"';
			return 0;
		} else {
			CHECK_MEMORY
				(*string)[*n] = '[';
			*n += 1;
			for(i=0; i<dimensions[0]; i++) {
				if (dump_data(string, n, max_n, data+(*strides)*i, nd-1, dimensions+1, strides+1, descr) < 0) return -1;
				CHECK_MEMORY
					if (i<dimensions[0]-1) {
						(*string)[*n] = ','; 
						(*string)[*n+1] = ' ';
						*n += 2;
					}
			}
			CHECK_MEMORY
				(*string)[*n] = ']'; *n += 1;
			return 0;
		}
	}
	
#undef CHECK_MEMORY
}

static PyObject * array_repr_builtin(PyArrayObject *self) {
	PyObject *ret;
	char *string;
	int n, max_n;
	
	max_n = NBYTES(self)*4*sizeof(char) + 7;
	
	if ((string = (char *)malloc(max_n)) == NULL) {
		PyErr_SetString(PyExc_MemoryError, "out of memory");
		return NULL;
	}
	
	n = 6;
	sprintf(string, "array(");
	
	if (dump_data(&string, &n, &max_n, self->data, self->nd, self->dimensions, 
		self->strides, self->descr) < 0) { free(string); return NULL; }
	
	sprintf(string+n, ", '%c')", self->descr->type);
	
	ret = PyString_FromStringAndSize(string, n+6);
	free(string);
	return ret;
}

static PyObject *PyArray_StrFunction=NULL;
static PyObject *PyArray_ReprFunction=NULL;

extern void PyArray_SetStringFunction(PyObject *op, int repr) {
	if (repr) {
		Py_XDECREF(PyArray_ReprFunction); /* Dispose of previous callback */
		Py_XINCREF(op); /* Add a reference to new callback */
		PyArray_ReprFunction = op; /* Remember new callback */
	} else {
		Py_XDECREF(PyArray_StrFunction); /* Dispose of previous callback */
		Py_XINCREF(op); /* Add a reference to new callback */
		PyArray_StrFunction = op; /* Remember new callback */
	}
}

static PyObject *array_repr(PyArrayObject *self) {
	PyObject *s, *arglist;
	
	if (PyArray_ReprFunction == NULL) {
		s = array_repr_builtin(self);
	} else {
		arglist = Py_BuildValue("(O)", self);
		s = PyEval_CallObject(PyArray_ReprFunction, arglist);
		Py_DECREF(arglist); 
	}
	return s;
}

static PyObject *array_str(PyArrayObject *self) {
	PyObject *s, *arglist;
	
	if (PyArray_StrFunction == NULL) {
		s = array_repr(self);
	} else {
		arglist = Py_BuildValue("(O)", self);
		s = PyEval_CallObject(PyArray_StrFunction, arglist);
		Py_DECREF(arglist); 
	}
	return s;
}

/* No longer implemented here.  Use the default in object.c
static int array_print(PyArrayObject *self, FILE *fp, int flags) {
	PyObject *s;
	int r=-1, n;

	s = array_str(self);

	if (s == NULL || !PyString_Check(s)) {
		Py_XDECREF(s);
		return -1;
	}
	n = PyString_Size(s);
	r = fwrite(PyString_AsString(s), sizeof(char), n, fp);
	Py_DECREF(s);
	return r==n ? 0 : -1;
}
*/

static void byte_swap_vector(void *p, int n, int size) {
	char *a, *b, c;
	
	switch(size) {
	case 2:
		for (a = (char*)p ; n > 0; n--, a += 1) {
			b = a + 1;
			c = *a; *a++ = *b; *b   = c;
		}
		break;
	case 4:
		for (a = (char*)p ; n > 0; n--, a += 2) {
			b = a + 3;
			c = *a; *a++ = *b; *b-- = c;
			c = *a; *a++ = *b; *b   = c;
		}
		break;
	case 8:
		for (a = (char*)p ; n > 0; n--, a += 4) {
			b = a + 7;
			c = *a; *a++ = *b; *b-- = c;
			c = *a; *a++ = *b; *b-- = c;
			c = *a; *a++ = *b; *b-- = c;
			c = *a; *a++ = *b; *b   = c;
		}
		break;
	default:
		break;
	}
}

static char doc_byteswapped[] = "m.byteswapped().  Swaps the bytes in the contents of array m.  Useful for reading data written by a machine with a different byte order.  Returns a new array with the byteswapped data.";

static PyObject * array_byteswap(PyArrayObject *self, PyObject *args) {
	PyArrayObject *ret;
	
	if (!PyArg_ParseTuple(args, "")) return NULL;
	
	if ((ret = (PyArrayObject *)PyArray_Copy(self)) == NULL) return NULL;
	
	if (self->descr->type_num < PyArray_CFLOAT) {
		byte_swap_vector(ret->data, SIZE(self), self->descr->elsize);
	} else {
		byte_swap_vector(ret->data, SIZE(self)*2, self->descr->elsize/2);
	}
	
	return (PyObject *)ret;
}

static char doc_tostring[] = "m.tostring().  Copy the data portion of the array to a python string and return that string.";

static PyObject *array_tostring(PyArrayObject *self, PyObject *args) {
	PyObject *so;
	PyArrayObject *mo;
	
	if (!PyArg_ParseTuple(args, "")) return NULL;
	
	if ((mo = PyArray_CONTIGUOUS(self)) ==NULL) return NULL;
	
	so = PyString_FromStringAndSize(mo->data, NBYTES(mo));
	
	Py_DECREF(mo);
	
	return so;
}


static PyObject *PyArray_ToList(PyObject *self) {
	PyObject *lp;
	PyObject *v;
	int sz, i;
	
	if (!PyArray_Check(self)) return self;
	
	if (((PyArrayObject *)self)->nd == 0) return ((PyArrayObject *)self)->descr->getitem(((PyArrayObject *)self)->data);
	
	sz = ((PyArrayObject *)self)->dimensions[0];
	lp = PyList_New(sz);
	
	for (i=0; i<sz; i++) {
        v=array_item((PyArrayObject *)self, i);
        PyList_SetItem(lp, i, PyArray_ToList(v));
        if (((PyArrayObject *)self)->nd>1){
            Py_DECREF(v);
        }
	}
	
	return lp;
}


static char doc_tolist[] = "m.tolist().  Copy the data portion of the array to a hierarchical python list and return that list.";

static PyObject *array_tolist(PyArrayObject *self, PyObject *args) {
	if (!PyArg_ParseTuple(args, "")) return NULL;
	if (self->nd <= 0) {
		PyErr_SetString(PyExc_ValueError, "Can't convert a 0d array to a list");
		return NULL;
	}
	
	return PyArray_ToList((PyObject *)self);
}

static char doc_toscalar[] = "m.toscalar().  Copy the first data point of the array to a Python scalar and return it.";

static PyObject *array_toscalar(PyArrayObject *self, PyObject *args) {
	if (!PyArg_ParseTuple(args, "")) return NULL;
	return self->descr->getitem(self->data);
}

static char doc_cast[] = "m.astype(t).  Cast array m to type t.  t can be either a string representing a typecode, or a python type object of type int, float, or complex.";

PyObject *array_cast(PyArrayObject *self, PyObject *args) {
	PyObject *op;
	char typecode;
	
	if (!PyArg_ParseTuple(args, "O", &op)) return NULL;
	
	if (PyString_Check(op) && (PyString_Size(op) == 1)) {
		return PyArray_Cast(self, PyString_AS_STRING((PyStringObject *)op)[0]);
	}
	
	if (PyType_Check(op)) {
		typecode = 'O';
		if (((PyTypeObject *)op) == &PyInt_Type) {
			typecode = PyArray_LONG;
		} 
		if (((PyTypeObject *)op) == &PyFloat_Type) {
			typecode = PyArray_DOUBLE;
		} 
		if (((PyTypeObject *)op) == &PyComplex_Type) {
			typecode = PyArray_CDOUBLE;
		} 
		return PyArray_Cast(self, typecode);
	}
	PyErr_SetString(PyExc_ValueError,
		"type must be either a 1-length string, or a python type object");
	return NULL;
}

static char doc_typecode[] = "m.typecode().  Return the single character code indicating the type of the elements of m.";

PyObject *array_typecode(PyArrayObject *self, PyObject *args) {
	if (!PyArg_ParseTuple(args, "")) return NULL;
	
	return PyString_FromStringAndSize(&(self->descr->type), 1);
}

static char doc_itemsize[] = "m.itemsize().  Return the size in bytes of a single element of the array m.";

PyObject *array_itemsize(PyArrayObject *self, PyObject *args) {
	if (!PyArg_ParseTuple(args, "")) return NULL;
	
	return PyInt_FromLong(self->descr->elsize);
}

static char doc_contiguous[] = "m.contiguous().  Return true if the memory used by the array m is contiguous.  A non-contiguous array can be converted to a contiguous one by m.copy().";

PyObject *array_contiguous(PyArrayObject *self, PyObject *args) {
	if (!PyArg_ParseTuple(args, "")) return NULL;
	
	return PyInt_FromLong(ISCONTIGUOUS(self));
}

/* Added savespace methods Jan. 2000 -- Travis Oliphant */
static char doc_savespace[] = "m.savespace(flag=1).  Change the savespace parameter of the array m to flag.";

PyObject *array_savespace(PyArrayObject *self, PyObject *args, PyObject *kws) 
{
        char flag=1;
	char *kwd[] = {"flag",NULL};

        if (!PyArg_ParseTupleAndKeywords(args,kws,"|b", kwd, &flag))
                return NULL;

	if (flag == 0)
	        self->flags &= ~SAVESPACE;
        else 
  	        self->flags |= SAVESPACE;
    
        Py_INCREF(Py_None);
	return Py_None;
}

static char doc_spacesaver[] = "m.spacesaver().  Return true if the spacesaving parameter is set for m.";

PyObject *array_spacesaver(PyArrayObject *self, PyObject *args) {
	if (!PyArg_ParseTuple(args, "")) return NULL;
	
	return PyInt_FromLong(PyArray_ISSPACESAVER(self) != 0);
}

	  
static char doc_copy[] = "m.__copy__().";

PyObject *array_copy(PyArrayObject *self, PyObject *args) {
	if (!PyArg_ParseTuple(args, "")) return NULL;
	
	return PyArray_Copy(self);
}


/* Changed by Travis Oliphant Jan, 2000 */
int array_compare(PyArrayObject *self, PyObject *other) {
        PyArrayObject *aother;
        PyObject *o1, *o2;
        int val, result;
        aother = (PyArrayObject *)other;
        if (self->nd != 0 || aother->nd != 0) {
	        PyErr_SetString(PyExc_TypeError,
		 "Comparison of multiarray objects other than rank-0 arrays is not implemented.");
	        return -1;
	}
        o1 = self->descr->getitem(self->data);
        o2 = aother->descr->getitem(aother->data);
        if (o1 == NULL || o2 == NULL) return -1;
        val = PyObject_Cmp(o1,o2,&result);
        Py_DECREF(o1);
        Py_DECREF(o2);
        if (val < 0) {
                PyErr_SetString(PyExc_TypeError, "objects cannot be compared.");
                return -1;
        }
        return result;
}

static struct PyMethodDef array_methods[] = {
  {"tostring",   (PyCFunction)array_tostring,	1, doc_tostring},
  {"tolist",   (PyCFunction)array_tolist,	1, doc_tolist},
  {"toscalar", (PyCFunction)array_toscalar,     METH_VARARGS, doc_toscalar},
  {"byteswapped",   (PyCFunction)array_byteswap,	1, doc_byteswapped},
  {"astype", (PyCFunction)array_cast, 1, doc_cast},
  {"typecode",   (PyCFunction)array_typecode, 1, doc_typecode},
  {"itemsize",   (PyCFunction)array_itemsize, 1, doc_itemsize},
  {"iscontiguous", (PyCFunction)array_contiguous, 1, doc_contiguous},
  {"savespace", (PyCFunction)array_savespace, METH_VARARGS | METH_KEYWORDS, doc_savespace},
  {"spacesaver", (PyCFunction)array_spacesaver, METH_VARARGS, doc_spacesaver},
  {"__copy__", (PyCFunction)array_copy, 1, doc_copy},
  {NULL,		NULL}		/* sentinel */
};

/* ---------- */
static PyObject *array_getattr(PyArrayObject *self, char *name) {
	PyArrayObject *ret;
	
	if (strcmp(name, "shape") == 0) {
		PyObject *s, *o;
		int i;
		
		if ((s=PyTuple_New(self->nd)) == NULL) return NULL;
		
		for(i=self->nd; --i >= 0;) {
			if ((o=PyInt_FromLong(self->dimensions[i])) == NULL) return NULL;
			if (PyTuple_SetItem(s,i,o) == -1) return NULL;
		}
		return s;
	}

	if (strcmp(name, "real") == 0) {
		if (self->descr->type_num == PyArray_CFLOAT || 
			self->descr->type_num == PyArray_CDOUBLE) {
			ret = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
				self->dimensions, 
				self->descr->type_num-2,
				self->data);
			if (ret == NULL) return NULL;
			memcpy(ret->strides, self->strides, sizeof(int)*ret->nd);
			ret->flags &= ~CONTIGUOUS;
			Py_INCREF(self);
			ret->base = (PyObject *)self;
			return (PyObject *)ret;
		} else {
			ret = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
				self->dimensions,
				self->descr->type_num,
				self->data);
			if (ret == NULL) return NULL;
			Py_INCREF(self);
			ret->base = (PyObject *)self;
		}
	}
	
	if ((strcmp(name, "imaginary") == 0) || (strcmp(name, "imag") == 0)) { 
		if (self->descr->type_num == PyArray_CFLOAT || 
			self->descr->type_num == PyArray_CDOUBLE) {
			ret = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
				self->dimensions, 
				self->descr->type_num-2,
				self->data+self->descr->elsize/2);
			if (ret == NULL) return NULL;
			memcpy(ret->strides, self->strides, sizeof(int)*ret->nd);
			ret->flags &= ~CONTIGUOUS;
			Py_INCREF(self);
			ret->base = (PyObject *)self;
			return (PyObject *)ret;
		} else {
			PyErr_SetString(PyExc_ValueError, "No imaginary part to real array");
			return NULL;
		}
	}
	
	if (strcmp(name, "flat") == 0) {
		int n;
		n = SIZE(self);
		if (!ISCONTIGUOUS(self)) {
			PyErr_SetString(PyExc_ValueError, "flattened indexing only available for contiguous array");
			return NULL;
		}
		ret = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(1, &n,
			self->descr,
			self->data);
		if (ret == NULL) return NULL;  
		Py_INCREF(self);
		ret->base = (PyObject *)self;
		return (PyObject *)ret;
	}


	return Py_FindMethod(array_methods, (PyObject *)self, name);
}

static int array_setattr(PyArrayObject *self, char *name, PyObject *op) {
	PyArrayObject *ap;
	int ret;
	
	if (strcmp(name, "shape") == 0) {
		/* This can be made more efficient by copying code from array_reshape if needed */
		if ((ap = (PyArrayObject *)PyArray_Reshape(self, op)) == NULL) return -1;
		if(self->flags & OWN_DIMENSIONS) free(self->dimensions);
		self->dimensions = ap->dimensions;
		if(self->flags & OWN_STRIDES) free(self->strides);
		self->strides = ap->strides;
		self->nd = ap->nd;
		self->flags &= ~(OWN_DIMENSIONS | OWN_STRIDES);
		self->flags |= ap->flags & (OWN_DIMENSIONS | OWN_STRIDES);
		ap->flags &= ~(OWN_DIMENSIONS | OWN_STRIDES);
		Py_DECREF(ap);
		return 0;
	}
	
	if (strcmp(name, "real") == 0) {
		if (self->descr->type_num == PyArray_CFLOAT || 
			self->descr->type_num == PyArray_CDOUBLE) {
			ap = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
				self->dimensions, 
				self->descr->type_num-2,
				self->data);
			if (ap == NULL) return -1;
			memcpy(ap->strides, self->strides, sizeof(int)*ap->nd);
			ap->flags &= ~CONTIGUOUS;
			ret = PyArray_CopyObject(ap, op);
			Py_DECREF(ap);
			return ret;
		} else {
			return PyArray_CopyObject(self, op);
		}
	}
	
	if ((strcmp(name, "imaginary") == 0) || (strcmp(name, "imag") == 0)) { 
		if (self->descr->type_num == PyArray_CFLOAT || 
			self->descr->type_num == PyArray_CDOUBLE) {
			ap = (PyArrayObject *)PyArray_FromDimsAndData(self->nd, 
				self->dimensions, 
				self->descr->type_num-2,
				self->data+self->descr->elsize/2);
			if (ap == NULL) return -1;
			memcpy(ap->strides, self->strides, sizeof(int)*ap->nd);
			ap->flags &= ~CONTIGUOUS;
			ret = PyArray_CopyObject(ap, op);
			Py_DECREF(ap);
			return ret;
		} else {
			PyErr_SetString(PyExc_ValueError, "No imaginary part to real array");
			return -1;
		}
	}


	PyErr_SetString(PyExc_AttributeError, "Attribute does not exist or cannot be set");
	return -1;
}

static char Arraytype__doc__[] = 
"A array object represents a multidimensional, homogeneous array of basic values.  It has the folowing data members, m.shape (the size of each dimension in the array), m.itemsize (the size (in bytes) of each element of the array), and m.typecode (a character representing the type of the matrices elements).  Matrices are sequence, mapping and numeric objects.  Sequence indexing is similar to lists, with single indices returning a reference that points to the old matrices data, and slices returning by copy.  A array is also allowed to be indexed by a sequence of sequences or ints.  Each member of the sequence indexes the corresponding dimension of the array.  If it is a sequence, it returns the elments along that dimension listed in the sequence, if a single number, it returns just that element (reducing the dimensionality of the returned array by 1.  Numeric operations operate on matrices in an element-wise fashion.";


PyTypeObject PyArray_Type = { 
	PyObject_HEAD_INIT(0)
	0,				/*ob_size*/
	"array",			/*tp_name*/
	sizeof(PyArrayObject),		/*tp_basicsize*/
	0,				/*tp_itemsize*/
		/* methods */
	(destructor)array_dealloc,	/*tp_dealloc*/
	(printfunc)NULL,		/*tp_print*/
	(getattrfunc)array_getattr,	/*tp_getattr*/
	(setattrfunc)array_setattr,	/*tp_setattr*/
	(cmpfunc)array_compare,		/*tp_compare*/
	(reprfunc)array_repr,		/*tp_repr*/
	&array_as_number,			/*tp_as_number*/
	&array_as_sequence,		/*tp_as_sequence*/
	&array_as_mapping, 		/*tp_as_mapping*/
	(hashfunc)0,		/*tp_hash*/
	(ternaryfunc)0,		/*tp_call*/
	(reprfunc)array_str,		/*tp_str*/
		
		/* Space for future expansion */
	0L,0L,0L,0L,
	Arraytype__doc__ /* Documentation string */
};


/* The rest of this code is to build the right kind of array from a python */
/* object. */

static int discover_depth(PyObject *s, int max, int stop_at_string) {
	int d=0;
	PyObject *e;
	
	if(max < 1) return -1;
	
	/*if(! PySequence_Check(s) || PyInstance_Check(s)) return 0;
	/*
	Relax this to allow class instances.  Dangerous? */
	if(! PySequence_Check(s) || PySequence_Length(s) < 0) {
			PyErr_Clear(); return 0;
	}
	if(PyString_Check(s)) return stop_at_string ? 0:1;
	if (PySequence_Length(s) == 0) return 1;
	
	if ((e=PySequence_GetItem(s,0)) == NULL) return -1;
	if(e!=s)
	{
		d=discover_depth(e,max-1, stop_at_string);
		if(d >= 0) d++;
	}
	Py_DECREF(e);
	return d;
}

static int discover_dimensions(PyObject *s, int nd, int *d) {
	PyObject *e;
	int r, n, i, n_lower;
	
	n=PyObject_Length(s);
	*d = n;
	if(*d < 0) return -1;
	if(nd <= 1) return 0;
	n_lower = 0;
	for(i=0; i<n; i++) {
		if ((e=PySequence_GetItem(s,i)) == NULL) return -1;
		r=discover_dimensions(e,nd-1,d+1);
		Py_DECREF(e);
		
		if (r == -1) return -1; 
		if (d[1] > n_lower) n_lower = d[1];
	}
	d[1] = n_lower;
	
	return 0;
}

#ifndef max
#define max(a,b) (a)>(b)?(a):(b);
#endif

static int array_objecttype(PyObject *op, int minimum_type, int savespaceflag) {
	int l;
	PyObject *ip;
	
	if (minimum_type == -1) return -1;
	
	if (PyArray_Check(op)) return max((int)((PyArrayObject *)op)->descr->type_num, minimum_type);
	
	if (PyInstance_Check(op)) {
		if (PyObject_HasAttrString(op, "__array__")) {
			PyObject *ap, *arglist;
			
			arglist = Py_BuildValue("()");
			ap = PyObject_GetAttrString(op, "__array__");
			ip = PyEval_CallObject(ap, arglist);
			Py_DECREF(arglist);
			
			return max(minimum_type, (int)((PyArrayObject *)ip)->descr->type_num);
		} else {
			if (PySequence_Length(op)<0)
				PyErr_Clear();
				return (int)PyArray_OBJECT;
		}
	}
	
	if (PyString_Check(op)) {
		return max(minimum_type, (int)PyArray_CHAR);
	}
	
	if (PySequence_Check(op)) {
		l = PyObject_Length(op);
		if (l == 0 && minimum_type == 0) minimum_type = (savespaceflag ? PyArray_SHORT : PyArray_LONG);
		while (--l >= 0) {
			ip = PySequence_GetItem(op, l);
			minimum_type = array_objecttype(ip, minimum_type, savespaceflag);
			Py_DECREF(ip);
		}
		return minimum_type;
	}
	
	if (PyInt_Check(op)) {
		return max(minimum_type, (int) (savespaceflag ? PyArray_SHORT: PyArray_LONG));
	} else {
		if (PyFloat_Check(op)) {
			return max(minimum_type, (int) (savespaceflag ? PyArray_FLOAT : PyArray_DOUBLE));
		} else {
			if (PyComplex_Check(op)) {
			        return max(minimum_type, (int) (savespaceflag ? PyArray_CFLOAT : PyArray_CDOUBLE));
			} else {
				return (int)PyArray_OBJECT;
			}
		}
	}
}

static int Assign_Array(PyArrayObject *self, PyObject *v) {
	PyObject *e;
	int l, r;
	
	if (!PySequence_Check(v)) {
		PyErr_SetString(PyExc_ValueError,"assignment from non-sequence");
		return -1;
	}
	
	l=PyObject_Length(v);
	if(l < 0) return -1;
	
	while(--l >= 0)
	{
		e=PySequence_GetItem(v,l);
		if (e == NULL) return -1;
		r = PySequence_SetItem((PyObject*)self,l,e);
		Py_DECREF(e);
		if(r == -1) return -1;
	}
	return 0;
}

static PyObject *
Array_FromSequence(PyObject *s, int type, int min_depth, int max_depth)
{
	PyArrayObject *r;
	int nd, *d;
	
	if (!PySequence_Check(s)) {
		PyErr_SetString(PyExc_ValueError,"expect source sequence");
		return NULL;
	}
	if (!((nd=discover_depth(s,99, type == PyArray_OBJECT || type == 'O')) > 0)) {
		PyErr_SetString(PyExc_ValueError, "invalid input sequence");
		return NULL;
	}
	if ((max_depth && nd > max_depth) || (min_depth && nd < min_depth)) {
		PyErr_SetString(PyExc_ValueError, "invalid number of dimensions");
		return NULL;
	}
	if ((d=(int *)malloc(nd*sizeof(int))) == NULL) {
		PyErr_SetString(PyExc_MemoryError,"out of memory");
	}
	if(discover_dimensions(s,nd,d) == -1)
	{
		free(d);
		return 0;
	}
	
	if (type == PyArray_CHAR && nd > 0 && d[nd-1] == 1) {
		nd = nd-1;
	}
	
	r=(PyArrayObject*)PyArray_FromDims(nd,d,type);
	free(d);
	if(! r) return NULL;
	if(Assign_Array(r,s) == -1)
	{
		Py_DECREF(r);
		return NULL;
	}
	return (PyObject*)r;
}

PyObject *PyArray_FromScalar(PyObject *op, int type) {
	PyArrayObject *ret;
	
	if ((ret = (PyArrayObject *)PyArray_FromDims(0, NULL, type)) == NULL)
		return NULL;
	
	ret->descr->setitem(op, ret->data);
	
	if (PyErr_Occurred()) {
		array_dealloc(ret);
		return NULL;
	} else {
		return (PyObject *)ret;
	}
}

PyObject *PyArray_Cast(PyArrayObject *mp, int type) {
	PyArrayObject *rp, *tmp;
	
	if (mp->descr->type_num == PyArray_OBJECT) {
		return PyArray_FromObject((PyObject *)mp, type, mp->nd, mp->nd);
	}
	
	if ((tmp = PyArray_CONTIGUOUS(mp)) == NULL) return NULL;
	
	rp = (PyArrayObject *)PyArray_FromDims(tmp->nd, tmp->dimensions, type);
	mp->descr->cast[(int)rp->descr->type_num](tmp->data, 1, rp->data, 1, SIZE(mp));
	
	Py_DECREF(tmp);
	return (PyObject *)rp;
}

#define ByCopy 1
#define Contiguous 2

PyObject *array_fromobject(PyObject *op_in, int type, int min_depth, int max_depth, int flags) {
        /* This is the main code to make a NumPy array from a Python Object.  It is 
           called from lot's of different places which is why there are so many checks.  
           The comments try to explain some of the checks. */

	PyObject *r, *op;
        char temp;
	int real_type;

	/* The type argument can have a (char-based) flag set which states that this array 
	   should be spacesaving */
	temp = type;
	temp &= ~SAVESPACEBIT;  
	real_type = temp;       /* Type without the spacesaving flag  */
	
	r = NULL;

	/* Code that returns the object to convert for a non multiarray input object 
           from the __array__ attribute of the object. */
	if (!PyArray_Check(op_in) && PyObject_HasAttrString(op_in, "__array__")) {
		PyObject *ap, *arglist;
		
		if (real_type == PyArray_NOTYPE) {
			arglist = Py_BuildValue("()");
		} else {
			arglist = Py_BuildValue("(c)", real_type);
		}
		ap = PyObject_GetAttrString(op_in, "__array__");
		op = PyEval_CallObject(ap, arglist);
		Py_DECREF(ap);
		Py_DECREF(arglist);
		
 		if (op == NULL) return NULL;
 	} else {
		op = op_in;
		Py_INCREF(op);
	}

	/* If the typecode is NOTYPE get the real_type from the object and 
           reset type with the SAVESPACEBIT */
	if (real_type == PyArray_NOTYPE) {
		real_type = array_objecttype(op, 0, type & SAVESPACEBIT);
		if (type & SAVESPACEBIT) {
		        type = real_type | SAVESPACEBIT;
		}
	}
	
	/* Is input object an array but not an object array and we don't want an object
           array output */
	if (PyArray_Check(op) && (((PyArrayObject *)op)->descr->type_num != PyArray_OBJECT || 
		real_type == PyArray_OBJECT || real_type == 'O')) {
	        /* Is the desired output type the same as the input array type */
		if ((((PyArrayObject *)op)->descr->type_num == real_type ||
			((PyArrayObject *)op)->descr->type == real_type)) {
		       /* Should the output be copied (1) because it's a flag, (2) because a contiguous
                          array is requested */
		        if ((flags & ByCopy) || (flags&Contiguous && !ISCONTIGUOUS((PyArrayObject *)op))) {
				r = PyArray_Copy((PyArrayObject *)op);
			} 
                        /* If no copy then just increase the reference count and return the input */
                        else {  
				Py_INCREF(op);
				r = op;
			}
		} 
                /* The desired output type is different than the input array type */
                else {
		        /* If the typecode is too big, (i.e. it is a character) make it a number. */
			if (real_type > PyArray_NTYPES) real_type = PyArray_DescrFromType(real_type)->type_num;
                        /* Cast to the desired type if we can do it safely.  If the SAVESPACE 
                           bit is set do it anyway.  Also cast if source is a rank-0 array to mimic
                           behavior with Python scalars */
			if (PyArray_CanCastSafely(((PyArrayObject *)op)->descr->type_num, real_type) || (type & SAVESPACEBIT) || (((PyArrayObject *)op)->nd == 0))
				r = PyArray_Cast((PyArrayObject *)op, real_type);
			else {
				PyErr_SetString(PyExc_TypeError, "Array can not be safely cast to required type");
				r = NULL;
			}
		}
	}
        /* The input object is either not an array or it is an object array, 
           or the desired output type is object. */ 
        else {
		r = Array_FromSequence(op, real_type, min_depth,max_depth);
		if (r == NULL && min_depth <= 0) {
				PyErr_Clear();
				r = PyArray_FromScalar(op, real_type);
		}
	}
	Py_DECREF(op);
	/* If we didn't succed return NULL */
	if (r == NULL) return NULL;
	
	if(!PyArray_Check(r)) {
		PyErr_SetString(PyExc_ValueError, "Internal error array_fromobject not producing an array");
		return NULL;
	}
	if (min_depth != 0 && ((PyArrayObject *)r)->nd < min_depth) {
		Py_DECREF(r);
		PyErr_SetString(PyExc_ValueError, 
			"Object of too small depth for desired array");
		return NULL;
	}
	if (max_depth != 0 && ((PyArrayObject *)r)->nd > max_depth) {
		Py_DECREF(r);
		PyErr_SetString(PyExc_ValueError, 
			"Object too deep for desired array");
		return NULL;
	}
	return r;
}


int PyArray_ObjectType(PyObject *op, int minimum_type) {
        return array_objecttype(op, minimum_type, 0);
}

PyObject *PyArray_FromObject(PyObject *op, int type, int min_depth, int max_depth) {
	return array_fromobject(op, type, min_depth, max_depth, 0);
}
PyObject *PyArray_ContiguousFromObject(PyObject *op, int type, int min_depth, int max_depth) {
	return array_fromobject(op, type, min_depth, max_depth, Contiguous);
}
PyObject *PyArray_CopyFromObject(PyObject *op, int type, int min_depth, int max_depth) {
	return array_fromobject(op, type, min_depth, max_depth, ByCopy);
}

extern int PyArray_CanCastSafely(int fromtype, int totype) {
	if (fromtype == totype) return 1;
	if (totype == PyArray_OBJECT) return 1;
	if (fromtype == PyArray_OBJECT) return 0;
	
	switch(fromtype) {
	case PyArray_CHAR:
		return 0;
	case PyArray_UBYTE:
		return (totype >= PyArray_SHORT);
	case PyArray_SBYTE:
	case PyArray_SHORT:
		return (totype > fromtype);
	case PyArray_INT:
	case PyArray_LONG:
		/*Shouldn't allow Longs to be cast to Ints, not safe on 64-bit machines!*/
		return (totype >= PyArray_INT) && (totype != PyArray_FLOAT); 
	case PyArray_FLOAT:
		return (totype > PyArray_FLOAT);
	case PyArray_DOUBLE:
		return (totype == PyArray_CDOUBLE);
	case PyArray_CFLOAT:
		return (totype == PyArray_CDOUBLE);
	case PyArray_CDOUBLE:
		return 0;
	default:
		return 0;
	}
}

extern int PyArray_As1D(PyObject **op, char **ptr, int *d1, int typecode) {
	PyArrayObject *ap;
	
	if ((ap = (PyArrayObject *)PyArray_ContiguousFromObject(*op, typecode, 1,
		1)) == NULL)
		return -1;
	
	*op = (PyObject *)ap;
	*ptr = ap->data;
	*d1 = ap->dimensions[0];
	return 0;
}

extern int PyArray_As2D(PyObject **op, char ***ptr, int *d1, int *d2, int typecode) {
	PyArrayObject *ap;
	int i, n, size;
	char **data;
	
	if ((ap = (PyArrayObject *)PyArray_ContiguousFromObject(*op, typecode, 2, 2)) == NULL)
		return -1;
	
	n = ap->dimensions[0];
	data = (char **)malloc(n*sizeof(char *));
	for(i=0; i<n; i++) {
		data[i] = ap->data + i*ap->strides[0];
	}
	*op = (PyObject *)ap;
	*ptr = data;
	*d1 = ap->dimensions[0];
	*d2 = ap->dimensions[1];
	return 0;
}


extern int PyArray_Free(PyObject *op, char *ptr) {
	PyArrayObject *ap = (PyArrayObject *)op;
	int i, n;
	
	if (ap->nd > 2) return -1;
	if (ap->nd == 3) {
		n = ap->dimensions[0];
		for (i=0; i<n; i++) {
			free(((char **)ptr)[i]);
		}
	}
	if (ap->nd >= 2) {
		free(ptr);
	}
	Py_DECREF(ap);
	return 0;
}

extern PyObject * PyArray_Reshape(PyArrayObject *self, PyObject *shape) {
	int i, n, s_original, i_unknown, s_known;
	int *dimensions;
	PyArrayObject *ret;
	
	if (!PyArray_ISCONTIGUOUS(self)) {
		PyErr_SetString(PyExc_ValueError, "reshape only works on contiguous matrices");
		return NULL;
	}
	
	if (PyArray_As1D(&shape, (char **)&dimensions, &n, PyArray_INT) == -1) return NULL;
	
	s_known = 1;
	i_unknown = -1;
	for(i=0; i<n; i++) {
		if (dimensions[i] < 0) {
			if (i_unknown == -1) {
				i_unknown = i;
			} else {
				PyErr_SetString(PyExc_ValueError, "can only specify one unknown dimension");
				goto fail;
			}
		} else {
			s_known *= dimensions[i];
		}
	}
	
	s_original = PyArray_SIZE(self);
	
	if (i_unknown >= 0) {
		if ((s_known == 0) || (s_original % s_known != 0)) {
			PyErr_SetString(PyExc_ValueError, "total size of new array must be unchanged");
			goto fail;
		}
		dimensions[i_unknown] = s_original/s_known;
	} else {
		if (s_original != s_known) {
			PyErr_SetString(PyExc_ValueError, "total size of new array must be unchanged");
			goto fail;
		}
	}
	if ((ret = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(n, dimensions, 
		self->descr,
		self->data)) == NULL)
		goto fail;
	
	Py_INCREF(self);
	ret->base = (PyObject *)self;
	
	PyArray_Free(shape, (char *)dimensions);
	return (PyObject *)ret;
	
fail:
	PyArray_Free(shape, (char *)dimensions);
	return NULL;
}

extern PyObject *PyArray_Take(PyObject *self0, PyObject *indices0, int axis) {
	PyArrayObject *self, *indices, *ret;
	int nd, shape[MAX_DIMS];
	int i, j, chunk, n, m, max_item, tmp;
	char *src, *dest;
	
	indices = ret = NULL;
	self = (PyArrayObject *)PyArray_ContiguousFromObject(self0, PyArray_NOTYPE, 1, 0);
	if (self == NULL) return NULL;
	
	if (axis < 0) axis = axis + self->nd;
	if ((axis < 0) || (axis >= self->nd)) {
		PyErr_SetString(PyExc_ValueError, "Invalid axis for this array");
		goto fail;
	}
	
	indices = (PyArrayObject *)PyArray_ContiguousFromObject(indices0, PyArray_LONG, 1, 0);
	if (indices == NULL) goto fail;
	
	n = m = chunk = 1;
	nd = self->nd + indices->nd - 1;
	for (i=0; i< nd; i++) {
		if (i < axis) {
			shape[i] = self->dimensions[i];
			n *= shape[i];
		} else {
			if (i < axis+indices->nd) {
				shape[i] = indices->dimensions[i-axis];
				m *= shape[i];
			} else {
				shape[i] = self->dimensions[i-indices->nd+1];
				chunk *= shape[i];
			}
		}
	}
	ret = (PyArrayObject *)PyArray_FromDims(nd, shape, self->descr->type_num);
	
	if (ret == NULL) goto fail;
	
	max_item = self->dimensions[axis];
	chunk = chunk * ret->descr->elsize;
	src = self->data;
	dest = ret->data;
	
	for(i=0; i<n; i++) {
		for(j=0; j<m; j++) {
			tmp = ((long *)(indices->data))[j];
			if (tmp < 0) tmp = tmp+max_item;
			if ((tmp < 0) || (tmp >= max_item)) {
				PyErr_SetString(PyExc_IndexError, "Index out of range for array");
				goto fail;
			}
			memcpy(dest, src+tmp*chunk, chunk);
			dest += chunk;
		}
		src += chunk*max_item;
	}
	
	PyArray_INCREF(ret);

	Py_XDECREF(indices);
	Py_XDECREF(self);

	return (PyObject *)ret;
	
	
fail:
	Py_XDECREF(ret);
	Py_XDECREF(indices);
	Py_XDECREF(self);
	return NULL;
}

extern PyObject *PyArray_Put(PyObject *self0, PyObject *indices0, 
                             PyObject* values0) {
	PyArrayObject  *indices, *values, *self;
	int i, j, chunk, ni, max_item, nv, tmp; 
	char *src, *dest;

        indices = NULL;
        values = NULL;

	if (!PyArray_Check(self0)) {
          PyErr_SetString(PyExc_ValueError, "put: first argument must be an array");
	  return NULL;
        }
        self = (PyArrayObject*) self0;
	if (!PyArray_ISCONTIGUOUS(self)) {
          PyErr_SetString(PyExc_ValueError, "put: first argument must be contiguous");
	  return NULL;
        }
        max_item = PyArray_SIZE(self);
	dest = self->data;
	chunk = self->descr->elsize;

	indices = (PyArrayObject *)PyArray_ContiguousFromObject(indices0, PyArray_LONG, 0, 0);
	if (indices == NULL) goto fail;
        ni = PyArray_SIZE(indices);

	values = (PyArrayObject *)PyArray_ContiguousFromObject(values0, 
                  self->descr->type, 0, 0);
	if (values == NULL) goto fail;
        nv = PyArray_SIZE(values);

	for(i=0; i<ni; i++) {
	   src = values->data + chunk * (i % nv);
 	   tmp = ((long *)(indices->data))[i];
	   if (tmp < 0) tmp = tmp+max_item;
  	   if ((tmp < 0) || (tmp >= max_item)) {
		PyErr_SetString(PyExc_IndexError, "Index out of range for array");
		goto fail;
	   }
	   memcpy(dest + tmp * chunk, src, chunk);
	}

	Py_XDECREF(values);
	Py_XDECREF(indices);
        Py_INCREF(Py_None);
        return Py_None;
	
fail:
	Py_XDECREF(indices);
	Py_XDECREF(values);
	return NULL;
}

extern PyObject *PyArray_PutMask(PyObject *self0, PyObject *mask0, 
                             PyObject* values0) {
	PyArrayObject  *mask, *values, *self;
	int i, j, chunk, ni, max_item, nv, tmp; 
	char *src, *dest;

        mask = NULL;
        values = NULL;

	if (!PyArray_Check(self0)) {
          PyErr_SetString(PyExc_ValueError, "putmask: first argument must be an array");
	  return NULL;
        }
        self = (PyArrayObject*) self0;
	if (!PyArray_ISCONTIGUOUS(self)) {
          PyErr_SetString(PyExc_ValueError, "putmask: first argument must be contiguous");
	  return NULL;
        }
        max_item = PyArray_SIZE(self);
	dest = self->data;
	chunk = self->descr->elsize;

	mask = (PyArrayObject *)PyArray_ContiguousFromObject(mask0, PyArray_LONG, 0, 0);
	if (mask == NULL) goto fail;
        ni = PyArray_SIZE(mask);
        if (ni != max_item) {
          PyErr_SetString(PyExc_ValueError, "putmask: mask and data must be the same size.");
          goto fail;
        }

	values = (PyArrayObject *)PyArray_ContiguousFromObject(values0, 
                  self->descr->type, 0, 0);
	if (values == NULL) goto fail;
        nv = PyArray_SIZE(values);

	for(i=0; i<ni; i++) {
	   src = values->data + chunk * (i % nv);
 	   tmp = ((long *)(mask->data))[i];
           if (tmp) {
	       memcpy(dest + i * chunk, src, chunk);
           }
	}

	Py_XDECREF(values);
	Py_XDECREF(mask);
        Py_INCREF(Py_None);
        return Py_None;
	
fail:
	Py_XDECREF(mask);
	Py_XDECREF(values);
	return NULL;
}


/* This conversion function can be used with the "O&" argument for
   PyArg_ParseTuple, as the O! form will not work with subclasses of
   PyArray_Type  -- courtesy of Travis Oliphant. */
    
extern int PyArray_Converter(PyObject *object, PyObject **address) {
  if (PyArray_Check(object)) {
    *address = object;
    return 1;
  }
  else {
    PyErr_SetString(PyExc_TypeError,"expected Array object in one of the arguments");
    return 0;
  }
}



syntax highlighted by Code2HTML, v. 0.9.1