C++ mergesortAsymmetric source code

Here's my current draft of mergesortAsymmetric(). But I'm suspicious there's something wrong with it somewhere, because when I tuned it for space it became suspiciously fast. So far, it seems a little faster, on average, than quickSortClassic(), for sorting integers. I don't see how that can be right! (It does more element moves and more comparisons than quicksort or standard binary mergesort - and I haven't fixed up either of the low-level sort routines to use pre-increment instead of post-increment). There must be a mistake buried in there somewhere, but I can't see what it is.

template <class T>void insertionSort(T *src, int count)
{
  //
  //Notes:  1.insertionSort is a "workhorse" sorting routine.  Many
  //      other sorting routines use it.  quickSort and Combsort
  //      may use it to "finish off the job".  Mergesorts may use it
  //      to create initial runs.
  //      2.insertionSort has a very low "set-up" cost, and for
  //      small enough N, is more efficient than quickSort.
  //      3.The following assumes that ((void*)src-sizeof(src[0]) >= 0),
  //          which may not be true if sizeof(src[0]) is large enough!
  //
  T v;
  T *stopAt  = src+count;  //pointer to location after last element to sort
  T *sorting;          //pointer to location after last element known
                //to be in sort order w.r.t. elements to its left
  T *scanning;        //used for finding where the element at *sorting
                //should be placed.

  for (sorting=src+1 ; sorting<stopAt; sorting++)
  {
    v=*sorting;
    for (scanning=sorting-1; src<=scanning && v<*scanning; scanning--)
    {
        scanning[1]=scanning[0];
    }
    scanning[1]=v;
  }
}

template <class T> int insertionSortExternal(T *input, int count, T *output)
{
  //note: it is assumed that there is no overlap between output and input
  //    (or in C notation, that this...
  //      !(input<=output && output<input+count) &&
  //      !(output<=input && input<output+count)
  //    ...holds true).
  //
  T *endOutput = output+count;
  T *sortedBoundary;
  T *scanBack;
  *output = *input++;        //copy over first element
  for (sortedBoundary=output+1;sortedBoundary<endOutput;sortedBoundary++)
  {
    for (scanBack=sortedBoundary-1;scanBack>=output;scanBack--)
      if (*scanBack<=*input) break; else scanBack[1]=*scanBack;
    scanBack[1]=*input++;
  }
  return 0;
}

template <class T> inline int indexOfFirstElementInBlockGreaterThan(T *block, int count, T &v)
{
  int loIndex = 0;
  int hiIndex = count;
  while (loIndex<hiIndex)
  {
    int tryIndex = ( loIndex + hiIndex ) / 2;
    if ( v < block[tryIndex]  )
      hiIndex  = tryIndex;
    else
      loIndex = tryIndex + 1;
  }
  return hiIndex;
}

//Date        By Change
//=========== == ======
//04-Jan-2010 JB Draft
//11-Jan-2010 JB Revisions to AsymmetricMergesorter, to reduce the space requirement 
//               from NR  (N=record count, R=record size) 
//               to   NPR (where P=ratio of smaller
//                         sublist to list size during each merge).
//
//               There are two sort routines, sortRight and sortLeft.
//               sortRight: 1. calls itself to sort the right-hand (n-np) records,
//                             back to the same location
//                          2. calls sortLeft to sort the left-hand np records into
//                             the np-record workarea
//                          3. merges the left-hand data from the workarea, with
//                             the right-hand data.
//               sortLeft:  1. sorts left hand of destination to source
//                          2. sorts right hand of destination to source
//                          3. merges left and right source to destination                          
//
//               There's an unexpected benefit: performance is considerably *better*
//               

template<class T> void mergeAsymmetricRadix2External(T* a, T* aStop, T* b, T* bStop, T* dest)
{
  //Assumes:  aStop-a is considerably less than bStop-b
  //Note:     We want to merge until we run out of elements in a.
  //          But, if elements in b will run out first, we can't do that safely.
  //

  if (aStop[-1]<=bStop[-1])//the easy case.  a runs out first
  {
    for (;a<aStop;++dest,++a)
    {
      for (;*b<*a;++dest,++b)
      {
        *dest=*b;
      }
      *dest=*a;
    }
    for (;b<bStop;++dest,++b)
    {
      *dest=*b;
    }
  }
  else //the harder case: b's will run out first.  but we to check for a's running out *first*
  {
    //
    //find a portion of the a's, left(a) that *will* run out first,
    //and merge left(a) and b until left(a) runs out.
    //
    T* aChop = a + indexOfFirstElementInBlockGreaterThan(a, aStop-a, bStop[-1]);
    for (;a<aChop;++dest,++a)
    {
      for (;*b<*a;++dest,++b)
      {
        *dest=*b;
      }
      *dest=*a;
    }
    for (;b<bStop;++dest,++b)
    {
      for (;*a<=*b;++dest,++a)
      {
        *dest=*a;
      }
      *dest = *b;
    }
    for (;a<aStop;++dest,++a)
    {
      *dest=*a;
    }
  }
}

template <class T> class AsymmetricMergesorter
{
private:
  T*  m_base;
  T*  m_workArea;
  int m_workCount;
  int m_count;
  int m_cutOff;
  int m_numerator;
  int m_denominator;
public:
  AsymmetricMergesorter(T* base, int count, int cutOff, int numerator, int denominator)
    : m_base(base)
    , m_count(count)
    , m_cutOff(cutOff)
    , m_numerator(numerator)
    , m_denominator(denominator)
  {
    m_workCount = (long)(long)m_count * m_numerator / m_denominator ;
    m_workArea  = new T [ m_workCount ]; //NPR space
  }
  ~AsymmetricMergesorter()
  {
    delete [] m_workArea;
  }
  void sortLeft(T *source, int count, T* dest, bool bSourceIsInput)
  {
    if (count<m_cutOff)
    {
      if (bSourceIsInput)
        insertionSortExternal(source, count, dest);
      else
        insertionSort(dest, count);
    }
    else
    {
      int leftCount = ( (long)(long)count * m_numerator / m_denominator );
      sortLeft( dest,             leftCount,       source,         !bSourceIsInput );
      sortLeft( source+leftCount, count-leftCount, dest+leftCount,  bSourceIsInput);
      mergeAsymmetricRadix2External( source, source+leftCount, dest+leftCount, dest+count, dest);
    }
  }
  void sortRight(T* source, int count)
  {
    if (count<m_cutOff)
    {  
      insertionSort(source, count);
    }
    else
    {
      int leftCount = ( (long)(long)count * m_numerator / m_denominator );    
      sortRight(source+leftCount, count-leftCount);
      sortLeft (source,           leftCount,        m_workArea,  true);
      mergeAsymmetricRadix2External(m_workArea, m_workArea+leftCount, source+leftCount, source+count, source);
    }
  }
  void sort()
  {
    sortRight(m_base, m_count);
  }
};

template<class T> void mergesortAsymmetricRadix2(T *a, int count,int cutOff=32,int numerator=1, int denominator=4)
{
  AsymmetricMergesorter<T> s(a, count, cutOff, numerator, denominator);
  s.sort();
}

Tuning Mergesorts for better branch prediction

Yesterday, I decided to play around with a forecasting binary mergesort.

Some background: way back in about 1999, when I was first experimenting with ternary mergesorts, I had noticed that the best "division" of the list into sublists was not a 1:1:1 ratio. On the Pentium II that I had then, the best "division" (for sorting 32-bitintegers, anyway) seemed to be 3:4:5. I filed that fact away but didn't think about it again until yesterday.

It occurred to me that, during a merge, it's possible to rig the merge so that, in the main loop, the smaller of the two lists runs out first (simple enough: if the smaller list runs out last, find where the last element in the larger list would go; which makes it easy to determine the latest element of the smaller list that goes before the last element in the larger list). The point is that, if the smaller of the lists consists of a proportion p, 0<p<0.5, of the n elements being merged, then the number of boundary checks will be np. Furthermore the actual element comparisons will be more predictable, because, on average, the element from the smaller list will be be less with probability p. There would be, I reasoned, some value of p below 0.5 for which performance would be better, because, asymptoticially, the sort would require

moves(m): -1/(plog2(p) + (1-p)log2(1-p)) .n.log2(n)
comparisons: O(n) less
boundary checks: pm

My guess was that the best value of p (for sorting integers, anyway) would be about 0.4, with about a 2% performance improvement (because I thought that boundary checks in a non-forecasting "stock" mergesort contribute about 24% of the total running time). But it wasn't so. The best value of p appears to be about 0.1875, with a 5% performance improvement. That shocked me, because plugging p=0.1875 into the formulae above yields...

moves: 1.436 n.log2(n) (contrast: the coefficient is 1, for p=0.5)
checks: 0.287 n.log2(n) (the coefficient is 0.5, for p=0.5)

The only plausible idea I could come up with to explain that was: the lower p, the more predictable the results of the comparisons (asymptotically, p is the branch misprediction rate).

Unfortunately, the best value for p depends on the element type.

I also looked at asymmetric quicksorts (e.g. take the 2nd of 4 elements from a sample, deliberately choosing a pivot value that won't divide the input evenly, but if there's any advantage at all - which I doubt - it is very small). I'll post the (single-threaded) source codeand some performance plots, for the asymmetric mergesort, in my next post. I also want to have another go at asymmetric ternary mergesort (it seems to me that the ratio 1:2:4 ought to do well on average).

I haven't given any thought yet to whether there are similar tweaks available for shellsort and combsort (trading some extra "work" for better branch prediction).

Sorting with less data movement than selection sort

Until recently, the wikipedia article on selection sort claimed it required less data movement than any other sorting algorithm. Luckily, that seems to have been just-about corrected; there's now an article (http://en.wikipedia.org/wiki/Cycle_sort) outlining a cycle sort algorithm that does "far fewer" array writes ("far fewer" is what Wikipedia says; it's about half as many, on average). I say, "just about corrected", because the cycle sort outlined there is still doing too many element moves (about as many as selection sort). The lines that go...

array[pos], item = item, array[pos]

still require three moves (because that's an exchange). If, however, you've got two item variables, item1 and item2, you can instead go...

array[pos], item2 = item1, array[pos]

half the time, and

array[pos], item1 = item2, array[pos]

the other half of the time, switching the sense of item1 and item2 each time (there's a very similar hack available in heapsort that reduces the number of index assignments). That change will reduce the average number of element moves by 1/3, to 2(n-f), where n is the number of elements in the array, and f is the number of elements already in their final position.

But that's still too many element moves. Applying a random permutation to n elements in place should require n+c-2*f element moves, where c is the number of cycles (including one-element cycles) - H(n) on average - and f is the number of elements already placed correctly (the number of one-element cycles), which yields: minimum 0, average n+H(n)-1, maximum 3n/2 moves. Even the "two item" version of cycle sort still requires on average almost twice as much data movement as it should.

I haven't been able to figure out an algorithm with a constant space overhead that achieves that. Maybe there isn't one. But I can see how to get closer to it, by maintaining more index variables. If you've got i index variables, the length of the current cycle can be reduced by i-1, at the cost of i element moves. The problem with that is handling duplicate values (it's bad enough in cycle sort but if you've got more index variables it gets much worse, because you have to skip over duplicates that will be moved into place when the cycle is "unraveled", but haven't actually been moved into place yet).

If the output isn't in the same location as the input, then n element moves will be sufficient (using a comparison counting sort, which will also be stable).

How to make Quicksort suck: use one-handed partitioning

While I was away (back in late 2009), somebody posted a Dual Pivot Quicksort, that advertised (and delivered) a 20% reduction in the number of exchanges, compared to Quicksort, as defined in wikipedia's article (http://en.wikipedia.org/wiki/Quicksort ). That sounds good, doesn't it? But...

Wikipedia insists on posting a quicksort with a partitioning routine that is downright lousy. The number of exchanges required during a quicksort of a randomly ordered array of n distinct elements shouldn't be ln(n)*n on average. It should be ln(n)*n/3. In other words, Dual Pivot Quicksort still does 12/5 = 2.4 times as many exchanges as are actually necessary.

Here's the pseudo-code for partition() in the Wiki article (verbatim) (which I think of as a "partition with one hand tied behind its back", or a one-handed partition):

  function partition(array, left, right, pivotIndex)
pivotValue := array[pivotIndex]
swap array[pivotIndex] and array[right] // Move pivot to end
storeIndex := left
for i  from  left to right - 1 // left ≤ i <>
    if array[i] ≤ pivotValue
        swap array[i] and array[storeIndex]
        storeIndex := storeIndex + 1
swap array[storeIndex] and array[right] // Move pivot to its final place
return storeIndex

That pseudocode is very pretty, very short, and very stupid. Better partitioning routines have been around for decades (google for Sedgewick Partition: certainly since the 1970s), and I wouldn't mind betting that the Quicksorts that beat Mergesort (for sorting arrays of value rather than reference types, anyway) actually use those. The pseudocode for the partition routine should look more like this (this is a "two-handed" version):

function partition(array, left, right, pivotIndex)
pivotValue := array[pivotIndex]
swap array[pivotIndex] and array[left] // Move pivot to middle
i = left+1
j = right
while i <= j and array[i] < pivotValue
    i = i + 1
while i < j
    while array[i] < pivotValue //search from left                    
        i = i + 1           
    while pivotValue < array[j] //search from right                    
        j = j - 1     
    swap array[i] and array[j]     
    i = i + 1     
    j = j - 1 
 if left < j
    swap array[left] and array[j]     
 return j 

The first while loop is there to handle the special case that there's no element in the remainder of the array that is not less than the pivot value. Otherwise it would be necessary to check that i is less than j when searching from the left. Equal elements get exchanged (that's so we don't have to check that i is less than j during the search from the right, and it also ensures that inputs containing large numbers of duplicates don't cause problems). The last if statement is necessary because there might be no other element greater than or equal to pivotValue.

It's more complicated, yes, but the extra complexity is worth it! The second algorithm "burns the candle at both ends", searching from the left for values larger than (or equal to) the pivot, and from the right for values smaller than (or equal to) the pivot, and exchanges each such pair that it finds. The number of exchanges is equal to the number of values that should be, but aren't yet, in the smaller of the two partitions.

Contrast this with the algorithm from Wikipedia's article, which performs one exchange for each value that should be in the left-hand partition. On average, to partition N elements, the "two-handed" version needs N/8 exchanges rather than the N/2 required by the "one-handed" version: four times fewer.

However, switching to "two-handed" partitioning will not reduce the expected number of exchanges for the entire quicksort by a factor of 4. It actually reduces it by a factor of only 3.

If you're ready to "invest" in two more index variables, you can switch to a routine that does "rotations" rather than exchanges, and will partition an array of N elements in, on average, N/4 moves. It's theoretically better (it reduces the number of moves by a third, on average), but it's much messier and only marginally faster.

Dzmitry Huba's parallel mergesort

Wow! No post in nearly eighteen month. My daughter took over the laptop. And the dog ate my homework.

Recently I spotted an interesting post, with C# code, at...

http://dzmitryhuba.blogspot.com/2010/10/parallel-merge-sort.html

...and the code was pretty damned good (nice how C# delegates make it all a lot cleaner; in my C++ code I have had to write a bunch of hideously ugly template functions to fake delegates; perhaps I should be coding in C# instead). A few minor, minor things:

• The main reasons that Dzmitry's parallel mergesort outperformed the parallel quicksort at http://msdn.microsoft.com/en-us/library/ff963551.aspx were probably that (a) the Partition() routine is an inefficient "one-handed" implementation [see next post], and (b) the parallel quicksort did not parallelize the partitioning step. I've got some code that does that (badly) lying around somewhere which I'll paste into a later post (if I can find it). If the partitioning is parallelized, and you're using "two-handed" partitioning, and you're sorting integers, parallel quicksort should win.
  
  
• The load balancing problem that Dzmitry mentioned is largely caused by the use of BinarySearch in his ParallelMerge(). The problem is that using the median of one of the ranges and binary searching for that value in the other range doesn't divide the merge evenly enough, particularly if the input was ordered or mostly ordered. If you are merging two ranges of size n, you need to search, instead, for an index, i, such that the first i elements of the left subrange will be merged with the first n-i elements of the right subrange. That will always handle an (n,n) merge by scheduling two same-sized merges: one of (i,n-i), one of (n-i,i). There's no risk of it scheduling an (n/2,n) and an (n/2,0) merge.
  
  Unfortunately I was lazy, and the code I've got that does that is wired in to the merge routines themselves (and furthermore it is generalized for splitting the load across P cores, where P is not necessarily a power of two, and takes P as a parameter, which makes it much harder to see what it is doing). I'll have to do a rewrite.
  
  

Zonesort mostly parallelized, but... it's a pyrrhic victory

Well... I have finally parallelized Zonesort (just about, I still need to do some extra work to deal with degenerate - that is, almost sorted or almost reverse-sorted input - cases properly). But there's a problem. It needs better load balancing. And I can't see a good way to do that. A single gantt chart illustrates the symptoms neatly.

This is a graph of CPU utilization (color, red is 100%, white is 0%) by time (across) by worker thread number (down), for a four thread Zonesort using 3-way merging (except for merging the last four zone lists), while it is sorting 9,395,322 distinct integers in about 0.32 seconds.

I've finally parallelized the last few merge operations (though I had to increase the space requirements by a factor of sqrt((M+1)/M) to achieve that), as well as the zone shuffling operation, but there's a load balancing problem earlier on in the sort, making - in this case, which is typical - the overall running time about 6% worse than it should be. If I could fix that it would still be considerably slower (by at least 10%) than a stock-standard 3-way mergesort on the same input. So I don't think I'll bother.

The other sorting routines I've ever parallelized never had significant per thread storage requirements, the way this one seems to. I can't share available zone lists between threads easily (and, no, serializing access to a global available zone list doesn't strike me as a good idea at all; that would mean altogether too much synchronization), so for a P-threaded sort I'm dividing the input into P pieces and sorting each piece independently. But some cores - and some threads - run faster than others, and my heroic (for which, read: stupid) assumption, that those P independent sorts would each run in about the same time, is clearly wrong.

I'm already using "barriers" in the last few merges and for the final zone shuffle (the time wasted by those barriers is visible in the gantt chart, and a barrier always means an idle thread), but at least I'm not doing any explicit thread control in the sort, and I'm not throwing away much extra memory (an important consideration when the whole point of this sorting routine is to save space compared to a standard merge sort) to co-ordinate threads running in parallel (so far I've been using up only about 200*P bytes on per-thread storage).

I'll post the parallel zoneSort once I've ironed out the bugs that stop it working for nearly-in-order or nearly-reverse-order inputs.

Zonesort - Yet another variant of mergesort

I've been on holiday in Thailand, so I haven't been posting. While I was in Thailand, in between looking at orchids, eating fantastic seafood, and getting short-changed everywhere I went, I worked a little on a variant of mergesort ("Zoned Mergesort") that requires a storage overhead of approximately

2 * sqrt( N.R.I.M.P )

where N is the number of records, R is the size of each record, I is the size of an integer, M is the order of merge, and P is the number of concurrent threads (mergesort overhead is usually something like N.R / 2 which is a lot more, for large enough N). It's an old idea of mine: I first wrote an implementation in 2004 (though only for M=2, P=1) but shelved it because it was considerably (typically 20%) slower than a standard mergesort. But it turned out that some minor rewrites to better exploit forecasting let it compete on even terms with the performance of standard mergesort, with P=1 and M=2. So I extended it to M=3, where it did even better. But I still haven't parallelized it properly (and that will take time, because it's a fairly complicated sorting routine).

Zoned Mergesort was derived, in a round-about way, from an idea of
Alexander Kronrod (see Knuth's Art of Programming, volume 3): merging dividing arrays of N records into approximately sqrt(N) "zones" of approximately sqrt(N) records. It works like ths:

1. Divide the input array into X zones of Y records each (where X is approximately sqrt(N*I/R/M/P))


2. Allocate storage for M*P*Y additional records.


3. Allocate storage for (M*P+X) integers. These will be used to maintain linked lists of sorted zones.


4. Sort each zone, using the additional zones as a scratch area, and put each zone's index into a one element-list


5. Start P separate lists, each of M zones, chosen from the M*P additional zones (these are "available zone" lists


6. Assemble linked lists of sorted zones by merging and relinking (I'll explain this step in more detail later)


7. When you get to one linked list of zones, permute the order of the zones so that the entire array is in order

The time overhead for steps 6 and 7 turns out not to be a serious problem. The linked list manipulation in step 6 has O(X) overhead, and can be pretty much ignored when N is large. The permutation in step 7 is rather more expensive but is still O(N) on average. In practice, when P=1 at least, for large N, the permutation cost is outweighed by the improved cache hit rate on writes (almost all writes end up being cache hits). The cache hit rate is better than that of a standard mergesort due to how the zone list merging works (which I'll outline for M=2):

1. Remove a destination zone from the current thread's list of available zones


2. Merge records from the zone at the head of the left list, and the zone at the head of the right list, into the destination zone, until either the destination zone is exhausted, or the left zone is exhausted, or the right zone is exhausted


3. If the destination zone is exhausted, add it to the tail of the output list, and remove a zone from the available zone list, to be the current destination zone


4. If the left hand zone is exhausted, move on the left zone list, so that the next zone in the list becomes the left input zone, and put the just-exhausted zone into the available zone list


5. As for step 4, but replacing "left" with "right"


6. If neither the left zone list or the right zone list is exhausted, go back to step 2


7. Copy records from the remaining zones in the input that hasn't been exhausted, until it runs out (allocating new destination zones and releasing exhausted input zones as above).

The cache benefit comes from the fact that, during most of the zone list merging phase, the current destination zone was almost always in use, recently, as an input zone. It turns out that you may have to read X.M-X+1 records before the first input zone is exhausted, which is why there need to be M "extra" zones in the initial available zone list for each thread.

There's some extra fiddling around to take account of the fact that the last zone will usually be smaller than the others (since N is unlikely to zero modulo Y), but nothing serious.

The big problem, though, is that parallelizing the last few zone list merges is painfully difficult. I've only got part of the way through it and I think it'll take me about six hours solid work to write even a clunky fully-parallel version. No fun at all!

Daily "It Ain't So" - wikipiedia's Mergesort article

For some reason blogger has decided this was published 29-Jun-2009. Which ain't so. It's actually 03-Jan-2011 as I type this!

Today's It-Ain't-So is for wikipedia's mergesort article (as it is today, not as it might be down the track).

• A "mergesort" using lists is not an in-place sort. Either the pointer is part of the record - in which case the list mergesort isn't sorting elements, it's modifying them, or the pointer is not part of the record, in which case mergesort isn't in-place.
  [that's how the ancient Greeks would have argued it, and they would have been right]
  
  
• It is not quite true that "mergesorting" lists is, in performance terms, no better than heapsort. It is true, if swapping values between nodes is not allowed (for a randomly-ordered input, as the sort proceeds the order of the addresses becomes increasingly jumbled up, leading to cache misses). But, if swapping values between nodes is allowed, and the input was initially of in-address-order nodes adjacent in memory, it's possible (and not very difficult) to rewrite the nodes in the sorted sublists after every, say, 8 levels of merging, avoiding the cache-miss problem, and beating heapsort comfortably, for large n. I've got some code for doing that buried somewhere (it's about ten years old and I haven't used it in about eight years, and I don't remember if it "cheated" and used additional storage for the rewrite. Maybe it did).
  
  
• If you can use "link fields" for mergesort (and it's still a mergesort?) then you can also use them for quicksort (and hence, quicksort can be stable).
  
  
• Mergesorting linked lists does not require an O(1) space overhead: it requires O(logn) space to record (at each of the O(logn) merging levels) which node is the head of the first sorted sublist, while the second sublist at the same level is being sorted, unless swapping values between nodes is allowed, or adjacent sorted sublists are "taped together" and you're happy to spend O(nlogn) time finding the start of the right hand sublist before each merge begins. I've never seen either of those techniques used; all the linked list mergesorts I've ever seen had an plog2(n) space overhead (on top of the np space overhead for the linkage fields), where p is the size of a link.
  
  
• While it's true that heapsort won't work on linked lists, there is a close analogue of heapsort that works on priority trees. As I understand it, heapsort was originally developed from a priority-tree algorithm.
  
  
• While it's technically true to say that a self-balancing binary search tree is better than an on-line mergesort, when "the received pieces are small compared to the sorted list", there's an invalid assumption wired into that statement. Why would there be one sorted list, rather than ~O(logm) sorted sublists, where m is the number of elements that have been "fed" to the on-line sort routine so far, with the (i+1t)h sublist at most half the length of the ith? That's the explicit stack used in a straight mergesort co-routine.

And now for the discussion page that hangs off the article...

• Where I said that Quicksort's real advantage over Mergesort comes down to the amount of record movement (July 2005), I was, well, wrong. I should have qualified that statement a lot more. I should have said that, a Quicksort with a two-handed partitioning subroutine will beat a (binary) Mergesort (even one tuned to do less boundary checks than Quicksort) if and only if element moves are significantly more expensive than element comparisons.
  

• Asturneresquire (August 2007) is wrong: Quicksort does not make significantly better use of the cache than Mergesort. It might at first seem that mergesort will put more of a load on the cache because it is accessing more memory, but it isn't that simple: mergesort is accessing each address fewer times on average. The cache (on, say, an I7) easily keeps up with both Mergesort and Quicksort. It isn't even "cracking a sweat".

Heapsort - Floyd's Trick - part 2 - Sorting Strings

An aside: I would have posted this one sooner, but I was running into an unpleasant problem due to a bug in a do-it-myself string class. Yes, the world has too many string classes - it probably already had one too many when the second one was written, perhaps even when the first was written - and mine was as buggy, or more so, than most. I know that writing my own string classes is stupid but I loathe std::string. It's soooo bloated, and with methods, not free functions. But still, that doesn't justify my writing my own string classes. Thinking it does, even temporarily, is a good example of the tu qoque phallacy. I have repented now! But the string class is already written. Someday, I'll compound my sin by posting some string classes too.

First, some preliminary comments about sorting strings. In later posts I shall be writing some posts about sort algorithms tailored specifically for sorting strings (e.g. Burstsort, in its various flavours, and B-Tree sorts), where prefix-extraction, tries, and the like, will come into play. But for now, the focus will be on Heapsorts.

Two-way mergesort is generally considered to be the reference record-comparison array sort for sorting strings. That's because, when sorting pointers to strings (which is often what happens in C++ and also, under the hood, in languages like Java), comparisons tend to be relatively expensive compared to record copies or exchanges. And fair enough too. String pointers may also be a good proxy for pointers to some more complex types and classes (when multiple members of those types/classes are being used as the sort key).

Asymptotically, a two-way top-down Mergesort with a randomly ordered input where all the values are distinct, requires N*log(N)/log(2) comparisons, on average (and the best case is closer to half that).

When sorting string pointers, with a single-threaded record comparison sort, the second most important consideration is minimizing the number of comparisons. Minimizing the number of cache misses is actually rather more important, and later on that consideration will lead to some bizarre sorting algorithms. But not yet.

I covered how Floyd's trick is actually implemented in an earlier post. As I remarked there, it reduces the asymptotic average number of comparisons required for a radix M heapsort from M.log(N)/log(M) to (M-1).log(N)/log(M). Using those formulae, and assuming that comparisons are what count, and using the performance, X, of mergesort as the reference, we can try to predict sort routine performance:

• For Radix-2 Heapsort with Floyd's trick: 1.0*X


• For classic Quicksort: 1/2/ln(2) ~= 0.72*X


• For Radix-2 Heapsort without Floyd's trick: 0.5*X


• For Radix-4 Heapsort with Floyd's trick: 2/3 ~= 0.67*X


• For Radix-4 Heapsort without Floyd's trick: 0.5*X

But it ain't so!





As you can see, from the two graphs, merely counting comparisons doesn't lead to a good prediction for the running time. Counting comparisons and cache misses works better. Even rough approximations of the cache miss count lead to reasonably good predictions for the running time.

• For most of the comparisons that take place during a mergesort, one of the records being compared was compared immediately before-hand. The access pattern in the array itself is predictable, and, asymptotically, half of the string look-ups will be cache hits.


• Quicksort likewise (though Quicksort does even better than I would have
  expected, and I don't have any explanation for that yet)


• String look-ups near enough to the top of the heap are cache hits. When the element being sifted down into the heap is being compared "on the way down" both string lookups are cache hits. Given that cache hits are what count, these comparisons are much cheaper than the others, low in the heap. That's the reason that Floyd's trick doesn't improve the performance of radix-2 heapsort much. On the other hand, the access pattern in the array itself is not predictable, so array element look-ups also result in cache misses.

Extra Horsepower

The graph shows the relative performance (for sorting an array of 10 million integers) for Quicksort, versus 2-way, 3-way, and 4-way mergesorts that parallelize the last (and only the last) merge operation at each level of merging.

Before parallelizing 3- and 4-way merges (yuu-uuk), I wondered... would it be worth parallelizing all the 2-way-merge operations, not just the last one at each level of merging (and using a StackDispatcher too)...? The answer turned out to be... no. It's very slightly slower (possibly because the extra synchronization overhead, plus the binary-search overhead for finding where to split up the inputs when dividing up the processing, outweighs the cache benefit, if any).

I also figured it might be a good idea to *check* that stability is preserved. Luckily that's pretty easy, with template sort routines. I declared a type, Horse, generated arrays of Horse (see below) (by setting the value members, randomly permuting, and then writing the starting index of each element into the satellite members), and checked that stability was preserved after sorting.

template<class V, class S> class ValueAndSatellite
{
public:
  V value;
  S satellite;
  ValueAndSatellite()         { }
  ValueAndSatellite(V v, S s) { value=v; satellite=s; }
  ~ValueAndSatellite()        { }
  bool operator<  (const ValueAndSatellite<V,S> &rhs) { return value<rhs.value; }
  bool operator<= (const ValueAndSatellite<V,S> &rhs) { return value<=rhs.value; }
  bool operator== (const ValueAndSatellite<V,S> &rhs) { return value==rhs.value; }
};
typedef ValueAndSatellite<int, int> Horse;

Stability was preserved, for the 2- and 3-way mergesorts, but not for the 4-way mergesort. I haven't yet figured out how the 4-way mergesort went wrong. I'll have to post a correction when I find out.

Incidentally, with an array of Horse (each Horse is twice the size of an int, but has approximately the same comparison cost), Quicksort beat two-way, three-way, and even four-way merging comfortably. Even if all of the value members of the Horses in the array were distinct, Quicksort was 11% faster than a four-way merge.






For 3- and 4-way merging I've opted not to worry so much about trying to "split evenly in the presence of large numbers of duplicates". It's too hard to think through all the cases. Instead, I... just hope the splits are about right ("on average, if there aren't too many duplicates, they'll be okay"). And, at least for sorting integers on 4 cores,

• my 3-way mergesort beats Quicksort comfortably. But not by the 8% I had expected. Only by about 5%.


• my 4-way mergesort beats Quicksort by about 9% (I had expected 8%). But then, it isn't stable (ouch!)

Here's my current code (less classes, functions, and so forth that have appeared in previous posts)...

template<class T> void mergeForecastBackwardsRadix2External(T *a, T* aStop, T* b, T* bStop, T* dest)
{
  if (aStop[1]<=bStop[1])
  {
    //the "b"s will run out before the "a"s.
    for (;b>bStop;*dest--=*b--)
    {
      for (;*b<*a;*dest--=*a--);
    }
    for (;a>aStop;*dest--=*a--);
  }
  else
  {
    for (;a>aStop;*dest--=*a--)
    {
      for (;*a<=*b;*dest--=*b--);
    }
    for (;b>bStop;*dest--=*b--);
  }
}


template <class T> inline int indexOfFirstElementInBlockGreaterOrEqualTo(T *block, int count, T &v)
{
  int loIndex = 0;
  int hiIndex = count;
  while (loIndex<hiIndex)
  {
    int tryIndex = ( loIndex + hiIndex ) / 2;
    if ( v <= block[tryIndex]  )
      hiIndex  = tryIndex;
    else
      loIndex = tryIndex + 1;
  }
  return hiIndex;
}

template <class T> inline int indexOfFirstElementInBlockGreaterThan(T *block, int count, T &v)
{
  int loIndex = 0;
  int hiIndex = count;
  while (loIndex<hiIndex)
  {
    int tryIndex = ( loIndex + hiIndex ) / 2;
    if ( v < block[tryIndex]  )
      hiIndex  = tryIndex;
    else
      loIndex = tryIndex + 1;
  }
  return hiIndex;
}

template <class T> void paranoidMergeRadix2External(T *a, int aCount, T* b, int bCount, T* dest, bool bLeftToRight, INotifiable *t)
{
  if (aCount==0)
  {
    for (;bCount>0;--bCount) *dest++=*b++;
  }
  else if (bCount==0)
  {
    for (;aCount>0;--aCount) *dest++=*a++;
  }
  else if (bLeftToRight)
  {
    mergeForecastRadix2External(a, a+aCount, b, b+bCount, dest);
  }
  else
  {
    mergeForecastBackwardsRadix2External(a+aCount-1, a-1, b+bCount-1, b-1, dest+aCount+bCount-1);
  }
  if (t!=NULL)
  {
    t->notify();
  }
}

template <class T> void parallelMergeRadix2External(T* a, int aCount, T* b, int bCount, T* dest, bool bLeftToRight, int processorCount, INotifiable *t)
{
  while (processorCount>1)
  {
    //Ick!  We need to split the merge into two parts.
    //We'll be kinda lazy.  We'll figure out where the median of the LHS would go in the RHS, and split that way.
    int middleLeft  = aCount/processorCount*(processorCount/2);
    int hiLeft = middleLeft+1+indexOfFirstElementInBlockGreaterThan(a+middleLeft+1, aCount-middleLeft-1, a[middleLeft]);  
      //index of first element >v, in a, or aCount if there is none
    int loRight   = indexOfFirstElementInBlockGreaterOrEqualTo(b, bCount, a[middleLeft]);
      //index of first element >=v, in b, or bCount if there is none

    //Normally: merge a[0..hiLeft-1] with b[0..loRight-1], so set... takeLeft=hiLeft.
    //But:      We want takeLeft + loRight to be roughly (aCount+bCount)/processorCount*(proessorCount/2), and takeLeft
    //      can be any value between loLeft and hiLeft.  

    int takeLeft     = hiLeft;
    int desiredSplit = (aCount+bCount)/processorCount*(processorCount/2);
    int loLeft       = -1;
    if (takeLeft + loRight > desiredSplit) //too many records will be going to the "left" merge operation
    {
      if (desiredSplit - loRight >= middleLeft)
        takeLeft = desiredSplit - loRight;    //since any value between middleLeft and hiLeft is fine for # records to take from a
      else
      {
        //otherwise: need to know the minimum number of records we could *safely* take from block "a".
        loLeft = indexOfFirstElementInBlockGreaterOrEqualTo(a, middleLeft, a[middleLeft]);
        if (desiredSplit - loRight > loLeft)
        {
          takeLeft = desiredSplit - loRight;
        }
        else
        {
          takeLeft = loLeft;
        }
      }
    }    

    parallelMergeRadix2External(a, takeLeft, b, loRight, dest, bLeftToRight, processorCount/2, t);
    a              += takeLeft;
    b              += loRight;
    aCount         -= takeLeft;
    bCount         -= loRight;
    dest           += takeLeft+loRight;
    processorCount -= processorCount/2;
  }
  gpDispatcher->queueOrDo(CallAction7(paranoidMergeRadix2External<T>, a, aCount, b, bCount, dest, bLeftToRight, t));
}

template <class T> void parallelMergesortSubArrayRadix2(T *src, T *dest, int elementCount
                            , int processorCount, int multiThreadingThreshold
                            , bool srcIsInput, bool isLastMerge, INotifiable *t)
{
  if (elementCount>multiThreadingThreshold)
  {  
    int halfCount = elementCount >> 1;
    IAction*     mergeAction;       //the merge itself
    INotifiable* pendingMerge;    //for notifying the count-down latch that prevents the merge running                

    if (isLastMerge)
    {
      INotifiable* notificationDemultiplexer = NewPendingNotification(processorCount, t);
      mergeAction = CallAction8(parallelMergeRadix2External<T>, src, halfCount, src+halfCount, elementCount-halfCount
              , dest, srcIsInput, processorCount, notificationDemultiplexer);
      pendingMerge = NewPendingAction(2, mergeAction);      

    }
    else 
    {
      if (srcIsInput)
      {
        mergeAction = CallAction5(mergeForecastRadix2External<T>,  src, src+halfCount, src+halfCount, src+elementCount, dest);
      }
      else
      {
        mergeAction = CallAction5(mergeForecastBackwardsRadix2External<T>, src+halfCount-1, src-1, src+elementCount-1, src+halfCount-1, dest+elementCount-1);
      }
      IAction* notifyingAction =  NewNotifyingAction(mergeAction,t);
      pendingMerge = NewPendingAction(2, notifyingAction);
    }        

    parallelMergesortSubArrayRadix2(dest, src, halfCount
      , processorCount, multiThreadingThreshold, !srcIsInput, false, pendingMerge);
    parallelMergesortSubArrayRadix2(dest+halfCount, src+halfCount, elementCount-halfCount
      , processorCount, multiThreadingThreshold, !srcIsInput, isLastMerge, pendingMerge);

  }
  else
  {
    gpDispatcher->queueOrDo( NewNotifyingAction( CallAction5(mergesortForecastRadix2External<T>, src, elementCount, dest, srcIsInput, 32), t));
  }
}

template <class T>   void parallelMergesort(T *base, int elementCount, int processorCount, int multiThreadingThreshold)
{
  T *workArea = new T[elementCount];
  CountDown t(1);
  gpDispatcher->queueOrDo(CallAction8(parallelMergesortSubArrayRadix2<T>,workArea, base, elementCount, processorCount, multiThreadingThreshold, false, 
    true /*set to false to turn off parallization of the last merge on each level*//, (INotifiable*)&t));
  t.wait();
  delete [] workArea;
}

template <class T> void paranoidMergeRadix3External(T *a, int aCount, T* b, int bCount, T* c, int cCount, T* dest, bool bLeftToRight, INotifiable *t)
{
  if (aCount==0) 
    paranoidMergeRadix2External(b, bCount, c, cCount, dest, bLeftToRight, NULL);
  else if (bCount==0)
    paranoidMergeRadix2External(a, aCount, c, cCount, dest, bLeftToRight, NULL);
  else if (cCount==0)
    paranoidMergeRadix2External(a, aCount, b, bCount, dest, bLeftToRight, NULL);
  else if (bLeftToRight)
    mergeRadix3FastForward(a, a+aCount, b, b+bCount, c, c+cCount, dest);
  else
    mergeRadix3FastBackward(a+aCount-1, a-1, b+bCount-1, b-1, c+cCount-1, c-1, dest+aCount+bCount+cCount-1);

  if (t!=NULL)
  {
    t->notify();
  }
}

template <class T> void parallelMergeRadix3External(T* a, int aCount, T* b, int bCount, T* c, int cCount, T* dest, bool bLeftToRight, int processorCount, INotifiable *t)
{
  while (processorCount > 1)
  {
    int aMiddle  = aCount/processorCount*(processorCount/2);
    int bMiddle  = indexOfFirstElementInBlockGreaterThan(b, bCount, a[aMiddle]);
    int cMiddle  = indexOfFirstElementInBlockGreaterThan(c, cCount, a[aMiddle]);
        aMiddle += indexOfFirstElementInBlockGreaterThan(a+aMiddle, aCount-aMiddle, a[aMiddle]);
    parallelMergeRadix3External(a, aMiddle, b, bMiddle, c, cMiddle, dest, bLeftToRight, processorCount/2, t);
    a += aMiddle; aCount -= aMiddle;
    b += bMiddle; bCount -= bMiddle;
    c += cMiddle; cCount -= cMiddle;
    dest += aMiddle + bMiddle + cMiddle;
    processorCount -= processorCount/2;
  }
  gpDispatcher->queueOrDo(CallAction9(paranoidMergeRadix3External<T>, a, aCount, b, bCount, c, cCount, dest, bLeftToRight, t));
}


template <class T> void parallelMergesortSubArrayRadix3(T *src, int elementCount, T *dest, bool isSourceInput,
                          bool isLastMerge, int processorCount, int multiThreadingThreshold,
                          INotifiable* t)
{
  if (elementCount>multiThreadingThreshold)
  {  
    int stepSize = elementCount / 3 ;
    int twoStep  = stepSize + stepSize;
    IAction *mergeAction;
    INotifiable *pendingMerge;

    if (isLastMerge)
    {
      INotifiable *notificationDemultiplexer = NewPendingNotification(processorCount, t);
      mergeAction = CallAction10(parallelMergeRadix3External<T>, src, stepSize, src+stepSize, stepSize, src+twoStep, elementCount-twoStep
              , dest, isSourceInput, processorCount, notificationDemultiplexer);
      pendingMerge = NewPendingAction(3, mergeAction);      
    }
    else 
    {
      if (isSourceInput)
      {
        mergeAction = CallAction7(voidMergeRadix3FastForward<T>,  src, src+stepSize, src+stepSize, src+stepSize+stepSize, src+stepSize+stepSize,
                  src+elementCount, dest);
      }
      else
      {
        mergeAction = CallAction7(voidMergeRadix3FastBackward<T>, src+stepSize-1, src-1, src+stepSize+stepSize-1, src+stepSize-1,
                      src+elementCount-1, src+stepSize+stepSize-1, dest+elementCount-1);
      }
      IAction* notifyingAction =  NewNotifyingAction(mergeAction,t);
      pendingMerge = NewPendingAction(3, notifyingAction);
    }

    parallelMergesortSubArrayRadix3(dest, stepSize, src
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix3(dest+stepSize, stepSize, src+stepSize
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix3(dest+stepSize+stepSize, elementCount-stepSize-stepSize, src+stepSize+stepSize
      , !isSourceInput , isLastMerge, processorCount, multiThreadingThreshold, pendingMerge);
  }
  else
  {
    gpDispatcher->queueOrDo( NewNotifyingAction( CallAction4(mergeSortExternalRadix3Fast<T>, src, elementCount, dest, isSourceInput), t));
  }
}


template <class T>   void parallelMergesortRadix3(T *base, int elementCount, int processorCount, int multiThreadingThreshold=50000)
{
  T *workArea = new T[elementCount];
  CountDown t(1);
  gpDispatcher->queueOrDo(CallAction8(parallelMergesortSubArrayRadix3<T>, workArea, elementCount, base, false, true, processorCount, multiThreadingThreshold, (INotifiable*)&t));
  t.wait();
  delete [] workArea;
}

template <class T> void paranoidMergeRadix4External(T *a, int aCount, T* b, int bCount, T* c, int cCount, T* d, int dCount, T* dest, bool bLeftToRight, INotifiable *t)
{
  if (aCount==0) 
    paranoidMergeRadix3External(b, bCount, c, cCount, d, dCount, dest, bLeftToRight, NULL);
  else if (bCount==0)
    paranoidMergeRadix3External(a, aCount, c, cCount, d, dCount, dest, bLeftToRight, NULL);
  else if (cCount==0)
    paranoidMergeRadix3External(a, aCount, b, bCount, d, dCount, dest, bLeftToRight, NULL);
  else if (dCount==0)
    paranoidMergeRadix3External(a, aCount, b, bCount, c, cCount, dest, bLeftToRight, NULL);
  else if (bLeftToRight)
    mergeRadix4FastForward(a, a+aCount, b, b+bCount, c, c+cCount, d, d+dCount, dest);
  else
    mergeRadix4FastBackward(a+aCount-1, a-1, b+bCount-1, b-1, c+cCount-1, c-1, d+dCount-1, d-1, dest+aCount+bCount+cCount+dCount-1);

  if (t!=NULL)
  {
    t->notify();
  }
}

template <class T> void parallelMergeRadix4External(
  T* a, int aCount, T* b, int bCount, T* c, int cCount, T* d, int dCount,
  T* dest, bool bLeftToRight, int processorCount, INotifiable *t)
{
  while (processorCount > 1)
  {
    int aMiddle  = aCount/processorCount*(processorCount/2);
    int bMiddle  = indexOfFirstElementInBlockGreaterThan(b, bCount, a[aMiddle]);
    int cMiddle  = indexOfFirstElementInBlockGreaterThan(c, cCount, a[aMiddle]);
    int dMiddle  = indexOfFirstElementInBlockGreaterThan(d, dCount, a[aMiddle]);
        aMiddle += indexOfFirstElementInBlockGreaterThan(a+aMiddle, aCount-aMiddle, a[aMiddle]);
    parallelMergeRadix4External(a, aMiddle, b, bMiddle, c, cMiddle, d, dMiddle, dest, bLeftToRight, processorCount/2, t);
    a += aMiddle; aCount -= aMiddle;
    b += bMiddle; bCount -= bMiddle;
    c += cMiddle; cCount -= cMiddle;
    d += dMiddle; dCount -= dMiddle;
    dest += aMiddle + bMiddle + cMiddle + dMiddle;
    processorCount -= processorCount/2;
  }
  gpDispatcher->queueOrDo(CallAction11(paranoidMergeRadix4External<T>, a, aCount, b, bCount, c, cCount, d, dCount, dest, bLeftToRight, t));
}

template <class T> void parallelMergesortSubArrayRadix4(T *src, int elementCount, T *dest, bool isSourceInput,
                          bool isLastMerge, int processorCount, int multiThreadingThreshold,
                          INotifiable* t)
{
  if (elementCount>multiThreadingThreshold)
  {  
    int stepSize  = elementCount >> 2 ;
    int twoStep   = stepSize + stepSize;
    int threeStep = twoStep  + stepSize;
    IAction *mergeAction;
    INotifiable *pendingMerge;

    if (isLastMerge)
    {
      INotifiable *notificationDemultiplexer = NewPendingNotification(processorCount, t);
      mergeAction = CallAction12(parallelMergeRadix4External<T>, src, stepSize, src+stepSize, stepSize
              , src+twoStep, stepSize, src+threeStep, elementCount-threeStep
              , dest, isSourceInput, processorCount, notificationDemultiplexer);
      pendingMerge = NewPendingAction(4, mergeAction);      
    }
    else 
    {
      if (isSourceInput)
      {
        mergeAction = CallAction9(voidMergeRadix4FastForward<T>,  src, src+stepSize, src+stepSize, src+twoStep, src+twoStep,
                  src+threeStep, src+threeStep, src+elementCount, dest);
      }
      else
      {
        mergeAction = CallAction9(voidMergeRadix4FastBackward<T>, src+stepSize-1, src-1, src+twoStep-1, src+stepSize-1,
                  src+threeStep-1, src+twoStep-1, src+elementCount-1, src+threeStep-1, dest+elementCount-1);
      }
      IAction* notifyingAction =  NewNotifyingAction(mergeAction,t);
      pendingMerge = NewPendingAction(4, notifyingAction);
    }

    parallelMergesortSubArrayRadix4(dest, stepSize, src
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix4(dest+stepSize, stepSize, src+stepSize
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix4(dest+twoStep, stepSize, src+twoStep
      , !isSourceInput , false,       processorCount, multiThreadingThreshold, pendingMerge);
    parallelMergesortSubArrayRadix4(dest+threeStep, elementCount-threeStep, src+threeStep
      , !isSourceInput , isLastMerge, processorCount, multiThreadingThreshold, pendingMerge);
  }
  else
  {
    gpDispatcher->queueOrDo( NewNotifyingAction( CallAction4(mergeSortExternalRadix4Fast<T>, src, elementCount, dest, isSourceInput), t));
  }
}

template <class T>   void parallelMergesortRadix4(T *base, int elementCount, int processorCount, int multiThreadingThreshold=50000)
{
  T *workArea = new T[elementCount];
  CountDown t(1);
  gpDispatcher->queueOrDo(CallAction8(parallelMergesortSubArrayRadix4<T>, workArea, elementCount, base, false, true, 
    processorCount, multiThreadingThreshold, (INotifiable*)&t));
  t.wait();
  delete [] workArea;
}