/***********************************************************************\
*                              sortbench.d                              *
*                                                                       *
*               Benchmark of sorting algorithms on arrays               *
*                                                                       *
*                           by Stewart Gordon                           *
*                             February 2007                             *
\***********************************************************************/
import std.random;
import std.date;
import std.stdio;
import std.stream;
import std.conv;
import std.path;
import std.cstream;

alias void function(uint[]) Algorithm;

template Data(uint len) {
	struct Data {
		uint[len][] random;
		uint[len] sorted;
		uint[len] reversed;

		void generate() {
			random.length = lists;

			for (int index = 0; index < len; index++) {
				foreach (inout uint[len] list; random) {
					list[index] = rand();
				}
				sorted[index] = index;
				reversed[index] = len - index;
			}
		}

		debug (100) {
			void print() {
				foreach (int iList, uint[len] list; random) {
					writefln("Random list %d:", iList);
					printList(list);
				}
				dout.writeLine("Sorted list:");
				printList(sorted);
				dout.writeLine("Reversed list:");
				printList(reversed);
			}
		}

		double benchRandom(Algorithm alg) {
			d_time start = getUTCtime();

			foreach (uint[len] list; random) {
				alg(list.dup);
				dout.write('.');
				dout.flush();
			}

			double result = cast(double) (getUTCtime() - start)
			  / TicksPerSecond / lists;

			writefln("\n%5d data: %f", len, result);
			if (csvOutput) csvOutput.writef(",%f", result);

			return result;
		}

		double benchSorted(Algorithm alg) {
			d_time start = getUTCtime();

			for (uint time = 0; time < repeats; time++) {
				/*	no need to dup sorted, as the before and after orders are
				 *	the same
				 */
				alg(sorted);
				dout.write('.');
				dout.flush();
			}

			double result = cast(double) (getUTCtime() - start)
			  / TicksPerSecond / repeats;

			writefln("\n%5d data: %f", len, result);
			if (csvOutput) csvOutput.writef(",%f", result);

			return result;
		}

		double benchReversed(Algorithm alg) {
			d_time start = getUTCtime();

			for (uint time = 0; time < repeats; time++) {
				alg(reversed.dup);
				dout.write('.');
				dout.flush();
			}

			double result = cast(double) (getUTCtime() - start)
			  / TicksPerSecond / repeats;

			writefln("\n%5d data: %f", len, result);
			if (csvOutput) csvOutput.writef(",%f", result);

			return result;
		}
	}
}

debug (100) {
	Data!(32) data;
} else {
	Data!(4096) data1;
	Data!(8192) data2;
	Data!(16384) data3;
}

uint lists = 64;
uint repeats = 32;
File csvOutput;


int main(char[][] args) {
	d_time start = getUTCtime();
	char[] csvName;

	foreach (arg; args[1..$]) {
		if (arg[0] == '-') arg = arg[1..$];
		if (arg.length < 1) {
			return showUsage();
		}
		switch (arg[0]) {
			case 'r': case 'R':
				repeats = toUint(arg[1..$]);
				break;

			case 'l': case 'L':
				lists = toUint(arg[1..$]);
				break;

			case 'o': case 'O':
				if (arg.length == 1) {
					csvName = "sortbench.csv";
				} else {
					csvName = defaultExt(arg[1..$], "csv");
				}
				break;

			default:
				return showUsage();
		}
	}

	debug (100) {
		data.generate();

		sortBench("Naive bubble sort", &naiveBubbleSort);
	} else {
		data1.generate();
		data2.generate();
		data3.generate();

		if (csvName.length != 0) {
			csvOutput = new File(csvName, FileMode.OutNew);
			csvOutput.writeLine(`"Algorithm",`
			  `"4096 random","8192 random","16384 random",`
			  `"4096 sorted","8192 sorted","16384 sorted",`
			  `"4096 reversed","8192 reversed","16384 reversed"`);
		}

		sortBench("Naive bubble sort", &naiveBubbleSort);
		sortBench("Flagged bubble sort", &flaggedBubbleSort);
		sortBench("Bubble sort with last swap memory", &lastSwapBubbleSort);
		sortBench("Bidirectional bubble sort", &bidirBubbleSort);
		sortBench("Comb sort", &combSort);

		sortBench("Selection sort", &selectionSort);
		sortBench("Shaker sort", &shakerSort);

		sortBench("Insertion sort", &insertionSort);
		sortBench("Buffered insertion sort", &bufferedInsertionSort);

		sortBench("Quick sort", &quickSort);
		sortBench("Median-of-three quick sort", &quickerSort);

		sortBench("Merge sort", &mergeSort);
		sortBench("Heap sort", &heapSort);
		sortBench("Smooth sort", &smoothSort);
		sortBench("Binary radix sort", &radixSort2);
		sortBench("Quaternary radix sort", &radixSort4);

		sortBench("Built-in sort", &builtInSort);
	}

	// finished - display elapsed time
	d_time elapsed = getUTCtime() - start;
	writefln("\nTotal time: %d:%02d.%03d", elapsed / TicksPerSecond / 60,
	  (elapsed / TicksPerSecond) % 60, elapsed % TicksPerSecond);

	return 0;
}


int showUsage() {
	derr.writeLine("Sort benchmark by Stewart Gordon
Usage: sortbench [options]
where options are
	-l<number>  number of lists of each length to sort (default 64)
	-r<number>  number of times to repeat sorted/reversed sort (default 32)
	-o          output a CSV file of the results as sortbench.csv
	-o<name>    output a CSV file of the results");
	return 1;
}


void printList(uint[] list) {
	foreach (uint datum; list) writef("%d ", datum);
	writefln();
}


void sortBench(char[] algName, Algorithm alg) {
	writefln();
	dout.writeLine(algName);
	foreach (char c; algName) dout.write('-');

	debug (100) {
		dout.writeLine("\nSorting random data");
		writefln("32 data: %f", data.benchRandom(alg));

		dout.writeLine("Sorting sorted data");
		writefln("32 data: %f", data.benchSorted(alg));

		dout.writeLine("Sorting reversed data");
		writefln("32 data: %f", data.benchReversed(alg));
	} else {
		if (csvOutput) csvOutput.writef(`"%s"`, algName);

		writefln("\nSorting random data");
		data1.benchRandom(alg);
		data2.benchRandom(alg);
		data3.benchRandom(alg);

		writefln("Sorting sorted data");
		data1.benchSorted(alg);
		data2.benchSorted(alg);
		data3.benchSorted(alg);

		writefln("Sorting reversed data");
		data1.benchReversed(alg);
		data2.benchReversed(alg);
		data3.benchReversed(alg);

		if (csvOutput) csvOutput.writefln();
	}
}



void checkSorted(uint[] data) {
	for (int index = 0; index < data.length - 1; index++) {
		assert (data[index] <= data[index + 1]);
	}
}


void swap(inout uint x, inout uint y) {
	uint temp = x;
	x = y;
	y = temp;
}


struct Portion {
	uint* bottom, top;
}


/+
/*	Stooge sort is a bottom-up recursive algorithm that works by recursively
 *	sorting the first two thirds, the last two thirds and then the first two
 *	thirds again.  I experimented with including it in the benchmark, but
 *	found it far too slow for the amount of data involved.
 *	Time complexity:
 *		O(n^2.7) (approx) all cases
 *	Space complexity:
 *		O(log n) all cases
 */
void stoogeSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	// the stack of portions to sort
	Portion[] stack = new Portion[1];
	stack[0].bottom = &data[0];
	stack[0].top = &data[$-1];

	// repeat while there is something left to do
	while (stack.length != 0) {
		Portion current = stack[$ - 1];

		with (current) {
			if (*bottom > *top) swap(*bottom, *top);

			if (top - bottom < 2) {
				stack.length = stack.length - 1;
			} else {
				uint subLen = (2 * (top - bottom) + 1) / 3;

				stack.length = stack.length + 2;

				stack[$-2].top = top;
				stack[$-2].bottom = top - subLen;
				stack[$-3].top = bottom + subLen;
				stack[$-1] = stack[$-3];
			}
		}
	}
}
+/


/*	The simplest implementation of bubble sort.  Passes decrease in length
 *	by one each time, and continue until the pass reaches length 2.
 *	Time complexity:
 *		O(n^2) all cases
 *	Working memory:
 *		O(1) all cases
 */
void naiveBubbleSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	for (int pass = data.length - 1; pass > 0; pass--) {
		for (int index = 0; index < pass; index++) {
			if (data[index] > data[index + 1]) {
				uint temp = data[index + 1];
				data[index + 1] = data[index];
				data[index] = temp;
			}
		}
	}
}


/*	A flag is used to note whether any swaps have taken place in the current
 *	pass.  The algorithm halts when a pass is made with no swaps.
 *	Time complexity:
 *		O(n^2) worst case
 *		O(n) best case
 *	Working memory:
 *		O(1) all cases
 */
void flaggedBubbleSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	for (int pass = data.length - 1; pass > 0; pass--) {
		bit anySwaps = false;

		for (int index = 0; index < pass; index++) {
			if (data[index] > data[index + 1]) {
				anySwaps = true;
				swap(data[index], data[index + 1]);
			}
		}
		if (!anySwaps) return;
	}
}


/*	A further optimisation of flagged bubble sort, in which the actual
 *	position of the last swap is remembered.  Subsequent passes are then
 *	shortened to this position.
 *	Time complexity:
 *		O(n^2) worst case
 *		O(n) best case
 *	Working memory:
 *		O(1) all cases
 */
void lastSwapBubbleSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	int stopAt = data.length - 1;

	while (stopAt != 0) {
		int lastSwap = 0;

		for (int index = 0; index < stopAt; index++) {
			if (data[index] > data[index + 1]) {
				lastSwap = index;
				swap(data[index], data[index + 1]);
			}
		}
		stopAt = lastSwap;
	}
}


/*	I thought of this variation of bubble sort some years ago, as an
 *	algorithm that would be suited to an application that I was developing.
 *	This was partly because the program loaded a window of records into
 *	memory at a time, and partly because the algorithm works especially well
 *	on data that is already almost sorted.  But it's intuitive enough that
 *	several people have come up with it independently.  I've also seen it
 *	called 'cocktail shaker sort' among other names.  This implementation
 *	uses the aforementioned last swap optimisation.
 *	Time complexity:
 *		O(n^2) worst case
 *		O(n) best case
 *	Working memory:
 *		O(1) all cases
 */
void bidirBubbleSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	int stopForward = data.length - 1;
	int stopBackward = -1;
	int lastSwap = 0;

	while (stopForward > stopBackward) {
		for (int index = stopBackward + 1; index < stopForward; index++) {
			if (data[index] > data[index + 1]) {
				lastSwap = index;
				swap(data[index], data[index + 1]);
			}
		}
		stopForward = lastSwap;

		for (int index = stopForward - 1; index > stopBackward; index--) {
			if (data[index] > data[index + 1]) {
				lastSwap = index;
				swap(data[index], data[index + 1]);
			}
		}
		stopBackward = lastSwap;
	}
}


/*	Another simplistic algorithm.  Though on average slightly faster than
 *	bubble sort thanks to fewer swaps, it unfortunately doesn't optimise
 *	well for already sorted or almost sorted lists.
 *	Time complexity:
 *		O(n^2) all cases
 *	Working memory:
 *		O(1) all cases
 */
void selectionSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	for(int pass = 0; pass < data.length - 1; pass++) {
		// set minimum to first element
		int minIndex = -1;	// index relative to slice
		uint min = data[pass];

		foreach (int index, uint datum; data[pass+1..data.length]) {
			if (datum < min) {
				minIndex = index;
				min = datum;
			}
		}

		if (minIndex >= 0) {
			data[minIndex + pass + 1] = data[pass];
			data[pass] = min;
		}
	}
}


/*	A double selection sort, picking out both the minimum and maximum
 *	elements in each pass.  Not to be confused with the bidirectional
 *	bubble sort, which I've also seen called 'shaker sort'.
 *	Time complexity:
 *		O(n^2) all cases
 *	Working memory:
 *		O(1) all cases
 */
void shakerSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	for(int pass = 0; pass < (data.length) / 2; pass++) {
		debug (100) {
			printList(data);
			getch();
		}
		// set minimum to first element
		int minIndex = -1, maxIndex = -1;	// index relative to slice
		uint min = data[pass], max = min;

		foreach (int index, uint datum; data[pass+1..data.length-pass]) {
			if (datum < min) {
				minIndex = index;
				min = datum;
			} else if (datum > max) {
				maxIndex = index;
				max = datum;
			}
		}

		debug (100) {
			writefln("pass %d, min %d, max %d", pass, minIndex, maxIndex);
		}

		if (minIndex >= 0) {
			data[minIndex + pass + 1] = data[pass];
			data[pass] = min;
		}

		if (maxIndex == -1) {
			// in which case it's been moved by the min swap
			data[minIndex + pass + 1] = data[data.length - pass - 1];
			data[data.length - pass - 1] = max;

		} else if (maxIndex < data.length - 2 * pass - 2) {
			data[maxIndex + pass + 1] = data[data.length - pass - 1];
			data[data.length - pass - 1] = max;
		}
	}
}


/*	The basic quicksort algorithm, in which the first element of each
 *	portion is chosen as the pivot.
 *	Time complexity:
 *		O(n^2) worst case
 *		O(n log n) average
 *	Working memory:
 *		O(log n) worst case
 *		O(1) best case
 */
void quickSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	// the stack of portions to sort
	Portion[] stack = new Portion[1];
	stack[0].bottom = data.ptr;
	stack[0].top = data.ptr + data.length;

	// repeat while there is something left to do
	while (stack.length != 0) {
		debug (100) {
			foreach (Portion p; stack) {
				writef("%d..%d  ", p.bottom - data.ptr, p.top - data.ptr);
			}
			putchar('\n');
			printList(data);
			getch();
		}
		Portion current = stack[$ - 1];
		// do the current pass
		uint pivot = *current.bottom;
		uint* lowUp = current.bottom;
		uint* highDown = current.top;

		while (lowUp != highDown) {
			// gap at lowUp - move highDown
			for (highDown--; lowUp != highDown && *highDown >= pivot;
			  highDown--) {}
			if (lowUp == highDown) break;
			*lowUp = *highDown;

			// gap at highDown - move lowUp
			for (lowUp++; lowUp != highDown && *lowUp <= pivot; lowUp++) {}
			*highDown = *lowUp;
		}

		*lowUp = pivot;

		// now update the stack
		int len = stack.length;

		if (current.top - current.bottom == 3
			  && lowUp != current.bottom + 1) {
			// just possibly a swap needed
			if (current.bottom[0] > current.bottom[1]) {
				swap(current.bottom[0], current.bottom[1]);
			} else if (current.bottom[1] > current.bottom[2]) {
				swap(current.bottom[1], current.bottom[2]);
			}

			// now it's sorted
			stack.length = len - 1;

		} else if (current.top - current.bottom < 4) {
			// portion sorted
			stack.length = len - 1;

		} else if (lowUp - current.bottom >= current.top - lowUp) {
			// gap towards top - push lower half first
			stack[len - 1].top = lowUp;
			if (current.top - lowUp >= 3) {
				stack.length = len + 1;
				stack[len].bottom = lowUp + 1;
				stack[len].top = current.top;
			}

		} else {
			// gap towards bottom - push upper half first
			stack[len - 1].bottom = lowUp + 1;
			if (lowUp - current.bottom >= 2) {
				stack.length = len + 1;
				stack[len].top = lowUp;
				stack[len].bottom = current.bottom;
			}
		}
	}
}


/*	A variation of quicksort, in which the median of the first three elements
 *	is chosen as the pivot.
 *	Time complexity:
 *		O(n^2) worst case
 *		O(n log n) average
 *	Working memory:
 *		O(log n) worst case
 *		O(1) best case
 */
void quickerSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	// the stack of portions to sort
	Portion[] stack = new Portion[1];
	stack[0].bottom = data.ptr;
	stack[0].top = data.ptr + data.length;

	// repeat while there is something left to do
	while (stack.length != 0) {
		debug (100) {
			printList(data);
			foreach (Portion p; stack) {
				writef("%d..%d  ", p.bottom - data.ptr, p.top - data.ptr);
			}
			writefln("\n");
			getch();
		}
		Portion current = stack[$ - 1];

		// do the current pass
		uint* lowUp = current.bottom;
		uint* highDown = current.top;
		uint pivot = *current.bottom;

		// find median of 3
		if (highDown - lowUp >= 3) {
			if (lowUp[0] < lowUp[1] && lowUp[0] < lowUp[2]) {
				// median is one of the other two
				if (lowUp[1] <= lowUp[2]) {
					// lowUp[1] is the median - take as pivot
					pivot = lowUp[1];
					lowUp[1] = lowUp[0];
				} else {
					// lowUp[2] is the median
					pivot = lowUp[2];
					lowUp[2] = lowUp[0];
				}
			} else if (lowUp[0] > lowUp[1] && lowUp[0] > lowUp[2]) {
				// median is one of the other two
				if (lowUp[1] >= lowUp[2]) {
					// lowUp[1] is the median - take as pivot
					pivot = lowUp[1];
					lowUp[1] = lowUp[0];
				} else {
					// lowUp[2] is the median
					pivot = lowUp[2];
					lowUp[2] = lowUp[0];
				}
			}
		}

		while (lowUp != highDown) {
			// gap at lowUp - move highDown
			for (highDown--; lowUp != highDown && *highDown >= pivot;
			  highDown--) {}
			if (lowUp == highDown) break;
			*lowUp = *highDown;

			// gap at highDown - move lowUp
			for (lowUp++; lowUp != highDown && *lowUp <= pivot; lowUp++) {}
			*highDown = *lowUp;
		}

		*lowUp = pivot;

		// now update the stack
		int len = stack.length;

		if (current.top - current.bottom == 3
			  && lowUp != current.bottom + 1) {
			// just possibly a swap needed
			if (current.bottom[0] > current.bottom[1]) {
				swap(current.bottom[0], current.bottom[1]);
			} else if (current.bottom[1] > current.bottom[2]) {
				swap(current.bottom[1], current.bottom[2]);
			}

			// now it's sorted
			stack.length = len - 1;

		} else if (current.top - current.bottom < 4) {
			// portion sorted
			stack.length = len - 1;

		} else if (lowUp - current.bottom >= current.top - lowUp) {
			// gap towards top - push lower half first
			stack[len - 1].top = lowUp;
			if (current.top - lowUp >= 3) {
				stack.length = len + 1;
				stack[len].bottom = lowUp + 1;
				stack[len].top = current.top;
			}

		} else {
			// gap towards bottom - push upper half first
			stack[len - 1].bottom = lowUp + 1;
			if (lowUp - current.bottom >= 2) {
				stack.length = len + 1;
				stack[len].top = lowUp;
				stack[len].bottom = current.bottom;
			}
		}
	}
}


/*	Buffered merge sort.  There's also an in-place merge sort, but that's
 *	O(n^2) at least in the form I've seen it.  This implementation assumes,
 *	as is true in the benchmarks, that data.length is a power of 2.
 *	Time complexity:
 *		O(n log n) all cases
 *	Working memory:
 *		O(n) all cases
 */
void mergeSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	int blockSize = 2;
	uint[] aux = new uint[data.length];

	void doMerge(uint[] fromFirst, uint[] fromSecond, uint[] to) {
		uint* inFirst = fromFirst.ptr, inSecond = fromSecond.ptr;
		uint* endFirst = inFirst + fromFirst.length;
		uint* endSecond = inSecond + fromSecond.length;

		foreach (int index, inout uint dest; to) {
			if (inFirst == endFirst) {
				to[index..to.length] = inSecond[0..to.length-index];
				break;
			} else if (inSecond == endSecond) {
				to[index..to.length] = inFirst[0..to.length-index];
				break;
			} else if (*inFirst <= *inSecond) {
				dest = *(inFirst++);
			} else {
				dest = *(inSecond++);
			}
		}
	}

	void doPass(uint[] from, uint[] to) {
		for (int index = 0; index < from.length; index += 2 * blockSize) {
			doMerge(from[index..index + blockSize],
				    from[index + blockSize..index + 2 * blockSize],
				    to[index..index + 2 * blockSize]);
		}
		blockSize *= 2;
	}

	// order pairs first
	for (int index = 0; index < data.length; index += 2) {
		if (data[index] <= data[index + 1]) {
			aux[index..index+2] = data[index..index+2];
		} else {
			aux[index] = data[index + 1];
			aux[index + 1] = data[index];
		}
	}

	assert (data.length >= 4);
	doPass(aux, data);

	while (blockSize < data.length) {
		doPass(data, aux);
		if (blockSize < data.length) {
			doPass(aux, data);
		} else {
			data[] = aux;
		}
	}
}


/*	Another bubble sort variation, this one makes passes in which the
 *	compared elements are spaced apart, by a smaller distance with each
 *	successive pass.  This is a surprisingly fast algorithm.
 *	Working memory:
 *		O(1) all cases
 */
void combSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	const float FACTOR = 0.7;

	for (int space = cast(int) (data.length * FACTOR);
		  space > 1; space *= FACTOR) {
		for (int index = 0; index < data.length - space; index++) {
			if (data[index] > data[index + space]) {
				swap(data[index], data[index + space]);
			}
		}
	}

	lastSwapBubbleSort(data);
}


/*	Insertion sorts are efficient on linked lists, but this is an array.
 *	This algorithm has quite a bit in common with bubble sort.
 *	Time complexity:
 *		O(n^2) worst case
 *		O(n) best case
 *	Working memory:
 *		O(1) all cases
 */
void insertionSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	foreach (int index, uint datum; data[1..data.length]) {
		debug (100) {
			printList(data);
			getch();
		}

		int pos;
		for (pos = index; pos >= 0 && data[pos] > datum; pos--) {
			data[pos + 1] = data[pos];
		}

		data[pos + 1] = datum;
	}
}


/*	This version uses a buffer to transfer the array contents.  It makes use
 *	of the binary search technique to find the insertion point.  It could be
 *	implemented using .dup instead of a pre-allocated buffer, thereby making
 *	O(n) merely the worst case working memory, but doing so slows things down
 *	noticeably.
 *	Time complexity:
 *		O(n^2) worst case
 *		O(n log n) best case
 *	Working memory:
 *		O(n) all cases
 */
void bufferedInsertionSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	uint[] buffer = new uint[data.length];

	// do first 2 elements straight off for efficiency
	if (data[0] > data[1]) swap(data[0], data[1]);

	foreach (int index, uint datum; data[1..data.length]) {
		// do binary search
		int highDown = index + 1;
		int lowUp = 0;

		while (lowUp != highDown) {
			int midpoint = (highDown + lowUp) / 2;

			if (data[midpoint] == datum) {			// matching item
				lowUp = midpoint + 1;
				break;
			} else if (data[midpoint] > datum) {	// datum goes before
				highDown = midpoint;
			} else {								// datum goes after
				lowUp = midpoint + 1;
			}
		}

		// found the space - do the insertion
		buffer[lowUp..index + 1] = data[lowUp..index + 1];
		data[lowUp + 1..index + 2] = buffer[lowUp..index + 1];
		//data[lowUp + 1..index + 2] = data[lowUp..index + 1].dup;
		data[lowUp] = datum;
	}
}


/*	A slight variation on the standard textbook heapsort, in which the first
 *	element is omitted from the initial building of the heap, and then
 *	swapped with the minimum element.  This is purely so that indexes into
 *	the heap form a nice, easy to process pattern.
 *	Time complexity:
 *		O(n log n) all cases
 *	Working memory:
 *		O(1) all cases
 */
void heapSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	// Bubble an element down to its position in the heap
	void bubbleDown(int start, int len) {
		for (int index = start; index < ((len + 1) >> 1);) {
			// find the max element in the triangle
			if (((index << 1) | 1) == len
				  || data[index] >= data[(index << 1) | 1]) {
				// it's not the next right, so try the left
				if (data[index] >= data[index << 1]) {
					// this subtree already heaped
					break;
				}
				// next left is the max
				//int temp = data[index];
				swap(data[index], data[index << 1]);
				//swap(index, index << 1);
				index <<= 1;

			} else {
				// next right > current, but is it the max?
				// take the greater of next left and next right
				int next = (index << 1) |
				  (data[index << 1] <= data[(index << 1) | 1]);
				swap(data[index], data[next]);
				//swap(index, next);
				index = next;
			}
		}
	}

	// first build the heap
	int halfLength = (data.length + 1) >> 1;

	// do it bottom-up, to take advantage of subtrees being already heaped
	for (int index = halfLength - 1; index >= 1; index--) {
		bubbleDown(index, data.length);
	}

	// now pick out the smallest element, and put it at the beginning
	int minIndex = 0;
	for (int index = halfLength; index < data.length; index++) {
		if (data[index] < data[minIndex]) minIndex = index;
	}
	if (minIndex != 0) {
		//swap(minIndex, 0);
		swap(data[minIndex], data[0]);

		// bubble up item moved from beginning
		for (int index = minIndex;
			  index > 1 && data[index] > data[index >> 1]; index >>= 1) {
			//swap(index, index >> 1);
			swap(data[index], data[index >> 1]);
		}
	}

	// now do the sort
	for (int leftLength = data.length - 1; leftLength >= 4; leftLength--) {
		//swap(1, leftLength);
		swap(data[1], data[leftLength]);
		bubbleDown(1, leftLength);
	}

	// only top 3 remain - finish it off
	swap(data[1], data[3]);
	if (data[1] > data[2]) swap(data[1], data[2]);
}


/*	An algorithm developed by Edsger Dijkstra based on the same principle as
 *	heapsort.  By using a postorder traversal of a heap-like structure, O(n)
 *	time complexity is achieved for already-sorted data.
 *	Time complexity:
 *		O(n log n) worst case
 *		O(n) best case
 *	Working memory:
 *		O(1) all cases
 */
void smoothSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	void up(inout int b, inout int c) {
		int temp = b + c + 1;
		c = b;
		b = temp;
	}

	void down(inout int b, inout int c) {
		int temp = b - c - 1;
		b = c;
		c = temp;
	}

	void sift(int r, int b, int c) {
		while (b >= 3) {
			int r2 = r - b + c;

			if (data[r2] < data[r - 1]) {
				r2 = r - 1;
				down(b, c);
			}

			if (data[r] >= data[r2]) {
				b = 1;
			} else {
				swap(data[r], data[r2]);
				r = r2;
				down(b, c);
			}
		}
	}

	void trinkle(int r, int p, int b, int c) {
		while (p > 0) {
			int r3;

			while (p % 2 == 0) {
				p /= 2;
				up(b, c);
			}
			r3 = r - b;

			if (p == 1 || data[r3] <= data[r]) {
				p = 0;
			} else {
				p--;
				if (b == 1) {
					swap(data[r], data[r3]);
					r = r3;
				} else if (b >= 3) {
					int r2 = r - b + c;
					if (data[r2] < data[r - 1]) {
						r2 = r - 1;
						down(b, c);
						p *= 2;
					}

					if (data[r3] >= data[r2]) {
						swap(data[r], data[r3]);
						r = r3;
					} else {
						swap(data[r], data[r2]);
						r = r2;
						down(b, c);
						p = 0;
					}
				}
			}
		}
		sift(r, b, c);
	}

	void semitrinkle(int r, int p, int b, int c) {
		int r1 = r - c;
		if (data[r1] > data[r]) {
			swap(data[r], data[r1]);
			trinkle(r1, p, b, c);
		}
	}

	void swap(inout int x, inout int y) {
		int temp = x;
		x = y;
		y = temp;
	}

	int q = 1, r, p = 1, b = 1, c = 1;

	while (q != data.length) {
		if (p % 8 == 3) {
			sift(r, b, c);
			p = (p+1)/4;
			up(b, c); up(b, c);
		} else if (p % 4 == 1) {
			if (q + c < data.length) {
				sift(r, b, c);
			} else {
				trinkle(r, p, b, c);
			}
			do {
				down(b, c);
				p *= 2;
			} while (b != 1);
			p++;
		}
		q++; r++;
	}
	trinkle(r, p, b, c);

	while (q != 1) {
		q--;
		if (b == 1) {
			r--; p--;
			while (p % 2 == 0) {
				p /= 2;
				up(b, c);
			}
		} else if (b >= 3) {
			p--;
			r += c - b;
			if (p != 0) semitrinkle(r, p, b, c);
			down(b, c);
			p = 2*p + 1;
			r += c;
			semitrinkle(r, p, b, c);
			down(b, c);
			p = 2*p + 1;
		}
	}
}


/*	Right-to-left binary radix sort.  Radix sorts are unusual among sorting
 *	algorithms in that it works by inspecting the data itself, and not just
 *	comparing the items with each other.  In this variant, values are stably
 *	sorted by each bit in turn, starting with the least significant bit.
 *	Time complexity:
 *		O(n) all cases
 *	Working memory:
 *		O(n) all cases
 */
void radixSort2(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	uint[][2] bucket;
	uint[] temp;

	// reserve space at the beginning to save overhead
	bucket[0].length = bucket[1].length = data.length;
	bucket[0].length = bucket[1].length = 0;

	for (uint sortBit = 1; sortBit != 0; sortBit <<= 1) {
		//debug writefln(sortBit);
		// sort into buckets
		foreach (uint value; data) bucket[!!(value & sortBit)] ~= value;

		// now join them together
		bucket[0] ~= bucket[1];

		/*	bucket[1] doesn't need its data at this time; however, to use it
		 *	as swapping space would mean losing one of the pre-allocated
		 *	buffers
		 */
		temp = data;
		data = bucket[0];
		bucket[0] = temp;
		bucket[0].length = bucket[1].length = 0;
	}

	/*	With successive passes, the data being sorted is passed between the
	 *	buffer passed in as data and the buffer initially allocated to
	 *	bucket[0].  This is done an even number of times, and so in the end,
	 *	it ends up in the original data buffer.
	 */
}


/*	Right-to-left quaternary radix sort.  Values are stably sorted by each
 *	crumb in turn, starting with the least significant crumb.
 *	Time complexity:
 *		O(n) all cases
 *	Working memory:
 *		O(n) all cases
 */
void radixSort4(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	uint[][4] bucket;
	uint[] temp;

	// reserve space at the beginning to save overhead
	bucket[0].length = data.length;
	bucket[0].length = 0;
	foreach (inout b; bucket[1..$]) {
		b.length = data.length / 2;
		b.length = 0;
	}

	for (uint sortCrumb = 3, bitIndex = 0; sortCrumb != 0;
		  sortCrumb <<= 2, bitIndex += 2) {
		debug (100) writefln(bitIndex);
		// sort into buckets
		foreach (uint value; data) {
			bucket[(value & sortCrumb) >> bitIndex] ~= value;
		}

		// now join them together
		foreach (b; bucket[1..$]) bucket[0] ~= b;

		temp = data;
		data = bucket[0];
		bucket[0] = temp;
		foreach (inout b; bucket) b.length = 0;
	}

	/*	With successive passes, the data being sorted is passed between the
	 *	buffer passed in as data and the buffer initially allocated to
	 *	bucket[0].  This is done an even number of times, and so in the end,
	 *	it ends up in the original data buffer.
	 */
}


//	As a control experiment, we can't ignore D's built in sort facility.
void builtInSort(uint[] data)
out {
	debug (100) printList(data);
	checkSorted(data);
}
body {
	debug (100) printList(data);
	data.sort;
}