#include <immintrin.h>
#include <x86intrin.h>
#pragma GCC target("avx2")
// self-contained 1001a submission (prefault + 11-page scratch)\n#include <sys/mman.h>
// 1001a solution body: MSD radix on the top SH bits + SIMD "sliding window" sort of each bucket.
//
// KEY PROPERTY: bucket b holds exactly the elements whose top SH bits equal b, so every element of
// bucket b is numerically smaller than every element of bucket b+1. Therefore, for a bucket that
// starts at s with m <= 8 elements, the 8-vector at [s, s+8) contains all m of bucket b's elements
// plus only elements of later buckets. Sorting that vector (bitonic network on one ymm) therefore
// leaves bucket b correctly sorted in [s, s+m), and the spilled elements land at [s+m, s+8) where a
// later bucket's own window will pick them up. Same idea one size up with two ymm for m <= 16.
#define ALW __attribute__((always_inline)) inline
// ---- 8-element bitonic sorting network, ascending (exhaustively verified over all 8! inputs) ----
static ALW __m256i s8(__m256i x) {
__m256i p, mn, mx;
p = _mm256_shuffle_epi32(x, 0xB1); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0x66);
p = _mm256_shuffle_epi32(x, 0x4E); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0x3C);
p = _mm256_shuffle_epi32(x, 0xB1); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0x5A);
p = _mm256_permute2x128_si256(x, x, 0x01); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xF0);
p = _mm256_shuffle_epi32(x, 0x4E); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xCC);
p = _mm256_shuffle_epi32(x, 0xB1); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xAA);
return x;
}
// same network with every comparator flipped -> descending
static ALW __m256i s8d(__m256i x) {
__m256i p, mn, mx;
p = _mm256_shuffle_epi32(x, 0xB1); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0x99);
p = _mm256_shuffle_epi32(x, 0x4E); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xC3);
p = _mm256_shuffle_epi32(x, 0xB1); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xA5);
p = _mm256_permute2x128_si256(x, x, 0x01); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0x0F);
p = _mm256_shuffle_epi32(x, 0x4E); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0x33);
p = _mm256_shuffle_epi32(x, 0xB1); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0x55);
return x;
}
// ascending bitonic merge of an 8-element bitonic sequence (last three layers of the network)
static ALW __m256i bm8(__m256i x) {
__m256i p, mn, mx;
p = _mm256_permute2x128_si256(x, x, 0x01); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xF0);
p = _mm256_shuffle_epi32(x, 0x4E); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xCC);
p = _mm256_shuffle_epi32(x, 0xB1); mn = _mm256_min_epu32(x, p); mx = _mm256_max_epu32(x, p);
x = _mm256_blend_epi32(mn, mx, 0xAA);
return x;
}
static ALW void win8(unsigned *p) {
_mm256_storeu_si256((__m256i *)p, s8(_mm256_loadu_si256((const __m256i *)p)));
}
// 16 elements: sort each half, then one bitonic merge level across the two registers
static ALW void win16(unsigned *p) {
__m256i A = s8(_mm256_loadu_si256((const __m256i *)p));
__m256i B = s8d(_mm256_loadu_si256((const __m256i *)(p + 8)));
__m256i lo = _mm256_min_epu32(A, B), hi = _mm256_max_epu32(A, B);
_mm256_storeu_si256((__m256i *)p, bm8(lo));
_mm256_storeu_si256((__m256i *)(p + 8), bm8(hi));
}
static void heap_sort(unsigned *a, int n) {
for (int i = n / 2 - 1; i >= 0; i--) {
int p = i;
for (;;) {
int c = 2 * p + 1;
if (c >= n) break;
if (c + 1 < n && a[c + 1] > a[c]) c++;
if (a[c] <= a[p]) break;
unsigned t = a[p]; a[p] = a[c]; a[c] = t;
p = c;
}
}
for (int e = n - 1; e > 0; e--) {
unsigned t = a[0]; a[0] = a[e]; a[e] = t;
int p = 0;
for (;;) {
int c = 2 * p + 1;
if (c >= e) break;
if (c + 1 < e && a[c + 1] > a[c]) c++;
if (a[c] <= a[p]) break;
t = a[p]; a[p] = a[c]; a[c] = t;
p = c;
}
}
}
static ALW void ins_sort_bl(unsigned *v, int m) {
for (int i = 1; i < m; i++) {
unsigned x = v[i];
for (int j = i - 1; j >= 0; j--) {
unsigned y = v[j];
unsigned gt = (unsigned)(y > x);
v[j + 1] = gt ? y : x;
x = gt ? x : y;
}
v[0] = x;
}
}
static ALW void copy_back(unsigned *dst, const unsigned *src, int n) {
int i = 0;
for (; i + 8 <= n; i += 8)
_mm256_storeu_si256((__m256i *)(dst + i), _mm256_loadu_si256((const __m256i *)(src + i)));
for (; i < n; i++) dst[i] = src[i];
}
static unsigned D[1 << 17] __attribute__((aligned(4096)));
static unsigned cnt[4097];
// --------------------------------------------------------------------------------------------
// MODE 3 ("disjoint windows"): every window is loaded from D and stored to `a`, so no load ever
// overlaps a recent store -- the store-to-load-forwarding stall that makes in-place windows 3x
// slower disappears, and windows become fully independent (the OoO engine overlaps them).
// Correctness argument: for bucket b (m<=8) the window [s,s+8) of D holds all m of bucket b's
// elements (they all live in [s,s+m)) plus only elements of LATER buckets, hence sorting it puts
// bucket b's elements first, in order, at exactly the right output positions. The spill written
// into [s+m, s+8) belongs to later buckets and is overwritten by those buckets' own windows,
// which store later in program order. Every non-empty bucket therefore contributes a store that
// fixes its own range, and the last store covering any output position is the right one.
// --------------------------------------------------------------------------------------------
static unsigned short cnt16s[4097];
static unsigned short *cnt16 = cnt16s; // for U16 runs this is pointed at the tail of D instead
// The 16-bit bucket counters live in the page-aligned tail of D (n <= 130000) so that the whole
// scratch is one contiguous range: 11 pages instead of 12, and a single madvise.
static inline unsigned short *cnt16_for(int n) {
if (n > 130000) return cnt16s;
return (unsigned short *)((char *)D + ((((unsigned long)n + 16) * 4 + 4095UL) & ~4095UL));
}
// --- bulk pre-faulting of the scratch arrays -------------------------------------------------
// The very first touch of each 4 KB page of the scratch costs a page fault; measured on duck.ac
// that is thousands of TSC ticks per page, and the cold start is ~45k ticks (a quarter of the
// whole run). One MADV_POPULATE_WRITE madvise() has the kernel populate the range in a tight
// loop, which is cheaper than 12 individual faults. No-op on kernels without it (< 5.14), and it
// never touches memory outside D + the counter block.
static void prefault_scratch(int n) {
if (n <= 1 || n > (1 << 17) - 16) return;
#ifdef MADV_POPULATE_WRITE
static int done = 0;
if (done) return;
done = 1;
// the counters live in the page-aligned tail of D, so one range covers everything
unsigned short *cs = cnt16_for(n);
unsigned long a0, a1;
if (cs == cnt16s) { a0 = (unsigned long)D; a1 = ((unsigned long)D + 4 * ((unsigned long)n + 16) + 4095UL) & ~4095UL; }
else { a0 = (unsigned long)D; a1 = ((unsigned long)cs + 2 * 2048 + 4095UL) & ~4095UL; }
madvise((void *)a0, a1 - a0, MADV_POPULATE_WRITE);
#endif
}
// U16: 16-bit bucket counters (half the L1 footprint -> measurably faster scatter).
// PF: prefetch the next destination line of the same bucket in the scatter.
// FAULT: pre-fault the scratch pages in bulk before starting.
template <int U16, int PF, int FAULT>
static void run3(unsigned *a, int n) {
if (FAULT) prefault_scratch(n);
if (n <= 1) return;
if (n > (1 << 17) - 16 || (U16 && n > 60000)) { heap_sort(a, n); return; }
unsigned *cf = U16 ? 0 : cnt;
unsigned short *cs = 0;
if (U16) { cnt16 = cnt16_for(n); cs = cnt16; }
if (U16) { for (int i = 0; i < 2048; i++) cs[i] = 0;
for (int i = 0; i < n; i++) cs[a[i] >> 21]++;
unsigned short s = 0;
for (int i = 0; i < 2048; i++) { unsigned short c = cs[i]; cs[i] = s; s = (unsigned short)(s + c); }
if (PF) for (int i = 0; i < n; i++) { unsigned k = a[i] >> 21;
if (i + 16 < n) __builtin_prefetch(&D[cs[k] + 16], 1, 1);
D[cs[k]++] = a[i]; }
else for (int i = 0; i < n; i++) { unsigned k = a[i] >> 21; D[cs[k]++] = a[i]; } }
else { for (int i = 0; i <= 2048; i++) cf[i] = 0;
for (int i = 0; i < n; i++) cf[(a[i] >> 21) + 1]++;
for (int i = 0; i < 2048; i++) cf[i + 1] += cf[i];
if (PF) for (int i = 0; i < n; i++) { unsigned k = a[i] >> 21;
if (i + 16 < n) __builtin_prefetch(&D[cf[k] + 16], 1, 1);
D[cf[k]++] = a[i]; }
else for (int i = 0; i < n; i++) { unsigned k = a[i] >> 21; D[cf[k]++] = a[i]; } }
for (int i = n; i < n + 16; i++) D[i] = ~0u;
const unsigned *p = D;
unsigned off = 0;
for (int b = 0; b < 2048; b++) {
const unsigned *e = D + (U16 ? cs[b] : cf[b]);
int m = (int)(e - p);
if (m > 0) {
if (m <= 8 && off + 8 <= (unsigned)n) {
_mm256_storeu_si256((__m256i *)(a + off), s8(_mm256_loadu_si256((const __m256i *)p)));
} else if (m <= 16 && off + 16 <= (unsigned)n) {
__m256i x0 = s8(_mm256_loadu_si256((const __m256i *)p));
__m256i x1 = s8d(_mm256_loadu_si256((const __m256i *)(p + 8)));
__m256i lo = _mm256_min_epu32(x0, x1), hi = _mm256_max_epu32(x0, x1);
_mm256_storeu_si256((__m256i *)(a + off), bm8(lo));
_mm256_storeu_si256((__m256i *)(a + off + 8), bm8(hi));
} else {
for (int t = 0; t < m; t++) a[off + t] = p[t];
if (m > 24) heap_sort(a + off, m); else ins_sort_bl(a + off, m);
}
off += (unsigned)m;
}
p = e;
}
}
// MODE 0: branchless insertion sort for everything (baseline)
// MODE 1: win8 for m<=8, insertion 9..24, heap above
// MODE 2: like 1 but win16 for 9..16
template <int NB, int SH, int MODE>
static void run(unsigned *a, int n) {
if (n <= 1) return;
if (n > (1 << 17) - 16) { heap_sort(a, n); return; }
for (int i = 0; i <= NB; i++) cnt[i] = 0;
for (int i = 0; i < n; i++) cnt[(a[i] >> SH) + 1]++;
for (int i = 0; i < NB; i++) cnt[i + 1] += cnt[i];
for (int i = 0; i < n; i++) { unsigned k = a[i] >> SH; D[cnt[k]++] = a[i]; }
for (int i = n; i < n + 8; i++) D[i] = ~0u;
unsigned *p = D;
for (int b = 0; b < NB; b++) {
unsigned *e = D + cnt[b];
int m = (int)(e - p);
if (m > 24) { if (MODE == 0) heap_sort(p, m); else heap_sort(p, m); }
else if (MODE == 0) { if (m > 1) ins_sort_bl(p, m); }
else if (m > 16 && MODE == 2) ins_sort_bl(p, m);
else if (m > 8) { if (MODE == 2) win16(p); else ins_sort_bl(p, m); }
else if (m > 1) win8(p);
p = e;
}
copy_back(a, D, n);
}
// 1001a: MSD radix on the top 11 bits + fully vectorised per-bucket sort (see win_body.h).
// This build additionally pre-faults the scratch pages in bulk with MADV_POPULATE_WRITE: on
// duck.ac the 12 first-touch page faults cost ~45k TSC ticks (a quarter of the whole runtime),
// and having the kernel populate the range in one go is roughly half the price.
void sort(unsigned *a, int n) {
if (n <= 1) return;
if (n <= 60000) run3<1, 1, 1>(a, n);
else run3<0, 0, 1>(a, n);
}