代码

#include <string.h>
typedef unsigned int u32;
typedef unsigned char u8;

static u32 tmp[(1u<<27) + (1u<<22)] __attribute__((aligned(64)));
static u32 cnt[512];
static u32 rbase[512];
static u32 pbase[512];
static u32 off[512];
static u8 pos[512];
static u32 staging[512 * 24] __attribute__((aligned(64)));
static u32 scratch[1u<<21] __attribute__((aligned(64)));
static u32 scratch2[8192] __attribute__((aligned(64)));
static u32 bc[1024];
static u32 bc2[256];
static u32 sstart[1024];
static u32 spos[1024];
static u32 scount[1024];

#define THRESH (1u<<18)

static void scatter_top(u32 *src, u32 *dst, int n) {
    int i, j;
    for (j = 0; j < 512; j++) cnt[j] = 0;
    int n8 = n & ~7;
    for (i = 0; i < n8; i += 8) {
        cnt[src[i]>>23]++; cnt[src[i+1]>>23]++; cnt[src[i+2]>>23]++; cnt[src[i+3]>>23]++;
        cnt[src[i+4]>>23]++; cnt[src[i+5]>>23]++; cnt[src[i+6]>>23]++; cnt[src[i+7]>>23]++;
    }
    for (; i < n; i++) cnt[src[i] >> 23]++;
    { u32 t = 0, p = 0; for (j = 0; j < 512; j++) { u32 c = cnt[j]; rbase[j] = t; pbase[j] = p; t += c; p = (p + c + 23) & ~23u; } }
    for (j = 0; j < 512; j++) { off[j] = pbase[j]; pos[j] = 0; }
#define S1(X) do { u32 x_ = (X); u32 b_ = x_ >> 23; u32 p_ = pos[b_]; \
    staging[b_ * 24 + p_] = x_; p_++; \
    if (p_ == 24) { u32 go_ = off[b_]; __builtin_memcpy(dst+go_, staging+b_*24, 96); off[b_] = go_ + 24; 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 < 512; j++) { u32 p = pos[j]; if (p > 0) memcpy(dst+off[j], staging+j*24, p*4); }
}

// sort16: sort by low 8 bits then (mask2) bits; fused
static void sort16(u32 *src, u32 *dst, int m, u32 mask2) {
    int i, j;
    for (j = 0; j < 256; j++) bc[j] = 0;
    int m8 = m & ~7;
    for (i = 0; i < m8; i += 8) {
        bc[src[i]&255]++; bc[src[i+1]&255]++; bc[src[i+2]&255]++; bc[src[i+3]&255]++;
        bc[src[i+4]&255]++; bc[src[i+5]&255]++; bc[src[i+6]&255]++; bc[src[i+7]&255]++;
    }
    for (; i < m; i++) bc[src[i] & 255]++;
    for (j = 1; j < 256; j++) bc[j] += bc[j-1];
    for (j = 0; j < 256; j++) bc2[j] = 0;
    int rm = m & ~15;
    for (i = m - 1; i >= rm; i--) { u32 x = src[i]; scratch2[--bc[x & 255]] = x; bc2[(x >> 8) & mask2]++; }
    for (i = rm - 1; i >= 0; i -= 8) {
        u32 x0=src[i],x1=src[i-1],x2=src[i-2],x3=src[i-3],x4=src[i-4],x5=src[i-5],x6=src[i-6],x7=src[i-7];
        scratch2[--bc[x0&255]]=x0; bc2[(x0>>8)&mask2]++;
        scratch2[--bc[x1&255]]=x1; bc2[(x1>>8)&mask2]++;
        scratch2[--bc[x2&255]]=x2; bc2[(x2>>8)&mask2]++;
        scratch2[--bc[x3&255]]=x3; bc2[(x3>>8)&mask2]++;
        scratch2[--bc[x4&255]]=x4; bc2[(x4>>8)&mask2]++;
        scratch2[--bc[x5&255]]=x5; bc2[(x5>>8)&mask2]++;
        scratch2[--bc[x6&255]]=x6; bc2[(x6>>8)&mask2]++;
        scratch2[--bc[x7&255]]=x7; bc2[(x7>>8)&mask2]++;
    }
    for (j = 1; j < 256; j++) bc2[j] += bc2[j-1];
    for (i = m - 1; i >= rm; i--) { u32 x = scratch2[i]; dst[--bc2[(x >> 8) & mask2]] = x; }
    for (i = rm - 1; i >= 0; i -= 8) {
        u32 x0=scratch2[i],x1=scratch2[i-1],x2=scratch2[i-2],x3=scratch2[i-3],x4=scratch2[i-4],x5=scratch2[i-5],x6=scratch2[i-6],x7=scratch2[i-7];
        dst[--bc2[(x0>>8)&mask2]]=x0;
        dst[--bc2[(x1>>8)&mask2]]=x1;
        dst[--bc2[(x2>>8)&mask2]]=x2;
        dst[--bc2[(x3>>8)&mask2]]=x3;
        dst[--bc2[(x4>>8)&mask2]]=x4;
        dst[--bc2[(x5>>8)&mask2]]=x5;
        dst[--bc2[(x6>>8)&mask2]]=x6;
        dst[--bc2[(x7>>8)&mask2]]=x7;
    }
}

// level-2 with `nb` sub-buckets (8-bit: nb=256 shift=15; 10-bit: nb=1024 shift=13), then sort16(rem)
static void sort_bucket_adapt(u32 *src, u32 *dst, int m, int nb, int shift, u32 mask2) {
    int i, j;
    int nm = nb - 1;
    for (j = 0; j < nb; j++) bc[j] = 0;
    int mb8 = m & ~7;
    for (i = 0; i < mb8; i += 8) {
        bc[(src[i]>>shift)&nm]++; bc[(src[i+1]>>shift)&nm]++; bc[(src[i+2]>>shift)&nm]++; bc[(src[i+3]>>shift)&nm]++;
        bc[(src[i+4]>>shift)&nm]++; bc[(src[i+5]>>shift)&nm]++; bc[(src[i+6]>>shift)&nm]++; bc[(src[i+7]>>shift)&nm]++;
    }
    for (; i < m; i++) bc[(src[i] >> shift) & nm]++;
    { u32 t = 0; for (j = 0; j < nb; 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],b0=(x0>>shift)&nm; scratch[spos[b0]++]=x0;
        u32 x1=src[i+1],b1=(x1>>shift)&nm; scratch[spos[b1]++]=x1;
        u32 x2=src[i+2],b2=(x2>>shift)&nm; scratch[spos[b2]++]=x2;
        u32 x3=src[i+3],b3=(x3>>shift)&nm; scratch[spos[b3]++]=x3;
        u32 x4=src[i+4],b4=(x4>>shift)&nm; scratch[spos[b4]++]=x4;
        u32 x5=src[i+5],b5=(x5>>shift)&nm; scratch[spos[b5]++]=x5;
        u32 x6=src[i+6],b6=(x6>>shift)&nm; scratch[spos[b6]++]=x6;
        u32 x7=src[i+7],b7=(x7>>shift)&nm; scratch[spos[b7]++]=x7;
    }
    for (; i < m; i++) { u32 x = src[i]; u32 b = (x >> shift) & nm; scratch[spos[b]++] = x; }
    for (j = 0; j < nb; j++) { int c = scount[j]; if (c > 0) sort16(scratch + sstart[j], dst + sstart[j], c, mask2); }
}

void sort(unsigned *aa, int n) {
    u32 *a = (u32*)aa;
    int j;
    scatter_top(a, tmp, n);
    for (j = 0; j < 512; j++) {
        int m = (int)cnt[j];
        if (m <= 0) continue;
        if (m <= THRESH) {
            sort_bucket_adapt(tmp + pbase[j], a + rbase[j], m, 256, 15, 127);
        } else {
            sort_bucket_adapt(tmp + pbase[j], a + rbase[j], m, 1024, 13, 31);
        }
    }
}

评测结果