/******************************************************************************
* File: BinaryQuicksort.hh
* Author: Keith Schwarz (htiek@cs.stanford.edu)
*
* An implementation of the "binary quicksort" algorithm, another name for MSD-
* radix sort in base 2.  The algorithm is a non-comparison sort for sorting
* integers that works by sorting the numbers one bit at a time, starting with
* the most significant digit and ending with the least-significant digit.
* Internally, the algorithm works by first sorting the numbers into ascending
* order by the first bit, then recursively sorting each half of the numbers by
* their second digit, etc.  For example, given the numbers 0-7 in scrambled
* order, like this:
*
*         3   1   4   0   5   7   2   6
*
* We would write these numbers in binary to get
*
*       011 001 100 000 101 111 010 110
*
* First, we sort the numbers by just their first digit.  For example:
*
*       011 001 000 010 100 101 111 110
*       ^               ^
*       |               |
*       |               +--- Numbers beginning with a 1
*       |
*       +------------------- Numbers beginning with a 0
*
*
* Notice that we now have the array split into two regions - one region
* consisting of values beginning with a 0 and one region consisting of values
* beginning with a 1.  We now recursively descend into these two regions and
* sort each one by their second digit.  This gives the following:
*
*       001 000 011 010 100 101 111 110
*
* From here we sort by the third digit:
*
*       000 001 010 011 100 101 110 111
*
* And the numbers are now in order.
*
* In order to convert this sketch of an algorithm into a concrete
* implementation, we will need to devise a means for sorting a range of
* numbers by just a single digit.  One option for doing this out-of-place is
* to build two queues, one for numbers with a 0 in that bit and one for
* numbers with a 1 in at that bit, then add each number in the range to the
* appropriate queue.  We can then dequeue the numbers from the 0 queue and the
* 1 queue, in order, and insert them back into the original array.  For
* example, given these numbers, which we want to sort by the most significant
* bit:
*
*       011 001 100 000 101 111 010 110
*
* We would build two queues and insert the numbers as follows:
*
*   Begins with 0: 011 001 000 010
*   Begins with 1: 100 101 111 110
*
* We would then concatenate these two queues and store the result back into
* the original array:
*
*      011 001 000 010 100 101 111 110
*
* While this works, it has the drawback that it requires O(n) auxiliary
* storage.  Fortunately, there is a faster algorithm for solving this problem
* that works in-place.  The intuition behind the algorithm is that we want to
* reorder the elements so that the array ends up being split into two regions-
* a region on the left consisting of all values that contain a 0 at the
* appropriate bit, and a region on the right consisting of all values that
* contain a 1 at the appropriate bit.  We can therefore build an algorithm
* that scans inward from the two sides of the array looking for numbers that
* violate this invariant, then exchange the numbers that are out of place.
* For example, suppose we have this sequence of 0s and 1s:
*
*        0  1  0  1  0  1  1  0  1
*
* We would march two pointers in from the left and right side looking for a
* point where elements are found that are out of place.  This would first find
* this pair:
*
*        0  1  0  1  0  1  1  0  1
*           ^                 ^
*
* We exchange those numbers, as shown here:
*
*        0  0  0  1  0  1  1  1  1
*           ^                 ^
*
* and then march inward to find this next pair:
*
*        0  0  0  1  0  1  1  1  1
*                 ^  ^
*
* which we swap to obtain
*
*        0  0  0  0  1  1  1  1  1
*                 ^  ^
*
* and at this point are done.  This algorithm is very similar to the partition
* step found in most quicksort implementations, and for this reason the MSD-
* radix sort is often called binary quicksort.
*
* To analyze the runtime of the algorithm, note that running the partitioning
* step on k integers takes O(k) time, since each element is visited at most
* O(1) times.  Initially, this means that we do O(n) work to partition the
* values by their first bit, then a total of O(n) work to partition each half
* by the second bit, etc.  More generally, if the largest possible value that
* can appear is U, there will be O(log U) partioning steps required, one for
* each of the bits in U.  This means that the total runtime is O(n log U).
* In cases where U = O(n), this is asympototically better than or comparable
* to comparison sorts like quicksort or heap sort.
*/

#ifndef BinaryQuicksort_Included
#define BinaryQuicksort_Included

#include <climits>   // For CHAR_BIT
#include <iterator>
#include <limits>
#include <algorithm> // For std::iter_swap, std::rotate, std::find_if

/**
* Function: BinaryQuicksort(RandomIterator begin, RandomIterator end);
* Usage: BinaryQuicksort(v.begin(), v.end());
* ------------------------------------------------------------------------
* Applies the binary quicksort algorithm to sort the specified list of
* numbers.  It is assumes that the iterators are traversing a list of
* integral types, and will not function properly otherwise.
*/

template <typename RandomIterator>
void BinaryQuicksort(RandomIterator begin, RandomIterator end);

/* * * * * Implementation Below This Point * * * * */

namespace binaryquicksort_detail {

/* Utility function to partition the elements of a range by moving all
* elements in the range having a 0 in a given position to the right and all
* elements in the range having a 1 in a given position to the left.  The
* function then returns an iterator to the beginning of the range that
* contains a 1.
*
* This algorithm works by having begin point one step past the end of the
* range of values known to be 0 and end point at the range of values known
* to be 1.  The endpoints are then marched inward until they collide (in
* which case we're done) or a pair of mismatched elements are found.
*/

template <typename RandomIterator>
RandomIterator PartitionAtBit(RandomIterator begin, RandomIterator end,
signed int bit) {
/* Typedef defining the type of the elements being traversed. */
typedef typename std::iterator_traits<RandomIterator>::value_type T;

/* Compute the bitmask we'll use to test whether the bit is set. */
const T bitmask = T(1) << bit;

/* Move these two together until they meet or we find two elements that
* are out of place.
*/

while (true) {
/* Find the first 1 after the 0s; it's either the end or we've just
* found the element that's out of place.
*/

while (begin < end && !(*begin & bitmask))
++ begin;

/* If the begin is now sitting atop the end, we're done and all of the
* remaining values are 1s.  We're therefore done.
*/

if (begin == end)
return begin;

/* Otherwise, the end is just before the elements containing 1s.  Start
* moving it over until we find a 0 or realize that we've already
* found everything.
*
* Since we need to back up the end iterator before we read it, this
* loop is written as a do ... while rather than as a while loop.
*/

do {
--end;
} while (begin < end && !!(*end & bitmask));

/* If the two are equal, we've found the crossover point and are done.
* We can hand back this element as the pivot point.
*/

if (begin == end)
return begin;

/* Otherwise, swap the two elements and repeat. */
std::iter_swap(begin, end);
}
}

/* Utility function which actually performs the binary quicksort algorithm,
* beginning with the specified bit.
*/

template <typename RandomIterator>
void BinaryQuicksortAtBit(RandomIterator begin, RandomIterator end,
signed int bit) {
/* Typedef defining the type of the elements being traversed. */
typedef typename std::iterator_traits<RandomIterator>::value_type T;

/* Borrowing an optimization technique from quicksort, we will have this
* function work iteratively and recursively.  To avoid having a large
* number of function calls made, we will iteratively process the larger
* of the two partitions we find, and will recursively process the other.
*/

/* If we've processed all the bits, or if the range has fewer than one
* element in it, we're done.
*/

while (bit >= 0 && std::distance(begin, end) > 1) {
/* Apply the partitioning step on this bit and get the start of the
* range of values containing the 1s.
*/

RandomIterator pivot = PartitionAtBit(begin, end, bit);

/* Drop the index of the bit we're processing; this will cause the next
* loop iteration to use the right bit and will make the recursive calls
* correct.
*/

--bit;

/* Determine which range is larger - the range holding the 0s or the
* range holding the 1s.  Based on which is smaller, recursively process
* one of the ranges.
*/

if (std::distance(begin, pivot) < std::distance(pivot, end)) {
/* There are fewer numbers beginning with 0; go recursively sort
* them.
*/

BinaryQuicksortAtBit(begin, pivot, bit);
begin = pivot;
} else {
/* There are fewer numbers beginning with 1; go recursively sort
* them.
*/

BinaryQuicksortAtBit(pivot, end, bit);
end = pivot;
}
}
}

/* If we are dealing with signed numbers, then negative numbers will
* incorrectly appear at the end of the range rather than the start, since
* the signed two's-complement representation will cause the sign bit to
* be set, making the negative values appear larger than positive values.
* This function applies a rotation to the final array to pull the negative
* values (if any) to the front.
*/

template <typename RandomIterator>
void RotateNegativeValues(RandomIterator begin, RandomIterator end) {
/* Typedef defining the type of the elements being traversed. */
typedef typename std::iterator_traits<RandomIterator>::value_type T;

/* Walk forward until we find a negative value.  If we find one, do a
* rotate to rectify the elements.
*/

for (RandomIterator itr = begin; itr != end; ++itr) {
/* If the value is negative, do a rotate starting here. */
if (*itr < T(0)) {
std::rotate(begin, itr, end);
return;
}
}
}
}

/* Actual implementation of binary quicksort. */
template <typename RandomIterator>
void BinaryQuicksort(RandomIterator begin, RandomIterator end) {
/* Typedef defining the type of the elements being traversed. */
typedef typename std::iterator_traits<RandomIterator>::value_type T;

/* Find out how many bits we need to process. */
const signed int kNumBits = (signed int)(CHAR_BIT * sizeof(T));

/* Run binary quicksort on the elements, starting with the MSD. */
binaryquicksort_detail::BinaryQuicksortAtBit(begin, end, kNumBits - 1);

/* If the numbers are signed, we need to do a rotate to pull all of the
* negative numbers to the front of the range, since otherwise (because
* their MSB is set) they'll be at the end instead of the front.
*/

if (std::numeric_limits<T>::is_signed)
binaryquicksort_detail::RotateNegativeValues(begin, end);
}

#endif