Update C-API to add features needed for numarray compatibility. Output argument added for several functions and clipmode argument added for a few others.

17 files changed, 1660 insertions, 469 deletions

diff --git a/numpy/core/code_generators/multiarray_api_order.txt b/numpy/core/code_generators/multiarray_api_order.txt
index 10bd1219e..cadab138c 100644
--- a/numpy/core/code_generators/multiarray_api_order.txt
+++ b/numpy/core/code_generators/multiarray_api_order.txt
@@ -1,6 +1,6 @@
 PyArray_Transpose
-PyArray_Take
-PyArray_Put
+PyArray_TakeOut
+PyArray_PutIn
 PyArray_PutMask
 PyArray_Repeat
 PyArray_Choose
@@ -73,3 +73,5 @@ PyArray_RegisterCanCast
 PyArray_InitArrFuncs
 PyArray_IntTupleFromIntp
 PyArray_TypeNumFromName
+PyArray_ClipmodeConverter
+PyArray_OutputConverter
diff --git a/numpy/core/defmatrix.py b/numpy/core/defmatrix.py
index f562cba06..5b47ce6e3 100644
--- a/numpy/core/defmatrix.py
+++ b/numpy/core/defmatrix.py
@@ -234,45 +234,45 @@ class matrix(N.ndarray):
             raise ValueError, "unsupported axis"
 
     # To preserve orientation of result...
-    def sum(self, axis=None, dtype=None):
+    def sum(self, axis=None, dtype=None, out=None):
         """Sum the matrix over the given axis.  If the axis is None, sum
         over all dimensions.  This preserves the orientation of the
         result as a row or column.
         """
-        return N.ndarray.sum(self, axis, dtype)._align(axis)
+        return N.ndarray.sum(self, axis, dtype, out=None)._align(axis)
 
-    def mean(self, axis=None):
-        return N.ndarray.mean(self, axis)._align(axis)
+    def mean(self, axis=None, out=None):
+        return N.ndarray.mean(self, axis, out)._align(axis)
 
-    def std(self, axis=None, dtype=None):
-        return N.ndarray.std(self, axis, dtype)._align(axis)
+    def std(self, axis=None, dtype=None, out=None):
+        return N.ndarray.std(self, axis, dtype, out)._align(axis)
 
-    def var(self, axis=None, dtype=None):
-        return N.ndarray.var(self, axis, dtype)._align(axis)
+    def var(self, axis=None, dtype=None, out=None):
+        return N.ndarray.var(self, axis, dtype, out)._align(axis)
     
-    def prod(self, axis=None, dtype=None):
-        return N.ndarray.prod(self, axis, dtype)._align(axis)
+    def prod(self, axis=None, dtype=None, out=None):
+        return N.ndarray.prod(self, axis, dtype, out)._align(axis)
     
-    def any(self, axis=None):
-        return N.ndarray.any(self, axis)._align(axis)        
+    def any(self, axis=None, out=None):
+        return N.ndarray.any(self, axis, out)._align(axis)        
 
-    def all(self, axis=None):
-        return N.ndarray.all(self, axis)._align(axis)
+    def all(self, axis=None, out=None):
+        return N.ndarray.all(self, axis, out)._align(axis)
         
-    def max(self, axis=None):
-        return N.ndarray.max(self, axis)._align(axis)
+    def max(self, axis=None, out=None):
+        return N.ndarray.max(self, axis, out)._align(axis)
 
     def argmax(self, axis=None):
         return N.ndarray.argmax(self, axis)._align(axis)
     
-    def min(self, axis=None):
-        return N.ndarray.min(self, axis)._align(axis)
+    def min(self, axis=None, out=None):
+        return N.ndarray.min(self, axis, out)._align(axis)
     
     def argmin(self, axis=None):
         return N.ndarray.argmin(self, axis)._align(axis)
     
-    def ptp(self, axis=None):
-        return N.ndarray.ptp(self, axis)._align(axis)
+    def ptp(self, axis=None, out=None):
+        return N.ndarray.ptp(self, axis, out)._align(axis)
 
     # Needed becase tolist method expects a[i]
     #  to have dimension a.ndim-1
diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py
index 5c21edd3b..56d6b86f3 100644
--- a/numpy/core/fromnumeric.py
+++ b/numpy/core/fromnumeric.py
@@ -43,12 +43,12 @@ def _wrapit(obj, method, *args, **kwds):
         result = wrap(result)
     return result
 
-def take(a, indices, axis=None):
+def take(a, indices, axis=None, out=None, mode='raise'):
     try:
         take = a.take
     except AttributeError:
-        return _wrapit(a, 'take', indices, axis)
-    return take(indices, axis)
+        return _wrapit(a, 'take', indices, axis, out, mode)
+    return take(indices, axis, out, mode)
 
 # not deprecated --- copy if necessary, view otherwise
 def reshape(a, newshape, order='C'):
@@ -61,12 +61,12 @@ def reshape(a, newshape, order='C'):
         return _wrapit(a, 'reshape', newshape, order=order)
     return reshape(newshape, order=order)
 
-def choose(a, choices):
+def choose(a, choices, out=None, mode='raise'):
     try:
         choose = a.choose
     except AttributeError:
-        return _wrapit(a, 'choose', choices)
-    return choose(choices)
+        return _wrapit(a, 'choose', choices, out=out, mode=mode)
+    return choose(choices, out=out, mode=mode)
 
 def repeat(a, repeats, axis=None):
     """repeat elements of a repeats times along axis
@@ -82,7 +82,7 @@ def repeat(a, repeats, axis=None):
         return _wrapit(a, 'repeat', repeats, axis)
     return repeat(repeats, axis)
 
-def put (a, ind, v):
+def put (a, ind, v, mode='raise'):
     """put(a, ind, v) results in a[n] = v[n] for all n in ind
        If v is shorter than mask it will be repeated as necessary.
        In particular v can be a scalar or length 1 array.
@@ -94,7 +94,7 @@ def put (a, ind, v):
            for i in ind: a.flat[i] = v[i]
        a must be a contiguous numpy array.
     """
-    return a.put(v,ind)
+    return a.put(v,ind, mode)
 
 def putmask (a, mask, v):
     """putmask(a, mask, v) results in a = v for all places mask is true.
@@ -124,22 +124,22 @@ def transpose(a, axes=None):
         return _wrapit(a, 'transpose', axes)
     return transpose(axes)
 
-def sort(a, axis=-1):
+def sort(a, axis=-1, kind='quicksort'):
     """sort(a,axis=-1) returns array with elements sorted along given axis.
     """
     a = asanyarray(a).copy()
-    a.sort(axis)
+    a.sort(axis, kind)
     return a
 
-def argsort(a, axis=-1):
+def argsort(a, axis=-1, kind='quicksort'):
     """argsort(a,axis=-1) return the indices into a of the sorted array
     along the given axis.
     """
     try:
         argsort = a.argsort
     except AttributeError:
-        return _wrapit(a, 'argsort', axis)
-    return argsort(axis)
+        return _wrapit(a, 'argsort', axis, kind)
+    return argsort(axis, kind)
 
 def argmax(a, axis=-1):
     """argmax(a,axis=-1) returns the indices to the maximum value of the
@@ -215,11 +215,11 @@ def diagonal(a, offset=0, axis1=0, axis2=1):
     """
     return asarray(a).diagonal(offset, axis1, axis2)
 
-def trace(a, offset=0, axis1=0, axis2=1, dtype=None):
+def trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):
     """trace(a,offset=0, axis1=0, axis2=1) returns the sum along diagonals
     (defined by the last two dimenions) of the array.
     """
-    return asarray(a).trace(offset, axis1, axis2, dtype)
+    return asarray(a).trace(offset, axis1, axis2, dtype, out)
 
 def ravel(m,order='C'):
     """ravel(m) returns a 1d array corresponding to all the elements of it's
@@ -250,15 +250,15 @@ def shape(a):
         result = asarray(a).shape
     return result
 
-def compress(condition, m, axis=-1):
+def compress(condition, m, axis=-1, out=None):
     """compress(condition, x, axis=-1) = those elements of x corresponding
     to those elements of condition that are "true".  condition must be the
     same size as the given dimension of x."""
     try:
         compress = m.compress
     except AttributeError:
-        return _wrapit(m, 'compress', condition, axis)
-    return compress(condition, axis)
+        return _wrapit(m, 'compress', condition, axis, out)
+    return compress(condition, axis, out)
 
 def clip(m, m_min, m_max):
     """clip(m, m_min, m_max) = every entry in m that is less than m_min is
@@ -271,7 +271,7 @@ def clip(m, m_min, m_max):
         return _wrapit(m, 'clip', m_min, m_max)
     return clip(m_min, m_max)
 
-def sum(x, axis=None, dtype=None):
+def sum(x, axis=None, dtype=None, out=None):
     """Sum the array over the given axis.  The optional dtype argument
     is the data type for intermediate calculations.
 
@@ -293,97 +293,100 @@ def sum(x, axis=None, dtype=None):
     array([1, 5])
     """
     if isinstance(x, _gentype):
-        return _sum_(x)
+        res = _sum_(x)
+        if out is not None:
+            out[...] = res
+            return out
     try:
         sum = x.sum
     except AttributeError:
-        return _wrapit(x, 'sum', axis, dtype)
-    return sum(axis, dtype)
+        return _wrapit(x, 'sum', axis, dtype, out)
+    return sum(axis, dtype, out)
 
-def product (x, axis=None, dtype=None):
+def product (x, axis=None, dtype=None, out=None):
     """Product of the array elements over the given axis."""
     try:
         prod = x.prod
     except AttributeError:
-        return _wrapit(x, 'prod', axis, dtype)
-    return prod(axis, dtype)
+        return _wrapit(x, 'prod', axis, dtype, out)
+    return prod(axis, dtype, out)
 
-def sometrue (x, axis=None):
+def sometrue (x, axis=None, out=None):
     """Perform a logical_or over the given axis."""
     try:
         any = x.any
     except AttributeError:
-        return _wrapit(x, 'any', axis)
-    return any(axis)
+        return _wrapit(x, 'any', axis, out)
+    return any(axis, out)
 
-def alltrue (x, axis=None):
+def alltrue (x, axis=None, out=None):
     """Perform a logical_and over the given axis."""
     try:
         all = x.all
     except AttributeError:
-        return _wrapit(x, 'all', axis)
-    return all(axis)
+        return _wrapit(x, 'all', axis, out)
+    return all(axis, out)
 
-def any(x,axis=None):
+def any(x,axis=None, out=None):
     """Return true if any elements of x are true:  
     """
     try:
         any = x.any
     except AttributeError:
-        return _wrapit(x, 'any', axis)
-    return any(axis)
+        return _wrapit(x, 'any', axis, out)
+    return any(axis, out)
 
-def all(x,axis=None):
+def all(x,axis=None, out=None):
     """Return true if all elements of x are true:  
     """
     try:
         all = x.all
     except AttributeError:
-        return _wrapit(x, 'all', axis)
-    return all(axis)
+        return _wrapit(x, 'all', axis, out)
+    return all(axis, out)
 
-def cumsum (x, axis=None, dtype=None):
+def cumsum (x, axis=None, dtype=None, out=None):
     """Sum the array over the given axis."""
     try:
         cumsum = x.cumsum
     except AttributeError:
-        return _wrapit(x, 'cumsum', axis, dtype)
-    return cumsum(axis, dtype)
+        return _wrapit(x, 'cumsum', axis, dtype, out)
+    return cumsum(axis, dtype, out)
 
-def cumproduct (x, axis=None, dtype=None):
+def cumproduct (x, axis=None, dtype=None, out=None):
     """Sum the array over the given axis."""
     try:
         cumprod = x.cumprod
     except AttributeError:
-        return _wrapit(x, 'cumprod', axis, dtype)
-    return cumprod(axis, dtype)
+        return _wrapit(x, 'cumprod', axis, dtype, out)
+    return cumprod(axis, dtype, out)
 
-def ptp(a, axis=None):
+def ptp(a, axis=None, out=None):
     """Return maximum - minimum along the the given dimension
     """
     try:
         ptp = a.ptp
     except AttributeError:
-        return _wrapit(a, 'ptp', axis)
-    return ptp(axis)
+        return _wrapit(a, 'ptp', axis, out)
+    return ptp(axis, out)
 
-def amax(a, axis=None):
+def amax(a, axis=None, out=None):
     """Return the maximum of 'a' along dimension axis.
     """
     try:
         max = a.max
     except AttributeError:
-        return _wrapit(a, 'max', axis)
-    return max(axis)
+        return _wrapit(a, 'max', axis, out)
+    return max(axis, out)
 
-def amin(a, axis=None):
+def amin(a, axis=None, out=None):
     """Return the minimum of a along dimension axis.
     """
     try:
         min = a.min
     except AttributeError:
-        return _wrapit(a, 'min', axis)
-    return min(axis)
+        return _wrapit(a, 'min', axis, out)
+    return min(axis, out)
 
 def alen(a):
     """Return the length of a Python object interpreted as an array
@@ -394,23 +397,23 @@ def alen(a):
     except TypeError:
         return len(array(a,ndmin=1))
 
-def prod(a, axis=None, dtype=None):
+def prod(a, axis=None, dtype=None, out=None):
     """Return the product of the elements along the given axis
     """
     try:
         prod = a.prod
     except AttributeError:
-        return _wrapit(a, 'prod', axis, dtype)
-    return prod(axis, dtype)
+        return _wrapit(a, 'prod', axis, dtype, out)
+    return prod(axis, dtype, out)
 
-def cumprod(a, axis=None, dtype=None):
+def cumprod(a, axis=None, dtype=None, out=None):
     """Return the cumulative product of the elments along the given axis
     """
     try:
         cumprod = a.cumprod
     except AttributeError:
-        return _wrapit(a, 'cumprod', axis, dtype)
-    return cumprod(axis, dtype)
+        return _wrapit(a, 'cumprod', axis, dtype, out)
+    return cumprod(axis, dtype, out)
 
 def ndim(a):
     try:
@@ -440,7 +443,7 @@ def size (a, axis=None):
         except AttributeError:
             return asarray(a).shape[axis]
 
-def round_(a, decimals=0):
+def round_(a, decimals=0, out=None):
     """Round 'a' to the given number of decimal places.  Rounding
     behaviour is equivalent to Python.
 
@@ -450,12 +453,12 @@ def round_(a, decimals=0):
     try:
         round = a.round
     except AttributeError:
-        return _wrapit(a, 'round', decimals)
-    return round(decimals)
+        return _wrapit(a, 'round', decimals, out)
+    return round(decimals, out)
 
 around = round_
 
-def mean(a, axis=None, dtype=None):
+def mean(a, axis=None, dtype=None, out=None):
     """mean(a, axis=None, dtype=None)
     Return the arithmetic mean.
     
@@ -466,10 +469,10 @@ def mean(a, axis=None, dtype=None):
     try:
         mean = a.mean
     except AttributeError:
-        return _wrapit(a, 'mean', axis, dtype)
-    return mean(axis, dtype)
+        return _wrapit(a, 'mean', axis, dtype, out)
+    return mean(axis, dtype, out)
 
-def std(a, axis=None, dtype=None):
+def std(a, axis=None, dtype=None, out=None):
     """std(sample, axis=None, dtype=None)
     Return the standard deviation, a measure of the spread of a distribution.
 
@@ -481,10 +484,10 @@ def std(a, axis=None, dtype=None):
     try:
         std = a.std
     except AttributeError:
-        return _wrapit(a, 'std', axis, dtype)
-    return std(axis, dtype)
+        return _wrapit(a, 'std', axis, dtype, out)
+    return std(axis, dtype, out)
 
-def var(a, axis=None, dtype=None):
+def var(a, axis=None, dtype=None, out=None):
     """var(sample, axis=None, dtype=None)
     Return the variance, a measure of the spread of a distribution.
 
@@ -496,5 +499,5 @@ def var(a, axis=None, dtype=None):
     try:
         var = a.var
     except AttributeError:
-        return _wrapit(a, 'var', axis, dtype)
-    return var(axis, dtype)
+        return _wrapit(a, 'var', axis, dtype, out)
+    return var(axis, dtype, out)
diff --git a/numpy/core/include/numpy/arrayobject.h b/numpy/core/include/numpy/arrayobject.h
index c70271d87..f4fc725d4 100644
--- a/numpy/core/include/numpy/arrayobject.h
+++ b/numpy/core/include/numpy/arrayobject.h
@@ -37,7 +37,7 @@ extern "C" CONFUSE_EMACS
 #define NPY_SUCCEED 1
 
         /* Helpful to distinguish what is installed */
-#define NPY_VERSION 0x01000000
+#define NPY_VERSION 0x01000001
 
 	/* Some platforms don't define bool, long long, or long double.
 	   Handle that here.
@@ -207,6 +207,12 @@ typedef enum {
 } NPY_ORDER;
 
 
+typedef enum {
+        NPY_CLIP=0,
+        NPY_WRAP=1,
+        NPY_RAISE=2
+} NPY_CLIPMODE;
+ 
 	/* Define bit-width array types and typedefs */
 
 #define NPY_MAX_INT8 127
@@ -1526,6 +1532,12 @@ typedef struct {
 #define PyArray_Cast(mp, type_num) \
 	PyArray_CastToType(mp, PyArray_DescrFromType(type_num), 0)
 
+#define PyArray_Take(ap, items, axis) \
+        PyArray_TakeOut(ap, items, axis, NULL, NPY_RAISE)
+
+#define PyArray_Put(ap, items, values) \
+        PyArray_PutIn(ap, items, values, NPY_RAISE)
+
 /* Compatibility with old Numeric stuff -- don't use in new code */
 
 #define PyArray_FromDimsAndData(nd, d, type, data) \
diff --git a/numpy/core/include/numpy/ufuncobject.h b/numpy/core/include/numpy/ufuncobject.h
index ab53c6afe..2446f5431 100644
--- a/numpy/core/include/numpy/ufuncobject.h
+++ b/numpy/core/include/numpy/ufuncobject.h
@@ -158,6 +158,7 @@ typedef struct {
         PyObject *decref;
 
         int obj;
+        int retbase;
 
 } PyUFuncReduceObject;
 
diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index 36dcb0ab6..d4537f435 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -11,14 +11,15 @@ __all__ = ['newaxis', 'ndarray', 'flatiter', 'ufunc',
            'array2string', 'get_printoptions', 'set_printoptions',
            'array_repr', 'array_str', 'set_string_function',
            'little_endian', 'require',
-           'fromiter',
+           'fromiter', 'array_equal', 'array_equiv',
            'indices', 'fromfunction',
            'load', 'loads', 'isscalar', 'binary_repr', 'base_repr',
            'ones', 'identity', 'allclose',
            'seterr', 'geterr', 'setbufsize', 'getbufsize',
            'seterrcall', 'geterrcall', 'flatnonzero',
            'Inf', 'inf', 'infty', 'Infinity',
-           'nan', 'NaN', 'False_', 'True_', 'bitwise_not']
+           'nan', 'NaN', 'False_', 'True_', 'bitwise_not',
+           'CLIP', 'RAISE', 'WRAP', 'MAXDIMS', 'BUFSIZE', 'ALLOW_THREADS']
 
 import sys
 import multiarray
@@ -29,6 +30,14 @@ from numerictypes import *
 
 bitwise_not = invert
 
+CLIP = multiarray.CLIP
+WRAP = multiarray.WRAP
+RAISE = multiarray.RAISE
+MAXDIMS = multiarray.MAXDIMS
+ALLOW_THREADS = multiarray.ALLOW_THREADS
+BUFSIZE = multiarray.BUFSIZE
+
+
 # from Fernando Perez's IPython
 def zeros_like(a):
     """Return an array of zeros of the shape and typecode of a.
@@ -475,6 +484,25 @@ def allclose (a, b, rtol=1.e-5, atol=1.e-8):
     return d.ravel().all()
 
 
+def array_equal(a1, a2):
+    try:
+        a1, a2 = asarray(a1), asarray(a2)
+    except:
+        return 0
+    if a1.shape != a2.shape:
+        return 0
+    return logical_and.reduce(equal(a1,a2).ravel())
+
+def array_equiv(a1, a2):
+    try:
+        a1, a2 = asarray(a1), asarray(a2)
+    except:
+        return 0
+    try:
+        return logical_and.reduce(equal(a1,a2).ravel())
+    except ValueError:
+        return 0
+
 
 _errdict = {"ignore":ERR_IGNORE,
             "warn":ERR_WARN,
diff --git a/numpy/core/src/arraymethods.c b/numpy/core/src/arraymethods.c
index cc885b90f..75f0d5b9b 100644
--- a/numpy/core/src/arraymethods.c
+++ b/numpy/core/src/arraymethods.c
@@ -10,15 +10,21 @@ array_take(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int dimension=MAX_DIMS;
 	PyObject *indices;
-	static char *kwlist[] = {"indices", "axis", NULL};
+        PyArrayObject *out=NULL;
+        NPY_CLIPMODE mode=NPY_RAISE;
+	static char *kwlist[] = {"indices", "axis", "out", "mode", NULL};
 	
 	dimension=0;
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|O&O&O&", kwlist, 
 					 &indices, PyArray_AxisConverter,
-					 &dimension))
+					 &dimension,
+                                         PyArray_OutputConverter,
+                                         &out,
+                                         PyArray_ClipmodeConverter,
+                                         &mode))
 		return NULL;
 	
-	return _ARET(PyArray_Take(self, indices, dimension));
+	return _ARET(PyArray_TakeOut(self, indices, dimension, out, mode));
 }
 
 static char doc_fill[] = "a.fill(value) places the scalar value at every "\
@@ -35,7 +41,7 @@ array_fill(PyArrayObject *self, PyObject *args)
 	return Py_None;
 }
 
-static char doc_put[] = "a.put(values, indices) sets a.flat[n] = v[n] "\
+static char doc_put[] = "a.put(values, indices, mode) sets a.flat[n] = v[n] "\
 	"for each n in indices. v can be scalar or shorter than indices, "\
 	"will repeat.";
 
@@ -43,12 +49,15 @@ static PyObject *
 array_put(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	PyObject *indices, *values;
-	static char *kwlist[] = {"values", "indices", NULL};
+        NPY_CLIPMODE mode=NPY_RAISE;
+	static char *kwlist[] = {"values", "indices", "mode", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO", kwlist,
-					 &values, &indices))
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO|O&", kwlist,
+					 &values, &indices,
+                                         PyArray_ClipmodeConverter,
+                                         &mode))
 		return NULL;
-	return PyArray_Put(self, values, indices);
+	return PyArray_PutIn(self, values, indices, mode);
 }
 
 static char doc_putmask[] = "a.putmask(values, mask) sets a.flat[n] = v[n] "\
@@ -185,14 +194,17 @@ static PyObject *
 array_max(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int axis=MAX_DIMS;
-	static char *kwlist[] = {"axis", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
 					 PyArray_AxisConverter,
-					 &axis))
+					 &axis,
+                                         PyArray_OutputConverter,
+                                         &out))
 		return NULL;	
 	
-	return PyArray_Max(self, axis);
+	return PyArray_Max(self, axis, out);
 }
 
 static char doc_ptp[] = "a.ptp(axis=None) a.max(axis)-a.min(axis)";
@@ -201,14 +213,17 @@ static PyObject *
 array_ptp(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int axis=MAX_DIMS;
-	static char *kwlist[] = {"axis", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
 					 PyArray_AxisConverter,
-					 &axis))
+					 &axis,
+                                         PyArray_OutputConverter,
+                                         &out))
 		return NULL;	
 		
-	return PyArray_Ptp(self, axis);
+	return PyArray_Ptp(self, axis, out);
 }
 
 
@@ -218,14 +233,17 @@ static PyObject *
 array_min(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int axis=MAX_DIMS;
-	static char *kwlist[] = {"axis", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
 					 PyArray_AxisConverter,
-					 &axis))
+					 &axis,
+                                         PyArray_OutputConverter,
+                                         &out))
 		return NULL;	
 	
-	return PyArray_Min(self, axis);
+	return PyArray_Min(self, axis, out);
 }
 
 static char doc_swapaxes[] = "a.swapaxes(axis1, axis2)  returns new view with axes swapped.";
@@ -451,7 +469,7 @@ array_tofile(PyArrayObject *self, PyObject *args, PyObject *kwds)
 	if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|ss", kwlist, 
                                          &file, &sep, &format)) return NULL;
 
-	if (PyString_Check(file)) {
+	if (PyString_Check(file) || PyUnicode_Check(file)) {
 		file = PyObject_CallFunction((PyObject *)&PyFile_Type,
 					     "Os", file, "wb");
 		if (file==NULL) return NULL;
@@ -610,7 +628,7 @@ array_copy(PyArrayObject *self, PyObject *args)
         return PyArray_NewCopy(self, fortran);
 }
 
-static char doc_resize[] = "self.resize(new_shape, refcheck=True, fortran=False).  "\
+static char doc_resize[] = "self.resize(new_shape, refcheck=True, order=False).  "\
 	"Change size and shape of self inplace.\n"\
 	"\n    Array must own its own memory and not be referenced by other " \
 	"arrays\n    Returns None.";
@@ -633,7 +651,7 @@ array_resize(PyArrayObject *self, PyObject *args, PyObject *kwds)
 				return NULL;
 			}
 		}
-		ref = PyDict_GetItemString(kwds, "fortran");
+		ref = PyDict_GetItemString(kwds, "order");
 		if (ref != NULL || 
 		    (PyArray_OrderConverter(ref, &fortran) == PY_FAIL))
 			return NULL;
@@ -678,7 +696,7 @@ array_repeat(PyArrayObject *self, PyObject *args, PyObject *kwds) {
 	return _ARET(PyArray_Repeat(self, repeats, axis));
 }
 
-static char doc_choose[] = "a.choose(b0, b1, ..., bn)\n"\
+static char doc_choose[] = "a.choose(b0, b1, ..., bn, out=None, mode='raise')\n"\
 	"\n"                                                            \
 	"Return an array that merges the b_i arrays together using 'a' as the index\n"
         "The b_i arrays and 'a' must all be broadcastable to the same shape.\n"
@@ -687,10 +705,12 @@ static char doc_choose[] = "a.choose(b0, b1, ..., bn)\n"\
         "an integer array with entries from 0 to n+1.";
 
 static PyObject *
-array_choose(PyArrayObject *self, PyObject *args) 
+array_choose(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	PyObject *choices;
 	int n;
+        PyArrayObject *out=NULL;
+        NPY_CLIPMODE clipmode=NPY_RAISE;
 	
 	n = PyTuple_Size(args);
 	if (n <= 1) {
@@ -700,8 +720,18 @@ array_choose(PyArrayObject *self, PyObject *args)
         else {
 		choices = args;
 	}
+        if (kwds && PyDict_Check(kwds)) {
+                if (PyArray_OutputConverter(PyDict_GetItemString(kwds, 
+                                                                 "out"),
+                                            &out) == PY_FAIL)
+                        return NULL;
+                if (PyArray_ClipmodeConverter(PyDict_GetItemString(kwds, 
+                                                                   "mode"), 
+                                              &clipmode) == PY_FAIL)
+                        return NULL;
+        }
 	
-	return _ARET(PyArray_Choose(self, choices));
+	return _ARET(PyArray_Choose(self, choices, out, clipmode));
 }
 
 static char doc_sort[] = "a.sort(axis=-1,kind='quicksort') sorts in place along axis.  Return is None and kind can be 'quicksort', 'mergesort', or 'heapsort'";
@@ -1254,14 +1284,17 @@ array_mean(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int axis=MAX_DIMS;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"axis", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "dtype", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist,
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist,
 					 PyArray_AxisConverter, 
 					 &axis, PyArray_DescrConverter2,
-					 &dtype)) return NULL;
+					 &dtype,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 
-	return PyArray_Mean(self, axis, _CHKTYPENUM(dtype));
+	return PyArray_Mean(self, axis, _CHKTYPENUM(dtype), out);
 }
 
 static char doc_sum[] = "a.sum(axis=None, dtype=None)\n\n"\
@@ -1293,32 +1326,38 @@ array_sum(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int axis=MAX_DIMS;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"axis", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "dtype", "out", NULL};
 
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist, 
 					 PyArray_AxisConverter, 
 					 &axis, PyArray_DescrConverter2,
-					 &dtype)) return NULL;
+					 &dtype,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 	
-	return PyArray_Sum(self, axis, _CHKTYPENUM(dtype));
+	return PyArray_Sum(self, axis, _CHKTYPENUM(dtype), out);
 }
 
 
-static char doc_cumsum[] = "a.cumsum(axis=None, dtype=None)";
+static char doc_cumsum[] = "a.cumsum(axis=None, dtype=None, out=None)";
 
 static PyObject *
 array_cumsum(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int axis=MAX_DIMS;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"axis", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "dtype", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist, 
 					 PyArray_AxisConverter, 
 					 &axis, PyArray_DescrConverter2,
-					 &dtype)) return NULL;
+					 &dtype,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 	
-	return PyArray_CumSum(self, axis, _CHKTYPENUM(dtype));
+	return PyArray_CumSum(self, axis, _CHKTYPENUM(dtype), out);
 }
 
 static char doc_prod[] = "a.prod(axis=None, dtype=None)";
@@ -1328,14 +1367,17 @@ array_prod(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int axis=MAX_DIMS;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"axis", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "dtype", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist, 
 					 PyArray_AxisConverter, 
 					 &axis, PyArray_DescrConverter2,
-					 &dtype)) return NULL;
+					 &dtype,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 	
-	return PyArray_Prod(self, axis, _CHKTYPENUM(dtype));
+	return PyArray_Prod(self, axis, _CHKTYPENUM(dtype), out);
 }
 
 
@@ -1346,31 +1388,37 @@ array_cumprod(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int axis=MAX_DIMS;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"axis", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "dtype", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist, 
 					 PyArray_AxisConverter, 
 					 &axis, PyArray_DescrConverter2,
-					 &dtype)) return NULL;
+					 &dtype,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 	
-	return PyArray_CumProd(self, axis, _CHKTYPENUM(dtype));
+	return PyArray_CumProd(self, axis, _CHKTYPENUM(dtype), out);
 }
 
 
-static char doc_any[] = "a.any(axis=None)";
+static char doc_any[] = "a.any(axis=None, out=None)";
 
 static PyObject *
 array_any(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int axis=MAX_DIMS;
-	static char *kwlist[] = {"axis", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
 					 PyArray_AxisConverter,
-					 &axis))
+					 &axis,
+                                         PyArray_OutputConverter,
+                                         &out))
 		return NULL;	
 	
-	return PyArray_Any(self, axis);
+	return PyArray_Any(self, axis, out);
 }
 
 static char doc_all[] = "a.all(axis=None)";
@@ -1379,17 +1427,20 @@ static PyObject *
 array_all(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int axis=MAX_DIMS;
-	static char *kwlist[] = {"axis", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
 					 PyArray_AxisConverter,
-					 &axis))
+					 &axis,
+                                         PyArray_OutputConverter,
+                                         &out))
 		return NULL;	
 
-	return PyArray_All(self, axis);
+	return PyArray_All(self, axis, out);
 }
 
-static char doc_stddev[] = "a.std(axis=None, dtype=None)\n"
+static char doc_stddev[] = "a.std(axis=None, dtype=None, out=None)\n"
 "Return the standard deviation, a measure of the spread of a distribution.\n"
 "\n"
 "The standard deviation is the square root of the average of the squared\n"
@@ -1402,14 +1453,17 @@ array_stddev(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int axis=MAX_DIMS;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"axis", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "dtype", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist, 
 					 PyArray_AxisConverter, 
 					 &axis, PyArray_DescrConverter2,
-					 &dtype)) return NULL;
+					 &dtype,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 	
-	return PyArray_Std(self, axis, _CHKTYPENUM(dtype), 0);
+	return PyArray_Std(self, axis, _CHKTYPENUM(dtype), out, 0);
 }
 
 static char doc_variance[] = "a.var(axis=None, dtype=None)";
@@ -1419,30 +1473,36 @@ array_variance(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int axis=MAX_DIMS;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"axis", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"axis", "dtype", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist, 
 					 PyArray_AxisConverter, 
 					 &axis, PyArray_DescrConverter2,
-					 &dtype)) return NULL;
+					 &dtype,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 	
-	return PyArray_Std(self, axis, _CHKTYPENUM(dtype), 1);
+	return PyArray_Std(self, axis, _CHKTYPENUM(dtype), out, 1);
 }
 
-static char doc_compress[] = "a.compress(condition=, axis=None)";
+static char doc_compress[] = "a.compress(condition=, axis=None, out=None)";
 
 static PyObject *
 array_compress(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int axis=MAX_DIMS;
 	PyObject *condition;	
-	static char *kwlist[] = {"condition", "axis", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"condition", "axis", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|O&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|O&O&", kwlist, 
 					 &condition, PyArray_AxisConverter,
-					 &axis)) return NULL;
+					 &axis,
+                                         PyArray_OutputConverter,
+                                         &out)) return NULL;
 
-	return _ARET(PyArray_Compress(self, condition, axis));
+	return _ARET(PyArray_Compress(self, condition, axis, out));
 }
 
 static char doc_nonzero[] = \
@@ -1465,7 +1525,7 @@ array_nonzero(PyArrayObject *self, PyObject *args)
 }
 
 
-static char doc_trace[] = "a.trace(offset=0, axis1=0, axis2=1, dtype=None)\n"\
+static char doc_trace[] = "a.trace(offset=0, axis1=0, axis2=1, dtype=None, out=None)\n"\
 	"return the sum along the offset diagonal of the arrays indicated\n" \
 	"axis1 and axis2.";
 
@@ -1474,15 +1534,17 @@ array_trace(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
 	int axis1=0, axis2=1, offset=0;
 	PyArray_Descr *dtype=NULL;
-	static char *kwlist[] = {"offset", "axis1", "axis2", "dtype", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"offset", "axis1", "axis2", "dtype", "out", NULL};
 	
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|iiiO&", kwlist, 
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|iiiO&O&", kwlist, 
 					 &offset, &axis1, &axis2,
-					 PyArray_DescrConverter2, &dtype))
+					 PyArray_DescrConverter2, &dtype,
+                                         PyArray_OutputConverter, &out))
 		return NULL;
 	
 	return _ARET(PyArray_Trace(self, offset, axis1, axis2, 
-				   _CHKTYPENUM(dtype)));
+				   _CHKTYPENUM(dtype), out));
 }
 
 #undef _CHKTYPENUM
@@ -1558,19 +1620,21 @@ array_ravel(PyArrayObject *self, PyObject *args)
 	return PyArray_Ravel(self, fortran);
 }
 
-static char doc_round[] = "a.round(decimals=0)";
+static char doc_round[] = "a.round(decimals=0, out=None)";
 
 static PyObject *
 array_round(PyArrayObject *self, PyObject *args, PyObject *kwds) 
 {
 	int decimals = 0;
-	static char *kwlist[] = {"decimals", NULL};
+        PyArrayObject *out=NULL;
+	static char *kwlist[] = {"decimals", "out", NULL};
 
-	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|i", kwlist,
-					 &decimals)) 
-		return NULL;
-	
-	return _ARET(PyArray_Round(self, decimals));
+	if (!PyArg_ParseTupleAndKeywords(args, kwds, "|iO&", kwlist,
+					 &decimals, PyArray_OutputConverter,
+                                         &out))
+            return NULL;
+        	
+	return _ARET(PyArray_Round(self, decimals, out));
 }
 
 
@@ -1702,7 +1766,7 @@ static PyMethodDef array_methods[] = {
 	{"repeat",	(PyCFunction)array_repeat, 
 	 METH_VARARGS|METH_KEYWORDS, doc_repeat},
 	{"choose",	(PyCFunction)array_choose, 
-	 METH_VARARGS, doc_choose},	
+	 METH_VARARGS|METH_KEYWORDS, doc_choose},	
 	{"sort",	(PyCFunction)array_sort, 
 	 METH_VARARGS|METH_KEYWORDS, doc_sort},
 	{"argsort",	(PyCFunction)array_argsort, 
diff --git a/numpy/core/src/arrayobject.c b/numpy/core/src/arrayobject.c
index f17c2098d..5b19dcdca 100644
--- a/numpy/core/src/arrayobject.c
+++ b/numpy/core/src/arrayobject.c
@@ -3218,55 +3218,70 @@ PyArray_GetNumericOps(void)
 }
 
 static PyObject *
+_get_keywords(int rtype, PyArrayObject *out)
+{
+        PyObject *kwds=NULL;
+        if (rtype != PyArray_NOTYPE || out != NULL) {
+                kwds = PyDict_New();
+                if (rtype != PyArray_NOTYPE) {
+                        PyArray_Descr *descr;
+                        descr = PyArray_DescrFromType(rtype);
+                        if (descr) {
+                                PyDict_SetItemString(kwds, "dtype", 
+                                                     (PyObject *)descr);
+                                Py_DECREF(descr);
+                        }
+                }
+                if (out != NULL) {
+                        PyDict_SetItemString(kwds, "out", 
+                                             (PyObject *)out);
+                }
+        }
+        return kwds;
+}
+
+static PyObject *
 PyArray_GenericReduceFunction(PyArrayObject *m1, PyObject *op, int axis,
-                              int rtype)
+                              int rtype, PyArrayObject *out)
 {
         PyObject *args, *ret=NULL, *meth;
+        PyObject *kwds;
         if (op == NULL) {
                 Py_INCREF(Py_NotImplemented);
                 return Py_NotImplemented;
         }
-        if (rtype == PyArray_NOTYPE)
-                args = Py_BuildValue("(Oi)", m1, axis);
-        else {
-                PyArray_Descr *descr;
-                descr = PyArray_DescrFromType(rtype);
-                args = Py_BuildValue("(Oic)", m1, axis, descr->type);
-                Py_DECREF(descr);
-        }
+        args = Py_BuildValue("(Oi)", m1, axis);
+        kwds = _get_keywords(rtype, out);
         meth = PyObject_GetAttrString(op, "reduce");
         if (meth && PyCallable_Check(meth)) {
-                ret = PyObject_Call(meth, args, NULL);
+                ret = PyObject_Call(meth, args, kwds);
         }
         Py_DECREF(args);
         Py_DECREF(meth);
+        Py_XDECREF(kwds);
         return ret;
 }
 
 
 static PyObject *
 PyArray_GenericAccumulateFunction(PyArrayObject *m1, PyObject *op, int axis,
-                                  int rtype)
+                                  int rtype, PyArrayObject *out)
 {
         PyObject *args, *ret=NULL, *meth;
+        PyObject *kwds;
         if (op == NULL) {
                 Py_INCREF(Py_NotImplemented);
                 return Py_NotImplemented;
         }
-        if (rtype == PyArray_NOTYPE)
-                args = Py_BuildValue("(Oi)", m1, axis);
-        else {
-                PyArray_Descr *descr;
-                descr = PyArray_DescrFromType(rtype);
-                args = Py_BuildValue("(Oic)", m1, axis, descr->type);
-                Py_DECREF(descr);
-        }
+        args = Py_BuildValue("(Oi)", m1, axis);
+        kwds = _get_keywords(rtype, out);
         meth = PyObject_GetAttrString(op, "accumulate");
         if (meth && PyCallable_Check(meth)) {
-                ret = PyObject_Call(meth, args, NULL);
+                ret = PyObject_Call(meth, args, kwds);
         }
         Py_DECREF(args);
         Py_DECREF(meth);
+        Py_XDECREF(kwds);
         return ret;
 }
 
diff --git a/numpy/core/src/multiarraymodule.c b/numpy/core/src/multiarraymodule.c
index c33f307f2..7f57c3632 100644
--- a/numpy/core/src/multiarraymodule.c
+++ b/numpy/core/src/multiarraymodule.c
@@ -213,22 +213,36 @@ power_of_ten(int n)
  Round
 */
 static PyObject *
-PyArray_Round(PyArrayObject *a, int decimals)
+PyArray_Round(PyArrayObject *a, int decimals, PyArrayObject *out)
 {
-	PyObject *f, *ret=NULL, *tmp;
+	PyObject *f, *ret=NULL, *tmp, *op1, *op2;
+        if (out && (!PyArray_SAMESHAPE(out, a) || 
+                    !PyArray_EquivTypes(a->descr, out->descr))) {
+                PyErr_SetString(PyExc_ValueError, 
+                                "output array must have the same shape"
+                                "and type");
+                return NULL;
+        }
 	if (PyArray_ISCOMPLEX(a)) {
 		PyObject *part;
 		PyObject *round_part;
 		PyObject *new;
 		int res;
-		new = PyArray_Copy(a);
-		if (new == NULL) return NULL;
+                if (out) {
+                        new = (PyObject *)out;
+                        Py_INCREF(new);
+                }
+		else {
+                        new = PyArray_Copy(a);
+                        if (new == NULL) return NULL;
+                }
 
 		/* new.real = a.real.round(decimals) */
 		part = PyObject_GetAttrString(new, "real");
 		if (part == NULL) {Py_DECREF(new); return NULL;}
 		part = PyArray_EnsureAnyArray(part);
-		round_part = PyArray_Round((PyArrayObject *)part, decimals);
+		round_part = PyArray_Round((PyArrayObject *)part, 
+                                           decimals, NULL);
 		Py_DECREF(part);
 		if (round_part == NULL) {Py_DECREF(new); return NULL;}
 		res = PyObject_SetAttrString(new, "real", round_part);
@@ -239,7 +253,8 @@ PyArray_Round(PyArrayObject *a, int decimals)
 		part = PyObject_GetAttrString(new, "imag");
 		if (part == NULL) {Py_DECREF(new); return NULL;}
 		part = PyArray_EnsureAnyArray(part); 
-		round_part = PyArray_Round((PyArrayObject *)part, decimals);
+		round_part = PyArray_Round((PyArrayObject *)part, 
+                                           decimals, NULL);
 		Py_DECREF(part);
 		if (round_part == NULL) {Py_DECREF(new); return NULL;}
 		res = PyObject_SetAttrString(new, "imag", round_part);
@@ -248,64 +263,50 @@ PyArray_Round(PyArrayObject *a, int decimals)
 		return new;
 	}
 	/* do the most common case first */
-	if (decimals == 0) {
+	if (decimals >= 0) {
 		if (PyArray_ISINTEGER(a)) {
-			Py_INCREF(a);
-			return (PyObject *)a;
-		}
-		return PyArray_GenericUnaryFunction((PyAO *)a, n_ops.rint);
-	}
-	if (decimals > 0) {
-		if (PyArray_ISINTEGER(a)) {
-			Py_INCREF(a);
-			return (PyObject *)a;
-		}
-		f = PyFloat_FromDouble(power_of_ten(decimals));
-		if (f==NULL) return NULL;
-		ret = PyNumber_Multiply((PyObject *)a, f);
-		if (ret==NULL) goto finish;
-		if (PyArray_IsScalar(ret, Generic)) {
-			/* array scalars cannot be modified inplace */
-			tmp = PyObject_CallFunction(n_ops.rint, "O", ret);
-			Py_DECREF(ret);
-			if (tmp == NULL) {ret=NULL; goto finish;}
-			ret = PyObject_CallFunction(n_ops.divide, "OO", 
-						    tmp, f);
-			Py_DECREF(tmp);
-		} else {
-			tmp = PyObject_CallFunction(n_ops.rint, "OO", ret, ret);
-			if (tmp == NULL) {Py_DECREF(ret); ret=NULL; goto finish;}
-			Py_DECREF(tmp);
-			tmp = PyObject_CallFunction(n_ops.divide, "OOO", ret, 
-						    f, ret);
-			if (tmp == NULL) {Py_DECREF(ret); ret=NULL; goto finish;}
-			Py_DECREF(tmp);
-		}
-	} 
-	else {
-		/* remaining case: decimals < 0 */
-		f = PyFloat_FromDouble(power_of_ten(-decimals));
-		if (f==NULL) return NULL;
-		ret = PyNumber_Divide((PyObject *)a, f);
-		if (ret==NULL) goto finish;
-		if (PyArray_IsScalar(ret, Generic)) {
-			/* array scalars cannot be modified inplace */
-			tmp = PyObject_CallFunction(n_ops.rint, "O", ret);
-			Py_DECREF(ret);
-			if (tmp == NULL) {ret=NULL; goto finish;}
-			ret = PyObject_CallFunction(n_ops.multiply, "OO", 
-						    tmp, f);
-			Py_DECREF(tmp);
-		} else {
-			tmp = PyObject_CallFunction(n_ops.rint, "OO", ret, ret);
-			if (tmp == NULL) {Py_DECREF(ret); ret=NULL; goto finish;}
-			Py_DECREF(tmp);
-			tmp = PyObject_CallFunction(n_ops.multiply, "OOO", ret, 
-						    f, ret);
-			if (tmp==NULL) {Py_DECREF(ret); ret=NULL; goto finish;}
-			Py_DECREF(tmp);
+                        if (out) {
+                                if (PyArray_CopyInto(out, a) < 0) return NULL;
+                                Py_INCREF(out);
+                                return (PyObject *)out;
+                        }
+                        else {
+                                Py_INCREF(a);
+                                return (PyObject *)a;
+                        }
 		}
-	}
+                if (decimals == 0) {
+                        if (out) {
+                                return PyObject_CallFunction(n_ops.rint, "OO", 
+                                                             a, out);
+                        }
+                        return PyObject_CallFunction(n_ops.rint, "O", a);
+                }
+                op1 = n_ops.multiply;
+                op2 = n_ops.divide;
+        }
+        else {
+                op1 = n_ops.divide;
+                op2 = n_ops.multiply;
+                decimals = -decimals;
+        }
+        if (!out) {
+                Py_INCREF(a->descr);
+                out = (PyArrayObject *)PyArray_Empty(a->nd, a->dimensions, 
+                                                     a->descr, 
+                                                     PyArray_ISFORTRAN(a));
+                if (out == NULL) return NULL;
+        }
+        f = PyFloat_FromDouble(power_of_ten(decimals));
+        if (f==NULL) return NULL;
+        ret = PyObject_CallFunction(op1, "OOO", a, f, out);
+        if (ret==NULL) goto finish;
+        tmp = PyObject_CallFunction(n_ops.rint, "OO", ret, ret);
+        if (tmp == NULL) {Py_DECREF(ret); ret=NULL; goto finish;}
+        Py_DECREF(tmp);
+        tmp = PyObject_CallFunction(op2, "OOO", ret, f, ret);
+        if (tmp == NULL) {Py_DECREF(ret); ret=NULL; goto finish;}
+        Py_DECREF(tmp);
 
  finish:
 	Py_DECREF(f);
@@ -594,7 +595,7 @@ PyArray_Squeeze(PyArrayObject *self)
  Mean
 */
 static PyObject *
-PyArray_Mean(PyArrayObject *self, int axis, int rtype)
+PyArray_Mean(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 {
 	PyObject *obj1=NULL, *obj2=NULL;
 	PyObject *new, *ret;
@@ -602,7 +603,7 @@ PyArray_Mean(PyArrayObject *self, int axis, int rtype)
 	if ((new = _check_axis(self, &axis, 0))==NULL) return NULL;
 
 	obj1 = PyArray_GenericReduceFunction((PyAO *)new, n_ops.add, axis,
-					     rtype);
+					     rtype, out);
 	obj2 = PyFloat_FromDouble((double) PyArray_DIM(new,axis));
         Py_DECREF(new);
 	if (obj1 == NULL || obj2 == NULL) {
@@ -610,10 +611,14 @@ PyArray_Mean(PyArrayObject *self, int axis, int rtype)
 		Py_XDECREF(obj2);
 		return NULL;
 	}
-
-	ret = PyNumber_Divide(obj1, obj2);
-	Py_DECREF(obj1);
-	Py_DECREF(obj2);
+        if (!out) {
+                ret = PyNumber_Divide(obj1, obj2);
+        }
+        else {
+                ret = PyObject_CallFunction(n_ops.divide, "OOO", out, obj2, out);
+        }
+        Py_DECREF(obj1);
+        Py_DECREF(obj2);
 	return ret;
 }
 
@@ -622,7 +627,8 @@ PyArray_Mean(PyArrayObject *self, int axis, int rtype)
  Std
 */
 static PyObject *
-PyArray_Std(PyArrayObject *self, int axis, int rtype, int variance)
+PyArray_Std(PyArrayObject *self, int axis, int rtype, PyArrayObject *out, 
+            int variance)
 {
 	PyObject *obj1=NULL, *obj2=NULL, *new=NULL;
 	PyObject *ret=NULL, *newshape=NULL;
@@ -630,9 +636,9 @@ PyArray_Std(PyArrayObject *self, int axis, int rtype, int variance)
 	intp val;
 
 	if ((new = _check_axis(self, &axis, 0))==NULL) return NULL;
-	
+
 	/* Compute and reshape mean */
-	obj1 = PyArray_EnsureArray(PyArray_Mean((PyAO *)new, axis, rtype));
+	obj1 = PyArray_EnsureArray(PyArray_Mean((PyAO *)new, axis, rtype, NULL));
 	if (obj1 == NULL) {Py_DECREF(new); return NULL;} 
 	n = PyArray_NDIM(new);
 	newshape = PyTuple_New(n);
@@ -660,7 +666,7 @@ PyArray_Std(PyArrayObject *self, int axis, int rtype, int variance)
 
 	/* Compute add.reduce(x*x,axis) */
 	obj1 = PyArray_GenericReduceFunction((PyAO *)obj2, n_ops.add,
-					     axis, rtype);
+					     axis, rtype, NULL);
 	Py_DECREF(obj2);
 	if (obj1 == NULL) {Py_DECREF(new); return NULL;}
 
@@ -686,6 +692,15 @@ PyArray_Std(PyArrayObject *self, int axis, int rtype, int variance)
 	if (obj1 == NULL) return NULL;
 	ret = PyArray_View((PyAO *)obj1, NULL, self->ob_type);
 	Py_DECREF(obj1);
+        if (out) {
+                if (PyArray_CopyAnyInto(out, (PyArrayObject *)ret) < 0) {
+                        Py_DECREF(ret);
+                        return NULL;
+                }
+                Py_DECREF(ret);
+                Py_INCREF(out);
+                return (PyObject *)out;
+        }
 	return ret;
 }
 
@@ -694,14 +709,14 @@ PyArray_Std(PyArrayObject *self, int axis, int rtype, int variance)
  Sum
 */
 static PyObject *
-PyArray_Sum(PyArrayObject *self, int axis, int rtype)
+PyArray_Sum(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 {
 	PyObject *new, *ret;
 
 	if ((new = _check_axis(self, &axis, 0))==NULL) return NULL;
 
 	ret = PyArray_GenericReduceFunction((PyAO *)new, n_ops.add, axis, 
-					    rtype);
+					    rtype, out);
 	Py_DECREF(new);
 	return ret;
 }
@@ -710,14 +725,14 @@ PyArray_Sum(PyArrayObject *self, int axis, int rtype)
  Prod
 */
 static PyObject *
-PyArray_Prod(PyArrayObject *self, int axis, int rtype)
+PyArray_Prod(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 {
 	PyObject *new, *ret;
 
 	if ((new = _check_axis(self, &axis, 0))==NULL) return NULL;
 
 	ret = PyArray_GenericReduceFunction((PyAO *)new, n_ops.multiply, axis,
-					    rtype);
+					    rtype, out);
 	Py_DECREF(new);
 	return ret;
 }
@@ -726,14 +741,14 @@ PyArray_Prod(PyArrayObject *self, int axis, int rtype)
  CumSum
 */
 static PyObject *
-PyArray_CumSum(PyArrayObject *self, int axis, int rtype)
+PyArray_CumSum(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 {
 	PyObject *new, *ret;
 
 	if ((new = _check_axis(self, &axis, 0))==NULL) return NULL;
 
 	ret = PyArray_GenericAccumulateFunction((PyAO *)new, n_ops.add, axis,
-						rtype);
+						rtype, out);
 	Py_DECREF(new);
 	return ret;
 }
@@ -742,7 +757,7 @@ PyArray_CumSum(PyArrayObject *self, int axis, int rtype)
  CumProd
 */
 static PyObject *
-PyArray_CumProd(PyArrayObject *self, int axis, int rtype)
+PyArray_CumProd(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 {
 	PyObject *new, *ret;
 
@@ -750,7 +765,7 @@ PyArray_CumProd(PyArrayObject *self, int axis, int rtype)
 
 	ret = PyArray_GenericAccumulateFunction((PyAO *)new, 
 						n_ops.multiply, axis,
-						rtype);
+						rtype, out);
 	Py_DECREF(new);
 	return ret;
 }
@@ -759,7 +774,7 @@ PyArray_CumProd(PyArrayObject *self, int axis, int rtype)
  Any
 */
 static PyObject *
-PyArray_Any(PyArrayObject *self, int axis)
+PyArray_Any(PyArrayObject *self, int axis, PyArrayObject *out)
 {
 	PyObject *new, *ret;
 
@@ -767,7 +782,7 @@ PyArray_Any(PyArrayObject *self, int axis)
 
 	ret = PyArray_GenericReduceFunction((PyAO *)new, 
 					    n_ops.logical_or, axis, 
-					    PyArray_BOOL);
+					    PyArray_BOOL, out);
 	Py_DECREF(new);
 	return ret;
 }
@@ -776,7 +791,7 @@ PyArray_Any(PyArrayObject *self, int axis)
  All
 */
 static PyObject *
-PyArray_All(PyArrayObject *self, int axis)
+PyArray_All(PyArrayObject *self, int axis, PyArrayObject *out)
 {
 	PyObject *new, *ret;
 
@@ -784,7 +799,7 @@ PyArray_All(PyArrayObject *self, int axis)
 
 	ret = PyArray_GenericReduceFunction((PyAO *)new, 
 					    n_ops.logical_and, axis, 
-					    PyArray_BOOL);
+					    PyArray_BOOL, out);
 	Py_DECREF(new);
 	return ret;
 }
@@ -794,7 +809,8 @@ PyArray_All(PyArrayObject *self, int axis)
  Compress
 */
 static PyObject *
-PyArray_Compress(PyArrayObject *self, PyObject *condition, int axis)
+PyArray_Compress(PyArrayObject *self, PyObject *condition, int axis, 
+                 PyArrayObject *out)
 {
         PyArrayObject *cond;
 	PyObject *res, *ret;
@@ -908,7 +924,7 @@ PyArray_Clip(PyArrayObject *self, PyObject *min, PyObject *max)
 
 	newtup = Py_BuildValue("(OOO)", (PyObject *)self, min, max);
 	if (newtup == NULL) {Py_DECREF(selector); return NULL;}
-	ret = PyArray_Choose((PyAO *)selector, newtup);
+	ret = PyArray_Choose((PyAO *)selector, newtup, NULL, NPY_RAISE);
 	Py_DECREF(selector);
 	Py_DECREF(newtup);
 	return ret;
@@ -961,13 +977,13 @@ PyArray_Conjugate(PyArrayObject *self)
 */
 static PyObject *
 PyArray_Trace(PyArrayObject *self, int offset, int axis1, int axis2, 
-int rtype)
+              int rtype, PyArrayObject *out)
 {
 	PyObject *diag=NULL, *ret=NULL;
 
 	diag = PyArray_Diagonal(self, offset, axis1, axis2);
 	if (diag == NULL) return NULL;
-	ret = PyArray_GenericReduceFunction((PyAO *)diag, n_ops.add, -1, rtype);
+	ret = PyArray_GenericReduceFunction((PyAO *)diag, n_ops.add, -1, rtype, out);
 	Py_DECREF(diag);
 	return ret;
 }
@@ -1762,19 +1778,20 @@ PyArray_ConvertToCommonType(PyObject *op, int *retn)
 }
 
 
-
 /*MULTIARRAY_API
 */
 static PyObject *
-PyArray_Choose(PyArrayObject *ip, PyObject *op)
+PyArray_Choose(PyArrayObject *ip, PyObject *op, PyArrayObject *ret, 
+               NPY_CLIPMODE clipmode)
 {
 	intp *sizes, offset;
-	int i,n,m,elsize;
+	int n, elsize;
+        intp i, m;
 	char *ret_data;
-	PyArrayObject **mps, *ap, *ret;
+	PyArrayObject **mps, *ap;
 	intp *self_data, mi;
+        int copyret=0;
 	ap = NULL;
-	ret = NULL;
 	       	
 	/* Convert all inputs to arrays of a common type */
 	mps = PyArray_ConvertToCommonType(op, &n);
@@ -1784,8 +1801,8 @@ PyArray_Choose(PyArrayObject *ip, PyObject *op)
 	if (sizes == NULL) goto fail;
 	
 	ap = (PyArrayObject *)PyArray_ContiguousFromAny((PyObject *)ip, 
-							   PyArray_INTP, 
-							   0, 0);
+                                                        PyArray_INTP, 
+                                                        0, 0);
 	if (ap == NULL) goto fail;
 	
 	/* Check the dimensions of the arrays */
@@ -1805,29 +1822,70 @@ PyArray_Choose(PyArrayObject *ip, PyObject *op)
 		sizes[i] = PyArray_NBYTES(mps[i]);
 	}
 	
-	Py_INCREF(mps[0]->descr);
-	ret = (PyArrayObject *)PyArray_NewFromDescr(ap->ob_type, 
-						    mps[0]->descr,
-						    ap->nd,
-						    ap->dimensions, 
-						    NULL, NULL, 0,
-						    (PyObject *)ap);
-	if (ret == NULL) goto fail;
-	
+        Py_INCREF(mps[0]->descr);
+        if (!ret) {
+                ret = (PyArrayObject *)PyArray_NewFromDescr(ap->ob_type, 
+                                                            mps[0]->descr,
+                                                            ap->nd,
+                                                            ap->dimensions, 
+                                                            NULL, NULL, 0,
+                                                            (PyObject *)ap);
+                if (ret == NULL) goto fail;
+        }
+        else {
+                PyArrayObject *obj;
+                int flags = NPY_CARRAY | NPY_UPDATEIFCOPY;
+
+                if (!PyArray_SAMESHAPE(ap, ret)) {
+                        PyErr_SetString(PyExc_TypeError, 
+                                        "invalid shape for output array.");
+                        ret = NULL;
+                        Py_DECREF(mps[0]->descr);
+                        goto fail;
+                }
+                if (clipmode == NPY_RAISE) {
+                        /* we need to make sure and get a copy
+                           so the input array is not changed
+                           before the error is called 
+                        */
+                        flags |= NPY_ENSURECOPY;
+                }
+                obj = (PyArrayObject *)PyArray_FromArray(ret, mps[0]->descr, 
+                                                         flags);
+                if (obj != ret) copyret = 1;
+                ret = obj;
+        }
+                
 	elsize = ret->descr->elsize;
 	m = PyArray_SIZE(ret);
 	self_data = (intp *)ap->data;
-	ret_data = ret->data;
-	
+	ret_data = ret->data;           
+
 	for (i=0; i<m; i++) {
 		mi = *self_data;
 		if (mi < 0 || mi >= n) {
-			PyErr_SetString(PyExc_ValueError, 
-					"invalid entry in choice array");
-			goto fail;
-		}
-		offset = i*elsize;
-		if (offset >= sizes[mi]) {offset = offset % sizes[mi]; }
+                        switch(clipmode) {
+                        case NPY_RAISE:
+                                PyErr_SetString(PyExc_ValueError, 
+                                                "invalid entry in choice "\
+                                                "array");
+                                goto fail;
+                        case NPY_WRAP:
+                                if (mi < 0) {
+                                        while(mi<0) mi += n;
+                                }
+                                else {
+                                        while(mi>=n) mi -= n;
+                                }
+                                break;
+                        case NPY_CLIP:
+                                if (mi < 0) mi=0;
+                                else if (mi>=n) mi=n-1;
+                                break;
+                        }
+                }
+                offset = i*elsize;
+                if (offset >= sizes[mi]) {offset = offset % sizes[mi]; }
 		memmove(ret_data, mps[mi]->data+offset, elsize);
 		ret_data += elsize; self_data++;
 	}
@@ -1837,7 +1895,13 @@ PyArray_Choose(PyArrayObject *ip, PyObject *op)
 	Py_DECREF(ap);
 	PyDataMem_FREE(mps);
 	_pya_free(sizes);
-
+        if (copyret) {
+                PyObject *obj;
+                obj = ret->base;
+                Py_INCREF(obj);
+                Py_DECREF(ret); 
+                ret = (PyArrayObject *)obj;
+        }
 	return (PyObject *)ret;
 	
  fail:
@@ -1845,7 +1909,7 @@ PyArray_Choose(PyArrayObject *ip, PyObject *op)
 	Py_XDECREF(ap);
 	PyDataMem_FREE(mps);
 	_pya_free(sizes);
-	Py_XDECREF(ret);
+        PyArray_XDECREF_ERR(ret);
 	return NULL;
 }
 
@@ -2984,7 +3048,7 @@ PyArray_ArgMin(PyArrayObject *ap, int axis)
  Max
 */
 static PyObject *
-PyArray_Max(PyArrayObject *ap, int axis)
+PyArray_Max(PyArrayObject *ap, int axis, PyArrayObject *out)
 {
 	PyArrayObject *arr;
 	PyObject *ret;
@@ -2992,7 +3056,7 @@ PyArray_Max(PyArrayObject *ap, int axis)
 	if ((arr=(PyArrayObject *)_check_axis(ap, &axis, 0))==NULL)
 		return NULL;
 	ret = PyArray_GenericReduceFunction(arr, n_ops.maximum, axis, 
-					    arr->descr->type_num);
+					    arr->descr->type_num, out);
 	Py_DECREF(arr);
 	return ret;	    
 }
@@ -3001,7 +3065,7 @@ PyArray_Max(PyArrayObject *ap, int axis)
  Min
 */
 static PyObject *
-PyArray_Min(PyArrayObject *ap, int axis)
+PyArray_Min(PyArrayObject *ap, int axis, PyArrayObject *out)
 {
 	PyArrayObject *arr;
 	PyObject *ret;
@@ -3009,7 +3073,7 @@ PyArray_Min(PyArrayObject *ap, int axis)
 	if ((arr=(PyArrayObject *)_check_axis(ap, &axis, 0))==NULL)
 		return NULL;
 	ret = PyArray_GenericReduceFunction(arr, n_ops.minimum, axis,
-					    arr->descr->type_num);
+					    arr->descr->type_num, out);
 	Py_DECREF(arr);
 	return ret;	    
 }
@@ -3018,7 +3082,7 @@ PyArray_Min(PyArrayObject *ap, int axis)
  Ptp
 */
 static PyObject *
-PyArray_Ptp(PyArrayObject *ap, int axis)
+PyArray_Ptp(PyArrayObject *ap, int axis, PyArrayObject *out)
 {
 	PyArrayObject *arr;
 	PyObject *ret;
@@ -3026,12 +3090,17 @@ PyArray_Ptp(PyArrayObject *ap, int axis)
 
 	if ((arr=(PyArrayObject *)_check_axis(ap, &axis, 0))==NULL)
 		return NULL;
-	obj1 = PyArray_Max(arr, axis);
+	obj1 = PyArray_Max(arr, axis, out);
 	if (obj1 == NULL) goto fail;
-	obj2 = PyArray_Min(arr, axis);
+	obj2 = PyArray_Min(arr, axis, NULL);
 	if (obj2 == NULL) goto fail;
 	Py_DECREF(arr);
-	ret = PyNumber_Subtract(obj1, obj2);
+        if (out) {
+                ret = PyObject_CallFunction(n_ops.subtract, "OOO", out, obj2, out);
+        }
+        else {
+                ret = PyNumber_Subtract(obj1, obj2);
+        }
 	Py_DECREF(obj1);
 	Py_DECREF(obj2);
 	return ret;
@@ -3132,12 +3201,14 @@ PyArray_ArgMax(PyArrayObject *op, int axis)
  Take
 */
 static PyObject *
-PyArray_Take(PyArrayObject *self0, PyObject *indices0, int axis)
+PyArray_TakeOut(PyArrayObject *self0, PyObject *indices0, int axis, 
+                PyArrayObject *ret, NPY_CLIPMODE clipmode)
 {
-        PyArrayObject *self, *indices, *ret;
+        PyArrayObject *self, *indices;
         intp nd, i, j, n, m, max_item, tmp, chunk;
 	intp shape[MAX_DIMS];
         char *src, *dest;
+        int copyret=0;
 	
         indices = ret = NULL;
 	self = (PyAO *)_check_axis(self0, &axis, CARRAY);
@@ -3165,45 +3236,111 @@ PyArray_Take(PyArrayObject *self0, PyObject *indices0, int axis)
                 }
         }
 	Py_INCREF(self->descr);
-        ret = (PyArrayObject *)PyArray_NewFromDescr(self->ob_type, 
-						    self->descr,
-						    nd, shape, 
-						    NULL, NULL, 0, 
-						    (PyObject *)self);
-	
-        if (ret == NULL) goto fail;
+        if (!ret) {
+                ret = (PyArrayObject *)PyArray_NewFromDescr(self->ob_type, 
+                                                            self->descr,
+                                                            nd, shape, 
+                                                            NULL, NULL, 0, 
+                                                            (PyObject *)self);
+                
+                if (ret == NULL) goto fail;
+        }
+        else {
+                PyArrayObject *obj;
+                int flags = NPY_CARRAY | NPY_UPDATEIFCOPY;
+
+                if ((ret->nd != nd) || 
+                    !PyArray_CompareLists(ret->dimensions, shape, nd)) {
+                        PyErr_SetString(PyExc_ValueError, 
+                                        "bad shape in output array");
+                        ret = NULL;
+                        Py_DECREF(self->descr);
+                        goto fail;
+                }
+
+                if (clipmode == NPY_RAISE) {
+                        /* we need to make sure and get a copy
+                           so the input array is not changed
+                           before the error is called 
+                        */
+                        flags |= NPY_ENSURECOPY;
+                }
+                obj = (PyArrayObject *)PyArray_FromArray(ret, self->descr,
+                                                         flags);
+                if (obj != ret) copyret = 1;
+                ret = obj;                
+        }
 	
         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 = ((intp *)(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;
+
+        switch(clipmode) {
+        case NPY_RAISE:
+                for(i=0; i<n; i++) {
+                        for(j=0; j<m; j++) {
+                                tmp = ((intp *)(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;
+                                }
+                                memmove(dest, src+tmp*chunk, chunk);
+                                dest += chunk;
+                        }
+                        src += chunk*max_item;
+                }
+                break;
+        case NPY_WRAP:
+                for(i=0; i<n; i++) {
+                        for(j=0; j<m; j++) {
+                                tmp = ((intp *)(indices->data))[j];
+                                if (tmp < 0) while (tmp < 0) tmp += max_item;
+                                else if (tmp >= max_item)
+                                        while (tmp >= max_item) 
+                                                tmp -= max_item;
+                                memmove(dest, src+tmp*chunk, chunk);
+                                dest += chunk;
+                        }
+                        src += chunk*max_item;
+                }
+                break;
+        case NPY_CLIP:
+                for(i=0; i<n; i++) {
+                        for(j=0; j<m; j++) {
+                                tmp = ((intp *)(indices->data))[j];
+                                if (tmp < 0) 
+                                        tmp = 0;
+                                else if (tmp >= max_item) 
+                                        tmp = max_item-1;
+                                memmove(dest, src+tmp*chunk, chunk);
+                                dest += chunk;
                         }
-                        memmove(dest, src+tmp*chunk, chunk);
-                        dest += chunk;
+                        src += chunk*max_item;
                 }
-                src += chunk*max_item;
+                break;
         }
-	
+
         PyArray_INCREF(ret);
 
         Py_XDECREF(indices);
         Py_XDECREF(self);
+        if (copyret) {
+                PyObject *obj;
+                obj = ret->base;
+                Py_INCREF(obj);
+                Py_DECREF(ret); 
+                ret = (PyArrayObject *)obj;
+        }
 
         return (PyObject *)ret;
 	
 	
  fail:
-        Py_XDECREF(ret);
+        PyArray_XDECREF_ERR(ret);
         Py_XDECREF(indices);
         Py_XDECREF(self);
         return NULL;
@@ -3213,28 +3350,40 @@ PyArray_Take(PyArrayObject *self0, PyObject *indices0, int axis)
  Put values into an array
 */
 static PyObject *
-PyArray_Put(PyArrayObject *self, PyObject* values0, PyObject *indices0) 
+PyArray_PutIn(PyArrayObject *self, PyObject* values0, PyObject *indices0,
+              NPY_CLIPMODE clipmode)
 {
         PyArrayObject  *indices, *values;
         int i, chunk, ni, max_item, nv, tmp, thistype; 
         char *src, *dest;
+        int copied = 0;
 
         indices = NULL;
         values = NULL;
 
         if (!PyArray_Check(self)) {
-                PyErr_SetString(PyExc_TypeError, "put: first argument must be an array");
+                PyErr_SetString(PyExc_TypeError, 
+                                "put: first argument must be an array");
                 return NULL;
         }
         if (!PyArray_ISCONTIGUOUS(self)) {
-                PyErr_SetString(PyExc_ValueError, "put: first argument must be contiguous");
-                return NULL;
+                PyArrayObject *obj;
+                int flags = NPY_CARRAY | NPY_UPDATEIFCOPY;
+                if (clipmode == NPY_RAISE) {
+                        flags |= NPY_ENSURECOPY;
+                }
+                Py_INCREF(self->descr);
+                obj = (PyArrayObject *)PyArray_FromArray(self, 
+                                                         self->descr, flags);
+                if (obj != self) copied = 1;
+                self = obj;
         }
         max_item = PyArray_SIZE(self);
         dest = self->data;
         chunk = self->descr->elsize;
 
-        indices = (PyArrayObject *)PyArray_ContiguousFromAny(indices0, PyArray_INTP, 0, 0);
+        indices = (PyArrayObject *)PyArray_ContiguousFromAny(indices0, 
+                                                             PyArray_INTP, 0, 0);
         if (indices == NULL) goto fail;
         ni = PyArray_SIZE(indices);
 
@@ -3244,44 +3393,107 @@ PyArray_Put(PyArrayObject *self, PyObject* values0, PyObject *indices0)
 						  DEFAULT | FORCECAST, NULL); 
         if (values == NULL) goto fail;
         nv = PyArray_SIZE(values);
-        if (nv > 0) { /* nv == 0 for a null array */
-                if (self->descr->hasobject) {   
+        if (nv <= 0) goto finish;
+        if (self->descr->hasobject) {   
+                switch(clipmode) {
+                case NPY_RAISE:
                         for(i=0; i<ni; i++) {
                                 src = values->data + chunk * (i % nv);
                                 tmp = ((intp *)(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");
+                                        PyErr_SetString(PyExc_IndexError, 
+                                                        "index out of " \
+                                                        "range for array");
                                         goto fail;
                                 }
-				PyArray_Item_INCREF(src, self->descr);
-				PyArray_Item_XDECREF(dest+tmp*chunk, self->descr);
+                                PyArray_Item_INCREF(src, self->descr);
+                                PyArray_Item_XDECREF(dest+tmp*chunk, self->descr);
                                 memmove(dest + tmp * chunk, src, chunk);
                         }
+                        break;
+                case NPY_WRAP:
+                        for(i=0; i<ni; i++) {
+                                src = values->data + chunk * (i % nv);
+                                tmp = ((intp *)(indices->data))[i];
+                                if (tmp < 0) while(tmp < 0) tmp+=max_item;
+                                else if (tmp >= max_item) 
+                                        while(tmp >= max_item)
+                                                tmp -= max_item;
+                                PyArray_Item_INCREF(src, self->descr);
+                                PyArray_Item_XDECREF(dest+tmp*chunk, self->descr);
+                                memmove(dest + tmp * chunk, src, chunk);
+                        }
+                        break;
+                case NPY_CLIP:
+                        for(i=0; i<ni; i++) {
+                                src = values->data + chunk * (i % nv);
+                                tmp = ((intp *)(indices->data))[i];
+                                if (tmp < 0) tmp = 0;
+                                else if (tmp >= max_item)
+                                        tmp = max_item - 1;
+                                PyArray_Item_INCREF(src, self->descr);
+                                PyArray_Item_XDECREF(dest+tmp*chunk, self->descr);
+                                memmove(dest + tmp * chunk, src, chunk);
+                        }
+                        break;
                 }
-                else {
+        }
+        else {
+                switch(clipmode) {
+                case NPY_RAISE:
                         for(i=0; i<ni; i++) {
                                 src = values->data + chunk * (i % nv);
                                 tmp = ((intp *)(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");
+                                        PyErr_SetString(PyExc_IndexError, 
+                                                        "index out of " \
+                                                        "range for array");
                                         goto fail;
                                 }
                                 memmove(dest + tmp * chunk, src, chunk);
                         }
+                        break;
+                case NPY_WRAP:
+                        for(i=0; i<ni; i++) {
+                                src = values->data + chunk * (i % nv);
+                                tmp = ((intp *)(indices->data))[i];
+                                if (tmp < 0) while(tmp < 0) tmp+=max_item;
+                                else if (tmp >= max_item) 
+                                        while(tmp >= max_item)
+                                                tmp -= max_item;
+                                memmove(dest + tmp * chunk, src, chunk);
+                        }
+                        break;
+                case NPY_CLIP:
+                        for(i=0; i<ni; i++) {
+                                src = values->data + chunk * (i % nv);
+                                tmp = ((intp *)(indices->data))[i];
+                                if (tmp < 0) tmp = 0;
+                                else if (tmp >= max_item)
+                                        tmp = max_item - 1;
+                        memmove(dest + tmp * chunk, src, chunk);
+                        }
+                        break;
                 }
-
         }
 
+ finish: 
         Py_XDECREF(values);
         Py_XDECREF(indices);
+        if (copied) {
+                Py_DECREF(self);
+        }
         Py_INCREF(Py_None);
         return Py_None;
 	
  fail:
         Py_XDECREF(indices);
         Py_XDECREF(values);
+        if (copied) {
+                PyArray_XDECREF_ERR(self);
+        }
         return NULL;
 }
 
@@ -3289,11 +3501,12 @@ PyArray_Put(PyArrayObject *self, PyObject* values0, PyObject *indices0)
  Put values into an array according to a mask.
 */
 static PyObject *
-PyArray_PutMask(PyArrayObject *self, PyObject* values0, PyObject* mask0) 
+PyArray_PutMask(PyArrayObject *self, PyObject* values0, PyObject* mask0)
 {
         PyArrayObject  *mask, *values;
         int i, chunk, ni, max_item, nv, tmp, thistype;
         char *src, *dest;
+        int copied=0;
 
         mask = NULL;
         values = NULL;
@@ -3305,9 +3518,13 @@ PyArray_PutMask(PyArrayObject *self, PyObject* values0, PyObject* mask0)
                 return NULL;
         }
         if (!PyArray_ISCONTIGUOUS(self)) {
-                PyErr_SetString(PyExc_ValueError, 
-				"putmask: first argument must be contiguous");
-                return NULL;
+                PyArrayObject *obj;
+                int flags = NPY_CARRAY | NPY_UPDATEIFCOPY;
+                Py_INCREF(self->descr);
+                obj = (PyArrayObject *)PyArray_FromArray(self, 
+                                                         self->descr, flags);
+                if (obj != self) copied = 1;
+                self = obj;
         }
 
         max_item = PyArray_SIZE(self);
@@ -3330,6 +3547,12 @@ PyArray_PutMask(PyArrayObject *self, PyObject* values0, PyObject* mask0)
 		PyArray_ContiguousFromAny(values0, thistype, 0, 0);
 	if (values == NULL) goto fail;
         nv = PyArray_SIZE(values);	 /* zero if null array */
+        if (nv <= 0) {
+                Py_XDECREF(values);
+                Py_XDECREF(mask);
+                Py_INCREF(Py_None);
+                return Py_None;
+        }
         if (nv > 0) {
                 if (self->descr->hasobject) {
                         for(i=0; i<ni; i++) {
@@ -3353,12 +3576,18 @@ PyArray_PutMask(PyArrayObject *self, PyObject* values0, PyObject* mask0)
 
         Py_XDECREF(values);
         Py_XDECREF(mask);
+        if (copied) {
+                Py_DECREF(self);
+        }
         Py_INCREF(Py_None);
         return Py_None;
 	
  fail:
         Py_XDECREF(mask);
         Py_XDECREF(values);
+        if (copied) {
+                PyArray_XDECREF_ERR(self);
+        }
         return NULL;
 }
 
@@ -3391,6 +3620,30 @@ PyArray_Converter(PyObject *object, PyObject **address)
 }
 
 /*MULTIARRAY_API
+ Useful to pass as converter function for O& processing in
+ PyArgs_ParseTuple for output arrays
+*/
+static int 
+PyArray_OutputConverter(PyObject *object, PyArrayObject **address) 
+{
+        if (object == NULL || object == Py_None) {
+                *address = NULL;
+                return PY_SUCCEED;
+        }
+        if (PyArray_Check(object)) {
+                *address = (PyArrayObject *)object;
+                return PY_SUCCEED;
+        }
+        else {
+                PyErr_SetString(PyExc_TypeError, 
+                                "output must be an array");
+                *address = NULL;
+                return PY_FAIL;
+        }
+}
+
+
+/*MULTIARRAY_API
  Convert an object to true / false
 */
 static int
@@ -3428,17 +3681,67 @@ PyArray_OrderConverter(PyObject *object, NPY_ORDER *val)
                 if (str[0] == 'C' || str[0] == 'c') {
                         *val = PyArray_CORDER; 
                 }
-                if (str[0] == 'F' || str[0] == 'f') {
+                else if (str[0] == 'F' || str[0] == 'f') {
                         *val = PyArray_FORTRANORDER;
                 }
-                if (str[0] == 'A' || str[0] == 'a') {
+                else if (str[0] == 'A' || str[0] == 'a') {
                         *val = PyArray_ANYORDER;
                 }
+                else {
+                        PyErr_SetString(PyExc_TypeError, 
+                                        "order not understood");
+                        return PY_FAIL;
+                }
+        }
+        return PY_SUCCEED;
+}
+
+/*MULTIARRAY_API
+ Convert an object to NPY_RAISE / NPY_CLIP / NPY_WRAP
+*/
+static int
+PyArray_ClipmodeConverter(PyObject *object, NPY_CLIPMODE *val)
+{
+        if (object == NULL || object == Py_None) {
+                *val = NPY_RAISE;
+        }
+        else if (PyString_Check(object)) {
+                char *str;
+                str = PyString_AS_STRING(object);
+                if (str[0] == 'C' || str[0] == 'c') {
+                        *val = NPY_CLIP; 
+                }
+                else if (str[0] == 'W' || str[0] == 'w') {
+                        *val = NPY_WRAP;
+                }
+                else if (str[0] == 'R' || str[0] == 'r') {
+                        *val = NPY_RAISE;
+                }
+                else {
+                        PyErr_SetString(PyExc_TypeError,
+                                        "clipmode not understood");
+                        return PY_FAIL;
+                }
         }
+        else {
+                int number;
+                number = PyInt_AsLong(object);
+                if (number == -1 && PyErr_Occurred()) goto fail;
+                if (number <= (int) NPY_RAISE &&
+                    number >= (int) NPY_CLIP)
+                        *val = (NPY_CLIPMODE) number;
+                else goto fail;
+        }        
         return PY_SUCCEED;
+
+ fail:
+        PyErr_SetString(PyExc_TypeError,
+                        "clipmode not understood");
+        return PY_FAIL;
 }
 
 
+
 /*MULTIARRAY_API
  Typestr converter
 */
@@ -5358,7 +5661,7 @@ array_fromfile(PyObject *ignored, PyObject *args, PyObject *keywds)
 
 	if (type == NULL) type = PyArray_DescrFromType(PyArray_DEFAULT);
 
-	if (PyString_Check(file)) {
+	if (PyString_Check(file) || PyUnicode_Check(file)) {
 		file = PyObject_CallFunction((PyObject *)&PyFile_Type,
 					     "Os", file, "rb");
 		if (file==NULL) return NULL;
@@ -5901,7 +6204,7 @@ PyArray_Where(PyObject *condition, PyObject *x, PyObject *y)
 	tup = Py_BuildValue("(OO)", y, x);
 	if (tup == NULL) {Py_DECREF(obj); return NULL;}
 
-	ret = PyArray_Choose((PyAO *)obj, tup);
+	ret = PyArray_Choose((PyAO *)obj, tup, NULL, NPY_RAISE);
 
 	Py_DECREF(obj);
 	Py_DECREF(tup);
@@ -6287,6 +6590,18 @@ PyMODINIT_FUNC initmultiarray(void) {
 	PyDict_SetItemString(d, "BUFSIZE", s);
 	Py_DECREF(s);
 
+	s = PyInt_FromLong(NPY_CLIP);
+	PyDict_SetItemString(d, "CLIP", s);
+	Py_DECREF(s);
+
+	s = PyInt_FromLong(NPY_RAISE);
+	PyDict_SetItemString(d, "RAISE", s);
+	Py_DECREF(s);
+
+	s = PyInt_FromLong(NPY_WRAP);
+	PyDict_SetItemString(d, "WRAP", s);
+	Py_DECREF(s);
+
 	s = PyInt_FromLong(NPY_MAXDIMS);
 	PyDict_SetItemString(d, "MAXDIMS", s);
 	Py_DECREF(s);
diff --git a/numpy/core/src/scalartypes.inc.src b/numpy/core/src/scalartypes.inc.src
index efd0816e4..1f92ee24e 100644
--- a/numpy/core/src/scalartypes.inc.src
+++ b/numpy/core/src/scalartypes.inc.src
@@ -1071,7 +1071,7 @@ gentype_wraparray(PyObject *scalar, PyObject *args)
 
 /**begin repeat
 
-#name=tolist, item, tostring, astype, copy, __deepcopy__, choose, searchsorted, view, swapaxes, conj, conjugate, nonzero, flatten, ravel, fill, transpose, newbyteorder#
+#name=tolist, item, tostring, astype, copy, __deepcopy__, searchsorted, view, swapaxes, conj, conjugate, nonzero, flatten, ravel, fill, transpose, newbyteorder#
 */
 
 static PyObject *
@@ -1130,7 +1130,7 @@ gentype_byteswap(PyObject *self, PyObject *args)
 
 /**begin repeat
 
-#name=take, getfield, put, putmask, repeat, tofile, mean, trace, diagonal, clip, std, var, sum, cumsum, prod, cumprod, compress, sort, argsort, round, argmax, argmin, max, min, ptp, any, all, resize, reshape#
+#name=take, getfield, put, putmask, repeat, tofile, mean, trace, diagonal, clip, std, var, sum, cumsum, prod, cumprod, compress, sort, argsort, round, argmax, argmin, max, min, ptp, any, all, resize, reshape, choose#
 */
 
 static PyObject *
@@ -1344,7 +1344,7 @@ static PyMethodDef gentype_methods[] = {
         {"repeat",      (PyCFunction)gentype_repeat,
          METH_VARARGS|METH_KEYWORDS, NULL},
         {"choose",      (PyCFunction)gentype_choose,
-         METH_VARARGS, NULL},
+         METH_VARARGS|METH_KEYWORDS, NULL},
         {"sort",        (PyCFunction)gentype_sort,
          METH_VARARGS, NULL},
         {"argsort",     (PyCFunction)gentype_argsort,
diff --git a/numpy/core/src/ufuncobject.c b/numpy/core/src/ufuncobject.c
index 573bff639..a2b711093 100644
--- a/numpy/core/src/ufuncobject.c
+++ b/numpy/core/src/ufuncobject.c
@@ -1755,25 +1755,27 @@ _create_reduce_copy(PyUFuncReduceObject *loop, PyArrayObject **arr, int rtype)
 }
 
 static PyUFuncReduceObject *
-construct_reduce(PyUFuncObject *self, PyArrayObject **arr, int axis, 
-		 int otype, int operation, intp ind_size, char *str)
+construct_reduce(PyUFuncObject *self, PyArrayObject **arr, PyArrayObject *out, 
+                 int axis, int otype, int operation, intp ind_size, char *str)
 {
         PyUFuncReduceObject *loop;
         PyArrayObject *idarr;
 	PyArrayObject *aar;
-        intp loop_i[MAX_DIMS];
+        intp loop_i[MAX_DIMS], outsize;
         int arg_types[3] = {otype, otype, otype};
 	PyArray_SCALARKIND scalars[3] = {PyArray_NOSCALAR, PyArray_NOSCALAR, 
 					 PyArray_NOSCALAR};
 	int i, j;
 	int nd = (*arr)->nd;
+        int flags;
 	/* Reduce type is the type requested of the input 
 	   during reduction */
         
         if ((loop = _pya_malloc(sizeof(PyUFuncReduceObject)))==NULL) {
                 PyErr_NoMemory(); return loop;
         }
-        
+
+        loop->retbase=0;
         loop->swap = 0;
 	loop->index = 0;
 	loop->ufunc = self;
@@ -1853,6 +1855,7 @@ construct_reduce(PyUFuncObject *self, PyArrayObject **arr, int axis,
         }
 	
         /* Construct return array */
+        flags = NPY_CARRAY | NPY_UPDATEIFCOPY | NPY_FORCECAST;
 	switch(operation) {
 	case UFUNC_REDUCE:
 		for (j=0, i=0; i<nd; i++) {
@@ -1860,23 +1863,38 @@ construct_reduce(PyUFuncObject *self, PyArrayObject **arr, int axis,
 				loop_i[j++] = (aar)->dimensions[i];
 			
 		}
-		loop->ret = (PyArrayObject *)				\
-			PyArray_New(aar->ob_type, aar->nd-1, loop_i, otype, 
-				    NULL, NULL, 0, 0, (PyObject *)aar);
+                if (out == NULL) {
+                        loop->ret = (PyArrayObject *)                   \
+                                PyArray_New(aar->ob_type, aar->nd-1, loop_i, 
+                                            otype, NULL, NULL, 0, 0, 
+                                            (PyObject *)aar);
+                }
+                else {
+                        outsize == PyArray_MultiplyList(loop_i, aar->nd-1);
+                }
 		break;
 	case UFUNC_ACCUMULATE:
-		loop->ret = (PyArrayObject *)				\
-			PyArray_New(aar->ob_type, aar->nd, aar->dimensions, 
-				    otype, NULL, NULL, 0, 0, (PyObject *)aar);
+                if (out == NULL) {
+                        loop->ret = (PyArrayObject *)                   \
+                                PyArray_New(aar->ob_type, aar->nd, aar->dimensions, 
+                                            otype, NULL, NULL, 0, 0, (PyObject *)aar);
+                }
+                else {
+                        outsize = PyArray_MultiplyList(aar->dimensions, aar->nd);
+                }
 		break;
 	case UFUNC_REDUCEAT:
 		memcpy(loop_i, aar->dimensions, nd*sizeof(intp));
 		/* Index is 1-d array */
 		loop_i[axis] = ind_size; 
-		loop->ret = (PyArrayObject *)\
-			PyArray_New(aar->ob_type, aar->nd, loop_i, otype,
-				    NULL, NULL, 0, 0, (PyObject *)aar);
-		if (loop->ret == NULL) goto fail;
+                if (out == NULL) {
+                        loop->ret = (PyArrayObject *)                   \
+                                PyArray_New(aar->ob_type, aar->nd, loop_i, otype,
+                                            NULL, NULL, 0, 0, (PyObject *)aar);
+                }
+                else {
+                        outsize = PyArray_MultiplyList(loop_i, aar->nd);
+                }
 		if (ind_size == 0) {
 			loop->meth = ZERODIM_REDUCELOOP;
 			return loop;
@@ -1885,6 +1903,19 @@ construct_reduce(PyUFuncObject *self, PyArrayObject **arr, int axis,
 			loop->meth = NOBUFFER_REDUCELOOP;
 		break;
 	}
+        if (out) {
+                if (PyArray_SIZE(out) != outsize) {
+                        PyErr_SetString(PyExc_ValueError, 
+                                        "wrong shape for output");
+                        goto fail;
+                }
+                loop->ret = (PyArrayObject *)                   \
+                        PyArray_FromArray(out, PyArray_DescrFromType(otype), 
+                                          flags);
+                if (loop->ret && loop->ret != out) {
+                        loop->retbase = 1;
+                }
+        }
         if (loop->ret == NULL) goto fail;
         loop->insize = aar->descr->elsize;
         loop->outsize = loop->ret->descr->elsize;
@@ -1982,7 +2013,8 @@ construct_reduce(PyUFuncObject *self, PyArrayObject **arr, int axis,
 */
 
 static PyObject *
-PyUFunc_Reduce(PyUFuncObject *self, PyArrayObject *arr, int axis, int otype)
+PyUFunc_Reduce(PyUFuncObject *self, PyArrayObject *arr, PyArrayObject *out, 
+               int axis, int otype)
 {
         PyArrayObject *ret=NULL;
         PyUFuncReduceObject *loop;
@@ -1991,7 +2023,7 @@ PyUFunc_Reduce(PyUFuncObject *self, PyArrayObject *arr, int axis, int otype)
         NPY_BEGIN_THREADS_DEF
         	
         /* Construct loop object */
-        loop = construct_reduce(self, &arr, axis, otype, UFUNC_REDUCE, 0,
+        loop = construct_reduce(self, &arr, out, axis, otype, UFUNC_REDUCE, 0,
 				"reduce");
 	if (!loop) return NULL;
 
@@ -2108,8 +2140,9 @@ PyUFunc_Reduce(PyUFuncObject *self, PyArrayObject *arr, int axis, int otype)
 
         NPY_LOOP_END_THREADS
 
-        ret = loop->ret;
 	/* Hang on to this reference -- will be decref'd with loop */
+        if (loop->retbase) ret = (PyArrayObject *)loop->ret->base;
+        else ret = loop->ret;
         Py_INCREF(ret);
         ufuncreduce_dealloc(loop);
         return (PyObject *)ret;
@@ -2123,8 +2156,8 @@ PyUFunc_Reduce(PyUFuncObject *self, PyArrayObject *arr, int axis, int otype)
 
 
 static PyObject *
-PyUFunc_Accumulate(PyUFuncObject *self, PyArrayObject *arr, int axis, 
-		   int otype)
+PyUFunc_Accumulate(PyUFuncObject *self, PyArrayObject *arr, PyArrayObject *out, 
+                   int axis, int otype)
 {
         PyArrayObject *ret=NULL;
         PyUFuncReduceObject *loop;
@@ -2133,7 +2166,7 @@ PyUFunc_Accumulate(PyUFuncObject *self, PyArrayObject *arr, int axis,
 	NPY_BEGIN_THREADS_DEF
         
         /* Construct loop object */
-        loop = construct_reduce(self, &arr, axis, otype, UFUNC_ACCUMULATE, 0,
+        loop = construct_reduce(self, &arr, out, axis, otype, UFUNC_ACCUMULATE, 0,
 				"accumulate");
 	if (!loop) return NULL;
 
@@ -2253,8 +2286,10 @@ PyUFunc_Accumulate(PyUFuncObject *self, PyArrayObject *arr, int axis,
         }
 
 	NPY_LOOP_END_THREADS
-        ret = loop->ret;
+
 	/* Hang on to this reference -- will be decref'd with loop */
+        if (loop->retbase) ret = (PyArrayObject *)loop->ret->base;
+        else ret = loop->ret;
         Py_INCREF(ret);
         ufuncreduce_dealloc(loop);
         return (PyObject *)ret;
@@ -2287,7 +2322,7 @@ output shape is based on the size of indices
 
 static PyObject *
 PyUFunc_Reduceat(PyUFuncObject *self, PyArrayObject *arr, PyArrayObject *ind, 
-                 int axis, int otype)
+                 PyArrayObject *out, int axis, int otype)
 {	
 	PyArrayObject *ret;
         PyUFuncReduceObject *loop;
@@ -2310,7 +2345,7 @@ PyUFunc_Reduceat(PyUFuncObject *self, PyArrayObject *arr, PyArrayObject *ind,
 	
 	ptr = (intp *)ind->data;
         /* Construct loop object */
-        loop = construct_reduce(self, &arr, axis, otype, UFUNC_REDUCEAT, nn,
+        loop = construct_reduce(self, &arr, out, axis, otype, UFUNC_REDUCEAT, nn,
 				"reduceat");
 	if (!loop) return NULL;
 
@@ -2404,9 +2439,10 @@ PyUFunc_Reduceat(PyUFuncObject *self, PyArrayObject *arr, PyArrayObject *ind,
 	}
 
 	NPY_LOOP_END_THREADS
-	
-        ret = loop->ret;
-	/* Hang on to this reference -- will be decref'd with loop */
+
+        /* Hang on to this reference -- will be decref'd with loop */	
+        if (loop->retbase) ret = (PyArrayObject *)loop->ret->base;
+        else ret = loop->ret;
         Py_INCREF(ret);
         ufuncreduce_dealloc(loop);
         return (PyObject *)ret;
@@ -2434,8 +2470,9 @@ PyUFunc_GenericReduction(PyUFuncObject *self, PyObject *args,
 	PyObject *obj_ind, *context; 
 	PyArrayObject *indices = NULL;
 	PyArray_Descr *otype=NULL;
-	static char *kwlist1[] = {"array", "axis", "dtype", NULL};
-	static char *kwlist2[] = {"array", "indices", "axis", "dtype", NULL}; 
+        PyArrayObject *out=NULL;
+	static char *kwlist1[] = {"array", "axis", "dtype", "out", NULL};
+	static char *kwlist2[] = {"array", "indices", "axis", "dtype", "out", NULL}; 
         static char *_reduce_type[] = {"reduce", "accumulate", \
 				       "reduceat", NULL};
 	if (self == NULL) {
@@ -2460,19 +2497,23 @@ PyUFunc_GenericReduction(PyUFuncObject *self, PyObject *args,
 	if (operation == UFUNC_REDUCEAT) {
 		PyArray_Descr *indtype;
 		indtype = PyArray_DescrFromType(PyArray_INTP);
-		if(!PyArg_ParseTupleAndKeywords(args, kwds, "OO|iO&", kwlist2, 
+		if(!PyArg_ParseTupleAndKeywords(args, kwds, "OO|iO&O&", kwlist2, 
 						&op, &obj_ind, &axis, 
 						PyArray_DescrConverter2, 
-						&otype)) return NULL;
+						&otype,
+                                                PyArray_OutputConverter,
+                                                &out)) return NULL;
                 indices = (PyArrayObject *)PyArray_FromAny(obj_ind, indtype, 
 							   1, 1, CARRAY, NULL);
                 if (indices == NULL) return NULL;
 	}
 	else {
-		if(!PyArg_ParseTupleAndKeywords(args, kwds, "O|iO&", kwlist1,
+		if(!PyArg_ParseTupleAndKeywords(args, kwds, "O|iO&O&", kwlist1,
 						&op, &axis, 
 						PyArray_DescrConverter2, 
-						&otype)) return NULL;
+						&otype,
+                                                PyArray_OutputConverter,
+                                                &out)) return NULL;
 	}
 	
 	/* Ensure input is an array */	
@@ -2511,7 +2552,14 @@ PyUFunc_GenericReduction(PyUFuncObject *self, PyObject *args,
 		return NULL;
 	}
 
-	/* Get default type to reduce over if not given */
+	/* If out is specified it determines otype unless otype
+           already specified.
+         */
+        if (otype == NULL && out != NULL) {
+                otype = out->descr;
+                Py_INCREF(otype);
+        }
+
         if (otype == NULL) {
 		/* For integer types --- makes sure at 
 		   least a long is used */
@@ -2533,15 +2581,15 @@ PyUFunc_GenericReduction(PyUFuncObject *self, PyObject *args,
 
         switch(operation) {
         case UFUNC_REDUCE:
-                ret = (PyArrayObject *)PyUFunc_Reduce(self, mp, axis, 
+                ret = (PyArrayObject *)PyUFunc_Reduce(self, mp, out, axis, 
                                                       otype->type_num);
 		break;
         case UFUNC_ACCUMULATE:
-                ret = (PyArrayObject *)PyUFunc_Accumulate(self, mp, axis, 
+                ret = (PyArrayObject *)PyUFunc_Accumulate(self, mp, out, axis, 
                                                           otype->type_num);
 		break;
         case UFUNC_REDUCEAT:
-                ret = (PyArrayObject *)PyUFunc_Reduceat(self, mp, indices, 
+                ret = (PyArrayObject *)PyUFunc_Reduceat(self, mp, indices, out, 
                                                         axis, otype->type_num);
                 Py_DECREF(indices);
 		break;
diff --git a/numpy/numarray/__init__.py b/numpy/numarray/__init__.py
index e48401433..6a902ddee 100644
--- a/numpy/numarray/__init__.py
+++ b/numpy/numarray/__init__.py
@@ -2,21 +2,24 @@ from util import *
 from numerictypes import *
 from functions import *
 from ufuncs import *
+from session import *
 
 import util
 import numerictypes
 import functions
 import ufuncs
 import compat
+import session
 
-__all__ = util.__all__
+__all__ = ['session', 'numerictypes']
+__all__ += util.__all__
 __all__ += numerictypes.__all__
 __all__ += functions.__all__
 __all__ += ufuncs.__all__
 __all__ += compat.__all__
+__all__ += session.__all__
 
 del util
-del numerictypes
 del functions
 del ufuncs
 del compat
diff --git a/numpy/numarray/compat.py b/numpy/numarray/compat.py
index c9fd3ef86..e0d13a7c2 100644
--- a/numpy/numarray/compat.py
+++ b/numpy/numarray/compat.py
@@ -1,6 +1,4 @@
 
-__all__ = ['NewAxis']
+__all__ = ['NewAxis', 'ArrayType']
 
-from numpy import newaxis
-
-NewAxis = newaxis
+from numpy import newaxis as NewAxis, ndarray as ArrayType
diff --git a/numpy/numarray/functions.py b/numpy/numarray/functions.py
index d549f938f..1f358385c 100644
--- a/numpy/numarray/functions.py
+++ b/numpy/numarray/functions.py
@@ -1,71 +1,80 @@
 
 # missing Numarray defined names (in from numarray import *)
-##__all__ = ['ArrayType', 'CLIP', 'ClassicUnpickler', 'Complex32_fromtype',
-##           'Complex64_fromtype', 'ComplexArray', 'EarlyEOFError', 'Error',
-##           'FileSeekWarning', 'MAX_ALIGN', 'MAX_INT_SIZE', 'MAX_LINE_WIDTH',
-##           'MathDomainError', 'NDArray', 'NewArray', 'NewAxis', 'NumArray',
-##           'NumError', 'NumOverflowError', 'PRECISION', 'Py2NumType',
+##__all__ = ['ClassicUnpickler', 'Complex32_fromtype',
+##           'Complex64_fromtype', 'ComplexArray', 'Error',
+##           'MAX_ALIGN', 'MAX_INT_SIZE', 'MAX_LINE_WIDTH',
+##           'NDArray', 'NewArray', 'NumArray',
+##           'NumError', 'PRECISION', 'Py2NumType',
 ##           'PyINT_TYPES', 'PyLevel2Type', 'PyNUMERIC_TYPES', 'PyREAL_TYPES',
-##           'RAISE', 'SLOPPY', 'STRICT', 'SUPPRESS_SMALL', 'SizeMismatchError',
-##           'SizeMismatchWarning', 'SuitableBuffer', 'USING_BLAS',
-##           'UnderflowError', 'UsesOpPriority', 'WARN', 'WRAP', 'all',
-##           'allclose', 'alltrue', 'and_', 'any', 'arange', 'argmax',
-##           'argmin', 'argsort', 'around', 'array2list', 'array_equal',
-##           'array_equiv', 'array_repr', 'array_str', 'arrayprint',
-##           'arrayrange', 'average', 'choose', 'clip',
-##           'codegenerator', 'compress', 'concatenate', 'conjugate',
-##           'copy', 'copy_reg', 'diagonal', 'divide_remainder',
-##           'dotblas', 'e', 'explicit_type', 'flush_caches', 'fromfile',
-##           'fromfunction', 'fromlist', 'fromstring', 'generic',
-##           'genericCoercions', 'genericPromotionExclusions', 'genericTypeRank',
-##           'getShape', 'getTypeObject', 'handleError', 'identity', 'indices',
-##           'info', 'innerproduct', 'inputarray', 'isBigEndian',
-##           'kroneckerproduct', 'lexsort', 'libnumarray', 'libnumeric',
-##           'load', 'make_ufuncs', 'math', 'memory',
-##           'numarrayall', 'numarraycore', 'numerictypes', 'numinclude',
-##           'operator', 'os', 'outerproduct', 'pi', 'put', 'putmask',
-##           'pythonTypeMap', 'pythonTypeRank', 'rank', 'repeat',
-##           'reshape', 'resize', 'round', 'safethread', 'save', 'scalarTypeMap',
-##           'scalarTypes', 'searchsorted', 'session', 'shape', 'sign', 'size',
-##           'sometrue', 'sort', 'swapaxes', 'sys', 'take', 'tcode',
-##           'tensormultiply', 'tname', 'trace', 'transpose', 'typeDict',
-##           'typecode', 'typecodes', 'typeconv', 'types', 'ufunc',
-##           'ufuncFactory', 'value', 'ieeemask', 'cumproduct', 'cumsum',
-##           'nonzero']
-
+##           'SUPPRESS_SMALL',
+##           'SuitableBuffer', 'USING_BLAS',
+##           'UsesOpPriority', 
+##           'codegenerator', 'generic', 'libnumarray', 'libnumeric',
+##           'make_ufuncs', 'memory',
+##           'numarrayall', 'numarraycore', 'numinclude', 'safethread',
+##           'typecode', 'typecodes', 'typeconv', 'ufunc', 'ufuncFactory',
+##           'ieeemask']
 
 __all__ = ['asarray', 'ones', 'zeros', 'array', 'where']
 __all__ += ['vdot', 'dot', 'matrixmultiply', 'ravel', 'indices',
-            'arange', 'concatenate']
+            'arange', 'concatenate', 'all', 'allclose', 'alltrue', 'and_',
+            'any', 'argmax', 'argmin', 'argsort', 'around', 'array_equal',
+            'array_equiv', 'arrayrange', 'array_str', 'array_repr',
+            'array2list', 'average', 'choose', 'CLIP', 'RAISE', 'WRAP',
+            'clip', 'compress', 'concatenate', 'copy', 'copy_reg',
+            'diagonal', 'divide_remainder', 'e', 'explicit_type', 'pi',
+            'flush_caches', 'fromfile', 'os', 'sys', 'STRICT',
+            'SLOPPY', 'WARN', 'EarlyEOFError', 'SizeMismatchError',
+            'SizeMismatchWarning', 'FileSeekWarning', 'fromstring',
+            'fromfunction', 'fromlist', 'getShape', 'getTypeObject',
+            'identity', 'indices', 'info', 'innerproduct', 'inputarray',
+            'isBigEndian', 'kroneckerproduct', 'lexsort', 'math',
+            'operator', 'outerproduct', 'put', 'putmask', 'rank',
+            'repeat', 'reshape', 'resize', 'round', 'searchsorted',
+            'shape', 'size', 'sometrue', 'sort', 'swapaxes', 'take',
+            'tcode', 'tname', 'tensormultiply', 'trace', 'transpose',
+            'types', 'value', 'cumsum', 'cumproduct', 'nonzero'
+            ]
 
-from numpy import dot as matrixmultiply, dot, vdot, ravel, concatenate
+import copy, copy_reg, types
+import os, sys, math, operator
 
-def array(sequence=None, typecode=None, copy=1, savespace=0,
-          type=None, shape=None, dtype=None):
-    dtype = type2dtype(typecode, type, dtype)
-    if sequence is None:
-        if shape is None:
-            return None
-        if dtype is None:
-            dtype = 'l'
-        return N.empty(shape, dtype)
-    arr = N.array(sequence, dtype, copy=copy)
-    if shape is not None:
-        arr.shape = shape
-    return arr
+from numpy import dot as matrixmultiply, dot, vdot, ravel, concatenate, all,\
+     allclose, any, argmax, argmin, around, argsort, array_equal, array_equiv,\
+     array_str, array_repr, average, CLIP, RAISE, WRAP, clip, concatenate, \
+     diagonal, e, pi, fromfunction, indices, inner as innerproduct, nonzero, \
+     outer as outerproduct, kron as kroneckerproduct, lexsort, putmask, rank, \
+     resize, searchsorted, shape, size, sort, swapaxes, trace, transpose
+import numpy as N
 
-def asarray(seq, type=None, typecode=None, dtype=None):
-    if seq is None:
-        return None
-    dtype = type2dtype(typecode, type, dtype)
-    return N.array(seq, dtype, copy=0)
+from numerictypes import typefrom
+
+isBigEndian = sys.byteorder != 'little'
+value = tcode = 'f'
+tname = 'Float32'
 
+#  If dtype is not None, then it is used
+#  If type is not None, then it is used
+#  If typecode is not None then it is used
+#  If use_default is True, then the default
+#   data-type is returned if all are None
+def type2dtype(typecode, type, dtype, use_default=True):
+    if dtype is None:
+        if type is None:
+            if use_default or typecode is not None:
+                dtype = N.dtype(typecode)
+        else:
+            dtype = N.dtype(type)
+    if use_default and dtype is None:
+        dtype = N.dtype(None)
+    return dtype
+    
 def ones(shape, type=None, typecode=None, dtype=None):
-    dtype = type2dtype(typecode, type, dtype)
+    dtype = type2dtype(typecode, type, dtype, 1)
     return N.ones(shape, dtype)
 
 def zeros(shape, type=None, typecode=None, dtype=None):
-    dtype = type2dtype(typecode, type, dtype)
+    dtype = type2dtype(typecode, type, dtype, 1)
     return N.zeros(shape, dtype)
 
 def where(condition, x=None, y=None, out=None):
@@ -83,5 +92,349 @@ def indices(shape, type=None):
 
 def arange(a1, a2=None, stride=1, type=None, shape=None,
            typecode=None, dtype=None):
-    dtype = type2dtype(typecode, type, dtype)
+    dtype = type2dtype(typecode, type, dtype, 0)
     return N.arange(a1, a2, stride, dtype)
+
+arrayrange = arange
+
+def alltrue(x, axis=0):
+    return N.alltrue(x, axis)
+
+def and_(a, b):
+    """Same as a & b
+    """
+    return a & b
+
+def divide_remainder(a, b):
+    a, b = asarray(a), asarray(b)
+    return (a/b,a%b)
+
+def around(array, digits=0, output=None):
+    ret = N.around(array, digits, output)
+    if output is None:
+        return ret
+    return
+    
+def array2list(arr):
+    return arr.tolist()
+
+
+def choose(selector, population, outarr=None, clipmode=RAISE):
+    a = N.asarray(selector)
+    ret = a.choose(population, out=outarr, mode=clipmode)
+    if outarr is None:
+        return ret
+    return
+
+def compress(condition, a, axis=0):
+    return N.compress(condition, a, axis)
+
+# only returns a view
+def explicit_type(a):
+    x = a.view()
+    return x
+
+# stub
+def flush_caches():
+    pass
+
+
+class EarlyEOFError(Exception):
+    "Raised in fromfile() if EOF unexpectedly occurs."
+    pass
+
+class SizeMismatchError(Exception):
+    "Raised in fromfile() if file size does not match shape."
+    pass
+
+class SizeMismatchWarning(Warning):
+    "Issued in fromfile() if file size does not match shape."
+    pass
+
+class FileSeekWarning(Warning):
+    "Issued in fromfile() if there is unused data and seek() fails"
+    pass
+
+
+STRICT, SLOPPY, WARN = range(3)
+
+_BLOCKSIZE=1024
+
+# taken and adapted directly from numarray
+def fromfile(infile, type=None, shape=None, sizing=STRICT,
+             typecode=None, dtype=None):
+    if isinstance(infile, (str, unicode)):
+        infile = open(infile, 'rb')
+    dtype = type2dtype(typecode, type, dtype, True)
+    if shape is None:
+        shape = (-1,)
+    if not isinstance(shape, tuple):
+        shape = (shape,)
+        
+    if (list(shape).count(-1)>1):
+        raise ValueError("At most one unspecified dimension in shape")
+
+    if -1 not in shape:
+        if sizing != STRICT:
+            raise ValueError("sizing must be STRICT if size complete")
+        arr = N.empty(shape, dtype)
+        bytesleft=arr.nbytes
+        bytesread=0
+        while(bytesleft > _BLOCKSIZE):
+            data = infile.read(_BLOCKSIZE)
+            if len(data) != _BLOCKSIZE:
+                raise EarlyEOFError("Unexpected EOF reading data for size complete array")
+            arr.data[bytesread:bytesread+_BLOCKSIZE]=data
+            bytesread += _BLOCKSIZE
+            bytesleft -= _BLOCKSIZE
+        if bytesleft > 0:
+            data = infile.read(bytesleft)
+            if len(data) != bytesleft:
+                raise EarlyEOFError("Unexpected EOF reading data for size complete array")
+            arr.data[bytesread:bytesread+bytesleft]=data
+        return arr
+
+
+    ##shape is incompletely specified
+    ##read until EOF
+    ##implementation 1: naively use memory blocks
+    ##problematic because memory allocation can be double what is
+    ##necessary (!)
+
+    ##the most common case, namely reading in data from an unchanging
+    ##file whose size may be determined before allocation, should be
+    ##quick -- only one allocation will be needed.
+    
+    recsize = dtype.itemsize * N.product([i for i in shape if i != -1])
+    blocksize = max(_BLOCKSIZE/recsize, 1)*recsize
+
+    ##try to estimate file size
+    try:
+        curpos=infile.tell()
+        infile.seek(0,2)
+        endpos=infile.tell()
+        infile.seek(curpos)
+    except (AttributeError, IOError):
+        initsize=blocksize
+    else:
+        initsize=max(1,(endpos-curpos)/recsize)*recsize
+
+    buf = N.newbuffer(initsize)
+
+    bytesread=0
+    while 1:
+        data=infile.read(blocksize)
+        if len(data) != blocksize: ##eof
+            break
+        ##do we have space?
+        if len(buf) < bytesread+blocksize:
+            buf=_resizebuf(buf,len(buf)+blocksize)
+            ## or rather a=resizebuf(a,2*len(a)) ?
+        assert len(buf) >= bytesread+blocksize
+        buf[bytesread:bytesread+blocksize]=data
+        bytesread += blocksize
+
+    if len(data) % recsize != 0:
+        if sizing == STRICT:
+            raise SizeMismatchError("Filesize does not match specified shape")
+        if sizing == WARN:
+            _warnings.warn("Filesize does not match specified shape",
+                           SizeMismatchWarning)
+        try:
+            infile.seek(-(len(data) % recsize),1)
+        except AttributeError:
+            _warnings.warn("Could not rewind (no seek support)",
+                           FileSeekWarning)
+        except IOError:
+            _warnings.warn("Could not rewind (IOError in seek)",
+                           FileSeekWarning)
+    datasize = (len(data)/recsize) * recsize
+    if len(buf) != bytesread+datasize:
+        buf=_resizebuf(buf,bytesread+datasize)
+    buf[bytesread:bytesread+datasize]=data[:datasize]
+    ##deduce shape from len(buf)
+    shape = list(shape)
+    uidx = shape.index(-1)
+    shape[uidx]=len(buf) / recsize
+
+    a = N.ndarray(shape=shape, dtype=type, buffer=buf)
+    if a.dtype.char == '?':
+        N.not_equal(a, 0, a)
+    return a
+    
+def fromstring(datastring, type=None, shape=None, typecode=None, dtype=None):
+    dtype = type2dtype(typecode, type, dtype, True)
+    if shape is None:
+        count = -1
+    else:
+        count = N.product(shape)*dtype.itemsize
+    res = N.fromstring(datastring, count=count)
+    if shape is not None:
+        res.shape = shape
+    return res                
+
+
+# check_overflow is ignored
+def fromlist(seq, type=None, shape=None, check_overflow=0, typecode=None, dtype=None):
+    dtype = type2dtype(typecode, type, dtype, False)
+    return N.array(seq, dtype)
+
+def array(sequence=None, typecode=None, copy=1, savespace=0,
+          type=None, shape=None, dtype=None):
+    dtype = type2dtype(typecode, type, dtype, 0)
+    if sequence is None:
+        if shape is None:
+            return None
+        if dtype is None:
+            dtype = 'l'
+        return N.empty(shape, dtype)
+    if isinstance(sequence, file):
+        return fromfile(sequence, dtype=dtype, shape=shape)
+    if isinstance(sequence, str):
+        return fromstring(sequence, dtype=dtype, shape=shape)
+    if isinstance(sequence, buffer):
+        arr = N.frombuffer(sequence, dtype=dtype)
+    else:
+        arr = N.array(sequence, dtype, copy=copy)
+    if shape is not None:
+        arr.shape = shape
+    return arr
+
+def asarray(seq, type=None, typecode=None, dtype=None):
+    if isinstance(seq, N.ndarray) and type is None and \
+           typecode is None and dtype is None:
+        return seq
+    return array(seq, type=type, typecode=typecode, copy=0, dtype=dtype)
+
+inputarray = asarray
+
+
+def getTypeObject(sequence, type):
+    if type is not None:
+        return type
+    try:
+        return typefrom(N.array(sequence))
+    except:
+        raise TypeError("Can't determine a reasonable type from sequence")
+
+def getShape(shape, *args):
+    try:
+        if shape is () and not args:
+            return ()
+        if len(args) > 0:
+            shape = (shape, ) + args
+        else:
+            shape = tuple(shape)
+        dummy = N.array(shape)
+        if not issubclass(dummy.dtype.type, N.integer):
+            raise TypeError
+        if len(dummy) > N.MAXDIMS:
+            raise TypeError        
+    except:
+        raise TypeError("Shape must be a sequence of integers")
+    return shape
+
+    
+def identity(n, type=None, typecode=None, dtype=None):
+    dtype = type2dtype(typecode, type, dtype, True)
+    return N.identity(n, dtype)
+
+def info(obj):
+    print "class: ", type(obj)
+    print "shape: ", obj.shape
+    print "strides: ", obj.strides
+    print "byteoffset: 0"
+    print "bytestride: ", obj.strides[0]
+    print "itemsize: ", obj.itemsize
+    print "aligned: ", obj.flags.isaligned
+    print "contiguous: ", obj.flags.contiguous
+    print "buffer: ", obj.data
+    print "data pointer:", obj._as_paramater_, "(DEBUG ONLY)"
+    print "byteorder: ",
+    endian = obj.dtype.byteorder
+    if endian in ['|','=']:
+        print sys.byteorder
+    elif endian == '>':
+        print "big"
+    else:
+        print "little"
+    print "byteswap: ", not obj.dtype.isnative
+    print "type: ", typefrom(obj)
+
+#clipmode is ignored if axis is not 0 and array is not 1d
+def put(array, indices, values, axis=0, clipmode=RAISE):
+    if not isinstance(array, N.ndarray):
+        raise TypeError("put only works on subclass of ndarray")
+    work = asarray(array)
+    if axis == 0:
+        if array.ndim == 1:
+            work.put(indices, values, clipmode)
+        else:
+            work[indices] = values
+    elif isinstance(axis, (int, long, N.integer)):
+        work = work.swapaxes(0, axis)
+        work[indices] = values
+        work = work.swapaxes(0, axis)
+    else:
+        def_axes = range(work.ndim)
+        for x in axis:
+            def_axes.remove(x)
+        axis = list(axis)+def_axes
+        work = work.transpose(axis)
+        work[indices] = values
+        work = work.transpose(axis)
+
+def repeat(array, repeats, axis=0):
+    return N.repeat(array, repeats, axis)
+
+
+def reshape(array, shape, *args):
+    if len(args) > 0:
+        shape = (shape,) + args
+    return N.reshape(array, shape)
+
+
+import warnings as _warnings
+def round(*args, **keys):
+    _warnings.warn("round() is deprecated. Switch to around()",
+                   DeprecationWarning)
+    return around(*args, **keys)
+
+def sometrue(array, axis=0):
+    return N.sometrue(array, axis)
+
+#clipmode is ignored if axis is not an integer
+def take(array, indices, axis=0, outarr=None, clipmode=RAISE):
+    array = N.asarray(array)
+    if isinstance(axis, (int, long, N.integer)):
+        res = array.take(indices, axis, outarr, clipmode)
+        if outarr is None:
+            return res
+        return
+    else:
+        def_axes = range(array.ndim)
+        for x in axis:
+            def_axes.remove(x)
+        axis = list(axis) + def_axes
+        work = array.transpose(axis)
+        res = work[indices]
+        if outarr is None:
+            return res
+        out[...] = res
+        return 
+    
+def tensormultiply(a1, a2):
+    a1, a2 = N.asarray(a1), N.asarray(a2)
+    if (a1.shape[-1] != a2.shape[0]):
+        raise ValueError("Unmatched dimensions")
+    shape = a1.shape[:-1] + a2.shape[1:]
+    return N.reshape(dot(N.reshape(a1, (-1, a1.shape[-1])),
+                         N.reshape(a2, (a2.shape[0],-1))),
+                     shape)
+
+def cumsum(a1, axis=0, out=None, type=None, dim=0):
+    return N.asarray(a1).cumsum(axis,dtype=type,out=out)
+
+def cumproduct(a1, axis=0, out=None, type=None, dim=0):
+    return N.asarray(a1).cumprod(axis,dtype=type,out=out)
+
diff --git a/numpy/numarray/numerictypes.py b/numpy/numarray/numerictypes.py
index 549fb8ca2..eadaaeee3 100644
--- a/numpy/numarray/numerictypes.py
+++ b/numpy/numarray/numerictypes.py
@@ -530,7 +530,6 @@ _scipy_dtypechar_inverse = {}
 for key,value in _scipy_dtypechar.items():
     _scipy_dtypechar_inverse[value] = key
 
-
 def typefrom(obj):
     return _scipy_dtypechar_inverse[obj.dtype.char]
 
diff --git a/numpy/numarray/session.py b/numpy/numarray/session.py
new file mode 100644
index 000000000..25520bd41
--- /dev/null
+++ b/numpy/numarray/session.py
@@ -0,0 +1,349 @@
+""" This module contains a "session saver" which saves the state of a
+NumPy session to a file.  At a later time, a different Python
+process can be started and the saved session can be restored using
+load().
+
+The session saver relies on the Python pickle protocol to save and
+restore objects.  Objects which are not themselves picklable (e.g.
+modules) can sometimes be saved by "proxy",  particularly when they
+are global constants of some kind.  If it's not known that proxying
+will work,  a warning is issued at save time.  If a proxy fails to
+reload properly (e.g. because it's not a global constant),  a warning
+is issued at reload time and that name is bound to a _ProxyFailure
+instance which tries to identify what should have been restored.
+
+First, some unfortunate (probably unnecessary) concessions to doctest
+to keep the test run free of warnings.
+
+>>> del _PROXY_ALLOWED
+>>> del copy
+>>> del __builtins__
+
+By default, save() stores every variable in the caller's namespace:
+
+>>> import numpy as na
+>>> a = na.arange(10)
+>>> save()
+
+Alternately,  save() can be passed a comma seperated string of variables:
+
+>>> save("a,na")
+
+Alternately,  save() can be passed a dictionary, typically one you already
+have lying around somewhere rather than created inline as shown here:
+
+>>> save(dictionary={"a":a,"na":na})
+
+If both variables and a dictionary are specified, the variables to be
+saved are taken from the dictionary.
+
+>>> save(variables="a,na",dictionary={"a":a,"na":na})
+
+Remove names from the session namespace
+
+>>> del a, na
+
+By default, load() restores every variable/object in the session file
+to the caller's namespace.
+
+>>> load()
+
+load() can be passed a comma seperated string of variables to be
+restored from the session file to the caller's namespace:
+
+>>> load("a,na")
+
+load() can also be passed a dictionary to *restore to*:
+
+>>> d = {}
+>>> load(dictionary=d)
+
+load can be passed both a list variables of variables to restore and a
+dictionary to restore to:
+
+>>> load(variables="a,na", dictionary=d)
+
+>>> na.all(a == na.arange(10))
+1
+>>> na.__name__
+'numpy'
+
+NOTE:  session saving is faked for modules using module proxy objects.
+Saved modules are re-imported at load time but any "state" in the module
+which is not restored by a simple import is lost.
+
+"""
+
+__all__ = ['load', 'save']
+
+import copy
+import sys
+import pickle
+
+SAVEFILE="session.dat"
+VERBOSE = False           # global import-time  override
+
+def _foo(): pass
+
+_PROXY_ALLOWED = (type(sys),  # module
+                  type(_foo), # function
+                  type(None)) # None
+
+def _update_proxy_types():
+    """Suppress warnings for known un-picklables with working proxies."""
+    pass
+
+def _unknown(_type):
+    """returns True iff _type isn't known as OK to proxy"""
+    return (_type is not None) and (_type not in _PROXY_ALLOWED)
+
+# caller() from the following article with one extra f_back added.
+# from http://www.python.org/search/hypermail/python-1994q1/0506.html
+# SUBJECT: import ( how to put a symbol into caller's namespace )
+# SENDER:  Steven D. Majewski (sdm7g@elvis.med.virginia.edu)
+# DATE:  Thu, 24 Mar 1994 15:38:53 -0500
+
+def _caller():
+    """caller() returns the frame object of the function's caller."""
+    try:
+        1 + '' # make an error happen
+    except: # and return the caller's caller's frame
+        return sys.exc_traceback.tb_frame.f_back.f_back.f_back
+
+def _callers_globals():
+    """callers_globals() returns the global dictionary of the caller."""
+    frame = _caller()
+    return frame.f_globals
+
+def _callers_modules():
+    """returns a list containing the names of all the modules in the caller's
+    global namespace."""
+    g = _callers_globals()
+    mods = []
+    for k,v in g.items():
+        if type(v) == type(sys):
+            mods.append(getattr(v,"__name__"))
+    return mods
+
+def _errout(*args):
+    for a in args:
+        print >>sys.stderr, a,
+    print >>sys.stderr
+
+def _verbose(*args):
+    if VERBOSE:
+        _errout(*args)
+
+class _ProxyingFailure:
+    """Object which is bound to a variable for a proxy pickle which failed to reload"""
+    def __init__(self, module, name, type=None):
+        self.module = module
+        self.name = name
+        self.type = type
+    def __repr__(self):
+        return "ProxyingFailure('%s','%s','%s')" % (self.module, self.name, self.type)
+
+class _ModuleProxy(object):
+    """Proxy object which fakes pickling a module"""
+    def __new__(_type, name, save=False):
+        if save:
+            _verbose("proxying module", name)
+            self = object.__new__(_type)
+            self.name = name
+        else:
+            _verbose("loading module proxy", name)
+            try:
+                self = _loadmodule(name)
+            except ImportError:
+                _errout("warning: module", name,"import failed.")
+        return self
+
+    def __getnewargs__(self):
+        return (self.name,)
+
+    def __getstate__(self):
+        return False
+
+def _loadmodule(module):
+    if not sys.modules.has_key(module):
+        modules = module.split(".")
+        s = ""
+        for i in range(len(modules)):
+            s = ".".join(modules[:i+1])
+            exec "import " + s
+    return sys.modules[module]
+
+class _ObjectProxy(object):
+    """Proxy object which fakes pickling an arbitrary object.  Only global
+    constants can really be proxied."""
+    def __new__(_type, module, name, _type2, save=False):
+        if save:
+            if _unknown(_type2):
+                _errout("warning: proxying object", module + "." + name,
+                        "of type", _type2, "because it wouldn't pickle...",
+                        "it may not reload later.")
+            else:
+                _verbose("proxying object", module, name)
+            self = object.__new__(_type)
+            self.module, self.name, self.type = module, name, str(_type2)
+        else:
+            _verbose("loading object proxy", module, name)
+            try:
+                m = _loadmodule(module)
+            except (ImportError, KeyError):
+                _errout("warning: loading object proxy", module + "." + name,
+                        "module import failed.")
+                return _ProxyingFailure(module,name,_type2)
+            try:
+                self = getattr(m, name)
+            except AttributeError:
+                _errout("warning: object proxy", module + "." + name,
+                        "wouldn't reload from", m)
+                return _ProxyingFailure(module,name,_type2)
+        return self
+
+    def __getnewargs__(self):
+        return (self.module, self.name, self.type)
+
+    def __getstate__(self):
+        return False
+
+
+class _SaveSession(object):    
+    """Tag object which marks the end of a save session and holds the
+    saved session variable names as a list of strings in the same
+    order as the session pickles."""
+    def __new__(_type, keys, save=False):
+        if save:
+            _verbose("saving session", keys)
+        else:
+            _verbose("loading session", keys)
+        self = object.__new__(_type)
+        self.keys = keys
+        return self
+
+    def __getnewargs__(self):
+        return (self.keys,)
+
+    def __getstate__(self):
+        return False
+
+class ObjectNotFound(RuntimeError):
+    pass
+
+def _locate(modules, object):
+    for mname in modules:
+        m = sys.modules[mname]
+        if m:
+            for k,v in m.__dict__.items():
+                if v is object:
+                    return m.__name__, k
+    else:
+        raise ObjectNotFound(k)
+
+def save(variables=None, file=SAVEFILE, dictionary=None, verbose=False):
+
+    """saves variables from a numpy session to a file.  Variables
+    which won't pickle are "proxied" if possible.
+
+    'variables'       a string of comma seperated variables: e.g. "a,b,c"
+                      Defaults to dictionary.keys().
+
+    'file'            a filename or file object for the session file.
+    
+    'dictionary'      the dictionary in which to look up the variables.
+                      Defaults to the caller's globals()
+    
+    'verbose'         print additional debug output when True.
+    """
+    
+    global VERBOSE
+    VERBOSE = verbose
+
+    _update_proxy_types()
+
+    if isinstance(file, str):
+        file = open(file, "wb")
+        
+    if dictionary is None:
+        dictionary = _callers_globals()
+
+    if variables is None:
+        keys = dictionary.keys()
+    else:
+        keys = variables.split(",")
+        
+    source_modules = _callers_modules() + sys.modules.keys()
+
+    p = pickle.Pickler(file, protocol=2)
+
+    _verbose("variables:",keys)
+    for k in keys:
+        v = dictionary[k]
+        _verbose("saving", k, type(v))
+        try:  # Try to write an ordinary pickle
+            p.dump(v)
+            _verbose("pickled", k)
+        except (pickle.PicklingError, TypeError, SystemError):
+            # Use proxies for stuff that won't pickle
+            if isinstance(v, type(sys)): # module
+                proxy = _ModuleProxy(v.__name__, save=True)
+            else:
+                try:
+                    module, name = _locate(source_modules, v)
+                except ObjectNotFound:
+                    _errout("warning: couldn't find object",k,
+                            "in any module... skipping.")
+                    continue
+                else:
+                    proxy = _ObjectProxy(module, name, type(v), save=True)
+            p.dump(proxy)
+    o = _SaveSession(keys, save=True)
+    p.dump(o)
+    file.close()
+
+def load(variables=None, file=SAVEFILE, dictionary=None, verbose=False):
+    
+    """load a numpy session from a file and store the specified
+    'variables' into 'dictionary'.
+
+    'variables'       a string of comma seperated variables: e.g. "a,b,c"
+                      Defaults to dictionary.keys().
+
+    'file'            a filename or file object for the session file.
+    
+    'dictionary'      the dictionary in which to look up the variables.
+                      Defaults to the caller's globals()
+
+    'verbose'         print additional debug output when True.
+    """
+    
+    global VERBOSE
+    VERBOSE = verbose
+            
+    if isinstance(file, str):
+        file = open(file, "rb")
+    if dictionary is None:
+        dictionary = _callers_globals()
+    values = []
+    p = pickle.Unpickler(file)
+    while 1:
+        o = p.load()
+        if isinstance(o, _SaveSession):
+            session = dict(zip(o.keys, values))
+            _verbose("updating dictionary with session variables.")
+            if variables is None:
+                keys = session.keys()
+            else:
+                keys = variables.split(",")
+            for k in keys:
+                dictionary[k] = session[k]
+            return None
+        else:
+            _verbose("unpickled object", str(o))
+            values.append(o)
+
+def test():
+    import doctest, numpy.numarray.session
+    return doctest.testmod(numpy.numarray.session)
+
diff --git a/numpy/numarray/ufuncs.py b/numpy/numarray/ufuncs.py
index 685b76587..3fb5671ce 100644
--- a/numpy/numarray/ufuncs.py
+++ b/numpy/numarray/ufuncs.py
@@ -8,7 +8,8 @@ __all__ = ['abs', 'absolute', 'add', 'arccos', 'arccosh', 'arcsin', 'arcsinh',
            'logical_or', 'logical_xor', 'lshift', 'maximum', 'minimum',
            'minus', 'multiply', 'negative', 'not_equal',
            'power', 'product', 'remainder', 'rshift', 'sin', 'sinh', 'sqrt',
-           'subtract', 'sum', 'tan', 'tanh', 'true_divide']
+           'subtract', 'sum', 'tan', 'tanh', 'true_divide',
+           'conjugate', 'sign']
 
 from numpy import absolute as abs, absolute, add, arccos, arccosh, arcsin, \
      arcsinh,  arctan, arctan2, arctanh, bitwise_and, invert as bitwise_not, \
@@ -18,4 +19,4 @@ from numpy import absolute as abs, absolute, add, arccos, arccosh, arcsin, \
      logical_not, logical_or, logical_xor, left_shift as lshift, \
      maximum, minimum, negative as minus, multiply, negative, \
      not_equal, power, product, remainder, right_shift as rshift, sin, \
-     sinh, sqrt, subtract, sum, tan, tanh, true_divide
+     sinh, sqrt, subtract, sum, tan, tanh, true_divide, conjugate, sign