3.0

# An LSD radix sort algorithm in IDL

Mark Alonzo
(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> print, 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:
```function radix_sort_r1, data, radix=radix
compile_opt idl2, logical_predicate

; default radix if not specified

sorted = data
mx = max(sorted)
factor = 1ull
repeat begin
rem = sorted / factor
h = histogram(digit, min=0, max=radix-1, binsize=1, reverse_indices=ri)
endrep until mx eq 0

return, sorted
end```
Test RADIX_SORT_R1 on the array used above:
```IDL> print, 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:
```function radix_sort_r2, data, radix=radix, index=index
compile_opt idl2, logical_predicate

; default radix if not specified

; index output is requested
if arg_present(index) then begin
index = lindgen(n_elements(data))
endif

sorted = data
mx = max(sorted)
factor = 1ull
repeat begin
rem = sorted / factor
h = histogram(digit, min=0, max=radix-1, binsize=1, reverse_indices=ri)
if arg_present(index) then index = index[ri[radix+1:*]]
endrep until mx eq 0

return, sorted
end```
Test on the array used above:
```IDL> void = radix_sort_r2([2,2,3,2,3,3], index=index)
IDL> print, index, format='(6(i0,x))'
0 1 3 2 4 5```
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 ints
case typename(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 = data
endcase```
And, at the end, insert:
```   case typename(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': begin
sorted = 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.
sorted xor= ((sorted and ishft(1ull,63)) ne 0) * (not ishft(1ull,63))
sorted += ishft(1ull, 63)
end
'FLOAT': begin
sorted = ulong(data,0,n_elements(data))
sorted xor= ((sorted and ishft(1ul,31)) ne 0) * (not ishft(1ul,31))
sorted += ishft(1ul, 31)
end```
And for post processing back into the floating point types:
```      'DOUBLE': begin
sorted -= ishft(1ull, 63)
sorted xor= ((sorted and ishft(1ull,63)) ne 0) * (not ishft(1ull,63))
sorted = double(sorted, 0, n_elements(data))
end
'FLOAT': begin
sorted -= ishft(1ul, 31)
sorted xor= ((sorted and ishft(1ul,31)) ne 0) * (not ishft(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) & print, systime(1)-t0
31.062000
IDL> t0 = systime(1) & radix_output = radix_sort(data, index=index) & print, systime(1)-t0
13.468000
IDL> print, 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) & print, 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) & print, systime(1)-t0
4.8600001
2.0630000
IDL> print, 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) & print, systime(1)-t0
10.750000
2.7970002
IDL> print, 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)
MOST RECENT

### U.S. DoD Partner Uses DTM Data to Mitigate Effects of Climate Change

12/2/2020

The Earth’s climate is rapidly changing and few areas on the planet have been left untouched.... more »

### Geospatial Answers as a Service | EAS 2020 Panel Discussion

11/18/2020

Join Stuart Blundell (Presagis) as he moderates this panel discussion about Geospatial Answers as a... more »