代码

// Based on the supplied Duck.ac v6 implementation and https://duck.ac/submission/48070
// Sorting networks: Bert Dobbelaere, https://bertdobbelaere.github.io/sorting_networks.html
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <cstring>
#include <immintrin.h>
#include <utility>

#if defined(__GNUC__)
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#define DUCK_INLINE inline __attribute__((always_inline))
#define DUCK_NOINLINE __attribute__((noinline))
#define DUCK_LAMBDA_INLINE __attribute__((always_inline))
#define DUCK_EXPECT(x, value) __builtin_expect((x), (value))
#else
#define DUCK_INLINE __forceinline
#define DUCK_NOINLINE __declspec(noinline)
#define DUCK_LAMBDA_INLINE
#define DUCK_EXPECT(x, value) (x)
#endif

#ifndef DUCK_SORT_NETWORK
#define DUCK_SORT_NETWORK 12
#endif
#ifndef DUCK_SORT_BLOCK_BYTES
#define DUCK_SORT_BLOCK_BYTES 128
#endif
#ifndef DUCK_SORT_RECORDS
#define DUCK_SORT_RECORDS 42
#endif
#ifndef DUCK_SORT_FUNCTION
#define DUCK_SORT_FUNCTION sort
#endif

namespace {
using U = unsigned;
using B = unsigned char;
using H = std::uint16_t;
using Z = std::uint64_t;
using V = __m256i;
using W = __m128i;
constexpr U Network = DUCK_SORT_NETWORK;
constexpr U BlockBytes = DUCK_SORT_BLOCK_BYTES;
constexpr U Records = DUCK_SORT_RECORDS;
constexpr std::size_t WorkspaceBytes = 384000000;
constexpr std::size_t ScratchOffset = 380000000;
constexpr U GuardBytes = 25344;
static_assert(sizeof(U) == 4 && sizeof(H) == 2, "32-bit unsigned required");
static_assert(Network == 8 || Network == 10 || Network == 12, "Unsupported network");
static_assert(BlockBytes == 64 || BlockBytes == 128 || BlockBytes == 256, "Unsupported block");
static_assert(Records > 0 && Records * 3 + 1 <= BlockBytes, "Invalid packed block");
static_assert((8192 + Records - 1) / Records * BlockBytes <= GuardBytes, "Insufficient guard");
B* workspace;

template<class F, std::size_t... I>
DUCK_INLINE void each(F&& f, std::index_sequence<I...>) {
    (f(std::integral_constant<std::size_t, I>{}), ...);
}
template<std::size_t N, class F>
DUCK_INLINE void repeat(F&& f) {
    each(std::forward<F>(f), std::make_index_sequence<N>{});
}
DUCK_INLINE U load16(const H* p) {
    H x;
    std::memcpy(&x, p, sizeof(x));
    return x;
}
DUCK_INLINE U load24(const B* p) {
    U x;
    std::memcpy(&x, p, sizeof(x));
    return x & 0xffffffu;
}
DUCK_INLINE U packed_at(const B* p, U i) {
    return load24(p + std::size_t(i / Records) * BlockBytes + i % Records * 3);
}

DUCK_NOINLINE void dense(const H* src, U n, U* dst, U base) {
    alignas(64) U counts[65536] = {};
    for (U i = 0; i < n; ++i) ++counts[load16(src + i)];
    for (U v = 0; v < 65536; ++v)
        for (U c = counts[v]; c; --c) *dst++ = base | v;
}

DUCK_INLINE void batch_sort(const B* src, B* first, B* tail) {
    V x[Network], a[8], b[8];
    repeat<Network>([&](auto j) DUCK_LAMBDA_INLINE {
        x[j] = _mm256_load_si256(reinterpret_cast<const V*>(src + 256 * j));
    });
#define CMP(A, B) do { V lo = _mm256_min_epu8(x[A], x[B]); x[B] = _mm256_max_epu8(x[A], x[B]); x[A] = lo; } while (0)
    if constexpr (Network == 8) {
        CMP(0,2); CMP(1,3); CMP(4,6); CMP(5,7);
        CMP(0,4); CMP(1,5); CMP(2,6); CMP(3,7);
        CMP(0,1); CMP(2,3); CMP(4,5); CMP(6,7);
        CMP(2,4); CMP(3,5);
        CMP(1,4); CMP(3,6);
        CMP(1,2); CMP(3,4); CMP(5,6);
    } else if constexpr (Network == 10) {
        CMP(0,8); CMP(1,9); CMP(2,7); CMP(3,5); CMP(4,6);
        CMP(0,2); CMP(1,4); CMP(5,8); CMP(7,9);
        CMP(0,3); CMP(2,4); CMP(5,7); CMP(6,9);
        CMP(0,1); CMP(3,6); CMP(8,9);
        CMP(1,5); CMP(2,3); CMP(4,8); CMP(6,7);
        CMP(1,2); CMP(3,5); CMP(4,6); CMP(7,8);
        CMP(2,3); CMP(4,5); CMP(6,7);
        CMP(3,4); CMP(5,6);
    } else {
        CMP(0,8); CMP(1,7); CMP(2,6); CMP(3,11); CMP(4,10); CMP(5,9);
        CMP(0,1); CMP(2,5); CMP(3,4); CMP(6,9); CMP(7,8); CMP(10,11);
        CMP(0,2); CMP(1,6); CMP(5,10); CMP(9,11);
        CMP(0,3); CMP(1,2); CMP(4,6); CMP(5,7); CMP(8,11); CMP(9,10);
        CMP(1,4); CMP(3,5); CMP(6,8); CMP(7,10);
        CMP(1,3); CMP(2,5); CMP(6,9); CMP(8,10);
        CMP(2,3); CMP(4,5); CMP(6,7); CMP(8,9);
        CMP(4,6); CMP(5,7);
        CMP(3,4); CMP(5,6); CMP(7,8);
    }
#undef CMP
    repeat<4>([&](auto i) DUCK_LAMBDA_INLINE {
        a[2*i] = _mm256_unpacklo_epi8(x[2*i], x[2*i+1]);
        a[2*i+1] = _mm256_unpackhi_epi8(x[2*i], x[2*i+1]);
    });
    repeat<2>([&](auto i) DUCK_LAMBDA_INLINE {
        repeat<2>([&](auto j) DUCK_LAMBDA_INLINE {
            b[4*i+2*j] = _mm256_unpacklo_epi16(a[4*i+j], a[4*i+j+2]);
            b[4*i+2*j+1] = _mm256_unpackhi_epi16(a[4*i+j], a[4*i+j+2]);
        });
    });
    repeat<4>([&](auto i) DUCK_LAMBDA_INLINE {
        V lo = _mm256_unpacklo_epi32(b[i], b[i+4]);
        V hi = _mm256_unpackhi_epi32(b[i], b[i+4]);
        _mm256_store_si256(reinterpret_cast<V*>(first + 32*i), _mm256_permute2x128_si256(lo, hi, 0x20));
        _mm256_store_si256(reinterpret_cast<V*>(first + 128 + 32*i), _mm256_permute2x128_si256(lo, hi, 0x31));
    });
    if constexpr (Network == 10) {
        V lo = _mm256_unpacklo_epi8(x[8], x[9]);
        V hi = _mm256_unpackhi_epi8(x[8], x[9]);
        _mm256_store_si256(reinterpret_cast<V*>(tail), _mm256_permute2x128_si256(lo, hi, 0x20));
        _mm256_store_si256(reinterpret_cast<V*>(tail + 32), _mm256_permute2x128_si256(lo, hi, 0x31));
    } else if constexpr (Network == 12) {
        a[0] = _mm256_unpacklo_epi8(x[8], x[9]);
        a[1] = _mm256_unpackhi_epi8(x[8], x[9]);
        a[2] = _mm256_unpacklo_epi8(x[10], x[11]);
        a[3] = _mm256_unpackhi_epi8(x[10], x[11]);
        b[0] = _mm256_unpacklo_epi16(a[0], a[2]);
        b[1] = _mm256_unpackhi_epi16(a[0], a[2]);
        b[2] = _mm256_unpacklo_epi16(a[1], a[3]);
        b[3] = _mm256_unpackhi_epi16(a[1], a[3]);
        repeat<2>([&](auto i) DUCK_LAMBDA_INLINE {
            _mm256_store_si256(reinterpret_cast<V*>(tail + 32*i), _mm256_permute2x128_si256(b[2*i], b[2*i+1], 0x20));
            _mm256_store_si256(reinterpret_cast<V*>(tail + 64 + 32*i), _mm256_permute2x128_si256(b[2*i], b[2*i+1], 0x31));
        });
    }
}

DUCK_INLINE W insert_byte(W x, U value) {
    return _mm_min_epu8(x, _mm_max_epu8(_mm_slli_si128(x, 1), _mm_set1_epi8(static_cast<char>(value))));
}
DUCK_INLINE void write8(W x, U* dst, U count, U base, U* end) {
    V y = _mm256_or_si256(_mm256_cvtepu8_epi32(x), _mm256_set1_epi32(base));
    if (end - dst >= 8) {
        _mm256_storeu_si256(reinterpret_cast<V*>(dst), y);
    } else {
        V mask = _mm256_cmpgt_epi32(_mm256_set1_epi32(count), _mm256_setr_epi32(0,1,2,3,4,5,6,7));
        _mm256_maskstore_epi32(reinterpret_cast<int*>(dst), mask, y);
    }
}

void leaf(const H* src, U n, U* dst, U base, U* end) {
    if (n < 64) {
        U copy[64];
        for (U i = 0; i < n; ++i) copy[i] = base | load16(src + i);
        std::sort(copy, copy + n);
        std::memcpy(dst, copy, n * sizeof(U));
        return;
    }
    if (n > 8192) {
        dense(src, n, dst, base);
        return;
    }
    alignas(64) static B columns[256 * 8192];
    alignas(64) U count[256];
    std::memset(columns, 255, 256 * Network);
    for (U i = 0; i < 256; ++i) count[i] = i;
    for (U i = 0; i < n; ++i) {
        U v = load16(src + i), h = v >> 8, pos = count[h];
        count[h] = pos + 256;
        columns[pos] = static_cast<B>(v);
    }
    U overflow = 0;
    for (U i = 0; i < 256; ++i) {
        count[i] >>= 8;
        overflow |= count[i] > 32;
    }
    if (overflow) {
        dense(src, n, dst, base);
        return;
    }
    for (U block = 0; block < 256; block += 32) {
        alignas(32) B first[256], tail[128];
        batch_sort(columns + block, first, tail);
        V prefix = _mm256_set1_epi32(base | (block << 8));
        for (U i = 0; i < 32; ++i) {
            U h = block + i, c = count[h], p = base | (h << 8);
            W x = _mm_loadl_epi64(reinterpret_cast<const W*>(first + 8*i));
            U extra = 0;
            if constexpr (Network > 8)
                std::memcpy(&extra, tail + (Network - 8)*i, Network - 8);
            if (DUCK_EXPECT(c <= Network, 1)) {
                V wide = _mm256_or_si256(_mm256_cvtepu8_epi32(x), prefix);
                if (DUCK_EXPECT(end - dst >= Network, 1)) {
#if defined(__GNUC__)
                    __asm__ volatile("vmovdqu {%1, %0|%0, %1}"
                        : "=m"(*reinterpret_cast<__m256i_u*>(dst)) : "x"(wide));
#else
                    _mm256_storeu_si256(reinterpret_cast<V*>(dst), wide);
#endif
                    if constexpr (Network > 8) {
                        W y = _mm_or_si128(_mm_cvtepu8_epi32(_mm_cvtsi32_si128(extra)), _mm256_castsi256_si128(prefix));
                        if constexpr (Network == 10)
                            _mm_storel_epi64(reinterpret_cast<W*>(dst + 8), y);
                        else
                            _mm_storeu_si128(reinterpret_cast<W*>(dst + 8), y);
                    }
                } else {
                    alignas(32) U copy[12];
                    _mm256_store_si256(reinterpret_cast<V*>(copy), wide);
                    if constexpr (Network > 8) {
                        W y = _mm_or_si128(_mm_cvtepu8_epi32(_mm_cvtsi32_si128(extra)), _mm256_castsi256_si128(prefix));
                        _mm_store_si128(reinterpret_cast<W*>(copy + 8), y);
                    }
                    std::memcpy(dst, copy, c * sizeof(U));
                }
            } else if (c <= 16) {
                constexpr Z sentinel = Network == 8 ? ~Z(0) : Network == 10 ? 0xffffffffffff0000ull : 0xffffffff00000000ull;
                x = _mm_unpacklo_epi64(x, _mm_cvtsi64_si128(static_cast<long long>(Z(extra) | sentinel)));
                for (U j = Network; j < c; ++j) x = insert_byte(x, columns[256*j + h]);
                write8(x, dst, 8, p, end);
                write8(_mm_srli_si128(x, 8), dst + 8, c - 8, p, end);
            } else {
                B copy[32];
                for (U j = 0; j < c; ++j) copy[j] = columns[256*j + h];
                std::sort(copy, copy + c);
                for (U j = 0; j < c; ++j) dst[j] = p | copy[j];
            }
            dst += c;
            prefix = _mm256_add_epi32(prefix, _mm256_set1_epi32(256));
        }
    }
}

DUCK_INLINE bool fits(B** p, B** end) {
    V bad = _mm256_setzero_si256();
    for (U k = 0; k < 256; k += 4)
        bad = _mm256_or_si256(bad, _mm256_cmpgt_epi64(
            _mm256_load_si256(reinterpret_cast<const V*>(p + k)),
            _mm256_load_si256(reinterpret_cast<const V*>(end + k))));
    return _mm256_testz_si256(bad, bad) != 0;
}

bool split(const U* a, U n, const U* capacity, B** start, U* count) {
    alignas(64) B cache[256 * BlockBytes] = {};
    alignas(64) B* p[256];
    alignas(64) B* end[256];
    U position[256];
    std::size_t offset = 0;
    for (U k = 0; k < 256; ++k) {
        start[k] = p[k] = workspace + offset;
        end[k] = p[k] + ((std::size_t(capacity[k]) + Records - 1) / Records) * BlockBytes;
        offset = std::size_t(end[k] - workspace) + GuardBytes;
        position[k] = k * BlockBytes;
    }
    auto push = [&](U x) DUCK_LAMBDA_INLINE {
        U k = x >> 24, t = position[k];
        std::memcpy(cache + t, &x, sizeof(x));
        t += 3;
        if (DUCK_EXPECT((t & (BlockBytes - 1)) == Records * 3, 0)) {
            t -= Records * 3;
            repeat<BlockBytes / 32>([&](auto j) DUCK_LAMBDA_INLINE {
                _mm256_stream_si256(reinterpret_cast<V*>(p[k] + 32*j),
                    _mm256_load_si256(reinterpret_cast<const V*>(cache + t + 32*j)));
            });
            p[k] += BlockBytes;
        }
        position[k] = t;
    };
    for (U i = 0; i < n;) {
        U stop = std::min(n, i + 8192);
        for (; i + 8 <= stop; i += 8)
            repeat<8>([&](auto j) DUCK_LAMBDA_INLINE { push(a[i + j]); });
        for (; i < stop; ++i) push(a[i]);
        if (!fits(p, end)) {
            _mm_sfence();
            return false;
        }
    }
    for (U k = 0; k < 256; ++k) {
        U r = position[k] & (BlockBytes - 1);
        count[k] = U((p[k] - start[k]) / BlockBytes) * Records + r / 3;
        if (r) {
            std::memcpy(p[k], cache + k * BlockBytes, BlockBytes);
            p[k] += BlockBytes;
        }
    }
    _mm_sfence();
    return fits(p, end);
}

void middle(const B* src, U n, U base, U* dst) {
    if (n < 4096) {
        for (U i = 0; i < n; ++i) dst[i] = base | packed_at(src, i);
        std::sort(dst, dst + n);
        return;
    }
    U count[256] = {}, begin[257], position[256];
    if (n > 600000 || n < 16384) {
        for (U i = 0; i < n; ++i) ++count[packed_at(src, i) >> 16];
        begin[0] = 0;
        for (U k = 0; k < 256; ++k) {
            position[k] = begin[k];
            begin[k + 1] = begin[k] + count[k];
        }
        for (U i = 0; i < n; ++i) {
            U x = packed_at(src, i);
            H low = static_cast<H>(x);
            std::memcpy(reinterpret_cast<B*>(dst) + 2 * position[x >> 16]++, &low, sizeof(low));
        }
        for (int k = 255; k >= 0; --k)
            if (count[k]) leaf(reinterpret_cast<const H*>(reinterpret_cast<const B*>(dst) + 2*begin[k]),
                count[k], dst + begin[k], base | (U(k) << 16), dst + begin[k + 1]);
        return;
    }
    B* tmp = workspace + ScratchOffset;
    B* p[256];
    U bytes = 0;
    for (U k = 0; k < 256; ++k) {
        begin[k] = bytes;
        p[k] = tmp + bytes;
        bytes += 64 * (((((n + 255) / 256 * 5 / 4 + 32) * 2 + 63) / 64) | 1);
    }
    begin[256] = bytes;
    for (;;) {
        U i = 0;
        const B* block = src;
        for (; i + Records <= n; i += Records, block += BlockBytes) {
            _mm_prefetch(reinterpret_cast<const char*>(block + BlockBytes), _MM_HINT_T0);
            constexpr U Unroll = Records % 12 == 0 ? 12 : Records % 7 == 0 ? 7 : 5;
            U j = 0;
#if defined(__GNUC__)
#pragma GCC unroll 1
#endif
            for (; j + Unroll <= Records; j += Unroll)
                repeat<Unroll>([&](auto q) DUCK_LAMBDA_INLINE {
                    U x = load24(block + 3 * (j + q));
                    H low = static_cast<H>(x);
                    std::memcpy(p[x >> 16], &low, sizeof(low));
                    p[x >> 16] += 2;
                });
            for (; j < Records; ++j) {
                U x = load24(block + 3*j);
                H low = static_cast<H>(x);
                std::memcpy(p[x >> 16], &low, sizeof(low));
                p[x >> 16] += 2;
            }
        }
        for (; i < n; ++i) {
            U x = packed_at(src, i);
            H low = static_cast<H>(x);
            std::memcpy(p[x >> 16], &low, sizeof(low));
            p[x >> 16] += 2;
        }
        bool overflow = false;
        for (U k = 0; k < 256; ++k) {
            count[k] = U(p[k] - (tmp + begin[k])) / 2;
            overflow |= p[k] > tmp + begin[k + 1];
        }
        if (!overflow) break;
        bytes = 0;
        for (U k = 0; k < 256; ++k) {
            begin[k] = bytes;
            p[k] = tmp + bytes;
            bytes += 2 * count[k];
        }
        begin[256] = bytes;
    }
    U* end = dst + n;
    for (U k = 0; k < 256; ++k) {
        if (count[k]) leaf(reinterpret_cast<const H*>(tmp + begin[k]), count[k], dst, base | (k << 16), end);
        dst += count[k];
    }
}
}

void DUCK_SORT_FUNCTION(unsigned* a, int n) {
    if (n < 2) return;
    if (n < 4096 || n > 100000000) {
        std::sort(a, a + n);
        return;
    }
    void* raw = std::malloc(WorkspaceBytes + 63);
    if (!raw) {
        std::sort(a, a + n);
        return;
    }
    workspace = reinterpret_cast<B*>((reinterpret_cast<std::uintptr_t>(raw) + 63) & ~std::uintptr_t(63));
    U capacity[256], count[256];
    B* start[256];
    for (U k = 0; k < 256; ++k) capacity[k] = U((Z(n) + 255) / 256) * 6 / 5 + 32;
    if (!split(a, U(n), capacity, start, count)) {
        std::memset(capacity, 0, sizeof(capacity));
        for (int i = 0; i < n; ++i) ++capacity[a[i] >> 24];
        split(a, U(n), capacity, start, count);
    }
    U offset = 0;
    for (U k = 0; k < 256; ++k) {
        if (count[k]) middle(start[k], count[k], k << 24, a + offset);
        offset += count[k];
    }
    std::free(raw);
}

评测结果