From 5aa7dc5daf821d1bdacde2fe23523a5406c70e8e Mon Sep 17 00:00:00 2001
From: Russ Cox <rsc@golang.org>
Date: Mon, 17 Nov 2008 11:51:34 -0800
Subject: [PATCH] adopt suggestions from Bentley and McIlroy (SP&E Nov 1993) to
 make qsort more robust:

	* use "ninther" to choose pivot.
	* use three-way partition to avoid quadratic
 	  behavior on all-one-value arrays.

also add tests suggested in that paper.

the immediate cause of the slowness we observed was
in fact none of these: the recursive call was sorting
data[0:m] instead of data[a:m].

also rename package to "sort" to match convention.

R=r,gri
DELTA=358  (255 added, 21 deleted, 82 changed)
OCL=19341
CL=19373
---
 src/lib/sort.go | 127 ++++++++++++++++++++-------
 test/sorting.go | 225 +++++++++++++++++++++++++++++++++++++++++-------
 2 files changed, 293 insertions(+), 59 deletions(-)

diff --git a/src/lib/sort.go b/src/lib/sort.go
index fb5f77f4711..381388223f2 100644
--- a/src/lib/sort.go
+++ b/src/lib/sort.go
@@ -2,7 +2,7 @@
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
-package Sort
+package sort
 
 export type SortInterface interface {
 	len() int;
@@ -10,43 +10,112 @@ export type SortInterface interface {
 	swap(i, j int);
 }
 
-
-func Pivot(data SortInterface, a, b int) int {
-	// if we have at least 10 elements, find a better median
-	// by selecting the median of 3 elements and putting it
-	// at position a
-	if b - a >= 10 {
-		m0 := (a + b) / 2;
-		m1 := a;
-		m2 := b - 1;
-		// bubble sort on 3 elements
-		if data.less(m1, m0) { data.swap(m1, m0); }
-		if data.less(m2, m1) { data.swap(m2, m1); }
-		if data.less(m1, m0) { data.swap(m1, m0); }
-		// "m0 <= m1 <= m2"
+func min(a, b int) int {
+	if a < b {
+		return a;
 	}
-	
-	m := a;
-	for i := a + 1; i < b; i++ {
-		if data.less(i, a) {
-			m++;
-			data.swap(i, m);
+	return b;
+}
+
+// Insertion sort
+func InsertionSort(data SortInterface, a, b int) {
+	for i := a+1; i < b; i++ {
+		for j := i; j > a && data.less(j, j-1); j-- {
+			data.swap(j, j-1);
 		}
 	}
-	data.swap(a, m);
-	
-	return m;
 }
 
+// Quicksort, following Bentley and McIlroy,
+// ``Engineering a Sort Function,'' SP&E November 1993.
 
-func Quicksort(data SortInterface, a, b int) {
-	if a + 1 < b {
-		m := Pivot(data, a, b);
-		Quicksort(data, 0, m);
-		Quicksort(data, m + 1, b);
+// Move the median of the three values data[a], data[b], data[c] into data[a].
+func MedianOfThree(data SortInterface, a, b, c int) {
+	m0 := b;
+	m1 := a;
+	m2 := c;
+
+	// bubble sort on 3 elements
+	if data.less(m1, m0) { data.swap(m1, m0); }
+	if data.less(m2, m1) { data.swap(m2, m1); }
+	if data.less(m1, m0) { data.swap(m1, m0); }
+	// now data[m0] <= data[m1] <= data[m2]
+}
+
+func SwapRange(data SortInterface, a, b, n int) {
+	for i := 0; i < n; i++ {
+		data.swap(a+i, b+i);
 	}
 }
 
+func Pivot(data SortInterface, lo, hi int) (midlo, midhi int) {
+	m := (lo+hi)/2;
+	if hi - lo > 40 {
+		// Tukey's ``Ninther,'' median of three medians of three.
+		s := (hi - lo) / 8;
+		MedianOfThree(data, lo, lo+s, lo+2*s);
+		MedianOfThree(data, m, m-s, m+s);
+		MedianOfThree(data, hi-1, hi-1-s, hi-1-2*s);
+	}
+	MedianOfThree(data, lo, m, hi-1);
+
+	// Invariants are:
+	//	data[lo] = pivot (set up by ChoosePivot)
+	//	data[lo <= i < a] = pivot
+	//	data[a <= i < b] < pivot
+	//	data[b <= i < c] is unexamined
+	//	data[c <= i < d] > pivot
+	//	data[d <= i < hi] = pivot
+	//
+	// Once b meets c, can swap the "= pivot" sections
+	// into the middle of the array.
+	pivot := lo;
+	a, b, c, d := lo+1, lo+1, hi, hi;
+	for b < c {
+		if data.less(b, pivot) {	// data[b] < pivot
+			b++;
+			continue;
+		}
+		if !data.less(pivot, b) {	// data[b] = pivot
+			data.swap(a, b);
+			a++;
+			b++;
+			continue;
+		}
+		if data.less(pivot, c-1) {	// data[c-1] > pivot
+			c--;
+			continue;
+		}
+		if !data.less(c-1, pivot) {	// data[c-1] = pivot
+			data.swap(c-1, d-1);
+			c--;
+			d--;
+			continue;
+		}
+		// data[b] > pivot; data[c-1] < pivot
+		data.swap(b, c-1);
+		b++;
+		c--;
+	}
+
+	n := min(b-a, a-lo);
+	SwapRange(data, lo, b-n, n);
+
+	n = min(hi-d, d-c);
+	SwapRange(data, c, hi-n, n);
+
+	return lo+b-a, hi-(d-c);
+}
+
+func Quicksort(data SortInterface, a, b int) {
+	if b - a > 7 {
+		mlo, mhi := Pivot(data, a, b);
+		Quicksort(data, a, mlo);
+		Quicksort(data, mhi, b);
+	} else if b - a > 1 {
+		InsertionSort(data, a, b);
+	}
+}
 
 export func Sort(data SortInterface) {
 	Quicksort(data, 0, data.len());
diff --git a/test/sorting.go b/test/sorting.go
index ae9dafb7519..ae278141f8f 100644
--- a/test/sorting.go
+++ b/test/sorting.go
@@ -6,13 +6,19 @@
 
 package main
 
-import Sort "sort"
+import (
+	"fmt";
+	"rand";
+	"sort";
+)
+
+func BentleyMcIlroyTests();
 
 func main() {
 	{	data := []int{74, 59, 238, -784, 9845, 959, 905, 0, 0, 42, 7586, -5467984, 7586};
-		a := Sort.IntArray{&data};
-		
-		Sort.Sort(&a);
+		a := sort.IntArray{&data};
+
+		sort.Sort(&a);
 
 		/*
 		for i := 0; i < len(data); i++ {
@@ -20,16 +26,16 @@ func main() {
 		}
 		print("\n");
 		*/
-		
-		if !Sort.IsSorted(&a) {
+
+		if !sort.IsSorted(&a) {
 			panic();
 		}
 	}
 
 	{	data := []float{74.3, 59.0, 238.2, -784.0, 2.3, 9845.768, -959.7485, 905, 7.8, 7.8};
-		a := Sort.FloatArray{&data};
-		
-		Sort.Sort(&a);
+		a := sort.FloatArray{&data};
+
+		sort.Sort(&a);
 
 		/*
 		for i := 0; i < len(data); i++ {
@@ -37,16 +43,16 @@ func main() {
 		}
 		print("\n");
 		*/
-		
-		if !Sort.IsSorted(&a) {
+
+		if !sort.IsSorted(&a) {
 			panic();
 		}
 	}
 
 	{	data := []string{"", "Hello", "foo", "bar", "foo", "f00", "%*&^*&^&", "***"};
-		a := Sort.StringArray{&data};
-		
-		Sort.Sort(&a);
+		a := sort.StringArray{&data};
+
+		sort.Sort(&a);
 
 		/*
 		for i := 0; i < len(data); i++ {
@@ -54,17 +60,17 @@ func main() {
 		}
 		print("\n");
 		*/
-		
-		if !Sort.IsSorted(&a) {
+
+		if !sort.IsSorted(&a) {
 			panic();
 		}
 	}
-	
+
 	// Same tests again, this time using the convenience wrappers
-	
+
 	{	data := []int{74, 59, 238, -784, 9845, 959, 905, 0, 0, 42, 7586, -5467984, 7586};
-		
-		Sort.SortInts(&data);
+
+		sort.SortInts(&data);
 
 		/*
 		for i := 0; i < len(data); i++ {
@@ -72,15 +78,15 @@ func main() {
 		}
 		print("\n");
 		*/
-		
-		if !Sort.IntsAreSorted(&data) {
+
+		if !sort.IntsAreSorted(&data) {
 			panic();
 		}
 	}
 
 	{	data := []float{74.3, 59.0, 238.2, -784.0, 2.3, 9845.768, -959.7485, 905, 7.8, 7.8};
-		
-		Sort.SortFloats(&data);
+
+		sort.SortFloats(&data);
 
 		/*
 		for i := 0; i < len(data); i++ {
@@ -88,15 +94,15 @@ func main() {
 		}
 		print("\n");
 		*/
-		
-		if !Sort.FloatsAreSorted(&data) {
+
+		if !sort.FloatsAreSorted(&data) {
 			panic();
 		}
 	}
 
 	{	data := []string{"", "Hello", "foo", "bar", "foo", "f00", "%*&^*&^&", "***"};
-		
-		Sort.SortStrings(&data);
+
+		sort.SortStrings(&data);
 
 		/*
 		for i := 0; i < len(data); i++ {
@@ -104,9 +110,168 @@ func main() {
 		}
 		print("\n");
 		*/
-		
-		if !Sort.StringsAreSorted(&data) {
+
+		if !sort.StringsAreSorted(&data) {
 			panic();
 		}
 	}
+
+	{
+		data := new([]int, 100000);
+		for i := 0; i < len(data); i++ {
+			data[i] = rand.rand() % 100;
+		}
+		if sort.IntsAreSorted(data) {
+			panic("terrible rand.rand");
+		}
+		sort.SortInts(data);
+		if !sort.IntsAreSorted(data) {
+			panic();
+		}
+	}
+
+	BentleyMcIlroyTests();
+}
+
+const (
+	Sawtooth = iota;
+	Rand;
+	Stagger;
+	Plateau;
+	Shuffle;
+	NDist;
+)
+
+const (
+	Copy = iota;
+	Reverse;
+	ReverseFirstHalf;
+	ReverseSecondHalf;
+	Sort;
+	Dither;
+	NMode;
+);
+
+type TestingData struct {
+	data *[]int;
+	maxswap int;	// number of swaps allowed
+	nswap int;
+}
+
+func (d *TestingData) len() int { return len(d.data); }
+func (d *TestingData) less(i, j int) bool { return d.data[i] < d.data[j]; }
+func (d *TestingData) swap(i, j int) {
+	if d.nswap >= d.maxswap {
+		panicln("used", d.nswap, "swaps sorting", len(d.data), "array");
+	}
+	d.nswap++;
+	d.data[i], d.data[j] = d.data[j], d.data[i];
+}
+
+func Lg(n int) int {
+	i := 0;
+	for 1<<uint(i) < n {
+		i++;
+	}
+	return i;
+}
+
+func Min(a, b int) int {
+	if a < b {
+		return a;
+	}
+	return b;
+}
+
+func SortIntsTest(mode int, data, x *[]int) {
+	switch mode {
+	case Copy:
+		for i := 0; i < len(data); i++ {
+			x[i] = data[i];
+		}
+	case Reverse:
+		for i := 0; i < len(data); i++ {
+			x[i] = data[len(data)-i-1];
+		}
+	case ReverseFirstHalf:
+		n := len(data)/2;
+		for i := 0; i < n; i++ {
+			x[i] = data[n-i-1];
+		}
+		for i := n; i < len(data); i++ {
+			x[i] = data[i];
+		}
+	case ReverseSecondHalf:
+		n := len(data)/2;
+		for i := 0; i < n; i++ {
+			x[i] = data[i];
+		}
+		for i := n; i < len(data); i++ {
+			x[i] = data[len(data)-(i-n)-1];
+		}
+	case Sort:
+		for i := 0; i < len(data); i++ {
+			x[i] = data[i];
+		}
+		// sort.SortInts is known to be correct
+		// because mode Sort runs after mode Copy.
+		sort.SortInts(x[0:len(data)]);
+	case Dither:
+		for i := 0; i < len(data); i++ {
+			x[i] = data[i] + i%5;
+		}
+	}
+	d := &TestingData{x[0:len(data)], len(data)*Lg(len(data))*12/10, 0};
+	sort.Sort(d);
+
+	// If we were testing C qsort, we'd have to make a copy
+	// of the array and sort it ourselves and then compare
+	// x against it, to ensure that qsort was only permuting
+	// the data, not (for example) overwriting it with zeros.
+	//
+	// In go, we don't have to be so paranoid: since the only
+	// mutating method sort.Sort can call is TestingData.swap,
+	// it suffices here just to check that the final array is sorted.
+	if !sort.IntsAreSorted(x[0:len(data)]) {
+		panicln("incorrect sort");
+	}
+}
+
+func BentleyMcIlroyTests() {
+	sizes := []int{100, 1023, 1024, 1025};
+	var x, tmp [1025]int;
+	for ni := 0; ni < len(sizes); ni++ {
+		n := sizes[ni];
+		for m := 1; m < 2*n; m *= 2 {
+			for dist := 0; dist < NDist; dist++ {
+				j := 0;
+				k := 1;
+				for i := 0; i < n; i++ {
+					switch dist {
+					case Sawtooth:
+						x[i] = i % m;
+					case Rand:
+						x[i] = rand.rand() % m;
+					case Stagger:
+						x[i] = (i*m + i) % n;
+					case Plateau:
+						x[i] = Min(i, m);
+					case Shuffle:
+						if rand.rand() % m != 0 {
+							j += 2;
+							x[i] = j;
+						} else {
+							k += 2;
+							x[i] = k;
+						}
+					}
+				}
+				data := (&x)[0:n];
+				for i := 0; i < NMode; i++ {
+					SortIntsTest(i, data, &tmp);
+				}
+			}
+		}
+	}
 }
+