代码

// wc2017b1 v13: COMPRESSED-INTERMEDIATE [16][16] MSD radix sort.
// Pass1: histogram of top-16 bits (scalar). Pass2: compressed scatter -> u16 scratch
// (stores only low 16 bits; top16 recovered from region). Pass3: per-bucket (6KB L1-resident)
// 2x8-bit LSD sort of the u16s + sequential NT reconstruct to a[].
// DRAM traffic ~3.2GB (vs v12's ~4.0GB).
#include <immintrin.h>
#include <string.h>
#include <emmintrin.h>
typedef unsigned int u32;
typedef unsigned short u16;
typedef unsigned char u8;

#define NB 65536

static u16 scratch16[200000000] __attribute__((aligned(64)));
static u32 cnt16[NB];
static u32 B[NB];
static u32 used[NB];
static u16 buf[4096] __attribute__((aligned(64)));
static u16 buf2[4096] __attribute__((aligned(64)));
static u32 bc[256];
static u32 bc2[256];

__attribute__((target("avx2")))
static inline void flush16_nt(u32 *dst, const u16 *src, int m, u32 top16) {
    int mi = 0;
    while ((((unsigned long)(dst + mi)) & 31) && mi < m) { dst[mi] = ((u32)src[mi]) | top16; mi++; }
    __m256i t = _mm256_set1_epi32((int)top16);
    for (; mi + 16 <= m; mi += 16) {
        __m128i lo = _mm_loadu_si128((const __m128i*)(src + mi));
        __m128i hi = _mm_loadu_si128((const __m128i*)(src + mi + 8));
        __m256i a0 = _mm256_cvtepu16_epi32(lo);
        __m256i a1 = _mm256_cvtepu16_epi32(hi);
        _mm256_stream_si256((__m256i*)(dst + mi + 0), _mm256_or_si256(a0, t));
        _mm256_stream_si256((__m256i*)(dst + mi + 8), _mm256_or_si256(a1, t));
    }
    for (; mi < m; mi++) dst[mi] = ((u32)src[mi]) | top16;
}

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

    // PASS 1: histogram of top-16 bits
    for (k = 0; k < NB; k++) cnt16[k] = 0;
    int n8 = n & ~3;
    for (i = 0; i < n8; i += 4) {
        cnt16[a[i+0]>>16]++;
        cnt16[a[i+1]>>16]++;
        cnt16[a[i+2]>>16]++;
        cnt16[a[i+3]>>16]++;
    }
    for (; i < n; i++) cnt16[a[i] >> 16]++;

    // prefix -> B[k] start offsets
    { u32 s = 0; for (k = 0; k < NB; k++) { u32 c = cnt16[k]; B[k] = s; s += c; } }

    // PASS 2: compressed scatter (store low 16 bits only)
    for (k = 0; k < NB; k++) used[k] = B[k];
    for (i = 0; i < n8; i += 4) {
        u32 x0=a[i+0]; u32 k0=x0>>16; scratch16[used[k0]++]=(u16)x0;
        u32 x1=a[i+1]; u32 k1=x1>>16; scratch16[used[k1]++]=(u16)x1;
        u32 x2=a[i+2]; u32 k2=x2>>16; scratch16[used[k2]++]=(u16)x2;
        u32 x3=a[i+3]; u32 k3=x3>>16; scratch16[used[k3]++]=(u16)x3;
    }
    for (; i < n; i++) { u32 x = a[i]; u32 kk = x >> 16; scratch16[used[kk]++] = (u16)x; }

    // PASS 3: per-bucket sort + sequential reconstruct
    u32 out = 0;
    for (k = 0; k < NB; k++) {
        u32 m = cnt16[k];
        if (m == 0) continue;
        u16 *region = scratch16 + B[k];
        // 2x8-bit LSD sort of the u16 region
        int j;
        for (j = 0; j < 256; j++) bc[j] = 0;
        int m8 = (int)m & ~3;
        for (i = 0; i < m8; i += 4) {
            bc[region[i+0]&255]++;
            bc[region[i+1]&255]++;
            bc[region[i+2]&255]++;
            bc[region[i+3]&255]++;
        }
        for (; i < (int)m; i++) bc[region[i] & 255]++;
        for (j = 1; j < 256; j++) bc[j] += bc[j-1];
        for (j = 0; j < 256; j++) bc2[j] = 0;
        int rm = (int)m & ~3;
        for (i = (int)m - 1; i >= rm; i--) { u16 x = region[i]; buf[--bc[x & 255]] = x; bc2[(x >> 8) & 255]++; }
        for (i = rm - 1; i >= 0; i -= 4) {
            u16 x0=region[i-0]; buf[--bc[x0&255]]=x0; bc2[(x0>>8)&255]++;
            u16 x1=region[i-1]; buf[--bc[x1&255]]=x1; bc2[(x1>>8)&255]++;
            u16 x2=region[i-2]; buf[--bc[x2&255]]=x2; bc2[(x2>>8)&255]++;
            u16 x3=region[i-3]; buf[--bc[x3&255]]=x3; bc2[(x3>>8)&255]++;
        }
        for (j = 1; j < 256; j++) bc2[j] += bc2[j-1];
        for (i = (int)m - 1; i >= rm; i--) { u16 x = buf[i]; buf2[--bc2[(x >> 8) & 255]] = x; }
        for (i = rm - 1; i >= 0; i -= 4) {
            u16 x0=buf[i-0]; buf2[--bc2[(x0>>8)&255]]=x0;
            u16 x1=buf[i-1]; buf2[--bc2[(x1>>8)&255]]=x1;
            u16 x2=buf[i-2]; buf2[--bc2[(x2>>8)&255]]=x2;
            u16 x3=buf[i-3]; buf2[--bc2[(x3>>8)&255]]=x3;
        }
        flush16_nt(a + out, buf2, (int)m, (u32)k << 16);
        out += m;
    }
}

评测结果