代码

// 1001 v15: [8][8][16] MSD + PACKED 3-byte compress (4 elements -> 3 u32 stores) + direct scatter.
// Pass 2: staging gathers 16 same-bucket elements; flush packs each 4 elements' low-24 bits into
//         3 u32 stores (shifts+ors, 0.75 stores/elem) written as 12 contiguous bytes to a packed
//         3-byte scratch s[3N]. Replaces v10/v14's 3 byte-stores/elem.
// Pass 3+4: identical to v10 (hist byte2 + reconstruct u32 -> a[], then 2x8-bit low-16 sort).
#include <string.h>
typedef unsigned int u32;
typedef unsigned char u8;
typedef unsigned short u16;

#ifndef N
#define N 100000000
#endif

#define SUBMAX 8192
#define W 16

static u8 s[3*N + 256*4] __attribute__((aligned(64)));
static u32 stg[256 * W] __attribute__((aligned(64)));
static u8  pos[256];
static u32 cur[256];       // byte cursor in s (advances 3/elem)
static u32 cnt[256];
static u32 B[257];
static u32 c2[256];
static u32 sstart[256];
static u32 scurs[256];
static u32 h[256];
static u32 h2[256];
static u32 tmpbuf[SUBMAX] __attribute__((aligned(64)));

// pack 4 elements (low-24) into 12 bytes at dst
static inline void pack4(u8 *dst, const u32 *e) {
    u32 w0 = (e[0] & 0xFFFFFF) | (e[1] << 24);
    u32 w1 = ((e[1] & 0xFFFFFF) >> 8) | ((e[2] & 0xFFFFFF) << 16);
    u32 w2 = ((e[2] & 0xFFFFFF) >> 16) | ((e[3] & 0xFFFFFF) << 8);
    __builtin_memcpy(dst,      &w0, 4);
    __builtin_memcpy(dst + 4,  &w1, 4);
    __builtin_memcpy(dst + 8,  &w2, 4);
}

// flush a full W-element staging block of bucket k
static inline void flush_full(int k) {
    u32 b = cur[k];
    const u32 *e = stg + k * W;
    pack4(s + b + 0,  e + 0);
    pack4(s + b + 12, e + 4);
    pack4(s + b + 24, e + 8);
    pack4(s + b + 36, e + 12);
    cur[k] = b + 48;
}

static inline void flush_partial(int k, int m) {
    u32 b = cur[k];
    for (int j = 0; j < m; j++) {
        u32 x = stg[k * W + j];
        s[b] = (u8)x; s[b+1] = (u8)(x >> 8); s[b+2] = (u8)(x >> 16);
        b += 3;
    }
    cur[k] = b;
}

// stable 2x8-bit LSD sort of low-16 of m u32 in base[] (top-16 fixed)
static void sort_low16(u32 *base, u32 m) {
    int i, j;
    if (m <= 1) return;
    for (j = 0; j < 256; j++) h[j] = 0;
    int m8 = (int)m & ~7;
    for (i = 0; i < m8; i += 8) {
        h[base[i]&255]++; h[base[i+1]&255]++; h[base[i+2]&255]++; h[base[i+3]&255]++;
        h[base[i+4]&255]++; h[base[i+5]&255]++; h[base[i+6]&255]++; h[base[i+7]&255]++;
    }
    for (; i < (int)m; i++) h[base[i] & 255]++;
    for (j = 1; j < 256; j++) h[j] += h[j-1];
    int rm = (int)m & ~7;
    for (i = (int)m - 1; i >= rm; i--) { u32 x = base[i]; tmpbuf[--h[x & 255]] = x; }
    for (i = rm - 1; i >= 0; i -= 8) {
        u32 x0=base[i],x1=base[i-1],x2=base[i-2],x3=base[i-3],x4=base[i-4],x5=base[i-5],x6=base[i-6],x7=base[i-7];
        tmpbuf[--h[x0&255]]=x0; tmpbuf[--h[x1&255]]=x1; tmpbuf[--h[x2&255]]=x2; tmpbuf[--h[x3&255]]=x3;
        tmpbuf[--h[x4&255]]=x4; tmpbuf[--h[x5&255]]=x5; tmpbuf[--h[x6&255]]=x6; tmpbuf[--h[x7&255]]=x7;
    }
    for (j = 0; j < 256; j++) h2[j] = 0;
    for (i = 0; i < m8; i += 8) {
        h2[(tmpbuf[i]>>8)&255]++; h2[(tmpbuf[i+1]>>8)&255]++; h2[(tmpbuf[i+2]>>8)&255]++; h2[(tmpbuf[i+3]>>8)&255]++;
        h2[(tmpbuf[i+4]>>8)&255]++; h2[(tmpbuf[i+5]>>8)&255]++; h2[(tmpbuf[i+6]>>8)&255]++; h2[(tmpbuf[i+7]>>8)&255]++;
    }
    for (; i < (int)m; i++) h2[(tmpbuf[i] >> 8) & 255]++;
    for (j = 1; j < 256; j++) h2[j] += h2[j-1];
    for (i = (int)m - 1; i >= rm; i--) { u32 x = tmpbuf[i]; base[--h2[(x >> 8) & 255]] = x; }
    for (i = rm - 1; i >= 0; i -= 8) {
        u32 x0=tmpbuf[i],x1=tmpbuf[i-1],x2=tmpbuf[i-2],x3=tmpbuf[i-3],x4=tmpbuf[i-4],x5=tmpbuf[i-5],x6=tmpbuf[i-6],x7=tmpbuf[i-7];
        base[--h2[(x0>>8)&255]]=x0; base[--h2[(x1>>8)&255]]=x1; base[--h2[(x2>>8)&255]]=x2; base[--h2[(x3>>8)&255]]=x3;
        base[--h2[(x4>>8)&255]]=x4; base[--h2[(x5>>8)&255]]=x5; base[--h2[(x6>>8)&255]]=x6; base[--h2[(x7>>8)&255]]=x7;
    }
}

void sort(unsigned *aa, int n) {
    u32 *a = (u32*)aa;
    int i, j;
    (void)n;

    // Pass 1: histogram top-8
    for (j = 0; j < 256; j++) cnt[j] = 0;
    int n8 = N & ~7;
    for (i = 0; i < n8; i += 8) {
        cnt[a[i]>>24]++; cnt[a[i+1]>>24]++; cnt[a[i+2]>>24]++; cnt[a[i+3]>>24]++;
        cnt[a[i+4]>>24]++; cnt[a[i+5]>>24]++; cnt[a[i+6]>>24]++; cnt[a[i+7]>>24]++;
    }
    for (; i < N; i++) cnt[a[i] >> 24]++;
    { u32 acc = 0; B[0] = 0; for (j = 0; j < 256; j++) { acc += cnt[j]; B[j+1] = acc; } }

    // Pass 2: packed compress via staging (16 elem/bucket) + 4-elem pack
    for (j = 0; j < 256; j++) { pos[j] = 0; cur[j] = B[j] * 3; }
    for (i = 0; i < n8; i += 8) {
        u32 x0=a[i+0],k0=x0>>24,p0=pos[k0]; stg[k0*W+p0]=x0; p0++; if(p0==W){flush_full(k0); p0=0;} pos[k0]=(u8)p0;
        u32 x1=a[i+1],k1=x1>>24,p1=pos[k1]; stg[k1*W+p1]=x1; p1++; if(p1==W){flush_full(k1); p1=0;} pos[k1]=(u8)p1;
        u32 x2=a[i+2],k2=x2>>24,p2=pos[k2]; stg[k2*W+p2]=x2; p2++; if(p2==W){flush_full(k2); p2=0;} pos[k2]=(u8)p2;
        u32 x3=a[i+3],k3=x3>>24,p3=pos[k3]; stg[k3*W+p3]=x3; p3++; if(p3==W){flush_full(k3); p3=0;} pos[k3]=(u8)p3;
        u32 x4=a[i+4],k4=x4>>24,p4=pos[k4]; stg[k4*W+p4]=x4; p4++; if(p4==W){flush_full(k4); p4=0;} pos[k4]=(u8)p4;
        u32 x5=a[i+5],k5=x5>>24,p5=pos[k5]; stg[k5*W+p5]=x5; p5++; if(p5==W){flush_full(k5); p5=0;} pos[k5]=(u8)p5;
        u32 x6=a[i+6],k6=x6>>24,p6=pos[k6]; stg[k6*W+p6]=x6; p6++; if(p6==W){flush_full(k6); p6=0;} pos[k6]=(u8)p6;
        u32 x7=a[i+7],k7=x7>>24,p7=pos[k7]; stg[k7*W+p7]=x7; p7++; if(p7==W){flush_full(k7); p7=0;} pos[k7]=(u8)p7;
    }
    for (; i < N; i++) {
        u32 x = a[i]; u32 k = x >> 24; u32 p = pos[k];
        stg[k*W + p] = x; p++;
        if (p == W) { flush_full(k); p = 0; }
        pos[k] = (u8)p;
    }
    for (j = 0; j < 256; j++) if (pos[j]) flush_partial(j, pos[j]);

    // Pass 3+4: per top-8 bucket
    for (int k = 0; k < 256; k++) {
        u32 m = cnt[k];
        if (m == 0) continue;
        u32 base3 = B[k] * 3;
        u32 baseA = B[k];
        int mb8 = (int)m & ~7;

        for (j = 0; j < 256; j++) c2[j] = 0;
        for (i = 0; i < mb8; i += 8) {
            c2[s[base3+3*i+2]]++; c2[s[base3+3*(i+1)+2]]++; c2[s[base3+3*(i+2)+2]]++; c2[s[base3+3*(i+3)+2]]++;
            c2[s[base3+3*(i+4)+2]]++; c2[s[base3+3*(i+5)+2]]++; c2[s[base3+3*(i+6)+2]]++; c2[s[base3+3*(i+7)+2]]++;
        }
        for (; i < (int)m; i++) c2[s[base3 + 3*i + 2]]++;
        { u32 acc = 0; for (j = 0; j < 256; j++) { sstart[j] = acc; scurs[j] = acc; acc += c2[j]; } }

        u32 kk = (u32)k << 24;
        for (i = 0; i < mb8; i += 8) {
            u32 b0 = base3 + 3*i;
            u32 g0 = s[b0+2]; u32 p0 = scurs[g0]++; a[baseA+p0] = kk | ((u32)g0<<16) | ((u32)s[b0+1]<<8) | (u32)s[b0];
            u32 b1 = b0 + 3;
            u32 g1 = s[b1+2]; u32 p1 = scurs[g1]++; a[baseA+p1] = kk | ((u32)g1<<16) | ((u32)s[b1+1]<<8) | (u32)s[b1];
            u32 b2 = b1 + 3;
            u32 g2 = s[b2+2]; u32 p2 = scurs[g2]++; a[baseA+p2] = kk | ((u32)g2<<16) | ((u32)s[b2+1]<<8) | (u32)s[b2];
            u32 b3 = b2 + 3;
            u32 g3 = s[b3+2]; u32 p3 = scurs[g3]++; a[baseA+p3] = kk | ((u32)g3<<16) | ((u32)s[b3+1]<<8) | (u32)s[b3];
            u32 b4 = b3 + 3;
            u32 g4 = s[b4+2]; u32 p4 = scurs[g4]++; a[baseA+p4] = kk | ((u32)g4<<16) | ((u32)s[b4+1]<<8) | (u32)s[b4];
            u32 b5 = b4 + 3;
            u32 g5 = s[b5+2]; u32 p5 = scurs[g5]++; a[baseA+p5] = kk | ((u32)g5<<16) | ((u32)s[b5+1]<<8) | (u32)s[b5];
            u32 b6 = b5 + 3;
            u32 g6 = s[b6+2]; u32 p6 = scurs[g6]++; a[baseA+p6] = kk | ((u32)g6<<16) | ((u32)s[b6+1]<<8) | (u32)s[b6];
            u32 b7 = b6 + 3;
            u32 g7 = s[b7+2]; u32 p7 = scurs[g7]++; a[baseA+p7] = kk | ((u32)g7<<16) | ((u32)s[b7+1]<<8) | (u32)s[b7];
        }
        for (; i < (int)m; i++) {
            u32 b = base3 + 3*i;
            u32 g = s[b+2];
            u32 p = scurs[g]++;
            a[baseA+p] = kk | ((u32)g<<16) | ((u32)s[b+1]<<8) | (u32)s[b];
        }

        for (j = 0; j < 256; j++) {
            u32 c = c2[j];
            if (c > 0) sort_low16(a + baseA + sstart[j], c);
        }
    }
}

评测结果