代码

/*
Duck.ac problem 1001 -- pdoom optimization, C++17.

Credits for code and optimization techniques used in this submission:

1. pdoom, submissions 92970 and 93376:
   Base radix pipeline, 24-bit workspace, bounded speculative capacities,
   bitmap/radix fallback for skewed data, and final-array boundary handling.
   https://duck.ac/submission/92970
   https://duck.ac/submission/93376

2. iMMIQ, submission 48133:
   The byte_sort namespace and supporting helpers are adapted from this code:
   byte-column storage; 32 concurrent 12-input SIMD sorting networks; transpose;
   vector count restoration; insertion for longer columns; and assembly output
   with overlapping stores and exact boundary handling. Also the separate
   byte-key/16-bit-payload loads used in the middle scatter.
   https://duck.ac/submission/48133

3. saffah_cc_v41_260924, submission 92971:
   Packed 3-byte records via overlapping 4-byte staging stores; 192-byte blocks
   with 320-byte staging stride; absolute byte cursors; post-increment boundary
   detection; and complete-cache-line non-temporal flushes. Reimplemented here
   with capacity checks before every flush and exact-count retry on overflow.
   https://duck.ac/submission/92971

4. saffah_cc_v41_260924, submission 93285:
   Write-touch the following input page before reading it, using OR-with-zero
   to preserve every value. Here it is scheduled by a 1024-element outer loop.
   The page preparation is part of the timed sort operation.
   https://duck.ac/submission/93285

5. Elo, submission 48186:
   Additional reference for independent byte-key and word-payload loads in
   the packed-24-bit middle scatter (also present in iMMIQ's implementation).
   https://duck.ac/submission/48186

6. Bert Dobbelaere, as credited by iMMIQ and Elo:
   The 12-input / 39-comparator / 9-layer sorting network.
   https://bertdobbelaere.github.io/sorting_networks.html#N12L39D9

The integration retains repeated values, exact fallback sorting, and bounded
output. No expected-answer table or input-specific output is used.
*/
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("avx2,bmi,bmi2,popcnt")
#include <immintrin.h>
#include <algorithm>
#include <cstring>
#include <cstdint>
#include <utility>
#ifndef BUFFER_BITS
#define BUFFER_BITS 6
#endif
#ifndef TOP_BITS
#define TOP_BITS 8
#endif
#ifndef UNROLL
#define UNROLL 1
#endif
static constexpr unsigned TB=TOP_BITS, K=1u<<TB, BS=1u<<BUFFER_BITS;
alignas(64) static unsigned char workspace[3*(106250000+K*BS)+16];
alignas(256) static unsigned char staging_bytes[K*320];
static unsigned byte_cursor[K];
alignas(64) static unsigned hist[K], starts[K], pos[K];

static void fallback(unsigned *s, unsigned *d, unsigned n) {
    constexpr unsigned Q=1u<<(16-TB);
    unsigned h0[256]={},h1[256]={},h2[Q]={};
    for(unsigned i=0;i<n;i++) {unsigned x=s[i]; ++h0[x&255]; ++h1[(x>>8)&255]; ++h2[(x>>16)&(Q-1)];}
    unsigned z=0;
    for(unsigned k=0;k<256;k++){unsigned t=h0[k];h0[k]=z;z+=t;}
    z=0;for(unsigned k=0;k<256;k++){unsigned t=h1[k];h1[k]=z;z+=t;}
    z=0;for(unsigned k=0;k<Q;k++){unsigned t=h2[k];h2[k]=z;z+=t;}
    for(unsigned i=0;i<n;i++){unsigned x=s[i];d[h0[x&255]++]=x;}
    for(unsigned i=0;i<n;i++){unsigned x=d[i];s[h1[(x>>8)&255]++]=x;}
    for(unsigned i=0;i<n;i++){unsigned x=s[i];d[h2[(x>>16)&(Q-1)]++]=x;}
}

static inline unsigned read24(const unsigned char *s){unsigned x;memcpy(&x,s,4);return x&0xffffff;}
static inline void write24(unsigned char *s,unsigned x){s[0]=x;s[1]=x>>8;s[2]=x>>16;}
static void fallback24(unsigned char *s,unsigned *d,unsigned n,unsigned prefix){
    unsigned h0[256]={},h1[256]={},h2[256]={};
    for(unsigned i=0;i<n;i++){unsigned x=read24(s+3*i);++h0[x&255];++h1[(x>>8)&255];++h2[x>>16];}
    unsigned z=0;for(unsigned k=0;k<256;k++){unsigned t=h0[k];h0[k]=z;z+=t;}
    z=0;for(unsigned k=0;k<256;k++){unsigned t=h1[k];h1[k]=z;z+=t;}
    z=0;for(unsigned k=0;k<256;k++){unsigned t=h2[k];h2[k]=z;z+=t;}
    for(unsigned i=0;i<n;i++){unsigned x=read24(s+3*i);d[h0[x&255]++]=x;}
    for(unsigned i=0;i<n;i++){unsigned x=d[i];write24(s+3*h1[(x>>8)&255]++,x);}
    for(unsigned i=0;i<n;i++){unsigned x=read24(s+3*i);d[h2[x>>16]++]=prefix|x;}
}
alignas(64) static uint64_t bitmap[1u<<(26-TB)];
alignas(64) static unsigned duplicates[8192], dup_unsorted[8192];
template<unsigned Bits,bool Bounded,bool WithDuplicates>
static unsigned emit(unsigned *d,unsigned n,unsigned prefix,unsigned nd){
    constexpr unsigned Mask=(1u<<Bits)-1;
    unsigned *out=d,*end=d+n;
    unsigned di=0;
    unsigned nextWord=nd ? ((duplicates[0]&Mask)>>6) : (1u<<(Bits-6));
    unsigned i=0;
    const unsigned words=(1u<<(Bits-6));
    #pragma GCC unroll 4
    for(;i<words && (!Bounded || end-out>=3);i++){
        uint64_t w=bitmap[i];bitmap[i]=0;
        unsigned base=prefix+(i<<6);
        if(WithDuplicates && i==nextWord){
            while(w){
                unsigned x=base+_tzcnt_u64(w);w=_blsr_u64(w);*out++=x;
                while(di<nd && duplicates[di]==x){*out++=x;++di;}
            }
            nextWord=di<nd ? ((duplicates[di]&Mask)>>6) : words;
            continue;
        }
        unsigned cnt=_mm_popcnt_u64(w);
        out[0]=base+_tzcnt_u64(w); w=_blsr_u64(w);
        out[1]=base+_tzcnt_u64(w); w=_blsr_u64(w);
        out[2]=base+_tzcnt_u64(w); w=_blsr_u64(w);
        unsigned j=3;
        while(w){out[j++]=base+_tzcnt_u64(w);w=_blsr_u64(w);}
        out+=cnt;
    }
    if constexpr(Bounded) for(;i<words;i++){
        uint64_t w=bitmap[i];bitmap[i]=0;unsigned base=prefix+(i<<6);
        while(w){unsigned x=base+_tzcnt_u64(w);*out++=x;w=_blsr_u64(w);while(di<nd&&duplicates[di]==x){*out++=x;++di;}}
    }

    return out-d;
}
static void small_original(unsigned char *s,unsigned *d,unsigned n,unsigned prefix,unsigned *array_end) {
    if(n<128){for(unsigned i=0;i<n;i++)d[i]=prefix|read24(s+3*i);std::sort(d,d+n);return;}
    constexpr unsigned Mask=(1u<<(32-TB))-1;
    unsigned nd=0;
    for(unsigned i=0;i<n;i++){
        unsigned x=read24(s+3*i);
        uint64_t old=bitmap[x>>6];bool existed;
        asm("btsq %2,%1" : "=@ccc"(existed), "+r"(old) : "r"(uint64_t(x)));
        bitmap[x>>6]=old;
        if(existed){
            if(nd==8192){memset(bitmap,0,sizeof(bitmap));fallback24(s,d,n,prefix);return;}
            dup_unsorted[nd++]=prefix|x;
        }
    }
    fallback(dup_unsorted,duplicates,nd);
    bool bounded=array_end-(d+n)<3;
    if(nd){if(bounded)emit<24,true,true>(d,n,prefix,nd);else emit<24,false,true>(d,n,prefix,nd);}
    else {if(bounded)emit<24,true,false>(d,n,prefix,nd);else emit<24,false,false>(d,n,prefix,nd);}

}

static void unpack_soa(unsigned char *s,unsigned n){
    for(unsigned i=0;i<n;i+=32){
        unsigned count=std::min(32u,n-i);unsigned char flat[96];
        const uint16_t *lo=(const uint16_t*)(s+3*i);const unsigned char *hi=s+3*i+64;
        for(unsigned j=0;j<count;j++)write24(flat+3*j,unsigned(lo[j])|(unsigned(hi[j])<<16));
        memcpy(s+3*i,flat,3*count);
    }
}
alignas(64) static uint16_t middle[1<<21];
static bool optimistic_unique;

namespace fastsort {
constexpr unsigned kRadix = 256, kRecordsPerBlock = 84, kBlockBytes = 256;
constexpr unsigned kLeafLimit = 8192;
constexpr int kSmallSort = 4096;
constexpr size_t kScratchOffset = 380'000'000, kWorkspaceBytes = 384'000'000;
using U32 = unsigned; using Byte = uint8_t; using U16 = uint16_t; using U64 = uint64_t;
using Vec256 = __m256i; using Vec128 = __m128i;
template <class T> inline T load(const void *src) { T value; std::memcpy(&value, src, sizeof value); return value; }
template <class F, size_t... I> [[gnu::always_inline]] inline void each(F &&f, std::index_sequence<I...>) {
  (f(std::integral_constant<size_t, I>{}), ...);
}
template <size_t N, class F> [[gnu::always_inline]] inline void repeat(F &&f) {
  each(std::forward<F>(f), std::make_index_sequence<N>{});
}
} // namespace fastsort

namespace byte_sort {
using namespace fastsort;
[[gnu::noinline]] inline void counting_sort(const U16 *src, U32 n, U32 *dst, U32 base) {
  alignas(64) U32 counts[65536] = {};
  for (U32 i = 0; i < n; ++i) ++counts[load<U16>(src + i)];
  for (U32 v = 0; v < 65536; ++v)
    for (U32 c = counts[v]; c; --c) *dst++ = base | v;
}
inline void expand8(Vec128 x, U32 *dst, U32 c, U32 base, U32 *end) {
  Vec256 y = _mm256_or_si256(_mm256_cvtepu8_epi32(x), _mm256_set1_epi32(base));
  if (end - dst >= 8) _mm256_storeu_si256((Vec256 *)dst, y);
  else _mm256_maskstore_epi32((int *)dst,
      _mm256_cmpgt_epi32(_mm256_set1_epi32(c), _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7)), y);
}
// Sort 32 independent columns, then transpose each into an 8-byte head and a 4-byte tail.
[[gnu::always_inline]] inline void sort12_columns(const Byte *src, Byte *out, Byte *tail) {
  Vec256 x[12], a[8], b[8];
  repeat<12>([&](auto j) __attribute__((always_inline)) {
    x[j] = _mm256_load_si256((const Vec256 *)(src + 256 * j));
  });
  // clang-format off
#define CMP(A, B) do { Vec256 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);
  // clang-format on
#undef CMP
#pragma GCC unroll 4
  for (U32 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 (U32 i = 0; i < 2; ++i) {
#pragma GCC unroll 2
    for (U32 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 (U32 i = 0; i < 4; ++i) {
    Vec256 lo = _mm256_unpacklo_epi32(b[i], b[i + 4]);
    Vec256 hi = _mm256_unpackhi_epi32(b[i], b[i + 4]);
    _mm256_store_si256((Vec256 *)(out + 32 * i), _mm256_permute2x128_si256(lo, hi, 0x20));
    _mm256_store_si256((Vec256 *)(out + 128 + 32 * i), _mm256_permute2x128_si256(lo, hi, 0x31));
  }
  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 (U32 i = 0; i < 2; ++i) {
    _mm256_store_si256((Vec256 *)(tail + 32 * i), _mm256_permute2x128_si256(b[2 * i], b[2 * i + 1], 0x20));
    _mm256_store_si256((Vec256 *)(tail + 64 + 32 * i), _mm256_permute2x128_si256(b[2 * i], b[2 * i + 1], 0x31));
  }
}
// Insert into a sorted 16-byte vector with at least one trailing 255 sentinel.
inline Vec128 insert_byte(Vec128 x, U32 byte) {
  return _mm_min_epu8(x, _mm_max_epu8(_mm_slli_si128(x, 1), _mm_set1_epi8(char(byte))));
}
// Encoded counter = bucket index + 256 * count. Recover counts and flag >12 / >32.
[[gnu::always_inline]] inline U32 restore_counts(U32 *count, U32 *exceptional) {
  U32 *end = count + 256; U32 over;
  const Vec256 limit12 = _mm256_set1_epi32(12), limit32 = _mm256_set1_epi32(32);
  __asm__ volatile("vpxor %%ymm4, %%ymm4, %%ymm4\n\t"
                   ".p2align 4\n\t1:\n\t"
                   ".irp reg,0,1,2,3; vmovdqa 32*\\reg(%[count]), %%ymm\\reg; .endr\n\t"
                   ".irp reg,0,1,2,3; vpsrld $8, %%ymm\\reg, %%ymm\\reg; .endr\n\t"
                   ".irp reg,0,1,2,3; vmovdqa %%ymm\\reg, 32*\\reg(%[count]); .endr\n\t"
                   "vpmaxud %%ymm1, %%ymm0, %%ymm0; vpmaxud %%ymm3, %%ymm2, %%ymm2\n\t"
                   "vpmaxud %%ymm2, %%ymm0, %%ymm0\n\t"
                   "vpcmpgtd %[limit12], %%ymm0, %%ymm1\n\t"
                   "vmovmskps %%ymm1, %%eax; movl %%eax, (%[flags])\n\t"
                   "vpmaxud %%ymm0, %%ymm4, %%ymm4\n\t"
                   "addq $128, %[count]; addq $4, %[flags]\n\t"
                   "cmpq %[end], %[count]; jb 1b\n\t"
                   "vpcmpgtd %[limit32], %%ymm4, %%ymm0; vmovmskps %%ymm0, %k[over]\n\t"
                   : [count] "+&r"(count), [flags] "+&r"(exceptional), [over] "=&r"(over)
                   : [end] "r"(end), [limit12] "x"(limit12), [limit32] "x"(limit32)
                   : "rax", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "cc", "memory");
  return over;
}
// All 32 counts must be <=12 and dst must have 384 writable elements.
// Overlapping 12-element stores are repaired by subsequent buckets in forward order.
[[gnu::always_inline]] inline U32 *emit32_buckets(const Byte *first, const Byte *second, const U32 *count,
                                              U32 *dst, Vec256 prefix) {
  const Vec256 step = _mm256_set1_epi32(256); const U32 *end = count + 32;
  __asm__ volatile(".p2align 4\n\t1:\n\t"
                   ".irp lane,0,1,2,3,4,5,6,7\n\t"
                   "movl 4*\\lane(%[count]), %%eax\n\t"
                   "vpmovzxbd 8*\\lane(%[first]), %%ymm1; vpmovzxbd 4*\\lane(%[second]), %%xmm2\n\t"
                   "vpor %[prefix], %%ymm1, %%ymm1; vpor %x[prefix], %%xmm2, %%xmm2\n\t"
                   "vmovdqu %%ymm1, (%[dst]); vmovdqu %%xmm2, 32(%[dst])\n\t"
                   "leaq (%[dst],%%rax,4), %[dst]\n\t"
                   "vpaddd %[step], %[prefix], %[prefix]\n\t"
                   ".endr\n\t"
                   "addq $64, %[first]; addq $32, %[second]; addq $32, %[count]\n\t"
                   "cmpq %[end], %[count]; jb 1b\n\t"
                   : [dst] "+&r"(dst), [prefix] "+&x"(prefix), [count] "+&r"(count), [first] "+&r"(first),
                     [second] "+&r"(second)
                   : [end] "r"(end), [step] "x"(step)
                   : "rax", "ymm1", "ymm2", "cc", "memory");
  return dst;
}
// end may extend past this leaf only when output cannot overwrite unread source.
inline void sort_leaf(const U16 *src, U32 n, U32 *dst, U32 base, U32 *end) {
  if (n < 64) {
    U32 copy[64];
    for (U32 i = 0; i < n; ++i) copy[i] = base | load<U16>(src + i);
    std::sort(copy, copy + n);
    std::memcpy(dst, copy, n * 4); return;
  }
  if (n > kLeafLimit) return counting_sort(src, n, dst, base);
  alignas(64) static Byte columns[kRadix * kLeafLimit];
  alignas(64) U32 count[256];
  // Only the first 12 rows need sentinels; later rows are read up to their real count.
  std::memset(columns, 255, 12 * kRadix);
  for (U32 i = 0; i < 256; ++i) count[i] = i;
  for (U32 i = 0; i < n; ++i) {
    U32 v = load<U16>(src + i), h = v >> 8, slot = count[h];
    count[h] = slot + 256; columns[slot] = Byte(v);
  }
  U32 exceptional[8];
  if (restore_counts(count, exceptional)) return counting_sort(src, n, dst, base);
  for (U32 block = 0; block < 256; block += 32) {
    alignas(32) Byte first[256], second[128];
    sort12_columns(columns + block, first, second);
    Vec256 prefix = _mm256_set1_epi32(base | (block << 8));
    if (exceptional[block >> 5] == 0 && end - dst >= 384) {
      dst = emit32_buckets(first, second, count + block, dst, prefix);
      continue;
    }
    for (U32 i = 0; i < 32; ++i) {
      U32 h = block + i, c = count[h], p = base | (h << 8);
      if (__builtin_expect(c <= 12, 1)) {
        Vec128 x = _mm_loadl_epi64((const Vec128 *)(first + 8 * i));
        U32 four = load<U32>(second + 4 * i);
        Vec256 wide = _mm256_or_si256(_mm256_cvtepu8_epi32(x), prefix);
        if (__builtin_expect(end - dst >= 12, 1)) {
          __asm__ volatile("vmovdqu {%1, %0|%0, %1}" : "=m"(*reinterpret_cast<__m256i_u *>(dst)) : "x"(wide));
          Vec128 y = _mm_or_si128(_mm_cvtepu8_epi32(_mm_cvtsi32_si128(four)), _mm256_castsi256_si128(prefix));
          _mm_storeu_si128((Vec128 *)(dst + 8), y);
        } else {
          alignas(32) U32 copy[16];
          _mm256_store_si256((Vec256 *)copy, wide);
          Vec128 y = _mm_or_si128(_mm_cvtepu8_epi32(_mm_cvtsi32_si128(four)), _mm256_castsi256_si128(prefix));
          _mm_store_si128((Vec128 *)(copy + 8), y);
          std::memcpy(dst, copy, c * 4);
        }
      } else if (c <= 16) {
        Vec128 x = _mm_loadl_epi64((const Vec128 *)(first + 8 * i));
        U32 four = load<U32>(second + 4 * i);
        x = _mm_unpacklo_epi64(x, _mm_cvtsi64_si128(uint64_t(four) | 0xffffffff00000000ull));
        for (U32 j = 12; j < c; ++j) x = insert_byte(x, columns[256 * j + h]);
        expand8(x, dst, 8, p, end);
        expand8(_mm_srli_si128(x, 8), dst + 8, c - 8, p, end);
      } else {
        Byte copy[32];
        for (U32 j = 0; j < c; ++j) copy[j] = columns[256 * j + h];
        std::sort(copy, copy + c);
        for (U32 j = 0; j < c; ++j) dst[j] = p | copy[j];
      }
      dst += c;
      prefix = _mm256_add_epi32(prefix, _mm256_set1_epi32(256));
    }
  }
}
} // namespace byte_sort
static void small(unsigned char *__restrict s,unsigned *__restrict d,unsigned n,unsigned prefix,unsigned *array_end){
    if(n<65536){small_original(s,d,n,prefix,array_end);return;}
    unsigned cap=(((n+255)/256*5/4+63)/64)*64+32;
    if(cap>4096){small_original(s,d,n,prefix,array_end);return;}
    unsigned mp[256];
    for(unsigned k=0;k<256;k++)mp[k]=k*cap;
    #pragma GCC unroll 12
    for(unsigned i=0;i<n;i++){
        const unsigned char *r=s+3*i;unsigned k=r[2];uint16_t value;memcpy(&value,r,2);
        middle[mp[k]++]=value;
    }
    for(unsigned k=0;k<256;k++)if(mp[k]-k*cap>cap){small_original(s,d,n,prefix,array_end);return;}
    for(unsigned k=0;k<256;k++){
        unsigned count=mp[k]-k*cap,pre=prefix|(k<<16);
        const uint16_t *src=middle+k*cap;
        byte_sort::sort_leaf(src,count,d,pre,array_end);
        d+=count;
    }
}

template<bool Reserved> static bool partition(unsigned *__restrict a,unsigned n,unsigned capacity=0){
    for(unsigned k=0;k<K;k++)byte_cursor[k]=k*320;
    for(unsigned page=0;page<n;page+=1024){
      if(page+1024<n)asm volatile("orl $0, %0" : "+m"(a[page+1024]));
      unsigned stop=std::min(n,page+1024);
      #pragma GCC unroll 8
      for(unsigned i=page;i<stop;i++){
        unsigned x=a[i],k=x>>24,o=byte_cursor[k]+3;
        memcpy(staging_bytes+o-3,&x,4);
        if(__builtin_expect((o&63)==0,0)){
            if(Reserved && pos[k]+64-starts[k]>capacity){_mm_sfence();return false;}
            unsigned char *d=workspace+3*pos[k];
            const unsigned char *src=staging_bytes+o-192;
            #pragma GCC unroll 6
            for(unsigned j=0;j<192;j+=32)
                _mm256_stream_si256((__m256i*)(d+j),_mm256_load_si256((const __m256i*)(src+j)));
            pos[k]+=64;o-=192;
        }
        byte_cursor[k]=o;
      }
    }
    for(unsigned k=0;k<K;k++)pos[k]+=(byte_cursor[k]-k*320)/3;
    _mm_sfence();
    if(Reserved)for(unsigned k=0;k<K;k++)if(pos[k]-starts[k]>capacity)return false;
    return true;
}

void sort(unsigned *a,int n){
    optimistic_unique=true;
    if(n<256){std::sort(a,a+n);return;}
    unsigned capacity=((unsigned(n)+K-1)/K*17/16+BS-1)&~(BS-1);
    for(unsigned k=0;k<K;k++)starts[k]=pos[k]=k*capacity;
    if(!partition<true>(a,n,capacity)){
        memset(hist,0,sizeof(hist));
        for(int i=0;i<n;i++)++hist[a[i]>>(32-TB)];
        unsigned z=0;
        for(unsigned k=0;k<K;k++){starts[k]=pos[k]=z;z=(z+hist[k]+BS-1)&~(BS-1);}
        partition<false>(a,n);
    }
    unsigned z=0;
    for(unsigned k=0;k<K;k++){
        unsigned size=pos[k]-starts[k];
        memcpy(workspace+3*(pos[k]&~(BS-1)),staging_bytes+k*320,3*(pos[k]&(BS-1)));
        small(workspace+3*starts[k],a+z,size,k<<24,a+n);
        z+=size;
    }
}

评测结果