CS Inside · Sorting

The sort that can't have the best case

  • Sorting

Introduction

You probably know this sorting algorithm. You probably learned it back in school, or during your first few computer science classes at university. You probably can implement it with your eyes closed. Still, I think there are a few interesting things about it you may never have thought about.

This sorting algorithm is selection sort.

When I was actively learning algorithms and data structures back in school and university, I loved exploring algorithms inside and out. Usually, it led me much further than just the implementation. Selection sort is no exception. I think we can find much more around these few lines of code than it may seem at first.

Yes, at first glance, there is not that much to say about it. The algorithm itself is relatively trivial. Its implementation is relatively easy and small, Its complexity is well known and easy to analyze. All of this is true, at least to a certain degree. Still, I think it is a very good algorithm to experiment with precisely because of its simplicity. It gives us a small and understandable playground where we can observe things that are usually hidden behind more complicated logic. In this post, we will take a look at one of these things - branch prediction.

So, selection sort...

Selection sort

Let's start with a relatively boring part. As I mentioned at the beginning of this post, I am sure that you know this algorithm, even though it is probably not something you have used in production (although you still can find 1, 2). So, just in case you forgot the details, let's take a look at how Donald Knuth defines it in his classic TAOCP:

Another important family of sorting techniques is based on the idea of repeated selection. The simplest selection method is perhaps the following:

  1. Find the smallest key; transfer the corresponding record to the output area; then replace the key by the value $$\infty$$ (which is assumed to be higher than any actual key).

  2. Repeat step 1. This time the second smallest key will be selected, since the smallest key has been replaced by $$\infty$$.

  3. Continue repeating step 1 until $$N$$ records have been selected.

Donald E. Knuth, The Art of Computer Programming, vol. 3, section 5.2.3, "Sorting by Selection"

Remember it? I suppose so.

Knuth starts with the simplest form of selection: repeatedly find the smallest element and move it to the output. In practice, however, we usually do not need a separate output area. At least I have never seen it in practice. Obviously, we can keep everything in the same array and, on every iteration, just place the smallest remaining element into its final position.

This gives us the selection sort implementation that most of us remember:

void selection_sort(void *base, size_t n, size_t size,
                    int (*comparator)(const void *, const void *)) {
    unsigned char *p = base;

    for (size_t i = 0; i + 1 < n; i++) {
        size_t current_min_idx = i;

        for (size_t j = i + 1; j < n; j++) {
            if (comparator(p + j * size, p + current_min_idx * size) < 0) {
                current_min_idx = j;
            }
        }

        if (current_min_idx != i) {
            swap_bytes(p + i * size, p + current_min_idx * size, size);
        }
    }
}

void swap_bytes(void *a, void *b, size_t size) {
    unsigned char *p = a;
    unsigned char *q = b;

    while (size--) {
        unsigned char tmp = *p;
        *p++ = *q;
        *q++ = tmp;
    }
}

Remembered? Good.

Here is a visualization of the same algorithm running on seven numbers:

5 2 9 1 7 3 8 0 1 2 3 4 5 6

Seven values in no particular order.

Now that the algorithm is fresh in our heads again, let’s start looking at it a little more carefully.

The complexity of this sorting algorithm is well known - $$O(n^2)$$, but it is also easy to analyze it to be sure. Using simple words, the outer for loop does $$n - 1$$ iterations. On each of these iterations, there is an inner for loop that does $$n - i - 1$$ iterations. i grows from $$0$$ to $$n - 2$$, so we have the following number of operations:

$$ (n - 1) + (n - 2) + (n - 3) + \ldots + (n - (n - 2) - 1) $$

or

$$ (n - 1) + (n - 2) + (n - 3) + \ldots + 1 $$

Since these terms form an arithmetic progression, we can easily calculate the sum, which is $$\frac{n(n - 1)}{2}$$. So, in terms of big-O it is really $$O(n^2)$$.

But let's return for a moment to the step right before the analysis. What did we try to count? I named it operations. What are these operations? In fact, we have counted the number of comparisons. Yes, there is no best, average, or worst case in terms of comparisons. Selection sort always does exactly $$\frac{n(n - 1)}{2}$$ comparisons, absolutely independent of the input order. Even an already sorted array does not help it. On the other hand, if we count the number of swaps, we will see that it is between $$0$$ and $$n - 1$$. This is already something that $$O(n^2)$$ does not tell us.

A comparison and a swap are two very different operations. There is no reason to assume that they have the same cost.

Actually, we do not even need to go far to see it. Look again at our implementation. A comparison is a call to comparator, while a swap is a call to swap_bytes, which moves the specified number of bytes one by one. You can say that it is not the most efficient implementation and we can move more than 1 byte at once. But that does not really matter here. If we sort strings, the comparator may need to inspect several characters before it can decide which string is smaller. If we sort large records, swapping two elements may require copying quite a lot of data anyway.

So, yes, selection sort cannot have a best case in comparisons because of its structure.

A traditional tutorial about selection sort would probably stop here. I think it is only the beginning. Yes, maybe you will not see this algorithm often in production, but it is so simple to understand and implement that it is nice to experiment with and learn some important lessons about code optimization. Let's see, for example, whether we can make it faster.

Can we make it faster?

OK, so we know that the performance of this sorting algorithm should not depend on the input. But can we come up with a faster implementation than the one above? For sure, there are traditional ways such as SIMD or parallel algorithms. But since we are experimenting, I think there could be something more interesting.

You can often hear that fewer branches make code faster. That can be true when branches are hard to predict, because a misprediction is expensive. But removing a branch is not automatically an optimization. It is always interesting to experiment to see how things behave in reality.

Let's see if we can get rid of this comparison:

if (comparator(p + j * size, p + current_min_idx * size) < 0) {
    current_min_idx = j;
}

The answer is both yes and no. We can not fully get rid of the comparison, because in any case we need to determine whether the current element is smaller than the minimum we have seen so far. But we can ask the compiler to look at it from a different angle. What we may be able to remove is the control-flow branch associated with updating the minimum.

NOTE

All my machines are x86-64, so here I will rely on the capabilities of this architecture.

There are at least two common ways the compiler can implement this condition:

  1. Conditional jump
  2. Conditional move

The difference between these approaches is that in the first case the compiler may generate something like this:

call comparator
test eax, eax
jns <skip next instruction>
mov current_min_idx, j

And in the second case, something like this:

call comparator
test eax, eax
cmovs current_min_idx, j

Which one is better? Need to measure.

But before trying to understand what is better and what is not, let's see what the compiler generates for our code. I tried both GNU GCC and Clang using only the -O2 -march=x86-64 flags in both cases.

For this condition:

if (comparator(p + j * size, p + current_min_idx * size) < 0) {
    current_min_idx = j;
}

The compilers produce something like this:

GNU GCC conditional move
  1. movrsi, r15
  2. imulrsi, r14
  3. addrsi, rbp
  4. movrdi, r12
  5. call[QWORD PTR [rsp]]
  6. testeax, eax
  7. cmovsr14, r13

Instructions run sequentially, no branch to predict.

Clang conditional jump
  1. movrsi, rbx
  2. imulrsi, r14
  3. addrsi, r12
  4. movrdi, r13
  5. callqword ptr [rsp + 80]
  6. movrcx, rbp
  7. testeax, eax
  8. js.LBB3_5

Instructions may jump, so the CPU predicts where.

So, GNU GCC already assumes that the cmov instruction will be better for this case. To compare both variants, we should force the compiler to use the cmov instruction in the first case and the jump instruction in the second.

There is very little point in rewriting the condition using the ternary operator. I could try, but writing the condition differently in C does not guarantee a cmov anyway. To avoid guessing compiler flags and relying on compiler heuristics in general, we can try to implement these small pieces explicitly using inline assembly. Here we run up against yet another superstition - "compiler knows better". Maybe it is true. But I think it is still worth trying. At least out of curiosity.

Let's take a look at my first attempt to rewrite the selection sort implementation to force the compiler to use the cmov instruction:

void selection_sort_cmov(void *base, size_t n, size_t size,
                         int (*comparator)(const void *, const void *)) {
    unsigned char *p = base;

    for (size_t i = 0; i + 1 < n; i++) {
        unsigned char *current_min = p + i * size;

        for (size_t j = i + 1; j < n; j++) {
            unsigned char *candidate = p + j * size;
            int res = comparator(candidate, current_min);

            __asm__(
                    "testl %[res], %[res]\n\t"
                    "cmovsq %[candidate], %[min]"
                    : [min] "+r" (current_min)
                    : [res] "r" (res),
                      [candidate] "r" (candidate)
                    : "cc"
                   );
        }

        if (current_min != p + i * size) {
            swap_bytes(p + i * size, current_min, size);
        }
    }
}

Now the compiler must emit cmov instead of the conditional jump. Did it get faster? Let's see.

What the measurements say

I ran the test with 32 thousand elements arranged in the array in a special way. Half of its elements are smaller than everything that stands before them, and these elements are put at random positions. This means that the swap is executed in about half of the cases, and there is no pattern in it that could be learned.

Here is what perf showed for the test program built with GNU GCC and with Clang:

Counter GCC jump GCC cmov Clang jump Clang cmov
cycles 7,566,060,605 6,038,949,794 3,004,740,343 5,760,465,017
instructions 9,874,308,573 9,739,620,541 3,151,747,340 6,159,000,872
branches 2,053,607,286 1,538,441,094 831,555,533 515,891,909
branch-misses 123,276,288 233,201 93,648,212 215,223
L1-dcache-loads 2,052,376,291 2,051,942,243 1,029,945,404 1,028,959,022
L1-dcache-load-misses 30,704,676 30,651,016 30,084,737 30,662,880
LLC-load-misses 10,560 9,879 15,702 30,145
time elapsed 2.189633685 s 1.741969967 s 0.874063903 s 1.680908270 s

Quite interesting, isn't it?

Under GNU GCC, the version with cmov looks better. But at the same time, the version with cmov is almost 2x slower if we build the program with Clang. We can see that using cmov, we clearly get fewer branch misses. The Intel® 64 and IA-32 Architectures Optimization Reference Manual says:

The largest challenge with mispredicted branches is finding the branch which caused them. Branch mispredictions incur penalty of about 20 cycles.

Twenty cycles per miss is not nothing, and probably it is enough to explain the difference between the two GCC builds. Even if this is a rough number, and the exact value depends on many different factors including the CPU microarchitecture, it says nothing about Clang, where the cmov build misses more than 400 times fewer branches and still loses. Something in the cmov version costs more than the mispredictions it saves. So what is the problem?

Knowing the exact time complexity of the selection sort algorithm, we can calculate the number of iterations our loops do. For 32 thousand elements it is:

iterations = n * (n - 1) / 2 = 32000 * 31999 / 2 = 511,984,000

Having this information, we can calculate how much one iteration costs. Let's do it and compare everything:

build time cycles/iter instructions/iter IPC branch-misses/iter L1 loads/iter
gcc jump 2.190 s 14.78 19.29 1.31 0.2408 4.01
gcc cmov 1.742 s 11.80 19.02 1.61 0.0005 4.01
clang jump 0.874 s 5.87 6.16 1.05 0.1829 2.01
clang cmov 1.681 s 11.25 12.03 1.07 0.0004 2.01

We can see that both cmov builds cost nearly the same per iteration - 11.80 and 11.25 cycles. This is very strange, because GCC's cmov loop runs 19.02 instructions per iteration and Clang's runs 12.03, so GCC executes about 58% more instructions. A loop that executes fewer instructions and finishes in almost the same time is telling us that it waits for something.

This leaves us with two questions: why do the jump builds lose so many cycles to wrong guesses, and what does the cmov build wait for?

Let's start with the first one.

Branch prediction in action

Let's take a look at what perf record will say for the Clang run without cmov:

16.21 :   400ea0:        addq    $0x2, %rdi
 2.24 :   400ea4:        movq    %rsi, %r8
 3.26 :   400ea7:        cmpq    %rbx, %rdi
 0.00 :   400eaa:        je      0x400e7c <main+0xdc>
12.18 :   400eac:        movl    (%r14,%rdi,4), %r9d
 0.91 :   400eb0:        movq    %rdi, %rsi
 4.27 :   400eb3:        cmpl    (%r14,%r8,4), %r9d
 2.04 :   400eb7:        jl     0x400ebc <main+0x11c>
13.07 :   400eb9:        movq    %r8, %rsi
23.59 :   400ebc:        movl    0x4(%r14,%rdi,4), %r8d
 3.20 :   400ec1:        cmpl    (%r14,%rsi,4), %r8d
 5.58 :   400ec5:        jge    0x400ea0 <main+0x100>
13.38 :   400ec7:        leaq   0x1(%rdi), %rsi
 0.07 :   400ecb:        jmp    0x400ea0 <main+0x100>

First of all, we can note that there is no call to int_cmp at all. Everything was inlined and the inner loop was unrolled, it processes two elements per iteration. Clang seems understood that the function pointer always points to the same function, put its body directly into the loop, and then noticed that the result of (x > y) - (x < y) is only used to answer "is it less than zero?". So the whole comparator collapsed into a single cmpl and a single conditional jump.

Besides that, we can note one more interesting thing. Almost all the most expensive or hot places are located immediately after the jumps. This can be the branch misprediction in action! We can re-run perf record with the -e branch-misses:pp flag, collecting samples where branch misprediction occurred, and that gives exactly these two comparisons:

 0.00 :   400ea0:        addq    $0x2, %rdi
 0.00 :   400ea4:        movq    %rsi, %r8
 0.00 :   400ea7:        cmpq    %rbx, %rdi
 0.05 :   400eaa:        je      0x400e7c <main+0xdc>
 0.00 :   400eac:        movl    (%r14,%rdi,4), %r9d
 0.00 :   400eb0:        movq    %rdi, %rsi
29.59 :   400eb3:        cmpl    (%r14,%r8,4), %r9d
22.26 :   400eb7:        jl      0x400ebc <main+0x11c>
 0.00 :   400eb9:        movq    %r8, %rsi
 0.00 :   400ebc:        movl    0x4(%r14,%rdi,4), %r8d
24.99 :   400ec1:        cmpl    (%r14,%rsi,4), %r8d
23.02 :   400ec5:        jge     0x400ea0 <main+0x100>
 0.00 :   400ec7:        leaq    0x1(%rdi), %rsi
 0.00 :   400ecb:        jmp     0x400ea0 <main+0x100>

These are exactly the places where the CPU did not guess the proper branch and "wasted" cycles. As predicted, the version without cmov spends time on recovering from the wrong guesses. Every wrong guess means that the CPU throws away the work it already started and does the same iteration again.

Since the CPU sometimes did not guess the proper continuation, at other times it did guess right, and those places run fast. If this theory is true, we can try to run the test with the array that is initially sorted in descending order and get better numbers. I ran such a test, and let's see what I got:

Counter 50% unsorted descending
cycles 2,924,064,866 779,178,169
instructions 3,128,083,743 3,205,408,576
branches 830,319,519 906,714,060
branch-misses 93,103,953 195,197
L1-dcache-loads 1,024,584,745 1,019,640,213
L1-dcache-load-misses 29,665,807 29,539,116
LLC-load-misses 1,010 701
time elapsed 0.843521867 s 0.225241245 s

The number of reads from memory did not change and the caches behaved exactly the same. The number of instructions grew a little, by 2.5%, because in the array sorted in descending order every element is a new minimum, so the update is executed every time. And the only counter that really significantly changed is branch-misses, from 93 million down to 195 thousand. The change of this single counter led to a 3.7x performance improvement!

So yes, branch prediction definitely does its job. But what about cmov? We reduced the branch misses, but at the same time we eliminated the branch prediction. Let's look at what perf record says about the version with the cmov instruction:

10.33 :   400f20:        movl    (%r8), %r9d
42.44 :   400f23:        cmpl    (%rsi), %r9d
 9.00 :   400f26:        setl    %r9b
 1.87 :   400f2a:        setg    %r10b
12.81 :   400f2e:        subb    %r9b, %r10b
 7.43 :   400f31:        movsbl  %r10b, %r9d
 7.20 :   400f35:        testl   %r9d, %r9d
 7.65 :   400f38:        cmovsq  %r8, %rsi
 1.20 :   400f3c:        addq    $0x4, %r8
 0.00 :   400f40:        decq    %rdi
 0.00 :   400f43:        cmpq    %rdi, %rcx
 0.02 :   400f46:        jne     0x400f20 <main+0x180>

Here we can see the completely different picture:

At every step of the inner loop there is always the same question - is the given value smaller than the smallest one it has seen so far? This cmpl instruction is exactly the answer to this question. If we annotate the assembly code, we will see that the instructions from 400f23 to 400f31 compute (x > y) - (x < y). To calculate this value, the CPU must know values of x and y. In our case, x is the current candidate and y is the current minimum we have seen so far.

The r8 and rsi registers contain the addresses of these two values. The address of x is in the r8 register and the address of y is in the rsi register. The important thing here is that the address in rsi is written by the cmovsq instruction at 400f38 at each previous iteration of the inner loop.

So, the cmov instruction by itself is not slow. The problem is that it makes every iteration of the inner loop wait for the previous one.

Each iteration of the inner loop does the same three things:

To read the next element from the array, the CPU does not need to know anything except the address of the array and offset to the element. In this case, nothing here depends on the previous iteration. Comparing two values does not depend on it either. But the address of the minimum value is written by the cmov instruction at each previous iteration of the inner loop. So the iterations go strictly one after another, and each of them waits for the one before it.

A modern CPU does not execute instructions one by one. It works more like a conveyor belt. While one instruction is still in progress, the next ones can already be fetched, decoded and sent to an execution unit. But this works only while the CPU knows what instructions should be executed next or at least it can try to guess. Because of the dependencies described above, the CPU does not have this information.

The version with the conditional jump does not wait at all. The CPU just guesses the result of the comparison and continues. The address of the new minimum can be only one of two values, and both of them are already in the registers, so the next comparison starts immediately and multiple of them can be in flight at the same time. Yes, sometimes the guess is wrong and all this work must be thrown away, but as we have seen above, the CPU guesses right much more often than it is wrong.

Now it should be clear that my "optimization" did not optimize anything, but actually made things worse. I still think it was interesting to experiment with.

Yes, the idea failed, and we did not end up with a faster version of the algorithm. But I can not say that the time was wasted. Often, in engineering, negative results can be just as valuable as successful ones. They tell us which assumptions were wrong and often teach us more than an optimization that happened to work on the first try. In this case, we saw branch prediction in action and not just as something described in a textbook, and that, I think, is very valuable by itself.

So although the optimization failed, I think the experiment did not.

Before you go

You may notice that we did not pay much attention to the code produced by GNU GCC. Interestingly, even the version with jump is definitely slower than the same version compiled with Clang. I found some factors why it is like this, but before we take a look at them, let's take another look at the original numbers:

Counter GCC jump Clang jump
time elapsed 2.189633685 s 0.874063903 s

I have compared the assembly output of both compilers and found three differences:

  1. GCC does not inline functions and makes a real call
  2. GCC does not unroll the loops
  3. GCC produces more instructions in the inner loop

The first two problems we can solve by adding inline to functions and passing the -funroll-loops compiler flag. After this, the execution time for the same number of elements becomes close to Clang:

 ./selection_sort branch 32000

real    0m0.971s
user    0m0.968s
sys	    0m0.000s

I still have not found a way to force GCC to reduce the number of instructions in the loop. If you experiment with it and figure out how to do it, please let me know on X.