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


#include "Python.h"
#include "numpy/noprefix.h"
#include "numpy/npy_math.h"

#include "npy_config.h"

#define NOT_USED NPY_UNUSED(unused)
#define PYA_QS_STACK 100
#define SMALL_QUICKSORT 15
#define SMALL_MERGESORT 20
#define SMALL_STRING 16

/*
 *****************************************************************************
 **                        SWAP MACROS                                      **
 *****************************************************************************
 */

/**begin repeat
 *
 * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
 *         LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE, CFLOAT,
 *         CDOUBLE,CLONGDOUBLE, INTP#
 * #type = npy_bool, npy_byte, npy_ubyte, npy_short, npy_ushort, npy_int,
 *         npy_uint, npy_long, npy_ulong, npy_longlong, npy_ulonglong,
 *         npy_float, npy_double, npy_longdouble, npy_cfloat, npy_cdouble,
 *         npy_clongdouble, npy_intp#
 */
#define @TYPE@_SWAP(a,b) {@type@ tmp = (b); (b)=(a); (a) = tmp;}

/**end repeat**/

/*
 *****************************************************************************
 **                        COMPARISON FUNCTIONS                             **
 *****************************************************************************
 */

/**begin repeat
 *
 * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
 *         LONGLONG, ULONGLONG#
 * #type = Bool, byte, ubyte, short, ushort, int, uint, long, ulong,
 *         longlong, ulonglong#
 */
NPY_INLINE static int
@TYPE@_LT(@type@ a, @type@ b)
{
    return a < b;
}
/**end repeat**/


/**begin repeat
 *
 * #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
 * #type = float, double, longdouble#
 */
NPY_INLINE static int
@TYPE@_LT(@type@ a, @type@ b)
{
    return a < b || (b != b && a == a);
}
/**end repeat**/


/*
 * For inline functions SUN recommends not using a return in the then part
 * of an if statement. It's a SUN compiler thing, so assign the return value
 * to a variable instead.
 */

/**begin repeat
 *
 * #TYPE = CFLOAT, CDOUBLE, CLONGDOUBLE#
 * #type = cfloat, cdouble, clongdouble#
 */
NPY_INLINE static int
@TYPE@_LT(@type@ a, @type@ b)
{
    int ret;

    if (a.real < b.real) {
        ret = a.imag == a.imag || b.imag != b.imag;
    }
    else if (a.real > b.real) {
        ret = b.imag != b.imag && a.imag == a.imag;
    }
    else if (a.real == b.real || (a.real != a.real && b.real != b.real)) {
        ret =  a.imag < b.imag || (b.imag != b.imag && a.imag == a.imag);
    }
    else {
        ret = b.real != b.real;
    }

    return ret;
}
/**end repeat**/


/* 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)
{
    memcpy(s1, s2, len);
}


NPY_INLINE static void
STRING_SWAP(char *s1, char *s2, size_t len)
{
    while(len--) {
        const char t = *s1;
        *s1++ = *s2;
        *s2++ = t;
    }
}


NPY_INLINE static int
STRING_LT(char *s1, char *s2, size_t len)
{
    const unsigned char *c1 = (unsigned char *)s1;
    const unsigned char *c2 = (unsigned char *)s2;
    size_t i;
    int ret = 0;

    for (i = 0; i < len; ++i) {
        if (c1[i] != c2[i]) {
            ret = c1[i] < c2[i];
            break;
        }
    }
    return ret;
}


NPY_INLINE static void
UNICODE_COPY(npy_ucs4 *s1, npy_ucs4 *s2, size_t len)
{
    while(len--) {
        *s1++ = *s2++;
    }
}


NPY_INLINE static void
UNICODE_SWAP(npy_ucs4 *s1, npy_ucs4 *s2, size_t len)
{
    while(len--) {
        const npy_ucs4 t = *s1;
        *s1++ = *s2;
        *s2++ = t;
    }
}


NPY_INLINE static int
UNICODE_LT(npy_ucs4 *s1, npy_ucs4 *s2, size_t len)
{
    size_t i;
    int ret = 0;

    for (i = 0; i < len; ++i) {
        if (s1[i] != s2[i]) {
            ret = s1[i] < s2[i];
            break;
        }
    }
    return ret;
}


/*
 *****************************************************************************
 **                            NUMERIC SORTS                                **
 *****************************************************************************
 */


/**begin repeat
 *
 * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
 *         LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE,
 *         CFLOAT, CDOUBLE, CLONGDOUBLE#
 * #type = Bool, byte, ubyte, short, ushort, int, uint, long, ulong,
 *         longlong, ulonglong, float, double, longdouble,
 *         cfloat, cdouble, clongdouble#
 */


static int
@TYPE@_quicksort(@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;
}

static int
@TYPE@_aquicksort(@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;
}


static int
@TYPE@_heapsort(@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;
}

static int
@TYPE@_aheapsort(@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
@TYPE@_mergesort0(@type@ *pl, @type@ *pr, @type@ *pw)
{
    @type@ vp, *pi, *pj, *pk, *pm;

    if (pr - pl > SMALL_MERGESORT) {
        /* merge sort */
        pm = pl + ((pr - pl) >> 1);
        @TYPE@_mergesort0(pl, pm, pw);
        @TYPE@_mergesort0(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;
        }
    }
}

static int
@TYPE@_mergesort(@type@ *start, npy_intp num, void *NOT_USED)
{
    @type@ *pl, *pr, *pw;

    pl = start;
    pr = pl + num;
    pw = (@type@ *) PyDataMem_NEW((num/2)*sizeof(@type@));
    if (!pw) {
        PyErr_NoMemory();
        return -1;
    }
    @TYPE@_mergesort0(pl, pr, pw);

    PyDataMem_FREE(pw);
    return 0;
}

static void
@TYPE@_amergesort0(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);
        @TYPE@_amergesort0(pl,pm-1,v,pw);
        @TYPE@_amergesort0(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;
        }
    }
}

static int
@TYPE@_amergesort(@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;
    }

    @TYPE@_amergesort0(pl, pr, v, pw);
    PyDimMem_FREE(pw);

    return 0;
}


/**end repeat**/

/*
 *****************************************************************************
 **                             STRING SORTS                                **
 *****************************************************************************
 */


/**begin repeat
 *
 * #TYPE = STRING, UNICODE#
 * #type = char, PyArray_UCS4#
 */

static void
@TYPE@_mergesort0(@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;
        @TYPE@_mergesort0(pl, pm, pw, vp, len);
        @TYPE@_mergesort0(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);
        }
    }
}

static int
@TYPE@_mergesort(@type@ *start, npy_intp num, PyArrayObject *arr)
{
    const size_t elsize = arr->descr->elsize;
    const size_t len = elsize / sizeof(@type@);
    @type@ *pl, *pr, *pw, *vp;
    int err = 0;

    pl = start;
    pr = pl + num*len;
    pw = (@type@ *) PyDataMem_NEW((num/2)*elsize);
    if (!pw) {
        PyErr_NoMemory();
        err = -1;
        goto fail_0;
    }
    vp = (@type@ *) PyDataMem_NEW(elsize);
    if (!vp) {
        PyErr_NoMemory();
        err = -1;
        goto fail_1;
    }
    @TYPE@_mergesort0(pl, pr, pw, vp, len);

    PyDataMem_FREE(vp);
fail_1:
    PyDataMem_FREE(pw);
fail_0:
    return err;
}

static int
@TYPE@_quicksort(@type@ *start, npy_intp num, PyArrayObject *arr)
{
    const size_t len = arr->descr->elsize/sizeof(@type@);
    @type@ *vp = malloc(arr->descr->elsize);
    @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;
}


static int
@TYPE@_heapsort(@type@ *start, npy_intp n, PyArrayObject *arr)
{
    size_t len = arr->descr->elsize/sizeof(@type@);
    @type@ *tmp = malloc(arr->descr->elsize);
    @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;
}


static int
@TYPE@_aheapsort(@type@ *v, npy_intp *tosort, npy_intp n, PyArrayObject *arr)
{
    size_t len = arr->descr->elsize/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;
}


static int
@TYPE@_aquicksort(@type@ *v, npy_intp* tosort, npy_intp num, PyArrayObject *arr)
{
    size_t len = arr->descr->elsize/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
@TYPE@_amergesort0(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);
        @TYPE@_amergesort0(pl,pm,v,pw,len);
        @TYPE@_amergesort0(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;
        }
    }
}


static int
@TYPE@_amergesort(@type@ *v, npy_intp *tosort, npy_intp num, PyArrayObject *arr)
{
    const size_t elsize = arr->descr->elsize;
    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;
    }
    @TYPE@_amergesort0(pl, pr, v, pw, len);

    PyDimMem_FREE(pw);
    return 0;
}
/**end repeat**/

static void
add_sortfuncs(void)
{
    PyArray_Descr *descr;

    /**begin repeat
     *
     * #TYPE = BOOL, BYTE, UBYTE, SHORT, USHORT, INT, UINT, LONG, ULONG,
     *         LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE,
     *         CFLOAT, CDOUBLE, CLONGDOUBLE, STRING, UNICODE#
     */
    descr = PyArray_DescrFromType(PyArray_@TYPE@);
    descr->f->sort[PyArray_QUICKSORT] =
        (PyArray_SortFunc *)@TYPE@_quicksort;
    descr->f->argsort[PyArray_QUICKSORT] =
        (PyArray_ArgSortFunc *)@TYPE@_aquicksort;
    descr->f->sort[PyArray_HEAPSORT] =
        (PyArray_SortFunc *)@TYPE@_heapsort;
    descr->f->argsort[PyArray_HEAPSORT] =
        (PyArray_ArgSortFunc *)@TYPE@_aheapsort;
    descr->f->sort[PyArray_MERGESORT] =
        (PyArray_SortFunc *)@TYPE@_mergesort;
    descr->f->argsort[PyArray_MERGESORT] =
        (PyArray_ArgSortFunc *)@TYPE@_amergesort;
    /**end repeat**/

}

static struct PyMethodDef methods[] = {
    {NULL, NULL, 0, NULL}
};


#if defined(NPY_PY3K)
static struct PyModuleDef moduledef = {
        PyModuleDef_HEAD_INIT,
        "_sort",
        NULL,
        -1,
        methods,
        NULL,
        NULL,
        NULL,
        NULL
};
#endif

/* Initialization function for the module */
#if defined(NPY_PY3K)
PyObject *PyInit__sort(void) {
    PyObject *m;
    m = PyModule_Create(&moduledef);
    if (!m) {
        return NULL;
    }
    import_array();
    add_sortfuncs();
    return m;
}
#else
PyMODINIT_FUNC
init_sort(void) {
    Py_InitModule("_sort", methods);

    import_array();
    add_sortfuncs();
}
#endif