Power-of-two base MSD radix sort for small arrays of unsigned 64-bit integers.

msd_radix.c

#include <stddef.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>

#define CEIL_DIV(n, d) (((n) / (d)) + (((n) % (d)) != 0))

#define N_DIGITS(lg_base, elem) CEIL_DIV(64 - __builtin_clzll(elem), lg_base)
#define DIGIT_AT(lg_base, elem, i_digit) (((elem) >> ((i_digit) * (lg_base))) & ((1lu << (lg_base)) - 1))

// Swap two uint64_t array pointers.
void swap(uint64_t **a, uint64_t **b) {
  uint64_t *c;
  c = *a;
  *a = *b;
  *b = c;
}

// Perform counting sort using i_digit (given lg_base) of src_elems and return in out_elems, with counts and bucket offsets.
void countingSort(int lg_base, uint64_t *src_elems, size_t n_elems, uint64_t *dst_elems, int i_digit, size_t *counts, size_t *offsets) {
  // Reset all counts.
  size_t n_counts = 1lu << (1 << lg_base);
  memset(counts, 0, n_counts * sizeof(size_t));

  // Set counts for current digit.
  for (size_t i_elem = 0; i_elem < n_elems; i_elem++) {
    counts[DIGIT_AT(lg_base, src_elems[i_elem], i_digit)]++;
  }

  // Set bucket offsets for current digit.
  size_t i_elem = 0;
  for (size_t i_count = 0; i_count < n_counts; i_count++) {
    offsets[i_count] = (i_elem += counts[i_count]);
  }

  // Copy elements from source array to destination array at corresponding offset.
  for (size_t i_elem = n_elems; i_elem > 0; i_elem--) {
    dst_elems[--offsets[DIGIT_AT(lg_base, src_elems[i_elem - 1], i_digit)]] = src_elems[i_elem - 1];
  }
}

// Recursively perform MSD radix sort starting with the i_digit (given lg_base) in src_elems and return in dst_elems.
void msdRadixSortRecursive(int lg_base, uint64_t *src_elems, size_t n_elems, uint64_t *dst_elems, int i_digit) {
  // Counts for each digit.
  size_t n_counts = 1lu << (1 << lg_base);
  size_t counts[n_counts];

  // Offsets for buckets.
  size_t offsets[n_counts];

  // Perform counting sort for current digit.
  countingSort(lg_base, src_elems, n_elems, dst_elems, i_digit, counts, offsets);

  // Swap source and destination arrays.
  swap(&src_elems, &dst_elems);

  // Only recurse if digits remain.
  if (i_digit > 0) {
    // Recursively apply radix sort for each remaining digit, if nonempty bucket.
    for (size_t i_count = 0; i_count < n_counts; i_count++) {
      if (counts[i_count]) {
        msdRadixSortRecursive(lg_base, src_elems + offsets[i_count], counts[i_count], dst_elems + offsets[i_count], i_digit - 1);
      }
    }
  }
}

// Use MSD radix sort (given lg_base) to sort in_elems with max value and return in out_elems.
void msdRadixSort(int lg_base, uint64_t *in_elems, size_t n_elems, uint64_t *out_elems, uint64_t max) {
  // Temporary element array.
  uint64_t tmp_elems[n_elems];

  // Number of rounds for radix sort (one for each digit).
  int n_digits = N_DIGITS(lg_base, max);

  // Initialize source and destination element arrays.
  uint64_t *src_elems = in_elems;
  uint64_t *dst_elems = (n_digits % 2) ? out_elems : tmp_elems;

  // Counts for each digit.
  size_t n_counts = 1lu << (1 << lg_base);
  size_t counts[n_counts];

  // Offsets for buckets.
  size_t offsets[n_counts];

  // Perform counting sort on most-significant digit.
  countingSort(lg_base, src_elems, n_elems, dst_elems, n_digits - 1, counts, offsets);

  // Swap source and destination arrays.
  src_elems = dst_elems;
  dst_elems = (n_digits % 2) ? tmp_elems : out_elems;

  // Only recurse if digits remain.
  if (n_digits > 1) {
    // Recursively apply radix sort for each remaining digit, if nonempty bucket.
    for (size_t i_count = 0; i_count < n_counts; i_count++) {
      if (counts[i_count]) {
        msdRadixSortRecursive(lg_base, src_elems + offsets[i_count], counts[i_count], dst_elems + offsets[i_count], n_digits - 2);
      }
    }
  }
}