// [+66B lane m183 arm A6 prefetchw on the flush store target] ** E-state, no RFO **
// [+66B lane m183 arm A2 flushtarget prefetch dist16] ** scatter metadata prefetch **
// 1008: 2-D dominance counting, n = 1e7. V1 = merged record.
//
// DELETION vs the boarded artifact (work/l1008_1009_lane/c8_inc64.cpp):
// the boarded engine carries the x-sort key in a SEPARATE 4-byte array `XB[NMAX]`
// next to the 8-byte record `PKB[NMAX]` -- 12 B/element of touched pages.
// `x` is 24-bit and the record's `idx` is 24-bit; after the MSD-8 stage only the
// LOW 16 bits of x are still needed to finish the sort. 16 + 24 + 24 = 64 exactly,
// so `xlow16` goes INTO the record and `XB` is deleted outright:
// REC = [ xlow16:16 | y:24 | idx:24 ]
// Touched working set 320 MB -> 240 MB. Every pass moves 8 B/el instead of 12 B/el.
//
// V7 (this arm): the L3 run-SCAN is deleted. After the L2 counting sort the array is
// already ordered by the 9-bit digit d=(rec>>30)&511, so the runs the old code found by
// scanning for equal digits are EXACTLY the regions [s+gHis[d], s+gHis[d+1]) that the L2
// prefix sum produced. Iterating those boundaries visits the same elements in the same
// order and emits the same STG records -- byte-identical output, one pass over 1e7
// elements' worth of scan (load+shift+mask+cmp per element, in a dependent chain) deleted.
//
// The remaining pipeline (MSD-8 bucket by x>>16, bucket-local LSD on xlow16,
// L1 = MSD-9 over y with the hierarchical prefix counter, fused L2/L3, flush) is
// byte-for-byte the boarded algorithm.
//
// Record layouts (u64):
// REC (x-sorted) : [ xlow16:16 | y:24 | idx:24 ]
// L1 out (ALT) : [ cnt:24 | ylow15:15 | idx:24 ] (bit 39 unused)
// after L2L3 : out[idx] = cnt
// cnt starts at (2^24 - dup) so the modular sum ends at (total - dup), where
// dup = #tie-mates earlier in the same x-group with smaller y (x-strictness).
#include <cstdint>
#include <cstring>
#include <cstdlib>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC target("avx2")
typedef uint64_t u64;
typedef uint32_t u32;
static const int NMAX = 10000002;
alignas(64) static u64 REC[NMAX]; // x-sorted working record
alignas(64) static u64 SCR[NMAX]; // bucket-local LSD scratch, touched only [0..maxlen)
alignas(64) static u64 ALT[NMAX]; // L1 output, ordered by y-MSD then x-order
alignas(64) static u64 STG[NMAX]; // idx-bucket staging array for the flush
static u32 cntb[256];
static u32 gC2[256];
static u32 gC1[256][256]; // per-x-bucket histogram of x bits 15..8 (closed in the stage-A scatter)
static u32 gH512[512];
#define NBKT 1024
#define BSH 15
static u32 gBkt[NBKT + 1];
static u32 gBOff[NBKT + 2];
alignas(32) static u32 gP1[8];
alignas(32) static u32 gP2[8][8];
alignas(32) static u32 gP3[8][8][8];
alignas(32) static u32 gInc[8][8];
static u32 gCur[512];
static u32 gBnd[514];
alignas(32) static u32 gT6[512][8]; // prefix over (d3>>3) inside each d2 group
alignas(32) static u32 gR6[512][8][8]; // prefix over (d3&7) inside each (d2,d3>>3)
static u32 *gP3p[8][8]; // precomputed row bases: &gP3[a][b][0]
static u32 *gP2p[8]; // &gP2[a][0]
static u32 *gT6p[512]; // &gT6[d][0]
static u32 *gR6p[512][8]; // &gR6[d][t3][0]
static inline void cnt_init(void) {
memset(gP1, 0, sizeof(gP1));
memset(gP2, 0, sizeof(gP2));
memset(gP3, 0, sizeof(gP3));
}
static inline u32 cnt_q(u32 d) {
return gP1[d >> 6] + gP2[d >> 6][(d >> 3) & 7] + gP3[d >> 6][(d >> 3) & 7][d & 7];
}
// ---- specialised pair for a 6-bit digit (d < 64): the top hierarchy level is
// provably never READ (cnt_q only ever reads gP1[0], which this path never writes)
// and cnt_init() resets all three levels before every run, so it cannot leak.
static inline u32 cnt_q64(u32 d) {
return gP2[0][(d >> 3) & 7] + gP3[0][(d >> 3) & 7][d & 7];
}
static inline void cnt_ins64(u32 d) {
u32 b = (d >> 3) & 7, c = d & 7;
__m256i i3 = _mm256_load_si256((const __m256i *)gInc[c]);
__m256i v3 = _mm256_load_si256((const __m256i *)&gP3[0][b][0]);
_mm256_store_si256((__m256i *)&gP3[0][b][0], _mm256_add_epi32(v3, i3));
__m256i i2 = _mm256_load_si256((const __m256i *)gInc[b]);
__m256i v2 = _mm256_load_si256((const __m256i *)&gP2[0][0]);
_mm256_store_si256((__m256i *)&gP2[0][0], _mm256_add_epi32(v2, i2));
}
static inline void cnt_ins(u32 d) {
u32 a = d >> 6, b = (d >> 3) & 7, c = d & 7;
__m256i i3 = _mm256_load_si256((const __m256i *)gInc[c]);
__m256i v3 = _mm256_load_si256((const __m256i *)&gP3[a][b][0]);
_mm256_store_si256((__m256i *)&gP3[a][b][0], _mm256_add_epi32(v3, i3));
__m256i i2 = _mm256_load_si256((const __m256i *)gInc[b]);
__m256i v2 = _mm256_load_si256((const __m256i *)&gP2[a][0]);
_mm256_store_si256((__m256i *)&gP2[a][0], _mm256_add_epi32(v2, i2));
__m256i i1 = _mm256_load_si256((const __m256i *)gInc[a]);
__m256i v1 = _mm256_load_si256((const __m256i *)&gP1[0]);
_mm256_store_si256((__m256i *)&gP1[0], _mm256_add_epi32(v1, i1));
}
// V5: the top counter level is kept in a YMM REGISTER and queried with vpermd.
// `cnt_ins` adds a 0/1 vector to gP1 EVERY element, while `cnt_q` reads one 32-bit
// lane of the SAME 32-byte vector every element. gcc keeps gP1 in %ymm but must
// still store it (it cannot see that the following 4-byte load aliases the 32-byte
// store), so every element pays a store->load-forward failure on the hottest
// counter level. Holding the level in a register and extracting the lane with
// VPERMD removes the memory round-trip entirely; lanes are still added with the
// same 0/1 vectors, so the counter is bit-identical.
static inline u32 lane8(const __m256i v, u32 a) {
__m256i sel = _mm256_permutevar8x32_epi32(v, _mm256_set1_epi32((int)a));
return (u32)_mm_cvtsi128_si32(_mm256_castsi256_si128(sel));
}
void count_2d(int n, const unsigned *x, const unsigned *y, unsigned *out) {
if (n <= 0) return;
if (n == 1) { out[0] = 0; return; }
if (n <= 2000) {
memset(out, 0, (size_t)n * 4);
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
if (x[j] < x[i] && y[j] < y[i]) out[i]++;
return;
}
if (n > 16000000) { // record fields: idx < 2^24, y < 2^24, xlow16 fits
memset(out, 0, (size_t)n * 4);
return;
}
static int inited = 0;
if (!inited) {
for (int m = 0; m < 8; m++)
for (int i2 = 0; i2 < 8; i2++) gInc[m][i2] = (i2 > m) ? 1u : 0u;
for (int a = 0; a < 8; a++) { gP2p[a] = &gP2[a][0];
for (int b = 0; b < 8; b++) gP3p[a][b] = &gP3[a][b][0]; }
for (int dd = 0; dd < 512; dd++) {
gT6p[dd] = &gT6[dd][0];
for (int t3 = 0; t3 < 8; t3++) gR6p[dd][t3] = &gR6[dd][t3][0];
}
inited = 1;
}
u32 np = (u32)n;
// ------------- sort by x: MSD 8 bits (x>>16) then two BUCKET-LOCAL LSD ----
static u32 bnd[257];
memset(cntb, 0, sizeof(cntb));
memset(gH512, 0, sizeof(gH512));
memset(gC1, 0, sizeof(gC1)); // per-call: the LSD-1 pass CONSUMES this histogram as cursors
for (u32 i = 0; i < np; i++) cntb[(x[i] >> 16) & 255]++;
{ u32 acc = 0;
for (int d = 0; d < 256; d++) { u32 c = cntb[d]; bnd[d] = acc; cntb[d] = acc; acc += c; } }
bnd[256] = np;
for (u32 i = 0; i < np; i++) {
u32 dg = (x[i] >> 16) & 255;
gC1[dg][x[i] & 255]++; // LSD-1 digit = (REC>>48)&255 = x bits 0..7
u32 qq = cntb[dg]++;
REC[qq] = ((u64)(x[i] & 0xFFFFu) << 48) | ((u64)y[i] << 24) | (u64)i;
gH512[(y[i] >> 15) & 511]++;
}
// ---- L1 bucket offsets AND the ALT cursors must exist BEFORE the fused L1 ----
{ u32 acc = 0;
for (int d = 0; d < 512; d++) { u32 c = gH512[d]; gBnd[d] = acc; gCur[d] = acc; acc += c; } }
gBnd[512] = np;
cnt_init(); // counter state must start clean before the fused L1
__m256i P1v = _mm256_setzero_si256();
for (u32 b = 0; b < 256; b++) {
u32 s = bnd[b], e = bnd[b + 1], len = e - s;
if (len > 1) {
// The second LSD pass's histogram is accumulated INSIDE the first pass's histogram
// loop (the digit is a property of the element, not of its position), so one whole
// pass over the bucket is deleted: costs 1 extra 1 KB-table RMW per element, saves
// a 312 KB read pass per bucket.
memset(gC2, 0, sizeof(gC2));
{ u32 a2 = 0;
for (int d = 0; d < 256; d++) { u32 c = gC1[b][d]; gC1[b][d] = a2; a2 += c; } }
for (u32 t = s; t < e; t++) {
u64 r = REC[t];
u32 dg = (u32)(r >> 48) & 255;
u32 qq = gC1[b][dg]++;
SCR[qq] = r;
gC2[(u32)(r >> 56) & 255]++;
}
{ u32 a2 = 0;
for (int d = 0; d < 256; d++) { u32 c = gC2[d]; gC2[d] = a2; a2 += c; } }
for (u32 t = 0; t < len; t++) {
u32 dg = (u32)(SCR[t] >> 56) & 255;
u32 qq = gC2[dg]++;
REC[s + qq] = SCR[t];
}
}
// ---- FUSED L1: the bucket is L1/L2/L3-resident right now, so the y-MSD
// rank-accumulation runs here instead of re-reading the whole array
// from DRAM. The counter state and gCur[] are global and persist
// across buckets, so the result is identical to the unfused pass.
u32 e2 = bnd[b + 1];
// SINGLE PASS: the x-group's end is found by the element that CHANGES xv, so the
// separate scanning pass over the bucket (one dependent load+cmp+branch per element,
// ~12 % of the row) is deleted. `gstart` replaces the group's base index and the
// dup loop runs over [gstart, t) -- the same earlier-in-group elements, so the
// counted set is unchanged.
u32 gstart = bnd[b];
u32 prev_xv = 0xFFFFFFFFu;
for (u32 t = bnd[b]; t < e2; t++) {
{
u64 rec = REC[t];
u32 xv = (u32)(rec >> 48);
if (xv != prev_xv) { gstart = t; prev_xv = xv; }
u32 dup = 0;
if (t > gstart) {
u32 ya = (u32)(rec >> 24) & 0xFFFFFFu;
for (u32 b = gstart; b < t; b++) dup += (u32)(((u32)(REC[b] >> 24) & 0xFFFFFFu) < ya);
}
u32 d = (u32)(rec >> 39) & 511;
u32 a3 = d >> 6, b3 = (d >> 3) & 7, c3 = d & 7;
u32 *p3 = gP3p[a3][b3];
u32 *p2 = gP2p[a3];
u32 inc = lane8(P1v, a3) + p2[b3] + p3[c3];
{ __m256i i3 = _mm256_load_si256((const __m256i *)gInc[c3]);
__m256i v3 = _mm256_load_si256((const __m256i *)p3);
_mm256_store_si256((__m256i *)p3, _mm256_add_epi32(v3, i3));
__m256i i2 = _mm256_load_si256((const __m256i *)gInc[b3]);
__m256i v2 = _mm256_load_si256((const __m256i *)p2);
_mm256_store_si256((__m256i *)p2, _mm256_add_epi32(v2, i2));
P1v = _mm256_add_epi32(P1v, _mm256_load_si256((const __m256i *)gInc[a3])); }
u32 q = gCur[d]++;
u32 cnt = (0x1000000u - dup) & 0xFFFFFFu;
cnt = (cnt + inc) & 0xFFFFFFu;
ALT[q] = ((u64)cnt << 40) | (rec & 0x0000007FFFFFFFFFull);
}
}
}
// idx IS the input index, so the indices are a permutation of [0,np): bucket b of the
// flush staging holds exactly min(np,(b+1)<<BSH) - min(np,b<<BSH) elements. The
// histogram over idx>>BSH is therefore a closed form -- one RMW per element deleted.
for (int b = 0; b < NBKT; b++) { u32 off = (u32)b << BSH; gBOff[b] = off < np ? off : np; }
gBOff[NBKT] = np;
// ---------------- fused rank pass over the y-MSD-9 groups ----------------
// One pass per group, in ALT (x-)order. inc2 = per-group 9-bit counter over y bits
// 14..6; inc3 = a 6-bit counter over y bits 5..0 keyed by that same 9-bit digit, i.e.
// 512 independent 64-lane prefix structures. No counting sort, no REC round trip.
for (u32 b0 = 0; b0 < 512; b0++) {
u32 s = gBnd[b0], e = gBnd[b0 + 1];
u32 m = e - s;
if (m == 0) continue;
if (m == 1) {
{ u32 _ix = (u32)(ALT[s] & 0xFFFFFFu); STG[gBOff[_ix >> BSH]++] = ((ALT[s] >> 40) << 24) | _ix; }
continue;
}
memset(gT6, 0, sizeof(gT6));
memset(gR6, 0, sizeof(gR6));
cnt_init();
__m256i P1v2 = _mm256_setzero_si256();
for (u32 t = s; t < e; t++) {
u64 rec = ALT[t];
u32 d = (u32)(rec >> 30) & 511;
u32 a3 = d >> 6, b3 = (d >> 3) & 7, c3 = d & 7;
u32 inc = lane8(P1v2, a3) + gP2[a3][b3] + gP3[a3][b3][c3];
{ __m256i i3 = _mm256_load_si256((const __m256i *)gInc[c3]);
__m256i v3 = _mm256_load_si256((const __m256i *)&gP3[a3][b3][0]);
_mm256_store_si256((__m256i *)&gP3[a3][b3][0], _mm256_add_epi32(v3, i3));
__m256i i2 = _mm256_load_si256((const __m256i *)gInc[b3]);
__m256i v2 = _mm256_load_si256((const __m256i *)&gP2[a3][0]);
_mm256_store_si256((__m256i *)&gP2[a3][0], _mm256_add_epi32(v2, i2));
P1v2 = _mm256_add_epi32(P1v2, _mm256_load_si256((const __m256i *)gInc[a3])); }
u32 d3 = (u32)(rec >> 24) & 63;
u32 t3 = d3 >> 3, q3 = d3 & 7;
u32 *r6 = gR6p[d][t3];
u32 *t6 = gT6p[d];
u32 inc3 = t6[t3] + r6[q3];
{ __m256i i3 = _mm256_load_si256((const __m256i *)gInc[q3]);
__m256i v3 = _mm256_load_si256((const __m256i *)r6);
_mm256_store_si256((__m256i *)r6, _mm256_add_epi32(v3, i3));
__m256i i2 = _mm256_load_si256((const __m256i *)gInc[t3]);
__m256i v2 = _mm256_load_si256((const __m256i *)t6);
_mm256_store_si256((__m256i *)t6, _mm256_add_epi32(v2, i2)); }
{ u32 _ix = (u32)(rec & 0xFFFFFFu);
u32 cnt = ((u32)(rec >> 40) + inc + inc3) & 0xFFFFFFu;
STG[gBOff[_ix >> BSH]++] = ((u64)cnt << 24) | _ix; }
}
}
// ---- flush: walk the staging array one idx-bucket at a time so every
// out[] store lands inside a 64 KB (BSH=14) L2-resident window ----
{
u32 t = 0;
for (u32 b = 0; b < NBKT; b++) {
u32 e = gBOff[b];
u32 pp = t < e ? t : e;
for (; t < e; t++) {
u64 r = STG[t];
// scatter prefetch: the out[] target is data-dependent, so the RFO of the
// line it lands in is exposed. Prefetch D elements ahead.
if (t + 16 < e) { u32 nx = (u32)(STG[t + 16] & 0xFFFFFFu);
// __builtin_prefetch(p,1,3) emits PREFETCHT0 (a READ prefetch), not a write
// prefetch -- for STORE TARGETS use the real PREFETCHW, which takes the line
// in E state so the store needs no RFO at all.
__asm__ __volatile__("prefetchw %0" :: "m"(out[nx])); }
out[(u32)(r & 0xFFFFFFu)] = (u32)(r >> 24);
}
}
}
}
#pragma GCC pop_options