Double-Pivot Quicksort in C

Question

I've added a sort function to my library of utility functions. I will need to extend its behaviour at a later time but for now it's just a replacement for the qsort() function.

I've implemented Yaroslavskiy's Double Pivot QuickSort which is fast and easy to implement (and has been selected as the sort of choice for Java7).

I've removed recursion via a fixed stack and randomized the choice of the pivot elements to avoid O(n^2) behaviour when the array is already sorted. I've also tested the difference of having a final insertion sort run on the entire array instead of many insertion sort runs on shorter arrays but I've seen no noticeable difference.

As for the stack depth I've tested that the average depth increase logarithmically (as expected). For example for 100 000 elements the average is 22, for 1 000 000 (10x) is 55. The chosen size should be more than enough for any realistic array.

The swap operation is dominant: any improvement will directly impact the overall speed. I tried copying one-byte-at-times and memcpy() but they were slower.

The random number generator utl_rnd() is surely not the best PRNG around but is fast and should be good enough for the purpose.

I'd be especially interested in any comment about speed. Currently this is just a little slower than gcc's qsort() (e.g. 1.6ms vs 1.2ms to sort an array of 10 million random integers).

Code is below (alternatively as github gist: https://gist.github.com/rdentato/8cf8007525feacc8df519bc91e98f20d):

/* Quick and dirty PRNG (xorshift) */
static uint32_t utl_rnd()
{
  static uint32_t rnd = 0;
    while (rnd == 0) rnd = (uint32_t)time(0);
    rnd ^= rnd << 13;
    rnd ^= rnd >> 17;
    rnd ^= rnd << 5;
    return rnd;
}

#define utl_dpqswap(a,b) do { if (a!=b) { \
                                  uint32_t sz = esz;\
                                  uint8_t  tmp8; \
                                  uint32_t tmp32; \
                                  uint8_t *pa = ((uint8_t *)a); \
                                  uint8_t *pb = ((uint8_t *)b); \
                                  while (sz >= 4) { \
                                    tmp32 = *(uint32_t *)pa; \
                                    *(uint32_t *)pa = *(uint32_t *)pb;\
                                    *(uint32_t *)pb = tmp32;\
                                    sz-=4; pa+=4; pb+=4;\
                                  }\
                                  switch (sz) {\
                                    case 3: tmp8=*pa; *pa=*pb; *pb=tmp8; pa--; pb--;\
                                    case 2: tmp8=*pa; *pa=*pb; *pb=tmp8; pa--; pb--;\
                                    case 1: tmp8=*pa; *pa=*pb; *pb=tmp8; pa--; pb--;\
                                  }\
                                }\
                           } while (0)

#define utl_dpqptr(k)    ((uint8_t *)base+(k)*esz)
#define utl_dpqpush(l,r) do {stack[stack_top][0]=(l); stack[stack_top][1]=(r); stack_top++; } while(0)
#define utl_dpqpop(l,r)  do {stack_top--; l=stack[stack_top][0]; r=stack[stack_top][1];} while(0)

/* Dropin replacement for qsort() using double pivot quicksort */
void utl_dpqsort(void *base, uint32_t nel, uint32_t esz, int (*cmp)(const void *, const void *))
{
  int32_t left,right;
  uint8_t *leftptr, *rightptr; 

  uint32_t L,K,G;

  int32_t stack[128][2]; // Enough for 2^31 max elements in the array
  int16_t stack_top = 0;

  utl_dpqpush(0,nel-1);
  while (stack_top > 0) {
    utl_dpqpop(left, right);
    if (left < right) {
      if ((right - left) <= 16) {  // Use insertion sort
        for (int32_t i = left+1; i<=right; i++) {
          rightptr = utl_dpqptr(i);
          leftptr = rightptr - esz;
          for (int32_t j=i; j>0 && cmp(leftptr, rightptr) > 0; j--) {
            utl_dpqswap(rightptr, leftptr);
            rightptr = leftptr;
            leftptr = rightptr - esz;
          }
        }
      }
      else {
        leftptr = utl_dpqptr(left);
        rightptr = utl_dpqptr(right);

        /* Randomize pivot to avoid worst case (already sorted array) */
        L = left + (utl_rnd() % (right-left));
        G = left + (utl_rnd() % (right-left));
        utl_dpqswap(utl_dpqptr(L),leftptr);
        utl_dpqswap(utl_dpqptr(G),rightptr);

        if (cmp(leftptr, rightptr) > 0) {
          utl_dpqswap(leftptr, rightptr);
        }
        L=left+1; K=L; G=right-1;
        while (K <= G) {
          if (cmp(utl_dpqptr(K), leftptr) < 0) {
            utl_dpqswap(utl_dpqptr(K), utl_dpqptr(L));
            L++;
          }
          else if (cmp(utl_dpqptr(K), rightptr) > 0) {
            while ((cmp(utl_dpqptr(G), rightptr) > 0) && (K<G)) 
              G--;

            utl_dpqswap(utl_dpqptr(K), utl_dpqptr(G));
            G--;
            if (cmp(utl_dpqptr(K), leftptr) < 0) {
              utl_dpqswap(utl_dpqptr(K), utl_dpqptr(L));
              L++;
            }
          }
          K++;
        }
        L--; G++;
        utl_dpqswap(leftptr, utl_dpqptr(L));
        utl_dpqswap(rightptr, utl_dpqptr(G));

        utl_dpqpush(G+1, right);
        utl_dpqpush(L+1, G-1);
        utl_dpqpush(left, L-1);
      }
    }
  }
}

Answer 1

Tips on general sorting performance

I recently implemented a sorting library (in c++) for the use in my novel suffix array construction algorithm (SACA) and did a lot of research for it. You asked on how to improve performance and as algorithm usally beats implementation I want to give you some tips on the quicksort algorithm in general:

1. Branch misses:

The paper "How Branch Mispredictions Affect Quicksort" - Kaligosi, Sanders showed that when sorting simple types (e.g. integers) branch misspredictions are a huge slow down. I've implemented and use the ideas from the really ground breaking paper "BlockQuicksort: How Branch Mispredictions don't affect Quicksort" - Edelkamp, Weiss. It's almost twice as fast as std::sort() on integer data.

Note that regular quicksort favors slightly skewed distributions to achieve better branch prediction. This is not the case for BlockQuicksort: It really favors good pivots (I use up to 65 elements to find the pivots, from an information theoric background median-of-O(sqrt(n)) would be optimal for BlockQuicksort).

2. Cache misses:

Cache misses are a huge problem when sorting pointers/indices. The amount of cache misses can be greatly reduced by increasing the number of pivots. I currently use 3 to balance between swaps and cache trashing. Note that this is still faster than vanilla std::sort() on integer data. The extrem form of this is called superscalar sample sort.

The paper "How Good is Multi-Pivot Quicksort?" - Aumueller, Dietzfelbinger, Klaue is very interessting on this topic. This also shows that 3 or 5 pivots result in a minimum number of scanned elements so if your access is slow you want to use 3 or 5 pivots.

Number of scanned elements in the whole sort on average:
  3 pivots: 1.385 * n * log2(n) 
  5 pivots: 1.379 * n * log2(n)

Number of copied elements in the whole sort on average:
  3 pivots: 1.108 * n * log2(n) 
  5 pivots: 1.248 * n * log2(n)

In my implementation I use a heuristic to choose between copying together key&value + BlockQuicksort and three-pivot quicksort when sorting pointers to integers but this is specially tuned towards my SACA. Generally copying + BlockQuicksort is favored for bigger lists.

3. Complexity

Sorting takes \$O(n * log(n))\$ time in general. But there are corner cases like almost sorted lists or lists with a lot of duplicates. My implementation of BlockQuicksort switches to three-way quicksort (<, = and > partitions with single pivot) when the sample I use to find the pivots has a lot of duplicates. This way performance of \$O(n * log(m))\$ can be achieved where m is the number of different values.

pdqsort nicely tackles not so random permutations. I greatly encourage you to look at the source code on github for info and benchmark data (it also uses BlockQuicksort). It is able to detect partially sorted subsequences and sometimes sort them in \$O(n)\$ time.

4. Indirections

The standard qsort() interface has a major drawback when compared to c++ std::sort(): It uses function pointers. Calling them almost always results in a branch missprediction and is a huge slow down. This bypasses the idea of BlockQuicksort to get rid of branches.

Sadly it's quite complex to get a general sorting routine in c without function pointers because of the lack of templates. So if you really need a high performance generic sorting routine in c you probably need to plug in some preprocessor tricks (one of the reasons I switched to c++).

5. Insertion sort

Insertion sort is also the finisher in my implementation. I use it as soon as the list is small enough (somewhere between 16-64 elements). This is simply because when sorting pointers this results in the values still being in cache from the last partition call.

Using it afterwards results in less branches but requires to load everything another time which isn't a big problem when sorting integers but rather heavy for pointers.

6. Sorting networks

I use sorting networks to find the pivots which is really fast but I have yet to find a way that outperforms insertion sort when sorting pointers. This is due to the compare-and-swap operation being way slower when the values to be compared differ from the values being swapped (conditional move + xor swap doesn't work).

7. Heap sort

Another performance improvement really was to switch to heap sort with deep recursion. I don't know why I ran into them but adding in heapsort made my sort around 3-4% faster in my tests.

8. Stack usage

I found using recursive calls to be faster than using a separated stack. So you should probably test both. Stackoverflow is not a problem if you limit recursion by switching to heap sort after a limited number of recursions. Sadly I don't understand exactly why but from looking at the assembly it seems the compiler had problems optimizing the separated stack.

---

Rather special info:

My library also differs from the standard std::sort() interface: the comparision function. When sorting pointers a comparision function has to load both args and compare them. That is not necessary with quicksort: The pivot can sometimes be stored in a register and only one element needs to be loaded. That's why I use an index() function which gets the value currently to be partitioned (e.g. the pointer) and returns a value which is used to actually be compared (e.g. a tuple to sort structs or the value pointed to). This also allows me to call the index function only once per partitioned element and compare it multiple times against the pivots. This way I treat L1 reads for register access.

// Example for sorting a struct data e.g. a points:
// with comparision:
std::sort(elements.begin(), elements.end(), [](auto a, auto b) {
   return std::tie(a.x, a.y) < std::tie(b.x, b.y);
   // this would be a simple bug:
   // return std::tie(a.x, a.y) < std::tie(b.y, b.x);
});

// with index:
sort_index(elements.begin(), elements.end(), [](auto a) {
   return std::tie(a.x, a.y);
});

This example shows that usage is less error prone. It doesn't match watch people learn about sorting and therefore I wouldn't recomment it in a general library. If a comparision function is really needed then one can implement an index function that returns a special type which overloads the comparision function.

---

Finally note that some of those tips are strictly tuned towards my application and I've did a lot of benchmarking and optimization to find the best solutions for my problem of suffix array construction. I hope you can use a lot of it for you because this post contains around 2-3 month of research and development. Feel free to ask for more info.

Answer 2

Suggestion for limiting stack depth

Currently, if your partitions happen to be the worst case partitions (think sorted array and you randomly pick the leftmost indices every time), you will overflow your stack.

You can limit your stack depth to a logarithmic limit even with worst case partitioning if you push your subpartitions correctly. If you always push your largest partitions before your smaller ones, you will handle the smaller partitions first and this will prevent your stack from growing past a logarithmic limit. (This is because the largest partition will be popped off the stack before being handled so when you handle it, you will not have grown the stack at all.)

For regular quicksort which partitions into two partitions, this technique limits the stack depth to \$\log_2 n\$. For your quicksort with three partitions, the limit should be \$2*\log_3 n\$.