# Finding the two largest of five small integers as quickly as possible

I use a variation of a 5-cross median filter on image data on a small embedded system, i.e.

    x
x x x
x


The algorithm is really simple: read 5 unsigned integer values, get the highest 2, do some calculations on those and write back the unsigned integer result.

What is nice is that the 5 integer input values are all in the range of 0-20. The calculated integer value are also in the 0-20 range!

Through profiling, I have figured out that getting the largest two numbers is the bottleneck so I want to speed this part up. What is the fastest way to perform this selection?

The current algorithm uses a 32 bit mask with 1 in the position given by the 5 numbers and a HW-supported CLZ function.
I should say that the CPU is a proprietary one, not available outside of my company. My compiler is GCC but tailor made for this CPU.

I have tried to figure out if I can use a lookup-table but I have failed to generate a key that I can use.

I have $21^5$ combinations for the input but order isn't important, i.e. [5,0,0,0,5] is the same as [5,5,0,0,0].

It happens that the hash-function below produces a perfect hash without collisions!

def hash(x):
h = 0
for i in x:
h = 33*h+i
return h


But the hash is huge and there is simply not enough memory to use that.

Is there a better algorithm that I can use? Is it possible to solve my problem using a lookup-table and generating a key?

• Which algorithm do you currently use? Seven integer comparisons are enough, is that too slow? Your hash already performs more operations. Are subsequent calls to the method related, e.g. does the central x move through through the matrix row-by-row? – Raphael Feb 24 '15 at 22:54
• The filter is convolved through the image row by row. I.e get the 5 values and do the calculations then move everything one step to the right and repeat. The hash was only an example. I benchmarked several sliding-window solutions to minimize reading of data but it all boils down to finding the highest 2 values. – Fredrik Pihl Feb 24 '15 at 23:03
• Most likely your algorithm, if implemented properly, would be bounded by memory access and not by computation. Using a hashtable would only increase the amount of memory accesses and slow things down. Please post your current code so we can see how it can be improved - I believe only micro-optimization is possible. The most I can think of is: maybe we can take advantage of the fact that 2 values are in common between neighbouring windows? – jkff Feb 25 '15 at 5:10
• @jkff Depending on matrix, cache sizes and (cache) mapping function, every value might only have have to be loaded once; most operations should run on registers or L1 cache then. Pipelining is another issue, though. – Raphael Feb 25 '15 at 7:33
• By the way, do you do this in parallel already? This seems particularly suitable for vector parallelisation or SIMD (e.g. on a GPU). That route would help much more than save a few percent per cell. – Raphael Feb 25 '15 at 7:38

In my other answer I suggest that conditional jumps might be the main impediment to efficiency. As a consequence, sorting networks come to mind: they are data agnostic, that is the same sequence of comparisons is executed no matter the input, with only the swaps being conditional.

Of course, sorting may be too much work; we only need the biggest two numbers. Lucky for us, selection networks have also been studied. Knuth tells us that finding the two smallest numbers out of five² can be done with $\hat{U}_2(5) = 6$ comparisons [1, 5.3.4 ex 19] (and at most as many swaps).

The network he gives in the solutions (rewritten to zero-based arrays) is

$\qquad\displaystyle [0:4]\,[1:4]\,[0:3]\,[1:3]\,[0:2]\,[1:2]$

which implements -- after adjusting the direction of the comparisons -- in pseudocode as

def selMax2(a : int[])
a.swap(0,4) if a < a
a.swap(1,4) if a < a
a.swap(0,3) if a < a
a.swap(1,3) if a < a
a.swap(0,2) if a < a
a.swap(1,2) if a < a
return (a, a)
end


Now, naive implementations still have conditional jumps (across the swap code). Depending on your machine you can cirumvent them with conditional instructions, though. x86 seems to be its usual mudpit self; ARM looks more promising since apparently most operations are conditional in themselves. If I understand the instructions correctly, the first swap translates to this, assuming our array values have been loaded to registers R0 through R4:

CMP     R0,R4
MOVLT   R5 = R0
MOVLT   R0 = R4
MOVLT   R4 = R6


Yes, yes, of course you can use XOR swapping with EOR.

I just hope your processor has this, or something similar. Of course, if you build the thing for this purpose, maybe you can get the network hard-wired on there?

This is probably (provably?) the best you can do in the classical realm, i.e. without making use of the limited domain and performing wicked intra-word magicks.

1. Sorting and Searching by Donald E. Knuth; The Art of Computer Programming Vol. 3 (2nd ed, 1998)
2. Note that this leaves the two selected elements unordered. Ordering them requires an extra comparison, that is $\hat{W}_2(5) = 7$ many in total [1, p234 Table 1].
• I'm accepting this. I received a lot of new ideas that I need to benchmark before moving on. Referring to Knuth always works for me :-) Thanks for your effort and time! – Fredrik Pihl Feb 26 '15 at 15:11
• @FredrikPihl Cool, please let us know how it turns out in the end! – Raphael Feb 26 '15 at 15:21
• I will! Reading Chapter 5.3.3 right now. Love the start of the it with references to Lewis Carroll and the tennis tournament :-) – Fredrik Pihl Feb 26 '15 at 15:23
• Depending on the instruction set, using 2*max(a,b) = a + b + abs(a-b) along with the selection network could be useful; it could be less costly than unpredictable conditional jumps (even without an intrinsic or conditional move for abs: gcc, at least for x86, generate a jumpless sequence which doesn't seem to be dependend on x86). Having a jumpless sequence is also useful when combined with SIMD or a GPU. – AProgrammer Feb 27 '15 at 14:14

Just so that it's on the table, here's a direct algorithm:

// Sort x1, x2
if x1 < x2
M1 = x2
m1 = x1
else
M1 = x1
m1 = x2
end

// Sort x3, x4
if x3 < x4
M2 = x4
m2 = x3
else
M2 = x3
m2 = x4
end

// Pick largest two
if M1 > M2
M3 = M1
if m1 > M2
m3 = m1
else
m3 = M2
end
else
M3 = M2
if m2 > M1
m3 = m2
else
m3 = M1
end
end

// Insert x4
if x4 > M3
m3 = M3
M3 = x4
else if x4 > m3
m3 = x4
end


By clever implementation of if ... else, one can get rid of some unconditional jumps a direct translation would have.

This is ugly but takes only

• five or six comparisons (i.e. conditional jumps),
• nine to ten assignments (with 11 variables, all in registers) and

In fact, six comparisons is optimal for this problem as Theorem S in section 5.3.3 of  shows; here we need $W_2(5)$.

This can not be expected to be fast on machines with pipelining, though; given they high percentage of conditional jumps, most time would probably be spent in stall.

Note that a simpler variant -- sort x1 and x2, then insert the other values subsequently -- takes four to seven comparisons and only five to six assignments. Since I expect jumps to be of higher cost here, I stuck with this one.

1. Sorting and Searching by Donald E. Knuth; The Art of Computer Programming Vol. 3 (2nd ed, 1998)
• I wonder what an optimizing compiler can do with these. – Raphael Feb 25 '15 at 8:07
• I'll implement this and benchmark it against the current CLZ-based solution. Thanks for your time! – Fredrik Pihl Feb 25 '15 at 8:09
• @FredrikPihl What was the result of your benchmarks? – Raphael Mar 17 '15 at 21:03
• SWAP based approach beats CLZ! On mobile now. Can post more data another time, on mobile now – Fredrik Pihl Mar 17 '15 at 21:34
• @FredrikPihl Cool! I'm happy the good old theory approach can (still) be of practical use. :) – Raphael Mar 17 '15 at 21:36

This could be a great application and test case for the Souper project. Souper is a superoptimizer -- a tool that takes a short sequence of code as input, and tries to optimize it as much as possible (tries to find an equivalent sequence of code that will be faster).

Souper is open source. You might try running Souper on your code snippet to see if it can do any better.

See also John Regehr's contest on writing fast code to sort 16 4-bit values; it's possible that some of the techniques there might be useful.

• I'd be interested in what this can do on the programs the OP's been trying. – Raphael Feb 27 '15 at 7:55

You can use a $21^3$ table that gets three integers and outputs the largest two. You can then use three table lookups:

T[T[T[441*a+21*b+c]*21+d]*21+e]

Similarly, using a $21^4$ table, you can reduce it to two table lookups, though it's not clear that this would be faster.

If you really want a small table, you can use two $21^2$ tables to "sort" two numbers, and then use a sorting network. According to Wikipedia, this requires at most 18 table lookups (9 comparators); you might be able to do with less since (1) you only want to know the two largest elements, and (2) for some comparator gates, you might only be interested in the maximum.

You can also use a single $21^2$ table. Implementing a sorting network then uses less memory accesses but more arithmetic. This way you get at most 9 table lookups.