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

lsd_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];
  }
}

// Use LSD radix sort (given lg_base) to sort in_elems with max value and return in out_elems.
void lsdRadixSort(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];

  // Bucket offsets.
  uint64_t offsets[n_counts];

  // Perform counting sort for least-significant digit.
  countingSort(lg_base, src_elems, n_elems, dst_elems, 0, counts, offsets);

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

  // Repeat for each remaining digit.
  for (int i_digit = 1; i_digit < n_digits; i_digit++) {
    // 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);
  }
}