Background⌗

The computational complexity of an optimal sorting algorithm based on comparisons is O(nlog⁡n)\mathcal O(n \log n). We can prove this by reducing the problem of sorting an array to finding the permutation vector that sorts it. There are n!n! possible permutations of the set {1,2,…,n}\{1, 2, \ldots, n\}. Each comparison can at best halve the search space of the sorting algorithm. Hence, the amount of comparisons done by an optimal comparison sort is at least log⁡2(n!)=O(nlog⁡n)\log_2 (n!) = \mathcal O(n \log n).

Indeed, one such algorithm - the heap sort - attains the O(nlog⁡n)\mathcal O(n \log n) bound and is used by many sorting routines under the bonnet. To avoid unnecessarily discussing the underlying methods which are already explained great detail by algorithm books and other online resources, I will not provide an abstract explanation of the algorithm, focusing on the implementation:

void heapsort(float * x, size_t n) {
  size_t i = n / 2, parent, child;
  for (float t;; x[parent] = t) {
    if (i > 0) t = x[--i];
    else if (--n == 0) return;
    else t = x[n], x[n] = x[0];
    parent = i;  child = i * 2 + 1;
    while (child < n) {
      if (child + 1 < n && x[child + 1] > x[child]) child++;
      if (x[child] > t) {
        x[parent] = x[child];
        parent = child;
        child = parent * 2 + 1;
      } else break;
    }
  }
}

The heap sort algorithm is pretty slow in practice. Another common algorithm - quicksort - doesn’t provide as lucrative guarantees, but generally runs faster than the alternatives:

void quicksort(float * x, size_t n) {
  if (n < 2) return;
  float p = x[n / 2], * l = x, * r = x + n - 1, t;
  while (l <= r) {
    if (*l < p) l++;
    else if (*r > p) r--;
    else t = *l, *l++ = *r, *r-- = t;
  }
  quicksort(x, r - x + 1);
  quicksort(l, x + n - l);
}

We can run a simple benchmark to demonstrate this:

int main(void) {
  srand(time(0));
  float * x = malloc(0x1000000 * sizeof(float));
  for (size_t i = 0; i < 0x1000000; i++) x[i] = (float) rand();
  #define TEST_SORT(f, tab, name) \
    volatile float * tab = malloc(0x1000000 * sizeof(float)); \
    { \
      memcpy(tab, x, 0x1000000 * sizeof(float)); \
      uint64_t start = __rdtsc(); \
      struct timespec t0, t1; \
      clock_gettime(CLOCK_MONOTONIC, &t0); \
      f(tab, 0x1000000); \
      uint64_t end = __rdtsc(); \
      clock_gettime(CLOCK_MONOTONIC, &t1); \
      double elapsed = t1.tv_sec - t0.tv_sec + 1e-9 * (t1.tv_nsec - t0.tv_nsec); \
      printf(name " sort: %ld cycles, %f seconds\n", end - start, elapsed); \
    }
  TEST_SORT(heapsort, a1, "heap");
  TEST_SORT(quicksort, a2, "quick");
}

On my machine, I get the following timings:

heap sort: 15121146414 cycles, 4.590787 seconds (901 cycles/element)
quick sort: 4824085530 cycles, 1.464595 seconds (287 cycles/element)

More efficient algorithms like TimSort would be able to take better advantage of what linearithmic comparison-based sorts have to offer, but we will cut this story short for now - its performance will still be significantly worse than the best we can do.

IEEE754 floats⌗

Assuming non-degenerate data (i.e. no NaNs, infinities, etc.) we can sort IEEE754 floats via naive memory comparisons. This motivates the following concise implementation of the radix sort algorithm, which attains linear time:

void radix_sort(float * x, size_t n) {
  uint32_t b1[0x100] = { 0 }, b2[0x100] = { 0 }, si[0x100] = { 0 };
  uint32_t * result = calloc(n, 4), * y = (uint32_t *) x, key, temp, i;
  for (size_t i = 0; i < n; i++) y[i] ^= 0x80000000;
  #define DO_BUCKET(expr) \
    i = 0; if (n % 2) ++b1[expr]; while (i < n) ++b1[expr], ++b2[expr]; \
    for (i = 0; i < 0x100; ++i) b1[i] += b2[i]; memset(b2, 0, 0x400);
  #define SCAN_BUCKET \
    si[0] = 0; for (int j = 1; j < 0x100; ++j) si[j] = si[j - 1] + b1[j - 1];
  #define PERMUTE1(expr) \
    for (int i = n-1; i >= 0; --i) key = expr, result[si[key] + --b1[key]] = y[i];
  #define PERMUTE2(expr) \
    for (int i = n-1; i >= 0; --i) key = expr, y[si[key] + --b1[key]] = result[i];
  DO_BUCKET(y[i++] & 0xFF); SCAN_BUCKET; PERMUTE1(y[i] & 0xFF);
  DO_BUCKET((result[i++] & 0xFF00) >> 8); SCAN_BUCKET; PERMUTE2((result[i] & 0xFF00) >> 8);
  DO_BUCKET((y[i++] & 0xFF0000) >> 16); SCAN_BUCKET; PERMUTE1((y[i] & 0xFF0000) >> 16);
  DO_BUCKET((result[i++] & 0xFF000000) >> 24); SCAN_BUCKET; PERMUTE2((result[i] & 0xFF000000) >> 24);
  free(result);
  for (i = 0; i < n; i++) y[i] ^= 0x80000000;
  for (key = 0, i = 0; i < n; i++) if (y[i] & 0x80000000) key++;
  for (i = 0; i < key / 2; i++) temp = y[i], y[i] = y[key - i - 1], y[key - i - 1] = temp;
}

An alternative implementation of the same idea involves being a bit more wasteful with stack space to avoid superfluous memory accesses and allow batching of operations:

void radix_sort2(float * x, size_t n) {
  uint32_t b[0x800] = { 0 }, i, * result = calloc(n, 4), * y = (uint32_t *) x, * b1, * b2, key;
  for (size_t i = 0; i < n; i++) y[i] ^= 0x80000000;
  #define DO_BUCKET(z, expr) { \
    b1 = b + z; b2 = b + 0x100 + z; \
    i = 0; if (n % 2) ++b1[expr]; while (i < n) ++b1[expr], ++b2[expr]; \
    for (i = 0; i < 0x100; ++i) b1[i] += b2[i]; }
  #define SCAN_BUCKET(z) { \
    b1 = b + z; b2 = b + 0x100 + z; \
    b2[0] = 0; for (int j = 1; j < 0x100; ++j) b2[j] = b2[j - 1] + b1[j - 1]; \
    for (int j = 0; j < 0x100; ++j) b1[j] += b2[j]; }
  DO_BUCKET(0x000, y[i++] & 0xFF); DO_BUCKET(0x200, (y[i++] & 0xFF00) >> 8);
  DO_BUCKET(0x400, (y[i++] & 0xFF0000) >> 16); DO_BUCKET(0x600, (y[i++] & 0xFF000000) >> 24);
  SCAN_BUCKET(0x000); SCAN_BUCKET(0x200); SCAN_BUCKET(0x400); SCAN_BUCKET(0x600);
  #define PERMUTE1(z, expr) { \
    b1 = b + z; for (int i = n - 1; i >= 0; --i) key = expr, result[--b1[key]] = y[i]; }
  #define PERMUTE2(z, expr) { \
    b1 = b + z; for (int i = n - 1; i >= 0; --i) key = expr, y[--b1[key]] = result[i]; }
  PERMUTE1(0x000, y[i] & 0xFF); PERMUTE2(0x200, (result[i] & 0xFF00) >> 8);
  PERMUTE1(0x400, (y[i] & 0xFF0000) >> 16); PERMUTE2(0x600, (result[i] & 0xFF000000) >> 24);
  for (i = 0; i < n; i++) y[i] ^= 0x80000000;
  for (key = 0, i = 0; i < n; i++) if (y[i] & 0x80000000) key++;
  for (i = 0; i < key / 2; i++) temp = y[i], y[i] = y[key - i - 1], y[key - i - 1] = temp;
}

The new timings (including the radix sorts) follow:

radix sort: 352613448 cycles, 0.107044 seconds
radix2 sort: 508113969 cycles, 0.154266 seconds
heap sort: 15164556891 cycles, 4.603966 seconds
quick sort: 4501389717 cycles, 1.366624 seconds

The radix sort implementation is approximately 12 times faster.