# An LSD radix sort algorithm in IDL

*(Note: I'm at the HDF and HDF-EOS Workshop XV this week, so today I have a guest post by Atle Borsholm, a Senior Consultant in the Exelis VIS Professional Services Group. Atle is an IDL master. I’m trying to persuade him to share some of the cool programs he’s written on this blog. Here’s one. –MP) * The IDL SORT function uses the quicksort (or qsort) algorithm. This is a memory-efficient sorting algorithm that handles any data type that supports comparison (e.g., operators like LT, GT, EQ). However, it can be slow for large input data. Wikipedia says the [proportional] speed of quicksort can range from *n*log(n)* to *n^2* (where *n* is the number of elements). My main motivation for implementing a least significant digit (LSD) radix sort algorithm in IDL was to see if I could get a speed improvement, although a bonus with this algorithm is that it’s a so-called “stable” sort algorithm, meaning that the relative order of equal entries is retained. This is unfortunately not the case with IDL’s SORT. For example:

IDL>sort([2,2,3,2,3,3]), format='(6(i0,x))' 1 3 0 4 5 2

The first three indices in the result [1,3,0] all correspond to the number 2, which in a stable sort would be returned as [0,1,3]. The same goes for the 3s, which should be returned as [2,4,5]. So, a stable sort would give a result of [0,1,3,2,4,5] for this array. An LSD radix sort is designed primarily for arrays of unsigned integer types of any length. The idea is to start sorting on the least significant digit of each array element. Since the digits have a fixed range, a histogram is used to place the entries into bins without regard for the other entries (except retaining the original order within each bin for stability). The algorithm continues by sorting on the next significant digit and so on until there are no more digits. Here’s the first pass:

functionradix_sort_r1, data, radix=radixcompile_optidl2, logical_predicate ; default radix if not specifiedif~keyword_set(radix)thenradix =256radix =long64(radix) sorted = data mx =max(sorted) factor =1ullrepeatbeginmx /= radix rem = sorted / factor digit = remmodradix factor *= radix h =histogram(digit, min=0, max=radix-1, binsize=1, reverse_indices=ri) sorted = sorted[ri[radix+1:*]]endrepuntilmxeq0return, sortedend

Test RADIX_SORT_R1 on the array used above:

IDL>radix_sort_r1([2,2,3,2,3,3]), format='(6(i0,x))' 2 2 2 3 3 3

Sorting an array in-place is nice, but in most cases we want the list of indices — the way IDL’s SORT function returns its result — so we can use the indices on other arrays. So, adding a few lines of code will also support the index as an output:

functionradix_sort_r2, data, radix=radix, index=indexcompile_optidl2, logical_predicate ; default radix if not specifiedif~keyword_set(radix)thenradix =256radix =long64(radix) ; index output is requestedifarg_present(index)thenbeginindex =lindgen(n_elements(data))endifsorted = data mx =max(sorted) factor =1ullrepeatbeginmx /= radix rem = sorted / factor digit = remmodradix factor *= radix h =histogram(digit, min=0, max=radix-1, binsize=1, reverse_indices=ri) sorted = sorted[ri[radix+1:*]]ifarg_present(index)thenindex = index[ri[radix+1:*]]endrepuntilmxeq0return, sortedend

Test on the array used above:

IDL> void =radix_sort_r2([2,2,3,2,3,3], index=index) IDL>

This is the result as expected from a stable sort algorithm. There’s a problem with data type support, though; the above sort only works for non-negative integers. There are some tricks that can be performed on the data to expand support to include negative integers. Basically, we can achieve this by converting to unsigned integer type and shifting by half the integer range:

; pre-processing to support signed intscasetypename(data)of'INT': sorted =uint(data) +ishft(1us,15) 'LONG': sorted =ulong(data) +ishft(1ul,31) 'LONG64': sorted =ulong64(data) +ishft(1ull,63)else: sorted = dataendcase

And, at the end, insert:

casetypename(data)of'INT': sorted =fix(sorted -ishft(1us,15)) 'LONG': sorted =long(sorted -ishft(1ul,31)) 'LONG64': sorted =long64(sorted -ishft(1ull,63))else:endcase

These changes allow negative integers to be sorted. But what about floating point numbers? There’s a more complicated trick that will allow support for floats, at least finite ones. The IEEE standard floating point numbers are ordered in a fairly logical bitwise fashion, so by “casting” the data from double to unsigned 64-bit integers we can manipulate the bits to achieve integers that retain the same order. Here’s the code I added to the CASE statement for pre-processing:

'DOUBLE':beginsorted =ulong64(data,0,n_elements(data)) ; IEEE does not use 2's complement for negative numbers, so flip the ; bits for the negative numbers with this trick. This may not be ; 100% IEEE compliant sorting, but it works well for finite numbers. sortedxor= ((sortedandishft(1ull,63))ne0) * (notishft(1ull,63)) sorted +=ishft(1ull,63)end'FLOAT':beginsorted =ulong(data,0,n_elements(data)) sortedxor= ((sortedandishft(1ul,31))ne0) * (notishft(1ul,31)) sorted +=ishft(1ul,31)end

And for post processing back into the floating point types:

'DOUBLE':beginsorted -=ishft(1ull,63) sortedxor= ((sortedandishft(1ull,63))ne0) * (notishft(1ull,63)) sorted =double(sorted,0,n_elements(data))end'FLOAT':beginsorted -=ishft(1ul,31) sortedxor= ((sortedandishft(1ul,31))ne0) * (notishft(1ul,31)) sorted =float(sorted,0,n_elements(data))end

Now that the LSD radix sort program is complete, we can test its performance:

IDL> data =randomu(1,50000000, /long) ; 50 million long integers IDL> t0 =systime(1) & sort_output =sort(data) &systime(1)-t0 31.062000 IDL> t0 =systime(1) & radix_output =radix_sort(data, index=index) &systime(1)-t0 13.468000 IDL>array_equal(data[sort_output], radix_output),array_equal(sort_output, index) 1 1

RADIX_SORT is more than twice as fast as SORT in this example. The results are identical, even the indices, because RANDOMU doesn’t produce duplicates. If you want to trade memory for speed, you can up the radix from 256 (the default) to 65536. Then RADIX_SORT is even faster:

IDL> t0=systime(1) & radix_output =radix_sort(data, index=index, radix=65536) &systime(1)-t0 8.9219999

Here’s a test with floating-point numbers:

IDL> data =randomn(1,10000000) ; 10 million floats IDL> t0=systime(1) & sort_output =sort(data) &systime(1)-t0 4.8600001 IDL> t0=systime(1) & radix_output =radix_sort(data, index=index, radix=65536) &systime(1)-t0 2.0630000 IDL>array_equal(data[sort_output], radix_output),array_equal(sort_output, index) 1 0

Note that the RADIX_SORT output may have duplicate indices, but the sorted outputs are the same. Here’s a final test using integers with lots of duplicates:

IDL> data =long(randomn(1,30000000)*300) ; 30 million longs with duplicates IDL> t0=systime(1) & sort_output =sort(data) &systime(1)-t0 10.750000 IDL> t0=systime(1) & radix_output =radix_sort(data, index=index, radix=65536) &systime(1)-t0 2.7970002 IDL>array_equal(data[sort_output], data[index]) 1

In this last example we got greater than 3x speed improvement with the added bonus of a stable sort. The final code listing can be downloaded here. *(Thanks, Atle! –MP)* * *