#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#include <bits/stdc++.h>
#include <immintrin.h>
#include <x86intrin.h>
#ifdef DUCK_SORT_TUNE_REPORT
unsigned long long duck_tune_choice[4]={};
unsigned long long duck_tune_cost[4][6]={};
#endif
// Four-pass compressed radix sort with on-machine kernel selection.
// Calibration consumes ordinary work once; there is no artificial delay and no
// repeated full-array sort. Layouts, counts and overflow recovery are unchanged.
// Default kernel choices use the former 16/16/21/32 scalar implementations.
namespace duck_sort_v10_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;
}
// Six layout-compatible implementations per phase. Every trial does useful work.
constexpr int Choices=6;
struct Tuner {
unsigned trial=0, selected=0, phase;
uint64_t samples[Choices][3] = {};
uint64_t pending_ticks=0; unsigned pending_count=0;
explicit Tuner(unsigned phase_id):phase(phase_id){}
bool training(unsigned count) const {return trial<Choices*3 && count>=4096;}
unsigned candidate() const {
if(trial>=Choices*3)return selected;
return (trial+2*(trial/Choices)+phase)%Choices;
}
void record(unsigned candidate,uint64_t ticks,unsigned count) {
pending_ticks+=ticks;pending_count+=count;
if(phase==0 && pending_count<262144)return;
samples[candidate][trial/Choices]=(pending_ticks*65536)/pending_count;
pending_ticks=0;pending_count=0;
if(++trial==Choices*3) {
uint64_t med[Choices];
for(unsigned j=0;j<Choices;++j) {
uint64_t a=samples[j][0],b=samples[j][1],c=samples[j][2];
med[j]=a>b ? (b>c?b:(a>c?c:a)) : (a>c?a:(b>c?c:b));
}
unsigned best=0;
for(unsigned j=1;j<Choices;++j)if(med[j]<med[best])best=j;
selected=(med[best]*100<med[0]*98)?best:0;
#ifdef DUCK_SORT_TUNE_REPORT
duck_tune_choice[phase]=selected;
for(unsigned j=0;j<Choices;++j)duck_tune_cost[phase][j]=med[j];
#endif
}
}
};
inline uint64_t stamp() {
_mm_lfence();
uint64_t t=__rdtsc();
_mm_lfence();
return t;
}
template<int U,int Hint,int Distance=64>
__attribute__((noinline)) void pass1_kernel(const unsigned* __restrict__ src,int n,
uint8_t** __restrict__ pp) {
int i=0;
for(;i+U<=n;i+=U) {
if constexpr (Distance!=0)
_mm_prefetch((const char*)((uintptr_t)src+size_t(i+Distance)*4),static_cast<_mm_hint>(Hint));
#pragma GCC unroll 64
for(int j=0;j<U;++j){unsigned v=src[i+j],key=v>>24;store3(pp[key],v);pp[key]+=3;}
}
for(;i<n;++i){unsigned v=src[i],key=v>>24;store3(pp[key],v);pp[key]+=3;}
}
__attribute__((noinline)) void asm1_byte(const unsigned* in,int n,uint8_t** pp){
const uint8_t* src=reinterpret_cast<const uint8_t*>(in);
const uint8_t* end=src+size_t(n&~15)*4;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetchnta 256(%[s])\n\t"
"movl 0(%[s]),%%eax\n\t"
"movzbl 3(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 4(%[s]),%%eax\n\t"
"movzbl 7(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 8(%[s]),%%eax\n\t"
"movzbl 11(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 12(%[s]),%%eax\n\t"
"movzbl 15(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 16(%[s]),%%eax\n\t"
"movzbl 19(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 20(%[s]),%%eax\n\t"
"movzbl 23(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 24(%[s]),%%eax\n\t"
"movzbl 27(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 28(%[s]),%%eax\n\t"
"movzbl 31(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 32(%[s]),%%eax\n\t"
"movzbl 35(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 36(%[s]),%%eax\n\t"
"movzbl 39(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 40(%[s]),%%eax\n\t"
"movzbl 43(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 44(%[s]),%%eax\n\t"
"movzbl 47(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 48(%[s]),%%eax\n\t"
"movzbl 51(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 52(%[s]),%%eax\n\t"
"movzbl 55(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 56(%[s]),%%eax\n\t"
"movzbl 59(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 60(%[s]),%%eax\n\t"
"movzbl 63(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"addq $64,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&15);++i){unsigned v;std::memcpy(&v,src+4*i,4);unsigned key=v>>24;store3(pp[key],v);pp[key]+=3;}
}
__attribute__((noinline)) void asm1_word(const unsigned* in,int n,uint8_t** pp){
const uint8_t* src=reinterpret_cast<const uint8_t*>(in);
const uint8_t* end=src+size_t(n&~15)*4;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetchnta 256(%[s])\n\t"
"movl 0(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 4(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 8(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 12(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 16(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 20(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 24(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 28(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 32(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 36(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 40(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 44(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 48(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 52(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 56(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"movl 60(%[s]),%%eax\n\t"
"movl %%eax,%%ecx\n\t"
"shrl $24,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $3,(%[p],%%rcx,8)\n\t"
"addq $64,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&15);++i){unsigned v;std::memcpy(&v,src+4*i,4);unsigned key=v>>24;store3(pp[key],v);pp[key]+=3;}
}
using P1=void(*)(const unsigned*,int,uint8_t**);
const P1 kernels1[Choices]={pass1_kernel<16,_MM_HINT_NTA>,pass1_kernel<8,_MM_HINT_NTA>,
pass1_kernel<32,_MM_HINT_NTA>,asm1_byte,asm1_word,pass1_kernel<16,_MM_HINT_T0>};
template<int U,int Load>
__attribute__((noinline)) void pass2_kernel(const uint8_t* __restrict__ src,int n,
uint8_t** __restrict__ pp,unsigned* __restrict__ c1,unsigned* __restrict__ c2) {
int k=0;
for(;k+U<=n;k+=U) {
_mm_prefetch((const char*)((uintptr_t)src+size_t(k+64)*3),_MM_HINT_T0);
#pragma GCC unroll 64
for(int j=0;j<U;++j) {
const uint8_t* s=src+(k+j)*3;
if constexpr(Load==0) {
unsigned v=load3(s),key=v&255;
++c1[(v>>8)&255];++c2[(v>>16)&255];store2(pp[key],v>>8);pp[key]+=2;
} else if constexpr(Load==1) {
unsigned v=load2(s+1),key=s[0];
++c1[v&255];++c2[v>>8];store2(pp[key],v);pp[key]+=2;
} else {
unsigned key=s[0];++c1[s[1]];++c2[s[2]];
store2(pp[key],load2(s+1));pp[key]+=2;
}
}
}
for(;k<n;++k) {
unsigned v=load3(src+3*k),key=v&255;
++c1[(v>>8)&255];++c2[(v>>16)&255];store2(pp[key],v>>8);pp[key]+=2;
}
}
__attribute__((noinline)) void asm2_word(const uint8_t* src,int n,uint8_t** pp,unsigned* c1,unsigned* c2){
const uint8_t* end=src+size_t(n&~7)*3;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetcht0 192(%[s])\n\t"
"movzwl 1(%[s]),%%eax\n\t"
"movzbl 0(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movzwl 4(%[s]),%%eax\n\t"
"movzbl 3(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movzwl 7(%[s]),%%eax\n\t"
"movzbl 6(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movzwl 10(%[s]),%%eax\n\t"
"movzbl 9(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movzwl 13(%[s]),%%eax\n\t"
"movzbl 12(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movzwl 16(%[s]),%%eax\n\t"
"movzbl 15(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movzwl 19(%[s]),%%eax\n\t"
"movzbl 18(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movzwl 22(%[s]),%%eax\n\t"
"movzbl 21(%[s]),%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"addq $24,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp), [c1] "r"(c1), [c2] "r"(c2)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&7);++i){unsigned v=load3(src+3*i),key=v&255;++c1[(v>>8)&255];++c2[(v>>16)&255];store2(pp[key],v>>8);pp[key]+=2;}
}
__attribute__((noinline)) void asm2_dword(const uint8_t* src,int n,uint8_t** pp,unsigned* c1,unsigned* c2){
const uint8_t* end=src+size_t(n&~7)*3;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetcht0 192(%[s])\n\t"
"movl 0(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movl 3(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movl 6(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movl 9(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movl 12(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movl 15(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movl 18(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"movl 21(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"shrl $8,%%eax\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzbl %%al,%%ecx\n\t"
"incl (%[c1],%%rcx,4)\n\t"
"movzbl %%ah,%%ecx\n\t"
"incl (%[c2],%%rcx,4)\n\t"
"addq $24,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp), [c1] "r"(c1), [c2] "r"(c2)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&7);++i){unsigned v=load3(src+3*i),key=v&255;++c1[(v>>8)&255];++c2[(v>>16)&255];store2(pp[key],v>>8);pp[key]+=2;}
}
using P2=void(*)(const uint8_t*,int,uint8_t**,unsigned*,unsigned*);
const P2 kernels2[Choices]={pass2_kernel<16,0>,pass2_kernel<8,0>,pass2_kernel<32,0>,
pass2_kernel<16,1>,asm2_word,asm2_dword};
inline uint16_t low_byte(uint16_t v,unsigned digit) {
unsigned r=v;
asm("movb %b1,%b0" : "+r"(r) : "r"(digit));
return uint16_t(r);
}
template<int U,bool ByteMove>
__attribute__((noinline)) void pass3_kernel(const uint8_t* __restrict__ input,
const unsigned* __restrict__ counts,const unsigned* __restrict__ starts,
uint8_t** __restrict__ pp) {
for(unsigned b0=0;b0<256;++b0) {
unsigned c=counts[b0];const uint8_t* src=input+starts[b0];unsigned k=0;
for(;k+U<=c;k+=U) {
_mm_prefetch((const char*)((uintptr_t)src+size_t(k+64)*2),_MM_HINT_T0);
#pragma GCC unroll 64
for(unsigned j=0;j<U;++j) {
unsigned v=load2(src+(k+j)*2),key=v&255;
uint16_t value;
if constexpr(ByteMove)value=low_byte(v,b0);
else value=uint16_t(b0|(v&0xff00));
store2(pp[key],value);pp[key]+=2;
}
}
for(;k<c;++k){unsigned v=load2(src+k*2),key=v&255;store2(pp[key],b0|(v&0xff00));pp[key]+=2;}
}
}
__attribute__((noinline)) void asm3_8(const uint8_t* input,const unsigned* counts,const unsigned* starts,uint8_t** pp){
for(unsigned b0=0;b0<256;++b0){
int n=counts[b0];const uint8_t* src=input+starts[b0];
const uint8_t* end=src+size_t(n&~7)*2;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetcht0 128(%[s])\n\t"
"movzwl 0(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 2(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 4(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 6(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 8(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 10(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 12(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 14(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"addq $16,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp), [digit] "r"(b0)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&7);++i){unsigned v=load2(src+2*i),key=v&255;store2(pp[key],b0|(v&0xff00));pp[key]+=2;}
}
}
__attribute__((noinline)) void asm3_16(const uint8_t* input,const unsigned* counts,const unsigned* starts,uint8_t** pp){
for(unsigned b0=0;b0<256;++b0){
int n=counts[b0];const uint8_t* src=input+starts[b0];
const uint8_t* end=src+size_t(n&~15)*2;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetcht0 128(%[s])\n\t"
"movzwl 0(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 2(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 4(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 6(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 8(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 10(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 12(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 14(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 16(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 18(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 20(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 22(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 24(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 26(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 28(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"movzwl 30(%[s]),%%eax\n\t"
"movzbl %%al,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"movb %b[digit],%%al\n\t"
"movw %%ax,(%%r10)\n\t"
"addq $2,(%[p],%%rcx,8)\n\t"
"addq $32,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp), [digit] "r"(b0)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&15);++i){unsigned v=load2(src+2*i),key=v&255;store2(pp[key],b0|(v&0xff00));pp[key]+=2;}
}
}
using P3=void(*)(const uint8_t*,const unsigned*,const unsigned*,uint8_t**);
const P3 kernels3[Choices]={pass3_kernel<21,false>,pass3_kernel<8,false>,pass3_kernel<16,false>,
pass3_kernel<32,false>,asm3_8,asm3_16};
template<int U,int Construct>
__attribute__((noinline)) void pass4_kernel(const uint8_t* __restrict__ input,
const unsigned* __restrict__ counts,const unsigned* __restrict__ starts,
unsigned** __restrict__ pp,unsigned high) {
for(unsigned b1=0;b1<256;++b1) {
unsigned c=counts[b1],common=high|(b1<<8);const uint8_t* src=input+starts[b1];unsigned k=0;
for(;k+U<=c;k+=U) {
_mm_prefetch((const char*)((uintptr_t)src+size_t(k+64)*2),_MM_HINT_T0);
#pragma GCC unroll 64
for(unsigned j=0;j<U;++j) {
unsigned v=load2(src+(k+j)*2),key=v>>8,value;
if constexpr(Construct==0)value=common|_pdep_u32(v,0x00ff00ff);
else if constexpr(Construct==1)value=common|(v&255)|((v&0xff00)<<8);
else {key=src[(k+j)*2+1];value=common|_pdep_u32(v,0x00ff00ff);}
*pp[key]++=value;
}
}
for(;k<c;++k){unsigned v=load2(src+k*2);*pp[v>>8]++=common|_pdep_u32(v,0x00ff00ff);}
}
}
__attribute__((noinline)) void asm4_8(const uint8_t* input,const unsigned* counts,const unsigned* starts,unsigned** pp,unsigned high){
for(unsigned b1=0;b1<256;++b1){
int n=counts[b1];const uint8_t* src=input+starts[b1];unsigned common=high|(b1<<8);
const uint8_t* end=src+size_t(n&~7)*2;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetcht0 128(%[s])\n\t"
"movzwl 0(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 2(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 4(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 6(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 8(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 10(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 12(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 14(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"addq $16,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp), [mask] "r"(0x00ff00ffu), [common] "r"(common)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&7);++i){unsigned v=load2(src+2*i);*pp[v>>8]++=common|_pdep_u32(v,0x00ff00ff);}
}
}
__attribute__((noinline)) void asm4_16(const uint8_t* input,const unsigned* counts,const unsigned* starts,unsigned** pp,unsigned high){
for(unsigned b1=0;b1<256;++b1){
int n=counts[b1];const uint8_t* src=input+starts[b1];unsigned common=high|(b1<<8);
const uint8_t* end=src+size_t(n&~15)*2;
if(src!=end) {
asm volatile(
".p2align 5\n\t"
"1:\n\t"
"prefetcht0 128(%[s])\n\t"
"movzwl 0(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 2(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 4(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 6(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 8(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 10(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 12(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 14(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 16(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 18(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 20(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 22(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 24(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 26(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 28(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"movzwl 30(%[s]),%%eax\n\t"
"movzbl %%ah,%%ecx\n\t"
"movq (%[p],%%rcx,8),%%r10\n\t"
"pdep %k[mask],%%eax,%%eax\n\t"
"orl %k[common],%%eax\n\t"
"movl %%eax,(%%r10)\n\t"
"addq $4,(%[p],%%rcx,8)\n\t"
"addq $32,%[s]\n\t"
"cmpq %[e],%[s]\n\t"
"jne 1b\n\t"
: [s] "+&r"(src)
: [e] "r"(end), [p] "r"(pp), [mask] "r"(0x00ff00ffu), [common] "r"(common)
: "rax","rcx","r10","memory","cc"
);
}
for(int i=0;i<(n&15);++i){unsigned v=load2(src+2*i);*pp[v>>8]++=common|_pdep_u32(v,0x00ff00ff);}
}
}
using P4=void(*)(const uint8_t*,const unsigned*,const unsigned*,unsigned**,unsigned);
const P4 kernels4[Choices]={pass4_kernel<32,0>,pass4_kernel<8,0>,pass4_kernel<16,0>,
pass4_kernel<64,0>,asm4_8,asm4_16};
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;
Tuner tuner1(0),tuner2(1),tuner3(2),tuner4(3);
// 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);
const unsigned count=end-i;
const bool measure=tuner1.training(count) && i>=2*TILE;
const unsigned choice=measure?tuner1.candidate():tuner1.selected;
uint64_t before=measure?stamp():0;
kernels1[choice](src+i,count,pp);
if(measure)tuner1.record(choice,stamp()-before,count);
i=end;
}
// 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;
// ---------------------------------------------------------
// 分段处理：遍历 B3 的每一个 Bucket
// ---------------------------------------------------------
uint8_t* a_u8 = (uint8_t*)a;
// 局部直方图缓存
uint cnt0[256];
uint cnt1[256];
uint cnt2[256];
uint ptr0[256]; // Pass 2 (Write a) pointers
uint ptr1[256]; // Pass 3 (Write b) pointers
uint ptr2[256]; // Pass 4 (Write a Final) pointers
uint32_t a_offset_start = 0; // index offset in a (Writing Final)
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;
}
// Pass 2 Output (Temporary in a)
// Prefer the unused suffix of a; the range is checked before scattering.
uint8_t* seg_a_temp = a_u8 + (a_offset_start * 4);
// -----------------------------------------------------
// Step 2: Use Sampling for Pointers & Init Counters
// -----------------------------------------------------
// Use sampling to estimate B0 upper bounds for Pass 2 (ptr0)
// Budget logic:
// Pass 2 writes 2 bytes per item into a space reserved for 4 bytes per item (Final Output array).
// Effectively, we have capacity for count * 2 items of size 2 bytes.
int budget_pass2 = count * 2;
int sample_size = 5000;
// get_population_upper_bounds will default to stride 3 because we pass stride=3?
// Wait, function signature is (uint8_t* A, int N, int budget, int stride, int sample_size)
// seg_b_in has data [B0, B1, B2]... So stride=3, offset=0 is B0.
auto bounds = get_population_upper_bounds<3>(seg_b_in, count, budget_pass2, sample_size);
// Calculate ptr0 (bucket start offsets in 'a')
uint32_t tmp = 0;
for(int k=0; k<256; k++) {
ptr0[k] = tmp;
tmp += bounds[k] * 2;
}
// Before scattering, reserve enough address range even when EVERY
// element enters the last bucket. The source b is not modified here.
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) {
// Earlier high-byte buckets are already final in a.
// All remaining raw values are still intact in b.
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;
}
// Initialize cnt1, cnt2 for exact counting during Pass 2
memset(cnt1, 0, sizeof(cnt1));
memset(cnt2, 0, sizeof(cnt2));
// cnt0 will be recovered from pointer progress after Pass 2
// -----------------------------------------------------
// Pass 2: LSD Step 1 (Key B0)
// Read b (3B: B0,B1,B2) -> Write a (2B: B1,B2)
// AND compute exact histograms for B1, B2
// -----------------------------------------------------
{
uint p[256];
uint8_t* pp[256];
for(int z=0;z<256;++z) pp[z]=seg_a_temp+ptr0[z];
uint8_t* src = seg_b_in;
uint8_t* dst = seg_a_temp;
const bool measure=tuner2.training(count);
const unsigned choice=measure?tuner2.candidate():tuner2.selected;
uint64_t before=measure?stamp():0;
kernels2[choice](src,count,pp,cnt1,cnt2);
if(measure)tuner2.record(choice,stamp()-before,count);
// Post-Pass 2: Recover exact cnt0 and compute ptr1, ptr2
bool retry = false;
uint32_t tmp1 = 0;
uint32_t tmp2 = 0;
for(int k=0; k<256; k++) {
// Recover cnt0 from pointer progress (p - ptr0) / 2
cnt0[k] = (pp[k] - (seg_a_temp+ptr0[k])) >> 1;
retry |= cnt0[k] > bounds[k];
ptr1[k] = tmp1;
tmp1 += cnt1[k] * 2 + 4; // Padding
ptr2[k] = tmp2;
tmp2 += cnt2[k]; // No Padding (Dense)
}
if (retry) {
// Recalculate ptr0 with exact counts
uint32_t tmp = 0;
for(int k=0; k<256; k++) {
ptr0[k] = tmp;
tmp += cnt0[k] * 2 + 4;
}
// Rerun Pass 2 (Distribution Only, no counting)
uint p_retry[256];
memcpy(p_retry, ptr0, sizeof(p_retry));
uint8_t* src = seg_b_in;
uint8_t* dst = seg_a_temp;
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; // B0
store2(dst + 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;
store2(dst + p_retry[key], val >> 8);
p_retry[key] += 2;
}
}
}
// -----------------------------------------------------
// Pass 3: LSD Step 2 (Key B1)
// Read a (2B: B1,B2) -> Write b (2B: B0,B2)
// Iterate B0 buckets to restore B0
// -----------------------------------------------------
{
// count>=1024, so 2*count+1024 fits in the consumed 3-byte input bucket.
uint8_t* dst_base = seg_b_in;
uint p[256];
uint8_t* pp[256];
for(int z=0;z<256;++z) pp[z]=dst_base+ptr1[z];
const bool measure=tuner3.training(count);
const unsigned choice=measure?tuner3.candidate():tuner3.selected;
uint64_t before=measure?stamp():0;
kernels3[choice](seg_a_temp,cnt0,ptr0,pp);
if(measure)tuner3.record(choice,stamp()-before,count);
}
// -----------------------------------------------------
// Pass 4: LSD Step 3 (Key B2) & Finalize
// Read b (2B: B0,B2) -> Write a (4B: Full)
// Iterate B1 buckets to restore B1
// -----------------------------------------------------
{
uint* dst_base = a + a_offset_start;
uint8_t* src_base = seg_b_in;
uint p[256];
unsigned* pp[256];
for(int z=0;z<256;++z)pp[z]=dst_base+ptr2[z];
uint32_t val_b3_shifted = i_b3 << 24;
const bool measure=tuner4.training(count);
const unsigned choice=measure?tuner4.candidate():tuner4.selected;
uint64_t before=measure?stamp():0;
kernels4[choice](src_base,cnt1,ptr1,pp,val_b3_shifted);
if(measure)tuner4.record(choice,stamp()-before,count);
}
// Update offsets
a_offset_start += count;
}
std::free(scratch);
std::free(b);
}
} // namespace duck_sort_v10_detail
void sort(unsigned* a,int n) {
static_assert(sizeof(unsigned)==4,"32-bit unsigned required");
if(n==100000000)duck_sort_v10_detail::sort_impl<100000000>(a,n);
else duck_sort_v10_detail::sort_impl<0>(a,n);
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 580.248 ms | 669 MB + 864 KB | Accepted | Score: 100 | 显示更多 |