// References:
// - Duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/109709: copied its public 1008e8 fused-sort/count engine with its original citations retained below. No separate license notice appears; use is credited under duck.ac's public-code terms.
// - Duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/109648: retained our AVX2 split pack; https://duck.ac/submission/109844: our 9/10/8 radix repartition; https://duck.ac/submission/109849: our 64-byte FMem alignment. These are our own code and ideas.
// Approach:
// - Keep the faster 9/10/8 radix and 64-byte FMem alignment from our prior submissions, then shift FMem by another 128 bytes to change the relative cache-set mapping of the two fused 12-byte scatter buffers while preserving alignment and output order.
// Purpose:
// - Experimentally price this distinct cache-set placement on the official judge after its 15-million-point Try It result showed a small improvement; verify exact output and timing.
// ===== REFERENCES =====
// [1] duck.ac 用户 **saffah_codex_6s_agg2**，提交 **#109648** <https://duck.ac/submission/109648>
// （3497.918345 ms）—— 本文件正文的来源。本发正文 = 其 #109648 正文**逐字节复制**，
// 仅改**一处编译期摆放/代码生成旋钮**（见下），输出语义零变化（逐位闸门复核）。存档 `ref2/rival_109648.cpp`。
// 合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」⇒ 站上提交按站点规则公开可见；
// 其正文无独立许可证声明，本发按站点规则完整署名取用。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号）：**#109669** <https://duck.ac/submission/109669>
// （3498.16186 ms，控制件）· **#109694** <https://duck.ac/submission/109694>（3473.797083，本席 H1 隔离族最好件）·
// **#109431** <https://duck.ac/submission/109431>（3536.651047，其 #109648 的父件）。
// ======================
// ===== 思路 =====
// 【本发 = 纯摆放/代码生成臂（`§2.19.1181/1201/1211`）—— 正文与 #109669 逐字节相同，只改一处旋钮。】
//
// 已定价的两条（本席，判题机）：
// · `#109669`（控制，`align-functions=4096` + sections + `count_2d aligned(4096)`）= **3498.161860**
// · `#109694`（**H1 隔离**：`align-functions=4096` → `64`）= **3473.797083** ⇒ **−24.365 ms**
// ★ 即 `x14f` 的 `H1`（在 G3 上加 `align-functions=4096`）在**非融合基座**上只值 +10.14 ms，
// 而在**融合件**上值 **+24.4 ms** ⇒ **摆放税随二进制体量放大**（可入库）。
// H1 族四臂（af64 / af64+align-loops / af256 / af64+count_2d128）读数 3473.80~3475.23
// ⇒ **彼此相差仅 1.4 ms ⇒ 已进入平坦好带**，需换轴。
//
// 本发单变量：`去掉 `unroll-loops`（保留 `O2` + `align-functions=64`）`。其余一字不动。
// **逐位闸门**（`ref/gcc9/g9.sh -O2 -std=c++17 -Dcount_2d=count_2d_arm` + `work/b18z_drv.cpp`）：
// n=1e7 ⇒ `FNV 6982038405141804923 / SUM 14324786532795827997` = 与存档基线**逐位相同** ✓
// 目的：判题机定价（试验性）。判据 = 严支 `0.99*T+1µs = 3462.940162`（better 现为 #109694 = 3473.797083，需 −10.856921 ms）。
// ================
// References:
// - Duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/109431: copied its accepted 12-byte fused-record engine and eight-point SIMD split; no separate license notice appeared in public code.
// - Our account saffah_codex_6s_agg2, https://duck.ac/submission/109645: retained our three contiguous AVX2 stores that replace split8's stack and scalar stores; our change has no external license.
// Approach:
// Form each permuted original index as block_start + the already loaded 0..7 permutation entry, eliminating a redundant vector permutation per eight input points.
// Purpose:
// Test whether the reduced split instruction count closes the remaining official 1008e8 timing gap.
// References:
// - Duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/109431: copied the accepted 12-byte fused-record sorter and eight-point SIMD split as the base; public source carried no separate license notice.
// Approach:
// Pack eight selected 12-byte records entirely in AVX2 registers and write three contiguous 32-byte vectors, removing the eight-record stack buffer and sixteen scalar record stores per split block.
// Purpose:
// Experimentally test whether cheaper split writes improve the official 1008e8 runtime without changing the record layout or output.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #109143 <https://duck.ac/submission/109143>
// （3639.883380 ms）—— 本文件正文的**来源**。本发正文 = 他的 #109143 正文**逐字节复制**
// （剥掉头部注释后 diff = 0 处；正文另存为 `work/g1_body_109143.cpp`）。
// 他在正文上方的 `#pragma GCC optimize("O2,unroll-loops")` 一并原样保留。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #109061 <https://duck.ac/submission/109061>
// （3699.969619 ms）—— **他的 #109143 就是在本件之上改出来的**（逐函数 md5：`SAME=10 / DIFF=1`，
// 只有 `count_2d` 不同，且剥头 diff 只有 **3 个 hunk**：新增 `gSplitPerm[256][8]` +
// `split_perm_init()` + `split8<upper>()`，以及两处 split 调用点改成 8 路调用）。
// 也就是说：他这一发 = 我们 #109061 的刀 B/D/E **原样** + 他自己的 `split8`。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #109123 <https://duck.ac/submission/109123>
// （3687.914900 ms）= #109061 + 我们的刀 K/L。★ **他的件里没有 K/L**：那两刀落在**标量** split
// 循环里，而他把该循环换成了 8 路 SIMD（标量只余每个 2^IDXSH 块 ≤7 个元素的收尾）
// ⇒ K/L 在新基座上已无载体（收尾走 ≤5341 次），故本发**不叠 K/L**，保持他的正文逐字节。
// [4] 正文内部保留的完整上游引用链（逐 hunk 继承、未删改）：
// #108844 → #104989 → #104945 → #104764 → #104366 / #104255 / #103599 / #103504 / #103254 / #102212 …
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」，即站上提交按站点规则公开；
// ===== 思路 =====
// 【本发 = `b18z` 融合刀（y 融进记录缓存行，见下）+ `x14f_` 席的**代码摆放 mover**】
//
// 进入时口径：`ok=0 RED` · `mine = 3626.628876 (#109409)` · `T = 3639.883380 (#109143)`
// · 严支 `0.99*T+1µs = 3603.485546` ⇒ 需再快 **23.143330 ms（0.638%）**（本发判据）。
//
// 【为什么要重测：`x14f_` 席的 12 发摆放带**更正了本席此前的结论**】
// 实测带 = **[−13.26, +32.53] ms（全幅 45.8 ms = 1.26%）**；其中**代码**摆放臂 6/7 落在同一条好带
// （−11~−13 ms，互不叠加 ⇒ 现役件原在坏位，被挪出 ⇒ **已兑现 −12.87 ms**，即 #109409 的 3626.629）；
// **数据**摆放臂 2/2 = **+31~33 ms**（热数组页对齐就贵）。
// ⇒ ★ 任何**动数组/容器**的结构性件（本融合刀正是）自带 **±13~33 ms 摆放税** ⇒ 本席 #109350 的
// "+0.252 ms 中性"**不足以判"无效果"**（它可能是"真 −32 ms 收益 + 32 ms 摆放税"）⇒ 本发按新规程重测。
//
// 【本发做的两件事】
// 1. **语义件（不变）**：`(X,Y)` 融成单一 **12 B 记录**（`u64 rec @+0 | u32 y @+8`）：y 的散写与记录同一条
// 缓存行 ⇒ 散写活跃行 **1024 → 512**（= 32 KB = 恰好 L1）；趟 3 仍把记录写给 RC、y 写给 (u32*)RB
// ⇒ 消费者 `attach_l1`/`levels` **一行未改**；内存近似中性（+63 MB：FW M*12 取代 RA+Y1，
// 乒乓伙伴 FR = 单块 RC‖RB 区的头 12*M 字节，M*8 是 64 的倍数 ⇒ 行对齐）。
// 2. **摆放 mover（本发新增，语义零变化）**：完全照 `x14f-G3/H1` 的配方 ——
// `#pragma GCC optimize("O2,unroll-loops,align-functions=4096")` +
// `sort_x_fused`/`levels` 加 `__attribute__((section(".text.b18z")))` +
// `count_2d` 加 `__attribute__((aligned(4096)))`。
//
// 【闸门（提交前已跑过）】判题机工具链 `ref/gcc9 -O2 -std=c++17`：`-c` rc=0、链接 rc=0 ✓
// · **逐位闸门**（`work/b18z_drv.cpp`，A/B 各自独立进程、同驱动、FNV64 over out[] + 加权和）：
// n=10000000 seed=88172645463325252 ⇒ FNV **6982038405141804923** / SUM **14324786532795827997**
// ⇒ 与现役/未 mover 件**逐位相同** ✓（摆放只动地址，不动语义 ✓）
// · 早期同语义件（无 mover）已实测：n=1e8 FNV 14022627980954509510 与现役件逐位相同 ✓
// ================
#pragma GCC optimize("O2,align-functions=64") // [e8y]
#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.
static u32 gP1[8];
static u32 gP2[8][8];
static u32 gP3[8][8][8];
static u32 gMsk[8][8] __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++) gMsk[m][i] = (i > m) ? 1u : 0u;
gMskInit = 1;
}
memset(gP1, 0, sizeof(gP1));
memset(gP2, 0, sizeof(gP2));
memset(gP3, 0, sizeof(gP3));
}
static inline u32 cnt_q(u32 d) {
u32 a = d >> 6, b = (d >> 3) & 7, c = d & 7;
return gP1[a] + gP2[a][b] + gP3[a][b][c];
}
/* [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 a = d >> 6, b = (d >> 3) & 7, c = d & 7;
STU_M((__m256i *)&gP3[a][b][0],
ADDM_M((const __m256i *)gMsk[c], LDU_M((const __m256i *)&gP3[a][b][0])));
STU_M((__m256i *)&gP2[a][0],
ADDM_M((const __m256i *)gMsk[b], LDU_M((const __m256i *)&gP2[a][0])));
STU_M((__m256i *)&gP1[0],
ADDM_M((const __m256i *)gMsk[a], LDU_M((const __m256i *)&gP1[0])));
}
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.
__attribute__((section(".text.b18z"))) 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]++;
}
}
// ------------------------------------------------------------------
// 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>
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;
for (size_t a = 0; a < m; a++) {
{
u64 rec = 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 = Yv[a];
u32 dup;
{
size_t ng = a - p;
__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.
alignas(32) static u16 gL3P1[16];
alignas(32) static u16 gL3P2[16][16];
alignas(32) static u16 gL3Inc[16][16];
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++) gL3Inc[m][i] = (i > m) ? 1u : 0u;
gL3MaskInit = 1;
}
memset(gL3P1, 0, sizeof(gL3P1));
memset(gL3P2, 0, sizeof(gL3P2));
}
static inline u32 l3_q(u32 d) {
return (u32)gL3P1[d >> 4] + gL3P2[d >> 4][d & 15u];
}
/* [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) {
u32 a = d >> 4, b = d & 15u;
STA_M((__m256i *)&gL3P1[0],
ADDW_M((const __m256i *)&gL3Inc[a][0], LDA_M((const __m256i *)&gL3P1[0])));
STA_M((__m256i *)&gL3P2[a][0],
ADDW_M((const __m256i *)&gL3Inc[b][0], LDA_M((const __m256i *)&gL3P2[a][0])));
}
__attribute__((section(".text.b18z"))) 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[q] = ((u64)cnt << 44) | ((u64)((cnt >> 20) & 0x7Fu) << 35) |
(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 >> 44) | (((u32)(rec >> 35) & 0x7Fu) << 20)));
} 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 >> 44) | (((u32)(rec >> 35) & 0x7Fu) << 20)) + 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 >> 44) | (((u32)(rec >> 35) & 0x7Fu) << 20)) + 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 + 128); FMem = (unsigned char *)((((uintptr_t)FMem + 63) & ~(uintptr_t)63) + 128);
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, 512 * 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;
size_t i=base;
for (; i+8<=hi; i+=8) a += split8<false>(x,y,i,K,FW+12*a);
for (; i<hi; 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;
}
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). */
for (size_t j = 0; j < a; j++) {
gCnt[(u32)(*(const u64*)(FW+12*j)) & 511u]++;
}
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)
sort_x_fused(FW, FR, RC, (u32 *)RB, lo);
// the low half's sorted x values are emitted by attach_l1 itself (Y1 is free after
// sort_x's pass 3); the standalone copy sweep is gone
memcpy(gCnt, gH1, sizeof(gH1));
attach_l1<false>(RC, (u32 *)RB, X1, (u64 *)FW, (u8 *)Y2, lo, (const u32 *)0, 0);
levels(RC, (u64 *)FW, U1, (u8 *)Y2, lo, out, RB, K); // low half -> staging in RB -> out
// ---- split B: high half -> RB (records) / Y2 (y shifted into [0,2^26))
{
size_t b = 0;
memset(gCnt, 0, 512 * 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 */
size_t i=base;
for (; i+8<=hi; i+=8) b += split8<true>(x,y,i,K,FW+12*b);
for (; i<hi; 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;
}
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. */
for (size_t j = 0; j < b; j++) {
gCnt[(u32)(*(const u64*)(FW+12*j)) & 511u]++;
}
}
// ---- high half: x sort (records -> RC, y -> RA reused as a u32 buffer)
sort_x_fused(FW, FR, RC, (u32 *)RB, hi);
memcpy(gCnt, gH1, sizeof(gH1));
// the fused sweep reads Yv = (u32*)RA and X = RC, so the side array can be Y2 (the
// low half's Cb, dead now) and staging RA (dead after attach_l1), exactly as before
attach_l1<true>(RC, (u32 *)RB, (u32 *)0, (u64 *)FW, (u8 *)Y2, hi, X1, lo);
levels(RC, (u64 *)FW, U1, (u8 *)Y2, hi, out, (u64 *)RB, (u32)N - K); // high half -> staging in RA -> out
}