代码

// Duck.ac 1001, C++17 / AVX2, static -O2.
// Architecture: fixed-expansion radix split + SIMD prefix/network tail.
//   split  : 256 top-byte buckets, 24-bit packed records, non-temporal block writes.
//   middle : 256 sub-buckets by bits 16..23 (scratch) or in place.
//   leaf   : 8-bit columns + 12-input 39-comparator network + shift/max/min insertion.
// Leaf/emit refinements follow https://duck.ac/submission/48078 (374 ms), which refines
// https://duck.ac/submission/48070.  12-input network from Bert Dobbelaere,
// https://bertdobbelaere.github.io/sorting_networks.html#N12L39D9
//
// Tunables (override with -D):
//   DUCK_BLOCK_BYTES  64 | 128 | 256   packed block size (cache array = 256*this bytes)
//   DUCK_RECORDS      21 |  42 |  84   records per block, must satisfy 3*R+1 <= BLOCK
//   DUCK_HUGE_PAGES   0 | 1            madvise(MADV_HUGEPAGE) on the workspace (Linux only)
//   DUCK_DIRECT_STORE 0 | 1            see below; default 1
// Default: direct store, 256-byte group, 84 records per group.
// -DDUCK_DIRECT_STORE=0 restores the block staging + non-temporal path.
// Exact overflow recovery is kept for every speculative bucket.
#include <algorithm>
#include <cstdint>
#include <cstring>
#include <immintrin.h>
#include <utility>
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#include <cstdlib>

#ifndef DUCK_BLOCK_BYTES
#define DUCK_BLOCK_BYTES 256
#endif
#ifndef DUCK_RECORDS
#define DUCK_RECORDS 84
#endif
#ifndef DUCK_HUGE_PAGES
#define DUCK_HUGE_PAGES 0
#endif
// Middle pass: how many packed blocks ahead to prefetch (0 disables).
#ifndef DUCK_PREFETCH_BLOCKS
#define DUCK_PREFETCH_BLOCKS 1
#endif
// Local diagnosis only: accumulate per-phase cycles and report at exit.
#ifndef DUCK_PHASE_TIMING
#define DUCK_PHASE_TIMING 0
#endif
// 0 = stage full blocks and write them with non-temporal stores (no read-for-ownership).
// 1 = write packed records straight into the output streams (fewer uops, but every
//     64-byte line is read from DRAM before it is written).
#ifndef DUCK_DIRECT_STORE
#define DUCK_DIRECT_STORE 0
#endif
#ifndef DUCK_SORT_FUNCTION
#define DUCK_SORT_FUNCTION sort
#endif

#if defined(__linux__) && DUCK_HUGE_PAGES
#include <sys/mman.h>
#endif

#if DUCK_PHASE_TIMING
#include <cstdio>
#if defined(__GNUC__)
#define DUCK_RDTSC() __builtin_ia32_rdtsc()
#else
#include <intrin.h>
#define DUCK_RDTSC() __rdtsc()
#endif
static unsigned long long duck_phase[3];
static bool duck_reported;
#define DUCK_TICK(v) unsigned long long v = DUCK_RDTSC()
#define DUCK_TOCK(i, v) duck_phase[i] += DUCK_RDTSC() - (v)
#define DUCK_REPORT() do { if (!duck_reported) { duck_reported = true; \
    std::fprintf(stderr, "[phase] split=%llu middle=%llu leaf=%llu\n", \
        duck_phase[0], duck_phase[1], duck_phase[2]); } } while (0)
#else
#define DUCK_TICK(v)
#define DUCK_TOCK(i, v)
#define DUCK_REPORT()
#endif

#if defined(__GNUC__)
#define DUCK_NOINLINE __attribute__((noinline))
#define DUCK_ALWAYS inline __attribute__((always_inline))
#define DUCK_LAMBDA_AI __attribute__((always_inline))
#define DUCK_BARRIER(v) __asm__("" : "+x"(v))
#define DUCK_EXPECT(x, v) __builtin_expect((x), (v))
#define DUCK_STORE_U(ptr, v) __asm__ volatile("vmovdqu {%1, %0|%0, %1}" \
    : "=m"(*reinterpret_cast<__m256i_u*>(ptr)) : "x"(v))
#else
#define DUCK_NOINLINE __declspec(noinline)
#define DUCK_ALWAYS __forceinline
#define DUCK_LAMBDA_AI
#define DUCK_BARRIER(v)
#define DUCK_EXPECT(x, v) (x)
#define DUCK_STORE_U(ptr, v) _mm256_storeu_si256(reinterpret_cast<__m256i*>(ptr), (v))
#endif

namespace fastsort {
template<class F, size_t... I>
DUCK_ALWAYS void each(F&& f, std::index_sequence<I...>) {
    (f(std::integral_constant<size_t, I>{}), ...);
}
template<size_t N, class F>
DUCK_ALWAYS void repeat(F&& f) {
    each(std::forward<F>(f), std::make_index_sequence<N>{});
}
}

namespace v2_tail {
using U = unsigned;
using B = std::uint8_t;
using H = std::uint16_t;
using V = __m256i;
using W = __m128i;
constexpr U BlockBytes = DUCK_BLOCK_BYTES;
constexpr U Records = DUCK_RECORDS;
static_assert(BlockBytes == 64 || BlockBytes == 128 || BlockBytes == 256 || BlockBytes == 512,
    "Unsupported block size");
static_assert(Records > 0 && Records * 3 + 1 <= BlockBytes, "Invalid packed block");
alignas(64) static B columns[256 * 8192];
// Middle-pass unroll: must divide Records, otherwise the last group would read
// past the packed block.
constexpr U Unroll = (Records % 12 == 0) ? 12 : (Records % 7 == 0) ? 7
    : (Records % 6 == 0) ? 6 : (Records % 5 == 0) ? 5
    : (Records % 4 == 0) ? 4 : (Records % 3 == 0) ? 3 : (Records % 2 == 0) ? 2 : 1;
DUCK_ALWAYS U load16(const H* p) {
    H v;
    std::memcpy(&v, p, 2);
    return v;
}
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_ALWAYS void write8(W x, U* dst, U base) {
    V y = _mm256_or_si256(_mm256_cvtepu8_epi32(x), _mm256_set1_epi32(base));
    DUCK_STORE_U(dst, y);
}
// One AVX2 byte lane per bucket: 12 ranks x 32 buckets, then a 32x12 byte transpose.
// 'first' receives ranks 0..7 (8 bytes per bucket), 'tail' receives ranks 8..11.
DUCK_ALWAYS void batch_sort(const B* src, B* first, B* tail) {
    V x[12], a[8], b[8];
    fastsort::repeat<12>([&](auto j) DUCK_LAMBDA_AI {
        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)
    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
#pragma GCC unroll 4
    for (U i = 0; i < 4; ++i) {
        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]);
    }
#pragma GCC unroll 2
    for (U i = 0; i < 2; ++i) {
#pragma GCC unroll 2
        for (U j = 0; j < 2; ++j) {
            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]);
        }
    }
#pragma GCC unroll 4
    for (U i = 0; i < 4; ++i) {
        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 + 64*i), lo);
        _mm256_store_si256(reinterpret_cast<V*>(first + 64*i + 32), hi);
    }
    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]);
#pragma GCC unroll 2
    for (U i = 0; i < 2; ++i) {
        _mm256_store_si256(reinterpret_cast<V*>(tail + 64*i), b[2*i]);
        _mm256_store_si256(reinterpret_cast<V*>(tail + 64*i + 32), b[2*i+1]);
    }
}
// x is sorted and has at least one trailing 255 sentinel.
DUCK_ALWAYS W insert_byte(W x, U byte) {
    return _mm_min_epu8(x, _mm_max_epu8(_mm_slli_si128(x, 1), _mm_set1_epi8(char(byte))));
}
DUCK_NOINLINE void slow_output(const B* first8, const B* last4, U h, U c, U* dst, U p) {
    if (c <= 16) {
        W x = _mm_loadl_epi64(reinterpret_cast<const W*>(first8));
        std::uint32_t four;
        std::memcpy(&four, last4, 4);
        x = _mm_unpacklo_epi64(x, _mm_cvtsi64_si128(std::uint64_t(four) | 0xffffffff00000000ull));
        U j = 12;
        do { x = insert_byte(x, columns[256*j + h]); } while (++j < c);
        write8(x, dst, p);
        write8(_mm_srli_si128(x, 8), dst + 8, p);
    } 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];
    }
}
// Sorts one sub-bucket of <= 8192 two-byte keys into dst.
// Writes in 16-element units, so it needs 16 writable padding elements after dst+n
// unless columns_merge routes through last_leaf.
DUCK_NOINLINE void columns_fast(const H* src, U n, U* dst, U base) {
    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 * 4);
        return;
    }
    if (n > 8192) { dense(src, n, dst, base); return; }
    alignas(64) U count[256];
    std::memset(columns, 255, 3072);
    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, at = count[h];
        count[h] = at + 256;
        columns[at] = B(v);
    }
    V over = _mm256_setzero_si256();
    for (U i = 0; i < 256; i += 8) {
        V c = _mm256_srli_epi32(_mm256_load_si256(reinterpret_cast<const V*>(count + i)), 8);
        _mm256_store_si256(reinterpret_cast<V*>(count + i), c);
        over = _mm256_or_si256(over, _mm256_cmpgt_epi32(c, _mm256_set1_epi32(32)));
    }
    if (!_mm256_testz_si256(over, over)) { dense(src, n, dst, base); return; }
    for (U block = 0; block < 256; block += 32) {
        alignas(32) B first[256], second[128];
        batch_sort(columns + block, first, second);
        V prefix = _mm256_set1_epi32(base | (block << 8));
        auto emit = [&](U i, const B* first8, const B* last4) DUCK_LAMBDA_AI {
            U h = block + i, c = count[h];
            if (DUCK_EXPECT(c <= 12, 1)) {
                V wide = _mm256_or_si256(_mm256_cvtepu8_epi32(_mm_loadl_epi64(reinterpret_cast<const W*>(first8))), prefix);
                DUCK_STORE_U(dst, wide);
                std::uint32_t four;
                std::memcpy(&four, last4, 4);
                W y = _mm_or_si128(_mm_cvtepu8_epi32(_mm_cvtsi32_si128(four)), _mm256_castsi256_si128(prefix));
                _mm_storeu_si128(reinterpret_cast<W*>(dst + 8), y);
            } else slow_output(first8, last4, h, c, dst, base | (h << 8));
            dst += c;
            prefix = _mm256_add_epi32(prefix, _mm256_set1_epi32(256));
            DUCK_BARRIER(prefix);
        };
#pragma GCC unroll 1
        for (U half = 0; half < 2; ++half) {
#pragma GCC unroll 1
            for (U group = 0; group < 4; ++group) {
                fastsort::repeat<4>([&](auto r) DUCK_LAMBDA_AI {
                    emit(16*half + 4*group + r,
                         first + 16*half + 64*group + 8*r + 16*(r/2),
                         second + 16*half + 32*group + 4*r);
                });
            }
        }
    }
}
// Buffered path for the last sub-buckets of a range and for in-place overlap.
DUCK_NOINLINE void last_leaf(const H* src, U n, U* dst, U base) {
    alignas(32) U copy[8192 + 16];
    columns_fast(src, n, copy, base);
    std::memcpy(dst, copy, n * 4);
}
template<bool Full = false>
DUCK_ALWAYS void columns_merge(const H* src, U n, U* dst, U base, U* end) {
    if (!Full && n >= 64 && n <= 8192 && end - dst < n + 16) last_leaf(src, n, dst, base);
    else columns_fast(src, n, dst, base);
}
}

namespace fastsort {
using U = unsigned;
using B = unsigned char;
using H = std::uint16_t;
using Z = std::uint64_t;
// Kept for the process lifetime and reused if sort() is called again.
static B* mem;
DUCK_ALWAYS U get(const B* p) {
    U x;
    std::memcpy(&x, p, 4);
    return x & 0xffffff;
}
// Stride between packed groups: padded blocks, or exactly Records*3 when the
// records are written straight into the stream.
constexpr U Stride = DUCK_DIRECT_STORE ? DUCK_RECORDS * 3 : DUCK_BLOCK_BYTES;
DUCK_ALWAYS U at(const B* p, U i) {
    return get(p + std::size_t(i / DUCK_RECORDS) * Stride + i % DUCK_RECORDS * 3);
}
DUCK_ALWAYS bool fits(B** p, B** e) {
    __m256i 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 __m256i*>(p + k)),
        _mm256_load_si256(reinterpret_cast<const __m256i*>(e + k))));
    return _mm256_testz_si256(bad, bad) != 0;
}
// First pass: bucket by the top byte and pack the low 24 bits.
// 25344 bytes of guard per bucket cover the largest advance possible between two
// fits() checks (8192 records = 24576 bytes) plus one partial block.
bool split(U* a, U n, U* cap, B** start, U* cnt) {
    constexpr U guard = 25344;
    constexpr U Block = DUCK_BLOCK_BYTES, Rec = DUCK_RECORDS;
    alignas(64) B cache[DUCK_DIRECT_STORE ? 1 : 256 * Block] = {};
    alignas(64) B* p[256];
    alignas(64) B* e[256];
#if !DUCK_DIRECT_STORE
    U pos[256];
#endif
    std::size_t off = 0;
    for (U k = 0; k < 256; ++k) {
        start[k] = p[k] = mem + off;
#if DUCK_DIRECT_STORE
        e[k] = p[k] + std::size_t(cap[k]) * 3;
#else
        e[k] = p[k] + ((std::size_t(cap[k]) + Rec - 1) / Rec) * Block;
        pos[k] = k * Block;
#endif
        off = std::size_t(e[k] - mem) + guard;
    }
#if DUCK_DIRECT_STORE
    auto push = [&](U x) DUCK_LAMBDA_AI {
        U k = x >> 24;
        B* q = p[k];
        std::memcpy(q, &x, 4);
        p[k] = q + 3;
    };
#else
    auto push = [&](U x) DUCK_LAMBDA_AI {
        U k = x >> 24, t = pos[k];
        std::memcpy(cache + t, &x, 4);
        t += 3;
        if (DUCK_EXPECT((t & (Block - 1)) == Rec * 3, 0)) {
            t -= Rec * 3;
            fastsort::repeat<DUCK_BLOCK_BYTES / 32>([&](auto j) DUCK_LAMBDA_AI {
                _mm256_stream_si256(reinterpret_cast<__m256i*>(p[k] + 32*j),
                    _mm256_load_si256(reinterpret_cast<const __m256i*>(cache + t + 32*j)));
            });
            p[k] += Block;
        }
        pos[k] = t;
    };
#endif
    for (U i = 0; i < n;) {
        U stop = std::min(n, i + 8192);
        for (; i + 8 <= stop; i += 8) fastsort::repeat<8>([&](auto j) DUCK_LAMBDA_AI { push(a[i + j]); });
        for (; i < stop; ++i) push(a[i]);
        if (!fits(p, e)) { _mm_sfence(); return false; }
    }
    for (U k = 0; k < 256; ++k) {
#if DUCK_DIRECT_STORE
        cnt[k] = U(p[k] - start[k]) / 3;
#else
        U r = pos[k] & (Block - 1);
        cnt[k] = U((p[k] - start[k]) / Block) * Rec + r / 3;
        if (r) {
            std::memcpy(p[k], cache + k * Block, Block);
            p[k] += Block;
        }
#endif
    }
    _mm_sfence();
    return fits(p, e);
}
void leaf(const B* s, U n, U base, U* d, U* end = nullptr) {
    DUCK_TICK(t);
    v2_tail::columns_merge(reinterpret_cast<const H*>(s), n, d, base, end ? end : d + n);
    DUCK_TOCK(2, t);
}
void middle(const B* s, U n, U base, U* d) {
    if (n < 4096) {
        for (U i = 0; i < n; ++i) d[i] = base | at(s, i);
        std::sort(d, d + n);
        return;
    }
    U cnt[256] = {}, begin[257], pos[256];
    if (n > 600000 || n < 16384) {
        for (U i = 0; i < n; ++i) ++cnt[at(s, i) >> 16];
        begin[0] = 0;
        for (U k = 0; k < 256; ++k) { pos[k] = begin[k]; begin[k + 1] = begin[k] + cnt[k]; }
        for (U i = 0; i < n; ++i) {
            U x = at(s, i);
            H lo = x;
            std::memcpy(reinterpret_cast<B*>(d) + 2*pos[x >> 16]++, &lo, 2);
        }
        // Reverse order: finalized leaves live above the not-yet-consumed keys.
        for (int k = 255; k >= 0; --k)
            if (cnt[k]) leaf(reinterpret_cast<const B*>(d) + 2*begin[k], cnt[k], base | (U(k) << 16), d + begin[k]);
        return;
    }
    B* tmp = mem + 380000000;
    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* b = s;
        for (; i + DUCK_RECORDS <= n; i += DUCK_RECORDS, b += Stride) {
#if DUCK_PREFETCH_BLOCKS > 0
            _mm_prefetch(reinterpret_cast<const char*>(b + DUCK_PREFETCH_BLOCKS * DUCK_BLOCK_BYTES), _MM_HINT_T0);
#endif
            U j = 0;
#pragma GCC unroll 1
            for (; j + v2_tail::Unroll <= DUCK_RECORDS; j += v2_tail::Unroll)
                fastsort::repeat<v2_tail::Unroll>([&](auto q) DUCK_LAMBDA_AI {
                    U x = get(b + 3*(j + q));
                    H lo = x;
                    std::memcpy(p[x >> 16], &lo, 2);
                    p[x >> 16] += 2;
                });
            for (; j < DUCK_RECORDS; ++j) {
                U x = get(b + 3*j);
                H lo = x;
                std::memcpy(p[x >> 16], &lo, 2);
                p[x >> 16] += 2;
            }
        }
        for (; i < n; ++i) {
            U x = at(s, i);
            H lo = x;
            std::memcpy(p[x >> 16], &lo, 2);
            p[x >> 16] += 2;
        }
        bool bad = false;
        for (U k = 0; k < 256; ++k) {
            cnt[k] = U(p[k] - (tmp + begin[k])) / 2;
            bad |= p[k] > tmp + begin[k+1];
        }
        if (!bad) break;
        bytes = 0;
        for (U k = 0; k < 256; ++k) { begin[k] = bytes; p[k] = tmp + bytes; bytes += 2*cnt[k]; }
        begin[256] = bytes;
    }
    // Source is independent scratch, so a leaf may run past its own end into the
    // not-yet-finalized suffix of this bucket.
    U* end = d + n;
    for (U k = 0; k < 256; ++k) {
        if (cnt[k]) leaf(tmp + begin[k], cnt[k], base | (k << 16), d, end);
        d += cnt[k];
    }
}
}

void DUCK_SORT_FUNCTION(unsigned* a, int n) {
    using namespace fastsort;
    if (n < 2) return;
    if (!mem) {
        void* raw = std::malloc(384000063);
        if (!raw) { std::sort(a, a + n); return; }
        mem = reinterpret_cast<B*>((reinterpret_cast<std::uintptr_t>(raw) + 63) & ~std::uintptr_t(63));
#if defined(__linux__) && DUCK_HUGE_PAGES
        ::madvise(mem, 384000000, MADV_HUGEPAGE);
#endif
    }
    U cap[256], cnt[256];
    B* start[256];
    for (U k = 0; k < 256; ++k) cap[k] = U((Z(n) + 255) / 256) * 6 / 5 + 32;
    DUCK_TICK(t0);
    if (!split(a, U(n), cap, start, cnt)) {
        std::memset(cap, 0, sizeof(cap));
        for (int i = 0; i < n; ++i) ++cap[a[i] >> 24];
        split(a, U(n), cap, start, cnt);
    }
    DUCK_TOCK(0, t0);
    DUCK_TICK(t1);
    U off = 0;
    for (U k = 0; k < 256; ++k) {
        if (cnt[k]) middle(start[k], cnt[k], k << 24, a + off);
        off += cnt[k];
    }
    DUCK_TOCK(1, t1);
    DUCK_REPORT();
}

评测结果