Compile Time Recursive Sort with Selection Algorithm

Question

Is it possible to write a compile time recursive sort function in c++ using the selection algorithm?

I want to sort array x from element istart to element iend in ascending order. Array x has N elements. The data in the input array x is known only at runtime, hence the data can only be sorted at runtime. However, I want to generate the C++ code, i.e. all recursive function calls to sort_asc() at compile time. Furthermore, I want to use this code in a CUDA device function. Since CUDA is a subset of C++ with a handful of extensions, I thought this was the right place to ask. Unfortunately, I don't think CUDA supports the constexpr keyword, neither Boost nor STL.

I came up with the following code for sorting in ascending order.

// sort in ascending order.
template< int istart, int N, int iend, int iend_min_istart >
inline void sort_asc
(
    float *x
)
{
    int   min_idx = istart;
    float min_val = x[ min_idx ];
    #pragma unroll
    for( int i=istart+1; i<N; i++ ){
        if( x[ i ] < min_val ){
            min_idx = i;
            min_val = x[ i ];
        }
    }
    swap( x[ istart ], x[ min_idx ] );
    sort_asc< istart+1, N, iend, iend-(istart+1) >( x );
}   

where iend_min_istart = iend - istart. If iend_min_istart < 0 then the recursion has finished, thus we can write the termination condition as:

// sort recursion termination condition.
template< int istart, int N, int iend > 
inline void sort_asc< -1 >
(
    float *x
)
{
    // do nothing.
}

The swap function is defined as:

void d_swap
( 
    float &a, 
    float &b 
)
{
    float c = a;
    a = b;
    b = c;
}

Then I'd call the sort function as:

void main( int argc, char *argv[] )
{
    float x[] = { 3, 4, 9, 2, 7 };   // 5 element array.
    sort_asc< 0, 5, 2, 2-0 >( x );   // sort from the 1st till the 3th element.
    float x_median = cost[ 2 ];      // the 3th element is the median
}

However, this code does not compile since c++ does not support partial template specialization for functions. Furthermore, I don't know how to write this in C++ meta programming. Is there any way to make this code work?

Answer 1

Using selection sort this time. (Variant using merge sort.) Without any guarantees, but passes a simple test:

// generate a sequence of integers as non-type template arguments
// (a basic meta-programming tool)
template<int... Is> struct seq {};
template<int N, int... Is> struct gen_seq : gen_seq<N-1, N-1, Is...> {};
template<int... Is> struct gen_seq<0, Is...> : seq<Is...> {};


// an array type that can be returned from a function
// and has `constexpr` accessors (as opposed to C++11's `std::array`)
template<class T, int N>
struct c_array
{
    T m[N];

    constexpr T const& operator[](int p) const
    {  return m[p];  }

    constexpr T const* begin() const { return m+0; }
    constexpr T const* end() const { return m+N; }
};


// return the index of the smallest element
template<class T, int Size>
constexpr int c_min_index(c_array<T, Size> const& arr, int offset, int cur)
{
    return Size == offset ? cur :
           c_min_index(arr, offset+1, arr[cur] < arr[offset] ? cur : offset);
}

Our goal is to be able to sort at compile-time. We can use several tools for this (macros, class templates, constexpr functions), but constexpr functions tend to be the simplest in many cases.

The function c_min_index above is an example for the C++11 restrictions on constexpr functions: It can only consist of a single return statement (plus static_assert and using), and at least for some possible arguments must yield a constant expression. This implies: No loops, no assignment. So we need to use recursion here, and pass the intermediary result to the next invocation (cur).

// copy the array but with the elements at `index0` and `index1` swapped
template<class T, int Size, int... Is>
constexpr c_array<T, Size> c_swap(c_array<T, Size> const& arr,
                                  int index0, int index1, seq<Is...>)
{
    return {{arr[Is == index0 ? index1 : Is == index1 ? index0 : Is]...}};
}

As we cannot modify the argument due to constexpr function restrictions, we need to return a copy of the whole array just to swap two elements.

The function expects the Is... to be a sequence of integers 0, 1, 2 .. Size-1, as produced by a gen_seq<Size>{} (in its base class). Those integers are used to access the elements of the array arr. Some expression like arr[Is]... would produce arr[0], arr[1], arr[2] .. arr[Size-1]. Here, we swap the index by applying a conditional to the current integer: a[0 == index0 ? index1 : 0 == index1 ? index0 : 0], ..

// the selection sort algorithm
template<class T, int Size>
constexpr c_array<T, Size> c_sel_sort(c_array<T, Size> const& arr, int cur = 0)
{
    return cur == Size ? arr :
           c_sel_sort( c_swap(arr, cur, c_min_index(arr, cur, cur),
                              gen_seq<Size>{}),
                       cur+1 );
}

Again, we need to replace loop with recursion (cur). The rest is a straight-forward selection sort implementation with an awkward syntax: Starting from the index cur, we find the minimal element in the remaining array, and swap it with the element at cur. Then, we re-run the selection sort on the remaining unsorted array (by incrementing cur).

#include <iostream>
#include <iterator>

int main()
{
    // homework: write a wrapper so that C-style arrays can be passed
    //           to an overload of `c_sel_sort` ;)
    constexpr c_array<float,10> f = {{4, 7, 9, 0, 6, 2, 3, 8, 1, 5}};
    constexpr auto sorted = c_sel_sort(f);
    for(auto const& e : sorted) std::cout << e << ", ";
}

Answer 2

I came up with the following meta code and it works.

template < class T >
__forceinline void swap
( 
    T &a, 
    T &b 
)
{
    T c = a;
    a = b;
    b = c;
}

template< int istart, int N, int iend, int iend_min_istart >
class SORT
{
    public:
        static __forceinline void sort( float *x, int *idx )
        {
            // sort code.
            int   min_idx = istart;
            float min_val = x[ min_idx ];
            for( int i=istart+1; i<N; i++ ){
                if( x[ i ] < min_val ){
                    min_idx = i;
                    min_val = x[ i ];
                }
            }
            swap(   x[ istart ],   x[ min_idx ] );
            swap( idx[ istart ], idx[ min_idx ] );
            SORT<istart+1, N, iend, iend-(istart+1)>::sort( x, idx );
        }
};

template< int istart, int N, int iend >
class SORT<istart,N,iend,-1>
{
    public:
        static __forceinline void sort( float *x, int *idx )
        {

        }
};

void main( int argc, char *argv[] )
{
    float arr[] = {1,4,2,7,5};
    int   idx[] = {0,1,2,3,4};
    SORT<0,5,3,2-0>::sort( arr, idx );
}