/*
   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.
*/

#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"

/* 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) {
		for(j=0; j<copies; j++) {
			memcpy(dest, src, elsize);
			dest += 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;
	
	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);
	for(i=0; i<n; i++, data++) Py_XINCREF(*data);
	
	if (!ISCONTIGUOUS(mp)) free(data);
	
	return 0;
}

extern int PyArray_XDECREF(PyArrayObject *mp) {
	int i, n;
	PyObject **data;
	
	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);
	for(i=0; i<n; i++, data++) Py_XDECREF(*data);
	
	if (!ISCONTIGUOUS(mp)) free(data);
	
	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); }
#ifndef NUMPY_NOATTRIBUTES
	if (ret != -1) {
		dest->attributes = src->attributes;
		Py_INCREF(dest->attributes);
	}
#endif
	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_FromDimsAndData(int nd, int *d, int type, char *data) {
	PyArrayObject *self;
	int i,sd;
	PyArray_Descr *descr;
	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 ((descr = PyArray_DescrFromType(type)) == NULL) 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;
#ifndef NUMPY_NOATTRIBUTES
	self->attributes = Py_None; Py_INCREF(Py_None);
#endif
	
	return (PyObject*)self;
	
fail:
	
	if (flags & OWN_DATA) free(data);
	if (dimensions != NULL) free(dimensions);
	if (strides != NULL) free(strides);
	return NULL;    
}

PyObject *PyArray_FromDims(int nd, int *d, int type) {
	return PyArray_FromDimsAndData(nd, d, type, NULL);
}


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); 
#ifndef NUMPY_NOATTRIBUTES
	Py_XDECREF(self->attributes);
#endif
	
	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;
#ifndef NUMPY_NOATTRIBUTES
	        r->attributes = Py_None; Py_INCREF(Py_None);
#endif
		Py_INCREF(self);
		return (PyObject*)r;
	} else {
		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) {
		op = array_item(mp, 0);
		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_FromDimsAndData(self->nd, self->dimensions, self->descr->type_num, data);
	self->dimensions[0] = l;
	
	if (!ISCONTIGUOUS(self)) r->flags &= ~CONTIGUOUS;
	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;
}

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_FromDimsAndData(nd, 
		dimensions,
		self->descr->type_num, 
		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;
	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; 
	
	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*/
};

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

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);
}


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*/
	0, /*nb_int*/
	0, /*nb_long*/
	0, /*nb_float*/
	0,		/*nb_oct*/
	0,		/*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 array_item((PyArrayObject *)self, 0);
	
	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_cast[] = "m.cast(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));
}

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

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

int array_compare(PyArrayObject *self, PyObject *other) {
	PyErr_SetString(PyExc_TypeError,
		"Comparison of multiarray objects is not implemented.");
	return -1;
}


static struct PyMethodDef array_methods[] = {
	{"tostring",   (PyCFunction)array_tostring,	1, doc_tostring},
	{"tolist",   (PyCFunction)array_tolist,	1, doc_tolist},
	{"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},
	{"__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_FromDimsAndData(1, &n,
			self->descr->type_num,
			self->data);
		if (ret == NULL) return NULL;  
		Py_INCREF(self);
		ret->base = (PyObject *)self;
		return (PyObject *)ret;
	}
#ifndef NUMPY_NOEXTRA
	if (strcmp(name, "attributes") == 0) {
	    Py_INCREF(self->attributes);
	    return self->attributes;
	}
#endif
	
	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;
		}
	}
#ifndef NUMPY_NOEXTRA	
	if (strcmp(name, "attributes") == 0) {
	    Py_INCREF(op);
	    Py_DECREF(self->attributes);
	    self->attributes = op;
	    return 0;
	}
#endif

	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

int PyArray_ObjectType(PyObject *op, int minimum_type) {
	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 = PyArray_LONG;
		while (--l >= 0) {
			ip = PySequence_GetItem(op, l);
			minimum_type = PyArray_ObjectType(ip, minimum_type);
			Py_DECREF(ip);
		}
		return minimum_type;
	}
	
	if (PyInt_Check(op)) {
		return max(minimum_type, (int)PyArray_LONG);
	} else {
		if (PyFloat_Check(op)) {
			return max(minimum_type, (int)PyArray_DOUBLE);
		} else {
			if (PyComplex_Check(op)) {
				return max(minimum_type, (int)PyArray_CDOUBLE);
			} else {
				return (int)PyArray_OBJECT;
			}
		}
	}
}

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;
}

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)) > 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) {
	PyObject *r, *op;
	
	r = NULL;
	
	if (!PyArray_Check(op_in) && PyObject_HasAttrString(op_in, "__array__")) {
		PyObject *ap, *arglist;
		
		if (type == PyArray_NOTYPE) {
			arglist = Py_BuildValue("()");
		} else {
			arglist = Py_BuildValue("(c)", 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 (type == PyArray_NOTYPE) {
		type = PyArray_ObjectType(op, 0);
	}
	
	if (PyArray_Check(op) && (((PyArrayObject *)op)->descr->type_num != PyArray_OBJECT || 
		type == PyArray_OBJECT || type == 'O')) {
		if ((((PyArrayObject *)op)->descr->type_num == type ||
			((PyArrayObject *)op)->descr->type == type)) {
			if ((flags & ByCopy) || (flags&Contiguous && !ISCONTIGUOUS((PyArrayObject *)op))) {
				r = PyArray_Copy((PyArrayObject *)op);
			} else {
				Py_INCREF(op);
				r = op;
			}
		} else {
			if (type > PyArray_NTYPES) type = PyArray_DescrFromType(type)->type_num;
			if (PyArray_CanCastSafely(((PyArrayObject *)op)->descr->type_num, type))
				r = PyArray_Cast((PyArrayObject *)op, type);
			else {
				PyErr_SetString(PyExc_TypeError, "Array can not be safely cast to required type");
				r = NULL;
			}
		}
	} else {
		r = Array_FromSequence(op, type, min_depth,max_depth);
		if (r == NULL && min_depth <= 0) {
				PyErr_Clear();
				r = PyArray_FromScalar(op, type);
		}
	}
	Py_DECREF(op);
	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;
}

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:
		/*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 *));
	size = ap->descr->elsize * ap->dimensions[1];
	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_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_FromDimsAndData(n, dimensions, 
		self->descr->type_num,
		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;
}

#define MAX_DIMS 40
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;
}