tr-sort.hpp

Go to the documentation of this file.

 
#ifndef COM_SAXBOPHONE_TR_SORT_HPP
#define COM_SAXBOPHONE_TR_SORT_HPP
 
#include <cmath>       // ceil, nextafter, pow
#include <cstddef>
#include <span>
#include <type_traits>
#include <vector>
 
 
namespace com::saxbophone::tr_sort {
    template <
        typename T,
        std::size_t Extent = std::dynamic_extent,
        typename Real = long double
    >
    void sort(std::span<T, Extent> data) {
        // don't sort data of length {0..1}
        if (data.size() < 2) {
            return;
        }
        // gather stats on first pass of data
        std::size_t size = data.size();
        T min, max;
        if constexpr (std::is_floating_point<T>::value) {
            min = +std::numeric_limits<T>::infinity();
            max = -std::numeric_limits<T>::infinity();
        } else {
            min = std::numeric_limits<T>::max();
            max = std::numeric_limits<T>::lowest();
        }
        // Real mean = 0.0;
        // Real mid = 0.0;
        // must be Real because otherwise gives out-of-range values for signed types
        Real range = 0;
        T previous = data[0];
        bool already_sorted = true;
        for (auto datum : data) {
            // mean += datum;
            if (datum < min) {
                min = datum;
            }
            // NOTE: >= to preserve stable sorting, so last max value gets put at end
            /*
             * Not actually an issue for this implementation because it only deals
             * with primitives, however if extended to handle general-purpose types
             * (using some casting), this will matter.
             */
            if (datum >= max) {
                max = datum;
            }
            if (datum < previous) {
                already_sorted = false;
            }
            previous = datum;
        }
        // short cut for already sorted
        if (already_sorted) {
            return;
        }
        // short cut for size 2
        if (data.size() == 2) {
            data[0] = min;
            data[data.size() - 1] = max;
            return;
        }
        if constexpr (std::is_floating_point<T>::value) {
            if (min == -std::numeric_limits<T>::infinity()) {
                min = std::numeric_limits<T>::lowest();
            }
            if (max == +std::numeric_limits<T>::infinity()) {
                max = std::numeric_limits<T>::max();
            }
        }
        // mean /= size;
        // mid = (min + max) / 2.0;
        range = (Real)max - (Real)min;
        // temporary storage for sorting -- vector of sub-vectors to store partial sorts
        std::vector<std::vector<T>> sorts(data.size() + 2);
        for (auto n : data) {
            if constexpr (std::is_floating_point<T>::value) {
                if (n == -std::numeric_limits<T>::infinity()) {
                    sorts[0].push_back(n);
                    continue;
                } else if (n == +std::numeric_limits<T>::infinity()) {
                    sorts[sorts.size() - 1].push_back(n);
                    continue;
                }
            }
            // calculated sort position
            Real raw_pos = std::ceil((((Real)n - min) / range) * (size - 1));
            std::size_t pos = std::isnan(raw_pos) ? 1u : (std::size_t)raw_pos;
            if (raw_pos < 0) {
                pos = 0;
            } else if (raw_pos > (sorts.size() - 1)) {
                pos = sorts.size() - 1;
            }
            pos += 1;
            sorts[pos].push_back(n);
        }
        // pull data out of sorted buckets, recursively sorting each before pulling
        std::size_t i = 0;
        for (auto& bucket : sorts) {
            // recursively sort any sort buckets that are larger than 1
            if (bucket.size() > 1) {
                sort<T, std::dynamic_extent, Real>(bucket);
            }
            for (T datum : bucket) {
                data[i] = datum;
                i++;
            }
        }
        return;
    }
}
 
#endif // include guard