Indirect bucket/radix sort

bs.c

#define N_CHUNKS (16)
#define CHUNK_SIZE (2)
#define CHUNK_MASK (~(~0u << CHUNK_SIZE))
#define SHIFT(C) (CHUNK_SIZE * (C))
#define N_BUCKETS (1 << CHUNK_SIZE)

// Return an integer s and array p that comprise a list of
// indices of array a that put it in ascending sorted order:
// a[s] <= a[next[s]] <= a[next[next[s]]]] <= ...
void ibs(unsigned *a, int n, int *s, int *next) {

#define ENQUEUE(I, K) do { \
  int k = K; \
  if (h[k] == -1) h[k] = next[I] = I; \
  else { \
    int tail = h[k]; \
    int head = next[tail]; \
    h[k] = next[tail] = I; \
    next[I] = head; \
  } \
} while (0)

#define ACCUMULATE_LIST(T) do { \
  T = -1; \
  for (int i = 0; i < N_BUCKETS; ++i) { \
    if (h[i] == -1) continue; \
    if (T == -1) T = h[i]; \
    else { \
      int tail = h[i]; \
      int head = next[tail]; \
      next[tail] = next[T]; \
      next[T] = head; \
      T = tail; \
    } \
    h[i] = -1; \
  } \
} while (0)

  if (n <= 0) {
    *s = -1;
    return;
  }
  int h[N_BUCKETS], list_tail;
  for (int i = 0; i < N_BUCKETS; ++i) h[i] = -1;
  for (int i = 0; i < n; ++i) ENQUEUE(i, a[i] & CHUNK_MASK);
  ACCUMULATE_LIST(list_tail);
  for (int shift = CHUNK_SIZE;
       shift < N_CHUNKS * CHUNK_SIZE;
       shift += CHUNK_SIZE) {
    int src_head = next[list_tail];
    int p = src_head;
    do {
      int p_next = next[p];
      ENQUEUE(p, (a[p] >> shift) & CHUNK_MASK);
      p = p_next;
    } while (p != src_head);
    ACCUMULATE_LIST(list_tail);
  }
  *s = next[list_tail]; next[list_tail] = -1;
}