libguf/src/guf_sort.h

/*
    is parametrized: yes
*/

#ifndef GUF_SORT_H
#define GUF_SORT_H
#include "guf_common.h"

typedef enum guf_sort_opt {
    GUF_SORT_ASCENDING =  0, 
    GUF_SORT_DESCENDING = 1
} guf_sort_opt;

#endif

#ifdef GUF_T

#ifndef GUF_FN_NAME_PREFIX
    #define GUF_FN_NAME_PREFIX GUF_CAT(GUF_T, _arr)
#endif

#if defined(GUF_SORT_IMPL_STATIC)
    #define GUF_SORT_KWRDS static
#else
    #define GUF_SORT_KWRDS
#endif

GUF_SORT_KWRDS GUF_T *GUF_CAT(GUF_FN_NAME_PREFIX, _insertion_sort)(GUF_T *arr, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b));
GUF_SORT_KWRDS GUF_T *GUF_CAT(GUF_FN_NAME_PREFIX, _merge_sort)(GUF_T *restrict arr, GUF_T *restrict arr_tmp, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b));
GUF_SORT_KWRDS GUF_T *GUF_CAT(GUF_FN_NAME_PREFIX, _qsort)(GUF_T *arr, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b));

GUF_SORT_KWRDS bool GUF_CAT(GUF_FN_NAME_PREFIX, _is_sorted)(GUF_T *arr, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b));

#if defined(GUF_SORT_IMPL) || defined(GUF_SORT_IMPL_STATIC) 

#define guf_before(a_ptr, b_ptr) (cmp ? (sort_opt == GUF_SORT_ASCENDING ? 1 : -1) * cmp(a_ptr, b_ptr) == -1 : (sort_opt == GUF_SORT_ASCENDING) ? *(a_ptr) < *(b_ptr) : *(a_ptr) > *(b_ptr))
#define guf_before_or_equal(a_ptr, b_ptr) (cmp ? (sort_opt == GUF_SORT_ASCENDING ? 1 : -1) * cmp(a_ptr, b_ptr) <= 0 : (sort_opt == GUF_SORT_ASCENDING) ? *(a_ptr) <= *(b_ptr) : *(a_ptr) >= *(b_ptr))

/*
    Insertion sort. 
    - stable: yes
    - time: worst O(n^2); average O(n^2); best O(n) (if arr is already sorted)
    - space: O(1)
*/
GUF_SORT_KWRDS GUF_T *GUF_CAT(GUF_FN_NAME_PREFIX, _insertion_sort)(GUF_T *arr, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b))
{
    GUF_ASSERT_RELEASE(arr);
    GUF_ASSERT_RELEASE(n >= 0);
    GUF_ASSERT_RELEASE(sort_opt == GUF_SORT_ASCENDING || sort_opt == GUF_SORT_DESCENDING);
    for (ptrdiff_t i = 1; i < n; ++i) { // Range [0, i) is sorted.
        ptrdiff_t j = i; 
        while (j > 0 && guf_before(&arr[j], &arr[j - 1])) {
            GUF_SWAP(GUF_T, arr[j], arr[j - 1]);
            j = j - 1;
        }
    }
    return arr;
}

/*
    Iterative bottom-up merge-sort: cf. https://en.wikipedia.org/wiki/Merge_sort#Bottom-up_implementation (last-retrieved: 2025-01-26)
    - stable: yes
    - time: O(n * log n) (worst, average, and best)
    - space: always O(n) (for arr_tmp, allocated and freed by the caller)
*/ 
GUF_SORT_KWRDS GUF_T *GUF_CAT(GUF_FN_NAME_PREFIX, _merge_sort)(GUF_T *restrict arr, GUF_T *restrict arr_tmp, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b))
{
    GUF_ASSERT_RELEASE(arr);
    GUF_ASSERT_RELEASE(n >= 0);
    GUF_ASSERT_RELEASE(sort_opt == GUF_SORT_ASCENDING || sort_opt == GUF_SORT_DESCENDING);

    GUF_T *in = arr;
    GUF_T *out = arr_tmp;
    const size_t arr_len = n;
    for (size_t len = 1; len < arr_len; len = (len * 2 > len) ? 2 * len : SIZE_MAX) { // Subarray len 1, 2, 4, 8, ...
        for (size_t i = 0; i < arr_len; i = ((i + 2 * len) > i) ? (i + 2 * len) : SIZE_MAX) { // For each pair of subarrays of length len:
            const size_t left_begin = i; // left subarray: [left_begin, right_begin)
            const size_t right_begin = GUF_MIN(i + len, arr_len), right_end = GUF_MIN(i + 2 * len, arr_len); // right subarray [right_begin, right_end)
            size_t left_idx = left_begin, right_idx = right_begin;
            for (size_t idx = left_begin; idx < right_end; ++idx) { // Merge the left and right subarrays into arr_tmp.
                if (left_idx < right_begin && (right_idx >= right_end || guf_before_or_equal(&in[left_idx], &in[right_idx]))) {
                    out[idx] = in[left_idx++];
                } else {
                    GUF_ASSERT(right_idx < right_end);
                    out[idx] = in[right_idx++];
                }
            }  
        }
        GUF_SWAP(GUF_T*, in, out); // Avoid copying memory by switching pointers.
    }

    if (in != arr) {
        memcpy(arr, in, sizeof(GUF_T) * arr_len);
    } 
    return arr;
}

/*
    Quicksort with "Median-of-3" partitioning (non-randomised) and minimal O(log n) stack usage as described in CLRS. 
    cf. Cormen; Leiserson; Riverst; Stein (2009). "II.7 Quicksort". In "Introduction to Algorithms" (3rd ed, pp. 170-190). MIT Press 
    cf. https://en.wikipedia.org/wiki/Quicksort#Choice_of_pivot
    cf. https://en.wikipedia.org/wiki/Quicksort#Optimizations
    - stable: no
    - time: worst O(n^2); average O(n * log n); best O(n * log n)
    - space: worst O(log n) (stack space used for the recursive function calls)
*/
static GUF_T *GUF_CAT(GUF_FN_NAME_PREFIX, _qsort_range)(GUF_T *arr, ptrdiff_t first_idx, ptrdiff_t last_idx, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b))
{
    GUF_ASSERT(arr);
    GUF_ASSERT(sort_opt == GUF_SORT_ASCENDING || sort_opt == GUF_SORT_DESCENDING);
        
    while (first_idx >= 0 && first_idx < last_idx) {
        // 1.) Partition using "Median-of-3 partitioning", cf. https://en.wikipedia.org/wiki/Quicksort#Choice_of_pivot (last-retrieved 2025-01-30)
        const ptrdiff_t mid_idx = first_idx + ((last_idx - first_idx) / 2); // Equivalent to (first_idx + last_idx) / 2, cf. https://research.google/blog/extra-extra-read-all-about-it-nearly-all-binary-searches-and-mergesorts-are-broken/ (last-retrieved 2025-01-24)
        GUF_ASSERT(mid_idx >= 0 && mid_idx <= last_idx);
        ptrdiff_t pivot_idx = last_idx;
        if (guf_before(&arr[mid_idx], &arr[first_idx])) { 
            GUF_SWAP(GUF_T, arr[first_idx], arr[mid_idx]);
        } 
        if (guf_before(&arr[last_idx], &arr[first_idx])) {
            GUF_SWAP(GUF_T, arr[first_idx], arr[last_idx]);
        }
        if (guf_before(&arr[mid_idx], &arr[last_idx])) {
            GUF_SWAP(GUF_T, arr[mid_idx], arr[last_idx]);
        }
        ptrdiff_t i = first_idx - 1; // i: end idx of the subarray with elements <= pivot
        for (ptrdiff_t j = first_idx; j < last_idx; ++j) { // j: end idx of the subarray with elements > pivot
            if (guf_before_or_equal(&arr[j], &arr[pivot_idx])) { 
                ++i;
                GUF_ASSERT(i >= 0 && i < last_idx);
                GUF_SWAP(GUF_T, arr[i], arr[j]);
            }
        }
        GUF_ASSERT(i + 1 >= 0 && i + 1 < last_idx);
        GUF_SWAP(GUF_T, arr[i + 1], arr[last_idx]); // The pivot element is always <= itself (or >= for descending sorts).
        pivot_idx = i + 1;

        // 2.) Sort the two partitions [first_idx, pivot_idx) and (pivot_idx, last_idx] recursively.
        /* 
            "To make sure at most O(log n) space is used, recur first into the smaller side of the partition, 
            then use a tail call to recur into the other, or update the parameters to no longer include the now sorted smaller side, 
            and iterate to sort the larger side.", cf. https://en.wikipedia.org/wiki/Quicksort#Optimizations (last-retrieved 2025-01-25)
            
            (Solution to exercise "II.7.4c. Stack depth for quicksort". In "Introduction to Algorithms" (3rd ed, p. 188))
        */
        if (pivot_idx <= mid_idx) { // a.) Left subarray is smaller or equal than the right subarray -> recur into the smaller left subarray first.
            GUF_CAT(GUF_FN_NAME_PREFIX, _qsort_range)(arr, first_idx, pivot_idx - 1, sort_opt, cmp); 
            first_idx = pivot_idx + 1;
        } else { // b.) Right subarray is smaller than the left subarray -> recur into the smaller right subarray first.
            GUF_CAT(GUF_FN_NAME_PREFIX, _qsort_range)(arr, pivot_idx + 1, last_idx, sort_opt, cmp); 
            last_idx = pivot_idx - 1;
        }
    }
    return arr;
}

GUF_SORT_KWRDS GUF_T *GUF_CAT(GUF_FN_NAME_PREFIX, _qsort)(GUF_T *arr, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b))
{
    GUF_ASSERT_RELEASE(arr);
    GUF_ASSERT_RELEASE(sort_opt == GUF_SORT_ASCENDING || sort_opt == GUF_SORT_DESCENDING);
    GUF_ASSERT_RELEASE(n >= 0);
    if (n <= 1) {
        return arr;
    } else if (n <= 4) {
        return GUF_CAT(GUF_FN_NAME_PREFIX, _insertion_sort)(arr, n, sort_opt, cmp);
    } else {
        return GUF_CAT(GUF_FN_NAME_PREFIX, _qsort_range)(arr, 0, n - 1, sort_opt, cmp);
    }
}

GUF_SORT_KWRDS bool GUF_CAT(GUF_FN_NAME_PREFIX, _is_sorted)(GUF_T *arr, ptrdiff_t n, guf_sort_opt sort_opt, int(*cmp)(const GUF_T *a, const GUF_T *b))
{
    GUF_ASSERT_RELEASE(arr);
    for (ptrdiff_t i = 0; i < n - 1; ++i) {
        if (!guf_before_or_equal(arr + i, arr + (i+1))) {
            return false;
        }
    }
    return true;
}

#undef guf_before
#undef guf_before_or_equal

#undef GUF_SORT_IMPL
#undef GUF_SORT_IMPL_STATIC
#endif /* end #ifdef GUF_IMPL */

#undef GUF_SORT_KWRDS
#undef GUF_T
#undef GUF_FN_NAME_PREFIX

#endif /* end  #ifdef GUF_T */