Поразрядная сортировка LSD (Radix Sort)

В последнее время публикуется множество статей посвященных различным алгоритмам сортировки и их сравнению, решил внести и свои пять копеек.

Хочу рассказать про свой излюбленный алгоритм для поразрядной сортировки LSD (least significant digit — сначала младший разряд) с подсчётом (Radix Sort). Классический алгоритм был несколько переосмыслен автором в сторону некоторых оптимизаций в пользу ускорения и простоты.

Итак, предложенная сортировка является устойчивой. Сортировать будем целые 32 битные числа. Для работы потребуется ~(n+4Кбайт) дополнительной памяти, что несколько расточительно, но позволяет добиться некоторого увеличения производительности.

В данной разновидности LSD не используются сравнения и обмены, алгоритм полностью линеен. Вычислительная сложность O(N).

Основной особенностью алгоритма является высокая эффективность для сильно перемешанных или случайных наборов данных. На почти отсортированных наборах имеет смысл применять другие алгоритмы, так как выигрыш будет не так значителен. Плохо работает на малых массивах, менее пары сотен элементов.

Сортировка происходит по месту, для экономии памяти.

//==================================================
// RADIX сортировка (авторская версия by rebuilder)
// Устойчивый алгоритм сортировки, полностью линеен.
// Вычислительная сложность О(n), доп память ~(n+4k)
//==================================================
procedure RSort(var m: array of Longword);
//-------------------------------------------------- 
   procedure Sort_step(var source, dest, offset: array of Longword; const num: Byte);
   var i,temp : Longword;
       k : Byte;
   begin
     for i := High(source) downto 0 do
     begin
       temp := source[i];
       k := temp SHR num;
       dec(offset[k]);
       dest[offset[k]] := temp;
     end;
   end;
//-------------------------------------------------- 
// Объявляем массив корзин первым, для выравнивания на стеке
var s : array[0..3] of array[0..255] of Longword;  
    i,k : longword;
    // Смещение байт внутри переменной k
    offset : array[0..3] of byte absolute k;       
    m_temp : array of Longword;
begin
  SetLength(m_temp, Length(m));                
  // Быстрая очистка корзин    
  FillChar(s[0], 256 * 4 * SizeOf(Longword), 0);   

  // Заполнение корзин
  for i := 0 to High(m) do                         
  begin
    k := m[i];
    Inc(s[0,offset[0]]);
    Inc(s[1,offset[1]]);
    Inc(s[2,offset[2]]);
    Inc(s[3,offset[3]]);
  end;

  // Пересчёт смещений для корзин
  for i := 1 to 255 do                             
  begin
    Inc(s[0,i], s[0,i-1]);
    Inc(s[1,i], s[1,i-1]);
    Inc(s[2,i], s[2,i-1]);
    Inc(s[3,i], s[3,i-1]);
  end;

  // Вызов сортировки по байтам от младших к старшим
  Sort_step(m, m_temp, s[0], 0);                   
  Sort_step(m_temp, m, s[1], 8);
  Sort_step(m, m_temp, s[2], 16);
  Sort_step(m_temp, m, s[3], 24);

  SetLength(m_temp, 0);                            
end;
//==================================================
...
SetLength(m, n);
for i := 0 to n - 1 do 
  m[i] := Random(65536 * 65536);  
...
RSort(m);  
...

Код написан на языке Pascal, но не составит труда портировать его на любой удобный Вам язык.

Последовательность выполнения состоит из двух этапов:

1. Для каждого блока (по восемь двоичных разрядов – 1 байт (оптимальная величина)), методом подсчёта, вычисляется его положение в новом массиве.
2. Последовательно для каждого блока (от младших разрядов к старшим) производится перемещение на рассчитанную ранее позицию.

Улучшения:

1. Для массива смещений используем выравнивание в памяти, а благодаря небольшому объёму он помещается в L1 – кэша процессора.
2. Массив смещений заполняется сразу для всех разрядов, что позволяет пройтись по массиву для подсчёта только однажды.
3. Расчёт позиций начинается не с головы массива, а с конца, это решает две задачи:
   
   • конец массива на первом проходе уже находится в «прогретом» кэше, что при больших массивах, даёт небольшое ускорение;
   • во-вторых нисходящий цикл к нулю короче на одну инструкцию ассемблера, на каждом шаге цикла, относительно восходящего цикла.
4. Для каждой итерации (из четырёх) не используется вложенный цикл, пусть так менее красиво, но экономятся ещё несколько инструкций процессора.

Из-за простоты код практически идентичен по скорости как 32, так и 64 битной сборке компилятора. При необходимости легко представить версию алгоритма для 16 и 64 разрядных чисел.

Сравнение алгоритма на случайной выборке с быстрой сортировкой на 64 битной платформе (в среднем, по десять проходов).



Предложения и замечания по поводу оптимизаций одобряются.

Спасибо.