root/trunk/sort.d

/**A comprehensive sorting library for statistical functions.  Each function
 * takes N arguments, which are arrays or array-like objects, sorts the first
 * and sorts the rest in lockstep.  For merge and insertion sort, if the last
 * argument is a ulong*, increments the dereference of this ulong* by the bubble
 * sort distance between the first argument and the sorted version of the first
 * argument.  This is useful for some statistical calculations.
 *
 * All sorting functions have the precondition that all parallel input arrays
 * must have the same length.
 *
 * Notes:
 *
 * Comparison functions must be written such that compFun(x, x) == false.
 * For example, "a < b" is good, "a <= b" is not.
 *
 * These functions are heavily optimized for sorting arrays of
 * ints and floats (by far the most common case when doing statistical
 * calculations).  In these cases, they can be several times faster than the
 * equivalent functions in std.algorithm.  Since sorting is extremely important
 * for non-parametric statistics, this results in important real-world
 * performance gains.  However, it comes at a price in terms of generality:
 *
 * 1.  They assume that what they are sorting is cheap to copy via normal
 *     assignment.
 * 2.  They don't work at all with general ranges, only arrays and maybe
 *     ranges very similar to arrays.
 * 3.  All tuning and micro-optimization is done with ints and floats, not
 *     classes, large structs, strings, etc.
 *
 * Examples:
 * ---
 * auto foo = [3, 1, 2, 4, 5].dup;
 * auto bar = [8, 6, 7, 5, 3].dup;
 * qsort(foo, bar);
 * assert(foo == [1, 2, 3, 4, 5]);
 * assert(bar == [6, 7, 8, 5, 3]);
 * auto baz = [1.0, 0, -1, -2, -3].dup;
 * mergeSort!("a > b")(bar, foo, baz);
 * assert(bar == [8, 7, 6, 5, 3]);
 * assert(foo == [3, 2, 1, 4, 5]);
 * assert(baz == [-1.0, 0, 1, -2, -3]);
 * ---
 *
 * Author:  David Simcha
 */
 /*
 * License:
 * Boost Software License - Version 1.0 - August 17th, 2003
 *
 * Permission is hereby granted, free of charge, to any person or organization
 * obtaining a copy of the software and accompanying documentation covered by
 * this license (the "Software") to use, reproduce, display, distribute,
 * execute, and transmit the Software, and to prepare derivative works of the
 * Software, and to permit third-parties to whom the Software is furnished to
 * do so, all subject to the following:
 *
 * The copyright notices in the Software and this entire statement, including
 * the above license grant, this restriction and the following disclaimer,
 * must be included in all copies of the Software, in whole or in part, and
 * all derivative works of the Software, unless such copies or derivative
 * works are solely in the form of machine-executable object code generated by
 * a source language processor.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 */

module dstats.sort;

import std.traits, std.algorithm, std.math, std.functional, std.math, std.typecons,
       std.typetuple, std.range, std.array, std.traits, std.string : whitespace;

import dstats.alloc;

version(unittest) {
    import std.stdio, std.random;

    void main (){
    }
}

class SortException : Exception {
    this(string msg) {
        super(msg);
    }
}

/* CTFE function.  Used in isSimpleComparison.*/
/*private*/ string removeWhitespace(string input) pure nothrow {
    string ret;
    foreach(elem; input) {
        bool shouldAppend = true;
        foreach(whiteChar; whitespace) {
            if(elem == whiteChar) {
                shouldAppend = false;
                break;
            }
        }

        if(shouldAppend) {
            ret ~= elem;
        }
    }
    return ret;
}

/* Conservatively tests whether the comparison function is simple enough that
 * we can get away with comparing floats as if they were ints.
 */
/*private*/ template isSimpleComparison(alias comp) {
    static if(!isSomeString!(typeof(comp))) {
        enum bool isSimpleComparison = false;
    } else {
        enum bool isSimpleComparison =
            removeWhitespace(comp) == "a<b" ||
            removeWhitespace(comp) == "a>b";
    }
}

/*private*/ bool intIsNaN(I)(I i) {
    static if(is(I == int) || is(I == uint)) {
        // IEEE 754 single precision float has a 23-bit significand stored in the
        // lowest order bits, followed by an 8-bit exponent.  A NaN is when the
        // exponent bits are all ones and the significand is nonzero.
        enum uint significandMask = 0b111_1111_1111_1111_1111_1111UL;
        enum uint exponentMask = 0b1111_1111UL << 23;
    } else static if(is(I == long) || is(I == ulong)) {
        // IEEE 754 double precision float has a 52-bit significand stored in the
        // lowest order bits, followed by an 11-bit exponent.  A NaN is when the
        // exponent bits are all ones and the significand is nonzero.
        enum ulong significandMask =
            0b1111_1111_1111_1111_1111_1111_1111_1111_1111_1111_1111_1111_1111UL;
        enum ulong exponentMask = 0b111_1111_1111UL << 52;
    } else {
        static assert(0);
    }

    return ((i & exponentMask) == exponentMask) && ((i & significandMask) != 0);
}

unittest {
    // Test on randomly generated integers punned to floats.  We expect that
    // about 1 in 256 will be NaNs.
    foreach(i; 0..10_000) {
        uint randInt = uniform(0U, uint.max);
        assert(std.math.isNaN(*(cast(float*) &randInt)) == intIsNaN(randInt));
    }

    // Test on randomly generated integers punned to doubles.  We expect that
    // about 1 in 2048 will be NaNs.
    foreach(i; 0..1_000_000) {
        ulong randInt = (cast(ulong) uniform(0U, uint.max) << 32) + uniform(0U, uint.max);
        assert(std.math.isNaN(*(cast(double*) &randInt)) == intIsNaN(randInt));
    }
}

/*private*/ T prepareForSorting(alias comp, T)(T arr)
if(!isFloatingPoint!(ElementType!T)) {
    return arr;
}

/* Check for NaNs and throw an exception if they're present.*/
 real[] prepareForSorting(alias comp, F)(F arr)
if(is(F == real[])) {
    foreach(elem; arr) {
        if(isNaN(elem)) {
            throw new SortException("Can't sort NaNs.");
        }
    }

    return arr;
}

/* Check for NaN and do some bit twiddling so that a float or double can be
 * compared as an integer.  This results in approximately a 40% speedup
 * compared to just sorting as floats.
 */
 auto prepareForSorting(alias comp, F)(F arr)
if(is(F == double[]) || is(F == float[])) {
    static if(is(F == double[])) {
        alias long Int;
        enum signMask = 1UL << 63;
    } else {
        alias int Int;
        enum signMask = 1U << 31;
    }

    Int[] intArr = cast(Int[]) arr;
    foreach(i, ref elem; intArr) {
        if(intIsNaN(elem)) {
            // Roll back the bit twiddling in case someone catches the
            // exception, so that they don't see corrupted values.
            postProcess!comp(intArr[0..i]);

            throw new SortException("Can't sort NaNs.");
        }

        static if(isSimpleComparison!comp) {
            if(elem & signMask) {
                // Negative.
                elem ^= signMask;
                elem = ~elem;
            }
        }
    }

    static if(isSimpleComparison!comp) {
        return intArr;
    } else {
        return arr;
    }
}

/*private*/ void postProcess(alias comp, T)(T arr)
if(!isSimpleComparison!comp || (!is(T == double[]) && !is(T == float[]))) {}

/* Undo bit twiddling from prepareForSorting() to get back original
 * floating point numbers.
 */
/*private*/ void postProcess(alias comp, F)(F arr)
if((is(F == double[]) || is(F == float[])) && isSimpleComparison!comp) {
    static if(is(F == double[])) {
        alias long Int;
        enum mask = 1UL << 63;
    } else {
        alias int Int;
        enum mask = 1U << 31;
    }

    Int[] useMe = cast(Int[]) arr;
    foreach(ref elem; useMe) {
        if(elem & mask) {
            elem = ~elem;
            elem ^= mask;
        }
    }
}

version(unittest) {
    static void testFloating(alias fun, F)() {
        F[] testL = new F[1_000];
        foreach(ref e; testL) {
            e = uniform(-1_000_000, 1_000_000);
        }
        auto testL2 = testL.dup;

        static if(__traits(isSame, fun, mergeSortTemp)) {
            auto temp1 = testL.dup;
            auto temp2 = testL.dup;
        }

        foreach(i; 0..200) {
            randomShuffle(zip(testL, testL2));
            uint len = uniform(0, 1_000);

            static if(__traits(isSame, fun, mergeSortTemp)) {
                fun!"a > b"(testL[0..len], testL2[0..len], temp1[0..len], temp2[0..len]);
            } else {
                fun!("a > b")(testL[0..len], testL2[0..len]);
            }

            assert(isSorted!("a > b")(testL[0..len]));
            assert(testL == testL2, fun.stringof ~ '\t' ~ F.stringof);
        }
    }
}

void rotateLeft(T)(T input)
if(isRandomAccessRange!(T)) {
    if(input.length < 2) return;
    ElementType!(T) temp = input[0];
    foreach(i; 1..input.length) {
        input[i-1] = input[i];
    }
    input[$-1] = temp;
}

void rotateRight(T)(T input)
if(isRandomAccessRange!(T)) {
    if(input.length < 2) return;
    ElementType!(T) temp = input[$-1];
    for(size_t i = input.length - 1; i > 0; i--) {
        input[i] = input[i-1];
    }
    input[0] = temp;
}

/* Returns the index, NOT the value, of the median of the first, middle, last
 * elements of data.*/
size_t medianOf3(alias compFun, T)(T[] data) {
    alias binaryFun!(compFun) comp;
    immutable size_t mid = data.length / 2;
    immutable uint result = ((cast(uint) (comp(data[0], data[mid]))) << 2) |
                            ((cast(uint) (comp(data[0], data[$ - 1]))) << 1) |
                            (cast(uint) (comp(data[mid], data[$ - 1])));

    assert(result != 2 && result != 5 && result < 8); // Cases 2, 5 can't happen.
    switch(result) {
        case 1:  // 001
        case 6:  // 110
            return data.length - 1;
        case 3:  // 011
        case 4:  // 100
            return 0;
        case 0:  // 000
        case 7:  // 111
            return mid;
        default:
            assert(0);
    }
    assert(0);
}

unittest {
    assert(medianOf3!("a < b")([1,2,3,4,5]) == 2);
    assert(medianOf3!("a < b")([1,2,5,4,3]) == 4);
    assert(medianOf3!("a < b")([3,2,1,4,5]) == 0);
    assert(medianOf3!("a < b")([5,2,3,4,1]) == 2);
    assert(medianOf3!("a < b")([5,2,1,4,3]) == 4);
    assert(medianOf3!("a < b")([3,2,5,4,1]) == 0);
}


/**Quick sort.  Unstable, O(N log N) time average, worst
 * case, O(log N) space, small constant term in time complexity.
 *
 * In this implementation, the following steps are taken to avoid the
 * O(N<sup>2</sup>) worst case of naive quick sorts:
 *
 * 1.  At each recursion, the median of the first, middle and last elements of
 *     the array is used as the pivot.
 *
 * 2.  To handle the case of few unique elements, the "Fit Pivot" technique
 *     previously decribed by Andrei Alexandrescu is used.  This allows
 *     reasonable performance with few unique elements, with zero overhead
 *     in other cases.
 *
 * 3.  After a much larger than expected amount of recursion has occured,
 *     this function transitions to a heap sort.  This guarantees an O(N log N)
 *     worst case.*/
T[0] qsort(alias compFun = "a < b", T...)(T data)
in {
    assert(data.length > 0);
    size_t len = data[0].length;
    foreach(array; data[1..$]) {
        assert(array.length == len);
    }
} body {
    if(data[0].length < 25) {
        // Skip computing logarithm rather than waiting until qsortImpl to
        // do this.
        return insertionSort!compFun(data);
    }

    // Determines the transition point to a heap sort.
    uint TTL = cast(uint) (log2(cast(real) data[0].length) * 2);

    auto toSort = prepareForSorting!compFun(data[0]);

    /* qsort() throws if an invalid comparison function is passed.  Even in
     * this case, the data should be post-processed so the bit twiddling
     * hacks for floats can be undone.
     */
    try {
        qsortImpl!(compFun)(toSort, data[1..$], TTL);
    } finally {
        postProcess!compFun(data[0]);
    }

    return data[0];
}

//TTL = time to live, before transitioning to heap sort.
void qsortImpl(alias compFun, T...)(T data, uint TTL) {
    alias binaryFun!(compFun) comp;
    if(data[0].length < 25) {
         insertionSortImpl!(compFun)(data);
         return;
    }
    if(TTL == 0) {
        heapSortImpl!(compFun)(data);
        return;
    }
    TTL--;

    {
        immutable size_t med3 = medianOf3!(comp)(data[0]);
        foreach(array; data) {
            auto temp = array[med3];
            array[med3] = array[$ - 1];
            array[$ - 1] = temp;
        }
    }

    T less, greater;
    size_t lessI = size_t.max, greaterI = data[0].length - 1;

    auto pivot = data[0][$ - 1];
    if(comp(pivot, pivot)) {
        throw new SortException
            ("Comparison function must be such that compFun(x, x) == false.");
    }

    while(true) {
        while(comp(data[0][++lessI], pivot)) {}
        while(greaterI > 0 && comp(pivot, data[0][--greaterI])) {}

        if(lessI < greaterI) {
            foreach(array; data) {
                auto temp = array[lessI];
                array[lessI] = array[greaterI];
                array[greaterI] = temp;
            }
        } else break;
    }

    foreach(ti, array; data) {
        auto temp = array[$ - 1];
        array[$ - 1] = array[lessI];
        array[lessI] = temp;
        less[ti] = array[0..min(lessI, greaterI + 1)];
        greater[ti] = array[lessI + 1..$];
    }
    // Allow tail recursion optimization for larger block.  This guarantees
    // that, given a reasonable amount of stack space, no stack overflow will
    // occur even in pathological cases.
    if(greater[0].length > less[0].length) {
        qsortImpl!(compFun)(less, TTL);
        qsortImpl!(compFun)(greater, TTL);
        return;
    } else {
        qsortImpl!(compFun)(greater, TTL);
        qsortImpl!(compFun)(less, TTL);
    }
}

unittest {
    {  // Test integer.
        uint[] test = new uint[1_000];
        foreach(ref e; test) {
            e = uniform(0, 100);
        }
        auto test2 = test.dup;
        foreach(i; 0..1_000) {
            randomShuffle(zip(test, test2));
            uint len = uniform(0, 1_000);
            qsort(test[0..len], test2[0..len]);
            assert(isSorted(test[0..len]));
            assert(test == test2);
        }
    }

    testFloating!(qsort, float)();
    testFloating!(qsort, double)();
    testFloating!(qsort, real)();

    auto nanArr = [double.nan, 1.0];
    try {
        qsort(nanArr);
        assert(0);
    } catch(SortException) {}
}

/* Keeps track of what array merge sort data is in.  This is a speed hack to
 * copy back and forth less.*/
/*private*/ enum {
    DATA,
    TEMP
}

/**Merge sort.  O(N log N) time, O(N) space, small constant.  Stable sort.
 * If last argument is a ulong* instead of an array-like type,
 * the dereference of the ulong* will be incremented by the bubble sort
 * distance between the input array and the sorted version.  This is useful
 * in some statistics functions such as Kendall's tau.*/
T[0] mergeSort(alias compFun = "a < b", T...)(T data)
in {
    assert(data.length > 0);
    size_t len = data[0].length;
    foreach(array; data[1..$]) {
        static if(!is(typeof(array) == ulong*))
            assert(array.length == len);
    }
} body {
    if(data[0].length < 65) {  //Avoid mem allocation.
        return insertionSortImpl!(compFun)(data);
    }
    static if(is(T[$ - 1] == ulong*)) {
        enum dl = data.length - 1;
        alias data[$ - 1] swapCount;
    } else {
        enum dl = data.length;
        alias TypeTuple!() swapCount; // Place holder.
    }

    auto keyArr = prepareForSorting!compFun(data[0]);
    auto toSort = TypeTuple!(keyArr, data[1..dl]);

    auto stateCache = TempAlloc.getState;
    typeof(toSort) temp;
    foreach(i, array; temp) {
        temp[i] = newStack!(typeof(temp[i][0]))(data[i].length, stateCache);
    }

    uint res = mergeSortImpl!(compFun)(toSort, temp, swapCount);
    if(res == TEMP) {
        foreach(ti, array; temp) {
            toSort[ti][0..$] = temp[ti][0..$];
        }
    }

    foreach(array; temp) {
        TempAlloc.free(stateCache);
    }

    postProcess!compFun(data[0]);
    return data[0];
}

unittest {
    uint[] test = new uint[1_000], stability = new uint[1_000];
    uint[] temp1 = new uint[1_000], temp2 = new uint[1_000];
    foreach(ref e; test) {
        e = uniform(0, 100);  //Lots of ties.
    }
    foreach(i; 0..100) {
        ulong mergeCount = 0, bubbleCount = 0;
        foreach(j, ref e; stability) {
            e = cast(uint) j;
        }
        randomShuffle(test);
        uint len = uniform(0, 1_000);
        // Testing bubble sort distance against bubble sort,
        // since bubble sort distance computed by bubble sort
        // is straightforward, unlikely to contain any subtle bugs.
        bubbleSort(test[0..len].dup, &bubbleCount);
        if(i & 1)  // Test both temp and non-temp branches.
            mergeSort(test[0..len], stability[0..len], &mergeCount);
        else
            mergeSortTemp(test[0..len], stability[0..len], temp1[0..len],
                          temp2[0..len], &mergeCount);
        assert(bubbleCount == mergeCount);
        assert(isSorted(test[0..len]));
        foreach(j; 1..len) {
            if(test[j - 1] == test[j]) {
                assert(stability[j - 1] < stability[j]);
            }
        }
    }
    // Test without swapCounts.
    foreach(i; 0..1000) {
        foreach(j, ref e; stability) {
            e = cast(uint) j;
        }
        randomShuffle(test);
        uint len = uniform(0, 1_000);
        if(i & 1)  // Test both temp and non-temp branches.
            mergeSort(test[0..len], stability[0..len]);
        else
            mergeSortTemp(test[0..len], stability[0..len], temp1[0..len],
                          temp2[0..len]);
        assert(isSorted(test[0..len]));
        foreach(j; 1..len) {
            if(test[j - 1] == test[j]) {
                assert(stability[j - 1] < stability[j]);
            }
        }
    }

    testFloating!(mergeSort, float)();
    testFloating!(mergeSort, double)();
    testFloating!(mergeSort, real)();

    testFloating!(mergeSortTemp, float)();
    testFloating!(mergeSortTemp, double)();
    testFloating!(mergeSortTemp, real)();
}

/**Merge sort, allowing caller to provide a temp variable.  This allows
 * recycling instead of repeated allocations.  If D is data, T is temp,
 * and U is a ulong* for calculating bubble sort distance, this can be called
 * as mergeSortTemp(D, D, D, T, T, T, U) or mergeSortTemp(D, D, D, T, T, T)
 * where each D has a T of corresponding type.
 *
 * Examples:
 * ---
 * int[] foo = [3, 1, 2, 4, 5].dup;
 * int[] temp = new uint[5];
 * mergeSortTemp!("a < b")(foo, temp);
 * assert(foo == [1, 2, 3, 4, 5]); // The contents of temp will be undefined.
 * foo = [3, 1, 2, 4, 5].dup;
 * real bar = [3.14L, 15.9, 26.5, 35.8, 97.9];
 * real temp2 = new real[5];
 * mergeSortTemp(foo, bar, temp, temp2);
 * assert(foo == [1, 2, 3, 4, 5]);
 * assert(bar == [15.9L, 26.5, 3.14, 35.8, 97.9]);
 * // The contents of both temp and temp2 will be undefined.
 * ---
 */
T[0] mergeSortTemp(alias compFun = "a < b", T...)(T data)
in {
    assert(data.length > 0);
    size_t len = data[0].length;
    foreach(array; data[1..$]) {
        static if(!is(typeof(array) == ulong*))
            assert(array.length == len);
    }
} body {
    static if(is(T[$ - 1] == ulong*)) {
        enum dl = data.length - 1;
    } else {
        enum dl = data.length;
    }

    auto keyArr = prepareForSorting!compFun(data[0]);
    auto keyTemp = cast(typeof(keyArr)) data[dl / 2];
    auto toSort = TypeTuple!(
        keyArr,
        data[1..dl / 2],
        keyTemp,
        data[dl / 2 + 1..$]
    );

    uint res = mergeSortImpl!(compFun)(toSort);

    if(res == TEMP) {
        foreach(ti, array; toSort[0..$ / 2]) {
            toSort[ti][0..$] = toSort[ti + dl / 2][0..$];
        }
    }

    postProcess!compFun(data[0]);
    return data[0];
}

/*private*/ uint mergeSortImpl(alias compFun = "a < b", T...)(T dataIn) {
    static if(is(T[$ - 1] == ulong*)) {
        alias dataIn[$ - 1] swapCount;
        alias dataIn[0..dataIn.length / 2] data;
        alias dataIn[dataIn.length / 2..$ - 1] temp;
    } else {  // Make empty dummy tuple.
        alias TypeTuple!() swapCount;
        alias dataIn[0..dataIn.length / 2] data;
        alias dataIn[dataIn.length / 2..$] temp;
    }

    if(data[0].length < 50) {
        insertionSortImpl!(compFun)(data, swapCount);
        return DATA;
    }
    size_t half = data[0].length / 2;
    typeof(data) left, right, tempLeft, tempRight;
    foreach(ti, array; data) {
        left[ti] = array[0..half];
        right[ti] = array[half..$];
        tempLeft[ti] = temp[ti][0..half];
        tempRight[ti] = temp[ti][half..$];
    }

    /* Implementation note:  The lloc, rloc stuff is a hack to avoid constantly
     * copying data back and forth between the data and temp arrays.
     * Instad of copying every time, I keep track of which array the last merge
     * went into, and only copy at the end or if the two sides ended up in
     * different arrays.*/
    uint lloc = mergeSortImpl!(compFun)(left, tempLeft, swapCount);
    uint rloc = mergeSortImpl!(compFun)(right, tempRight, swapCount);
    if(lloc == DATA && rloc == TEMP) {
        foreach(ti, array; tempLeft) {
            array[] = left[ti][];
        }
        lloc = TEMP;
    } else if(lloc == TEMP && rloc == DATA) {
        foreach(ti, array; tempRight) {
            array[] = right[ti][];
        }
    }
    if(lloc == DATA) {
        merge!(compFun)(left, right, temp, swapCount);
        return TEMP;
    } else {
        merge!(compFun)(tempLeft, tempRight, data, swapCount);
        return DATA;
    }
}

/*private*/ void merge(alias compFun, T...)(T data) {
    alias binaryFun!(compFun) comp;

    static if(is(T[$ - 1] == ulong*)) {
        enum dl = data.length - 1;  //Length after removing swapCount;
        alias data[$ - 1] swapCount;
    } else {
        enum dl = data.length;
    }

    static assert(dl % 3 == 0);
    alias data[0..dl / 3] left;
    alias  data[dl / 3..dl * 2 / 3] right;
    alias data[dl * 2 / 3..dl] result;
    static assert(left.length == right.length && right.length == result.length);
    size_t i = 0, l = 0, r = 0;
    while(l < left[0].length && r < right[0].length) {
        if(comp(right[0][r], left[0][l])) {

            static if(is(T[$ - 1] == ulong*)) {
                *swapCount += left[0].length - l;
            }

            foreach(ti, array; result) {
                result[ti][i] = right[ti][r];
            }
            r++;
        } else {
            foreach(ti, array; result) {
                result[ti][i] = left[ti][l];
            }
            l++;
        }
        i++;
    }
    if(right[0].length > r) {
        foreach(ti, array; result) {
            result[ti][i..$] = right[ti][r..$];
        }
    } else {
        foreach(ti, array; result) {
            result[ti][i..$] = left[ti][l..$];
        }
    }
}

/**In-place merge sort, based on C++ STL's stable_sort().  O(N log<sup>2</sup> N)
 * time complexity, O(1) space complexity, stable.  Much slower than plain
 * old mergeSort(), so only use it if you really need the O(1) space.*/
T[0] mergeSortInPlace(alias compFun = "a < b", T...)(T data)
in {
    assert(data.length > 0);
    size_t len = data[0].length;
    foreach(array; data[1..$]) {
        assert(array.length == len);
    }
} body {
    auto toSort = prepareForSorting!compFun(data[0]);
    mergeSortInPlaceImpl!compFun(toSort, data[1..$]);
    postProcess!compFun(data[0]);
    return data[0];
}

/*private*/ T[0] mergeSortInPlaceImpl(alias compFun, T...)(T data) {
    if (data[0].length <= 100)
        return insertionSortImpl!(compFun)(data);

    T left, right;
    foreach(ti, array; data) {
        left[ti] = array[0..$ / 2];
        right[ti] = array[$ / 2..$];
    }

    mergeSortInPlace!(compFun, T)(right);
    mergeSortInPlace!(compFun, T)(left);
    mergeInPlace!(compFun)(data, data[0].length / 2);
    return data[0];
}

unittest {
    uint[] test = new uint[1_000], stability = new uint[1_000];
    foreach(ref e; test) {
        e = uniform(0, 100);  //Lots of ties.
    }
    uint[] test2 = test.dup;
    foreach(i; 0..1000) {
        foreach(j, ref e; stability) {
            e = cast(uint) j;
        }
        randomShuffle(zip(test, test2));
        uint len = uniform(0, 1_000);
        mergeSortInPlace(test[0..len], test2[0..len], stability[0..len]);
        assert(isSorted(test[0..len]));
        assert(test == test2);
        foreach(j; 1..len) {
            if(test[j - 1] == test[j]) {
                assert(stability[j - 1] < stability[j]);
            }
        }
    }

    testFloating!(mergeSortInPlace, float)();
    testFloating!(mergeSortInPlace, double)();
    testFloating!(mergeSortInPlace, real)();
}

// Loosely based on C++ STL's __merge_without_buffer().
/*private*/ void mergeInPlace(alias compFun = "a < b", T...)(T data, size_t middle) {
    static size_t largestLess(alias compFun, T)(T[] data, T value) {
        alias binaryFun!(compFun) comp;
        size_t len = data.length, first, last = data.length, half, middle;

        while (len > 0) {
            half = len / 2;
            middle = first + half;
            if (comp(data[middle], value)) {
                first = middle + 1;
                len = len - half - 1;
            } else
                len = half;
        }
        return first;
    }

    static size_t smallestGr(alias compFun, T)(T[] data, T value) {
        alias binaryFun!(compFun) comp;
        size_t len = data.length, first, last = data.length, half, middle;

        while (len > 0) {
            half = len / 2;
            middle = first + half;
            if (comp(value, data[middle]))
                len = half;
            else {
                first = middle + 1;
                len = len - half - 1;
            }
        }
        return first;
    }


    alias binaryFun!(compFun) comp;
    if (data[0].length < 2 || middle == 0 || middle == data[0].length)
        return;
    if (data[0].length == 2) {
        if(comp(data[0][1], data[0][0])) {
            foreach(array; data) {
                auto temp = array[0];
                array[0] = array[1];
                array[1] = temp;
            }
        }
        return;
    }

    size_t half1, half2, firstCut, secondCut;

    if (middle > data[0].length - middle) {
        half1 = middle / 2;
        auto pivot = data[0][half1];
        half2 = largestLess!(compFun)(data[0][middle..$], pivot);
    } else {
        half2 = (data[0].length - middle) / 2;
        auto pivot = data[0][half2 + middle];
        half1 = smallestGr!(compFun)(data[0][0..middle], pivot);
    }

    foreach(array; data) {
        bringToFront(array[half1..middle], array[middle..middle + half2]);
    }
    size_t newMiddle = half1 + half2;

    T left, right;
    foreach(ti, array; data) {
        left[ti] = array[0..newMiddle];
        right[ti] = array[newMiddle..$];
    }

    mergeInPlace!(compFun, T)(left, half1);
    mergeInPlace!(compFun, T)(right, half2 + middle - newMiddle);
}


/**Heap sort.  Unstable, O(N log N) time average and worst case, O(1) space,
 * large constant term in time complexity.*/
T[0] heapSort(alias compFun = "a < b", T...)(T data)
in {
    assert(data.length > 0);
    size_t len = data[0].length;
    foreach(array; data[1..$]) {
        assert(array.length == len);
    }
} body {
    auto toSort = prepareForSorting!compFun(data[0]);
    heapSortImpl!compFun(toSort, data[1..$]);
    postProcess!compFun(data[0]);
    return data[0];
}

/*private*/ T[0] heapSortImpl(alias compFun, T...)(T input) {
    // Heap sort has such a huge constant that insertion sort's faster for N <
    // 100 (for reals; even larger for smaller types).
    if(input[0].length <= 100) {
        return insertionSortImpl!(compFun)(input);
    }

    alias binaryFun!(compFun) comp;
    if(input[0].length < 2) return input[0];
    makeMultiHeap!(compFun)(input);
    for(size_t end = input[0].length - 1; end > 0; end--) {
        foreach(ti, ia; input) {
            auto temp = ia[end];
            ia[end] = ia[0];
            ia[0] = temp;
        }
        multiSiftDown!(compFun)(input, 0, end);
    }
    return input[0];
}

unittest {
    uint[] test = new uint[1_000];
    foreach(ref e; test) {
        e = uniform(0, 100_000);
    }
    auto test2 = test.dup;
    foreach(i; 0..1_000) {
        randomShuffle(zip(test, test2));
        uint len = uniform(0, 1_000);
        heapSort(test[0..len], test2[0..len]);
        assert(isSorted(test[0..len]));
        assert(test == test2);
    }

    testFloating!(heapSort, float)();
    testFloating!(heapSort, double)();
    testFloating!(heapSort, real)();
}

void makeMultiHeap(alias compFun = "a < b", T...)(T input) {
    if(input[0].length < 2)
        return;
    alias binaryFun!(compFun) comp;
    for(sizediff_t start = (input[0].length - 1) / 2; start >= 0; start--) {
        multiSiftDown!(compFun)(input, start, input[0].length);
    }
}

void multiSiftDown(alias compFun = "a < b", T...)
     (T input, size_t root, size_t end) {
    alias binaryFun!(compFun) comp;
    alias input[0] a;
    while(root * 2 + 1 < end) {
        size_t child = root * 2 + 1;
        if(child + 1 < end && comp(a[child], a[child + 1])) {
            child++;
        }
        if(comp(a[root], a[child])) {
            foreach(ia; input) {
                auto temp = ia[root];
                ia[root] = ia[child];
                ia[child] = temp;
            }
            root = child;
        }
        else return;
    }
}

/**Insertion sort.  O(N<sup>2</sup>) time worst, average case, O(1) space, VERY
 * small constant, which is why it's useful for sorting small subarrays in
 * divide and conquer algorithms.  If last argument is a ulong*, increments
 * the dereference of this argument by the bubble sort distance between the
 * input array and the sorted version of the input.*/
T[0] insertionSort(alias compFun = "a < b", T...)(T data)
in {
    assert(data.length > 0);
    size_t len = data[0].length;
    foreach(array; data[1..$]) {
        static if(!is(typeof(array) == ulong*))
            assert(array.length == len);
    }
} body {
    auto toSort = prepareForSorting!compFun(data[0]);
    insertionSortImpl!compFun(toSort, data[1..$]);
    postProcess!compFun(data[0]);
    return data[0];
}

private template IndexType(T) {
    alias typeof(T.init[0]) IndexType;
}

/*private*/ T[0] insertionSortImpl(alias compFun, T...)(T data) {
    alias binaryFun!(compFun) comp;
    static if(is(T[$ - 1] == ulong*)) {
        enum dl = data.length - 1;
        alias data[$ - 1] swapCount;
    } else {
        enum dl = data.length;
    }

    alias data[0] keyArray;
    if(keyArray.length < 2) {
        return keyArray;
    }

    // Yes, I measured this, caching this value is actually faster on DMD.
    immutable maxJ = keyArray.length - 1;
    for(size_t i = keyArray.length - 2; i != size_t.max; --i) {
        size_t j = i;

        Tuple!(staticMap!(IndexType, typeof(data[0..dl]))) temp = void;
        foreach(ti, Type; typeof(data[0..dl])) {
            static if(hasElaborateAssign!Type) {
                emplace(&(temp.field[ti]), data[ti][i]);
            } else {
                temp.field[ti] = data[ti][i];
            }
        }

        for(; j < maxJ && comp(keyArray[j + 1], temp.field[0]); ++j) {
            // It's faster to do all copying here than to call rotateLeft()
            // later, probably due to better ILP.
            foreach(array; data[0..dl]) {
                array[j] = array[j + 1];
            }
        }

        foreach(ti, Unused; typeof(temp.field)) {
            data[ti][j] = temp.field[ti];
        }

        static if(is(typeof(swapCount))) {
            *swapCount += (j - i);  //Increment swapCount variable.
        }
    }

    return keyArray;
}

unittest {
    uint[] test = new uint[100], stability = new uint[100];
    foreach(ref e; test) {
        e = uniform(0, 100);  //Lots of ties.
    }
    foreach(i; 0..1_000) {
        ulong insertCount = 0, bubbleCount = 0;
        foreach(j, ref e; stability) {
            e = cast(uint) j;
        }
        randomShuffle(test);
        uint len = uniform(0, 100);
        // Testing bubble sort distance against bubble sort,
        // since bubble sort distance computed by bubble sort
        // is straightforward, unlikely to contain any subtle bugs.
        bubbleSort(test[0..len].dup, &bubbleCount);
        insertionSort(test[0..len], stability[0..len], &insertCount);
        assert(bubbleCount == insertCount);
        assert(isSorted(test[0..len]));
        foreach(j; 1..len) {
            if(test[j - 1] == test[j]) {
                assert(stability[j - 1] < stability[j]);
            }
        }
    }
}

// Kept around only because it's easy to implement, and therefore good for
// testing more complex sort functions against.  Especially useful for bubble
// sort distance, since it's straightforward with a bubble sort, and not with
// a merge sort or insertion sort.
version(unittest) {
    T[0] bubbleSort(alias compFun = "a < b", T...)(T data) {
        alias binaryFun!(compFun) comp;
        static if(is(T[$ - 1] == ulong*))
            enum dl = data.length - 1;
        else enum dl = data.length;
        if(data[0].length < 2)
            return data[0];
        bool swapExecuted;
        foreach(i; 0..data[0].length) {
            swapExecuted = false;
            foreach(j; 1..data[0].length) {
                if(comp(data[0][j], data[0][j - 1])) {
                    swapExecuted = true;
                    static if(is(T[$ - 1] == ulong*))
                        (*(data[$-1]))++;
                    foreach(array; data[0..dl])
                        swap(array[j-1], array[j]);
                }
            }
            if(!swapExecuted) return data[0];
        }
        return data[0];
    }
}

unittest {
    //Sanity check for bubble sort distance.
    uint[] test = [4, 5, 3, 2, 1];
    ulong dist = 0;
    bubbleSort(test, &dist);
    assert(dist == 9);
    dist = 0;
    test = [6, 1, 2, 4, 5, 3];
    bubbleSort(test, &dist);
    assert(dist == 7);
}

/**Returns the kth largest/smallest element (depending on compFun, 0-indexed)
 * in the input array in O(N) time.  Allocates memory, does not modify input
 * array.*/
T quickSelect(alias compFun = "a < b", T)(T[] data, sizediff_t k) {
    auto dataDup = data.tempdup;
    scope(exit) TempAlloc.free;
    return partitionK!(compFun)(dataDup, k);
}

/**Partitions the input data according to compFun, such that position k contains
 * the kth largest/smallest element according to compFun.  For all elements e
 * with indices < k, !compFun(data[k], e) is guaranteed to be true.  For all
 * elements e with indices > k, !compFun(e, data[k]) is guaranteed to be true.
 * For example, if compFun is "a < b", all elements with indices < k will be
 * <= data[k], and all elements with indices larger than k will be >= k.
 * Reorders any additional input arrays in lockstep.
 *
 * Examples:
 * ---
 * auto foo = [3, 1, 5, 4, 2].dup;
 * auto secondSmallest = partitionK(foo, 1);
 * assert(secondSmallest == 2);
 * foreach(elem; foo[0..1]) {
 *     assert(elem <= foo[1]);
 * }
 * foreach(elem; foo[2..$]) {
 *     assert(elem >= foo[1]);
 * }
 * ---
 *
 * Returns:  The kth element of the array.
 */
ElementType!(T[0]) partitionK(alias compFun = "a < b", T...)(T data, ptrdiff_t k)
in {
    assert(data.length > 0);
    size_t len = data[0].length;
    foreach(array; data[1..$]) {
        assert(array.length == len);
    }
} body {
    // Don't use the float-to-int trick because it's actually slower here
    // because the main part of the algorithm is O(N), not O(N log N).
    return partitionKImpl!compFun(data, k);
}

/*private*/ ElementType!(T[0]) partitionKImpl(alias compFun, T...)(T data, ptrdiff_t k) {
    alias binaryFun!(compFun) comp;

    {
        immutable size_t med3 = medianOf3!(comp)(data[0]);
        foreach(array; data) {
            auto temp = array[med3];
            array[med3] = array[$ - 1];
            array[$ - 1] = temp;
        }
    }

    ptrdiff_t lessI = -1, greaterI = data[0].length - 1;
    auto pivot = data[0][$ - 1];
    while(true) {
        while(comp(data[0][++lessI], pivot)) {}
        while(greaterI > 0 && comp(pivot, data[0][--greaterI])) {}

        if(lessI < greaterI) {
            foreach(array; data) {
                auto temp = array[lessI];
                array[lessI] = array[greaterI];
                array[greaterI] = temp;
            }
        } else break;
    }
    foreach(array; data) {
        auto temp = array[lessI];
        array[lessI] = array[$ - 1];
        array[$ - 1] = temp;
    }

    if((greaterI < k && lessI >= k) || lessI == k) {
        return data[0][k];
    } else if(lessI < k) {
        foreach(ti, array; data) {
            data[ti] = array[lessI + 1..$];
        }
        return partitionK!(compFun, T)(data, k - lessI - 1);
    } else {
        foreach(ti, array; data) {
            data[ti] = array[0..min(greaterI + 1, lessI)];
        }
        return partitionK!(compFun, T)(data, k);
    }
}

template ArrayElemType(T : T[]) {
    alias T ArrayElemType;
}

unittest {
    enum n = 1000;
    uint[] test = new uint[n];
    uint[] test2 = new uint[n];
    uint[] lockstep = new uint[n];
    foreach(ref e; test) {
        e = uniform(0, 1000);
    }
    foreach(i; 0..1_000) {
        test2[] = test[];
        lockstep[] = test[];
        uint len = uniform(0, n - 1) + 1;
        qsort!("a > b")(test2[0..len]);
        int k = uniform(0, len);
        auto qsRes = partitionK!("a > b")(test[0..len], lockstep[0..len], k);
        assert(qsRes == test2[k]);
        foreach(elem; test[0..k]) {
            assert(elem >= test[k]);
        }
        foreach(elem; test[k + 1..len]) {
            assert(elem <= test[k]);
        }
        assert(test == lockstep);
    }
}

/**Given a set of data points entered through the put function, this output range
 * maintains the invariant that the top N according to compFun will be
 * contained in the data structure.  Uses a heap internally, O(log N) insertion
 * time.  Good for finding the largest/smallest N elements of a very large
 * dataset that cannot be sorted quickly in its entirety, and may not even fit
 * in memory. If less than N datapoints have been entered, all are contained in
 * the structure.
 *
 * Examples:
 * ---
 * Random gen;
 * gen.seed(unpredictableSeed);
 * uint[] nums = seq(0U, 100U);
 * auto less = TopN!(uint, "a < b")(10);
 * auto more = TopN!(uint, "a > b")(10);
 * randomShuffle(nums, gen);
 * foreach(n; nums) {
 *     less.put(n);
 *     more.put(n);
 * }
 *  assert(less.getSorted == [0U, 1,2,3,4,5,6,7,8,9]);
 *  assert(more.getSorted == [99U, 98, 97, 96, 95, 94, 93, 92, 91, 90]);
 *  ---
 */
struct TopN(T, alias compFun = "a > b") {
private:
    alias binaryFun!(compFun) comp;
    uint n;
    uint nAdded;

    T[] nodes;
public:
    /** The variable ntop controls how many elements are retained.*/
    this(uint ntop) {
        n = ntop;
        nodes = new T[n];
    }

    /** Insert an element into the topN struct.*/
    void put(T elem) {
        if(nAdded < n) {
            nodes[nAdded] = elem;
            if(nAdded == n - 1) {
                makeMultiHeap!(comp)(nodes);
            }
            nAdded++;
        } else if(nAdded >= n) {
             if(comp(elem, nodes[0])) {
                nodes[0] = elem;
                multiSiftDown!(comp)(nodes, 0, nodes.length);
            }
        }
    }

    /**Get the elements currently in the struct.  Returns a reference to
     * internal state, elements will be in an arbitrary order.  Cheap.*/
    T[] getElements() {
        return nodes[0..min(n, nAdded)];
    }

    /**Returns the elements sorted by compFun.  The array returned is a
     * duplicate of the input array.  Not cheap.*/
    T[] getSorted() {
        return qsort!(comp)(nodes[0..min(n, nAdded)].dup);
    }
}

unittest {
    alias TopN!(uint, "a < b") TopNLess;
    alias TopN!(uint, "a > b") TopNGreater;
    Random gen;
    gen.seed(unpredictableSeed);
    uint[] nums = new uint[100];
    foreach(i, ref n; nums) {
        n = cast(uint) i;
    }
    foreach(i; 0..100) {
        auto less = TopNLess(10);
        auto more = TopNGreater(10);
        randomShuffle(nums, gen);
        foreach(n; nums) {
            less.put(n);
            more.put(n);
        }
        assert(less.getSorted == [0U, 1,2,3,4,5,6,7,8,9]);
        assert(more.getSorted == [99U, 98, 97, 96, 95, 94, 93, 92, 91, 90]);
    }
    foreach(i; 0..100) {
        auto less = TopNLess(10);
        auto more = TopNGreater(10);
        randomShuffle(nums, gen);
        foreach(n; nums[0..5]) {
            less.put(n);
            more.put(n);
        }
        assert(less.getSorted == qsort!("a < b")(nums[0..5]));
        assert(more.getSorted == qsort!("a > b")(nums[0..5]));
    }
}

// Verify that there are no TempAlloc memory leaks anywhere in the code covered
// by the unittest.  This should always be the last unittest of the module.
unittest {
    auto TAState = TempAlloc.getState;
    assert(TAState.used == 0);
    assert(TAState.nblocks < 2);
}