代码

// 1001 v13: [8][8][16] MSD + u16 sidecar, but the per-sub-bucket low-16 sort uses a
// 65536-bit bitmap (values are DISTINCT = permutation). Set bits for each low-16 value,
// iterate the bitmap ascending = sorted output. Replaces the 2x8-bit count+scatter sort16.
#include <immintrin.h>
#include <string.h>
#include <emmintrin.h>
typedef unsigned int u32;
typedef unsigned short u16;
typedef unsigned char u8;
typedef unsigned long long u64;

#define NB 256
#define LINE 128

static u32 tmp[100000000 + NB * LINE] __attribute__((aligned(64)));
static u32 staging[NB * LINE] __attribute__((aligned(64)));
static u32 cnt[NB];
static u32 rbase[NB];
static u32 pbase[NB];
static u32 off[NB];
static u8  pos[NB];
static u16 scratch[1048576] __attribute__((aligned(64)));
static u32 bc[256];
static u32 sstart[256];
static u32 spos[256];
static u32 scount[256];
static u64 bitmap[1024] __attribute__((aligned(64)));

__attribute__((target("avx2")))
static inline void flush_nt(u32 *dst, const u32 *src) {
    __m256i v0 = _mm256_loadu_si256((const __m256i*)(src + 0));
    __m256i v1 = _mm256_loadu_si256((const __m256i*)(src + 8));
    __m256i v2 = _mm256_loadu_si256((const __m256i*)(src + 16));
    __m256i v3 = _mm256_loadu_si256((const __m256i*)(src + 24));
    __m256i v4 = _mm256_loadu_si256((const __m256i*)(src + 32));
    __m256i v5 = _mm256_loadu_si256((const __m256i*)(src + 40));
    __m256i v6 = _mm256_loadu_si256((const __m256i*)(src + 48));
    __m256i v7 = _mm256_loadu_si256((const __m256i*)(src + 56));
    __m256i v8 = _mm256_loadu_si256((const __m256i*)(src + 64));
    __m256i v9 = _mm256_loadu_si256((const __m256i*)(src + 72));
    __m256i v10 = _mm256_loadu_si256((const __m256i*)(src + 80));
    __m256i v11 = _mm256_loadu_si256((const __m256i*)(src + 88));
    __m256i v12 = _mm256_loadu_si256((const __m256i*)(src + 96));
    __m256i v13 = _mm256_loadu_si256((const __m256i*)(src + 104));
    __m256i v14 = _mm256_loadu_si256((const __m256i*)(src + 112));
    __m256i v15 = _mm256_loadu_si256((const __m256i*)(src + 120));
    _mm256_stream_si256((__m256i*)(dst + 0), v0);
    _mm256_stream_si256((__m256i*)(dst + 8), v1);
    _mm256_stream_si256((__m256i*)(dst + 16), v2);
    _mm256_stream_si256((__m256i*)(dst + 24), v3);
    _mm256_stream_si256((__m256i*)(dst + 32), v4);
    _mm256_stream_si256((__m256i*)(dst + 40), v5);
    _mm256_stream_si256((__m256i*)(dst + 48), v6);
    _mm256_stream_si256((__m256i*)(dst + 56), v7);
    _mm256_stream_si256((__m256i*)(dst + 64), v8);
    _mm256_stream_si256((__m256i*)(dst + 72), v9);
    _mm256_stream_si256((__m256i*)(dst + 80), v10);
    _mm256_stream_si256((__m256i*)(dst + 88), v11);
    _mm256_stream_si256((__m256i*)(dst + 96), v12);
    _mm256_stream_si256((__m256i*)(dst + 104), v13);
    _mm256_stream_si256((__m256i*)(dst + 112), v14);
    _mm256_stream_si256((__m256i*)(dst + 120), v15);
}

static void scatter_top(u32 *src, u32 *dst, int n) {
    int i, j;
    for (j = 0; j < NB; j++) cnt[j] = 0;
    int n8 = n & ~7;
    for (i = 0; i < n8; i += 8) {
        cnt[src[i]>>24]++; cnt[src[i+1]>>24]++; cnt[src[i+2]>>24]++; cnt[src[i+3]>>24]++;
        cnt[src[i+4]>>24]++; cnt[src[i+5]>>24]++; cnt[src[i+6]>>24]++; cnt[src[i+7]>>24]++;
    }
    for (; i < n; i++) cnt[src[i] >> 24]++;
    { u32 t = 0, p = 0; for (j = 0; j < NB; j++) { u32 c = cnt[j]; rbase[j] = t; pbase[j] = p; t += c; p = (p + c + 127) & ~127u; } }
    for (j = 0; j < NB; j++) { off[j] = pbase[j]; pos[j] = 0; }
#define S1(X) do { u32 x_ = (X); u32 b_ = x_ >> 24; u32 p_ = pos[b_]; \
    staging[b_ * LINE + p_] = x_; p_++; \
    if (p_ == LINE) { u32 go_ = off[b_]; flush_nt(dst + go_, staging + b_ * LINE); off[b_] = go_ + LINE; p_ = 0; } \
    pos[b_] = (u8)p_; } while (0)
    int n8s = n & ~7;
    for (i = 0; i < n8s; i += 8) {
        S1(src[i]); S1(src[i+1]); S1(src[i+2]); S1(src[i+3]);
        S1(src[i+4]); S1(src[i+5]); S1(src[i+6]); S1(src[i+7]);
    }
    for (; i < n; i++) S1(src[i]);
#undef S1
    for (j = 0; j < NB; j++) { u32 p = pos[j]; if (p > 0) __builtin_memcpy(dst + off[j], staging + j * LINE, p * 4); }
}

// sort a sub-bucket of m distinct low-16 values (src = u16) and write u32 = top16|x to dst ascending.
static void sort_sub(u16 *src, u32 *dst, int m, u32 top16) {
    memset(bitmap, 0, sizeof(bitmap));
    int i;
    int m8 = m & ~7;
    for (i = 0; i < m8; i += 8) {
        bitmap[src[i]>>6] |= 1ull << (src[i]&63);
        bitmap[src[i+1]>>6] |= 1ull << (src[i+1]&63);
        bitmap[src[i+2]>>6] |= 1ull << (src[i+2]&63);
        bitmap[src[i+3]>>6] |= 1ull << (src[i+3]&63);
        bitmap[src[i+4]>>6] |= 1ull << (src[i+4]&63);
        bitmap[src[i+5]>>6] |= 1ull << (src[i+5]&63);
        bitmap[src[i+6]>>6] |= 1ull << (src[i+6]&63);
        bitmap[src[i+7]>>6] |= 1ull << (src[i+7]&63);
    }
    for (; i < m; i++) bitmap[src[i]>>6] |= 1ull << (src[i]&63);
    u32 r = 0;
    for (int w = 0; w < 1024; w++) {
        u64 word = bitmap[w];
        u32 base = (u32)(w * 64);
        while (word) {
            int b = __builtin_ctzll(word);
            dst[r++] = top16 | (base + (u32)b);
            word &= word - 1;
        }
    }
}

static void sort_bucket(u32 *src, u32 *dst, int m, int bj) {
    int i, j;
    for (j = 0; j < 256; j++) bc[j] = 0;
    int mb8 = m & ~7;
    for (i = 0; i < mb8; i += 8) {
        bc[(src[i]>>16)&255]++; bc[(src[i+1]>>16)&255]++; bc[(src[i+2]>>16)&255]++; bc[(src[i+3]>>16)&255]++;
        bc[(src[i+4]>>16)&255]++; bc[(src[i+5]>>16)&255]++; bc[(src[i+6]>>16)&255]++; bc[(src[i+7]>>16)&255]++;
    }
    for (; i < m; i++) bc[(src[i] >> 16) & 255]++;
    { u32 t = 0; for (j = 0; j < 256; j++) { u32 c = bc[j]; scount[j] = c; sstart[j] = t; spos[j] = t; t += c; } }

    for (i = 0; i < mb8; i += 8) {
        u32 x0=src[i+0],b0=(x0>>16)&255; scratch[spos[b0]++]=(u16)x0;
        u32 x1=src[i+1],b1=(x1>>16)&255; scratch[spos[b1]++]=(u16)x1;
        u32 x2=src[i+2],b2=(x2>>16)&255; scratch[spos[b2]++]=(u16)x2;
        u32 x3=src[i+3],b3=(x3>>16)&255; scratch[spos[b3]++]=(u16)x3;
        u32 x4=src[i+4],b4=(x4>>16)&255; scratch[spos[b4]++]=(u16)x4;
        u32 x5=src[i+5],b5=(x5>>16)&255; scratch[spos[b5]++]=(u16)x5;
        u32 x6=src[i+6],b6=(x6>>16)&255; scratch[spos[b6]++]=(u16)x6;
        u32 x7=src[i+7],b7=(x7>>16)&255; scratch[spos[b7]++]=(u16)x7;
    }
    for (; i < m; i++) { u32 x = src[i]; u32 b = (x >> 16) & 255; scratch[spos[b]++] = (u16)x; }
    for (j = 0; j < 256; j++) { int c = scount[j]; if (c > 0) sort_sub(scratch + sstart[j], dst + sstart[j], c, ((u32)((bj<<8)|j))<<16); }
}

void sort(unsigned *aa, int n) {
    u32 *a = (u32*)aa;
    int j;
    scatter_top(a, tmp, n);
    for (j = 0; j < NB; j++) {
        int m = (int)cnt[j];
        if (m > 0) sort_bucket(tmp + pbase[j], a + rbase[j], m, j);
    }
}

评测结果