// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #121867 <https://duck.ac/submission/121867>
// （3118.385642 ms）—— **本文件正文的来源**：本发 = 其正文 + 两处改动
// （① `bins[]` 计数打包 `<<35`→`<<37`；② 首行加 `#pragma GCC optimize("O2,unroll-loops")`），见下"思路"节。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #121429 <https://duck.ac/submission/121429>
// （3162.657490 ms，Accepted）—— 改动①（`[n08c-B]`，`<<37`）的原始出处；#121867 抄件时把它改成了 `<<35`，
// 本发改回我方在册形态。其自带的上游引用链（#121248 / #121281 / #121153 / #120937 / #120790 / #111438 /
// #110329 / #109849 / #109709 / duck.ac 题面）在正文头注中原样保留 ✓
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」，即站上提交按站点规则公开；
// 本文件逐条标注原作者账号、原提交地址与用途 ✓
// ======================
// ===== 思路 =====
// 【本发 = #121867 正文 + 2 处：`bins[]` 打包位 `<<35`→`<<37`、首行 `#pragma GCC optimize("O2,unroll-loops")`】
//
// §一 改动①（`[n08d-B37]`，逐位等价）：`bins[]` 只作中转，其内部位布局不参与输出。
// 写 `(cnt<<37)|(rec&0x7FFFFFFFF)`、读 `(u32)(rec>>37)`；`cnt<2^27` ⇒ 与 `<<35` 语义逐位相同。
// 仪器（cachegrind 判题机几何，n=3e6）：Ir **−9,000,000（−0.88%）**，其中 `bins[q]=` 那一行 27M→21M。
// ★ 实测板面：本件（`<<37`）**#121963 = 3099.098951 ms** ↔ 纯抄件（`<<35`）**#121935 = 3118.555120 ms**
// ⇒ **−19.456 ms（−0.624%）** ⇒ 本族兑现率 **ΔIr −0.88% → Δtime −0.62%（比值 0.71）** ✓ 建议入库。
//
// §二 改动②（`[n08d-U]`，逐位等价，但**必须用闸门复核**）：首行加 `#pragma GCC optimize("O2,unroll-loops")`。
// ★ 依据 = **本件（#121867 系）上的 cachegrind 三栏实测**（判题机同款 `ref/gcc9 -O2 -std=c++17 -static`）：
// | 量 | 无 pragma | +O2,unroll-loops | Δ |
// | Ir | 1,010,407,683 | **968,649,538** | **−41,758,145（−4.13%）** |
// | Dr | 194,575,893 | 194,657,952 | +0.04% |
// | Dw | 88,351,608 | **85,362,653** | −3.38% |
// | I1mr/ILmr（取指） | 882 / 882 | 937 / 937 | ≈0（无取指爆炸 ✓） |
// 其中受益最大的单行 = `sort_x_local_att` 的 pass-1 x 直方图循环 `h0[...]++`：
// 21.0M → 12.75M Ir（−2.75/元素，2 路展开）✓
// ★★ 本条**更正 e8y_ 席的时代结论**（notes L2115-2127："去 `unroll-loops` = −18.341 ms"）：那条是在
// **融合前的 #109694/#109709 代码**上量的；本席在同一判题机口径下对**当代件**复测 ⇒ **反号**
// ⇒ 与 e8y 自己的律"展开轴是逐题/逐件的，必须实测"一致 ✓（**codegen 轴必须随结构性改动重扫**）
// ★ 诚实边界：Ir/Dw 是确定性账，但它**不是墙钟**（本题有"Ir 降而墙钟不降"的在册反例）⇒ 本发以板面定生死。
//
// §三 闸门（判题机同款 gcc 9.3 + `work/b18z_drv.cpp`，双尺度，本文件全文）：
// n=1e7 ⇒ FNV 6982038405141804923 / SUM 14324786532795827997 ✓ 与归档基线**逐位相同**
// n=1e8 ⇒ FNV 14022627980954509510 / SUM 8456426522220901567 ✓ 同样逐位相同（`unroll-loops` 不改语义 ✓）
// ================
#pragma GCC optimize("O2,unroll-loops")
#include <cstdint>
#include <cstring>
#include <cstdlib>
#include <sys/mman.h>
#include <immintrin.h>
#define LDU_M(p) ({ __m256i _v; __asm__("vmovdqu %1, %0":"=x"(_v):"m"(*(const __m256i*)(p))); _v; })
#define LDA_M(p) ({ __m256i _v; __asm__("vmovdqa %1, %0":"=x"(_v):"m"(*(const __m256i*)(p))); _v; })
#define STU_M(p, v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqu %1, %0":"+m"(*(__m256i*)(p)):"x"(_vv)); } while(0)
#define STA_M(p, v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqa %1, %0":"+m"(*(__m256i*)(p)):"x"(_vv)); } while(0)
#ifndef NO_AVX2
#pragma GCC push_options
#pragma GCC target("avx2,bmi2")
#endif
typedef uint8_t u8; /* [e8c-B] narrow side-array element type */
typedef uint16_t u16;
typedef uint32_t u32;
typedef uint64_t u64;
static const u64 M27 = 0x7FFFFFFull;
static const u64 LOW53 = 0x1FFFFFFFFFFFFFull; // idx:27 | y:26
static u64 *RA, *RB, *RC;
static unsigned char *FW, *FR, *FMem; // [b18z_fuse] fused 12 B (rec u64 | y u32) containers
static u32 *Y1, *Y2, *X1;
static u16 *U1;
static size_t gCap;
static const unsigned *gX, *gY;
static unsigned *gOut;
__attribute__((used)) alignas(64) static char m4e_pad[4160]; // [M4E] layout-only
static u32 gCnt[4096];
static u32 gH1[512]; // L1-digit histogram (y >> 17), built inside sort_x's pass 3
// from the value that pass is already moving to Yb; the
// fused attach_l1 scatter needs its prefix as the write
// cursor, so attach_l1 itself cannot build it
static u32 gNxt[1024]; // histogram of the NEXT radix pass, built during the current scatter
static u32 gBndL2[520]; // L2 bucket boundaries inside the current L1 bucket
static u32 gBndL1[520];
static size_t gOutLim = 0;
static u32 gCur[4096];
// ---- staging of the random out[] stores -------------------------------------------
// out[idx] = cnt with idx an essentially random permutation: 1e8 scattered 4-byte
// stores, each one a read-for-ownership fetch of a 64-byte line that is then written
// for 4 bytes. Measured on the judge at density 1 over 64 MB: 40.5 ticks per random
// store, versus 6.05 ticks when the same stores are issued bucket-by-bucket inside a
// 256 KB window. So: histogram idx >> IDXSH during the L1 pass (idx is already in the
// record there), scatter (idx << 27 | cnt) sequentially into a staging array, and replay
// staging into out[] one bucket at a time at the very end. The set of addresses written
// to out[] is exactly the same; only their order changes.
#define IDXSH 17
#define IDXNB (1024u)
static u32 gIdxCnt[IDXNB];
static u32 gIdxBase[IDXNB + 1];
static u32 gIdxCur[IDXNB];
#define PUT_STG(stg, idxv, cntv) do { u32 _i = (u32)(idxv); \
(stg)[gIdxCur[_i >> IDXSH]++] = ((u64)_i << 27) | (u32)(cntv); } while (0)
// ---- branchless 8x8x8 cumulative digit counter: query = 3 loads, insert = 3
// masked SIMD suffix adds over 8-lane u32; 2 KB of state, all L1 resident.
/* [e8y-bd] ONE base pointer for the whole 8x8x8 counter. The three digit tables
are laid out consecutively in gPB, so every address in cnt_q / cnt_ins is
`base + small displacement (+ index*scale)` off a single register. The old
form kept gP1/gP2/gP3 in three separate statics: the compiler had to hold
three base addresses (r11/r13/r15/rdx in the -O2 output) live across the whole
inner loop, and the spill traffic showed up as ~8 stack reloads per element.
Index map (all digit pairs collapse into the flat index -- a*8+b == d>>3 and
a*64+b*8+c == d for d < 512):
query : gPB[d>>6] + gPB[GP2 + ((d>>3)&63)] + gPB[GP3 + d]
insert: 8-lane groups start at gPB[GP3 + (d & 0x1F8)], gPB[GP2 + ((d&0x1C0)>>3)],
gPB[GP1]; the lane mask is gPB[GMK + 8*c] (32 B rows). */
#define GP1 0
#define GP2 8
#define GP3 72
#define GMK 584
static u32 gPB[8 + 64 + 512 + 64] __attribute__((aligned(32)));
static int gMskInit = 0;
static inline void cnt_init(void) {
if (!gMskInit) {
for (u32 m = 0; m < 8; m++)
for (u32 i = 0; i < 8; i++) gPB[GMK + m * 8 + i] = (i > m) ? 1u : 0u;
gMskInit = 1;
}
memset(gPB, 0, (GP3 + 512) * sizeof(u32));
}
static inline u32 cnt_q(u32 d) {
return gPB[d >> 6] + gPB[GP2 + ((d >> 3) & 63u)] + gPB[GP3 + (d & 511u)];
}
/* [SM] the mask operand is folded into the add: `vpaddd ymm, ymm, m256`.
VEX memory operands must be 32 B aligned; gMsk is alignas(32) and its rows are
32 B, so this is structurally safe (and it removes one instruction per level). */
#define ADDM_M(p, v) ({ __m256i _r; __asm__("vpaddd %1, %2, %0" : "=x"(_r) : \
"m"(*(const __m256i *)(p)), "x"(v)); _r; })
static inline void cnt_ins(u32 d) {
u32 *p3 = &gPB[GP3 + (d & 0x1F8u)];
u32 *p2 = &gPB[GP2 + ((d & 0x1C0u) >> 3)];
STU_M((__m256i *)p3,
ADDM_M((const __m256i *)&gPB[GMK + ((d & 7u) << 3)], LDU_M((const __m256i *)p3)));
STU_M((__m256i *)p2,
ADDM_M((const __m256i *)&gPB[GMK + (((d >> 3) & 7u) << 3)], LDU_M((const __m256i *)p2)));
STU_M((__m256i *)&gPB[GP1],
ADDM_M((const __m256i *)&gPB[GMK + ((d >> 6) << 3)], LDU_M((const __m256i *)&gPB[GP1])));
}
static void *bigalloc(size_t bytes) {
void *p = mmap(0, bytes, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
if (p != MAP_FAILED) return p; // mmap fails on JudgeDuck -> malloc
return malloc(bytes);
}
// ------------------------------------------------------------------
// 3-pass LSD radix sort by x of [spare:11|idx:27|x:27] records, permuting the
// parallel y array with the identical permutation; the result ends in X1 / Yb.
//
// ARM: the histogram of each pass is DERIVED DURING THE PREVIOUS PASS'S SCATTER,
// where the record is already in a register -- the three standalone histogram
// sweeps over the records are gone. Pass 1's histogram must be supplied by the
// caller, which has x[i] in hand while it is already building the record.
static void sort_x_fused(unsigned char *F0, unsigned char *F1, u64 *rec_out, u32 *y_out, size_t m) {
// [b18z_fuse] 3-pass LSD radix sort by x, with the y value carried INSIDE the record:
// each element's 12 B unit (u64 rec at +0, u32 y at +8) is one cache line's worth of
// scatter traffic instead of two. The scatter's concurrently-active lines drop from
// 512 records + 512 y (64 KB > 32 KB L1) to 512 (32 KB), and passes 1/2 no longer
// touch a separate y array at all. Pass 3 emits the consumer's containers unchanged
// (records -> rec_out, y -> y_out in this half's shifted y).
// ---- pass 1: gCnt must already hold the histogram of (rec & 511) ----
gCur[0] = 0;
for (u32 b = 1; b < 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
memset(gNxt, 0, sizeof(gNxt));
for (size_t i = 0; i < m; i++) {
const unsigned char *s = F0 + 12 * i;
u64 rec = *(const u64 *)s;
u32 yv = *(const u32 *)(s + 8);
u32 q = gCur[(u32)rec & 511u]++;
unsigned char *d = F1 + 12 * q;
*(u64 *)d = rec;
*(u32 *)(d + 8) = yv;
gNxt[(u32)(rec >> 9) & 1023u]++;
}
// ---- pass 2: ten-bit digit (x >> 9), build pass 3's eight-bit histogram ----
memcpy(gCnt, gNxt, sizeof(gNxt));
gCur[0] = 0;
for (u32 b = 1; b < 1024; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
memset(gNxt, 0, sizeof(gNxt));
for (size_t i = 0; i < m; i++) {
const unsigned char *s = F1 + 12 * i;
u64 rec = *(const u64 *)s;
u32 yv = *(const u32 *)(s + 8);
u32 q = gCur[(u32)(rec >> 9) & 1023u]++;
unsigned char *d = F0 + 12 * q;
*(u64 *)d = rec;
*(u32 *)(d + 8) = yv;
gNxt[(u32)(rec >> 19) & 255u]++;
}
// ---- pass 3: eight-bit digit (x >> 19) ----
memcpy(gCnt, gNxt, sizeof(gNxt));
gCur[0] = 0;
for (u32 b = 1; b < 256; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
/* pass 3 also builds the L1-digit histogram (y >> 17) of this half. A histogram only
needs the MULTISET of values, so counting here -- where the y value is already in
hand -- is bit-identical to counting it in a separate sweep over the compacted array. */
for (size_t i = 0; i < m; i++) {
const unsigned char *s = F0 + 12 * i;
u64 rec = *(const u64 *)s;
u32 yv = *(const u32 *)(s + 8);
u32 q = gCur[(u32)(rec >> 19) & 255u]++;
rec_out[q] = rec;
y_out[q] = yv;
gH1[(yv >> 17) & 511u]++;
}
}
// ======================= [n08a-att] fused consumer ==========================
/* `attach_l1`'s AOS path, specialised to one already-sorted bucket. `B` holds
`len` 12-byte records (rec u64 | y u32) in ascending x order; `abase` is the
bucket's absolute position in the half, so Xout[] indices stay identical.
Everything that used to be a loop-carried local of the whole-array sweep is
either passed in (xl/xrun: the monotone cross-term merge cursor) or reset per
bucket (prevx / the tie-group start / the rolling yhist lanes: an equal-x
group cannot straddle a bucket, since the bucket digit is the top of x). */
template<bool preseed>
static inline void attach_bucket(const u64 *B, size_t abase, size_t len, u32 xbase,
u32 *Xout, u64 *Dst, u8 *Cb,
const u32 *Xlow, size_t &xl) {
const u64 *it=B,*end=B+len;
u32 *xo=nullptr;
if constexpr(!preseed)xo=Xout+abase;
u32 ng=0;
__m128i prevx=_mm_set1_epi32(-1);
const __m128i xmask=_mm_set1_epi32(255);
__m256i yhist=_mm256_setzero_si256();
for(;it<end;++it){
const u64 rec=*it;
const __m128i xv=_mm_and_si128(_mm_cvtsi64_si128(rec),xmask);
const u32 same=(u32)_mm_cvtsi128_si32(_mm_cmpeq_epi32(xv,prevx));
ng=(ng+1u)&same;prevx=xv;
const u32 xvv=((u32)rec&255u)|xbase;
const u32 ia=(u32)((rec>>8)&M27);
const u32 ya=(u32)(rec>>35);
const __m256i yb=_mm256_set1_epi32((int)ya);
const unsigned hit=(unsigned)_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(yb,yhist)));
u32 dup=(u32)__builtin_popcount(_bzhi_u32(hit,ng));
if(__builtin_expect(ng>8,0)){
dup=0;
for(const u64 *p=it-ng;p<it;++p)dup+=((u32)(*p>>35)<ya);
}
yhist=_mm256_alignr_epi8(yhist,_mm256_permute2x128_si256(yb,yhist,0x20),12);
if constexpr(preseed){
const __m256i xv8=_mm256_set1_epi32((int)xvv);
const __m256i low8=LDU_M(Xlow+xl);
const unsigned hits=(unsigned)_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(xv8,low8)));
xl+=(unsigned)__builtin_ctz((~hits&255u)|256u);
while(Xlow[xl]<xvv)++xl;
}
u32 cnt=preseed?(u32)xl:0u;
cnt=(cnt-dup)&0x7ffffffu;
if constexpr(!preseed)*xo++=xvv;
const u32 d1=(ya>>17)&511u;
cnt=(cnt+cnt_q(d1))&0x7ffffffu;
cnt_ins(d1);
const u32 qq=gCur[d1]++;
Dst[qq]=((u64)cnt<<44)|(((u64)ya<<27)&0x00000fffffffffffull)|(u64)ia;
Cb[qq]=(u8)(cnt>>20);
}
}
/* scratch for the local sort's two bucket-local scatter buffers. Separate from
Y2 because Cb is (u8*)Y2 in both halves and the consumer now interleaves with
the sort. 4 MiB is SPLIT in two (pass-1b scratch and the compact pass-2 output), so a bucket
may hold ATT_GRAN/2/12 = 174762 records; the official n=1e8 gives ~24400 and the
guard below falls back to the old path if a half is skewed. */
#define ATT_GRAN (4u << 20)
static unsigned char gAttScr[ATT_GRAN] __attribute__((aligned(64)));
template<bool preseed>
static bool sort_x_local_att(unsigned char*F0,unsigned char*F1,u64*rec_out,u32*y_out,size_t m,
u32*Xout,u64*Dst,u8*Cb,const u32*Xlow,size_t lo){
static u32 bnd[2049],h0[256],h1[256],p0[256],p1[256];
u32 mx=0;for(unsigned b=0;b<2048;b++)if(gCnt[b]>mx)mx=gCnt[b];
if((size_t)mx*12>(ATT_GRAN/2)){
memset(gCnt,0,sizeof(gCnt));
for(size_t i=0;i<m;i++)gCnt[(u32)*(const u64*)(F0+12*i)&511u]++;
sort_x_fused(F0,F1,rec_out,y_out,m);return false;
}
bnd[0]=0;for(unsigned b=0;b<2048;b++){bnd[b+1]=bnd[b]+gCnt[b];gCur[b]=bnd[b];}
for(size_t i=0;i<m;i++){
const unsigned char*p=F0+12*i;u64 v=*(const u64*)p;u32 y=*(const u32*)(p+8);
unsigned char*q=F1+12*gCur[((u32)v>>16)&2047u]++;
*(u64*)q=v;*(u32*)(q+8)=y;
gH1[(y>>17)&511u]++; /* [n08a-att] moved here from pass 2: must be complete
before the first bucket is handed to the consumer */
}
/* the consumer's L1 scatter cursors, now that gH1 is final */
memcpy(gCnt, gH1, sizeof(gH1));
gCur[0] = 0;
for (u32 b = 1; b <= 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
for (u32 b = 0; b <= 512; b++) gBndL1[b] = gCur[b];
cnt_init();
size_t xl = 0;
if (preseed) ((u32 *)Xlow)[lo] = 0x7FFFFFFFu; /* merge sentinel (signed compare) */
unsigned char*temp=gAttScr,*temp2=gAttScr+(ATT_GRAN/2);
for(unsigned b=0;b<2048;b++){
u32 start=bnd[b],end=bnd[b+1],len=end-start;
if(!len)continue;
memset(h0,0,sizeof(h0));memset(h1,0,sizeof(h1));
for(u32 j=start;j<end;j++)h0[(u32)*(const u64*)(F1+12*j)&255u]++;
u32 ac=0;for(unsigned k=0;k<256;k++){p0[k]=ac;ac+=h0[k];}
for(u32 j=start;j<end;j++){
const unsigned char*p=F1+12*j;u64 v=*(const u64*)p;u32 y=*(const u32*)(p+8);
unsigned char*q=temp+12*p0[(u32)v&255u]++;
*(u64*)q=v;*(u32*)(q+8)=y;h1[((u32)v>>8)&255u]++;
}
/* [n08a-att] 0-based compact second-pass output: the sorted bucket goes to the
other half of the scratch instead of back into F0's slice, so (a) the consumer
reads it out of L2 (293 KB) rather than from a 600 MB array, and (b) F0 stops
being written at all -- which is what frees FW to serve as Dst below. */
ac=0;for(unsigned k=0;k<256;k++){p1[k]=ac;ac+=h1[k];}
u64 *packed=(u64*)temp2;
for(u32 j=0;j<len;j++){
const unsigned char*p=temp+12*j;u64 v=*(const u64*)p;u32 y=*(const u32*)(p+8);
const u32 mid=((u32)v>>8)&255u;
packed[p1[mid]++]=((u64)y<<35)|(((v>>27)&M27)<<8)|(v&255u);
}
u32 consumed=0;
for(u32 mid=0;mid<256;++mid){
const u32 stop=p1[mid];
if(stop>consumed)attach_bucket<preseed>(packed+consumed,start+consumed,stop-consumed,
(b<<16)|(mid<<8),Xout,Dst,Cb,Xlow,xl);
consumed=stop;
}
}
return true;
}
// ============================================================================
// ------------------------------------------------------------------
// FUSED attach + L1. y is already in sequence order (Yv); fold the x-tie correction
// into the 27-bit count, build the idx-bucket histogram for the staging pass, and
// scatter straight into the L1 layout -- the record never makes the round trip through
// X that attach()+levels() used to pay for. `preseed` = the count already holds the
// half's cross term (high half). `Dst`/`Cb` must be arrays this sweep does not read.
__attribute__((noinline)) static u32 dup_slow(const u32 *Yv, size_t p, size_t a, u32 ya) {
u32 d = 0;
for (size_t b = p; b < a; b++) d += (Yv[b] < ya);
return d;
}
template<bool preseed,bool AOS=false>
static void attach_l1(u64 *X, const u32 *Yv, u32 *Xout, u64 *Dst, u8 *Cb,
size_t m, const u32 *Xlow, size_t lo) {
// gCnt already holds the L1-digit histogram (built by the split sweep, copied in by
// the caller): its prefix IS the scatter cursor, so it must be final before the loop.
gCur[0] = 0;
for (u32 b = 1; b <= 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
for (u32 b = 0; b <= 512; b++) gBndL1[b] = gCur[b];
cnt_init();
/* cross term: #{low: x_j < x_i} by one monotone merge against the low half's
x-sorted array -- the SAME sweep the caller used to do separately, but here
the elements already arrive in exactly that order, so xl/xrun stay in
registers and the record never has to be patched. */
if (preseed) ((u32 *)Xlow)[lo] = 0x7FFFFFFFu; /* merge sentinel (signed compare) */
size_t xl = 0;
u32 xrun = 0;
/* The tie group is found by comparing with the PREVIOUS element instead of scanning
forward for the group's end: the array is x-sorted, so equal-x elements are
contiguous and `p` (the group start) only has to be updated when x changes. The
old form re-loaded X[a+1] for every element (one extra 8 B load per element over
1e8 elements, plus the outer while bookkeeping) purely to discover that the group
had length 1, which is the overwhelmingly common case for random x. */
size_t p = 0;
u32 prevx = 0xFFFFFFFFu;
__m256i yhist=_mm256_setzero_si256();
for (size_t a = 0; a < m; a++) {
{
u64 rec = AOS ? *(const u64*)((const unsigned char*)X+12*a) : X[a];
u32 xvv0 = (u32)rec & (u32)M27;
if (xvv0 != prevx) p = a;
prevx = xvv0;
u32 ia = (u32)((rec >> 27) & M27);
u32 xvv = (u32)rec & (u32)M27;
u32 ya = AOS ? *(const u32*)((const unsigned char*)X+12*a+8) : Yv[a];
u32 dup;
{
size_t ng = a - p;
if(AOS){
__m256i yb=_mm256_set1_epi32((int)ya);
unsigned hit=(unsigned)_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(yb,yhist)));
dup=(u32)__builtin_popcount(_bzhi_u32(hit,(unsigned)ng));
if(__builtin_expect(ng>8,0)){
dup=0;for(size_t j=p;j<a;j++)dup+=(*(const u32*)((const unsigned char*)X+12*j+8)<ya);
}
yhist=_mm256_blend_epi32(_mm256_permutevar8x32_epi32(yhist,_mm256_setr_epi32(7,0,1,2,3,4,5,6)),yb,1);
}else{
__m256i yv8 = LDU_M(Yv + p);
unsigned hit = (unsigned)_mm256_movemask_ps(_mm256_castsi256_ps(
_mm256_cmpgt_epi32(_mm256_set1_epi32((int)ya), yv8)));
dup = (u32)__builtin_popcount(_bzhi_u32(hit, (unsigned)ng));
}
}
if (preseed) {
const __m256i xv8 = _mm256_set1_epi32((int)xvv);
const __m256i low8 = LDU_M(Xlow + xl);
const unsigned hit = (unsigned)_mm256_movemask_ps(
_mm256_castsi256_ps(_mm256_cmpgt_epi32(xv8, low8)));
const unsigned adv = (unsigned)__builtin_ctz((~hit & 255u) | 256u);
xl += adv;
xrun += adv;
while (Xlow[xl] < xvv) { xl++; xrun++; }
}
u32 cnt = preseed ? xrun : 0u;
cnt = (cnt - dup) & 0x7FFFFFFu;
if (!preseed) Xout[a] = xvv;
u32 d1 = (ya >> 17) & 511u;
cnt = (cnt + cnt_q(d1)) & 0x7FFFFFFu;
cnt_ins(d1);
u32 qq = gCur[d1]++;
Dst[qq] = ((u64)cnt << 44) | (((u64)ya << 27) & 0x00000FFFFFFFFFFFull) | (u64)ia;
Cb[qq] = (u8)(cnt >> 20); /* [e8c-B] cnt < 2^27 => (cnt>>20) < 128 fits a byte */
}
}
}
// ------------------------------------------------------------------
// Two-pass MSD rank accumulation over the 26-bit y field (digits 9/9/8), with the
// last two digits resolved inside the L1 bucket (L1 cache resident) and out[] stored
// through the idx-bucket staging array. Replaces sol_v6's four-pass L1/L2/L3/L4 chain.
//
// L1 itself is done by attach_l1() -- see there. This routine starts from its output:
// incoming: X1[q] = [cntlo:20 @44..63 | d2:9 @35..43 | d3:8 @27..34 | idx:27 @0..26]
// Cb[q] = cnt >> 20
// working layout: same, with the count re-accumulated for digits d2 then d3
// side[q] = cnt >> 20 (low byte of the u16 side array)
// count <= n-1 < 2^27 => 20 + 8 bits is enough.
/* [e8y-bd] same one-base-pointer treatment as gPB above: 16x16 counter in one
array (a*16+b == d for d < 256). */
#define L3P1 0
#define L3P2 16
#define L3I 272
alignas(32) static u16 gLB[16 + 256 + 256];
static int gL3MaskInit = 0;
static inline void l3_init(void) {
if (!gL3MaskInit) {
for (u32 m = 0; m < 16; m++)
for (u32 i = 0; i < 16; i++) gLB[L3I + m * 16 + i] = (i > m) ? 1u : 0u;
gL3MaskInit = 1;
}
memset(gLB, 0, (L3P2 + 256) * sizeof(u16));
}
static inline u32 l3_q(u32 d) {
return (u32)gLB[L3P1 + (d >> 4)] + gLB[L3P2 + (d & 255u)];
}
/* [SM] same fold as cnt_ins: gL3Inc rows are 32 B and alignas(32). */
#define ADDW_M(p, v) ({ __m256i _r; __asm__("vpaddw %1, %2, %0" : "=x"(_r) : \
"m"(*(const __m256i *)(p)), "x"(v)); _r; })
static inline void l3_ins(u32 d) {
u16 *p2 = &gLB[L3P2 + (d & 0xF0u)];
STA_M((__m256i *)&gLB[L3P1],
ADDW_M((const __m256i *)&gLB[L3I + (d & 0xF0u)], LDA_M((const __m256i *)&gLB[L3P1])));
STA_M((__m256i *)p2,
ADDW_M((const __m256i *)&gLB[L3I + ((d & 15u) << 4)], LDA_M((const __m256i *)p2)));
}
static void levels(u64 *X0, u64 *X1, u16 *Ca, u8 *Cb, size_t m, unsigned *out, u64 *stg, u32 ylim) {
(void)m; // every bucket is delimited by gBndL1
// ARM: gCnt/gBndL1 were already consumed by attach_l1; only the idx-bucket prefix
// (whose histogram attach_l1 built) is still needed, and staging is the replay target.
gIdxBase[0] = 0;
for (u32 b = 0; b < IDXNB; b++) gIdxBase[b + 1] = gIdxBase[b] + gIdxCnt[b];
memcpy(gIdxCur, gIdxBase, sizeof(gIdxCur));
// ---------------- L2 + L3 fused, one L1 bucket at a time --------------
for (u32 b0 = 0; b0 < 512; b0++) {
size_t s = gBndL1[b0], e = gBndL1[b0 + 1];
if (e - s <= 1) {
if (e > s) PUT_STG(stg, (u32)(X1[s] & M27), ((u32)(X1[s] >> 44) | ((u32)Cb[s] << 20)) & 0x7FFFFFFu);
continue;
}
// Direct stable L2 buckets reuse Ca/U1, which the current engine otherwise
// leaves untouched after packing count high bits into X0 records.
size_t span = (size_t)ylim - ((size_t)b0 << 17);
if (span > (1u << 17)) span = 1u << 17;
const size_t active = (span + 255u) >> 8;
const size_t cap = (e - s + active - 1u) / active + 128u;
u64 *bins = (u64 *)Ca;
memset(gCnt, 0, 512 * sizeof(u32));
cnt_init();
for (size_t t = s; t < e; t++) {
u64 rec = X1[t];
u32 d2 = (u32)(rec >> 35) & 511u;
u32 inc = cnt_q(d2);
cnt_ins(d2);
const size_t q = (size_t)d2 * cap + gCnt[d2]++;
u32 cnt = (((u32)(rec >> 44) | ((u32)Cb[t] << 20)) + inc) & 0x7FFFFFFu;
// bins: index[0..26], final-y[27..34], contiguous count[35..61].
bins[q] = ((u64)cnt << 37) | (rec & 0x00000007FFFFFFFFull);
}
for (u32 d2 = 0; d2 < 512; d2++) {
const size_t nn = gCnt[d2];
if (!nn) continue;
const u64 *src = bins + (size_t)d2 * cap;
if (nn == 1) {
u64 rec = src[0];
PUT_STG(stg, (u32)(rec & M27), ((u32)(rec >> 37)));
} else if (nn <= 65535) {
l3_init();
for (size_t a = 0; a < nn; a++) {
u64 rec = src[a];
u32 d3 = (u32)(rec >> 27) & 255u;
u32 inc = l3_q(d3);
l3_ins(d3);
u32 cnt = (((u32)(rec >> 37)) + inc) & 0x7FFFFFFu;
PUT_STG(stg, (u32)(rec & M27), cnt);
}
} else {
cnt_init();
for (size_t a = 0; a < nn; a++) {
u64 rec = src[a];
u32 d3 = (u32)(rec >> 27) & 255u;
u32 inc = cnt_q(d3);
cnt_ins(d3);
u32 cnt = (((u32)(rec >> 37)) + inc) & 0x7FFFFFFu;
PUT_STG(stg, (u32)(rec & M27), cnt);
}
}
}
}
// ---------------- replay staging into out[], one idx bucket at a time ----
// 256 KB window per bucket: the scattered out[] stores become, per bucket, a
// sequential sweep of a window that stays in L1/L2.
for (u32 b = 0; b < IDXNB; b++) {
const u64 *p = stg + gIdxBase[b];
size_t nn = (size_t)(gIdxBase[b + 1] - gIdxBase[b]);
for (size_t t = 0; t < nn; t++) {
u64 v = p[t];
out[(size_t)(v >> 27)] = (u32)(v & M27); /* [e8c-D] 64-bit index: no zero-extend */
}
}
}
// ------------------------------------------------------------------
alignas(32) static u32 gSplitPerm[256][8];
static void split_perm_init() {
for (u32 mask=0;mask<256;mask++) {
u32 t=0;
for (u32 j=0;j<8;j++) if (mask & (1u<<j)) gSplitPerm[mask][t++]=j;
for (u32 j=0;j<8;j++) if (!(mask & (1u<<j))) gSplitPerm[mask][t++]=j;
}
}
template<bool upper>
static inline __attribute__((always_inline)) u32 split8(const unsigned *x,const unsigned *y,
size_t i,u32 K,unsigned char *fout) {
const __m256i xv=_mm256_loadu_si256((const __m256i*)(x+i));
const __m256i yv=_mm256_loadu_si256((const __m256i*)(y+i));
const __m256i kv=_mm256_set1_epi32((int)K);
const __m256i dy=upper ? _mm256_sub_epi32(yv,kv):yv;
const u32 mask=upper ? ((~(u32)_mm256_movemask_ps(_mm256_castsi256_ps(dy))) & 255u)
: (u32)_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(kv,yv)));
const __m256i perm=_mm256_load_si256((const __m256i*)gSplitPerm[mask]);
const __m256i xp=_mm256_permutevar8x32_epi32(xv,perm);
const __m256i yp=_mm256_permutevar8x32_epi32(dy,perm);
const __m256i ip=_mm256_add_epi32(_mm256_set1_epi32((int)i),perm);
const __m256i il=_mm256_cvtepu32_epi64(_mm256_castsi256_si128(ip));
const __m256i ih=_mm256_cvtepu32_epi64(_mm256_extracti128_si256(ip,1));
const __m256i xl=_mm256_cvtepu32_epi64(_mm256_castsi256_si128(xp));
const __m256i xh=_mm256_cvtepu32_epi64(_mm256_extracti128_si256(xp,1));
const __m256i rl=_mm256_or_si256(_mm256_slli_epi64(il,27),xl);
const __m256i rh=_mm256_or_si256(_mm256_slli_epi64(ih,27),xh);
const __m256i v0r=_mm256_permutevar8x32_epi32(rl,_mm256_setr_epi32(0,1,0,2,3,0,4,5));
const __m256i v0y=_mm256_permutevar8x32_epi32(yp,_mm256_setr_epi32(0,0,0,0,0,1,0,0));
const __m256i v0=_mm256_blend_epi32(v0r,v0y,0x24);
const __m256i v1l=_mm256_permutevar8x32_epi32(rl,_mm256_setr_epi32(0,6,7,0,0,0,0,0));
const __m256i v1h=_mm256_permutevar8x32_epi32(rh,_mm256_setr_epi32(0,0,0,0,0,1,0,2));
const __m256i v1y=_mm256_permutevar8x32_epi32(yp,_mm256_setr_epi32(2,0,0,3,0,0,4,0));
const __m256i v1=_mm256_blend_epi32(_mm256_blend_epi32(v1l,v1h,0xB0),v1y,0x49);
const __m256i v2r=_mm256_permutevar8x32_epi32(rh,_mm256_setr_epi32(3,0,4,5,0,6,7,0));
const __m256i v2y=_mm256_permutevar8x32_epi32(yp,_mm256_setr_epi32(0,5,0,0,6,0,0,7));
const __m256i v2=_mm256_blend_epi32(v2r,v2y,0x92);
_mm256_storeu_si256((__m256i*)(fout+0),v0);
_mm256_storeu_si256((__m256i*)(fout+32),v1);
_mm256_storeu_si256((__m256i*)(fout+64),v2);
return (u32)__builtin_popcount(mask);
}
__attribute__((aligned(4096))) void count_2d(int n, const unsigned *x, const unsigned *y, unsigned *out) { // [b18z-fusemv]
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;
}
size_t N = (size_t)n;
size_t M = N / 2 + 70000; // slack for the cut granularity
if (gCap < M) {
if (FW) munmap(FW, gCap * 12);
if (FMem) munmap(FMem, gCap * 16);
if (X1) munmap(X1, gCap * 4);
if (Y2) munmap(Y2, gCap * 2);
if (U1) munmap(U1, gCap * 2);
/* [b18z_fuse] footprint-neutral re-plumbing: the fused 12 B/pair array (FW) replaces
the old (records + y) input pair RA||Y1, and the sort's ping-pong partner (FR) is
the head of a single RC||RB region -- exactly the bytes pass 3 would overwrite with
the consumer's containers anyway. After the sort FW is dead and is reused as the
consumer's Dst (RA), so no extra memory is needed. M*8 is a multiple of 64, so RB
(and therefore every fused record) keeps line alignment. */
FW = (unsigned char *)bigalloc(M * 12);
FMem = (unsigned char *)bigalloc(M * 16 + 63); FMem = (unsigned char *)(((uintptr_t)FMem + 63) & ~(uintptr_t)63);
X1 = (u32 *)bigalloc(M * 4);
Y2 = (u32 *)bigalloc(M * 2);
U1 = (u16 *)bigalloc(M * 2);
gCap = M;
if (!FW || !FMem || !X1 || !Y2 || !U1) { gCap = 0; return; }
FR = FMem; RA = (u64 *)FW; RC = (u64 *)FMem; RB = (u64 *)(FMem + M * 8);
}
gX = x; gY = y; gOut = out; gOutLim = N;
// Random coordinates are uniform on [0,n); the midpoint cut gives
// two halves within the preallocated slack on the official instance.
// The first split below counts the exact low-half size as it writes.
u32 K = (u32)(((N / 2 + 65535u) >> 16) << 16);
if (K == 0) K = 1u << 16;
size_t lo = 0, hi = 0;
split_perm_init();
// ---- split A: low half -> RA (records) / Y1 (y, in the same order)
{
size_t a = 0;
memset(gCnt, 0, 2048 * sizeof(u32)); // pass-1 histogram for sort_x below
memset(gH1, 0, sizeof(gH1)); // L1-digit histogram for the fused sweep
/* ---- idx>>IDXSH histogram, counted HERE instead of in attach_l1.
attach_l1 had `gIdxCnt[ia >> IDXSH]++` inside its scatter, where ia is the
original element index and the scatter runs in x-sorted order, i.e. ia is a
random permutation: two memory ops per element on a random address. But this
loop walks i = 0..N-1 in order and i IS that index, so one register count per
2^IDXSH block does it with no extra traversal and no random access. The value
is bit-identical: attach_l1 used ia = (rec>>27) & M27 with the record's index
field set to i. The block check is hoisted out of the inner loop so the
per-element cost is a single register increment. */
memset(gIdxCnt, 0, sizeof(gIdxCnt));
for (size_t base = 0; base < N; base += ((size_t)1 << IDXSH)) {
size_t hi = base + ((size_t)1 << IDXSH);
if (hi > N) hi = N;
/* [e8c-E] blk_n is just this block's increment of the cursor `a' -- keeping it
as a second register accumulator cost one extra `add' per element across
both halves. Collapsing it to `a - a0' is exact and moves that add out of
the innermost loop (one per 2^IDXSH elements instead of one per element). */
const size_t a0 = a;
for(size_t sub=base;sub<hi;sub+=1024u) {
const size_t stop=(sub+1024u<hi)?sub+1024u:hi;
const size_t produced=a;
size_t i=sub;
for(;i+8<=stop;i+=8)a+=split8<false>(x,y,i,K,FW+12*a);
for(;i<stop;++i){
u32 yv=y[i],xv=x[i];u32 take=(u32)(yv<K);
*(u64*)(FW+12*a)=((u64)i<<27)|(u64)xv;*(u32*)(FW+12*a+8)=yv;a+=take;
}
#pragma GCC unroll 8
for(size_t j=produced;j<a;++j)
++gCnt[((u32)*(const u64*)(FW+12*j)>>16)&2047u];
}
gIdxCnt[(u32)(((size_t)base & M27) >> IDXSH)] = (u32)(a - a0);
}
/* the x pass-1 histogram is over the TAKEN set only, so it is built here from the
compacted array instead of inside the scan: the scan then touches the 2 KB table
0 times instead of once per element scanned, and this replacement pass runs over
lo ~= N/2 elements (1 load + 1 RMW each) instead of N (1 RMW each).
The y-side (gH1) histogram used to live in this sweep too; it is now built inside
sort_x's pass 3, where the y value is already being loaded (so that sweep pays
1 load + 1 RMW instead of 2 loads + 2 RMWs per element). */
lo = a;
hi = N - a;
}
// ---- low half: x sort (records -> RC, y -> RB; Y2 must stay free for the L1 side
// array, which the fused sweep below writes while it reads the y array)
/* [n08a-att] Dst must NOT be RC here: F1 = FR = RC = FMem are the SAME bytes, and
the fused consumer writes Dst while later buckets still read their F1 slice.
RB (FMem + M*8) is free until levels' staging, so Dst := RB and staging moves to
FW, which the consumer has finished reading by then. The fallback keeps the
original RC/FW split, so `lev` still has to know which path ran. */
bool aoslow=sort_x_local_att<false>(FW, FR, RC, (u32 *)RB, lo, X1, (u64 *)FW, (u8*)Y2, nullptr, 0);
if(!aoslow){ /* skewed half: the old two-step path */
memcpy(gCnt, gH1, sizeof(gH1));
attach_l1<false>(RC, (u32 *)RB, X1, (u64 *)FW, (u8 *)Y2, lo, (const u32 *)0, 0);
}
/* both paths now leave the consumer's output in FW (see the Dst note above) */
levels(RC, (u64 *)FW, U1, (u8 *)Y2, lo, out, RB, K);
// ---- split B: high half -> RB (records) / Y2 (y shifted into [0,2^26))
{
size_t b = 0;
memset(gCnt, 0, 2048 * sizeof(u32)); // pass-1 histogram for sort_x below
memset(gH1, 0, sizeof(gH1)); // L1-digit histogram for the fused sweep
/* ---- see split A: same per-block register accumulator, high half. */
memset(gIdxCnt, 0, sizeof(gIdxCnt));
for (size_t base = 0; base < N; base += ((size_t)1 << IDXSH)) {
size_t hi = base + ((size_t)1 << IDXSH);
if (hi > N) hi = N;
const size_t b0 = b; /* [e8c-E] see split A */
for(size_t sub=base;sub<hi;sub+=1024u) {
const size_t stop=(sub+1024u<hi)?sub+1024u:hi;
const size_t produced=b;
size_t i=sub;
for(;i+8<=stop;i+=8)b+=split8<true>(x,y,i,K,FW+12*b);
for(;i<stop;++i){
u32 yv=y[i],xv=x[i];u32 take=(u32)(yv>=K);
*(u64*)(FW+12*b)=((u64)i<<27)|(u64)xv;*(u32*)(FW+12*b+8)=yv-K;b+=take;
}
#pragma GCC unroll 8
for(size_t j=produced;j<b;++j)
++gCnt[((u32)*(const u64*)(FW+12*j)>>16)&2047u];
}
gIdxCnt[(u32)(((size_t)base & M27) >> IDXSH)] = (u32)(b - b0);
}
/* see split A: same x pass-1 histogram; the y-side histogram is built in sort_x. */
}
// ---- high half: x sort (records -> RC, y -> RA reused as a u32 buffer)
bool aoshigh=sort_x_local_att<true>(FW, FR, RC, (u32 *)RB, hi, nullptr, (u64 *)FW, (u8*)Y2, X1, lo);
if(!aoshigh){
memcpy(gCnt, gH1, sizeof(gH1));
attach_l1<true>(RC, (u32 *)RB, (u32 *)0, (u64 *)FW, (u8 *)Y2, hi, X1, lo);
}
levels(RC, (u64 *)FW, U1, (u8 *)Y2, hi, out, RB, (u32)N - K);
}