gravity_sort_benchmark.cpp

#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <fstream>
#include <atomic>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <sstream>
#include <random>
#include <string>
#include <thread>
#include <vector>

using Clock = std::chrono::steady_clock;
using ms    = std::chrono::duration<double, std::milli>;

struct Timer {
    Clock::time_point start;
    Timer() : start(Clock::now()) {}
    double elapsed_ms() const {
        return std::chrono::duration_cast<ms>(Clock::now() - start).count();
    }
};

struct TrialStats {
    double mean_ms   = 0;
    double median_ms = 0;
    double min_ms    = 0;
    double max_ms    = 0;
    double stddev_ms = 0;
    bool   correct   = true;
    int    reps      = 0;
};

// ---------------------------------------------------------------------------
// Synthetic dataset generators
// ---------------------------------------------------------------------------

static std::vector<std::uint32_t> generate_random(std::size_t n,
    std::uint32_t max_value = std::numeric_limits<std::uint32_t>::max())
{
    std::mt19937_64 rng(123456789ULL + n);
    std::uniform_int_distribution<std::uint32_t> dist(0, max_value);
    std::vector<std::uint32_t> v;
    v.reserve(n);
    for (std::size_t i = 0; i < n; ++i) v.push_back(dist(rng));
    return v;
}

static std::vector<std::uint32_t> generate_clustered(std::size_t n,
    std::size_t centers = 32, std::uint32_t spread = 128)
{
    std::mt19937_64 rng(987654321ULL + n);
    std::uniform_int_distribution<std::uint32_t> center_dist(0,
        std::numeric_limits<std::uint32_t>::max());
    std::uniform_int_distribution<std::size_t> pick_center(0, centers - 1);
    std::uniform_int_distribution<int> offset_dist(
        -static_cast<int>(spread), static_cast<int>(spread));

    std::vector<std::uint32_t> cv;
    cv.reserve(centers);
    for (std::size_t i = 0; i < centers; ++i) cv.push_back(center_dist(rng));

    std::vector<std::uint32_t> v;
    v.reserve(n);
    for (std::size_t i = 0; i < n; ++i) {
        long long x = static_cast<long long>(cv[pick_center(rng)])
                    + static_cast<long long>(offset_dist(rng));
        if (x < 0) x = 0;
        if (x > static_cast<long long>(std::numeric_limits<std::uint32_t>::max()))
            x = std::numeric_limits<std::uint32_t>::max();
        v.push_back(static_cast<std::uint32_t>(x));
    }
    return v;
}

static std::vector<std::uint32_t> generate_nearly_sorted(std::size_t n,
    std::size_t swaps = 64)
{
    std::mt19937_64 rng(4242424242ULL + n);
    std::vector<std::uint32_t> v(n);
    std::iota(v.begin(), v.end(), 0u);
    if (n == 0) return v;
    std::uniform_int_distribution<std::size_t> dist(0, n - 1);
    for (std::size_t i = 0; i < std::min(swaps, n / 2); ++i) {
        std::swap(v[dist(rng)], v[dist(rng)]);
    }
    return v;
}

static std::vector<std::uint32_t> generate_reverse_sorted(std::size_t n) {
    std::vector<std::uint32_t> v(n);
    for (std::size_t i = 0; i < n; ++i)
        v[i] = static_cast<std::uint32_t>(n - 1 - i);
    return v;
}

// Power-law: most values small, a few very large — common in real-world data
// (web traffic, file sizes, word frequencies, genomic hotspots)
static std::vector<std::uint32_t> generate_power_law(std::size_t n,
    double exponent = 1.5)
{
    std::mt19937_64 rng(1111111111ULL + n);
    std::uniform_real_distribution<double> unit(0.0, 1.0);
    std::vector<std::uint32_t> v;
    v.reserve(n);
    const double max32 = static_cast<double>(std::numeric_limits<std::uint32_t>::max());
    for (std::size_t i = 0; i < n; ++i) {
        // Inverse-CDF of a power-law: x = (1 - u)^(-1/(alpha-1))
        double u   = unit(rng);
        double raw = std::pow(1.0 - u + 1e-12, -1.0 / (exponent - 1.0));
        // Normalise into uint32 range; raw starts at 1 and grows without bound.
        // We clip at a reasonable multiple so ~99% of values fit spread across range.
        double clipped = std::min(raw / 1000.0, 1.0) * max32;
        v.push_back(static_cast<std::uint32_t>(clipped));
    }
    return v;
}

// Pipe-organ: 0,1,2,...,n/2-1,n/2-1,...,1,0 — stress-tests adaptive sorts
static std::vector<std::uint32_t> generate_pipe_organ(std::size_t n) {
    std::vector<std::uint32_t> v(n);
    std::size_t half = n / 2;
    for (std::size_t i = 0; i < half; ++i)
        v[i] = static_cast<std::uint32_t>(i);
    for (std::size_t i = half; i < n; ++i)
        v[i] = static_cast<std::uint32_t>(n - 1 - i);
    return v;
}

// All-equal: worst-case for algorithms that don't short-circuit on equal keys
static std::vector<std::uint32_t> generate_all_equal(std::size_t n) {
    return std::vector<std::uint32_t>(n, 42u);
}

// Already sorted: best-case for adaptive sorts, measures overhead baseline
static std::vector<std::uint32_t> generate_already_sorted(std::size_t n) {
    std::vector<std::uint32_t> v(n);
    std::iota(v.begin(), v.end(), 0u);
    return v;
}

// Bimodal: two well-separated clusters with a clean gap between them
// Perfect stress-test for Shock Wave Sort / Immiscible Fluids
static std::vector<std::uint32_t> generate_bimodal(std::size_t n) {
    std::mt19937_64 rng(7777777777ULL + n);
    std::uniform_int_distribution<std::uint32_t> lo(0,       100000);
    std::uniform_int_distribution<std::uint32_t> hi(3000000, 3100000);
    std::vector<std::uint32_t> v;
    v.reserve(n);
    for (std::size_t i = 0; i < n; ++i)
        v.push_back((i % 2 == 0) ? lo(rng) : hi(rng));
    // Shuffle to prevent order-based shortcuts
    std::shuffle(v.begin(), v.end(), rng);
    return v;
}

// Small range: all values in [0, 65535] — stresses algorithms that should detect
// the effective bit-width of the data (Isotope Separation Sort target)
static std::vector<std::uint32_t> generate_small_range(std::size_t n) {
    std::mt19937_64 rng(3333333333ULL + n);
    std::uniform_int_distribution<std::uint32_t> dist(0, 65535);
    std::vector<std::uint32_t> v;
    v.reserve(n);
    for (std::size_t i = 0; i < n; ++i) v.push_back(dist(rng));
    return v;
}

// Many sorted runs: K contiguous sorted segments concatenated — nearly-sorted
// variant that has more structure than random swaps; good for Phase Transition
static std::vector<std::uint32_t> generate_many_runs(std::size_t n,
    std::size_t num_runs = 128)
{
    std::mt19937_64 rng(5555555555ULL + n);
    std::uniform_int_distribution<std::uint32_t> dist(0,
        std::numeric_limits<std::uint32_t>::max());
    std::vector<std::uint32_t> v;
    v.reserve(n);
    std::size_t run_len = std::max<std::size_t>(1, n / num_runs);
    for (std::size_t r = 0; r < num_runs; ++r) {
        std::size_t start = v.size();
        std::size_t end   = std::min(start + run_len, n);
        while (v.size() < end) v.push_back(dist(rng));
        std::sort(v.begin() + static_cast<std::ptrdiff_t>(start), v.end());
    }
    while (v.size() < n) v.push_back(dist(rng));
    return v;
}

// High duplicates: only K distinct values across all n elements
// Pushes sorting to do mostly equality comparisons
static std::vector<std::uint32_t> generate_high_duplicates(std::size_t n,
    std::size_t distinct = 16)
{
    std::mt19937_64 rng(9999999999ULL + n);
    std::uniform_int_distribution<std::size_t> pick(0, distinct - 1);
    // Fixed distinct values spread across uint32 range
    std::vector<std::uint32_t> keys(distinct);
    for (std::size_t i = 0; i < distinct; ++i)
        keys[i] = static_cast<std::uint32_t>((static_cast<std::uint64_t>(i) *
                  std::numeric_limits<std::uint32_t>::max()) / (distinct - 1));
    std::vector<std::uint32_t> v;
    v.reserve(n);
    for (std::size_t i = 0; i < n; ++i) v.push_back(keys[pick(rng)]);
    return v;
}

// ---------------------------------------------------------------------------
// Real-file dataset loaders
// Reads raw bytes from a file and reinterprets every 4-byte block as a
// uint32 (little-endian host order).  This is fast and format-agnostic:
// DNA, XML, source code — everything becomes an integer sort workload.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> load_file_as_uint32(
    const std::string& path, std::size_t max_elements = 0)
{
    std::ifstream f(path, std::ios::binary | std::ios::ate);
    if (!f) return {};
    std::streamsize bytes = f.tellg();
    if (bytes <= 0) return {};
    f.seekg(0);

    std::size_t n_elems = static_cast<std::size_t>(bytes) / sizeof(std::uint32_t);
    if (max_elements > 0 && n_elems > max_elements) n_elems = max_elements;

    std::vector<std::uint32_t> v(n_elems);
    f.read(reinterpret_cast<char*>(v.data()), n_elems * sizeof(std::uint32_t));
    return v;
}

// ---------------------------------------------------------------------------
// Shared helpers
// ---------------------------------------------------------------------------

static void insertion_sort(std::vector<std::uint32_t>& a) {
    for (std::size_t i = 1; i < a.size(); ++i) {
        std::uint32_t x = a[i];
        std::size_t j = i;
        while (j > 0 && a[j - 1] > x) { a[j] = a[j - 1]; --j; }
        a[j] = x;
    }
}

// Sort each bucket in parallel using a work-stealing counter.
static void parallel_sort_buckets(std::vector<std::vector<std::uint32_t>>& buckets,
                                  unsigned threads)
{
    std::atomic<std::size_t> next{0};
    std::vector<std::thread> pool;
    pool.reserve(threads);
    auto worker = [&]() {
        for (;;) {
            std::size_t b = next.fetch_add(1, std::memory_order_relaxed);
            if (b >= buckets.size()) break;
            auto& bkt = buckets[b];
            if (bkt.size() <= 32) insertion_sort(bkt);
            else                  std::sort(bkt.begin(), bkt.end());
        }
    };
    for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker);
    for (auto& th : pool) th.join();
}

// ---------------------------------------------------------------------------
// Algorithm 1 (existing): Adaptive gravity bucket sort
//   - Samples the data to build quantile-based splitters
//   - O(log B) per element (binary search into splitter table)
//   - Good when the distribution is truly unknown and non-monotone
// ---------------------------------------------------------------------------

static std::uint32_t percentile_from_sample(std::vector<std::uint32_t>& s, double p) {
    if (s.empty()) return 0;
    std::sort(s.begin(), s.end());
    std::size_t i = static_cast<std::size_t>(p * (s.size() - 1));
    return s[std::min(i, s.size() - 1)];
}

static std::vector<std::uint32_t> choose_splitters(
    const std::vector<std::uint32_t>& arr, std::size_t B)
{
    if (arr.empty() || B <= 1) return {};
    std::size_t ss = std::min<std::size_t>(arr.size(), std::max<std::size_t>(1024, B * 32));
    std::vector<std::uint32_t> sample;
    sample.reserve(ss);
    if (ss == arr.size()) {
        sample = arr;
    } else {
        std::mt19937_64 rng(7777777ULL + arr.size());
        std::uniform_int_distribution<std::size_t> dist(0, arr.size() - 1);
        for (std::size_t i = 0; i < ss; ++i) sample.push_back(arr[dist(rng)]);
    }
    std::sort(sample.begin(), sample.end());
    std::vector<std::uint32_t> sp;
    sp.reserve(B - 1);
    for (std::size_t i = 1; i < B; ++i) {
        std::size_t idx = (i * sample.size()) / B;
        sp.push_back(sample[std::min(idx, sample.size() - 1)]);
    }
    sp.erase(std::unique(sp.begin(), sp.end()), sp.end());
    return sp;
}

static std::vector<std::uint32_t> gravity_bucket_sort_parallel(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    std::uint32_t minv = *std::min_element(arr.begin(), arr.end());
    std::uint32_t maxv = *std::max_element(arr.begin(), arr.end());
    if (minv == maxv) return arr;

    double range = static_cast<double>(maxv) - static_cast<double>(minv);
    double iqr_ratio = 0.0;
    {
        std::size_t ss = std::min<std::size_t>(arr.size(), 4096);
        std::vector<std::uint32_t> s;
        s.reserve(ss);
        if (ss == arr.size()) { s = arr; }
        else {
            std::mt19937_64 rng(1357911ULL + arr.size());
            std::uniform_int_distribution<std::size_t> d(0, arr.size() - 1);
            for (std::size_t i = 0; i < ss; ++i) s.push_back(arr[d(rng)]);
        }
        std::uint32_t q1 = percentile_from_sample(s, 0.25);
        std::uint32_t q3 = percentile_from_sample(s, 0.75);
        iqr_ratio = (range > 0) ? ((double)(q3 - q1) / (range + 1.0)) : 0.0;
    }

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)arr.size()) * 1.25);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    if (iqr_ratio < 0.05)      B = std::max<std::size_t>(32, B / 2);
    else if (iqr_ratio > 0.30) B = std::min<std::size_t>(4096, B * 2);
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, arr.size() / 8));

    auto splitters = choose_splitters(arr, B);
    if (splitters.empty()) return arr;
    std::size_t bucket_total = splitters.size() + 1;

    unsigned threads = threads_hint == 0 ? 4u : threads_hint;
    threads = std::max(1u, std::min<unsigned>(threads,
                            static_cast<unsigned>(arr.size())));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(bucket_total));

    {
        std::vector<std::thread> pool;
        pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t begin = (arr.size() * tid) / threads;
            std::size_t end   = (arr.size() * (tid + 1)) / threads;
            for (std::size_t i = begin; i < end; ++i) {
                std::uint32_t x = arr[i];
                auto it = std::upper_bound(splitters.begin(), splitters.end(), x);
                local[tid][static_cast<std::size_t>(it - splitters.begin())].push_back(x);
            }
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(bucket_total);
    for (std::size_t b = 0; b < bucket_total; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets, threads);

    std::vector<std::uint32_t> out;
    out.reserve(arr.size());
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

// ---------------------------------------------------------------------------
// Algorithm 2 (NEW): Gravity physics sort — fall-time log bucket mapping
//
// Physical model realised in code:
//   A particle of value x falls at speed proportional to x (heavier = faster,
//   like terminal velocity in a fluid).  Its landing time is H/x, which is
//   inversely proportional to x.  Equal time intervals between bucket
//   boundaries therefore correspond to *logarithmically spaced* value
//   boundaries, because:
//
//       t = H / x  →  x = H / t  →  log(x) is linear in log(t)
//
// Practical payoff:
//   • O(1) bucket assignment per element — one log + multiply, no binary search.
//   • No sampling pass needed — only min and max are required.
//   • Naturally handles power-law / right-skewed data (common in real files):
//     linear buckets pile everything into the first few; log buckets spread it
//     evenly so each bucket's local sort is O((n/B) log(n/B)) not O(n log n).
//   • Falls back to linear spacing when range is small (avoids log(0) edge case).
// ---------------------------------------------------------------------------

static std::vector<std::uint32_t> gravity_physics_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    // Single parallel pass: find min and max.
    std::uint32_t minv = *std::min_element(arr.begin(), arr.end());
    std::uint32_t maxv = *std::max_element(arr.begin(), arr.end());
    if (minv == maxv) return arr;

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)arr.size()) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, arr.size() / 4));

    // Precompute the log-space scale factor.
    // bucket(x) = floor( log1p(x - min) / log1p(max - min) * B )
    // log1p avoids log(0) when x == min (gives bucket 0) or when range is 1.
    const double log_range    = std::log1p(static_cast<double>(maxv - minv));
    const double inv_log_range = (log_range > 0.0) ? (1.0 / log_range) : 1.0;
    const double scale         = inv_log_range * static_cast<double>(B);

    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        double pos = std::log1p(static_cast<double>(x - minv)) * scale;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    unsigned threads = threads_hint == 0 ? 4u : threads_hint;
    threads = std::max(1u, std::min<unsigned>(threads,
                            static_cast<unsigned>(arr.size())));

    // Each thread scatters its slice into thread-local buckets to avoid locking.
    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));

    {
        std::vector<std::thread> pool;
        pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t begin = (arr.size() * tid) / threads;
            std::size_t end   = (arr.size() * (tid + 1)) / threads;
            for (std::size_t i = begin; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    // Merge thread-local buckets.
    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets, threads);

    std::vector<std::uint32_t> out;
    out.reserve(arr.size());
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

// ---------------------------------------------------------------------------
// Algorithm 3 (existing): Parallel 16-bit radix sort
// ---------------------------------------------------------------------------

static std::vector<std::uint32_t> parallel_radix_sort_16(
    const std::vector<std::uint32_t>& input,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (input.size() <= 1) return input;

    std::vector<std::uint32_t> a = input;
    std::vector<std::uint32_t> out(a.size());

    unsigned threads = threads_hint == 0 ? 4u : threads_hint;
    threads = std::max(1u, std::min<unsigned>(threads,
                            static_cast<unsigned>(a.size())));

    const std::uint32_t base = 1u << 16;
    const std::uint32_t mask = base - 1u;
    std::uint32_t maxv = *std::max_element(a.begin(), a.end());

    std::size_t passes = 0;
    for (std::uint32_t x = maxv; x > 0; x >>= 16) ++passes;
    if (passes == 0) passes = 1;

    std::vector<std::vector<std::uint32_t>> local_counts(threads,
        std::vector<std::uint32_t>(base, 0));
    std::vector<std::vector<std::uint32_t>> thread_starts(threads,
        std::vector<std::uint32_t>(base, 0));
    std::vector<std::uint32_t> global_counts(base, 0);

    for (std::size_t pass = 0; pass < passes; ++pass) {
        std::size_t shift = pass * 16;
        std::fill(global_counts.begin(), global_counts.end(), 0);
        for (unsigned t = 0; t < threads; ++t)
            std::fill(local_counts[t].begin(), local_counts[t].end(), 0);

        {
            std::vector<std::thread> pool;
            pool.reserve(threads);
            auto count_worker = [&](unsigned tid) {
                std::size_t begin = (a.size() * tid) / threads;
                std::size_t end   = (a.size() * (tid + 1)) / threads;
                auto& lc = local_counts[tid];
                for (std::size_t i = begin; i < end; ++i)
                    ++lc[(a[i] >> shift) & mask];
            };
            for (unsigned t = 0; t < threads; ++t) pool.emplace_back(count_worker, t);
            for (auto& th : pool) th.join();
        }

        for (std::uint32_t d = 0; d < base; ++d) {
            std::uint32_t sum = 0;
            for (unsigned t = 0; t < threads; ++t) sum += local_counts[t][d];
            global_counts[d] = sum;
        }

        std::uint32_t total = 0;
        std::vector<std::uint32_t> prefix(base, 0);
        for (std::uint32_t d = 0; d < base; ++d) { prefix[d] = total; total += global_counts[d]; }

        for (std::uint32_t d = 0; d < base; ++d) {
            std::uint32_t running = prefix[d];
            for (unsigned t = 0; t < threads; ++t) {
                thread_starts[t][d] = running;
                running += local_counts[t][d];
            }
        }

        {
            std::vector<std::thread> pool;
            pool.reserve(threads);
            auto scatter_worker = [&](unsigned tid) {
                std::size_t begin = (a.size() * tid) / threads;
                std::size_t end   = (a.size() * (tid + 1)) / threads;
                auto starts = thread_starts[tid];
                for (std::size_t i = begin; i < end; ++i) {
                    std::uint32_t x = a[i];
                    std::uint32_t d = (x >> shift) & mask;
                    out[starts[d]++] = x;
                }
            };
            for (unsigned t = 0; t < threads; ++t) pool.emplace_back(scatter_worker, t);
            for (auto& th : pool) th.join();
        }

        a.swap(out);
    }
    return a;
}

// ---------------------------------------------------------------------------
// Dead-Weight Sort (Prefix Compression Radix)
// Physics: leading zero bits are drag that doesn't affect ordering but wastes
// radix passes. Sort by (x - min) using digit_bits = ceil(effective_bits/2),
// so count arrays fit in L1 when the value range is small.
//
// For range <= 2^16: 1 pass of 16 bits (same as parallel_radix)
// For range in (2^16, 2^20]: 2 passes of 10 bits  (1024-entry count, L1-friendly)
// For range in (2^20, 2^24]: 2 passes of 12 bits  (4096-entry count, L1-friendly)
// For range in (2^24, 2^32]: 2 passes of 16 bits  (same as parallel_radix)
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> dead_weight_sort(
    const std::vector<std::uint32_t>& input,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    const std::size_t n = input.size();
    if (n <= 1) return input;

    unsigned threads = threads_hint == 0 ? 4u : threads_hint;
    threads = std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(n)));

    std::uint32_t minv = input[0], maxv = input[0];
    for (auto x : input) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
    }

    const std::uint32_t range = maxv - minv;
    if (range == 0) return input;  // all equal — already sorted

    // Compute effective bit-width: bits in range
    int effective_bits = 1;
    if (range > 0) effective_bits = 32 - __builtin_clz(range);

    // Choose digit width: smallest that covers effective_bits in at most 2 passes
    int digit_bits;
    int passes;
    if (effective_bits <= 16) {
        digit_bits = 16; passes = 1;
    } else if (effective_bits <= 20) {
        digit_bits = 10; passes = 2;
    } else if (effective_bits <= 24) {
        digit_bits = 12; passes = 2;
    } else {
        digit_bits = 16; passes = 2;
    }

    const std::uint32_t base = 1u << digit_bits;
    const std::uint32_t mask = base - 1u;

    // Work on shifted values: sort (x - minv) so only effective_bits matter
    std::vector<std::uint32_t> a(n), b(n);
    for (std::size_t i = 0; i < n; ++i) a[i] = input[i] - minv;

    std::vector<std::vector<std::uint32_t>> local_counts(threads,
        std::vector<std::uint32_t>(base, 0));
    std::vector<std::vector<std::uint32_t>> thread_starts(threads,
        std::vector<std::uint32_t>(base, 0));
    std::vector<std::uint32_t> global_counts(base, 0);

    for (int pass = 0; pass < passes; ++pass) {
        int shift = pass * digit_bits;
        std::fill(global_counts.begin(), global_counts.end(), 0);
        for (unsigned t = 0; t < threads; ++t)
            std::fill(local_counts[t].begin(), local_counts[t].end(), 0);

        {
            std::vector<std::thread> pool; pool.reserve(threads);
            auto count_worker = [&](unsigned tid) {
                std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
                auto& lc = local_counts[tid];
                for (std::size_t i = beg; i < end; ++i) ++lc[(a[i] >> shift) & mask];
            };
            for (unsigned t = 0; t < threads; ++t) pool.emplace_back(count_worker, t);
            for (auto& th : pool) th.join();
        }

        for (std::uint32_t d = 0; d < base; ++d) {
            std::uint32_t sum = 0;
            for (unsigned t = 0; t < threads; ++t) sum += local_counts[t][d];
            global_counts[d] = sum;
        }

        std::uint32_t total = 0;
        std::vector<std::uint32_t> prefix(base, 0);
        for (std::uint32_t d = 0; d < base; ++d) { prefix[d] = total; total += global_counts[d]; }

        for (std::uint32_t d = 0; d < base; ++d) {
            std::uint32_t running = prefix[d];
            for (unsigned t = 0; t < threads; ++t) {
                thread_starts[t][d] = running;
                running += local_counts[t][d];
            }
        }

        {
            std::vector<std::thread> pool; pool.reserve(threads);
            auto scatter_worker = [&](unsigned tid) {
                std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
                auto starts = thread_starts[tid];
                for (std::size_t i = beg; i < end; ++i) {
                    std::uint32_t x = a[i];
                    std::uint32_t d = (x >> shift) & mask;
                    b[starts[d]++] = x;
                }
            };
            for (unsigned t = 0; t < threads; ++t) pool.emplace_back(scatter_worker, t);
            for (auto& th : pool) th.join();
        }

        a.swap(b);
    }

    // Convert back to original values: add minv
    for (std::size_t i = 0; i < n; ++i) a[i] += minv;
    return a;
}

static std::vector<std::uint32_t> dead_weight_wrapper(
    const std::vector<std::uint32_t>& in) { return dead_weight_sort(in); }

// ---------------------------------------------------------------------------
// Wrappers (function-pointer compatible)
// ---------------------------------------------------------------------------

static std::vector<std::uint32_t> std_sort_wrapper(
    const std::vector<std::uint32_t>& input)
{
    std::vector<std::uint32_t> a = input;
    std::sort(a.begin(), a.end());
    return a;
}

static std::vector<std::uint32_t> std_stable_sort_wrapper(
    const std::vector<std::uint32_t>& input)
{
    std::vector<std::uint32_t> a = input;
    std::stable_sort(a.begin(), a.end());
    return a;
}

// ---------------------------------------------------------------------------
// Flash Sort (Neubert 1998) — parallel version using LINEAR bucket spacing
// Formula: bucket = floor( (B-1) * (x - min) / (max - min) )
// Same parallel structure as gravity_physics_sort for fair comparison.
// The ONLY difference is linear vs log bucket assignment.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> flash_sort_linear(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    std::uint32_t minv = *std::min_element(arr.begin(), arr.end());
    std::uint32_t maxv = *std::max_element(arr.begin(), arr.end());
    if (minv == maxv) return arr;

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)arr.size()) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, arr.size() / 4));

    // Neubert's linear formula (parallel version)
    const double range  = static_cast<double>(maxv) - static_cast<double>(minv);
    const double scale  = static_cast<double>(B - 1) / range;

    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        return static_cast<std::size_t>((static_cast<double>(x - minv)) * scale);
    };

    unsigned threads = threads_hint == 0 ? 4u : threads_hint;
    threads = std::max(1u, std::min<unsigned>(threads,
                            static_cast<unsigned>(arr.size())));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));

    {
        std::vector<std::thread> pool;
        pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t begin = (arr.size() * tid) / threads;
            std::size_t end   = (arr.size() * (tid + 1)) / threads;
            for (std::size_t i = begin; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets, threads);

    std::vector<std::uint32_t> out;
    out.reserve(arr.size());
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> flash_sort_wrapper(
    const std::vector<std::uint32_t>& input) { return flash_sort_linear(input); }

// ---------------------------------------------------------------------------
// Algorithm NEW: Adaptive Box-Cox Sort
//
// Physical story: GPS assumes particles fall through a power-law fluid (heavy
// drag, log spacing). Flash Sort assumes vacuum (linear spacing). Real data
// sits somewhere in between — the drag exponent equals the distribution shape.
//
// One extra statistic: compute mean alongside min/max (single pass, no
// sampling). The mean's position within [min,max] reveals the skew:
//
//   lambda = 2 * (mean - min) / (max - min)
//
//   lambda ~= 1  →  mean is centred  →  uniform data  →  linear buckets (Flash Sort)
//   lambda ~= 0  →  mean hugs left   →  power-law data →  log buckets  (GPS)
//   0 < lambda < 1  →  intermediate skew  →  Box-Cox buckets
//
// Bucket formula (Box-Cox with offset to avoid log(0)):
//   bucket(x) = floor( ((x-min+1)^lambda - 1) / ((max-min+1)^lambda - 1) * B )
//
// This is the first no-sampling distribution sort that adapts its spacing
// exponent from a single-pass moment estimate (mean).
// ---------------------------------------------------------------------------

static std::vector<std::uint32_t> adaptive_box_cox_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    // Single pass: min, max, and mean together — O(n), no sampling.
    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0;
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        sum += static_cast<double>(x);
    }
    if (minv == maxv) return arr;

    double mean  = sum / static_cast<double>(arr.size());
    double range = static_cast<double>(maxv) - static_cast<double>(minv);

    // Estimate lambda: normalised mean position tells us the distribution shape.
    double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));   // clamp away from singularities

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)arr.size()) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, arr.size() / 4));

    // Precompute scale factor based on chosen regime.
    // Regime thresholds chosen to avoid numerical instability near 0 and 1.
    const bool use_log    = (lam < 0.10);   // near GPS
    const bool use_linear = (lam > 0.90);   // near Flash Sort

    const double log_range    = std::log1p(range);
    const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(B - 1) / range;
    // Box-Cox denominator: (range+1)^lambda - 1
    const double denom_bc = std::pow(range + 1.0, lam) - 1.0;
    const double inv_bc   = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        double d = static_cast<double>(x - minv);
        double pos;
        if (use_log) {
            pos = std::log1p(d) * inv_log;
        } else if (use_linear) {
            pos = d * linear_scale;
        } else {
            pos = (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        }
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    unsigned threads = threads_hint == 0 ? 4u : threads_hint;
    threads = std::max(1u, std::min<unsigned>(threads,
                            static_cast<unsigned>(arr.size())));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));

    {
        std::vector<std::thread> pool;
        pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t begin = (arr.size() * tid) / threads;
            std::size_t end   = (arr.size() * (tid + 1)) / threads;
            for (std::size_t i = begin; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets, threads);

    std::vector<std::uint32_t> out;
    out.reserve(arr.size());
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> adaptive_box_cox_wrapper(
    const std::vector<std::uint32_t>& input) { return adaptive_box_cox_sort(input); }

// ---------------------------------------------------------------------------
// Idea 1: Variance-Adaptive Tank Depth
// Physics: tank depth needed = f(spread of particle sizes). Low σ → shallow
// tank (few buckets). High σ → deep tank (many buckets).
// Adds sum_sq to the single pass → σ → scale B by sigma_norm.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> variance_adaptive_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0, sum_sq = 0.0;
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        double xd = static_cast<double>(x);
        sum    += xd;
        sum_sq += xd * xd;
    }
    if (minv == maxv) return arr;

    const double nd    = static_cast<double>(arr.size());
    double mean        = sum / nd;
    double range       = static_cast<double>(maxv) - static_cast<double>(minv);
    double variance    = std::max(0.0, sum_sq / nd - mean * mean);
    double sigma_norm  = (range > 0.0) ? std::min(1.0, std::sqrt(variance) / (range / 2.0)) : 1.0;

    double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));

    // Variance-scaled tank depth: low spread → fewer buckets, high spread → more
    double b_scale = std::max(0.4, std::min(2.0, sigma_norm * 1.5 + 0.3));
    std::size_t B = static_cast<std::size_t>(std::sqrt((double)arr.size()) * 2.0 * b_scale);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, arr.size() / 4));

    const bool use_log    = (lam < 0.10);
    const bool use_linear = (lam > 0.90);
    const double log_range    = std::log1p(range);
    const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(B - 1) / range;
    const double denom_bc     = std::pow(range + 1.0, lam) - 1.0;
    const double inv_bc       = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        double d = static_cast<double>(x - minv);
        double pos = use_log ? std::log1p(d) * inv_log
                   : use_linear ? d * linear_scale
                   : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(arr.size())));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t begin = (arr.size() * tid) / threads;
            std::size_t end   = (arr.size() * (tid + 1)) / threads;
            for (std::size_t i = begin; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }
    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }
    parallel_sort_buckets(buckets, threads);
    std::vector<std::uint32_t> out; out.reserve(arr.size());
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}
static std::vector<std::uint32_t> variance_adaptive_wrapper(
    const std::vector<std::uint32_t>& in) { return variance_adaptive_sort(in); }

// ---------------------------------------------------------------------------
// Idea 2: Pressure Redistribution (Two-Pass Equalization)
// Physics: overpressured tank regions push particles sideways into adjacent
// zones. After the first GPS scatter, if any bucket holds >2*(n/B) elements,
// recompute value boundaries from cumulative occupancy and re-scatter.
// Zero overhead on well-behaved data.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> pressure_equalized_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    std::uint32_t minv = arr[0], maxv = arr[0];
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
    }
    if (minv == maxv) return arr;

    const std::size_t n = arr.size();
    std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

    const double log_range = std::log1p(static_cast<double>(maxv) - static_cast<double>(minv));
    const double inv_log   = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;

    auto gps_bucket = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        std::size_t b = static_cast<std::size_t>(
            std::log1p(static_cast<double>(x - minv)) * inv_log);
        return b < B ? b : B - 1;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    // --- Pass 1: GPS scatter ---
    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local[tid][gps_bucket(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    // --- Pressure check ---
    const std::size_t thresh = std::max<std::size_t>(2, 2 * n / B);
    bool overpressured = false;
    for (std::size_t b = 0; b < B && !overpressured; ++b)
        if (buckets[b].size() > thresh) overpressured = true;

    if (!overpressured) {
        // Well-behaved data: no second pass needed
        parallel_sort_buckets(buckets, threads);
        std::vector<std::uint32_t> out; out.reserve(n);
        for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
        return out;
    }

    // --- Pressure relief: equalized value thresholds from GPS occupancy ---
    // Cumulative GPS bucket sizes
    std::vector<std::size_t> cum(B + 1, 0);
    for (std::size_t b = 0; b < B; ++b) cum[b + 1] = cum[b] + buckets[b].size();

    // GPS inverse: bucket boundary b → lower value of that bucket
    auto gps_inv = [&](std::size_t b) -> std::uint32_t {
        if (b == 0) return minv;
        if (b >= B) return maxv;
        double v = std::expm1(static_cast<double>(b) / static_cast<double>(B) * log_range);
        std::uint64_t vi = static_cast<std::uint64_t>(v);
        std::uint64_t rng = static_cast<std::uint64_t>(maxv - minv);
        if (vi > rng) vi = rng;
        return minv + static_cast<std::uint32_t>(vi);
    };

    // Build B equalized value thresholds
    std::vector<std::uint32_t> val_thresh(B + 1);
    val_thresh[0] = minv;
    val_thresh[B] = maxv + 1;   // sentinel (exclusive upper bound)
    for (std::size_t k = 1; k < B; ++k) {
        std::size_t target = k * n / B;
        // GPS bucket whose cumulative total first exceeds target
        std::size_t gb = static_cast<std::size_t>(
            std::lower_bound(cum.begin(), cum.end(), target + 1) - cum.begin());
        if (gb > B) gb = B;
        val_thresh[k] = gps_inv(gb);
        // Enforce strict monotonicity
        if (val_thresh[k] <= val_thresh[k - 1])
            val_thresh[k] = val_thresh[k - 1] + 1;
        if (val_thresh[k] > maxv) val_thresh[k] = maxv;
    }

    // New bucket function via binary search on equalized thresholds
    auto eq_bucket = [&](std::uint32_t x) -> std::size_t {
        auto it = std::upper_bound(val_thresh.begin(), val_thresh.end() - 1, x);
        std::size_t b = static_cast<std::size_t>(it - val_thresh.begin());
        if (b > 0) --b;
        return b < B ? b : B - 1;
    };

    // Flatten back to a single array
    std::vector<std::uint32_t> flat; flat.reserve(n);
    for (auto& bkt : buckets) flat.insert(flat.end(), bkt.begin(), bkt.end());

    // --- Pass 2: re-scatter with equalized boundaries ---
    std::vector<std::vector<std::vector<std::uint32_t>>> local2(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker2 = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local2[tid][eq_bucket(flat[i])].push_back(flat[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker2, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets2(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local2[t][b].size();
        buckets2[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local2[t][b];
            buckets2[b].insert(buckets2[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets2, threads);
    std::vector<std::uint32_t> out; out.reserve(n);
    for (auto& bkt : buckets2) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> pressure_equalized_wrapper(
    const std::vector<std::uint32_t>& in) { return pressure_equalized_sort(in); }

// ---------------------------------------------------------------------------
// Idea 3: Layered / Stratified Fluid Sort
// Physics: a stratified ocean has density layers — particles falling through
// encounter different drag at different depths. After GPS scatter, each bucket
// that is still large and internally spread is sorted recursively with a fresh
// local λ estimated from the sub-bucket's own mean. Bounded to 2 recursion
// levels to keep total work O(n).
// ---------------------------------------------------------------------------

static void layered_sort_bucket(std::vector<std::uint32_t>& bkt, int depth)
{
    const std::size_t sz = bkt.size();
    if (sz <= 128 || depth >= 2) {
        if (sz <= 32) insertion_sort(bkt);
        else          std::sort(bkt.begin(), bkt.end());
        return;
    }

    // Single pass: local min, max, mean
    std::uint32_t lmin = bkt[0], lmax = bkt[0];
    double lsum = 0.0;
    for (auto x : bkt) {
        if (x < lmin) lmin = x;
        if (x > lmax) lmax = x;
        lsum += static_cast<double>(x);
    }
    if (lmin == lmax) return;   // all equal — already "sorted"

    const double lrange = static_cast<double>(lmax) - static_cast<double>(lmin);
    const double lmean  = lsum / static_cast<double>(sz);

    // Estimate local λ from mean position (same physics as Adaptive Box-Cox)
    double lam = 2.0 * (lmean - static_cast<double>(lmin)) / lrange;
    lam = std::max(0.02, std::min(1.98, lam));

    const std::size_t sub_B = std::max<std::size_t>(8,
        static_cast<std::size_t>(std::sqrt((double)sz)));

    const bool use_log    = (lam < 0.10);
    const bool use_linear = (lam > 0.90);
    const double log_range    = std::log1p(lrange);
    const double inv_log      = (log_range > 0.0)
                                 ? (static_cast<double>(sub_B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(sub_B - 1) / lrange;
    const double denom_bc     = std::pow(lrange + 1.0, lam) - 1.0;
    const double inv_bc       = (denom_bc > 1e-15)
                                 ? (static_cast<double>(sub_B) / denom_bc) : 1.0;

    auto sub_bucket = [&](std::uint32_t x) -> std::size_t {
        if (x <= lmin) return 0;
        if (x >= lmax) return sub_B - 1;
        double d = static_cast<double>(x - lmin);
        double pos = use_log    ? std::log1p(d) * inv_log
                   : use_linear ? d * linear_scale
                   : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < sub_B ? b : sub_B - 1;
    };

    std::vector<std::vector<std::uint32_t>> subs(sub_B);
    for (auto x : bkt) subs[sub_bucket(x)].push_back(x);

    // Recurse
    for (auto& s : subs) layered_sort_bucket(s, depth + 1);

    // Gather back
    bkt.clear();
    for (auto& s : subs) bkt.insert(bkt.end(), s.begin(), s.end());
}

static std::vector<std::uint32_t> layered_fluid_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    // GPS first pass (identical to gravity_physics_sort)
    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0;
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        sum += static_cast<double>(x);
    }
    if (minv == maxv) return arr;

    const std::size_t n = arr.size();
    const double mean  = sum / static_cast<double>(n);
    const double range = static_cast<double>(maxv) - static_cast<double>(minv);
    double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

    const bool use_log    = (lam < 0.10);
    const bool use_linear = (lam > 0.90);
    const double log_range    = std::log1p(range);
    const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(B - 1) / range;
    const double denom_bc     = std::pow(range + 1.0, lam) - 1.0;
    const double inv_bc       = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        double d = static_cast<double>(x - minv);
        double pos = use_log    ? std::log1p(d) * inv_log
                   : use_linear ? d * linear_scale
                   : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    // Layered sort: recurse into large buckets with fresh local λ
    {
        std::atomic<std::size_t> next{0};
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&]() {
            for (;;) {
                std::size_t b = next.fetch_add(1, std::memory_order_relaxed);
                if (b >= B) break;
                layered_sort_bucket(buckets[b], 0);
            }
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker);
        for (auto& th : pool) th.join();
    }

    std::vector<std::uint32_t> out; out.reserve(n);
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> layered_fluid_wrapper(
    const std::vector<std::uint32_t>& in) { return layered_fluid_sort(in); }

// ---------------------------------------------------------------------------
// Idea 4: Immiscible Fluids / Two-Population Split
// Physics: oil and water separate regardless of individual droplet sizes.
// A bimodal distribution has two populations that don't overlap — split at
// the detected gap, sort each half independently, merge by concatenation.
//
// Detection: |mean − mid| / range > 0.15  (mean is pulled away from center
// when one population dominates). Split at the mean, sort each half with
// adaptive_box_cox_sort, concatenate (no extra comparisons needed).
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> immiscible_fluids_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0;
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        sum += static_cast<double>(x);
    }
    if (minv == maxv) return arr;

    const double range = static_cast<double>(maxv) - static_cast<double>(minv);
    const double mean  = sum / static_cast<double>(arr.size());
    const double mid   = static_cast<double>(minv) + range * 0.5;

    // Detect skew: mean far from centre → bimodal or strongly skewed distribution
    // Use the split only when it meaningfully separates two populations
    const double skew_ratio = std::abs(mean - mid) / range;
    if (skew_ratio < 0.15) {
        // Well-centred distribution: no benefit from splitting, use AGS directly
        return adaptive_box_cox_sort(arr, threads_hint);
    }

    // Split at the mean (natural boundary between the two populations)
    const std::uint32_t split = static_cast<std::uint32_t>(mean);

    std::vector<std::uint32_t> lower, upper;
    lower.reserve(arr.size() / 2);
    upper.reserve(arr.size() / 2);
    for (auto x : arr) {
        if (x <= split) lower.push_back(x);
        else            upper.push_back(x);
    }

    // Sort each population independently with its own λ estimate
    auto lo_sorted = (lower.size() <= 1) ? lower : adaptive_box_cox_sort(lower, threads_hint);
    auto hi_sorted = (upper.size() <= 1) ? upper : adaptive_box_cox_sort(upper, threads_hint);

    // Merge by concatenation (gap guarantees the two sorted halves are in order)
    std::vector<std::uint32_t> out;
    out.reserve(arr.size());
    out.insert(out.end(), lo_sorted.begin(), lo_sorted.end());
    out.insert(out.end(), hi_sorted.begin(), hi_sorted.end());
    return out;
}

static std::vector<std::uint32_t> immiscible_fluids_wrapper(
    const std::vector<std::uint32_t>& in) { return immiscible_fluids_sort(in); }

// ---------------------------------------------------------------------------
// Idea 5: Sedimentation Rate Prediction (Single-Pass λ Estimation)
// Physics: a scientist watches the first few particles settle and predicts
// where the rest will land without waiting for all of them.
//
// Implementation: scan first k=5% of elements for min/max/mean → estimate λ
// → scatter ALL n elements in a single subsequent pass. Saves the second full
// statistics scan. Outliers (values outside the sample's min/max) are clamped
// to boundary buckets. Works best for distributions where the extremes are
// well-represented in the sample (uniform, clustered, nearly sorted).
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> sedimentation_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    const std::size_t n = arr.size();
    // Sample: first k% of the array
    const std::size_t k = std::max<std::size_t>(64, n / 20);   // ~5%, min 64

    std::uint32_t sminv = arr[0], smaxv = arr[0];
    double ssum = 0.0;
    for (std::size_t i = 0; i < k; ++i) {
        std::uint32_t x = arr[i];
        if (x < sminv) sminv = x;
        if (x > smaxv) smaxv = x;
        ssum += static_cast<double>(x);
    }
    if (sminv == smaxv) {
        // Sample is degenerate — fallback: full scan
        for (std::size_t i = k; i < n; ++i) {
            if (arr[i] < sminv) sminv = arr[i];
            if (arr[i] > smaxv) smaxv = arr[i];
        }
        if (sminv == smaxv) return arr;
        ssum = 0.0;
        for (auto x : arr) ssum += static_cast<double>(x);
    }

    const double range = static_cast<double>(smaxv) - static_cast<double>(sminv);
    const double smean = ssum / static_cast<double>(k);

    double lam = 2.0 * (smean - static_cast<double>(sminv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

    const bool use_log    = (lam < 0.10);
    const bool use_linear = (lam > 0.90);
    const double log_range    = std::log1p(range);
    const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(B - 1) / range;
    const double denom_bc     = std::pow(range + 1.0, lam) - 1.0;
    const double inv_bc       = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

    // Bucket function: outliers clamped to boundary buckets
    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= sminv) return 0;
        if (x >= smaxv) return B - 1;
        double d = static_cast<double>(x - sminv);
        double pos = use_log    ? std::log1p(d) * inv_log
                   : use_linear ? d * linear_scale
                   : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    // Single scatter pass over ALL n elements using the sample-estimated parameters
    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets, threads);

    std::vector<std::uint32_t> out; out.reserve(n);
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> sedimentation_wrapper(
    const std::vector<std::uint32_t>& in) { return sedimentation_sort(in); }

// ---------------------------------------------------------------------------
// Idea 6: Reynolds Number Switching (Small-n vs Large-n Regime)
// Physics: at low Reynolds number, drag is linear (Stokes flow, simple);
// at high Reynolds number, turbulent drag justifies a complex model.
// The "Reynolds number" of sorting = the n where GPS overhead becomes worth it.
//
// Below the threshold, thread setup + bucket allocation + log() calls cost
// more than std::sort on the flat array. Above it, GPS dominates.
// This version also applies the AGS formula (best distribution-sort) for
// the large-n regime — giving a two-regime adaptive algorithm.
// Threshold calibrated empirically: std::sort beats GPS below ~8192 elements.
// ---------------------------------------------------------------------------
static constexpr std::size_t REYNOLDS_THRESHOLD = 8192;

static std::vector<std::uint32_t> reynolds_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    const std::size_t n = arr.size();

    // Laminar regime (low Re): std::sort wins — no bucket overhead
    if (n < REYNOLDS_THRESHOLD) {
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());
        return out;
    }

    // Turbulent regime (high Re): AGS gives best distribution-adaptive result
    return adaptive_box_cox_sort(arr, threads_hint);
}

static std::vector<std::uint32_t> reynolds_wrapper(
    const std::vector<std::uint32_t>& in) { return reynolds_sort(in); }

// ---------------------------------------------------------------------------
// Idea 7: Magnetic Levitation Stable Points (Histogram-Equalized Boundaries)
// Physics: in magnetic levitation, stable positions are where gravitational
// force exactly balances magnetic repulsion — tuned by the field, not the
// particle. Here, bucket boundaries are "levitated" to quantile positions so
// each bucket receives exactly n/B elements.
//
// Two-pass: (1) GPS scatter into 256 coarse histogram bins; (2) compute
// cumulative quantile boundaries and scatter into B equalized buckets.
// The GPS first pass is better-calibrated for skewed data than linear bins,
// so fewer coarse bins are needed than a naive histogram sort.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> magnetic_levitation_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    std::uint32_t minv = arr[0], maxv = arr[0];
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
    }
    if (minv == maxv) return arr;

    const std::size_t n = arr.size();
    const std::size_t COARSE = 256;   // histogram bins (GPS-spaced)

    const double log_range   = std::log1p(static_cast<double>(maxv) - static_cast<double>(minv));
    const double inv_log_c   = (log_range > 0.0) ? (static_cast<double>(COARSE) / log_range) : 1.0;

    auto coarse_bucket = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return COARSE - 1;
        std::size_t b = static_cast<std::size_t>(
            std::log1p(static_cast<double>(x - minv)) * inv_log_c);
        return b < COARSE ? b : COARSE - 1;
    };

    // Pass 1: count occupancy of 256 coarse GPS bins
    std::vector<std::size_t> hist(COARSE, 0);
    for (auto x : arr) ++hist[coarse_bucket(x)];

    // Cumulative counts → find value threshold at each n/B quantile
    std::vector<std::size_t> cum(COARSE + 1, 0);
    for (std::size_t i = 0; i < COARSE; ++i) cum[i + 1] = cum[i] + hist[i];

    // Final bucket count (same as GPS)
    std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

    // GPS inverse for coarse bins: coarse bin c → lower value threshold
    auto coarse_inv = [&](std::size_t c) -> std::uint32_t {
        if (c == 0) return minv;
        if (c >= COARSE) return maxv;
        double v = std::expm1(static_cast<double>(c) / static_cast<double>(COARSE) * log_range);
        std::uint64_t vi = static_cast<std::uint64_t>(v);
        std::uint64_t rng = static_cast<std::uint64_t>(maxv - minv);
        if (vi > rng) vi = rng;
        return minv + static_cast<std::uint32_t>(vi);
    };

    // Build B equalized value thresholds from coarse histogram
    std::vector<std::uint32_t> val_thresh(B + 1);
    val_thresh[0] = minv;
    val_thresh[B] = maxv + 1;
    for (std::size_t k = 1; k < B; ++k) {
        std::size_t target = k * n / B;
        std::size_t c = static_cast<std::size_t>(
            std::lower_bound(cum.begin(), cum.end(), target + 1) - cum.begin());
        if (c > COARSE) c = COARSE;
        val_thresh[k] = coarse_inv(c);
        if (val_thresh[k] <= val_thresh[k - 1])
            val_thresh[k] = val_thresh[k - 1] + 1;
        if (val_thresh[k] > maxv) val_thresh[k] = maxv;
    }

    // Pass 2: scatter all elements into equalized buckets (binary search)
    auto eq_bucket = [&](std::uint32_t x) -> std::size_t {
        auto it = std::upper_bound(val_thresh.begin(), val_thresh.end() - 1, x);
        std::size_t b = static_cast<std::size_t>(it - val_thresh.begin());
        if (b > 0) --b;
        return b < B ? b : B - 1;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local[tid][eq_bucket(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets, threads);
    std::vector<std::uint32_t> out; out.reserve(n);
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> magnetic_levitation_wrapper(
    const std::vector<std::uint32_t>& in) { return magnetic_levitation_sort(in); }

// ---------------------------------------------------------------------------
// Idea 8: Buoyancy-Aware Reverse Sort
// Physics: particles less dense than the fluid float upward. For reverse-
// sorted input, the "lightest" (smallest) elements are at the bottom; they
// need to rise. Flipping the bucket formula sorts them naturally without
// an explicit reversal step.
//
// Detection: if arr[0] > arr[n-1] and mean > mid-range
//   → near-reverse-sorted → invert bucket assignment.
// Expected win on skewed reverse-sorted data where GPS normally piles
// all heavy particles into one bucket.
// Note: for uniformly distributed reverse-sorted data the flip is neutral
// (same bucket balance either way). The real win is for skewed distributions.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> buoyancy_aware_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    if (arr.size() <= 1) return arr;

    const std::size_t n = arr.size();
    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0;
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        sum += static_cast<double>(x);
    }
    if (minv == maxv) return arr;

    const double range = static_cast<double>(maxv) - static_cast<double>(minv);
    const double mean  = sum / static_cast<double>(n);
    const double mid   = static_cast<double>(minv) + range * 0.5;

    // Detect reverse-sorted / buoyant regime:
    // (a) First element > last element (decreasing trend)
    // (b) Mean above centre (distribution heavy toward high values, suggesting
    //     particles that would rise in GPS's log-spaced fluid)
    const bool reverse_trend  = (arr[0] > arr[n - 1]);
    const bool heavy_high     = (mean > mid);
    const bool use_flip       = (reverse_trend && heavy_high);

    // Box-Cox λ from mean
    double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

    const bool use_log    = (lam < 0.10);
    const bool use_linear = (lam > 0.90);
    const double log_range    = std::log1p(range);
    const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(B - 1) / range;
    const double denom_bc     = std::pow(range + 1.0, lam) - 1.0;
    const double inv_bc       = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

    auto bucket_of_fwd = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        double d = static_cast<double>(x - minv);
        double pos = use_log    ? std::log1p(d) * inv_log
                   : use_linear ? d * linear_scale
                   : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    // Buoyancy flip: invert bucket index so heavy (large) values float to bucket 0
    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        std::size_t b = bucket_of_fwd(x);
        return use_flip ? (B - 1 - b) : b;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    parallel_sort_buckets(buckets, threads);

    std::vector<std::uint32_t> out; out.reserve(n);
    if (use_flip) {
        // Flipped: bucket 0 has largest values, gather in reverse order
        for (std::size_t b = B; b-- > 0; )
            out.insert(out.end(), buckets[b].begin(), buckets[b].end());
    } else {
        for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    }
    return out;
}

static std::vector<std::uint32_t> buoyancy_aware_wrapper(
    const std::vector<std::uint32_t>& in) { return buoyancy_aware_sort(in); }

// ---------------------------------------------------------------------------
// COMBINED: Gravity Combined Sort
// Integrates three validated wins into one algorithm:
//   1. Reynolds Switching  — n < 8192 → std::sort (avoids GPS thread overhead)
//   2. Adaptive Box-Cox    — single-pass λ estimate, Box-Cox scatter
//   3. Layered Fluid       — recursive local-λ re-sort only for overloaded
//                            buckets (> 3×expected), zero overhead otherwise
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> gravity_combined_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    const std::size_t n = arr.size();
    if (n <= 1) return arr;

    // ── Reynolds: small n → std::sort wins outright ──────────────────────
    if (n < REYNOLDS_THRESHOLD) {
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());
        return out;
    }

    // ── AGS: single pass — min, max, mean ────────────────────────────────
    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0;
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        sum += static_cast<double>(x);
    }
    if (minv == maxv) return arr;

    const double mean  = sum / static_cast<double>(n);
    const double range = static_cast<double>(maxv) - static_cast<double>(minv);
    double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

    const bool use_log    = (lam < 0.10);
    const bool use_linear = (lam > 0.90);
    const double log_range    = std::log1p(range);
    const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(B - 1) / range;
    const double denom_bc     = std::pow(range + 1.0, lam) - 1.0;
    const double inv_bc       = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        double d = static_cast<double>(x - minv);
        double pos = use_log    ? std::log1p(d) * inv_log
                   : use_linear ? d * linear_scale
                   : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    // ── Turbulent rebalancing (O(B) scan): merge adjacent lopsided pairs
    //    before sorting. Only triggers when one bucket > 3× its neighbor.
    //    Cost on well-balanced data: O(B) = O(√n) — effectively free. ────
    for (std::size_t b = 0; b + 1 < B; ++b) {
        std::size_t sa = buckets[b].size(), sb = buckets[b + 1].size();
        if (sa == 0 || sb == 0) continue;
        if (sa > 3 * sb || sb > 3 * sa) {
            buckets[b].insert(buckets[b].end(),
                              buckets[b + 1].begin(), buckets[b + 1].end());
            buckets[b + 1].clear();
        }
    }

    // ── Layered: overloaded buckets get recursive sub-GPS; others get
    //    std::sort or insertion_sort as before. Threshold = 3×expected size
    //    so well-balanced distributions pay zero extra cost. ────────────────
    const std::size_t layered_thresh = std::max<std::size_t>(128, 3 * n / B);
    {
        std::atomic<std::size_t> next{0};
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&]() {
            for (;;) {
                std::size_t b = next.fetch_add(1, std::memory_order_relaxed);
                if (b >= B) break;
                auto& bkt = buckets[b];
                if      (bkt.size() > layered_thresh) layered_sort_bucket(bkt, 0);
                else if (bkt.size() > 32)             std::sort(bkt.begin(), bkt.end());
                else                                  insertion_sort(bkt);
            }
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker);
        for (auto& th : pool) th.join();
    }

    std::vector<std::uint32_t> out; out.reserve(n);
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> gravity_combined_wrapper(
    const std::vector<std::uint32_t>& in) { return gravity_combined_sort(in); }

// ---------------------------------------------------------------------------
// New Idea 1: Phase Transition Sort
// Physics: matter transitions between solid (sorted), liquid (nearly sorted),
// and gas (random) phases. The optimal sort strategy differs per phase.
//
// Key insight: probe pairs spaced n/64 apart — not contiguous windows.
// Gapped pairs detect GLOBAL order (arr[0] vs arr[n/64], arr[n/64] vs arr[n/32],
// etc.) so multi-run and pipe-organ data show descents that local windows miss.
//
// Solid  (0 desc): verify full array O(n) then return
// Liquid (few desc/asc): std::sort — fast for nearly/reverse sorted
// Gas    (mixed):  gravity_combined_sort — full adaptive pipeline
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> phase_transition_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    const std::size_t n = arr.size();
    if (n <= 1) return arr;

    // 64 gapped pairs: arr[k * gap] vs arr[k * gap + gap]
    // gap = n/65 so all 64 pairs fit within [0, n-1]; global array coverage
    const std::size_t samples = std::min<std::size_t>(64, n - 1);
    const std::size_t gap     = std::max<std::size_t>(1, n / (samples + 1));
    int asc = 0, desc = 0;
    for (std::size_t k = 0; k < samples; ++k) {
        std::size_t i = k * gap, j = i + gap;
        if (j >= n) break;
        if      (arr[i] < arr[j]) ++asc;
        else if (arr[i] > arr[j]) ++desc;
    }
    const int total = asc + desc;
    // If all sampled pairs are equal (ties), do a safe full verify
    if (total == 0) {
        bool sorted = true;
        for (std::size_t i = 0; i + 1 < n; ++i) {
            if (arr[i] > arr[i + 1]) { sorted = false; break; }
        }
        if (sorted) return arr;
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());
        return out;
    }

    // SOLID: no descents in any gapped pair → globally ordered
    if (desc == 0) {
        // Fast full verify; nearly free on truly sorted data
        bool sorted = true;
        for (std::size_t i = 0; i + 1 < n; ++i) {
            if (arr[i] > arr[i + 1]) { sorted = false; break; }
        }
        if (sorted) return arr;
        // A few local inversions — std::sort with timsort is fast here
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());
        return out;
    }

    // REVERSE: no ascending gapped pairs → globally descending
    if (asc == 0) {
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());   // timsort detects reverse in O(n)
        return out;
    }

    // LIQUID: strongly directional (< 15% minority direction)
    // Catches near-sorted and near-reverse, NOT many-runs (which shows ~50/50)
    if (desc * 7 < total || asc * 7 < total) {
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());
        return out;
    }

    // GAS: mixed order — choose based on effective bit-width of value range.
    // dead_weight_sort is dramatically faster when range fits in ≤24 bits
    // (count arrays of ≤4096 entries fit in L1 cache). For wider ranges
    // (effective_bits > 24), dead_weight degenerates to parallel_radix with
    // 65536-entry count arrays — use gravity_combined instead.
    {
        std::uint32_t minv = arr[0], maxv = arr[0];
        for (auto x : arr) {
            if (x < minv) minv = x;
            if (x > maxv) maxv = x;
        }
        const std::uint32_t range = maxv - minv;
        int effective_bits = (range > 0) ? (32 - __builtin_clz(range)) : 1;
        // dead_weight wins when digit_bits < 16 (L1-resident count arrays).
        // digit_bits < 16 only when effective_bits > 16 (else 1-pass 16-bit radix,
        // same count array size as standard radix — gravity_combined is better).
        if (effective_bits > 16 && effective_bits <= 24) {
            return dead_weight_sort(arr, threads_hint);
        }
    }
    return gravity_combined_sort(arr, threads_hint);
}

static std::vector<std::uint32_t> phase_transition_wrapper(
    const std::vector<std::uint32_t>& in) { return phase_transition_sort(in); }

// ---------------------------------------------------------------------------
// New Idea 2: Turbulent Diffusion Sort
// Physics: after GPS scatter, turbulent eddies mix locally imbalanced regions.
// Two adjacent buckets that are lopsided (one bucket > 3× the other) get merged
// and re-sorted as a unit. No second scatter pass; only touches bad pairs.
// Overhead: O(B) scan + O(k log k) for each merged pair — never O(n).
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> turbulent_diffusion_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    const std::size_t n = arr.size();
    if (n <= 1) return arr;

    // --- Full gravity_combined_sort scatter phase (copy from combined) ---
    if (n < REYNOLDS_THRESHOLD) {
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());
        return out;
    }

    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0;
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        sum += static_cast<double>(x);
    }
    if (minv == maxv) return arr;

    const double mean  = sum / static_cast<double>(n);
    const double range = static_cast<double>(maxv) - static_cast<double>(minv);
    double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));

    std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
    B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
    B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

    const bool use_log    = (lam < 0.10);
    const bool use_linear = (lam > 0.90);
    const double log_range    = std::log1p(range);
    const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
    const double linear_scale = static_cast<double>(B - 1) / range;
    const double denom_bc     = std::pow(range + 1.0, lam) - 1.0;
    const double inv_bc       = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

    auto bucket_of = [&](std::uint32_t x) -> std::size_t {
        if (x <= minv) return 0;
        if (x >= maxv) return B - 1;
        double d = static_cast<double>(x - minv);
        double pos = use_log    ? std::log1p(d) * inv_log
                   : use_linear ? d * linear_scale
                   : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
        std::size_t b = static_cast<std::size_t>(pos);
        return b < B ? b : B - 1;
    };

    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(B));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i)
                local[tid][bucket_of(arr[i])].push_back(arr[i]);
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> buckets(B);
    for (std::size_t b = 0; b < B; ++b) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
        buckets[b].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][b];
            buckets[b].insert(buckets[b].end(), src.begin(), src.end());
        }
    }

    // --- Turbulent rebalancing: merge adjacent lopsided pairs ---
    // Two adjacent buckets are lopsided if one is > 3× the other.
    // Merge them into bucket b (clear b+1) so the combined sort is cheaper.
    const std::size_t lopsided_ratio = 3;
    for (std::size_t b = 0; b + 1 < B; ++b) {
        std::size_t sa = buckets[b].size(), sb = buckets[b + 1].size();
        if (sa == 0 || sb == 0) continue;   // empty bucket — already trivial
        if (sa > lopsided_ratio * sb || sb > lopsided_ratio * sa) {
            // Merge b+1 into b, clear b+1
            buckets[b].insert(buckets[b].end(),
                              buckets[b + 1].begin(), buckets[b + 1].end());
            buckets[b + 1].clear();
        }
    }

    // Sort (layered threshold on overloaded buckets)
    const std::size_t layered_thresh = std::max<std::size_t>(128, 3 * n / B);
    {
        std::atomic<std::size_t> next{0};
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&]() {
            for (;;) {
                std::size_t b = next.fetch_add(1, std::memory_order_relaxed);
                if (b >= B) break;
                auto& bkt = buckets[b];
                if      (bkt.size() > layered_thresh) layered_sort_bucket(bkt, 0);
                else if (bkt.size() > 32)             std::sort(bkt.begin(), bkt.end());
                else                                  insertion_sort(bkt);
            }
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker);
        for (auto& th : pool) th.join();
    }

    std::vector<std::uint32_t> out; out.reserve(n);
    for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
    return out;
}

static std::vector<std::uint32_t> turbulent_diffusion_wrapper(
    const std::vector<std::uint32_t>& in) { return turbulent_diffusion_sort(in); }

// ---------------------------------------------------------------------------
// New Idea 3: Shock Wave Sort
// Physics: in supersonic gas, a shock front forms at a large pressure gap.
// Analog: histogram bins where hist[b]>0 and hist[b+1]==0 are gap boundaries.
// Use these gaps as bucket splitters instead of uniform Box-Cox spacing.
// Falls through to standard Box-Cox when no gaps are detected.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> shock_wave_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    const std::size_t n = arr.size();
    if (n <= 1) return arr;

    if (n < REYNOLDS_THRESHOLD) {
        std::vector<std::uint32_t> out = arr;
        std::sort(out.begin(), out.end());
        return out;
    }

    // Single-pass: collect min, max, mean AND 256-bin histogram
    std::uint32_t minv = arr[0], maxv = arr[0];
    double sum = 0.0;
    std::size_t hist[256] = {};
    for (auto x : arr) {
        if (x < minv) minv = x;
        if (x > maxv) maxv = x;
        sum += static_cast<double>(x);
        // Map value to [0,255] bin using linear scaling
        // We'll re-map properly after we know range; accumulate raw x for now
        // Defer actual histogram to second pass over arr (still O(n))
    }
    if (minv == maxv) return arr;

    const double range = static_cast<double>(maxv) - static_cast<double>(minv);
    const double mean  = sum / static_cast<double>(n);

    // Build the 256-bin histogram in one more O(n) pass
    const double hist_scale = 255.0 / range;
    for (auto x : arr) {
        std::size_t bin = static_cast<std::size_t>(
            static_cast<double>(x - minv) * hist_scale);
        if (bin > 255) bin = 255;
        ++hist[bin];
    }

    // Detect shock fronts: positions where hist[b]>0 and hist[b+1]==0
    // Collect gap bin boundaries (value just below the empty zone)
    std::vector<std::uint32_t> splitters;
    splitters.push_back(minv);  // always start here
    for (int b = 0; b < 255; ++b) {
        if (hist[b] > 0 && hist[b + 1] == 0) {
            // Find the last real value in bin b: bin b covers
            // values in [ minv + b/hist_scale, minv + (b+1)/hist_scale )
            // Use the upper edge as a splitter
            double edge = static_cast<double>(minv)
                        + static_cast<double>(b + 1) / hist_scale;
            std::uint32_t spl = static_cast<std::uint32_t>(std::ceil(edge));
            if (spl > minv && spl < maxv)
                splitters.push_back(spl);
        }
    }
    splitters.push_back(maxv + 1);  // sentinel past end

    const std::size_t S = splitters.size();  // number of split points including sentinels
    // S-1 = number of "natural" gap-separated segments
    // Only use gap-splitting when we found enough gaps to matter
    // If we only found 1 gap (uniform data), fall through to Box-Cox scatter
    const std::size_t num_segments = S - 1;

    // Decide: if no useful gaps found, use standard combined sort
    if (num_segments < 2) {
        // No gaps detected — fall through to gravity_combined_sort scatter
        double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
        lam = std::max(0.02, std::min(1.98, lam));

        std::size_t B = static_cast<std::size_t>(std::sqrt((double)n) * 2.0);
        B = std::max<std::size_t>(32, std::min<std::size_t>(4096, B));
        B = std::min<std::size_t>(B, std::max<std::size_t>(32, n / 4));

        const bool use_log    = (lam < 0.10);
        const bool use_linear = (lam > 0.90);
        const double log_range    = std::log1p(range);
        const double inv_log      = (log_range > 0.0) ? (static_cast<double>(B) / log_range) : 1.0;
        const double linear_scale = static_cast<double>(B - 1) / range;
        const double denom_bc     = std::pow(range + 1.0, lam) - 1.0;
        const double inv_bc       = (denom_bc > 1e-15) ? (static_cast<double>(B) / denom_bc) : 1.0;

        auto bucket_of = [&](std::uint32_t x) -> std::size_t {
            if (x <= minv) return 0;
            if (x >= maxv) return B - 1;
            double d = static_cast<double>(x - minv);
            double pos = use_log    ? std::log1p(d) * inv_log
                       : use_linear ? d * linear_scale
                       : (std::pow(d + 1.0, lam) - 1.0) * inv_bc;
            std::size_t b = static_cast<std::size_t>(pos);
            return b < B ? b : B - 1;
        };

        unsigned threads = std::max(1u, std::min<unsigned>(
            threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

        std::vector<std::vector<std::vector<std::uint32_t>>> local(
            threads, std::vector<std::vector<std::uint32_t>>(B));
        {
            std::vector<std::thread> pool; pool.reserve(threads);
            auto worker = [&](unsigned tid) {
                std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
                for (std::size_t i = beg; i < end; ++i)
                    local[tid][bucket_of(arr[i])].push_back(arr[i]);
            };
            for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
            for (auto& th : pool) th.join();
        }

        std::vector<std::vector<std::uint32_t>> buckets(B);
        for (std::size_t b = 0; b < B; ++b) {
            std::size_t total = 0;
            for (unsigned t = 0; t < threads; ++t) total += local[t][b].size();
            buckets[b].reserve(total);
            for (unsigned t = 0; t < threads; ++t) {
                auto& src = local[t][b];
                buckets[b].insert(buckets[b].end(), src.begin(), src.end());
            }
        }
        const std::size_t layered_thresh = std::max<std::size_t>(128, 3 * n / B);
        {
            std::atomic<std::size_t> next{0};
            std::vector<std::thread> pool; pool.reserve(threads);
            auto worker = [&]() {
                for (;;) {
                    std::size_t b = next.fetch_add(1, std::memory_order_relaxed);
                    if (b >= B) break;
                    auto& bkt = buckets[b];
                    if      (bkt.size() > layered_thresh) layered_sort_bucket(bkt, 0);
                    else if (bkt.size() > 32)             std::sort(bkt.begin(), bkt.end());
                    else                                  insertion_sort(bkt);
                }
            };
            for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker);
            for (auto& th : pool) th.join();
        }
        std::vector<std::uint32_t> out; out.reserve(n);
        for (auto& bkt : buckets) out.insert(out.end(), bkt.begin(), bkt.end());
        return out;
    }

    // --- Gap-detected path: scatter into num_segments natural segments ---
    // Each segment is bounded by [splitters[s], splitters[s+1])
    // Within each segment, sort independently with std::sort or layered.
    unsigned threads = std::max(1u, std::min<unsigned>(
        threads_hint == 0 ? 4u : threads_hint, static_cast<unsigned>(n)));

    std::vector<std::vector<std::vector<std::uint32_t>>> local(
        threads, std::vector<std::vector<std::uint32_t>>(num_segments));
    {
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&](unsigned tid) {
            std::size_t beg = (n * tid) / threads, end = (n * (tid+1)) / threads;
            for (std::size_t i = beg; i < end; ++i) {
                std::uint32_t x = arr[i];
                // Binary search into splitters to find segment
                std::size_t lo = 0, hi = num_segments - 1;
                while (lo < hi) {
                    std::size_t mid = (lo + hi) / 2;
                    if (x < splitters[mid + 1]) hi = mid;
                    else lo = mid + 1;
                }
                local[tid][lo].push_back(x);
            }
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    std::vector<std::vector<std::uint32_t>> segs(num_segments);
    for (std::size_t s = 0; s < num_segments; ++s) {
        std::size_t total = 0;
        for (unsigned t = 0; t < threads; ++t) total += local[t][s].size();
        segs[s].reserve(total);
        for (unsigned t = 0; t < threads; ++t) {
            auto& src = local[t][s];
            segs[s].insert(segs[s].end(), src.begin(), src.end());
        }
    }

    // Sort each segment (layered if large, std::sort otherwise)
    const std::size_t layered_thresh = std::max<std::size_t>(128, 3 * n / std::max(num_segments, (std::size_t)32));
    {
        std::atomic<std::size_t> next{0};
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&]() {
            for (;;) {
                std::size_t s = next.fetch_add(1, std::memory_order_relaxed);
                if (s >= num_segments) break;
                auto& seg = segs[s];
                if      (seg.size() > layered_thresh) layered_sort_bucket(seg, 0);
                else if (seg.size() > 32)             std::sort(seg.begin(), seg.end());
                else                                  insertion_sort(seg);
            }
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker);
        for (auto& th : pool) th.join();
    }

    std::vector<std::uint32_t> out; out.reserve(n);
    for (auto& seg : segs) out.insert(out.end(), seg.begin(), seg.end());
    return out;
}

static std::vector<std::uint32_t> shock_wave_wrapper(
    const std::vector<std::uint32_t>& in) { return shock_wave_sort(in); }

// ---------------------------------------------------------------------------
// New Idea 4: Osmosis Sort (GPS-Pivoted Quicksort)
// Physics: a semi-permeable membrane sits at the chemical equilibrium point —
// the distributional median. Box-Cox inverse gives the *value* at the exact
// median bucket without sampling. This is a recursive partition sort where
// every pivot is the physics-estimated median.
// Box-Cox inverse for median: m = minv + (denom_bc/2 + 1)^(1/lam) - 1
// Falls back to std::sort for small partitions (< REYNOLDS_THRESHOLD).
// ---------------------------------------------------------------------------
static void osmosis_recurse(std::uint32_t* data, std::size_t n, int depth)
{
    if (n <= 1) return;
    if (n < 1024 || depth >= 24) {
        std::sort(data, data + n);
        return;
    }

    // Single pass: min, max, mean
    std::uint32_t minv = data[0], maxv = data[0];
    double sum = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        if (data[i] < minv) minv = data[i];
        if (data[i] > maxv) maxv = data[i];
        sum += static_cast<double>(data[i]);
    }
    if (minv == maxv) return;  // all equal

    const double range = static_cast<double>(maxv) - static_cast<double>(minv);
    const double mean  = sum / static_cast<double>(n);

    // Estimate lambda from mean position
    double lam = 2.0 * (mean - static_cast<double>(minv)) / range;
    lam = std::max(0.02, std::min(1.98, lam));

    // Compute Box-Cox inverse at the median bucket (B/2)
    // Forward: pos = (d+1)^lam - 1, d = x - minv
    // Inverse: d = (pos+1)^(1/lam) - 1
    // Median pos = denom_bc / 2 where denom_bc = (range+1)^lam - 1
    std::uint32_t pivot;
    if (lam < 0.10) {
        // log-space: median at exp(log(range+1)/2) - 1
        double d = std::exp(std::log(range + 1.0) * 0.5) - 1.0;
        pivot = minv + static_cast<std::uint32_t>(d);
    } else {
        double denom_bc = std::pow(range + 1.0, lam) - 1.0;
        if (denom_bc <= 1e-15) {
            std::sort(data, data + n);
            return;
        }
        double median_pos = denom_bc * 0.5;
        double d = std::pow(median_pos + 1.0, 1.0 / lam) - 1.0;
        if (d < 0.0) d = 0.0;
        if (d > range) d = range;
        pivot = minv + static_cast<std::uint32_t>(d);
    }

    // In-place partition around pivot (Lomuto-style, stable enough)
    std::size_t lo = 0;
    for (std::size_t i = 0; i < n; ++i)
        if (data[i] <= pivot) std::swap(data[lo++], data[i]);

    // Avoid degenerate case where pivot is too extreme
    if (lo == 0 || lo == n) {
        // Pivot didn't split — fall back to std::sort to avoid infinite recursion
        std::sort(data, data + n);
        return;
    }

    osmosis_recurse(data,      lo,     depth + 1);
    osmosis_recurse(data + lo, n - lo, depth + 1);
}

static std::vector<std::uint32_t> osmosis_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned /*threads_hint*/ = std::thread::hardware_concurrency())
{
    std::vector<std::uint32_t> out = arr;
    osmosis_recurse(out.data(), out.size(), 0);
    return out;
}

static std::vector<std::uint32_t> osmosis_wrapper(
    const std::vector<std::uint32_t>& in) { return osmosis_sort(in); }

// ---------------------------------------------------------------------------
// New Idea 5: Isotope Separation Sort (Entropy-Maximising Radix)
// Physics: choose the centrifuge speed (bit position) that maximises separation.
// For each level, count 1-bits at all 32 bit positions in O(n). Choose the
// bit position whose 1-count is closest to n/2 — most balanced partition.
// Falls back to std::sort for small arrays (< 64).
// ---------------------------------------------------------------------------
static void isotope_recurse(std::uint32_t* data, std::size_t n,
                            std::uint32_t mask_avail, int depth)
{
    if (n < 64 || mask_avail == 0 || depth >= 32) {
        std::sort(data, data + n);
        return;
    }

    // Count 1-bits at each available bit position
    std::size_t best_bit = 32;
    std::size_t best_dist = n + 1;
    for (int b = 31; b >= 0; --b) {
        std::uint32_t bitmask = 1u << b;
        if (!(mask_avail & bitmask)) continue;
        std::size_t ones = 0;
        for (std::size_t i = 0; i < n; ++i)
            if (data[i] & bitmask) ++ones;
        std::size_t dist = (ones <= n / 2) ? (n / 2 - ones) : (ones - n / 2);
        if (dist < best_dist) {
            best_dist = dist;
            best_bit  = b;
        }
        // If perfectly balanced, stop looking
        if (dist == 0) break;
    }

    if (best_bit == 32 || best_dist == n / 2) {
        // All remaining bits are uniform → no useful split
        std::sort(data, data + n);
        return;
    }

    std::uint32_t split_mask = 1u << best_bit;
    // Partition: 0-bit elements first, 1-bit elements second
    std::size_t lo = 0;
    for (std::size_t i = 0; i < n; ++i)
        if (!(data[i] & split_mask)) std::swap(data[lo++], data[i]);

    std::uint32_t next_mask = mask_avail & ~split_mask;
    isotope_recurse(data,      lo,     next_mask, depth + 1);
    isotope_recurse(data + lo, n - lo, next_mask, depth + 1);
}

static std::vector<std::uint32_t> isotope_separation_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned /*threads_hint*/ = std::thread::hardware_concurrency())
{
    std::vector<std::uint32_t> out = arr;
    if (out.size() > 1)
        isotope_recurse(out.data(), out.size(), ~0u, 0);
    return out;
}

static std::vector<std::uint32_t> isotope_separation_wrapper(
    const std::vector<std::uint32_t>& in) { return isotope_separation_sort(in); }

// ---------------------------------------------------------------------------
// New Idea 6: Convection Cell Sort (Cache-Tiled Scatter)
// Physics: Rayleigh-Bénard cells have width proportional to fluid depth
// (L2 cache depth). Each tile fits in L2, scatter is warm throughout.
// Algorithm: divide into T=65536-element tiles (L2-sized for uint32_t at 256KB),
// sort each tile with gravity_combined_sort (warm L2 throughout), then
// k-way tournament merge of the n/T sorted tiles.
// ---------------------------------------------------------------------------
static std::vector<std::uint32_t> convection_cell_sort(
    const std::vector<std::uint32_t>& arr,
    unsigned threads_hint = std::thread::hardware_concurrency())
{
    const std::size_t n = arr.size();
    if (n <= 1) return arr;

    // Tile size: L2 ~ 256KB for uint32_t → 65536 elements per tile
    // At smaller n (< 2 tiles), just use gravity_combined
    const std::size_t TILE = 65536;

    unsigned threads = threads_hint == 0 ? 4u : threads_hint;
    threads = std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(n)));

    if (n <= TILE) {
        // Fits in one tile — no benefit from tiling
        return gravity_combined_sort(arr, threads_hint);
    }

    // Split into tiles and sort each independently (parallel)
    std::size_t num_tiles = (n + TILE - 1) / TILE;
    std::vector<std::vector<std::uint32_t>> tiles(num_tiles);

    {
        std::atomic<std::size_t> next{0};
        std::vector<std::thread> pool; pool.reserve(threads);
        auto worker = [&]() {
            for (;;) {
                std::size_t t = next.fetch_add(1, std::memory_order_relaxed);
                if (t >= num_tiles) break;
                std::size_t beg = t * TILE;
                std::size_t end = std::min(beg + TILE, n);
                std::vector<std::uint32_t> tile(arr.begin() + beg, arr.begin() + end);
                // Sort the tile — it fits in L2, so scatter stays warm
                tiles[t] = gravity_combined_sort(tile, 1u);  // 1 thread per tile (nested)
            }
        };
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker);
        for (auto& th : pool) th.join();
    }

    // k-way tournament merge using a min-heap
    // Each entry: (value, tile_index, element_index)
    struct Entry {
        std::uint32_t val;
        std::size_t tile;
        std::size_t idx;
        bool operator>(const Entry& o) const { return val > o.val; }
    };

    std::vector<std::uint32_t> out; out.reserve(n);
    std::vector<Entry> heap;
    heap.reserve(num_tiles);

    for (std::size_t t = 0; t < num_tiles; ++t) {
        if (!tiles[t].empty())
            heap.push_back({tiles[t][0], t, 0});
    }
    std::make_heap(heap.begin(), heap.end(), std::greater<Entry>{});

    while (!heap.empty()) {
        std::pop_heap(heap.begin(), heap.end(), std::greater<Entry>{});
        auto e = heap.back(); heap.pop_back();
        out.push_back(e.val);
        std::size_t next_idx = e.idx + 1;
        if (next_idx < tiles[e.tile].size()) {
            heap.push_back({tiles[e.tile][next_idx], e.tile, next_idx});
            std::push_heap(heap.begin(), heap.end(), std::greater<Entry>{});
        }
    }
    return out;
}

static std::vector<std::uint32_t> convection_cell_wrapper(
    const std::vector<std::uint32_t>& in) { return convection_cell_sort(in); }

static std::vector<std::uint32_t> gravity_bucket_sort_wrapper(
    const std::vector<std::uint32_t>& input) { return gravity_bucket_sort_parallel(input); }

static std::vector<std::uint32_t> gravity_physics_sort_wrapper(
    const std::vector<std::uint32_t>& input) { return gravity_physics_sort(input); }

static std::vector<std::uint32_t> parallel_radix_sort_wrapper(
    const std::vector<std::uint32_t>& input) { return parallel_radix_sort_16(input); }

// ---------------------------------------------------------------------------
// Benchmark harness
// ---------------------------------------------------------------------------

template <typename Fn>
static TrialStats bench_full(Fn&& fn, const std::vector<std::uint32_t>& input,
    const std::vector<std::uint32_t>& expected, int repeats = 5, int warmup = 1)
{
    // Warmup runs (not recorded)
    for (int w = 0; w < warmup; ++w) fn(input);

    std::vector<double> times;
    times.reserve(repeats);
    TrialStats s;
    s.reps = repeats;
    for (int i = 0; i < repeats; ++i) {
        Timer t;
        auto out = fn(input);
        times.push_back(t.elapsed_ms());
        // Correctness: exact match AND permutation check on first run
        if (i == 0) {
            if (out != expected) {
                s.correct = false;
            } else {
                // Permutation check: sorted output must match sorted input
                auto sorted_out = out;
                std::sort(sorted_out.begin(), sorted_out.end());
                auto sorted_in  = input;
                std::sort(sorted_in.begin(), sorted_in.end());
                s.correct = (sorted_out == sorted_in);
            }
        }
    }
    std::sort(times.begin(), times.end());
    s.median_ms = times[times.size() / 2];
    s.min_ms    = times.front();
    s.max_ms    = times.back();
    double sum  = 0;
    for (double x : times) sum += x;
    s.mean_ms   = sum / repeats;
    double var  = 0;
    for (double x : times) var += (x - s.mean_ms) * (x - s.mean_ms);
    s.stddev_ms = (repeats > 1) ? std::sqrt(var / (repeats - 1)) : 0.0;
    return s;
}

// Backward-compat shim
template <typename Fn>
static double bench_median(Fn&& fn, const std::vector<std::uint32_t>& input,
    const std::vector<std::uint32_t>& expected, int repeats = 3)
{
    return bench_full(fn, input, expected, repeats, 0).median_ms;
}

static bool gpu_sort_available() {
#ifdef USE_CUDA
    return true;
#else
    return false;
#endif
}

static std::vector<std::uint32_t> gpu_sort_placeholder(
    const std::vector<std::uint32_t>& input)
{
    (void)input;
    throw std::runtime_error(
        "GPU sort not wired in. Compile with your CUDA/OpenCL adapter.");
}

// ---------------------------------------------------------------------------
// Fuzz tests (Step 9)
// ---------------------------------------------------------------------------

struct FuzzResult {
    std::string algo_name;
    std::string case_name;
    bool        passed;
    std::string error_msg;
};

static std::vector<FuzzResult> run_fuzz_tests(
    const std::vector<std::pair<std::string,
        std::vector<std::uint32_t>(*)(const std::vector<std::uint32_t>&)>>& algs)
{
    std::vector<std::pair<std::string, std::vector<std::uint32_t>>> cases;

    // Edge cases
    cases.push_back({"empty",      {}});
    cases.push_back({"single",     {42u}});
    cases.push_back({"two_asc",    {1u, 2u}});
    cases.push_back({"two_desc",   {2u, 1u}});
    cases.push_back({"two_equal",  {5u, 5u}});
    cases.push_back({"all_equal_10", std::vector<std::uint32_t>(10, 7u)});
    cases.push_back({"two_distinct", {0u, std::numeric_limits<std::uint32_t>::max()}});
    cases.push_back({"extreme_vals", {0u, 1u, std::numeric_limits<std::uint32_t>::max()-1, std::numeric_limits<std::uint32_t>::max()}});

    // Sawtooth
    {
        std::vector<std::uint32_t> v;
        for (std::uint32_t i = 0; i < 30; ++i) v.push_back(i % 5);
        cases.push_back({"sawtooth_30", v});
    }

    // Already sorted / reverse
    {
        std::vector<std::uint32_t> v(20); std::iota(v.begin(), v.end(), 0u);
        cases.push_back({"sorted_20", v});
    }
    {
        std::vector<std::uint32_t> v(20); std::iota(v.begin(), v.end(), 0u); std::reverse(v.begin(), v.end());
        cases.push_back({"reverse_20", v});
    }

    // Random n=3..50 x 20 seeds
    std::mt19937 rng(42u);
    for (int seed = 0; seed < 20; ++seed) {
        std::size_t sz = 3 + rng() % 48;
        std::vector<std::uint32_t> v(sz);
        for (auto& x : v) x = rng() % 10000;
        cases.push_back({"rand_" + std::to_string(sz) + "_s" + std::to_string(seed), v});
    }

    std::vector<FuzzResult> results;
    for (const auto& [case_name, input] : cases) {
        auto expected = input;
        std::sort(expected.begin(), expected.end());
        auto expected_multiset = expected; // already sorted = sorted input

        for (const auto& [alg_name, fn] : algs) {
            FuzzResult fr;
            fr.algo_name  = alg_name;
            fr.case_name  = case_name;
            try {
                auto out = fn(input);
                auto out_sorted = out;
                std::sort(out_sorted.begin(), out_sorted.end());
                if (out == expected && out_sorted == expected_multiset) {
                    fr.passed = true;
                } else if (out_sorted != expected_multiset) {
                    fr.passed = false;
                    fr.error_msg = "permutation mismatch (lost/duped elements)";
                } else {
                    fr.passed = false;
                    fr.error_msg = "not sorted correctly";
                }
            } catch (const std::exception& e) {
                fr.passed    = false;
                fr.error_msg = std::string("exception: ") + e.what();
            }
            results.push_back(fr);
        }
    }
    return results;
}

// ---------------------------------------------------------------------------
// JSON helpers (Step 2)
// ---------------------------------------------------------------------------

// Minimal JSON string escaper
static std::string json_str(const std::string& s) {
    std::string out = "\"";
    for (char c : s) {
        if (c == '"') out += "\\\"";
        else if (c == '\\') out += "\\\\";
        else if (c == '\n') out += "\\n";
        else out += c;
    }
    out += "\"";
    return out;
}

static std::string build_info_json() {
    std::ostringstream o;
    o << "{";
#if defined(__clang__)
    o << "\"compiler\":" << json_str(std::string("clang ") + __clang_version__) << ",";
#elif defined(__GNUC__)
    o << "\"compiler\":" << json_str(std::string("gcc ") + std::to_string(__GNUC__) + "." + std::to_string(__GNUC_MINOR__)) << ",";
#else
    o << "\"compiler\":\"unknown\",";
#endif
    o << "\"opt_level\":\"O2\",";
    o << "\"cpp_std\":" << __cplusplus << ",";
#if defined(__APPLE__)
    o << "\"os\":\"macOS\",";
#elif defined(__linux__)
    o << "\"os\":\"Linux\",";
#else
    o << "\"os\":\"unknown\",";
#endif
    o << "\"hw_threads\":" << std::thread::hardware_concurrency() << ",";
    o << "\"compile_time\":" << json_str(std::string(__DATE__) + " " + __TIME__);
    o << "}";
    return o.str();
}

struct JsonAlgResult {
    std::string name;
    TrialStats  stats;
    double      speedup_vs_baseline;
};

struct JsonDataset {
    std::string name;
    std::size_t n;
    std::vector<JsonAlgResult> algs;
};

static void output_json(
    std::size_t n, int reps_used, int warmup_used,
    unsigned long long seed,
    const std::vector<JsonDataset>& datasets,
    const std::vector<FuzzResult>& fuzz_results = {})
{
    // Get timestamp
    auto now = std::chrono::system_clock::now();
    std::time_t t = std::chrono::system_clock::to_time_t(now);
    char tsbuf[32];
    std::strftime(tsbuf, sizeof(tsbuf), "%Y-%m-%dT%H:%M:%S", std::gmtime(&t));

    std::cout << "{\n";
    std::cout << "  \"timestamp\":" << json_str(tsbuf) << ",\n";
    std::cout << "  \"seed\":" << seed << ",\n";
    std::cout << "  \"n\":" << n << ",\n";
    std::cout << "  \"element_type\":\"uint32\",\n";
    std::cout << "  \"reps\":" << reps_used << ",\n";
    std::cout << "  \"warmup\":" << warmup_used << ",\n";
    std::cout << "  \"build\":" << build_info_json() << ",\n";
    std::cout << "  \"datasets\":[\n";
    for (std::size_t di = 0; di < datasets.size(); ++di) {
        const auto& ds = datasets[di];
        std::cout << "    {\"name\":" << json_str(ds.name)
                  << ",\"n\":" << ds.n << ",\"algorithms\":[\n";
        for (std::size_t ai = 0; ai < ds.algs.size(); ++ai) {
            const auto& a = ds.algs[ai];
            const auto& s = a.stats;
            std::cout << "      {\"name\":" << json_str(a.name)
                      << ",\"mean_ms\":"    << std::fixed << std::setprecision(4) << s.mean_ms
                      << ",\"median_ms\":"  << s.median_ms
                      << ",\"min_ms\":"     << s.min_ms
                      << ",\"max_ms\":"     << s.max_ms
                      << ",\"stddev_ms\":"  << s.stddev_ms
                      << ",\"correct\":"    << (s.correct ? "true" : "false")
                      << ",\"speedup_vs_baseline\":" << std::setprecision(4) << a.speedup_vs_baseline
                      << "}";
            if (ai + 1 < ds.algs.size()) std::cout << ",";
            std::cout << "\n";
        }
        std::cout << "    ]}";
        if (di + 1 < datasets.size()) std::cout << ",";
        std::cout << "\n";
    }
    std::cout << "  ],\n";
    if (fuzz_results.empty()) {
        std::cout << "  \"fuzz_results\":null\n";
    } else {
        std::cout << "  \"fuzz_results\":[\n";
        for (std::size_t fi = 0; fi < fuzz_results.size(); ++fi) {
            const auto& fr = fuzz_results[fi];
            std::cout << "    {\"algo\":" << json_str(fr.algo_name)
                      << ",\"case\":" << json_str(fr.case_name)
                      << ",\"passed\":" << (fr.passed ? "true" : "false");
            if (!fr.error_msg.empty())
                std::cout << ",\"error\":" << json_str(fr.error_msg);
            std::cout << "}";
            if (fi + 1 < fuzz_results.size()) std::cout << ",";
            std::cout << "\n";
        }
        std::cout << "  ]\n";
    }
    std::cout << "}\n";
}

// ---------------------------------------------------------------------------
// main
// ---------------------------------------------------------------------------

int main(int argc, char** argv) {
    std::size_t n = 1000000;  // default 1M
    if (argc >= 2) {
        std::string a1(argv[1]);
        if (a1 != "--json" && a1 != "--fuzz" && a1 != "--reps" &&
            a1 != "--warmup" && a1 != "--type" && a1.substr(0,2) != "--")
            n = static_cast<std::size_t>(std::stoull(a1));
    }

    // Parse optional flags
    int  reps    = 5;
    int  warmup  = 1;
    bool json_mode  = false;
    bool fuzz_mode  = false;
    std::string elem_type = "uint32";
    std::string dataset_dir = "datasets/";

    for (int i = 1; i < argc; ++i) {
        std::string a(argv[i]);
        if (a == "--json")   { json_mode = true; }
        else if (a == "--fuzz")   { fuzz_mode = true; }
        else if (a == "--reps"  && i+1 < argc) { reps    = std::stoi(argv[++i]); }
        else if (a == "--warmup" && i+1 < argc) { warmup  = std::stoi(argv[++i]); }
        else if (a == "--type"  && i+1 < argc) { elem_type = argv[++i]; }
        else if (a == "--dir"   && i+1 < argc) { dataset_dir = argv[++i]; }
    }

    // float type support: if --type float, run float variant then exit
    if (elem_type == "float") {
        // Generate float datasets
        auto gen_rand_f = [](std::size_t n) {
            std::mt19937_64 rng(123456789ULL + n);
            std::uniform_real_distribution<float> d(0.0f, 1.0f);
            std::vector<float> v(n);
            for (auto& x : v) x = d(rng);
            return v;
        };
        auto gen_power_f = [](std::size_t n) {
            std::mt19937_64 rng(1111111111ULL + n);
            std::uniform_real_distribution<double> d(0.0, 1.0);
            std::vector<float> v;
            v.reserve(n);
            for (std::size_t i = 0; i < n; ++i) {
                double u = std::max(d(rng), 1e-9);
                v.push_back(static_cast<float>(std::min(std::pow(u, -0.7) * 0.05, 1.0)));
            }
            return v;
        };

        struct FDataset { std::string name; std::vector<float> data; };
        std::vector<FDataset> fdatasets = {
            {"Float: Uniform random", gen_rand_f(n)},
            {"Float: Power-law",      gen_power_f(n)},
        };

        if (!json_mode) std::cout << "\n[Float element type]\n";
        for (const auto& fds : fdatasets) {
            auto expected_f = fds.data;
            std::sort(expected_f.begin(), expected_f.end());

            if (!json_mode) {
                std::cout << "Dataset: " << fds.name << "  (n=" << fds.data.size() << ")\n";
                std::cout << std::left  << std::setw(36) << "algorithm"
                          << std::right << std::setw(12) << "median ms" << "\n";
                std::cout << std::string(48, '-') << "\n";
            }

            // std::sort
            auto f_std = [](const std::vector<float>& a) { auto b=a; std::sort(b.begin(),b.end()); return b; };
            // gravity physics sort (float)
            auto f_grav = [](const std::vector<float>& arr) -> std::vector<float> {
                if (arr.size() <= 1) return arr;
                float mn = *std::min_element(arr.begin(), arr.end());
                float mx = *std::max_element(arr.begin(), arr.end());
                if (mn == mx) return arr;
                std::size_t B = std::max<std::size_t>(32, std::min<std::size_t>(4096,
                    static_cast<std::size_t>(std::sqrt((double)arr.size()) * 2.0)));
                double log_range = std::log1p(static_cast<double>(mx - mn));
                double scale = (log_range > 0) ? (static_cast<double>(B) / log_range) : 1.0;
                auto bucket_of = [&](float x) -> std::size_t {
                    if (x <= mn) return 0;
                    std::size_t b = static_cast<std::size_t>(std::log1p(static_cast<double>(x-mn))*scale);
                    return b < B ? b : B-1;
                };
                std::vector<std::vector<float>> buckets(B);
                for (float x : arr) buckets[bucket_of(x)].push_back(x);
                std::vector<float> out; out.reserve(arr.size());
                for (auto& bkt : buckets) { std::sort(bkt.begin(), bkt.end()); out.insert(out.end(), bkt.begin(), bkt.end()); }
                return out;
            };

            for (auto& [name, fn] : std::vector<std::pair<std::string, std::function<std::vector<float>(const std::vector<float>&)>>>{
                    {"std::sort (float)", f_std},
                    {"gravity physics sort (float)", f_grav},
                }) {
                // warmup
                fn(fds.data);
                std::vector<double> times;
                for (int i = 0; i < reps; ++i) {
                    Timer t; auto out = fn(fds.data);
                    times.push_back(t.elapsed_ms());
                    (void)out;
                }
                std::sort(times.begin(), times.end());
                if (!json_mode) {
                    std::cout << std::left  << std::setw(36) << name
                              << std::right << std::setw(12) << std::fixed
                              << std::setprecision(3) << times[times.size()/2] << "\n";
                }
            }
            if (!json_mode) std::cout << "\n";
        }
        return 0;
    }

    if (!json_mode) std::cout << "Gravity sort benchmark  (n = " << n << ")\n\n";

    // Build fuzz algs list (same as benchmark algs, excluding GPU)
    using SortFnPtr = std::vector<std::uint32_t>(*)(const std::vector<std::uint32_t>&);
    std::vector<std::pair<std::string, SortFnPtr>> fuzz_algs = {
        {"std::sort",                    &std_sort_wrapper},
        {"std::stable_sort",             &std_stable_sort_wrapper},
        {"flash sort (linear)",          &flash_sort_wrapper},
        {"parallel radix sort",          &parallel_radix_sort_wrapper},
        {"gravity bucket sort (sampled)",&gravity_bucket_sort_wrapper},
        {"gravity physics sort (log)",   &gravity_physics_sort_wrapper},
        {"adaptive box-cox sort",        &adaptive_box_cox_wrapper},
        {"variance-adaptive sort",       &variance_adaptive_wrapper},
        {"pressure-equalized sort",      &pressure_equalized_wrapper},
        {"layered fluid sort",           &layered_fluid_wrapper},
        {"immiscible fluids sort",       &immiscible_fluids_wrapper},
        {"sedimentation sort",           &sedimentation_wrapper},
        {"reynolds sort",                &reynolds_wrapper},
        {"magnetic levitation sort",     &magnetic_levitation_wrapper},
        {"buoyancy-aware sort",          &buoyancy_aware_wrapper},
        {"gravity combined sort",        &gravity_combined_wrapper},
        {"phase transition sort",        &phase_transition_wrapper},
        {"turbulent diffusion sort",     &turbulent_diffusion_wrapper},
        {"shock wave sort",              &shock_wave_wrapper},
        {"osmosis sort",                 &osmosis_wrapper},
        {"isotope separation sort",      &isotope_separation_wrapper},
        {"convection cell sort",         &convection_cell_wrapper},
        {"dead weight sort",             &dead_weight_wrapper},
    };

    if (fuzz_mode) {
        auto fuzz_results = run_fuzz_tests(fuzz_algs);
        if (json_mode) {
            output_json(n, reps, warmup, 123456789ULL, {}, fuzz_results);
        } else {
            int total = 0, passed = 0;
            for (const auto& fr : fuzz_results) {
                ++total;
                if (fr.passed) ++passed;
                else std::cout << "FAIL  " << fr.algo_name << "  case=" << fr.case_name
                               << "  " << fr.error_msg << "\n";
            }
            std::cout << "Fuzz: " << passed << "/" << total << " passed\n";
        }
        return 0;
    }

    struct Dataset {
        std::string name;
        std::vector<std::uint32_t> data;
    };

    std::vector<Dataset> datasets;

    // --- Synthetic ---
    datasets.push_back({"Uniform random",        generate_random(n)});
    datasets.push_back({"Clustered (32 groups)", generate_clustered(n)});
    datasets.push_back({"Nearly sorted",         generate_nearly_sorted(n)});
    datasets.push_back({"Reverse sorted",        generate_reverse_sorted(n)});
    datasets.push_back({"Power-law (skewed)",    generate_power_law(n)});
    datasets.push_back({"Pipe-organ",            generate_pipe_organ(n)});
    datasets.push_back({"All equal",             generate_all_equal(n)});
    datasets.push_back({"Already sorted",        generate_already_sorted(n)});
    datasets.push_back({"Bimodal (2 clusters)",  generate_bimodal(n)});
    datasets.push_back({"Small range [0,65535]", generate_small_range(n)});
    datasets.push_back({"Many runs (128)",       generate_many_runs(n)});
    datasets.push_back({"High duplicates (16k)", generate_high_duplicates(n)});
    // Dead-weight test case: values in high sub-range [2^20, 2^21-1]
    // Only 20 bits vary; existing radix uses maxv-based pass count (2 passes of 16 bits)
    // but dead-weight sort uses 2 passes of 10 bits (1024-entry count arrays, L1-friendly)
    {
        std::mt19937 rng2(54321u + static_cast<unsigned>(n));
        std::uniform_int_distribution<std::uint32_t> dist(1u << 20, (1u << 21) - 1u);
        std::vector<std::uint32_t> v(n); for (auto& x : v) x = dist(rng2);
        datasets.push_back({"High-range [2^20,2^21)", std::move(v)});
    }

    // --- Real files (Pizza & Chili corpus, loaded as raw uint32 words) ---
    // Decompress with: gunzip datasets/pizzachili/*.gz
    // Then re-run — the loader silently skips files that are not present.
    struct FileDataset { std::string name, path; };
    for (auto& fd : std::vector<FileDataset>{
            {"PizzaChili: DNA",    dataset_dir + "pizzachili/dna.200MB"},
            {"PizzaChili: XML",    dataset_dir + "pizzachili/dblp.xml.200MB"},
            {"PizzaChili: Source", dataset_dir + "pizzachili/sources.200MB"},
    }) {
        auto v = load_file_as_uint32(fd.path, n);
        if (!v.empty()) datasets.push_back({fd.name, std::move(v)});
    }

    using SortFn = std::vector<std::uint32_t>(*)(const std::vector<std::uint32_t>&);

    struct Alg { std::string name; SortFn fn; bool enabled; };
    std::vector<Alg> algs = {
        {"std::sort (baseline)",          &std_sort_wrapper,             true},
        {"std::stable_sort",              &std_stable_sort_wrapper,      true},
        {"flash sort (linear)",           &flash_sort_wrapper,           true},
        {"parallel radix sort",           &parallel_radix_sort_wrapper,  true},
        {"gravity bucket sort (sampled)", &gravity_bucket_sort_wrapper,  true},
        {"gravity physics sort (log)",    &gravity_physics_sort_wrapper, true},
        {"adaptive box-cox sort",         &adaptive_box_cox_wrapper,     true},
        {"variance-adaptive sort",        &variance_adaptive_wrapper,    true},
        {"pressure-equalized sort",       &pressure_equalized_wrapper,   true},
        {"layered fluid sort",            &layered_fluid_wrapper,        true},
        {"immiscible fluids sort",        &immiscible_fluids_wrapper,    true},
        {"sedimentation sort",            &sedimentation_wrapper,        true},
        {"reynolds sort",                 &reynolds_wrapper,             true},
        {"magnetic levitation sort",      &magnetic_levitation_wrapper,  true},
        {"buoyancy-aware sort",           &buoyancy_aware_wrapper,       true},
        {"gravity combined sort",         &gravity_combined_wrapper,     true},
        {"phase transition sort",         &phase_transition_wrapper,     true},
        {"turbulent diffusion sort",      &turbulent_diffusion_wrapper,  true},
        {"shock wave sort",               &shock_wave_wrapper,           true},
        {"osmosis sort",                  &osmosis_wrapper,              true},
        {"isotope separation sort",       &isotope_separation_wrapper,   false},
        {"convection cell sort",          &convection_cell_wrapper,      true},
        {"dead weight sort",              &dead_weight_wrapper,          true},
        {"GPU sort",                      &gpu_sort_placeholder,         gpu_sort_available()},
    };

    std::vector<JsonDataset> json_datasets;

    for (const auto& ds : datasets) {
        if (ds.data.empty()) continue;

        // Pre-compute expected result once per dataset.
        auto expected = std_sort_wrapper(ds.data);

        // Baseline timing for speedup column.
        double baseline_ms = 0.0;

        if (!json_mode) {
            std::cout << "Dataset: " << ds.name
                      << "  (n=" << ds.data.size() << ")\n";
            std::cout << std::left  << std::setw(36) << "algorithm"
                      << std::right << std::setw(12) << "median ms"
                      << std::right << std::setw(10) << "speedup"
                      << "\n";
            std::cout << std::string(58, '-') << "\n";
        }

        JsonDataset jds;
        jds.name = ds.name;
        jds.n    = ds.data.size();

        for (const auto& alg : algs) {
            if (!alg.enabled) {
                if (!json_mode) {
                    std::cout << std::left << std::setw(36) << alg.name
                              << std::right << std::setw(12) << "n/a"
                              << std::right << std::setw(10) << "n/a" << "\n";
                }
                continue;
            }
            try {
                auto s = bench_full(alg.fn, ds.data, expected, reps, warmup);
                double speedup = (baseline_ms > 0.0) ? (baseline_ms / s.median_ms) : 1.0;
                if (baseline_ms <= 0.0) baseline_ms = s.median_ms;

                if (!json_mode) {
                    std::string speedup_str;
                    char buf[32];
                    std::snprintf(buf, sizeof(buf), "%.2fx", speedup);
                    speedup_str = buf;
                    std::cout << std::left  << std::setw(36) << alg.name
                              << std::right << std::setw(12) << std::fixed
                              << std::setprecision(3) << s.median_ms
                              << std::right << std::setw(10) << speedup_str;
                    if (!s.correct) std::cout << "  [WRONG]";
                    std::cout << "\n";
                }

                jds.algs.push_back({alg.name, s, speedup});
            } catch (const std::exception& e) {
                if (!json_mode) {
                    std::cout << std::left  << std::setw(36) << alg.name
                              << std::right << std::setw(22)
                              << (std::string("err: ") + e.what()) << "\n";
                }
                TrialStats err_s; err_s.correct = false;
                jds.algs.push_back({alg.name, err_s, 0.0});
            }
        }
        if (!json_mode) std::cout << "\n";
        json_datasets.push_back(std::move(jds));
    }

    if (json_mode) {
        output_json(n, reps, warmup, 123456789ULL, json_datasets);
    } else {
        std::cout << "Note: Pizza&Chili datasets require decompressing the .gz files first:\n"
                  << "  gunzip " << dataset_dir << "pizzachili/*.gz\n"
                  << "Note: GPU slot requires a CUDA/OpenCL implementation wired into "
                     "gpu_sort_placeholder().\n";
    }
    return 0;
}