Idea: ska_sort based nth_element

[JobLeonard]

So I came across a series of blog posts about finding the median value of an array (which, really, is just a special version of the selection algorithm). The last post contains a very simple and very fast technique: use counting sort, but after the counting pass, don't even bother swapping. Just look at which bucket falls in the middle of the array. This must be the median. His implementation of this "counting median" is pretty fast:

Then he laments:

Of course, we get a speed-up because we exploit a special case, where the number of different values is small compared to the length of the list. But hey, why not?

Oh... but is it really such a limited special case? ;)

I bet you have already guessed what I have in mind: take ska_sort, but only recurse on the bucket with the range that contains the median value (it also has to pass an off-set for where the middle is in the sub-range). If you do that, it actually meets all the requirements of std::nth_element:

Rearranges the elements in the range [first,last), in such a way that the element at the nth position is the element that would be in that position in a sorted sequence.

The other elements are left without any specific order, except that none of the elements preceding nth are greater than it, and none of the elements following it are less.

With or without that last optimization this should be blazing fast, even for long radix keys: for uniformly distributed values, the number of elements can be expected to shrink by 256 each time the function recurses. The worst-case input would be an input where the bucket size shrinks a minimal amount (I guess just by one element) - that would just degenerate into ska_sort, so still perform pretty well.

Since most implementations of nth_element use quick_select, a "radix select" like this should probably beat it easily. Even better, it would require relatively little changes to the existing ska_sort code!

I also thought of an even simpler and likely faster method: instead of swapping everything, only swap the values in the median range. To make this easier, swap them to the front. In pseudo-code:

unsigned int i = 0;
// skip all the values of the median range already at the front
while(array[i] == median_range){ i++; }
unsigned int j = i++;
while(i < array_size && j < median_range_size){
  while(array[i] == median_val && j < median_range_size){
    swap(j++, i++);
  }
  i++;
}

Then we recurse over the median range, with an off-set to indicate where the median falls in the sub-range. And like the counting median approach above, we can skip the swapping on the deepest level of recursion, and only do the counting part.

This should be even faster than the previous version, because it reduces the nr of swaps to the bare minimum. The obvious downside is that this would completely scramble the array order.