I've been investigating an instance that I ran into while
processing some data, in which multiSort had apparently stalled
whereas two successive sort calls processed the data quickly.
After reducing the code and 20 million points, I've reduced the
problem to this test case:
enum N=1000000;
N.iota.retro.chain(N.only).array.sort();
Apparently, this simple case exhibits worst-case (quadratic)
complexity in our sort implementation. The .retro isn't
necessary - this also exhibits quadratic complexity (athough
the total number of predicate calls is halved), because the
first quickSort cycle will flip the order of elements during
partitioning:
N.iota.chain(0.only).array.sort();
Closer inspection shows that the problem is manifested due to
how the pivot is chosen. Currently, std.algorithm's getPivot
uses the median-of-three method (median value of first, middle
and last element[1]). However, in this case, the median element
will always be the first element - this remains true in all
quickSort iterations.
Here is an illustrated example with N=9. Note that at the same
time as choosing the pivot, getPivot also orders the candidates
(first, middle and last range elements) according to the sort
predicate. After calling getPivot, quickSortImpl moves the
pivot to the end of the range by swapping it with the last
element.
getPivot(0..10)
8,7,6,5,4,3,2,1,0,9 <- getPivot - before swap
9,7,6,5,4,8,2,1,0,3 <- getPivot - after swap
9,7,6,5,4,3,2,1,0,8 <- quickSortImpl - after swap
9,8,6,5,4,3,2,1,0,7 <- quickSortImpl - after partition
getPivot(2..10)
6,5,4,3,2,1,0,7 <- getPivot - before swap
7,5,4,3,6,1,0,2 <- getPivot - after swap
7,5,4,3,2,1,0,6 <- quickSortImpl - after swap
7,6,4,3,2,1,0,5 <- quickSortImpl - after partition
(...)
The algorithm gets stuck on "slopes" in sorted data, which are
not uncommon in real-world data (my data resembled a 1D
heightmap). Although the median-of-three way of picking a pivot
is recommended by some sources (notably Sedgewick), perhaps a
better method would be better suitable for Phobos:
* Many sources recommend using a random element as a pivot.
According to [2], "Randomized quicksort, for any input, it
requires only O(n log n) expected time (averaged over all
choices of pivots)". Also, if it is not possible to predict the
pivot choice, it is impossible to craft worst-case input, which
is a plus from a security point[3]. However, I'm not sure if
making the behavior of std.algorithm's sort nondeterministic is
desirable.
* It looks like many STL implementations use Introsort[4] - a
hybrid sorting algorithm which starts with QuickSort, but
switches to HeapSort if the recursion limit exceeds a threshold.
* I found a recent article[5] which describes an optimal
algorithm of picking a pivot, however it is behind a paywall.
[1]: Apparently the term is also used to refer to choosing
three points *randomly*, instead of first/middle/last:
http://stackoverflow.com/a/164183/21501
[2]:
http://en.wikipedia.org/wiki/Quicksort#Analysis_of_Randomized_quicksort
[3]:
https://www.google.com/search?q=quicksort+worst+case+denial+of+service
[4]: http://en.wikipedia.org/wiki/Introsort
[5]: https://news.ycombinator.com/item?id=6629117
