代码

#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#include <bits/stdc++.h>
#include <immintrin.h>

// V9: high-byte partition, next-high-byte partition, exact low-16-bit sorting.
// Global partition and its recovery checks are retained from the V4 baseline.
// The low-16-bit path uses two 8 KiB bitmaps, not four-pass LSD scattering.
// No timing code, background work, threads, I/O or input-distribution assumptions.
// Requires GCC-compatible x86-64 with the ISA options named above.
namespace duck_sort_v9_detail {
template <int STRIDE>
std::array<int, 256> get_population_upper_bounds(uint8_t* A, int N, int budget, int sample_size) {
    std::array<int, 256> results;
    results.fill(0);

    size_t n = (size_t)sample_size;
    
    // 1. Calculate Cache Line Alignment Info
    uintptr_t start_addr = (uintptr_t)A;
    // Align up to next 64-byte boundary
    uintptr_t aligned_start = (start_addr + 63) & ~63ULL;
    // Align down end address
    uintptr_t end_addr_exclusive = start_addr + (size_t)N * STRIDE;
    uintptr_t aligned_end = end_addr_exclusive & ~63ULL;

    if (aligned_end <= aligned_start) {
        // Not enough data for aligned sampling, fallback to full scan
         for (int i = 0; i < N; ++i) results[A[(size_t)i * STRIDE]]++;
         return results;
    }

    size_t num_lines = (aligned_end - aligned_start) / 64;
    
    // Safety check: if no cache lines available, fallback to full scan
    if (num_lines == 0) {
        for (int i = 0; i < N; ++i) results[A[(size_t)i * STRIDE]]++;
        return results;
    }
    
    // Determine the offset pattern for the FIRST aligned block
    // We need (aligned_start + off) % STRIDE == start_addr % STRIDE
    size_t diff = aligned_start - start_addr;
    int base_offset = (STRIDE - (diff % STRIDE)) % STRIDE;

    // 2. Adjust Sample Size to be in terms of cache lines
    // Average items per line
    int items_per_line_approx = 64 / STRIDE;
    size_t lines_to_sample = (n + items_per_line_approx - 1) / items_per_line_approx;
    
    // Cap at available lines
    if (lines_to_sample > num_lines) lines_to_sample = num_lines;
    
    // Recalculate actual n for statistics
    // (This is an approximation if stride=3 because different lines have different counts, 
    // but for large N it converges)
    // For Stride=4, count is always 16.
    // For Stride=3, count is 21 or 22 (avg 21.33).
    // optimizing: just counting actually sampled items is better, 
    // but user code expects 'n' to be passed to math formulas.
    // We will count exact sampled items in the loop.
    
    // 3. Sparse Sampling of Cache Lines
    std::array<int, 256> sample_counts;
    sample_counts.fill(0);
    size_t actual_sampled_count = 0;
    
    static std::mt19937 gen;
    const uint32_t mod_blocks = (uint32_t)num_lines;
    const uint64_t mu = ((unsigned __int128)1 << 64) / mod_blocks;

    for (size_t i = 0; i < lines_to_sample; ++i) {
        // Random Block Index
        uint32_t x = gen();
        uint64_t q = ((unsigned __int128)x * mu) >> 64;
        uint32_t blk_idx = x - q * mod_blocks;
        if (blk_idx >= mod_blocks) blk_idx -= mod_blocks;

        uint8_t* p_line = (uint8_t*)(aligned_start + (size_t)blk_idx * 64);

        // Calculate offset for this specific block
        // Block addr changes by 64. 64 % 3 = 1. 64 % 4 = 0.
        // offset_new = (offset_old - delta_addr) % STRIDE
        // delta_addr = blk_idx * 64
        int current_offset;
        if constexpr (STRIDE == 4) {
             current_offset = base_offset;
        } else {
             // STRIDE == 3
             // shift = (blk_idx) % 3
             // off = (base - shift) % 3
             int shift = blk_idx % 3;
             current_offset = base_offset - shift;
             if (current_offset < 0) current_offset += 3;
        }

        // Fetch fixed number of items per cache line
        // Safe max index check:
        // Stride 4: offset max 3. count 16. max idx = 3 + 15*4 = 63 < 64.
        // Stride 3: offset max 2. count 21. max idx = 2 + 20*3 = 62 < 64.
        const int ITEMS = 64 / STRIDE;
        
        #pragma GCC unroll 21
        for (int k = 0; k < ITEMS; ++k) {
            sample_counts[p_line[current_offset + k * STRIDE]]++;
        }
        actual_sampled_count += ITEMS;
    }
    
    n = actual_sampled_count;

    // 2. 二分查找最优 Z 值
    // 目标：找到最大的 Z，使得 Sum(UpperBounds(Z)) <= Budget
    double low_z = 0.0;
    double high_z = 10.0;
    double best_z = 0.0;
    
    double n_double = (double)n;
    double N_double = (double)N;
    double fpc = (double)(N - n) / (double)(N - 1);
    if (fpc < 0) fpc = 0; // Safety

    // 预计算 p_hat 以加速循环
    std::array<double, 256> p_hats;
    for(int i=0; i<256; ++i) p_hats[i] = sample_counts[i] / n_double;

    for (int iter = 0; iter < 20; ++iter) {
        double mid_z = (low_z + high_z) * 0.5;
        double z2 = mid_z * mid_z;
        double div_factor = 1.0 / (1.0 + z2 / n_double);
        
        long long current_sum = 0;
        
        for (int i = 0; i < 256; ++i) {
            double p_hat = p_hats[i];
            
            // Wilson Score Interval
            double term1 = p_hat + z2 / (2.0 * n_double);
            double variance_term = (p_hat * (1.0 - p_hat) / n_double) * fpc;
            if (variance_term < 0) variance_term = 0;
            double term2 = mid_z * std::sqrt(variance_term + z2 / (4.0 * n_double * n_double));
            double p_upper = (term1 + term2) * div_factor;
            
            int limit = (int)std::ceil(N_double * p_upper);
            current_sum += limit;
        }

        if (current_sum <= budget) {
            best_z = mid_z;
            low_z = mid_z;
        } else {
            high_z = mid_z;
        }
    }

    // 3. 使用最佳 Z 生成最终结果
    double z = best_z;
    double z2 = z * z;
    double div_factor = 1.0 / (1.0 + z2 / n_double);
    
    for (int i = 0; i < 256; ++i) {
        double p_hat = p_hats[i];
        double term1 = p_hat + z2 / (2.0 * n_double);
        double variance_term = (p_hat * (1.0 - p_hat) / n_double) * fpc;
        if (variance_term < 0) variance_term = 0;
        double term2 = z * std::sqrt(variance_term + z2 / (4.0 * n_double * n_double));
        double p_upper = (term1 + term2) * div_factor;
        
        int limit = (int)std::ceil(N_double * p_upper);
        if (limit > N) limit = N;
        results[i] = limit;
    }

    return results;
}

using namespace std;

const int n = 1e8;
const int PREFETCH_DIST = 64; // 元素个数：Pass1(256B), Pass2(192B), Pass3/4(128B)

// 辅助函数：向地址 p 写入 3 字节 (利用 uint32 覆盖写，需保证 buffer 有 padding)
// Input val: [B0, B1, B2, X] (Little Endian) -> Writes B0, B1, B2
inline void store3(uint8_t* __restrict__ p, uint32_t val) {
    std::memcpy(p, &val, 4);
}

// 辅助函数：向地址 p 写入 2 字节
inline void store2(uint8_t* __restrict__ p, uint16_t val) {
    std::memcpy(p, &val, 2);
}


inline uint32_t load3(const uint8_t* p) {
    uint32_t v; std::memcpy(&v,p,4); return v;
}
inline uint16_t load2(const uint8_t* p) {
    uint16_t v; std::memcpy(&v,p,2); return v;
}
// Verify that a full tile can be scattered without leaving the allocation.
// Bucket overlap is detected after the pass; the original source stays intact.
inline bool top_tile_fits(uint8_t* const* p, const uint8_t* limit) {
    const __m256i bound=_mm256_set1_epi64x(reinterpret_cast<intptr_t>(limit));
    __m256i bad=_mm256_setzero_si256();
    for (int k=0;k<256;k+=4) {
        const __m256i q=_mm256_loadu_si256(reinterpret_cast<const __m256i*>(p+k));
        bad=_mm256_or_si256(bad,_mm256_cmpgt_epi64(q,bound));
    }
    return _mm256_movemask_epi8(bad)==0;
}


// Sort a uint16_t multiset. Two bit planes encode multiplicities 1 and 2.
// Third and later copies are recorded separately; many repetitions use counting.
struct Low16Workspace {
    alignas(64) uint64_t seen[1024];
    alignas(64) uint64_t twice[1024];
    unsigned extra[65];
    unsigned* dense = nullptr;
    ~Low16Workspace(){std::free(dense);}
};
__attribute__((always_inline)) inline bool insert_bit(uint64_t& word, unsigned bit) {
    bool present;
    uint64_t x=word;
    __asm__("btsq %2, %0" : "+r"(x), "=@ccc"(present) : "r"(uint64_t(bit)));
    word=x;
    return present;
}
inline void low16_dense(const uint8_t* src, int c, unsigned* dst,
                        unsigned prefix, Low16Workspace& ws) {
    if(!ws.dense) ws.dense=static_cast<unsigned*>(std::malloc(65536*sizeof(unsigned)));
    if(!ws.dense){
        for(int i=0;i<c;++i)dst[i]=prefix|load2(src+2*size_t(i));
        std::sort(dst,dst+c);return;
    }
    std::memset(ws.dense,0,65536*sizeof(unsigned));
    for(int i=0;i<c;++i)++ws.dense[load2(src+2*size_t(i))];
    for(unsigned x=0;x<65536;++x){
        unsigned num=ws.dense[x];
        for(unsigned t=0;t<num;++t)*dst++=prefix|x;
    }
}
__attribute__((noinline))
void low16_bitmap(const uint8_t* src, int c, unsigned* dst,
                  unsigned prefix, Low16Workspace& ws) {
    if(c<32){
        for(int i=0;i<c;++i)dst[i]=prefix|load2(src+2*size_t(i));
        std::sort(dst,dst+c);return;
    }
    if(c>8192){low16_dense(src,c,dst,prefix,ws);return;}
    std::memset(ws.seen,0,sizeof(ws.seen));
    std::memset(ws.twice,0,sizeof(ws.twice));
    unsigned used=0;
    for(int i=0;i<c;++i) {
        unsigned v=load2(src+2*size_t(i));
        if(__builtin_expect(insert_bit(ws.seen[v>>6],v),0)) {
            if(__builtin_expect(insert_bit(ws.twice[v>>6],v),0)) {
                if(used==64){low16_dense(src,c,dst,prefix,ws);return;}
                ws.extra[used++]=v;
            }
        }
    }
    if(used>1)std::sort(ws.extra,ws.extra+used);
    ws.extra[used]=65536;
    unsigned ex=0;
    unsigned* end=dst+c;
    for(unsigned word=0;word<1024;++word) {
        uint64_t bits=ws.seen[word];
        unsigned base=prefix|(word<<6);
        if(!ws.twice[word] && end-dst>=4) {
            unsigned num=__builtin_popcountll(bits);
            unsigned v0=_tzcnt_u64(bits);bits=_blsr_u64(bits);
            unsigned v1=_tzcnt_u64(bits);bits=_blsr_u64(bits);
            unsigned v2=_tzcnt_u64(bits);bits=_blsr_u64(bits);
            unsigned v3=_tzcnt_u64(bits);bits=_blsr_u64(bits);
            __m128i values=_mm_setr_epi32(v0,v1,v2,v3);
            values=_mm_add_epi32(values,_mm_set1_epi32(base));
            _mm_storeu_si128(reinterpret_cast<__m128i*>(dst),values);
            dst+=std::min(num,4u);
            while(bits){unsigned x=_tzcnt_u64(bits);bits=_blsr_u64(bits);*dst++=base|x;}
            continue;
        }
        if(!ws.twice[word]){
            while(bits){unsigned x=_tzcnt_u64(bits);bits=_blsr_u64(bits);*dst++=base|x;}
        }else{
            uint64_t twice=ws.twice[word];
            while(bits) {
                unsigned x=_tzcnt_u64(bits);bits=_blsr_u64(bits);
                unsigned val=base|x;
                *dst++=val;
                if((twice>>x)&1){
                    *dst++=val;
                    unsigned low=(word<<6)|x;
                    while(ws.extra[ex]==low){*dst++=val;++ex;}
                }
            }
        }
    }
}

template<int FixedN>
void sort_impl(uint* a, int __n) {
    const int n=FixedN ? FixedN : __n;
    if(n<=1)return;
    if(n<4096){std::sort(a,a+n);return;}


    // ---------------------------------------------------------
    // Pass 1: Global MSD (Partition by B3)
    // Read: a (4 bytes) -> Write: b (3 bytes: [B0, B1, B2])
    // ---------------------------------------------------------
    
    uint cnt_global[256];
    // memset(cnt_global, 0, sizeof(cnt_global)); // No longer needed beforehand
    
    // 1.1 统计 B3 (Sampling & Upper Bounds)
    // Budget set to n * 1.47 (47% over-provisioning)
    int budget = (int)(n * 1.47);
    int sample_size = 20000;
    // A 是 (uint8_t*)a + 3 (B3 byte), stride = 4
    auto bounds = get_population_upper_bounds<4>((uint8_t*)a + 3, n, budget, sample_size);
    
    // 1.2 计算 B3 Offset (Bytes in b)
    // 增加 4 字节 Padding 以安全使用 store3
    uint ptr_global[256];
    uint32_t offset_b3 = 0;
    for (int i = 0; i < 256; i++) {
        ptr_global[i] = offset_b3;
        offset_b3 += bounds[i] * 3; // Use Upper Bound
    }
    
    // 申请 b 数组 
    constexpr int TILE=16384;
    const size_t main_bytes=size_t(budget)*3;
    uint8_t* b=static_cast<uint8_t*>(std::malloc(main_bytes+3*TILE+4096));
    if(!b){std::sort(a,a+n);return;}
    bool incomplete=false;
    // 1.3 执行 Pass 1 分发
    {
        uint* __restrict__ src = a;
        uint8_t* __restrict__ dst = b;
        uint p[256];
        uint8_t* pp[256];
        for(int z=0;z<256;++z)pp[z]=b+ptr_global[z];

        int i=0;
        for(;i<n;) {
            if (__builtin_expect(!top_tile_fits(pp,b+main_bytes),0)) {
                incomplete=true;break;
            }
            const int end=std::min(n,i+TILE);
            for(;i+16<=end;i+=16) {
                _mm_prefetch(reinterpret_cast<const char*>(
                    reinterpret_cast<uintptr_t>(src)+size_t(i+PREFETCH_DIST)*4),_MM_HINT_NTA);
                #pragma GCC unroll 16
                for(int j=0;j<16;++j) {
                    const unsigned v=src[i+j], k=v>>24;
                    uint8_t* q=pp[k];
                    store3(q,v); pp[k]=q+3;
                }
            }
            for(;i<end;++i){const unsigned v=src[i],k=v>>24;store3(pp[k],v);pp[k]+=3;}
        }
        
        // Reconstruct exact counts from pointer progress
        for(int k=0; k<256; ++k) {
            cnt_global[k] = (pp[k] - (b+ptr_global[k])) / 3;
        }
    }
    bool retry_global=incomplete;
    if(incomplete) {
        std::memset(cnt_global,0,sizeof(cnt_global));
        for(int i=0;i<n;++i)++cnt_global[a[i]>>24];
    }
    for(int k=0;k<256;++k)
        retry_global |= cnt_global[k] && cnt_global[k]>=static_cast<unsigned>(bounds[k]);
    if(retry_global) {
        unsigned off=0;
        uint8_t* pp[256];
        for(int k=0;k<256;++k) {
            ptr_global[k]=off;pp[k]=b+off;
            off+=3*cnt_global[k]+4; // Allow the fourth byte of store3.
        }
        int i=0;
        for(;i+16<=n;i+=16) {
            #pragma GCC unroll 16
            for(int j=0;j<16;++j) {
                unsigned v=a[i+j],k=v>>24;store3(pp[k],v);pp[k]+=3;
            }
        }
        for(;i<n;++i){unsigned v=a[i],k=v>>24;store3(pp[k],v);pp[k]+=3;}
    }

    uint8_t* scratch=nullptr;
    size_t scratch_capacity=0;
    Low16Workspace ws;
    unsigned emitted=0;
    for(unsigned h=0;h<256;++h) {
        unsigned count=cnt_global[h];
        if(!count)continue;
        uint8_t* src=b+ptr_global[h];
        unsigned* dst=a+emitted;
        if(count<1024){
            for(unsigned j=0;j<count;++j)dst[j]=(h<<24)|(load3(src+size_t(j)*3)&0xffffffu);
            std::sort(dst,dst+count);emitted+=count;continue;
        }
        auto cap=get_population_upper_bounds<3>(src+2,int(count),int(count*2),5000);
        unsigned start[256],cnt[256];
        size_t bytes=0;
        for(unsigned k=0;k<256;++k){start[k]=unsigned(bytes);bytes+=size_t(cap[k])*2;}
        size_t need=std::max(bytes+size_t(count)*2+8,size_t(count)*2+1032);
        uint8_t* temp;
        // The whole current output range and temp are disjoint. Unlike the
        // old four-pass variant, output is now final directly after Pass 2.
        size_t free_bytes=size_t(n-emitted-count)*4;
        if(need+64<=free_bytes){
            uintptr_t end_addr=reinterpret_cast<uintptr_t>(a+n);
            temp=reinterpret_cast<uint8_t*>((end_addr-need)&~uintptr_t(63));
        }else{
            if(need>scratch_capacity){
                uint8_t* q=static_cast<uint8_t*>(std::malloc(need));
                if(!q){
                    unsigned t=emitted;
                    for(unsigned k=h;k<256;++k)
                        for(unsigned j=0;j<cnt_global[k];++j)
                            a[t++]=(k<<24)|(load3(b+ptr_global[k]+size_t(j)*3)&0xffffffu);
                    std::sort(a+emitted,a+n);
                    std::free(scratch);std::free(b);return;
                }
                std::free(scratch);scratch=q;scratch_capacity=need;
            }
            temp=scratch;
        }
        uint8_t* p[256];
        for(unsigned k=0;k<256;++k)p[k]=temp+start[k];
        unsigned i=0;
        for(;i+16<=count;i+=16){
            _mm_prefetch(reinterpret_cast<const char*>(reinterpret_cast<uintptr_t>(src)+size_t(i+64)*3),_MM_HINT_T0);
            #pragma GCC unroll 16
            for(unsigned j=0;j<16;++j){
                unsigned v=load3(src+size_t(i+j)*3), key=(v>>16)&255;
                store2(p[key],uint16_t(v));p[key]+=2;
            }
        }
        for(;i<count;++i){unsigned v=load3(src+size_t(i)*3),key=(v>>16)&255;store2(p[key],uint16_t(v));p[key]+=2;}
        bool retry=false;
        for(unsigned k=0;k<256;++k){cnt[k]=unsigned(p[k]-(temp+start[k]))/2;retry|=cnt[k]>unsigned(cap[k]);}
        if(retry){
            unsigned off=0;
            for(unsigned k=0;k<256;++k){start[k]=off;p[k]=temp+off;off+=cnt[k]*2;}
            for(unsigned i=0;i<count;++i){unsigned v=load3(src+size_t(i)*3),key=(v>>16)&255;store2(p[key],uint16_t(v));p[key]+=2;}
        }
        for(unsigned k=0;k<256;++k){
            if(cnt[k])low16_bitmap(temp+start[k],int(cnt[k]),dst,(h<<24)|(k<<16),ws);
            dst+=cnt[k];
        }
        emitted+=count;
    }
    std::free(scratch);std::free(b);
}
} // namespace duck_sort_v9_detail
void sort(unsigned* a,int n) {
    static_assert(sizeof(unsigned)==4,"32-bit unsigned required");
    if(n==100000000)duck_sort_v9_detail::sort_impl<100000000>(a,n);
    else duck_sort_v9_detail::sort_impl<0>(a,n);
}

评测结果