6
\$\begingroup\$

I am trying to optimize the following (now C) radix sort code for use in my game engine library. The basis for this code was inspired by this youtube video on a C++ implementation called skasort that used templating to sort STL containers. The current first pass code sorts in-place in most-significant-byte (MSB) order and is designed to handle [u]int[8/16/32/64]_t types as well as float/f32, double/f64 and char * strings of variable lengths. Here is the raw C code (header-only) I wrote that my cython library is wrapping:

#include <stdint.h>
#include <stdbool.h>

#define RADIX 256
#define RADIX_THRESHOLD 128

#define CMP_CHECK(a, b, cmp_type) \
cmp_type a_i = *(cmp_type *)a; \
cmp_type b_i = *(cmp_type *)b; \
if(a_i < b_i) \
{ \
    return -1; \
} \
else if(a_i > b_i) \
{ \
    return 1; \
} \
else \
{ \
    return 0; \
} \

#define CMP_FUNC(a, b, cmp_type, cmp_suffix) \
int cmp_func_##cmp_suffix(const void *a, const void *b) \
{ \
    CMP_CHECK(a, b, cmp_type) \
} \

typedef enum RadixSortType
{
    RADIX_SORT_TYPE_U8,
    RADIX_SORT_TYPE_I8,
    RADIX_SORT_TYPE_U16,
    RADIX_SORT_TYPE_I16,
    RADIX_SORT_TYPE_U32,
    RADIX_SORT_TYPE_I32,
    RADIX_SORT_TYPE_U64,
    RADIX_SORT_TYPE_I64,
    RADIX_SORT_TYPE_F32,
    RADIX_SORT_TYPE_F64,
    RADIX_SORT_TYPE_STR,
} RadixSortType;

typedef int (* CmpFuncC)(const void *, const void *);
typedef uint8_t (* RadixKeyFuncC)(void *item, size_t byte_offset);

typedef struct RadixPartitionTableC
{
    size_t counts[RADIX];
    size_t prefix_sums[RADIX + 1];
    size_t shifted_sums[RADIX + 1];
} RadixPartitionTableC;

void table_clear(RadixPartitionTableC *table)
{
    memset(table->counts, 0, sizeof(size_t) * RADIX);
    memset(table->prefix_sums, 0, sizeof(size_t) * (RADIX + 1));
    memset(table->shifted_sums, 0, sizeof(size_t) * (RADIX + 1));
}

void item_swap(void *items, size_t a, size_t b, size_t size)
{
    uint8_t *a_ptr = (uint8_t *)items + (a * size);
    uint8_t *b_ptr = (uint8_t *)items + (b * size);
    uint8_t tmp;
    for(size_t i = 0; i < size; i++)
    {
        tmp = a_ptr[i];
        a_ptr[i] = b_ptr[i];
        b_ptr[i] = tmp;
    }
}

CMP_FUNC(a, b, uint8_t, u8)
CMP_FUNC(a, b, int8_t, i8)
CMP_FUNC(a, b, uint16_t, u16)
CMP_FUNC(a, b, int16_t, i16)
CMP_FUNC(a, b, uint32_t, u32)
CMP_FUNC(a, b, int32_t, i32)
CMP_FUNC(a, b, uint64_t, u64)
CMP_FUNC(a, b, int64_t, i64)
CMP_FUNC(a, b, float, f32)
CMP_FUNC(a, b, double, f64)
int cmp_func_str(const void *a, const void *b)
{
    return strcmp(*(char **)a, *(char **)b);
}

uint8_t radix_key_func_u8(void *item, size_t byte_offset)
{
    return ((uint8_t *)item + byte_offset)[0];
}

uint8_t radix_key_func_i8(void *item, size_t byte_offset)
{
    return ((uint8_t *)item + byte_offset)[0] + 128;
}

uint8_t radix_key_func_u16(void *item, size_t byte_offset)
{
    return ((uint8_t *)item + byte_offset)[0];
}

uint8_t radix_key_func_i16(void *item, size_t byte_offset)
{
    uint16_t value = ((uint16_t *)item)[0] + 32768;
    return ((uint8_t *)&value + byte_offset)[0];
}

uint8_t radix_key_func_u32(void *item, size_t byte_offset)
{
    return ((uint8_t *)item + byte_offset)[0];
}

uint8_t radix_key_func_i32(void *item, size_t byte_offset)
{
    uint32_t value = ((int32_t *)item)[0] + (int32_t)2147483648;
    return ((uint8_t *)&value + byte_offset)[0];
}

uint8_t radix_key_func_u64(void *item, size_t byte_offset)
{
    return ((uint8_t *)item + byte_offset)[0];
}

uint8_t radix_key_func_i64(void *item, size_t byte_offset)
{
    uint64_t value = ((int64_t *)item)[0] + (int64_t)9223372036854775808;
    return ((uint8_t *)&value + byte_offset)[0];
}

uint8_t radix_key_func_f32(void *item, size_t byte_offset)
{
    uint32_t value = ((uint32_t *)item)[0];
    uint32_t mask = -(int32_t)(value >> (uint32_t)31) | (uint32_t)2147483648;
    uint32_t shifted_value = value ^ mask;
    return ((uint8_t *)(&shifted_value) + byte_offset)[0];
}

uint8_t radix_key_func_f64(void *item, size_t byte_offset)
{
    uint64_t value = ((uint64_t *)item)[0];
    uint64_t mask = -(int64_t)(value >> (uint64_t)63) | (uint64_t)9223372036854775808;
    uint64_t shifted_value = value ^ mask;
    return ((uint8_t *)(&shifted_value) + byte_offset)[0];
}

uint8_t radix_key_func_str(void *item, size_t byte_offset)
{
    return (((uint8_t **)item)[0] + byte_offset)[0];
}

void c_radix_sort(void *items, size_t item_size, 
        size_t start, size_t end, size_t start_offset, 
        RadixSortType type_, RadixKeyFuncC key_func)
{
    size_t byte_offset;
    size_t num_bytes;
    bool flip_byte_order = true;
    bool use_null_term = false;
    CmpFuncC cmp_func;
    
    switch(type_)
    {
        case RADIX_SORT_TYPE_U8:
            byte_offset = 0;
            num_bytes = sizeof(uint8_t);
            key_func = radix_key_func_u8;
            cmp_func = cmp_func_u8;
            break;
        case RADIX_SORT_TYPE_I8:
            byte_offset = 0;
            num_bytes = sizeof(int8_t);
            key_func = radix_key_func_i8;
            cmp_func = cmp_func_i8;
            break;
        case RADIX_SORT_TYPE_U16:
            byte_offset = 1;
            num_bytes = sizeof(uint16_t);
            key_func = radix_key_func_u16;
            cmp_func = cmp_func_u16;
            break;
        case RADIX_SORT_TYPE_I16:
            byte_offset = 1;
            num_bytes = sizeof(int16_t);
            key_func = radix_key_func_i16;
            cmp_func = cmp_func_i16;
            break;
        case RADIX_SORT_TYPE_U32:
            byte_offset = 3;
            num_bytes = sizeof(uint32_t);
            key_func = radix_key_func_u32;
            cmp_func = cmp_func_u32;
            break;
        case RADIX_SORT_TYPE_I32:
            byte_offset = 3;
            num_bytes = sizeof(int32_t);
            key_func = radix_key_func_i32;
            cmp_func = cmp_func_i32;
            break;
        case RADIX_SORT_TYPE_U64:
            byte_offset = 7;
            num_bytes = sizeof(uint64_t);
            key_func = radix_key_func_u64;
            cmp_func = cmp_func_u64;
            break;
        case RADIX_SORT_TYPE_I64:
            byte_offset = 7;
            num_bytes = sizeof(int64_t);
            key_func = radix_key_func_i64;
            cmp_func = cmp_func_i64;
            break;
        case RADIX_SORT_TYPE_F32:
            byte_offset = 3;
            num_bytes = sizeof(float);
            key_func = radix_key_func_f32;
            cmp_func = cmp_func_f32;
            break;
        case RADIX_SORT_TYPE_F64:
            byte_offset = 7;
            num_bytes = sizeof(double);
            key_func = radix_key_func_f64;
            cmp_func = cmp_func_f64;
            break;
        case RADIX_SORT_TYPE_STR:
            byte_offset = 0;
            num_bytes = sizeof(char *);
            flip_byte_order = false;
            use_null_term = true;
            key_func = radix_key_func_str;
            cmp_func = cmp_func_str;
            break;
    }
    c_radix_sort_byte(
        items, item_size, start, end, start_offset,
        byte_offset, num_bytes, flip_byte_order, 
        use_null_term, key_func, cmp_func
    );
}

inline void c_radix_sort_byte(void *items, size_t item_size, 
        size_t start, size_t end, size_t start_offset, 
        size_t byte_offset, size_t num_bytes, bool flip_byte_order, 
        bool use_null_term, RadixKeyFuncC key_func, CmpFuncC cmp_func)
{
    void *item_ptr;
    uint8_t item;
    size_t total;
    size_t prefix_sum;
    size_t shifted_sum;
    size_t a, b;
    RadixPartitionTableC table;
    size_t count;
    size_t table_start;
    size_t table_end;
    int8_t byte_order = flip_byte_order ? -1 : 1;
    
    table_clear(&table);
    for(size_t i = start; i < end; i++)
    {
        item_ptr = ((uint8_t *)items) + (i * item_size);
        item = key_func(item_ptr, byte_offset);
        table.counts[item] += 1;
    }
    
    if(use_null_term)
    {
        table.counts[0] = 0;
    }
    
    total = 0;
    for(size_t i = 0; i < RADIX; i++)
    {
        total += table.counts[i];
        table.prefix_sums[i + 1] += total;
        table.shifted_sums[i] = table.prefix_sums[i + 1];
    }
    
    for(size_t i = start; i < end; i++)
    {
        while(true)
        {
            item_ptr = ((uint8_t *)items) + (i * item_size);
            item = key_func(item_ptr, byte_offset);
            prefix_sum = table.prefix_sums[item];
            shifted_sum = table.shifted_sums[item];
            if(prefix_sum == shifted_sum)
            {
                break;
            }
            a = i;
            b = start + prefix_sum;
            item_swap(items, a, b, item_size);
            table.prefix_sums[item] += 1;
        }
    }

    if(!flip_byte_order && byte_offset == num_bytes && !use_null_term)
    {
        return;
    }
    else if(flip_byte_order && byte_offset == 0 && !use_null_term)
    {
        return;
    }
    else
    {
        total = 0;
        for(size_t i = 0; i < RADIX; i++)
        {
            count = table.counts[i];
            table_start = start + total;
            table_end = start + total + count;
            total += count;
            if(count >= RADIX_THRESHOLD)
            {
                c_radix_sort_byte(
                    items, item_size, table_start, table_end, 
                    start_offset, byte_offset + byte_order, num_bytes, 
                    flip_byte_order, use_null_term, key_func, cmp_func
                );
            }
            else if(count > 1)
            {
                qsort(
                    items + (item_size * table_start) + start_offset, 
                    table_end - table_start, item_size, 
                    cmp_func
                );
            }
        }
    }
}

The benchmarking test I have is written in cython and runs radix sorts on each of the supported types (with the char * test being on random alphanumeric strings ranging from 20-50 chars in length), on a custom vector-like container from n=2^1 to n=2^24 items. I can add the benchmarking code if needed; in the meantime, the code for the game engine repo can be found here. The performance of this code is slower than I would like, ranging from 5x faster in the simple uint8_t case to about 0.9x the speed of the radix sort in the worst cases where the sort code needs to be recursively called (the 64-bit types and the variable length strings). I have the following questions:

  1. The code swaps items in-place, leading to a bunch of memcpy calls when items need to be swapped. Would pointer swapping and sorting out of place provide a meaningful performance benefit? This would come at the cost of allocating additional heap memory, which I was striving to avoid based on that above lecture.
  2. I have heard about "loop unrolling" being a possible technique to improve performance. How would I go about implementing this? Would SIMD be needed to do this? Any resources on this would be helpful.
  3. I would like to generalize this to a series of passes based on struct members to do more complex sorts based on these primitive types (e.g. sorting players by distance from a target, alphabetizing user names, etc.). I would appreciate any ideas for a convenient, user-friendly API that could handle arbitrary structs.
  4. The code is essentially already C-code. Would rewriting in C and making a cython wrapper provide any meaningful performance improvement to help the compiler optimize further? The code is already being compiled with reasonable compiler arguments in gcc ( -std=c11, -O3, -ffast-math, -march=native). This point is moot now that the code has been rewritten in C for both convenience (macros to deduplicate code) and reviewability (cython is a niche language). The code is still being compiled with the same flags and the performance is essentially unchanged.

As requested by @Reinderien, this is the complete gcc command that is invoked when the radix sort is compiled by gcc:

building 'pyorama.algs.radix_sort' extension
gcc -Wno-unused-result -Wsign-compare -DNDEBUG -march=x86-64 -mtune=generic -O2 -pipe -fwrapv -D__USE_MINGW_ANSI_STDIO=1 -D_WIN32_WINNT=0x0601 -DNDEBUG -march=x86-64 -mtune=generic -O2 -pipe -fwrapv -D__USE_MINGW_ANSI_STDIO=1 -D_WIN32_WINNT=0x0601 -DNDEBUG -DCIMGUI_DEFINE_ENUMS_AND_STRUCTS=True -DCGLTF_IMPLEMENTATION=True -DSTB_IMAGE_IMPLEMENTATION=True -DSTBI_FAILURE_USERMSG=True -DCGLM_CLIPSPACE_INCLUDE_ALL=True -DCGLM_ALL_UNALIGNED=True -I./pyorama/algs -I./pyorama/libs/include -IC:/msys64/mingw64/include/python3.9 -c ./pyorama/algs/radix_sort.c -o build/temp.mingw_x86_64-3.9/./pyorama/algs/radix_sort.o -w -std=c11 -O3 -ffast-math -march=native

And here is a link to the results (.csv) as well as a graph of the data:

enter image description here

\$\endgroup\$
5
  • 2
    \$\begingroup\$ What compiler, what flags? Do you have profiling results to add to the question? \$\endgroup\$
    – Reinderien
    Commented Feb 28, 2022 at 2:20
  • \$\begingroup\$ @Reinderien The compiler is gcc version 10.3.0 (Rev5, Built by MSYS2 project). The compiler flags are a bit trickier. As mentioned in the question, I use -std=c11, -03, -ffast-math, and -march-native. I will add the exact gcc build command that is generated by my setup.py (make file equivalent) as well as the benchmarking data from my cython test to the question shortly. \$\endgroup\$ Commented Feb 28, 2022 at 8:06
  • \$\begingroup\$ @Reinderien Added a graph of the benchmark data as well as a raw csv. Basically sorting arrays of random data of various types from size n = 2^1 to n = 2^24. \$\endgroup\$ Commented Feb 28, 2022 at 8:59
  • \$\begingroup\$ Your gcc invocation doesn't make much sense. It has repeated flags, some that even conflict with each other (see march). \$\endgroup\$
    – Reinderien
    Commented Feb 28, 2022 at 19:03
  • \$\begingroup\$ That is true, but it is autogenerated when compiling a cython extension. I can specify some custom flags, which go at the end and should override the prior flags though. Is there a certain set of flags that you would recommend trying, @Reinderien? \$\endgroup\$ Commented Feb 28, 2022 at 19:17

1 Answer 1

1
\$\begingroup\$

Some general review:

We're missing an include of <string.h> for memset().

There are some constants that aren't obviously correct, such as:

uint64_t value = ((int64_t *)item)[0] + (int64_t)9223372036854775808;

It's clearer to use -INT64_MIN there, and similarly in the other functions. If we're on a 2's complement machine, we might generate faster code with ((uint64_t*)item)[0] ^ ((uint64_t)1 << 63). But I'd want to read the compiler's output to be certain.

This expression looks weird:

((uint8_t *)item + byte_offset)[0];

Surely that's the same as ((uint8_t *)item)[byte_offset]?


Performance aspects:

I think that if you want to increase performance, the next thing to do is to actually build separate sort functions for the different types, rather than indirecting through the key-functions.

You could do that by writing macros, but I think it will be easier to repeatedly compile a single source file, using preprocessor to substitute the varying parts (either by multiple includes from one translation unit, or as individual translation units - you can use Make to orchestrate that quite effectively).

I wouldn't bother trying to unroll loops by hand - these ones look amenable to intelligent unrolling by your compiler. Consider experimenting with the code-generation options that affect loop unrolling, of course.

\$\endgroup\$
1
  • \$\begingroup\$ Wow it has been a while since I looked at the code in the question! I appreciate the feedback, especially for the performance aspects. \$\endgroup\$ Commented Sep 10, 2022 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.