Merge pull request #339 from charris/generic-sorts

Add generic versions of heapsort and mergesort so that those sorts will be available to all numpy types. The generic versions are probably not as fast as the type specific versions, but they will always be there.

12 files changed, 1427 insertions, 875 deletions

diff --git a/numpy/core/SConscript b/numpy/core/SConscript
index 46087610c..952544a25 100644
--- a/numpy/core/SConscript
+++ b/numpy/core/SConscript
@@ -429,7 +429,9 @@ env.Install('$distutils_installdir/lib/npy-pkg-config', mlib_ini)
 env.Install('$distutils_installdir/lib/npy-pkg-config', npymath_ini)
 
 # npysort core lib
-npysort_src = [env.GenerateFromTemplate(pjoin('src', 'npysort', 'sort.c.src'))]
+npysort_src = [env.GenerateFromTemplate(pjoin('src', 'npysort', 'quicksort.c.src')),
+               env.GenerateFromTemplate(pjoin('src', 'npysort', 'mergesort.c.src')),
+               env.GenerateFromTemplate(pjoin('src', 'npysort', 'heapsort.c.src'))]
 env.StaticExtLibrary("npysort", npysort_src)
 env.Prepend(LIBS=["npysort"])
 env.Prepend(LIBPATH=["."])
diff --git a/numpy/core/bento.info b/numpy/core/bento.info
index 04401f991..bf8ae3a81 100644
--- a/numpy/core/bento.info
+++ b/numpy/core/bento.info
@@ -10,7 +10,9 @@ Library:
             src/npymath/halffloat.c
     CompiledLibrary: npysort
         Sources:
-            src/npysort/sort.c.src
+            src/npysort/quicksort.c.src,
+            src/npysort/mergesort.c.src,
+            src/npysort/heapsort.c.src
     Extension: multiarray
         Sources:
             src/multiarray/multiarraymodule_onefile.c
diff --git a/numpy/core/setup.py b/numpy/core/setup.py
index 2d36f5ee1..a1000ae0b 100644
--- a/numpy/core/setup.py
+++ b/numpy/core/setup.py
@@ -671,7 +671,10 @@ def configuration(parent_package='',top_path=None):
 
     # This library is created for the build but it is not installed
     config.add_library('npysort',
-            sources = [join('src', 'npysort', 'sort.c.src')])
+            sources = [join('src', 'npysort', 'quicksort.c.src'),
+                       join('src', 'npysort', 'mergesort.c.src'),
+                       join('src', 'npysort', 'heapsort.c.src')])
+
 
     #######################################################################
     #                        multiarray module                            #
diff --git a/numpy/core/src/multiarray/item_selection.c b/numpy/core/src/multiarray/item_selection.c
index 977580d70..fe28f2482 100644
--- a/numpy/core/src/multiarray/item_selection.c
+++ b/numpy/core/src/multiarray/item_selection.c
@@ -19,6 +19,7 @@
 #include "lowlevel_strided_loops.h"
 
 #include "item_selection.h"
+#include "npy_sort.h"
 
 /*NUMPY_API
  * Take
@@ -890,7 +891,7 @@ _new_argsort(PyArrayObject *op, int axis, NPY_SORTKIND which)
 static PyArrayObject *global_obj;
 
 static int
-qsortCompare (const void *a, const void *b)
+sortCompare (const void *a, const void *b)
 {
     return PyArray_DESCR(global_obj)->f->compare(a,b,global_obj);
 }
@@ -954,7 +955,9 @@ PyArray_Sort(PyArrayObject *op, int axis, NPY_SORTKIND which)
     PyArrayObject *ap = NULL, *store_arr = NULL;
     char *ip;
     int i, n, m, elsize, orign;
-    int axis_orig=axis;
+    int res = 0;
+    int axis_orig = axis;
+    int (*sort)(void *, size_t, size_t, npy_comparator);
 
     n = PyArray_NDIM(op);
     if ((n == 0) || (PyArray_SIZE(op) == 1)) {
@@ -975,13 +978,30 @@ PyArray_Sort(PyArrayObject *op, int axis, NPY_SORTKIND which)
     if (PyArray_DESCR(op)->f->sort[which] != NULL) {
         return _new_sort(op, axis, which);
     }
-    if ((which != NPY_QUICKSORT)
-        || PyArray_DESCR(op)->f->compare == NULL) {
+
+    if (PyArray_DESCR(op)->f->compare == NULL) {
         PyErr_SetString(PyExc_TypeError,
-                        "desired sort not supported for this type");
+                "type does not have compare function");
         return -1;
     }
 
+    switch (which) {
+        case NPY_QUICKSORT :
+            sort = npy_quicksort;
+            break;
+        case NPY_HEAPSORT :
+            sort = npy_heapsort;
+            break;
+        case NPY_MERGESORT :
+            sort = npy_mergesort;
+            break;
+        default:
+            PyErr_SetString(PyExc_TypeError,
+                    "requested sort kind is not supported");
+            goto fail;
+    }
+
+
     SWAPAXES2(op);
 
     ap = (PyArrayObject *)PyArray_FromAny((PyObject *)op,
@@ -1001,13 +1021,26 @@ PyArray_Sort(PyArrayObject *op, int axis, NPY_SORTKIND which)
     store_arr = global_obj;
     global_obj = ap;
     for (ip = PyArray_DATA(ap), i = 0; i < n; i++, ip += elsize*m) {
-        qsort(ip, m, elsize, qsortCompare);
+        res = sort(ip, m, elsize, sortCompare);
+        if (res < 0) {
+            break;
+        }
     }
     global_obj = store_arr;
 
     if (PyErr_Occurred()) {
         goto fail;
     }
+    else if (res == -NPY_ENOMEM) {
+        PyErr_NoMemory();
+        goto fail;
+    }
+    else if (res == -NPY_ECOMP) {
+        PyErr_SetString(PyExc_TypeError,
+                "sort comparison failed");
+        goto fail;
+    }
+
 
  finish:
     Py_DECREF(ap);  /* Should update op if needed */
@@ -1045,6 +1078,8 @@ PyArray_ArgSort(PyArrayObject *op, int axis, NPY_SORTKIND which)
     npy_intp i, j, n, m, orign;
     int argsort_elsize;
     char *store_ptr;
+    int res = 0;
+    int (*sort)(void *, size_t, size_t, npy_comparator);
 
     n = PyArray_NDIM(op);
     if ((n == 0) || (PyArray_SIZE(op) == 1)) {
@@ -1071,14 +1106,32 @@ PyArray_ArgSort(PyArrayObject *op, int axis, NPY_SORTKIND which)
         return (PyObject *)ret;
     }
 
-    if ((which != NPY_QUICKSORT) || PyArray_DESCR(op2)->f->compare == NULL) {
+    if (PyArray_DESCR(op2)->f->compare == NULL) {
         PyErr_SetString(PyExc_TypeError,
-                        "requested sort not available for type");
+                "type does not have compare function");
         Py_DECREF(op2);
         op = NULL;
         goto fail;
     }
 
+    switch (which) {
+        case NPY_QUICKSORT :
+            sort = npy_quicksort;
+            break;
+        case NPY_HEAPSORT :
+            sort = npy_heapsort;
+            break;
+        case NPY_MERGESORT :
+            sort = npy_mergesort;
+            break;
+        default:
+            PyErr_SetString(PyExc_TypeError,
+                    "requested sort kind is not supported");
+            Py_DECREF(op2);
+            op = NULL;
+            goto fail;
+    }
+
     /* ap will contain the reference to op2 */
     SWAPAXES(ap, op2);
     op = (PyArrayObject *)PyArray_ContiguousFromAny((PyObject *)ap,
@@ -1109,11 +1162,27 @@ PyArray_ArgSort(PyArrayObject *op, int axis, NPY_SORTKIND which)
         for (j = 0; j < m; j++) {
             ip[j] = j;
         }
-        qsort((char *)ip, m, sizeof(npy_intp), argsort_static_compare);
+        res = sort((char *)ip, m, sizeof(npy_intp), argsort_static_compare);
+        if (res < 0) {
+            break;
+        }
     }
     global_data = store_ptr;
     global_obj = store;
 
+    if (PyErr_Occurred()) {
+        goto fail;
+    }
+    else if (res == -NPY_ENOMEM) {
+        PyErr_NoMemory();
+        goto fail;
+    }
+    else if (res == -NPY_ECOMP) {
+        PyErr_SetString(PyExc_TypeError,
+                "sort comparison failed");
+        goto fail;
+    }
+
  finish:
     Py_DECREF(op);
     SWAPBACK(op, ret);
@@ -1123,7 +1192,6 @@ PyArray_ArgSort(PyArrayObject *op, int axis, NPY_SORTKIND which)
     Py_XDECREF(op);
     Py_XDECREF(ret);
     return NULL;
-
 }
 
 
diff --git a/numpy/core/src/npysort/heapsort.c.src b/numpy/core/src/npysort/heapsort.c.src
new file mode 100644
index 000000000..6a8c79f7a
--- /dev/null
+++ b/numpy/core/src/npysort/heapsort.c.src
@@ -0,0 +1,341 @@
+/* -*- c -*- */
+
+/*
+ * The purpose of this module is to add faster sort functions
+ * that are type-specific.  This is done by altering the
+ * function table for the builtin descriptors.
+ *
+ * These sorting functions are copied almost directly from numarray
+ * with a few modifications (complex comparisons compare the imaginary
+ * part if the real parts are equal, for example), and the names
+ * are changed.
+ *
+ * The original sorting code is due to Charles R. Harris who wrote
+ * it for numarray.
+ */
+
+/*
+ * Quick sort is usually the fastest, but the worst case scenario can
+ * be slower than the merge and heap sorts.  The merge sort requires
+ * extra memory and so for large arrays may not be useful.
+ *
+ * The merge sort is *stable*, meaning that equal components
+ * are unmoved from their entry versions, so it can be used to
+ * implement lexigraphic sorting on multiple keys.
+ *
+ * The heap sort is included for completeness.
+ */
+
+#define NPY_NO_DEPRECATED_API NPY_API_VERSION
+
+#include <stdlib.h>
+#include "npy_sort.h"
+#include "npysort_common.h"
+
+#define NOT_USED NPY_UNUSED(unused)
+#define PYA_QS_STACK 100
+#define SMALL_QUICKSORT 15
+#define SMALL_MERGESORT 20
+#define SMALL_STRING 16
+
+
+/*
+ *****************************************************************************
+ **                            NUMERIC SORTS                                **
+ *****************************************************************************
+ */
+
+
+/**begin repeat
+ *
+ * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
+ *         LONGLONG, ULONGLONG, HALF, FLOAT, DOUBLE, LONGDOUBLE,
+ *         CFLOAT, CDOUBLE, CLONGDOUBLE#
+ * #suff = bool, byte, ubyte, short, ushort, int, uint, long, ulong,
+ *         longlong, ulonglong, half, float, double, longdouble,
+ *         cfloat, cdouble, clongdouble#
+ * #type = npy_bool, npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int,
+ *         npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong,
+ *         npy_ushort, npy_float, npy_double, npy_longdouble, npy_cfloat,
+ *         npy_cdouble, npy_clongdouble#
+ */
+
+int
+heapsort_@suff@(@type@ *start, npy_intp n, void *NOT_USED)
+{
+    @type@ tmp, *a;
+    npy_intp i,j,l;
+
+    /* The array needs to be offset by one for heapsort indexing */
+    a = start - 1;
+
+    for (l = n>>1; l > 0; --l) {
+        tmp = a[l];
+        for (i = l, j = l<<1; j <= n;) {
+            if (j < n && @TYPE@_LT(a[j], a[j+1])) {
+                j += 1;
+            }
+            if (@TYPE@_LT(tmp, a[j])) {
+                a[i] = a[j];
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        a[i] = tmp;
+    }
+
+    for (; n > 1;) {
+        tmp = a[n];
+        a[n] = a[1];
+        n -= 1;
+        for (i = 1, j = 2; j <= n;) {
+            if (j < n && @TYPE@_LT(a[j], a[j+1])) {
+                j++;
+            }
+            if (@TYPE@_LT(tmp, a[j])) {
+                a[i] = a[j];
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        a[i] = tmp;
+    }
+
+    return 0;
+}
+
+
+int
+aheapsort_@suff@(@type@ *v, npy_intp *tosort, npy_intp n, void *NOT_USED)
+{
+    npy_intp *a, i,j,l, tmp;
+    /* The arrays need to be offset by one for heapsort indexing */
+    a = tosort - 1;
+
+    for (l = n>>1; l > 0; --l) {
+        tmp = a[l];
+        for (i = l, j = l<<1; j <= n;) {
+            if (j < n && @TYPE@_LT(v[a[j]], v[a[j+1]])) {
+                j += 1;
+            }
+            if (@TYPE@_LT(v[tmp], v[a[j]])) {
+                a[i] = a[j];
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        a[i] = tmp;
+    }
+
+    for (; n > 1;) {
+        tmp = a[n];
+        a[n] = a[1];
+        n -= 1;
+        for (i = 1, j = 2; j <= n;) {
+            if (j < n && @TYPE@_LT(v[a[j]], v[a[j+1]])) {
+                j++;
+            }
+            if (@TYPE@_LT(v[tmp], v[a[j]])) {
+                a[i] = a[j];
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        a[i] = tmp;
+    }
+
+    return 0;
+}
+
+/**end repeat**/
+
+
+/*
+ *****************************************************************************
+ **                             STRING SORTS                                **
+ *****************************************************************************
+ */
+
+
+/**begin repeat
+ *
+ * #TYPE = STRING, UNICODE#
+ * #suff = string, unicode#
+ * #type = npy_char, npy_ucs4#
+ */
+
+int
+heapsort_@suff@(@type@ *start, npy_intp n, PyArrayObject *arr)
+{
+    size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
+    @type@ *tmp = malloc(PyArray_ITEMSIZE(arr));
+    @type@ *a = start - len;
+    npy_intp i, j, l;
+
+    for (l = n>>1; l > 0; --l) {
+        @TYPE@_COPY(tmp, a + l*len, len);
+        for (i = l, j = l<<1; j <= n;) {
+            if (j < n && @TYPE@_LT(a + j*len, a + (j+1)*len, len))
+                j += 1;
+            if (@TYPE@_LT(tmp, a + j*len, len)) {
+                @TYPE@_COPY(a + i*len, a + j*len, len);
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        @TYPE@_COPY(a + i*len, tmp, len);
+    }
+
+    for (; n > 1;) {
+        @TYPE@_COPY(tmp, a + n*len, len);
+        @TYPE@_COPY(a + n*len, a + len, len);
+        n -= 1;
+        for (i = 1, j = 2; j <= n;) {
+            if (j < n && @TYPE@_LT(a + j*len, a + (j+1)*len, len))
+                j++;
+            if (@TYPE@_LT(tmp, a + j*len, len)) {
+                @TYPE@_COPY(a + i*len, a + j*len, len);
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        @TYPE@_COPY(a + i*len, tmp, len);
+    }
+
+    free(tmp);
+    return 0;
+}
+
+
+int
+aheapsort_@suff@(@type@ *v, npy_intp *tosort, npy_intp n, PyArrayObject *arr)
+{
+    size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
+    npy_intp *a, i,j,l, tmp;
+
+    /* The array needs to be offset by one for heapsort indexing */
+    a = tosort - 1;
+
+    for (l = n>>1; l > 0; --l) {
+        tmp = a[l];
+        for (i = l, j = l<<1; j <= n;) {
+            if (j < n && @TYPE@_LT(v + a[j]*len, v + a[j+1]*len, len))
+                j += 1;
+            if (@TYPE@_LT(v + tmp*len, v + a[j]*len, len)) {
+                a[i] = a[j];
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        a[i] = tmp;
+    }
+
+    for (; n > 1;) {
+        tmp = a[n];
+        a[n] = a[1];
+        n -= 1;
+        for (i = 1, j = 2; j <= n;) {
+            if (j < n && @TYPE@_LT(v + a[j]*len, v + a[j+1]*len, len))
+                j++;
+            if (@TYPE@_LT(v + tmp*len, v + a[j]*len, len)) {
+                a[i] = a[j];
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        a[i] = tmp;
+    }
+
+    return 0;
+}
+
+/**end repeat**/
+
+
+/*
+ *****************************************************************************
+ **                             GENERIC SORT                                **
+ *****************************************************************************
+ */
+
+
+/*
+ * This sort has almost the same signature as libc qsort and is intended to
+ * provide a heapsort for array types that don't have type specific sorts.
+ * The difference in the signature is an error return, as it might be the
+ * case that a memory allocation fails.
+ */
+int
+npy_heapsort(void *base, size_t num, size_t size, npy_comparator cmp)
+{
+    char *tmp = (char *) malloc(size);
+    char *a = (char *) base - size;
+    size_t i, j, l;
+
+    if (tmp == NULL) {
+        return -NPY_ENOMEM;
+    }
+
+    for (l = num >> 1; l > 0; --l) {
+        GENERIC_COPY(tmp, a + l*size, size);
+        for (i = l, j = l << 1; j <= num;) {
+            if (j < num && GENERIC_LT(a + j*size, a + (j+1)*size, cmp))
+                j += 1;
+            if (GENERIC_LT(tmp, a + j*size, cmp)) {
+                GENERIC_COPY(a + i*size, a + j*size, size);
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        GENERIC_COPY(a + i*size, tmp, size);
+    }
+
+    for (; num > 1;) {
+        GENERIC_COPY(tmp, a + num*size, size);
+        GENERIC_COPY(a + num*size, a + size, size);
+        num -= 1;
+        for (i = 1, j = 2; j <= num;) {
+            if (j < num && GENERIC_LT(a + j*size, a + (j+1)*size, cmp))
+                j++;
+            if (GENERIC_LT(tmp, a + j*size, cmp)) {
+                GENERIC_COPY(a + i*size, a + j*size, size);
+                i = j;
+                j += j;
+            }
+            else {
+                break;
+            }
+        }
+        GENERIC_COPY(a + i*size, tmp, size);
+    }
+
+    free(tmp);
+    return 0;
+}
diff --git a/numpy/core/src/npysort/mergesort.c.src b/numpy/core/src/npysort/mergesort.c.src
new file mode 100644
index 000000000..1003b95d1
--- /dev/null
+++ b/numpy/core/src/npysort/mergesort.c.src
@@ -0,0 +1,433 @@
+/* -*- c -*- */
+
+/*
+ * The purpose of this module is to add faster sort functions
+ * that are type-specific.  This is done by altering the
+ * function table for the builtin descriptors.
+ *
+ * These sorting functions are copied almost directly from numarray
+ * with a few modifications (complex comparisons compare the imaginary
+ * part if the real parts are equal, for example), and the names
+ * are changed.
+ *
+ * The original sorting code is due to Charles R. Harris who wrote
+ * it for numarray.
+ */
+
+/*
+ * Quick sort is usually the fastest, but the worst case scenario can
+ * be slower than the merge and heap sorts.  The merge sort requires
+ * extra memory and so for large arrays may not be useful.
+ *
+ * The merge sort is *stable*, meaning that equal components
+ * are unmoved from their entry versions, so it can be used to
+ * implement lexigraphic sorting on multiple keys.
+ *
+ * The heap sort is included for completeness.
+ */
+
+#define NPY_NO_DEPRECATED_API NPY_API_VERSION
+
+#include <stdlib.h>
+#include "npy_sort.h"
+#include "npysort_common.h"
+
+#define NOT_USED NPY_UNUSED(unused)
+#define PYA_QS_STACK 100
+#define SMALL_QUICKSORT 15
+#define SMALL_MERGESORT 20
+#define SMALL_STRING 16
+
+
+/*
+ *****************************************************************************
+ **                            NUMERIC SORTS                                **
+ *****************************************************************************
+ */
+
+
+/**begin repeat
+ *
+ * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
+ *         LONGLONG, ULONGLONG, HALF, FLOAT, DOUBLE, LONGDOUBLE,
+ *         CFLOAT, CDOUBLE, CLONGDOUBLE#
+ * #suff = bool, byte, ubyte, short, ushort, int, uint, long, ulong,
+ *         longlong, ulonglong, half, float, double, longdouble,
+ *         cfloat, cdouble, clongdouble#
+ * #type = npy_bool, npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int,
+ *         npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong,
+ *         npy_ushort, npy_float, npy_double, npy_longdouble, npy_cfloat,
+ *         npy_cdouble, npy_clongdouble#
+ */
+
+static void
+mergesort0_@suff@(@type@ *pl, @type@ *pr, @type@ *pw)
+{
+    @type@ vp, *pi, *pj, *pk, *pm;
+
+    if (pr - pl > SMALL_MERGESORT) {
+        /* merge sort */
+        pm = pl + ((pr - pl) >> 1);
+        mergesort0_@suff@(pl, pm, pw);
+        mergesort0_@suff@(pm, pr, pw);
+        for (pi = pw, pj = pl; pj < pm;) {
+            *pi++ = *pj++;
+        }
+        pi = pw + (pm - pl);
+        pj = pw;
+        pk = pl;
+        while (pj < pi && pm < pr) {
+            if (@TYPE@_LT(*pm, *pj)) {
+                *pk++ = *pm++;
+            }
+            else {
+                *pk++ = *pj++;
+            }
+        }
+        while(pj < pi) {
+            *pk++ = *pj++;
+        }
+    }
+    else {
+        /* insertion sort */
+        for (pi = pl + 1; pi < pr; ++pi) {
+            vp = *pi;
+            pj = pi;
+            pk = pi - 1;
+            while (pj > pl && @TYPE@_LT(vp, *pk)) {
+                *pj-- = *pk--;
+            }
+            *pj = vp;
+        }
+    }
+}
+
+
+int
+mergesort_@suff@(@type@ *start, npy_intp num, void *NOT_USED)
+{
+    @type@ *pl, *pr, *pw;
+
+    pl = start;
+    pr = pl + num;
+    pw = (@type@ *) malloc((num/2) * sizeof(@type@));
+    if (pw == NULL) {
+        PyErr_NoMemory();
+        return -NPY_ENOMEM;
+    }
+    mergesort0_@suff@(pl, pr, pw);
+
+    free(pw);
+    return 0;
+}
+
+
+static void
+amergesort0_@suff@(npy_intp *pl, npy_intp *pr, @type@ *v, npy_intp *pw)
+{
+    @type@ vp;
+    npy_intp vi, *pi, *pj, *pk, *pm;
+
+    if (pr - pl > SMALL_MERGESORT) {
+        /* merge sort */
+        pm = pl + ((pr - pl) >> 1);
+        amergesort0_@suff@(pl, pm, v, pw);
+        amergesort0_@suff@(pm, pr, v, pw);
+        for (pi = pw, pj = pl; pj < pm;) {
+            *pi++ = *pj++;
+        }
+        pi = pw + (pm - pl);
+        pj = pw;
+        pk = pl;
+        while (pj < pi && pm < pr) {
+            if (@TYPE@_LT(v[*pm], v[*pj])) {
+                *pk++ = *pm++;
+            }
+            else {
+                *pk++ = *pj++;
+            }
+        }
+        while(pj < pi) {
+            *pk++ = *pj++;
+        }
+    }
+    else {
+        /* insertion sort */
+        for (pi = pl + 1; pi < pr; ++pi) {
+            vi = *pi;
+            vp = v[vi];
+            pj = pi;
+            pk = pi - 1;
+            while (pj > pl && @TYPE@_LT(vp, v[*pk])) {
+                *pj-- = *pk--;
+            }
+            *pj = vi;
+        }
+    }
+}
+
+
+int
+amergesort_@suff@(@type@ *v, npy_intp *tosort, npy_intp num, void *NOT_USED)
+{
+    npy_intp *pl, *pr, *pw;
+
+    pl = tosort;
+    pr = pl + num;
+    pw = (npy_intp *) malloc((num/2) * sizeof(npy_intp));
+    if (pw == NULL) {
+        PyErr_NoMemory();
+        return -NPY_ENOMEM;
+    }
+    amergesort0_@suff@(pl, pr, v, pw);
+    free(pw);
+
+    return 0;
+}
+
+/**end repeat**/
+
+
+/*
+ *****************************************************************************
+ **                             STRING SORTS                                **
+ *****************************************************************************
+ */
+
+
+/**begin repeat
+ *
+ * #TYPE = STRING, UNICODE#
+ * #suff = string, unicode#
+ * #type = npy_char, npy_ucs4#
+ */
+
+static void
+mergesort0_@suff@(@type@ *pl, @type@ *pr, @type@ *pw, @type@ *vp, size_t len)
+{
+    @type@ *pi, *pj, *pk, *pm;
+
+    if ((size_t)(pr - pl) > SMALL_MERGESORT*len) {
+        /* merge sort */
+        pm = pl + (((pr - pl)/len) >> 1)*len;
+        mergesort0_@suff@(pl, pm, pw, vp, len);
+        mergesort0_@suff@(pm, pr, pw, vp, len);
+        @TYPE@_COPY(pw, pl, pm - pl);
+        pi = pw + (pm - pl);
+        pj = pw;
+        pk = pl;
+        while (pj < pi && pm < pr) {
+            if (@TYPE@_LT(pm, pj, len)) {
+                @TYPE@_COPY(pk, pm, len);
+                pm += len;
+                pk += len;
+            }
+            else {
+                @TYPE@_COPY(pk, pj, len);
+                pj += len;
+                pk += len;
+            }
+        }
+        @TYPE@_COPY(pk, pj, pi - pj);
+    }
+    else {
+        /* insertion sort */
+        for (pi = pl + len; pi < pr; pi += len) {
+            @TYPE@_COPY(vp, pi, len);
+            pj = pi;
+            pk = pi - len;
+            while (pj > pl && @TYPE@_LT(vp, pk, len)) {
+                @TYPE@_COPY(pj, pk, len);
+                pj -= len;
+                pk -= len;
+            }
+            @TYPE@_COPY(pj, vp, len);
+        }
+    }
+}
+
+
+int
+mergesort_@suff@(@type@ *start, npy_intp num, PyArrayObject *arr)
+{
+    size_t elsize = PyArray_ITEMSIZE(arr);
+    size_t len = elsize / sizeof(@type@);
+    @type@ *pl, *pr, *pw, *vp;
+    int err = 0;
+
+    pl = start;
+    pr = pl + num*len;
+    pw = (@type@ *) malloc((num/2) * elsize);
+    if (pw == NULL) {
+        PyErr_NoMemory();
+        err = -NPY_ENOMEM;
+        goto fail_0;
+    }
+    vp = (@type@ *) malloc(elsize);
+    if (vp == NULL) {
+        PyErr_NoMemory();
+        err = -NPY_ENOMEM;
+        goto fail_1;
+    }
+    mergesort0_@suff@(pl, pr, pw, vp, len);
+
+    free(vp);
+fail_1:
+    free(pw);
+fail_0:
+    return err;
+}
+
+
+static void
+amergesort0_@suff@(npy_intp *pl, npy_intp *pr, @type@ *v, npy_intp *pw, size_t len)
+{
+    @type@ *vp;
+    npy_intp vi, *pi, *pj, *pk, *pm;
+
+    if (pr - pl > SMALL_MERGESORT) {
+        /* merge sort */
+        pm = pl + ((pr - pl) >> 1);
+        amergesort0_@suff@(pl, pm, v, pw, len);
+        amergesort0_@suff@(pm, pr, v, pw, len);
+        for (pi = pw, pj = pl; pj < pm;) {
+            *pi++ = *pj++;
+        }
+        pi = pw + (pm - pl);
+        pj = pw;
+        pk = pl;
+        while (pj < pi && pm < pr) {
+            if (@TYPE@_LT(v + (*pm)*len, v + (*pj)*len, len)) {
+                *pk++ = *pm++;
+            }
+            else {
+                *pk++ = *pj++;
+            }
+        }
+        while (pj < pi) {
+            *pk++ = *pj++;
+        }
+    }
+    else {
+        /* insertion sort */
+        for (pi = pl + 1; pi < pr; ++pi) {
+            vi = *pi;
+            vp = v + vi*len;
+            pj = pi;
+            pk = pi - 1;
+            while (pj > pl && @TYPE@_LT(vp, v + (*pk)*len, len)) {
+                *pj-- = *pk--;
+            }
+            *pj = vi;
+        }
+    }
+}
+
+
+int
+amergesort_@suff@(@type@ *v, npy_intp *tosort, npy_intp num, PyArrayObject *arr)
+{
+    size_t elsize = PyArray_ITEMSIZE(arr);
+    size_t len = elsize / sizeof(@type@);
+    npy_intp *pl, *pr, *pw;
+
+    pl = tosort;
+    pr = pl + num;
+    pw = (npy_intp *) malloc((num/2) * sizeof(npy_intp));
+    if (pw == NULL) {
+        PyErr_NoMemory();
+        return -NPY_ENOMEM;
+    }
+    amergesort0_@suff@(pl, pr, v, pw, len);
+    free(pw);
+
+    return 0;
+}
+
+/**end repeat**/
+
+
+/*
+ *****************************************************************************
+ **                             GENERIC SORT                                **
+ *****************************************************************************
+ */
+
+
+static void
+npy_mergesort0(char *pl, char *pr, char *pw, char *vp, size_t size, npy_comparator cmp)
+{
+    char *pi, *pj, *pk, *pm;
+
+    if ((size_t)(pr - pl) > SMALL_MERGESORT*size) {
+        /* merge sort */
+        pm = pl + (((pr - pl)/size) >> 1)*size;
+        npy_mergesort0(pl, pm, pw, vp, size, cmp);
+        npy_mergesort0(pm, pr, pw, vp, size, cmp);
+        GENERIC_COPY(pw, pl, pm - pl);
+        pi = pw + (pm - pl);
+        pj = pw;
+        pk = pl;
+        while (pj < pi && pm < pr) {
+            if (GENERIC_LT(pm, pj, cmp)) {
+                GENERIC_COPY(pk, pm, size);
+                pm += size;
+                pk += size;
+            }
+            else {
+                GENERIC_COPY(pk, pj, size);
+                pj += size;
+                pk += size;
+            }
+        }
+        GENERIC_COPY(pk, pj, pi - pj);
+    }
+    else {
+        /* insertion sort */
+        for (pi = pl + size; pi < pr; pi += size) {
+            GENERIC_COPY(vp, pi, size);
+            pj = pi;
+            pk = pi - size;
+            while (pj > pl && GENERIC_LT(vp, pk, cmp)) {
+                GENERIC_COPY(pj, pk, size);
+                pj -= size;
+                pk -= size;
+            }
+            GENERIC_COPY(pj, vp, size);
+        }
+    }
+}
+
+
+/*
+ * This sort has almost the same signature as libc qsort and is intended to
+ * provide a mergesort for array types that don't have type specific sorts.
+ * The difference in the signature is an error return, as it might be the
+ * case that a memory allocation fails.
+ */
+int
+npy_mergesort(void *base, size_t num, size_t size, npy_comparator cmp)
+{
+    char *pl, *pr, *pw, *vp;
+    int err = 0;
+
+    pl = (char *)base;
+    pr = pl + num*size;
+    pw = (char *) malloc((num/2) * size);
+    if (pw == NULL) {
+        err = -NPY_ENOMEM;
+        goto fail_0;
+    }
+    vp = (char *) malloc(size);
+    if (vp == NULL) {
+        err = -NPY_ENOMEM;
+        goto fail_1;
+    }
+    npy_mergesort0(pl, pr, pw, vp, size, cmp);
+
+    free(vp);
+fail_1:
+    free(pw);
+fail_0:
+    return err;
+}
diff --git a/numpy/core/src/npysort/npysort_common.h b/numpy/core/src/npysort/npysort_common.h
index bc1d278b0..318a80222 100644
--- a/numpy/core/src/npysort/npysort_common.h
+++ b/numpy/core/src/npysort/npysort_common.h
@@ -1,15 +1,9 @@
 #ifndef __NPY_SORT_COMMON_H__
 #define __NPY_SORT_COMMON_H__
 
-/*
-#include "Python.h"
-#include "numpy/npy_common.h"
-#include "numpy/npy_math.h"
-#include "npy_config.h"
-*/
+#include <stdlib.h>
 #include <numpy/ndarraytypes.h>
 
-
 /*
  *****************************************************************************
  **                        SWAP MACROS                                      **
@@ -155,14 +149,17 @@ npy_half_lt_nonan(npy_half h1, npy_half h2)
     if (h1&0x8000u) {
         if (h2&0x8000u) {
             return (h1&0x7fffu) > (h2&0x7fffu);
-        } else {
+        }
+        else {
             /* Signed zeros are equal, have to check for it */
             return (h1 != 0x8000u) || (h2 != 0x0000u);
         }
-    } else {
+    }
+    else {
         if (h2&0x8000u) {
             return 0;
-        } else {
+        }
+        else {
             return (h1&0x7fffu) < (h2&0x7fffu);
         }
     }
@@ -176,7 +173,8 @@ HALF_LT(npy_half a, npy_half b)
 
     if (npy_half_isnan(b)) {
         ret = !npy_half_isnan(a);
-    } else {
+    }
+    else {
         ret = !npy_half_isnan(a) && npy_half_lt_nonan(a, b);
     }
 
@@ -254,14 +252,6 @@ CLONGDOUBLE_LT(npy_clongdouble a, npy_clongdouble b)
 }
 
 
-/* The PyObject functions are stubs for later use */
-NPY_INLINE static int
-PyObject_LT(PyObject *pa, PyObject *pb)
-{
-    return 0;
-}
-
-
 NPY_INLINE static void
 STRING_COPY(char *s1, char *s2, size_t len)
 {
@@ -333,4 +323,29 @@ UNICODE_LT(npy_ucs4 *s1, npy_ucs4 *s2, size_t len)
     return ret;
 }
 
+
+NPY_INLINE static void
+GENERIC_COPY(char *a, char *b, size_t len)
+{
+    memcpy(a, b, len);
+}
+
+
+NPY_INLINE static void
+GENERIC_SWAP(char *a, char *b, size_t len)
+{
+    while(len--) {
+        const char t = *a;
+        *a++ = *b;
+        *b++ = t;
+    }
+}
+
+
+NPY_INLINE static int
+GENERIC_LT(char *a, char *b, int (*cmp)(const void *, const void *))
+{
+    return cmp(a, b) < 0;
+}
+
 #endif
diff --git a/numpy/core/src/npysort/quicksort.c.src b/numpy/core/src/npysort/quicksort.c.src
new file mode 100644
index 000000000..2225887b1
--- /dev/null
+++ b/numpy/core/src/npysort/quicksort.c.src
@@ -0,0 +1,363 @@
+/* -*- c -*- */
+
+/*
+ * The purpose of this module is to add faster sort functions
+ * that are type-specific.  This is done by altering the
+ * function table for the builtin descriptors.
+ *
+ * These sorting functions are copied almost directly from numarray
+ * with a few modifications (complex comparisons compare the imaginary
+ * part if the real parts are equal, for example), and the names
+ * are changed.
+ *
+ * The original sorting code is due to Charles R. Harris who wrote
+ * it for numarray.
+ */
+
+/*
+ * Quick sort is usually the fastest, but the worst case scenario can
+ * be slower than the merge and heap sorts.  The merge sort requires
+ * extra memory and so for large arrays may not be useful.
+ *
+ * The merge sort is *stable*, meaning that equal components
+ * are unmoved from their entry versions, so it can be used to
+ * implement lexigraphic sorting on multiple keys.
+ *
+ * The heap sort is included for completeness.
+ */
+
+#define NPY_NO_DEPRECATED_API NPY_API_VERSION
+
+#include <stdlib.h>
+#include "npy_sort.h"
+#include "npysort_common.h"
+
+#define NOT_USED NPY_UNUSED(unused)
+#define PYA_QS_STACK 100
+#define SMALL_QUICKSORT 15
+#define SMALL_MERGESORT 20
+#define SMALL_STRING 16
+
+
+/*
+ *****************************************************************************
+ **                            NUMERIC SORTS                                **
+ *****************************************************************************
+ */
+
+
+/**begin repeat
+ *
+ * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
+ *         LONGLONG, ULONGLONG, HALF, FLOAT, DOUBLE, LONGDOUBLE,
+ *         CFLOAT, CDOUBLE, CLONGDOUBLE#
+ * #suff = bool, byte, ubyte, short, ushort, int, uint, long, ulong,
+ *         longlong, ulonglong, half, float, double, longdouble,
+ *         cfloat, cdouble, clongdouble#
+ * #type = npy_bool, npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int,
+ *         npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong,
+ *         npy_ushort, npy_float, npy_double, npy_longdouble, npy_cfloat,
+ *         npy_cdouble, npy_clongdouble#
+ */
+
+int
+quicksort_@suff@(@type@ *start, npy_intp num, void *NOT_USED)
+{
+    @type@ vp;
+    @type@ *pl = start;
+    @type@ *pr = start + num - 1;
+    @type@ *stack[PYA_QS_STACK];
+    @type@ **sptr = stack;
+    @type@ *pm, *pi, *pj, *pk;
+
+    for (;;) {
+        while ((pr - pl) > SMALL_QUICKSORT) {
+            /* quicksort partition */
+            pm = pl + ((pr - pl) >> 1);
+            if (@TYPE@_LT(*pm, *pl)) @TYPE@_SWAP(*pm, *pl);
+            if (@TYPE@_LT(*pr, *pm)) @TYPE@_SWAP(*pr, *pm);
+            if (@TYPE@_LT(*pm, *pl)) @TYPE@_SWAP(*pm, *pl);
+            vp = *pm;
+            pi = pl;
+            pj = pr - 1;
+            @TYPE@_SWAP(*pm, *pj);
+            for (;;) {
+                do ++pi; while (@TYPE@_LT(*pi, vp));
+                do --pj; while (@TYPE@_LT(vp, *pj));
+                if (pi >= pj) {
+                    break;
+                }
+                @TYPE@_SWAP(*pi,*pj);
+            }
+            pk = pr - 1;
+            @TYPE@_SWAP(*pi, *pk);
+            /* push largest partition on stack */
+            if (pi - pl < pr - pi) {
+                *sptr++ = pi + 1;
+                *sptr++ = pr;
+                pr = pi - 1;
+            }
+            else {
+                *sptr++ = pl;
+                *sptr++ = pi - 1;
+                pl = pi + 1;
+            }
+        }
+
+        /* insertion sort */
+        for (pi = pl + 1; pi <= pr; ++pi) {
+            vp = *pi;
+            pj = pi;
+            pk = pi - 1;
+            while (pj > pl && @TYPE@_LT(vp, *pk)) {
+                *pj-- = *pk--;
+            }
+            *pj = vp;
+        }
+        if (sptr == stack) {
+            break;
+        }
+        pr = *(--sptr);
+        pl = *(--sptr);
+    }
+
+    return 0;
+}
+
+
+int
+aquicksort_@suff@(@type@ *v, npy_intp* tosort, npy_intp num, void *NOT_USED)
+{
+    @type@ vp;
+    npy_intp *pl = tosort;
+    npy_intp *pr = tosort + num - 1;
+    npy_intp *stack[PYA_QS_STACK];
+    npy_intp **sptr = stack;
+    npy_intp *pm, *pi, *pj, *pk, vi;
+
+    for (;;) {
+        while ((pr - pl) > SMALL_QUICKSORT) {
+            /* quicksort partition */
+            pm = pl + ((pr - pl) >> 1);
+            if (@TYPE@_LT(v[*pm],v[*pl])) INTP_SWAP(*pm, *pl);
+            if (@TYPE@_LT(v[*pr],v[*pm])) INTP_SWAP(*pr, *pm);
+            if (@TYPE@_LT(v[*pm],v[*pl])) INTP_SWAP(*pm, *pl);
+            vp = v[*pm];
+            pi = pl;
+            pj = pr - 1;
+            INTP_SWAP(*pm, *pj);
+            for (;;) {
+                do ++pi; while (@TYPE@_LT(v[*pi], vp));
+                do --pj; while (@TYPE@_LT(vp, v[*pj]));
+                if (pi >= pj) {
+                    break;
+                }
+                INTP_SWAP(*pi, *pj);
+            }
+            pk = pr - 1;
+            INTP_SWAP(*pi,*pk);
+            /* push largest partition on stack */
+            if (pi - pl < pr - pi) {
+                *sptr++ = pi + 1;
+                *sptr++ = pr;
+                pr = pi - 1;
+            }
+            else {
+                *sptr++ = pl;
+                *sptr++ = pi - 1;
+                pl = pi + 1;
+            }
+        }
+
+        /* insertion sort */
+        for (pi = pl + 1; pi <= pr; ++pi) {
+            vi = *pi;
+            vp = v[vi];
+            pj = pi;
+            pk = pi - 1;
+            while (pj > pl && @TYPE@_LT(vp, v[*pk])) {
+                *pj-- = *pk--;
+            }
+            *pj = vi;
+        }
+        if (sptr == stack) {
+            break;
+        }
+        pr = *(--sptr);
+        pl = *(--sptr);
+    }
+
+    return 0;
+}
+
+/**end repeat**/
+
+
+/*
+ *****************************************************************************
+ **                             STRING SORTS                                **
+ *****************************************************************************
+ */
+
+
+/**begin repeat
+ *
+ * #TYPE = STRING, UNICODE#
+ * #suff = string, unicode#
+ * #type = npy_char, npy_ucs4#
+ */
+
+int
+quicksort_@suff@(@type@ *start, npy_intp num, PyArrayObject *arr)
+{
+    const size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
+    @type@ *vp = malloc(PyArray_ITEMSIZE(arr));
+    @type@ *pl = start;
+    @type@ *pr = start + (num - 1)*len;
+    @type@ *stack[PYA_QS_STACK], **sptr = stack, *pm, *pi, *pj, *pk;
+
+    for (;;) {
+        while ((size_t)(pr - pl) > SMALL_QUICKSORT*len) {
+            /* quicksort partition */
+            pm = pl + (((pr - pl)/len) >> 1)*len;
+            if (@TYPE@_LT(pm, pl, len)) @TYPE@_SWAP(pm, pl, len);
+            if (@TYPE@_LT(pr, pm, len)) @TYPE@_SWAP(pr, pm, len);
+            if (@TYPE@_LT(pm, pl, len)) @TYPE@_SWAP(pm, pl, len);
+            @TYPE@_COPY(vp, pm, len);
+            pi = pl;
+            pj = pr - len;
+            @TYPE@_SWAP(pm, pj, len);
+            for (;;) {
+                do pi += len; while (@TYPE@_LT(pi, vp, len));
+                do pj -= len; while (@TYPE@_LT(vp, pj, len));
+                if (pi >= pj) {
+                    break;
+                }
+                @TYPE@_SWAP(pi, pj, len);
+            }
+            pk = pr - len;
+            @TYPE@_SWAP(pi, pk, len);
+            /* push largest partition on stack */
+            if (pi - pl < pr - pi) {
+                *sptr++ = pi + len;
+                *sptr++ = pr;
+                pr = pi - len;
+            }
+            else {
+                *sptr++ = pl;
+                *sptr++ = pi - len;
+                pl = pi + len;
+            }
+        }
+
+        /* insertion sort */
+        for (pi = pl + len; pi <= pr; pi += len) {
+            @TYPE@_COPY(vp, pi, len);
+            pj = pi;
+            pk = pi - len;
+            while (pj > pl && @TYPE@_LT(vp, pk, len)) {
+                @TYPE@_COPY(pj, pk, len);
+                pj -= len;
+                pk -= len;
+            }
+            @TYPE@_COPY(pj, vp, len);
+        }
+        if (sptr == stack) {
+            break;
+        }
+        pr = *(--sptr);
+        pl = *(--sptr);
+    }
+
+    free(vp);
+    return 0;
+}
+
+
+int
+aquicksort_@suff@(@type@ *v, npy_intp* tosort, npy_intp num, PyArrayObject *arr)
+{
+    size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
+    @type@ *vp;
+    npy_intp *pl = tosort;
+    npy_intp *pr = tosort + num - 1;
+    npy_intp *stack[PYA_QS_STACK];
+    npy_intp **sptr=stack;
+    npy_intp *pm, *pi, *pj, *pk, vi;
+
+    for (;;) {
+        while ((pr - pl) > SMALL_QUICKSORT) {
+            /* quicksort partition */
+            pm = pl + ((pr - pl) >> 1);
+            if (@TYPE@_LT(v + (*pm)*len, v + (*pl)*len, len)) INTP_SWAP(*pm, *pl);
+            if (@TYPE@_LT(v + (*pr)*len, v + (*pm)*len, len)) INTP_SWAP(*pr, *pm);
+            if (@TYPE@_LT(v + (*pm)*len, v + (*pl)*len, len)) INTP_SWAP(*pm, *pl);
+            vp = v + (*pm)*len;
+            pi = pl;
+            pj = pr - 1;
+            INTP_SWAP(*pm,*pj);
+            for (;;) {
+                do ++pi; while (@TYPE@_LT(v + (*pi)*len, vp, len));
+                do --pj; while (@TYPE@_LT(vp, v + (*pj)*len, len));
+                if (pi >= pj) {
+                    break;
+                }
+                INTP_SWAP(*pi,*pj);
+            }
+            pk = pr - 1;
+            INTP_SWAP(*pi,*pk);
+            /* push largest partition on stack */
+            if (pi - pl < pr - pi) {
+                *sptr++ = pi + 1;
+                *sptr++ = pr;
+                pr = pi - 1;
+            }
+            else {
+                *sptr++ = pl;
+                *sptr++ = pi - 1;
+                pl = pi + 1;
+            }
+        }
+
+        /* insertion sort */
+        for (pi = pl + 1; pi <= pr; ++pi) {
+            vi = *pi;
+            vp = v + vi*len;
+            pj = pi;
+            pk = pi - 1;
+            while (pj > pl && @TYPE@_LT(vp, v + (*pk)*len, len)) {
+                *pj-- = *pk--;
+            }
+            *pj = vi;
+        }
+        if (sptr == stack) {
+            break;
+        }
+        pr = *(--sptr);
+        pl = *(--sptr);
+    }
+
+    return 0;
+}
+
+/**end repeat**/
+
+
+/*
+ *****************************************************************************
+ **                             GENERIC SORT                                **
+ *****************************************************************************
+ */
+
+
+/*
+ * This sort has almost the same signature as libc qsort and is intended to
+ * supply an error return for compatibility with the other generic sort
+ * kinds.
+ */
+int
+npy_quicksort(void *base, size_t num, size_t size, npy_comparator cmp)
+{
+    qsort(base, num, size, cmp);
+    return 0;
+}
diff --git a/numpy/core/src/npysort/sort.c.src b/numpy/core/src/npysort/sort.c.src
deleted file mode 100644
index 514131748..000000000
--- a/numpy/core/src/npysort/sort.c.src
+++ /dev/null
@@ -1,800 +0,0 @@
-/* -*- c -*- */
-
-/*
- * The purpose of this module is to add faster sort functions
- * that are type-specific.  This is done by altering the
- * function table for the builtin descriptors.
- *
- * These sorting functions are copied almost directly from numarray
- * with a few modifications (complex comparisons compare the imaginary
- * part if the real parts are equal, for example), and the names
- * are changed.
- *
- * The original sorting code is due to Charles R. Harris who wrote
- * it for numarray.
- */
-
-/*
- * Quick sort is usually the fastest, but the worst case scenario can
- * be slower than the merge and heap sorts.  The merge sort requires
- * extra memory and so for large arrays may not be useful.
- *
- * The merge sort is *stable*, meaning that equal components
- * are unmoved from their entry versions, so it can be used to
- * implement lexigraphic sorting on multiple keys.
- *
- * The heap sort is included for completeness.
- */
-
-#define NPY_NO_DEPRECATED_API NPY_API_VERSION
-
-#include <stdlib.h>
-#include "npy_sort.h"
-#include "npysort_common.h"
-
-#define NOT_USED NPY_UNUSED(unused)
-#define PYA_QS_STACK 100
-#define SMALL_QUICKSORT 15
-#define SMALL_MERGESORT 20
-#define SMALL_STRING 16
-
-/* These should be changed to PyDataMem_NEW/FREE if npysort is moved
- *  to the multiarray directory.
- */
-#define NpySortArray_malloc(size) ((char *)malloc(size))
-#define NpySortArray_free(ptr)  free(ptr)
-
-
-/*
- *****************************************************************************
- **                            NUMERIC SORTS                                **
- *****************************************************************************
- */
-
-
-/**begin repeat
- *
- * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
- *         LONGLONG, ULONGLONG, HALF, FLOAT, DOUBLE, LONGDOUBLE,
- *         CFLOAT, CDOUBLE, CLONGDOUBLE#
- * #suff = bool, byte, ubyte, short, ushort, int, uint, long, ulong,
- *         longlong, ulonglong, half, float, double, longdouble,
- *         cfloat, cdouble, clongdouble#
- * #type = npy_bool, npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int,
- *         npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong,
- *         npy_ushort, npy_float, npy_double, npy_longdouble, npy_cfloat,
- *         npy_cdouble, npy_clongdouble#
- */
-
-int
-quicksort_@suff@(@type@ *start, npy_intp num, void *NOT_USED)
-{
-    @type@ *pl = start;
-    @type@ *pr = start + num - 1;
-    @type@ vp;
-    @type@ *stack[PYA_QS_STACK], **sptr = stack, *pm, *pi, *pj, *pk;
-
-    for (;;) {
-        while ((pr - pl) > SMALL_QUICKSORT) {
-            /* quicksort partition */
-            pm = pl + ((pr - pl) >> 1);
-            if (@TYPE@_LT(*pm, *pl)) @TYPE@_SWAP(*pm, *pl);
-            if (@TYPE@_LT(*pr, *pm)) @TYPE@_SWAP(*pr, *pm);
-            if (@TYPE@_LT(*pm, *pl)) @TYPE@_SWAP(*pm, *pl);
-            vp = *pm;
-            pi = pl;
-            pj = pr - 1;
-            @TYPE@_SWAP(*pm, *pj);
-            for (;;) {
-                do ++pi; while (@TYPE@_LT(*pi, vp));
-                do --pj; while (@TYPE@_LT(vp, *pj));
-                if (pi >= pj) {
-                    break;
-                }
-                @TYPE@_SWAP(*pi,*pj);
-            }
-            pk = pr - 1;
-            @TYPE@_SWAP(*pi, *pk);
-            /* push largest partition on stack */
-            if (pi - pl < pr - pi) {
-                *sptr++ = pi + 1;
-                *sptr++ = pr;
-                pr = pi - 1;
-            }
-            else {
-                *sptr++ = pl;
-                *sptr++ = pi - 1;
-                pl = pi + 1;
-            }
-        }
-
-        /* insertion sort */
-        for (pi = pl + 1; pi <= pr; ++pi) {
-            vp = *pi;
-            pj = pi;
-            pk = pi - 1;
-            while (pj > pl && @TYPE@_LT(vp, *pk)) {
-                *pj-- = *pk--;
-            }
-            *pj = vp;
-        }
-        if (sptr == stack) {
-            break;
-        }
-        pr = *(--sptr);
-        pl = *(--sptr);
-    }
-
-    return 0;
-}
-
-int
-aquicksort_@suff@(@type@ *v, npy_intp* tosort, npy_intp num, void *NOT_USED)
-{
-    @type@ vp;
-    npy_intp *pl, *pr;
-    npy_intp *stack[PYA_QS_STACK], **sptr=stack, *pm, *pi, *pj, *pk, vi;
-
-    pl = tosort;
-    pr = tosort + num - 1;
-
-    for (;;) {
-        while ((pr - pl) > SMALL_QUICKSORT) {
-            /* quicksort partition */
-            pm = pl + ((pr - pl) >> 1);
-            if (@TYPE@_LT(v[*pm],v[*pl])) INTP_SWAP(*pm,*pl);
-            if (@TYPE@_LT(v[*pr],v[*pm])) INTP_SWAP(*pr,*pm);
-            if (@TYPE@_LT(v[*pm],v[*pl])) INTP_SWAP(*pm,*pl);
-            vp = v[*pm];
-            pi = pl;
-            pj = pr - 1;
-            INTP_SWAP(*pm,*pj);
-            for (;;) {
-                do ++pi; while (@TYPE@_LT(v[*pi],vp));
-                do --pj; while (@TYPE@_LT(vp,v[*pj]));
-                if (pi >= pj) {
-                    break;
-                }
-                INTP_SWAP(*pi,*pj);
-            }
-            pk = pr - 1;
-            INTP_SWAP(*pi,*pk);
-            /* push largest partition on stack */
-            if (pi - pl < pr - pi) {
-                *sptr++ = pi + 1;
-                *sptr++ = pr;
-                pr = pi - 1;
-            }
-            else {
-                *sptr++ = pl;
-                *sptr++ = pi - 1;
-                pl = pi + 1;
-            }
-        }
-
-        /* insertion sort */
-        for (pi = pl + 1; pi <= pr; ++pi) {
-            vi = *pi;
-            vp = v[vi];
-            pj = pi;
-            pk = pi - 1;
-            while (pj > pl && @TYPE@_LT(vp, v[*pk])) {
-                *pj-- = *pk--;
-            }
-            *pj = vi;
-        }
-        if (sptr == stack) {
-            break;
-        }
-        pr = *(--sptr);
-        pl = *(--sptr);
-    }
-
-    return 0;
-}
-
-
-int
-heapsort_@suff@(@type@ *start, npy_intp n, void *NOT_USED)
-{
-    @type@ tmp, *a;
-    npy_intp i,j,l;
-
-    /* The array needs to be offset by one for heapsort indexing */
-    a = start - 1;
-
-    for (l = n>>1; l > 0; --l) {
-        tmp = a[l];
-        for (i = l, j = l<<1; j <= n;) {
-            if (j < n && @TYPE@_LT(a[j], a[j+1])) {
-                j += 1;
-            }
-            if (@TYPE@_LT(tmp, a[j])) {
-                a[i] = a[j];
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        a[i] = tmp;
-    }
-
-    for (; n > 1;) {
-        tmp = a[n];
-        a[n] = a[1];
-        n -= 1;
-        for (i = 1, j = 2; j <= n;) {
-            if (j < n && @TYPE@_LT(a[j], a[j+1])) {
-                j++;
-            }
-            if (@TYPE@_LT(tmp, a[j])) {
-                a[i] = a[j];
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        a[i] = tmp;
-    }
-
-    return 0;
-}
-
-int
-aheapsort_@suff@(@type@ *v, npy_intp *tosort, npy_intp n, void *NOT_USED)
-{
-    npy_intp *a, i,j,l, tmp;
-    /* The arrays need to be offset by one for heapsort indexing */
-    a = tosort - 1;
-
-    for (l = n>>1; l > 0; --l) {
-        tmp = a[l];
-        for (i = l, j = l<<1; j <= n;) {
-            if (j < n && @TYPE@_LT(v[a[j]], v[a[j+1]])) {
-                j += 1;
-            }
-            if (@TYPE@_LT(v[tmp], v[a[j]])) {
-                a[i] = a[j];
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        a[i] = tmp;
-    }
-
-    for (; n > 1;) {
-        tmp = a[n];
-        a[n] = a[1];
-        n -= 1;
-        for (i = 1, j = 2; j <= n;) {
-            if (j < n && @TYPE@_LT(v[a[j]], v[a[j+1]])) {
-                j++;
-            }
-            if (@TYPE@_LT(v[tmp], v[a[j]])) {
-                a[i] = a[j];
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        a[i] = tmp;
-    }
-
-    return 0;
-}
-
-static void
-mergesort0_@suff@(@type@ *pl, @type@ *pr, @type@ *pw)
-{
-    @type@ vp, *pi, *pj, *pk, *pm;
-
-    if (pr - pl > SMALL_MERGESORT) {
-        /* merge sort */
-        pm = pl + ((pr - pl) >> 1);
-        mergesort0_@suff@(pl, pm, pw);
-        mergesort0_@suff@(pm, pr, pw);
-        for (pi = pw, pj = pl; pj < pm;) {
-            *pi++ = *pj++;
-        }
-        pj = pw;
-        pk = pl;
-        while (pj < pi && pm < pr) {
-            if (@TYPE@_LT(*pm, *pj)) {
-                *pk = *pm++;
-            }
-            else {
-                *pk = *pj++;
-            }
-            pk++;
-        }
-        while(pj < pi) {
-            *pk++ = *pj++;
-        }
-    }
-    else {
-        /* insertion sort */
-        for (pi = pl + 1; pi < pr; ++pi) {
-            vp = *pi;
-            pj = pi;
-            pk = pi -1;
-            while (pj > pl && @TYPE@_LT(vp, *pk)) {
-                *pj-- = *pk--;
-            }
-            *pj = vp;
-        }
-    }
-}
-
-int
-mergesort_@suff@(@type@ *start, npy_intp num, void *NOT_USED)
-{
-    @type@ *pl, *pr, *pw;
-
-    pl = start;
-    pr = pl + num;
-    pw = (@type@ *) NpySortArray_malloc((num/2)*sizeof(@type@));
-    if (!pw) {
-        PyErr_NoMemory();
-        return -1;
-    }
-    mergesort0_@suff@(pl, pr, pw);
-
-    NpySortArray_free(pw);
-    return 0;
-}
-
-static void
-amergesort0_@suff@(npy_intp *pl, npy_intp *pr, @type@ *v, npy_intp *pw)
-{
-    @type@ vp;
-    npy_intp vi, *pi, *pj, *pk, *pm;
-
-    if (pr - pl > SMALL_MERGESORT) {
-        /* merge sort */
-        pm = pl + ((pr - pl + 1)>>1);
-        amergesort0_@suff@(pl, pm-1, v, pw);
-        amergesort0_@suff@(pm, pr, v, pw);
-        for (pi = pw, pj = pl; pj < pm; ++pi, ++pj) {
-            *pi = *pj;
-        }
-        for (pk = pw, pm = pl; pk < pi && pj <= pr; ++pm) {
-            if (@TYPE@_LT(v[*pj], v[*pk])) {
-                *pm = *pj;
-                ++pj;
-            }
-            else {
-                *pm = *pk;
-                ++pk;
-            }
-        }
-        for (; pk < pi; ++pm, ++pk) {
-            *pm = *pk;
-        }
-    }
-    else {
-        /* insertion sort */
-        for (pi = pl + 1; pi <= pr; ++pi) {
-            vi = *pi;
-            vp = v[vi];
-            for (pj = pi, pk = pi - 1; pj > pl && @TYPE@_LT(vp, v[*pk]); --pj, --pk) {
-                *pj = *pk;
-            }
-            *pj = vi;
-        }
-    }
-}
-
-int
-amergesort_@suff@(@type@ *v, npy_intp *tosort, npy_intp num, void *NOT_USED)
-{
-    npy_intp *pl, *pr, *pw;
-
-    pl = tosort; pr = pl + num - 1;
-    pw = PyDimMem_NEW((1+num/2));
-
-    if (!pw) {
-        PyErr_NoMemory();
-        return -1;
-    }
-
-    amergesort0_@suff@(pl, pr, v, pw);
-    PyDimMem_FREE(pw);
-
-    return 0;
-}
-
-
-/**end repeat**/
-
-/*
- *****************************************************************************
- **                             STRING SORTS                                **
- *****************************************************************************
- */
-
-
-/**begin repeat
- *
- * #TYPE = STRING, UNICODE#
- * #suff = string, unicode#
- * #type = npy_char, npy_ucs4#
- */
-
-static void
-mergesort0_@suff@(@type@ *pl, @type@ *pr, @type@ *pw, @type@ *vp, size_t len)
-{
-    @type@ *pi, *pj, *pk, *pm;
-
-    if ((size_t)(pr - pl) > SMALL_MERGESORT*len) {
-        /* merge sort */
-        pm = pl + (((pr - pl)/len) >> 1)*len;
-        mergesort0_@suff@(pl, pm, pw, vp, len);
-        mergesort0_@suff@(pm, pr, pw, vp, len);
-        @TYPE@_COPY(pw, pl, pm - pl);
-        pi = pw + (pm - pl);
-        pj = pw;
-        pk = pl;
-        while (pj < pi && pm < pr) {
-            if (@TYPE@_LT(pm, pj, len)) {
-                @TYPE@_COPY(pk, pm, len);
-                pm += len;
-            }
-            else {
-                @TYPE@_COPY(pk, pj, len);
-                pj += len;
-            }
-            pk += len;
-        }
-        @TYPE@_COPY(pk, pj, pi - pj);
-    }
-    else {
-        /* insertion sort */
-        for (pi = pl + len; pi < pr; pi += len) {
-            @TYPE@_COPY(vp, pi, len);
-            pj = pi;
-            pk = pi - len;
-            while (pj > pl && @TYPE@_LT(vp, pk, len)) {
-                @TYPE@_COPY(pj, pk, len);
-                pj -= len;
-                pk -= len;
-            }
-            @TYPE@_COPY(pj, vp, len);
-        }
-    }
-}
-
-int
-mergesort_@suff@(@type@ *start, npy_intp num, PyArrayObject *arr)
-{
-    const size_t elsize = PyArray_ITEMSIZE(arr);
-    const size_t len = elsize / sizeof(@type@);
-    @type@ *pl, *pr, *pw, *vp;
-    int err = 0;
-
-    pl = start;
-    pr = pl + num*len;
-    pw = (@type@ *) NpySortArray_malloc((num/2)*elsize);
-    if (!pw) {
-        PyErr_NoMemory();
-        err = -1;
-        goto fail_0;
-    }
-    vp = (@type@ *) NpySortArray_malloc(elsize);
-    if (!vp) {
-        PyErr_NoMemory();
-        err = -1;
-        goto fail_1;
-    }
-    mergesort0_@suff@(pl, pr, pw, vp, len);
-
-    NpySortArray_free(vp);
-fail_1:
-    NpySortArray_free(pw);
-fail_0:
-    return err;
-}
-
-int
-quicksort_@suff@(@type@ *start, npy_intp num, PyArrayObject *arr)
-{
-    const size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
-    @type@ *vp = malloc(PyArray_ITEMSIZE(arr));
-    @type@ *pl = start;
-    @type@ *pr = start + (num - 1)*len;
-    @type@ *stack[PYA_QS_STACK], **sptr = stack, *pm, *pi, *pj, *pk;
-
-    for (;;) {
-        while ((size_t)(pr - pl) > SMALL_QUICKSORT*len) {
-            /* quicksort partition */
-            pm = pl + (((pr - pl)/len) >> 1)*len;
-            if (@TYPE@_LT(pm, pl, len)) @TYPE@_SWAP(pm, pl, len);
-            if (@TYPE@_LT(pr, pm, len)) @TYPE@_SWAP(pr, pm, len);
-            if (@TYPE@_LT(pm, pl, len)) @TYPE@_SWAP(pm, pl, len);
-            @TYPE@_COPY(vp, pm, len);
-            pi = pl;
-            pj = pr - len;
-            @TYPE@_SWAP(pm, pj, len);
-            for (;;) {
-                do pi += len; while (@TYPE@_LT(pi, vp, len));
-                do pj -= len; while (@TYPE@_LT(vp, pj, len));
-                if (pi >= pj) {
-                    break;
-                }
-                @TYPE@_SWAP(pi, pj, len);
-            }
-            pk = pr - len;
-            @TYPE@_SWAP(pi, pk, len);
-            /* push largest partition on stack */
-            if (pi - pl < pr - pi) {
-                *sptr++ = pi + len;
-                *sptr++ = pr;
-                pr = pi - len;
-            }
-            else {
-                *sptr++ = pl;
-                *sptr++ = pi - len;
-                pl = pi + len;
-            }
-        }
-
-        /* insertion sort */
-        for (pi = pl + len; pi <= pr; pi += len) {
-            @TYPE@_COPY(vp, pi, len);
-            pj = pi;
-            pk = pi - len;
-            while (pj > pl && @TYPE@_LT(vp, pk, len)) {
-                @TYPE@_COPY(pj, pk, len);
-                pj -= len;
-                pk -= len;
-            }
-            @TYPE@_COPY(pj, vp, len);
-        }
-        if (sptr == stack) {
-            break;
-        }
-        pr = *(--sptr);
-        pl = *(--sptr);
-    }
-
-    free(vp);
-    return 0;
-}
-
-
-int
-heapsort_@suff@(@type@ *start, npy_intp n, PyArrayObject *arr)
-{
-    size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
-    @type@ *tmp = malloc(PyArray_ITEMSIZE(arr));
-    @type@ *a = start - len;
-    npy_intp i,j,l;
-
-    for (l = n>>1; l > 0; --l) {
-        @TYPE@_COPY(tmp, a + l*len, len);
-        for (i = l, j = l<<1; j <= n;) {
-            if (j < n && @TYPE@_LT(a + j*len, a + (j+1)*len, len))
-                j += 1;
-            if (@TYPE@_LT(tmp, a + j*len, len)) {
-                @TYPE@_COPY(a + i*len, a + j*len, len);
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        @TYPE@_COPY(a + i*len, tmp, len);
-    }
-
-    for (; n > 1;) {
-        @TYPE@_COPY(tmp, a + n*len, len);
-        @TYPE@_COPY(a + n*len, a + len, len);
-        n -= 1;
-        for (i = 1, j = 2; j <= n;) {
-            if (j < n && @TYPE@_LT(a + j*len, a + (j+1)*len, len))
-                j++;
-            if (@TYPE@_LT(tmp, a + j*len, len)) {
-                @TYPE@_COPY(a + i*len, a + j*len, len);
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        @TYPE@_COPY(a + i*len, tmp, len);
-    }
-
-    free(tmp);
-    return 0;
-}
-
-
-int
-aheapsort_@suff@(@type@ *v, npy_intp *tosort, npy_intp n, PyArrayObject *arr)
-{
-    size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
-    npy_intp *a, i,j,l, tmp;
-
-    /* The array needs to be offset by one for heapsort indexing */
-    a = tosort - 1;
-
-    for (l = n>>1; l > 0; --l) {
-        tmp = a[l];
-        for (i = l, j = l<<1; j <= n;) {
-            if (j < n && @TYPE@_LT(v + a[j]*len, v + a[j+1]*len, len))
-                j += 1;
-            if (@TYPE@_LT(v + tmp*len, v + a[j]*len, len)) {
-                a[i] = a[j];
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        a[i] = tmp;
-    }
-
-    for (; n > 1;) {
-        tmp = a[n];
-        a[n] = a[1];
-        n -= 1;
-        for (i = 1, j = 2; j <= n;) {
-            if (j < n && @TYPE@_LT(v + a[j]*len, v + a[j+1]*len, len))
-                j++;
-            if (@TYPE@_LT(v + tmp*len, v + a[j]*len, len)) {
-                a[i] = a[j];
-                i = j;
-                j += j;
-            }
-            else {
-                break;
-            }
-        }
-        a[i] = tmp;
-    }
-
-    return 0;
-}
-
-
-int
-aquicksort_@suff@(@type@ *v, npy_intp* tosort, npy_intp num, PyArrayObject *arr)
-{
-    size_t len = PyArray_ITEMSIZE(arr)/sizeof(@type@);
-    @type@ *vp;
-    npy_intp *pl = tosort;
-    npy_intp *pr = tosort + num - 1;
-    npy_intp *stack[PYA_QS_STACK];
-    npy_intp **sptr=stack;
-    npy_intp *pm, *pi, *pj, *pk, vi;
-
-    for (;;) {
-        while ((pr - pl) > SMALL_QUICKSORT) {
-            /* quicksort partition */
-            pm = pl + ((pr - pl) >> 1);
-            if (@TYPE@_LT(v + (*pm)*len, v + (*pl)*len, len)) INTP_SWAP(*pm, *pl);
-            if (@TYPE@_LT(v + (*pr)*len, v + (*pm)*len, len)) INTP_SWAP(*pr, *pm);
-            if (@TYPE@_LT(v + (*pm)*len, v + (*pl)*len, len)) INTP_SWAP(*pm, *pl);
-            vp = v + (*pm)*len;
-            pi = pl;
-            pj = pr - 1;
-            INTP_SWAP(*pm,*pj);
-            for (;;) {
-                do ++pi; while (@TYPE@_LT(v + (*pi)*len, vp, len));
-                do --pj; while (@TYPE@_LT(vp, v + (*pj)*len, len));
-                if (pi >= pj) {
-                    break;
-                }
-                INTP_SWAP(*pi,*pj);
-            }
-            pk = pr - 1;
-            INTP_SWAP(*pi,*pk);
-            /* push largest partition on stack */
-            if (pi - pl < pr - pi) {
-                *sptr++ = pi + 1;
-                *sptr++ = pr;
-                pr = pi - 1;
-            }
-            else {
-                *sptr++ = pl;
-                *sptr++ = pi - 1;
-                pl = pi + 1;
-            }
-        }
-
-        /* insertion sort */
-        for (pi = pl + 1; pi <= pr; ++pi) {
-            vi = *pi;
-            vp = v + vi*len;
-            pj = pi;
-            pk = pi - 1;
-            while (pj > pl && @TYPE@_LT(vp, v + (*pk)*len, len)) {
-                *pj-- = *pk--;
-            }
-            *pj = vi;
-        }
-        if (sptr == stack) {
-            break;
-        }
-        pr = *(--sptr);
-        pl = *(--sptr);
-    }
-
-    return 0;
-}
-
-
-static void
-amergesort0_@suff@(npy_intp *pl, npy_intp *pr, @type@ *v, npy_intp *pw, int len)
-{
-    @type@ *vp;
-    npy_intp vi, *pi, *pj, *pk, *pm;
-
-    if (pr - pl > SMALL_MERGESORT) {
-        /* merge sort */
-        pm = pl + ((pr - pl) >> 1);
-        amergesort0_@suff@(pl, pm, v, pw, len);
-        amergesort0_@suff@(pm, pr, v, pw, len);
-        for (pi = pw, pj = pl; pj < pm;) {
-            *pi++ = *pj++;
-        }
-        pj = pw;
-        pk = pl;
-        while (pj < pi && pm < pr) {
-            if (@TYPE@_LT(v + (*pm)*len, v + (*pj)*len, len)) {
-                *pk = *pm++;
-            } else {
-                *pk = *pj++;
-            }
-            pk++;
-        }
-        while (pj < pi) {
-            *pk++ = *pj++;
-        }
-    } else {
-        /* insertion sort */
-        for (pi = pl + 1; pi < pr; ++pi) {
-            vi = *pi;
-            vp = v + vi*len;
-            pj = pi;
-            pk = pi -1;
-            while (pj > pl && @TYPE@_LT(vp, v + (*pk)*len, len)) {
-                *pj-- = *pk--;
-            }
-            *pj = vi;
-        }
-    }
-}
-
-
-int
-amergesort_@suff@(@type@ *v, npy_intp *tosort, npy_intp num, PyArrayObject *arr)
-{
-    const size_t elsize = PyArray_ITEMSIZE(arr);
-    const size_t len = elsize / sizeof(@type@);
-    npy_intp *pl, *pr, *pw;
-
-    pl = tosort;
-    pr = pl + num;
-    pw = PyDimMem_NEW(num/2);
-    if (!pw) {
-        PyErr_NoMemory();
-        return -1;
-    }
-    amergesort0_@suff@(pl, pr, v, pw, len);
-
-    PyDimMem_FREE(pw);
-    return 0;
-}
-/**end repeat**/
diff --git a/numpy/core/src/private/npy_sort.h b/numpy/core/src/private/npy_sort.h
index 94d831f86..120e7efc2 100644
--- a/numpy/core/src/private/npy_sort.h
+++ b/numpy/core/src/private/npy_sort.h
@@ -6,6 +6,10 @@
 #include <numpy/npy_common.h>
 #include <numpy/ndarraytypes.h>
 
+#define NPY_ENOMEM 1
+#define NPY_ECOMP 2
+
+typedef int (*npy_comparator)(const void *, const void *);
 
 int quicksort_bool(npy_bool *vec, npy_intp cnt, void *null);
 int heapsort_bool(npy_bool *vec, npy_intp cnt, void *null);
@@ -166,4 +170,9 @@ int aquicksort_unicode(npy_ucs4 *vec, npy_intp *ind, npy_intp cnt, PyArrayObject
 int aheapsort_unicode(npy_ucs4 *vec, npy_intp *ind, npy_intp cnt, PyArrayObject *arr);
 int amergesort_unicode(npy_ucs4 *vec, npy_intp *ind, npy_intp cnt, PyArrayObject *arr);
 
+
+int npy_quicksort(void *base, size_t num, size_t size, npy_comparator cmp);
+int npy_heapsort(void *base, size_t num, size_t size, npy_comparator cmp);
+int npy_mergesort(void *base, size_t num, size_t size, npy_comparator cmp);
+
 #endif
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index 73abc81b9..99861748c 100644
--- a/numpy/core/tests/test_multiarray.py
+++ b/numpy/core/tests/test_multiarray.py
@@ -491,7 +491,7 @@ class TestMethods(TestCase):
         # of sorted items must be greater than ~50 to check the actual
         # algorithm because quick and merge sort fall over to insertion
         # sort for small arrays.
-        a = np.arange(100)
+        a = np.arange(101)
         b = a[::-1].copy()
         for kind in ['q','m','h'] :
             msg = "scalar sort, kind=%s" % kind
@@ -527,7 +527,7 @@ class TestMethods(TestCase):
 
         # test string sorts.
         s = 'aaaaaaaa'
-        a = np.array([s + chr(i) for i in range(100)])
+        a = np.array([s + chr(i) for i in range(101)])
         b = a[::-1].copy()
         for kind in ['q', 'm', 'h'] :
             msg = "string sort, kind=%s" % kind
@@ -538,9 +538,9 @@ class TestMethods(TestCase):
             c.sort(kind=kind)
             assert_equal(c, a, msg)
 
-        # test unicode sort.
+        # test unicode sorts.
         s = 'aaaaaaaa'
-        a = np.array([s + chr(i) for i in range(100)], dtype=np.unicode)
+        a = np.array([s + chr(i) for i in range(101)], dtype=np.unicode)
         b = a[::-1].copy()
         for kind in ['q', 'm', 'h'] :
             msg = "unicode sort, kind=%s" % kind
@@ -551,7 +551,55 @@ class TestMethods(TestCase):
             c.sort(kind=kind)
             assert_equal(c, a, msg)
 
-        # todo, check object array sorts.
+        # test object array sorts.
+        a = np.empty((101,), dtype=np.object)
+        a[:] = range(101)
+        b = a[::-1]
+        for kind in ['q', 'h', 'm'] :
+            msg = "object sort, kind=%s" % kind
+            c = a.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
+            c = b.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
+
+        # test record array sorts.
+        dt = np.dtype([('f',float),('i',int)])
+        a = array([(i,i) for i in range(101)], dtype = dt)
+        b = a[::-1]
+        for kind in ['q', 'h', 'm'] :
+            msg = "object sort, kind=%s" % kind
+            c = a.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
+            c = b.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
+
+        # test datetime64 sorts.
+        a = np.arange(0, 101, dtype='datetime64[D]')
+        b = a[::-1]
+        for kind in ['q', 'h', 'm'] :
+            msg = "datetime64 sort, kind=%s" % kind
+            c = a.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
+            c = b.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
+
+        # test timedelta64 sorts.
+        a = np.arange(0, 101, dtype='timedelta64[D]')
+        b = a[::-1]
+        for kind in ['q', 'h', 'm'] :
+            msg = "timedelta64 sort, kind=%s" % kind
+            c = a.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
+            c = b.copy();
+            c.sort(kind=kind)
+            assert_equal(c, a, msg)
 
         # check axis handling. This should be the same for all type
         # specific sorts, so we only check it for one type and one kind
@@ -567,10 +615,6 @@ class TestMethods(TestCase):
         d = a.copy()
         d.sort()
         assert_equal(d, c, "test sort with default axis")
-        # using None is known fail at this point
-        # d = a.copy()
-        # d.sort(axis=None)
-        #assert_equal(d, c, "test sort with axis=None")
 
 
     def test_sort_order(self):
@@ -613,7 +657,7 @@ class TestMethods(TestCase):
         # of sorted items must be greater than ~50 to check the actual
         # algorithm because quick and merge sort fall over to insertion
         # sort for small arrays.
-        a = np.arange(100)
+        a = np.arange(101)
         b = a[::-1].copy()
         for kind in ['q','m','h'] :
             msg = "scalar argsort, kind=%s" % kind
@@ -637,10 +681,10 @@ class TestMethods(TestCase):
 
         # test string argsorts.
         s = 'aaaaaaaa'
-        a = np.array([s + chr(i) for i in range(100)])
+        a = np.array([s + chr(i) for i in range(101)])
         b = a[::-1].copy()
-        r = arange(100)
-        rr = r[::-1].copy()
+        r = np.arange(101)
+        rr = r[::-1]
         for kind in ['q', 'm', 'h'] :
             msg = "string argsort, kind=%s" % kind
             assert_equal(a.copy().argsort(kind=kind), r, msg)
@@ -648,16 +692,57 @@ class TestMethods(TestCase):
 
         # test unicode argsorts.
         s = 'aaaaaaaa'
-        a = np.array([s + chr(i) for i in range(100)], dtype=np.unicode)
-        b = a[::-1].copy()
-        r = arange(100)
-        rr = r[::-1].copy()
+        a = np.array([s + chr(i) for i in range(101)], dtype=np.unicode)
+        b = a[::-1]
+        r = np.arange(101)
+        rr = r[::-1]
         for kind in ['q', 'm', 'h'] :
             msg = "unicode argsort, kind=%s" % kind
             assert_equal(a.copy().argsort(kind=kind), r, msg)
             assert_equal(b.copy().argsort(kind=kind), rr, msg)
 
-        # todo, check object array argsorts.
+        # test object array argsorts.
+        a = np.empty((101,), dtype=np.object)
+        a[:] = range(101)
+        b = a[::-1]
+        r = np.arange(101)
+        rr = r[::-1]
+        for kind in ['q', 'm', 'h'] :
+            msg = "object argsort, kind=%s" % kind
+            assert_equal(a.copy().argsort(kind=kind), r, msg)
+            assert_equal(b.copy().argsort(kind=kind), rr, msg)
+
+        # test structured array argsorts.
+        dt = np.dtype([('f',float),('i',int)])
+        a = array([(i,i) for i in range(101)], dtype = dt)
+        b = a[::-1]
+        r = np.arange(101)
+        rr = r[::-1]
+        for kind in ['q', 'm', 'h'] :
+            msg = "structured array argsort, kind=%s" % kind
+            assert_equal(a.copy().argsort(kind=kind), r, msg)
+            assert_equal(b.copy().argsort(kind=kind), rr, msg)
+
+        # test datetime64 argsorts.
+        a = np.arange(0, 101, dtype='datetime64[D]')
+        b = a[::-1]
+        r = np.arange(101)
+        rr = r[::-1]
+        for kind in ['q', 'h', 'm'] :
+            msg = "datetime64 argsort, kind=%s" % kind
+            assert_equal(a.copy().argsort(kind=kind), r, msg)
+            assert_equal(b.copy().argsort(kind=kind), rr, msg)
+
+        # test timedelta64 argsorts.
+        a = np.arange(0, 101, dtype='timedelta64[D]')
+        b = a[::-1]
+        r = np.arange(101)
+        rr = r[::-1]
+        for kind in ['q', 'h', 'm'] :
+            msg = "timedelta64 argsort, kind=%s" % kind
+            assert_equal(a.copy().argsort(kind=kind), r, msg)
+            assert_equal(b.copy().argsort(kind=kind), rr, msg)
+
 
         # check axis handling. This should be the same for all type
         # specific argsorts, so we only check it for one type and one kind
@@ -933,7 +1018,7 @@ class TestMethods(TestCase):
         ro_diag.flags.writeable = False
         assert_equal(collect_warning_types(getattr, ro_diag,
                                            "__array_struct__"), [])
-        
+
 
     def test_ravel(self):
         a = np.array([[0,1],[2,3]])
@@ -1093,7 +1178,7 @@ class TestFancyIndexing(TestCase):
         x = xorig.copy()
         x[m3] = 10
         assert_array_equal(x, array([[1,10,3,4],[5,6,7,8]]))
-                           
+
 
 class TestStringCompare(TestCase):
     def test_string(self):
@@ -1701,7 +1786,7 @@ class TestFlat(TestCase):
         self.b = a[::2,::2]
         self.a0 = a0
         self.b0 = a0[::2,::2]
-        
+
     def test_contiguous(self):
         testpassed = False
         try:
@@ -1735,7 +1820,7 @@ class TestFlat(TestCase):
         assert d.flags.updateifcopy is False
         assert e.flags.updateifcopy is False
         assert f.flags.updateifcopy is True
-        assert f.base is self.b0    
+        assert f.base is self.b0
 
 class TestResize(TestCase):
     def test_basic(self):
diff --git a/numpy/lib/tests/test_arraysetops.py b/numpy/lib/tests/test_arraysetops.py
index e40c155a4..b0d2ca7c3 100644
--- a/numpy/lib/tests/test_arraysetops.py
+++ b/numpy/lib/tests/test_arraysetops.py
@@ -8,31 +8,62 @@ from numpy.lib.arraysetops import *
 
 import warnings
 
-class TestAso(TestCase):
-    def test_unique( self ):
-        a = np.array( [5, 7, 1, 2, 1, 5, 7] )
-
-        ec = np.array( [1, 2, 5, 7] )
-        c = unique( a )
-        assert_array_equal( c, ec )
-
-        vals, indices = unique( a, return_index=True )
+class TestSetOps(TestCase):
 
 
-        ed = np.array( [2, 3, 0, 1] )
-        assert_array_equal(vals, ec)
-        assert_array_equal(indices, ed)
-
-        vals, ind0, ind1 = unique( a, return_index=True,
-                                     return_inverse=True )
-
+    def test_unique( self ):
 
-        ee = np.array( [2, 3, 0, 1, 0, 2, 3] )
-        assert_array_equal(vals, ec)
-        assert_array_equal(ind0, ed)
-        assert_array_equal(ind1, ee)
+        def check_all(a, b, i1, i2, dt):
+            msg = "check values failed for type '%s'" % dt
+            v = unique(a)
+            assert_array_equal(v, b, msg)
+
+            msg = "check indexes failed for type '%s'" % dt
+            v, j = unique(a, 1, 0)
+            assert_array_equal(v, b, msg)
+            assert_array_equal(j, i1, msg)
+
+            msg = "check reverse indexes failed for type '%s'" % dt
+            v, j = unique(a, 0, 1)
+            assert_array_equal(v, b, msg)
+            assert_array_equal(j, i2, msg)
+
+            msg = "check with all indexes failed for type '%s'" % dt
+            v, j1, j2 = unique(a, 1, 1)
+            assert_array_equal(v, b, msg)
+            assert_array_equal(j1, i1, msg)
+            assert_array_equal(j2, i2, msg)
+
+        a = [5, 7, 1, 2, 1, 5, 7]*10
+        b = [1, 2, 5, 7]
+        i1 = [2, 3, 0, 1]
+        i2 = [2, 3, 0, 1, 0, 2, 3]*10
+
+        # test for numeric arrays
+        types = []
+        types.extend(np.typecodes['AllInteger'])
+        types.extend(np.typecodes['AllFloat'])
+        types.append('datetime64[D]')
+        types.append('timedelta64[D]')
+        for dt in types:
+            aa = np.array(a, dt)
+            bb = np.array(b, dt)
+            check_all(aa, bb, i1, i2, dt)
+
+        # test for object arrays
+        dt = 'O'
+        aa = np.empty(len(a), dt)
+        aa[:] = a
+        bb = np.empty(len(b), dt)
+        bb[:] = b
+        check_all(aa, bb, i1, i2, dt)
+
+        # test for structured arrays
+        dt = [('', 'i'), ('', 'i')]
+        aa = np.array(zip(a,a), dt)
+        bb = np.array(zip(b,b), dt)
+        check_all(aa, bb, i1, i2, dt)
 
-        assert_array_equal([], unique([]))
 
     def test_intersect1d( self ):
         # unique inputs