// Written in the D programming language.

/*
 *  Copyright (C) 2008 by David Simcha.
 *
 *  This software is provided 'as-is', without any express or implied
 *  warranty. In no event will the authors be held liable for any damages
 *  arising from the use of this software.
 *
 *  Permission is granted to anyone to use this software for any purpose,
 *  including commercial applications, and to alter it and redistribute it
 *  freely, subject to the following restrictions:
 *
 *  o  The origin of this software must not be misrepresented; you must not
 *     claim that you wrote the original software. If you use this software
 *     in a product, an acknowledgment in the product documentation would be
 *     appreciated but is not required.
 *  o  Altered source versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  o  This notice may not be removed or altered from any source
 *     distribution.
 */

import std.traits, std.algorithm;

version(unittest) {
    import std.stdio, std.random, std.math;
    void main () {}
}

/**Returns a new T[], skipping initialization.  Useful when fast memory
 * allocation is necessary, and the initial values are definitely
 * not going to be used.*/

T[] newVoid(T)(size_t length) {
    T* ptr = cast(T*) std.gc.malloc(length * T.sizeof);
    return ptr[0..length];
}

///Helper function for stableSelectionSort.
void rotateRight(T)(T[] input) {
    if(input.length < 2) return;
    T temp = input[$-1];
    for(size_t i = input.length - 1; i > 0; i--) {
        input[i] = input[i-1];
    }
    input[0] = temp;
}

/**Unstable quick multisort.  Sorts N arrays by first array.
 * Faster than stable version, and sorts in place. */

T[0] qsort(T...)(T data) {  //Faster than stable version.
    if(data[0].length < 8) {
         selectionSort(data);
         return data[0];
    }
    size_t middle = data[0].length / 2;
    foreach(array; data) swap(array[middle], array[$-1]);  //Avoid worst case if already sorted.
    T less, greater;
    size_t swapIndex = 0;
    foreach(i; 0..data[0].length - 1) {
        if(data[0][i] < data[0][$-1]) {
            foreach(array; data) swap(array[i], array[swapIndex]);
            swapIndex++;
        }
    }
    foreach(i, array; data) {
        swap(array[$-1], array[swapIndex]);
        less[i] = array[0..swapIndex];
        greater[i] = array[swapIndex + 1..$];
    }
    qsort(less);
    qsort(greater);
    return data[0];  //Make this function like "normal" qsort if data.length is one.
}

/**Stable version of qsort() function.  Does not sort in place, and
 * is about 15% slower than unstable version, so use only when necessary.
 * Also, performance is O(N^2) on pre-sorted arrays, because choosing an
 * intelligent pivot is difficult or impossible in a stable qsort.
 * stableQsort creates temp variables and forwards all data to
 * stableQsortBackEnd.  However, stableQsortBackEnd is callable directly,
 * in case the user wants to recycle pre-allocated temp space. */

T[0] stableQsort(T...)(T data) {
    typeof(data) temp;
    foreach(i, array; temp) {
        temp[i] = newVoid!(typeof(data[i][0]))(data[i].length);
    }
    stableQsortBackEnd(data, temp);
    foreach(i, array; temp) {
        delete temp[i];  //Do ths manually since GC sometimmes fails on large arrays.
    }
    return data[0];
}

T[0] stableQsortBackEnd(T...)(T dataIn) {  //Assumes first half of input is data, second half is temp.
    auto data = dataIn[0..dataIn.length/2];
    auto temp = dataIn[dataIn.length/2..$];
    assert(data[0].ptr != temp[0].ptr);
    if(data[0].length < 8) {
         stableSelectionSort(data);
         return data[0];
    }
    typeof(data) less, greater, lessTemp, greaterTemp;
    size_t index = 0;
    foreach(i; 0..data[0].length - 1) {
        if(data[0][i] <= data[0][$-1]) {
            foreach(ti, array; data) temp[ti][index] = array[i];
            index++;
        }
    }
    foreach(ti, array; data) temp[ti][index] = array[$-1];
    size_t mid = index;
    index++;
    foreach(i; 0..data[0].length - 1) {
        if(data[0][i] > data[0][$-1]) {
            foreach(ti, array; data) temp[ti][index] = array[i];
            index++;
        }
    }
    foreach(i, array; temp) {
        data[i][0..$] = array[0..$];
    }
    foreach(i, array; data) {
        less[i] = array[0..mid];
        greater[i] = array[mid + 1..$];
        lessTemp[i] = temp[i][0..mid];
        greaterTemp[i] = temp[i][mid + 1..$];
    }
    stableQsortBackEnd(less, lessTemp);
    stableQsortBackEnd(greater, greaterTemp);
    return data[0];
}

/**Merge multisort.  Takes in N arrays and sorts them by the
 * first array.  If last argument is a ulong*, increments this by
 * the number of swaps that would have been necessary to bubble sort
 * the data.  This is useful for calculating statistical functions like
 * Kendall correlation.  Slower than stableQsort on average, use only
 * if you need to count the swaps for some statistical function,
 * or if you need O(N log N) worst case performance. Like stableQsort(),
 * front end function creates temp variables and forwards them to back
 * end function.*/

T[0] mergeSort(T...)(T data) {
    static if(!isArray!(typeof(data[$-1]))) {  //This static if block accounts for the possibility of a swapCount pointer being passed.
        invariant uint dl = data.length - 1;
        ulong* swapCount = data[$-1];
    } else {
        invariant uint dl = data.length;
        ulong dummy = void;  //This is a dummy variable.  It doesn't need to be initialized.
        ulong* swapCount = &dummy;  //Makes things a lot easier if I use a dummy swapcount variable instaed of writing a whole bunch of static ifs.
    }
    typeof(data[0..dl]) temp;
    foreach(i, array; temp) temp[i] = newVoid!(typeof(data[i][0]))(data[i].length);
    T[0] output = mergeSortBackEnd(data[0..dl], temp, swapCount);
    foreach(ti, array; temp) delete temp[ti];
    return output;
}

T[0] mergeSortBackEnd(T...)(T dataIn) {
    auto data = dataIn[0..dataIn.length/2];
    auto temp = dataIn[dataIn.length/2..$-1];
    ulong* swapCount = dataIn[$-1];
    if(data[0].length < 2) {
        return stableSelectionSort(data, swapCount);
    }
    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..$];
    }
    mergeSortBackEnd(left, tempLeft, swapCount);
    mergeSortBackEnd(right, tempRight, swapCount);
    merge(left, right, temp, swapCount);
    foreach(ti, array; temp) {
        data[ti][0..$] = temp[ti][0..$];
    }
    return data[0];
}

///Merge function for merge sort.
void merge(T...)(T data) {
    ulong* swapCount = data[$-1];
    invariant uint dl = data.length - 1;  //Real length after removing swapCount;
    static assert(dl % 3 == 0);  //After taking out swapcount, number of arguments should be divisible by three.
    auto left = data[0..dl/3];
    auto right = data[dl/3..dl*2/3];
    auto result = data[dl*2/3..dl];
    static assert(left.length == right.length && right.length == result.length);
    size_t i = result[0].length - 1;
    //writef(left, "\t", right, "\t", " >> " );
    while(left[0].length && right[0].length) {
        if(right[0][$ - 1] >= left[0][$ - 1]) {
            foreach(ti, array; result) {
                result[ti][i] = right[ti][$-1];
                right[ti].length = right[ti].length - 1;
            }
        } else {
            *swapCount += right[0].length;
            foreach(ti, array; result) {
                result[ti][i] = left[ti][$-1];
                left[ti].length = left[ti].length - 1;
            }
        }
        i--;
    }
    if(right[0].length) {
        foreach(ti, array; result) {
            result[ti][0..right[ti].length] = right[ti];
        }
    } else {
        foreach(ti, array; result) {
            result[ti][0..left[ti].length] = left[ti];
        }
    }
   // writefln(result);
}

/**Bubble multisort function.  Useful only because it's easy to implement
 * and test other sorting functions against.  Counts swaps if last
 * argument is a ulong*. */

T[0] bubbleSort(T...)(T data) {
    //If last argument is a pointer to an integer type, increments this integer as a swapCount variable.  If it is anything
    //other than an array or a pointer to an integer, undefined behavior.
    static if(!isArray!(typeof(data[$-1]))) invariant uint dl = data.length - 1;
    else invariant uint dl = data.length;
    if(data[0].length < 2) return data[0];
    uint swapExecuted = void;
    foreach(i; 0..data[0].length) {
        swapExecuted = 0;
        foreach(j; 1..data[0].length) {
            if(data[0][j-1] > data[0][j]) {
                swapExecuted = 1;
                static if(!isArray!(typeof(data[$-1]))) (*(data[$-1]))++;
                foreach(array; data[0..dl]) swap(array[j-1], array[j]);
            }
        }
        if(!swapExecuted) return data[0];
    }
    return data[0];
}

/**Unstable selection sort.  Used for small subarrays in qsort. */

T[0] selectionSort(T...)(T data) {  //Used for small subarrays in qsort.
    if(data[0].length < 2) return data[0];
    foreach(i; 0..data[0].length-1) {
        size_t minRow = i;
        foreach(j; i+1..data[0].length) {
            if(data[0][j] < data[0][minRow]) {
                minRow = j;
            }
        }
        foreach(array; data) swap(array[i], array[minRow]);
    }
    return data[0];
}

/**Stable selection sort.  Used for small subarrays in stableQsort and mergeSort.
 * If last argument is a ulong*, keeps track of number of swaps that would have been
 * necessary if this were a bubble sort.*/

T[0] stableSelectionSort(T...)(T data) {  //Used for small subarrays in stableQsort, mergeSort.
    static if(!isArray!(typeof(data[$-1]))) invariant uint dl = data.length - 1;
    else invariant uint dl = data.length;
    if(data[0].length < 2) {return data[0];}
    foreach(i; 0..data[0].length-1) {
        size_t minRow = i;
        foreach(j; i+1..data[0].length) {
            if(data[0][j] < data[0][minRow]) {
                minRow = j;
            }
        }
        foreach(array; data[0..dl]) rotateRight(array[i..minRow+1]);
        static if(!isArray!(typeof(data[$-1]))) {
            (*(data[$-1])) += minRow - i;  //Incement swapCount variable.
        }
    }
    return data[0];
}


unittest {  //Parts of this test the sort functions against each other by checking for agreement.  Assumes they won't ALL be wrong in the same way.
    Random gen;
    gen.seed(unpredictableSeed);
    uint[] qsTest = new uint[1000];
    foreach(i, ref m; qsTest) m = uniform(gen, 0, 10000);
    foreach(i; 0..1000) {
        ulong bsCount, msCount, msInPlaceCount, ssCount;
        size_t lowerBound = uniform(gen, 0, qsTest.length);  //Need to test a bunch of different sizes, especially for merge sort.
        size_t upperBound = uniform(gen, 0, qsTest.length);
        if(upperBound < lowerBound) swap(upperBound, lowerBound);
        randomShuffle(qsTest, gen);
        uint[] msTest = qsTest.dup,  bsTest = qsTest.dup, ssTest = qsTest.dup, sssTest = qsTest.dup;
        qsTest[lowerBound..upperBound].qsort();
        msTest[lowerBound..upperBound].mergeSort(&msCount);
        sssTest[lowerBound..upperBound].stableSelectionSort(&ssCount);
        ssTest[lowerBound..upperBound].selectionSort();
        bsTest[lowerBound..upperBound].bubbleSort(&bsCount);
        assert(isSorted(qsTest[lowerBound..upperBound]));
        assert(isSorted(bsTest[lowerBound..upperBound]));
        assert(isSorted(ssTest[lowerBound..upperBound]));
        assert(isSorted(sssTest[lowerBound..upperBound]));
        assert(isSorted(msTest[lowerBound..upperBound]));
        assert(bsCount == msCount);
        assert(msCount == ssCount);
    }
    //Tesing for stability of multisorts:
    uint[] bOne = new uint[1000], bTwo = new uint[1000];
    foreach(i, ref m; bOne) m = uniform(gen, 0, 100);  //Small range virtually guarantees some ties.
    foreach(i, ref m; bTwo) m = uniform(gen, 0, 100);
    foreach(i; 0..1000) {  //Test stable sorts.
        size_t lowerBound = uniform(gen, 0, bOne.length);
        size_t upperBound = uniform(gen, 0, bOne.length);
        if(upperBound < lowerBound) swap(upperBound, lowerBound);
        randomShuffle(bOne, gen);
        randomShuffle(bTwo, gen);
        uint[] qOne = bOne.dup, qTwo = bTwo.dup, sOne = bOne.dup, sTwo = bTwo.dup, oO = bOne.dup, tO = bTwo.dup,
            mOne = bOne.dup, mTwo = bTwo.dup;
        bubbleSort(bOne, bTwo);
        stableSelectionSort(sOne, sTwo);
        mergeSort(mOne, mTwo);
        stableQsort(qOne, qTwo);
        assert(isSorted(sOne));
        assert(isSorted(bOne));
        assert(isSorted(qOne));
        assert(isSorted(mOne));
        assert(sTwo == bTwo);
        assert(bTwo == qTwo);
        assert(qTwo == mTwo);
    }
    writefln("Passed sort test.");
}