代码

// C++17 / x86-64 AVX2 + BMI2. Single thread; final output is in a.
// V6: keep the V5 global compressed partition. In local buckets with >=65536
// elements, replace the low-byte sampler by checked speculative bucket slots,
// and recover B1 counts from third-pass write cursors instead of counting B1
// for every element in pass 2. B2 is still counted in pass 2.
// Slot pitch is an odd multiple of 64 bytes. Speculation affects performance,
// not correctness: original local input is retained for exact retry.
// Auxiliary allocation/layout are otherwise inherited from V5.
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#include <bits/stdc++.h>
#include <immintrin.h>

#ifndef DUCK_PREFETCH_BYTES
#define DUCK_PREFETCH_BYTES 32
#endif

namespace duck_sort_v6_detail {

__attribute__((always_inline)) inline void prefetch_destination(const void* p) {
    _mm_prefetch(reinterpret_cast<const char*>(
        reinterpret_cast<std::uintptr_t>(p) + DUCK_PREFETCH_BYTES), _MM_HINT_T1);
}

// Original sampled capacity allocator; used globally and for small local buckets.
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;

    uintptr_t start_addr = (uintptr_t)A;
    
    uintptr_t aligned_start = (start_addr + 63) & ~63ULL;
    
    uintptr_t end_addr_exclusive = start_addr + (size_t)N * STRIDE;
    uintptr_t aligned_end = end_addr_exclusive & ~63ULL;

    if (aligned_end <= aligned_start) {
        
         for (int i = 0; i < N; ++i) results[A[(size_t)i * STRIDE]]++;
         return results;
    }

    size_t num_lines = (aligned_end - aligned_start) / 64;

    if (num_lines == 0) {
        for (int i = 0; i < N; ++i) results[A[(size_t)i * STRIDE]]++;
        return results;
    }

    size_t diff = aligned_start - start_addr;
    int base_offset = (STRIDE - (diff % STRIDE)) % STRIDE;

    int items_per_line_approx = 64 / STRIDE;
    size_t lines_to_sample = (n + items_per_line_approx - 1) / items_per_line_approx;

    if (lines_to_sample > num_lines) lines_to_sample = num_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) {
        
        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);

        int current_offset;
        if constexpr (STRIDE == 4) {
             current_offset = base_offset;
        } else {

             int shift = blk_idx % 3;
             current_offset = base_offset - shift;
             if (current_offset < 0) current_offset += 3;
        }

        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;

    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; 

    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];

            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;
        }
    }

    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 PREFETCH_DIST = 64; // Input elements; destination distance is in bytes. 

inline void store3(uint8_t* __restrict__ p, uint32_t val) {
    std::memcpy(p, &val, 4);
}

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;
}

// Caller reserves space for every write in the next tile/group.
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;
}

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: B3 partition; preserve a until the partition is complete.
    uint cnt_global[256];

    int budget = (int)(n * 1.47);
    int sample_size = 20000;
    
    auto bounds = get_population_upper_bounds<4>((uint8_t*)a + 3, n, budget, sample_size);

    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; 
    }

    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;
    
    {
        uint* __restrict__ src = a;

        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];
                    prefetch_destination(q);
                    store3(q,v); pp[k]=q+3;
                }
            }
            for(;i<end;++i){const unsigned v=src[i],k=v>>24;prefetch_destination(pp[k]);store3(pp[k],v);pp[k]+=3;}
        }

        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; 
        }
        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;prefetch_destination(pp[k]);store3(pp[k],v);pp[k]+=3;
            }
        }
        for(;i<n;++i){unsigned v=a[i],k=v>>24;prefetch_destination(pp[k]);store3(pp[k],v);pp[k]+=3;}
    }
    
    uint8_t* scratch=nullptr;
    size_t scratch_capacity=0;

    uint8_t* a_u8 = (uint8_t*)a;

    uint cnt0[256];
    uint cnt1[256];
    uint cnt2[256];
    uint ptr0[256]; 
    uint ptr1[256]; 
    uint ptr2[256]; 

    uint32_t a_offset_start = 0; 

    for (int i_b3 = 0; i_b3 < 256; i_b3++) {
        
        int count = cnt_global[i_b3];
        if (count == 0) continue;

        uint8_t* seg_b_in = b + ptr_global[i_b3];
        if(count<1024) {
            uint* out=a+a_offset_start;
            for(int j=0;j<count;++j)out[j]=(unsigned(i_b3)<<24)|(load3(seg_b_in+j*3)&0xffffffu);
            std::sort(out,out+count);a_offset_start+=count;continue;
        }

        uint8_t* seg_a_temp = a_u8 + (a_offset_start * 4);

        int budget_pass2 = count * 2;
        int sample_size = 5000;

        // Low-byte estimated slots. Exact counts are recovered after scatter.
        std::array<int,256> bounds;
        if(count>=65536) {
            unsigned need=(5ull*count+1023)/1024;
            unsigned cap=(((need+31)/32)|1u)*32;
            bounds.fill(cap);
        } else bounds=get_population_upper_bounds<3>(seg_b_in,count,budget_pass2,sample_size);

        uint32_t tmp = 0;
        for(int k=0; k<256; k++) { 
            ptr0[k] = tmp; 
            tmp += bounds[k] * 2;
        }

        const size_t need=std::max(size_t(tmp)+size_t(count)*2+8,
                                   size_t(count)*2+1032);
        if(need>size_t(n-a_offset_start)*4) {
            if(need>scratch_capacity) {
                uint8_t* next=static_cast<uint8_t*>(std::malloc(need));
                if(!next) {

                    unsigned t=a_offset_start;
                    for(int h=i_b3;h<256;++h)
                        for(unsigned j=0;j<cnt_global[h];++j)
                            a[t++]=(unsigned(h)<<24)|(load3(b+ptr_global[h]+3*j)&0xffffffu);
                    std::sort(a+a_offset_start,a+n);
                    std::free(scratch);std::free(b);return;
                }
                std::free(scratch);scratch=next;scratch_capacity=need;
            }
            seg_a_temp=scratch;
        }

        const bool lazy1=count>=65536;
        // Pass 2: B0 -> [B1,B2]. Fast path counts only B2.
        
        memset(cnt1, 0, sizeof(cnt1));
        memset(cnt2, 0, sizeof(cnt2));

        {

            uint8_t* pp[256];
            for(int z=0;z<256;++z) pp[z]=seg_a_temp+ptr0[z];
            uint8_t* src = seg_b_in;

            if(lazy1) {
            int k = 0;
            for (; k <= count - 16; k += 16) {
                _mm_prefetch(reinterpret_cast<const char*>(reinterpret_cast<uintptr_t>(src)+size_t(k+PREFETCH_DIST)*3), _MM_HINT_T0);
                #pragma GCC unroll 16
                for (int j = 0; j < 16; j++) {
                    uint32_t val = load3(src + (k + j) * 3);
                    uint8_t key = val & 0xFF; 

                    cnt2[(val >> 16) & 0xFF]++;

                    prefetch_destination(pp[key]);
                    store2(pp[key],val>>8);pp[key]+=2;
                }
            }
            for (; k < count; k++) {
                uint32_t val = load3(src + k * 3);
                uint8_t key = val & 0xFF;

                cnt2[(val >> 16) & 0xFF]++;
                
                prefetch_destination(pp[key]);
                    store2(pp[key],val>>8);pp[key]+=2;
            }

            } else {
            int k = 0;
            for (; k <= count - 16; k += 16) {
                _mm_prefetch(reinterpret_cast<const char*>(reinterpret_cast<uintptr_t>(src)+size_t(k+PREFETCH_DIST)*3), _MM_HINT_T0);
                #pragma GCC unroll 16
                for (int j = 0; j < 16; j++) {
                    uint32_t val = load3(src + (k + j) * 3);
                    uint8_t key = val & 0xFF; 

                    cnt1[(val >> 8) & 0xFF]++;
                    cnt2[(val >> 16) & 0xFF]++;

                    prefetch_destination(pp[key]);
                    store2(pp[key],val>>8);pp[key]+=2;
                }
            }
            for (; k < count; k++) {
                uint32_t val = load3(src + k * 3);
                uint8_t key = val & 0xFF;
                
                cnt1[(val >> 8) & 0xFF]++;
                cnt2[(val >> 16) & 0xFF]++;
                
                prefetch_destination(pp[key]);
                    store2(pp[key],val>>8);pp[key]+=2;
            }

            }
            
            bool retry = false;
            uint32_t tmp1 = 0;
            uint32_t tmp2 = 0;
            for(int k=0; k<256; k++) {
                
                cnt0[k] = (pp[k] - (seg_a_temp+ptr0[k])) >> 1; 
                retry |= cnt0[k] > bounds[k];
                
                ptr1[k] = tmp1; 
                tmp1 += cnt1[k] * 2 + 4; 
                
                ptr2[k] = tmp2; 
                tmp2 += cnt2[k];         
            }
            if (retry) {
                
                uint32_t tmp = 0;
                for(int k=0; k<256; k++) {
                    ptr0[k] = tmp;
                    tmp += cnt0[k] * 2 + 4; 
                }

                uint p_retry[256];
                memcpy(p_retry, ptr0, sizeof(p_retry));
                uint8_t* src = seg_b_in;

                int k = 0;
                for (; k <= count - 16; k += 16) {
                    _mm_prefetch(reinterpret_cast<const char*>(reinterpret_cast<uintptr_t>(src)+size_t(k+PREFETCH_DIST)*3), _MM_HINT_T0);
                    #pragma GCC unroll 16
                    for (int j = 0; j < 16; j++) {
                        uint32_t val = load3(src + (k + j) * 3);
                        uint8_t key = val & 0xFF; 
                        prefetch_destination(seg_a_temp + p_retry[key]);
                        store2(seg_a_temp + p_retry[key], val >> 8);
                        p_retry[key] += 2;
                    }
                }
                for (; k < count; k++) {
                    uint32_t val = load3(src + k * 3);
                    uint8_t key = val & 0xFF;
                    prefetch_destination(seg_a_temp + p_retry[key]);
                        store2(seg_a_temp + p_retry[key], val >> 8);
                    p_retry[key] += 2;
                }
            }
        }

        {
            
            uint8_t* dst_base=seg_b_in;
            // Pass 3: B1 -> [B0,B2], traversing B0 buckets in order.
            if(lazy1) {
                const unsigned cap=((((5ull*count+1023)/1024+31)/32)|1u)*32;
                for(unsigned h=0;h<256;++h)ptr1[h]=2*cap*h;
            }
            uint8_t* pp[256];
            bool incomplete=false;
            for(int h=0;h<256;++h)pp[h]=dst_base+ptr1[h];
            for(unsigned b0=0;b0<256;++b0) {
                const int c=cnt0[b0];
                if(lazy1 && !top_tile_fits(pp,dst_base+3*count-2*c)) {
                    incomplete=true; break;
                }
                const uint8_t* src=seg_a_temp+ptr0[b0];
                int k=0;
                for(;k+21<=c;k+=21) {
                    _mm_prefetch((const char*)(uintptr_t(src)+(k+64)*2),_MM_HINT_T0);
                    #pragma GCC unroll 21
                    for(int j=0;j<21;++j) {
                        unsigned v=load2(src+(k+j)*2),key=v&255;
                        prefetch_destination(pp[key]);
                        store2(pp[key],b0|(v&0xff00));pp[key]+=2;
                    }
                }
                for(;k<c;++k) {
                    unsigned v=load2(src+k*2),key=v&255;
                    prefetch_destination(pp[key]);store2(pp[key],b0|(v&0xff00));pp[key]+=2;
                }
            }
            if(lazy1) {
                bool retry=incomplete;
                if(incomplete) {
                    std::memset(cnt1,0,sizeof(cnt1));
                    for(int h=0;h<256;++h) {
                        const uint8_t* src=seg_a_temp+ptr0[h];
                        for(unsigned j=0;j<cnt0[h];++j)++cnt1[load2(src+2*j)&255];
                    }
                } else {
                    for(int h=0;h<256;++h) {
                        cnt1[h]=(pp[h]-(dst_base+ptr1[h]))/2;
                        if(h<255 && pp[h]>dst_base+ptr1[h+1])retry=true;
                    }
                }
                if(retry) {
                    unsigned offset=0;
                    for(int h=0;h<256;++h){ptr1[h]=offset;pp[h]=dst_base+offset;offset+=cnt1[h]*2;}
                    for(unsigned b0=0;b0<256;++b0){
                        const uint8_t* src=seg_a_temp+ptr0[b0];
                        const int c=cnt0[b0];
                        int k=0;
                        for(;k+21<=c;k+=21) {
                            #pragma GCC unroll 21
                            for(int j=0;j<21;++j) {
                                unsigned v=load2(src+(k+j)*2),key=v&255;
                                prefetch_destination(pp[key]);store2(pp[key],b0|(v&0xff00));pp[key]+=2;
                            }
                        }
                        for(;k<c;++k) {
                            unsigned v=load2(src+k*2),key=v&255;
                            prefetch_destination(pp[key]);store2(pp[key],b0|(v&0xff00));pp[key]+=2;
                        }
                    }
                }
            }
        }

        {
            
            // Pass 4: B2 -> full unsigned, traversing B1 buckets in order.
            uint* dst_base = a + a_offset_start;
            uint8_t* src_base = seg_b_in; 

            unsigned* pp[256];
            for(int z=0;z<256;++z)pp[z]=dst_base+ptr2[z];

            uint32_t val_b3_shifted = i_b3 << 24;

            for (int b1 = 0; b1 < 256; b1++) {
                int c = cnt1[b1];
                if (c == 0) continue;

                uint8_t* src = src_base + ptr1[b1];

                uint32_t common_bits = val_b3_shifted | (b1 << 8);

                int k = 0;
                for (; k <= c - 32; k += 32) {
                    _mm_prefetch(reinterpret_cast<const char*>(reinterpret_cast<uintptr_t>(src)+size_t(k+PREFETCH_DIST)*2), _MM_HINT_T0);
                    #pragma GCC unroll 32
                    for (int j = 0; j < 32; j++) {
                        
                        uint32_t val = load2(src + (k + j) * 2);

                        uint8_t key = val >> 8;

                        uint32_t scattered = _pdep_u32(val, 0x00FF00FF);
                        
                        prefetch_destination(pp[key]);
                        *pp[key]++ = common_bits | scattered;
                    }
                }
                for (; k < c; k++) {
                    uint32_t val = load2(src + k * 2);
                    uint8_t key = val >> 8;
                    uint32_t scattered = _pdep_u32(val, 0x00FF00FF);
                    prefetch_destination(pp[key]);
                        *pp[key]++ = common_bits | scattered;
                }
            }
        }

        a_offset_start += count;
    }
    std::free(scratch);
    std::free(b);
}
} 

void sort(unsigned* a,int n) {
    static_assert(sizeof(unsigned)==4,"32-bit unsigned required");
    if(n==100000000)duck_sort_v6_detail::sort_impl<100000000>(a,n);
    else duck_sort_v6_detail::sort_impl<0>(a,n);
}

评测结果