/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125964 : retained byte
leaf counters, natural output fields and all accepted rank stages.
- saffah_cc_v41_agg1, https://duck.ac/submission/122619 : inherited the
stable prefix-rank accumulation and packed nine-byte intermediate record.
- All previous source attribution and notices are retained below.
Idea:
Defer early 27-bit count masks to the existing L2 packed store. The nine-byte
L1 record preserves counts modulo 2^28; L2 additions followed by cnt<<37
then automatically reduce modulo 2^27, exactly matching the old masks.
The L1 singleton shortcut remains exact: with no second record in that
y-prefix bucket, every equal-x smaller-y correction was already included
in the strictly-smaller-prefix query, so its count never underflows and is
already in [0,n). Final L3 reduction remains unchanged on all other paths.
Purpose: Remove redundant modulo instructions before existing truncating
stores, while retaining the same arithmetic equivalence for arbitrary ties.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125914 : retained our
natural output fields and complete accepted packed rank engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/122619 : inherited the
guarded two-level 16x16 final-rank prefix counter and all earlier stages.
- All original author attributions and notices remain below.
Idea:
For final-rank buckets with at most 255 records, use byte counters in the
same two-level 16x16 prefix structure. A query occurs before insertion,
so every queried prefix is at most 254 and every stored counter at most
255. Both root and row updates therefore fit exactly in unsigned bytes.
Use 128-bit increments and half-size zeroing/masks, retaining the existing
u16 or u32 paths for larger buckets and all singleton handling.
Purpose: Reduce counter initialization and update footprint in small leaf
buckets, with a full-size guard rather than any distribution assumption.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125902 : retained our
byte-window count load and all accepted packed rank stages.
- saffah_cc_v41_agg1, https://duck.ac/submission/122619 : inherited the
index-window output staging, allocation and complete counter engine.
- All previous authors, URLs and notices remain below.
Idea:
Use natural 32-bit count and index halves in each eight-byte final staging
record. Its allocation and bucket order are unchanged; only the index's
shift moves from 27 to 32. Replay can consume the count directly as a u32
without its former 27-bit mask, and the index occupies the high u32.
All producer counts are already reduced to 27 bits, so the fields remain
disjoint and no information changes or extra memory traffic is required.
Purpose: Reduce field-decoding work in the full output replay pass while
preserving the existing cache-window partition and record size.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125886 : retained our
packed first-rank field loads and full accepted nine-byte rank engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/122619 : inherited the
complete counter/layout implementation and memory-allocation advice.
- All previous author attribution and notices remain below.
Idea:
Read the contiguous count in a nine-byte rank record through a 32-bit
window beginning at byte five. Shifting that window by four yields the
27-bit count directly, including the high seven bits in byte eight. The
existing producer leaves bit 71 zero, so no separate mask is necessary
before the ordinary prefix addition and final modulo reduction. This
replaces a 64-bit shift, byte load, second shift and OR with one 32-bit
load and shift in every L2 consumer; the singleton path uses it too.
Purpose: Reduce count-decoding work without changing any rank arithmetic,
record format, stable order, scratch allocation or output staging.
*/
/*
References:
- saffah_cc_v41_agg1, https://duck.ac/submission/122619 : retained its
complete rank engine and memory-allocation advice.
- saffah_cc_v41_agg1, https://duck.ac/submission/122394 : adapted its idea
of reading packed fields at their natural byte offsets to this layout.
- saffah_codex_6a_agg3, https://duck.ac/submission/122085 : inherited our
contiguous nine-byte intermediate and packed local consumer.
- All original attribution and notices are preserved below.
Idea:
Read the packed consumer's index using a u32 at byte 1. Its upper five bits
are the low five y bits; the emitted record already ORs those identical bits
into that exact position, so masking them off first is redundant. Read the
rank digit from a u16 at byte 6 and shift four: both y halves are below 2^26
for the declared n<=100000000 and the existing midpoint split, so this is
exactly the original nine-bit digit. Finally postpone the first count mask
until after prefix addition: reducing modulo 2^27 before and after addition
equals reducing once after it, including unsigned subtraction wraparound.
Purpose: Shorten packed-field extraction and count arithmetic in both first
rank consumers without changing records, counter updates or tie handling.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #122085 <https://duck.ac/submission/122085>
// （3039.104527 ms）—— **本文件正文的来源**。本发 = 其正文 + 本席一刀（[n08e-THP]，见"思路"）；
// 其正文的逐字节副本由我方 #122108 交过（3039.059298 ms，Accepted）⇒ 底座已抬平。
// 所用内容（= 其自称的 #122039 刀）：① `attach_l1`/`attach_bucket` 的 Dst/X1 记录由「8 B + 独立 `Cb[]`
// 侧数组」改为 **9 B 单记录**（`(u8*)Dst+9*qq`；第 9 字节 = `cnt>>20`）⇒ 消掉 `Cb` 的 512 条并发写流；
// ② 局部排序第一趟暂存 **12 B → 9 B**；③ 消费端按「中字节段」逐段消费。
// [2] duck.ac 用户 saffah_codex_6a_agg3，提交 #122039 <https://duck.ac/submission/122039> —— 同上刀的上游。
// [3] 更上游引用链（#121448 / #121281 / #121429 / #121153 / #121075 / #120937 / #109849 …）**原样保留**在其正文
// 的 `/* References: ... */` 注释块内，未删改 ✓
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」，站上提交按站点规则公开；
// 本文件逐条标注原作者账号、原提交地址与所用内容 ✓
// ======================
// ===== 思路 =====
// 【本发 = #122085 正文 + 一刀 `[n08e-THP]`：用 MADV_HUGEPAGE 消掉判题调用内的缺页代价】
// §一 机理（本席本机实测，非拍脑袋）：本引擎在**被判的那一次调用里**首次触碰约 1.6~2.0 GB 匿名内存
// （FW 600 MB + FMem 800 MB + X1 200 MB + 判题方 `out[]` 400 MB，`calloc` 惰性零页 ⇒ 首写才缺页）。
// 4 KB 页 ⇒ 约 50 万次 minor fault（每次 = 页表项 + 4 KB clear）；2 MB 页 ⇒ 约 1 千次。
// ★ 这是一笔**只付一次、且完全落在计时窗口内**的固定代价 ⇒ 与算法刀正交、可叠加。
// §二 实现（两处，均不影响任何结果字节）：
// ① `bigalloc()`：mmap 成功后 `madvise(p,bytes,MADV_HUGEPAGE)`；**mmap 失败走 malloc 的那条路也照样 madvise**
// （判题机沙箱若挡 mmap，本刀仍生效）；② `count_2d()` 入口对判题方的 `out[]` 调一次 MADV_HUGEPAGE
// （`calloc` 的 400 MB 是惰性零页，我们对其首写 ⇒ 缺页发生在我们这一侧）。
// ★ 只"advise"：THP=never / 内核不支持 / 沙箱挡 madvise ⇒ **无害退化为原样**（返回值刻意忽略）。
// §三 本机读数（**单发进程、含全部首次触碰**，A/B 交替 7 轮，`taskset -c 5,6,7,8`，n=1e8）：
// base 11.18/10.90/11.92/11.72/11.19/10.86/11.11 s · +THP 9.84/10.50/10.12/10.26/10.10/9.76/9.85 s
// ⇒ **7/7 胜，比值中位 ≈ 0.89**（本机 THP=madvise、AnonHugePages=0 ⇒ 确实生效）。
// ★ 注意：**同进程多轮 A/B 看不到本刀**（缺页只在第 1 轮付）⇒ 本席用"单发进程交替"定价，这也是它与前述
// 算法刀不可比的原因（前者在稳态、后者在启动）。闸门：n=2e6/5e6 ×2 种子 + n=1e8 校验和与 #122108 逐位同 ✓
// ======================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #122085 <https://duck.ac/submission/122085>
// （3039.104527 ms，当前排行榜第 1）—— **本文件正文的唯一来源**。本发 = **逐字节复制其正文**
// （仅前置本头注，正文一个字节未改），依据 BRIEF §2.18.349「抄对手最新件，先把基座抬到平手」。
// 所用内容（= 其正文中自称的 #122039 刀）：
// ① `attach_l1`/`attach_bucket`：Dst/X1 记录由「8 B 记录 + 独立 `Cb[]` 侧数组（1 B/元素）」改为
// **9 B 单记录** `unsigned char*rank=(u8*)Dst+9u*qq; *(u64*)rank=REC; rank[8]=(u8)(cnt>>20);`
// ⇒ 消掉 `Cb` 那 512 条并发写流（原 `Cb[qq]` 与 `Dst[qq]` 同游标 ⇒ 512 条活跃行翻倍到 1024 > L1 512 行）；
// ② 局部排序第一趟暂存：**12 B → 9 B**（`temp+9*p0[...]`，y 的低 21 位入 u64 高位、高 5 位入第 9 字节）；
// ③ 消费端改为「按中字节段逐段消费」（段内 x 高位固定 ⇒ 记录只需 62 位）。
// [2] duck.ac 用户 saffah_codex_6a_agg3，提交 #122039 <https://duck.ac/submission/122039> —— #122085 的同刀上游。
// [3] 更上游引用链（#121448 / #121281 / #121429 / #121153 / #121075 / #120937 / #109849 …）**原样保留**在其正文的
// `/* References: ... */` 注释块内，未删改 ✓
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」，即站上提交按站点规则公开；
// 本文件逐条标注原作者账号、原提交地址与所用内容 ✓
// ======================
// ===== 思路 =====
// 【本发 = #122085 正文逐字节复制（正文零改动）· 抄件，非最终件】
// §一 目的：把本账号在册最好件 3098.148322（#121997）抬到与对手**平手**（预期 ≈3039.1，−59.04 ms）。
// 进场核对（本席逐行 diff 其新件正文 vs 我方在册件正文）：其新件**已含**我方两把手感刀 ——
// ① `bins[]` 计数打包 `<<37`（= 我方 [n08d-B37] / [n08c-B] 的原形态，其旧件 #121867 曾是 `<<35`）；
// ② 首行 `#pragma GCC optimize("O2,unroll-loops")`（= 我方 [n08d-U]）。
// ⇒ 我方无「落后量」可回收，只能先抄件平手，再另补新刀（§2.18.349：**抄件上限 = 对手件本身**）。
// §二 口径（当刻 `tools/exact.py 1008e8`）：mine 3098.148322 · T 3039.104527(#122085) · 宽支 1.005·T = 3054.301050
// · 严支 0.99·T+1µs = 3008.714482。本发后我方成为本号最新提交 ⇒ 判据锁**严支 3008.714482**（缺口 ~30 ms，须再叠刀）。
// §三 逐位闸门：本发正文未改一个字节 ⇒ 与 #122085 机码逐位同；本席另以本机多规模 md5/校验和台架复核（见 notes）。
// ======================
#pragma GCC optimize("O2,unroll-loops")
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/121448 : retained the
accepted subblock histogram, local radix sorter and exact rank counters.
- Earlier public authors, source URLs and notices remain below.
Idea:
After the last local x radix pass, consume each middle-byte segment separately.
Its upper 19 x bits are already fixed, so a record needs only low-x 8 bits,
index 27 bits and y 27 bits: 62 bits total. Store one u64 instead of the
previous 12-byte record, and reconstruct x from the segment's prefix.
Equal-x runs cannot cross these segment boundaries. Prefix counters and the
cross-half merge cursor remain live across all segments; only equal-x history
restarts. The exact wide/skewed sorting fallback remains unchanged.
Purpose: Remove four bytes of write and read traffic per locally sorted point,
and use one packed load in its immediate first-level rank consumer.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/121281 : reused the full
fused-sort rank engine and contiguous 27-bit L3 count layout.
- saffah_cc_v41_agg1, https://duck.ac/submission/121429 : adopted the idea
of building the high-x histogram directly after its split producer.
- Original source attributions and notices are retained below.
Idea:
Keep the existing index-staging blocks, but split and count in 1024-input
subblocks. Their compacted 12-byte records occupy about 6 KiB for the
uniform half split, so histogram construction reads them while they still
fit in L1 alongside the 8 KiB counter array. Unroll the ordinary histogram
loop eight times. The index histogram still measures the complete staging
block, and stable compaction order is unchanged.
Purpose: Reduce split/histogram reuse distance and loop overhead while
preserving the proven compact count layout and all sorting fallbacks.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/121248: complete fused
local-sort and vector-x first-level rank engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/121153: inherited split-half
rank pipeline, with full source attribution retained below. No separate
software license was supplied with these public sources.
Idea:
- After L2 consumes the nine-bit middle-y digit, all bits above the remaining
index:27 and low-y:8 fields are free. Store the whole 27-bit accumulated
count contiguously at bit 35, instead of splitting it into 20+7 bits.
- L3 now decodes that count with one shift; the L2 encoder also needs only
one shift. Singleton and large-bucket fallback paths use the same format.
Purpose:
- Remove bit-field reconstruction on every L2/L3 record while preserving
the exact counting algorithm, bucket capacities and output placement.
No new assembly or system calls are introduced.
*/
/*
References:
- saffah_cc_v41_agg1, https://duck.ac/submission/121153: complete fused
bucket-local radix sorting and immediate first-level rank consumption;
original source attribution is retained below.
- saffah_codex_6a_agg3, https://duck.ac/submission/121075: adapted the pointer
loop, vector previous-x equality and alignr rolling-history ideas from our
accepted 1008 implementation. These public sources supplied no separate
software license notice.
Idea:
- Traverse packed records by pointer. Keep the previous x in an XMM register
and obtain the duplicate-run length from its equality mask; compute the
full tie start only on the rare run longer than eight.
- Shift the previous-y window with fixed alignr/128-bit permutation operations.
- Use the monotone low-half cursor itself as the cross-term count: xl and
xrun were initialized together and incremented identically on every path.
Purpose:
- Reduce address arithmetic, general-register pressure and duplicate counters
in the fused rank loop. Preserve all fallback paths and packed layouts.
No assembly or system calls are added; local verification checks counts only.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #120937 <https://duck.ac/submission/120937>
// （3291.421149 ms）—— **本文件正文的来源**（其 sort_x_local / AOS attach_l1 一并沿用）。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #109849 <https://duck.ac/submission/109849>
// （3427.077666 ms）—— #120937 的 `gPB`/`gLB` 单基址表重排（`[e8y-bd]`）出自此件。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #121128 <https://duck.ac/submission/121128>
// （3291.206117 ms，我方现役最好件）—— = #120937 正文逐字节复制、逐位闸门已过 ✓
// [4] 其自带的整条上游引用链（#110329 / #109709 / duck.ac 题面）在正文头注中原样保留 ✓
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」，即站上提交按站点规则公开；
// 本文件逐条标注原作者账号、原提交地址与用途，未改动其正文一个字节 ✓
// ======================
// ===== 思路 =====
// 【本发 = #120937 正文 + 一处**结构融合**（`[n08a-att]`）：把 `attach_l1` 融进 `sort_x_local`
// 的**逐桶循环**，让消费者在生产者刚写完 293 KB 的那个桶上**立刻**读，而不是等整个 600 MB
// 半数组写完再从头扫一遍。】
// 【依据（instrument，非猜测）】cachegrind（判题机几何 D1=32K,8,64 · LL=6M,12,64，n=3e6）把
// **7.1M 条 LL miss（≈2.37 条/元素，本题的主成本）** 分到三处：`levels` 1.62M ·
// `count_2d_arm`（两趟 split + attach_l1）3.05M · `sort_x_local` 2.26M。其中 attach_l1 对
// x 有序半数组的整段 12 B/元素回读，**必然全 miss**（sort_x_local 最后一次写桶 0 已是 600 MB 之前）。
// 而 attach_l1 的消费次序**恰好就是** sort_x_local 的生产次序（都按 x 升序＝按桶升序）。
// 【改动（逐位等价，可证：输出只依赖 (x,y) 的多重集与计数，不依赖读写时序）】
// ① gH1（L1 数字直方图）从局部排序的第 3 趟搬到**第 1 趟**（直方图只认多重集）——这是本融合
// 此前做不成的原因：消费者要的游标前缀必须在**第一个桶被消费之前**就完整；
// ② 两趟桶内排序的散写目标从 `F0+12*start`／`Y2` 搬到一块专属 4 MiB 暂存（一半给第 2 趟
// 输入、一半给**紧凑的**第 2 趟输出）—— 因为 `Cb=(u8*)Y2` 在两半里都用，且原布局
// `FR==RC==FMem` 的复用**只在"attach 严格晚于 sort"时才安全**，融合后必须解耦；
// ③ 于是 F0(=FW) 在 pass A 之后**整个空闲**，正好当消费者输出 `Dst`（4N ≤ 6N ✓），
// levels 的 X1 两路都改读 FW，staging 仍走 RB。
// 净效果：**每半省掉 600 MB 的 F0 写 + 600 MB 的回读**（全 DRAM → L2 命中）。
// 【闸门（判题机同款 gcc 9.3 + `work/b18z_drv.cpp`，双尺度，逐位）】n=1e7 ⇒
// FNV 6982038405141804923 / SUM 14324786532795827997；n=1e8 ⇒ FNV 14022627980954509510 /
// SUM 8456426522220901567 —— 与归档基线逐位相同 ✓；A/B（同进程交替、n=1e8、7 轮）⇒ 中位
// 比值 0.907（**−9.3%**，7 轮里 6 轮 <1）· 每轮 `diff=0`（两实现输出逐元素相同）✓
// 【形状闸】`tools/gcc93check.py`：gcc 9.3 编译通过、无 ≥3 级诊断；vinserti128=4 / vmovdqu-xmm=4
// —— 与基座**完全一致**（未引入新的 32 B 载入拆分）✓
// ===== 本次思路（n08h_ 席，2026-10-03）=====
// 【本发 = 在册最好件 #122143 正文 + 一刀 [n08h-V9]（本席原创）】
// §一 定位（判题机同款 gcc9 -O2 + 判题机 probe 消融，非本机猜测）：本引擎每元素每半搬运约 144 B，
// 其中**局部排序链**（FW 12 B → F1 12 B → temp 9 B → packed 8 B）独占约四成 ⇒ 属"换表示"族
// （本队今日多题实证：换表示按 DRAM 行 1:1 兑现）⇒ 唯一有量的杠杆。
// §二 刀 [n08h-V9]：局部排序 pass 1（散播进 2048 路 x 桶）的输出**不需要存 x 的高 11 位** —— 桶号
// (x>>16)&2047 已隐含它们。原实现照抄 split 的 12 B 记录（u64 = i:27|x:27，u32 = y），再由桶内
// 第一趟 pass 把它重打成 9 B；本发让 pass 1 **直接产出那个 9 B 记录**（u64 = [y 低 21 位:43..63 |
// i:16..42 | x 低 16 位:0..15]，第 9 字节 = y>>21），桶内第一趟随即退化为 **9 B 原样复制**
// （原 12 B 读 + 重打包 5 条 ALU 全删）。语义/LSD 各位提取/最终记录格式**逐位不变**。
// §三 闸门：判题机同款 gcc9 -O2，n ∈ {2e6, 3e6, 6e6, 1.2e7} × seed ∈ {999, 12345} 共 8 组，
// out[] 校验和与在册最好件**逐位相同**；判题机 probe（自定义测试 n=1.2e7 / 1.6e7，各 3 发）FNV 同值 ✓
// §四 定价（判题机 probe，同批 3 发取 min）：n=1.2e7 **575.32 → 570.83（−4.49 ms，−0.78%）**；
// n=1.6e7 **774.87 → 768.55（−6.27 ms，−0.81%）** ⇒ 优势随 n **不衰减**（按仪器律做了 n 趋势）。
// 本机 Skylake-X 墙钟噪声 ≥1.7% ⇒ 未用于定价（仅作否决）。
// §五 同批判负（未采纳）：`out[]` 重放二级分桶（16 KB L1 窗，−0.39%：只降 L1 未中、不降 DRAM 流量）·
// 非时序存储 movnti+sfence（+20.7%，半行 WC 冲刷反噬）· staging 记录 8→6 B（−0.13%，6 B 步长破对齐）。
// ================#pragma GCC optimize("O2,align-functions=64") // [e8y]
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #120937 <https://duck.ac/submission/120937>
// （3291.421149 ms，本题场上最快）—— **本文件正文的来源**。本发 = 其正文 + 见下"思路"节。
// 其正文本身又是「saffah_cc_v41_agg1 #111438（我方件）正文 + sort_x_local 一级 MSD 局部基数」
// 演化而来，整条上游引用链（#120790 / #110329 / #109849 / #109709 / duck.ac 题面）在正文头注原样保留 ✓
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #111438 <https://duck.ac/submission/111438>
// （3360.593443 ms，我方现役最好件）—— 正文里 `gPB`/`gLB` 单基址表重排（`[e8y-bd]`）出自我方此件 ✓
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」，即站上提交按站点规则公开；
// 本文件逐条标注原作者账号、原提交地址与用途，未改动其正文一个字节 ✓
// ======================
// ===== 思路 =====
// 【本发 = 对手当代件 #120937 正文**逐字节复制**（我方旧件 #111438 = 3360.593443 ms 已明显落后）】
// 目的：把「我方现役最好件」从 3360.593443 抬到对手当代水平（≈3291.421149），作为下一步
// 「另补一刀」的基座 —— 全队律「抄件上限 = 与对手件平手 ⇒ 必须另补一刀」✓
// 【逐位闸门（判题机同款 gcc 9.3 + `work/b18z_drv.cpp`，双尺度）】n=1e7 ⇒ FNV 6982038405141804923 /
// SUM 14324786532795827997；n=1e8 ⇒ FNV 14022627980954509510 / SUM 8456426522220901567
// —— 与本题归档基线**逐位相同** ✓（复制无一处改动 ⇒ 输出不可能变 ✓）
// ======================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/111438:
* retained the split-half rank engine and inherited citations.
* saffah_codex_6a_agg3, https://duck.ac/submission/120790:
* retained the accepted local-radix/AoS first-level implementation.
* Public source displayed no separate software license notice.
* Idea: Use 2048 high-x buckets with local 8/8-bit passes, changing
* the working-set balance between random partition streams and local scratch.
* Purpose: Official geometry experiment motivated by scaling-dependent sort costs.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/111438:
* retained the accepted split-half rank engine and inherited citations.
* saffah_codex_6s_agg2, https://duck.ac/submission/110329:
* adapted high-x MSD partition plus bucket-local radix sorting, and rolling
* previous-y vectors for exact equal-x correction.
* Public sources displayed no separate software license notices.
* Idea: Replace three full-array LSD passes by a high-x partition followed by
* two bucket-local 8/7-bit radix passes, reusing dormant Y2 memory as scratch.
* Keep sorted 12-byte records and let first-level counting consume them directly,
* eliminating a global layout-conversion pass. Overlarge buckets use the original
* sorter; every scratch bound is checked before partitioning begins.
* Purpose: Official experiment reducing full-array traffic without increasing peak memory.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #109849 <https://duck.ac/submission/109849>
// （3427.077666 ms，本题场上最快）—— **本文件正文的来源**。本发 = 其正文 + 一处
// **纯表布局重排**（`[e8y-bd]`，见下"思路"节；语义逐位不变，判题机接口/输出零变化）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #109709 <https://duck.ac/submission/109709>
// （3455.456222 ms，我方现役最好件）—— 他的 #109849 就是在我方件上改出来的：剥注释 diff 仅 2 处
// （① 基数位段 9/9/9 -> **10/8/8**（桶 512->1024/256/256）② `FMem = bigalloc(M*16 + 63)` + **64 B 对齐**）✓
// [3] 其自带的整条上游引用链在本文件正文头注中原样保留（duck.ac 题面 + 前面各代提交）✓
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」，即站上提交按站点规则公开；
// 本文件逐条标注原作者账号、原提交地址与用途，未改动其正文一个字节 ✓
// ======================
// ===== 思路 =====
// 【本发 = 对手当代件 #109849 正文 + 一处 **纯表布局重排**（`[e8y-bd]`）：把两套数字计数器
// 的三张静态表各并成"一张表 + 一个基址寄存器"。零语义、零输出变化；只动地址计算。】
// 【依据（判题机同款 gcc 9.3 -O2 的反汇编）】`attach_l1` 内层环 = **78 条指令/元素**，其中
// ~10 条是纯表地址计算（`leaq _ZL3gP2/_ZL3gP3(%rip)` 在**环内**被反复重载）+ ~8 条栈重载。
// 根因：`gP1/gP2/gP3/gMsk` 四组静态表各要一个基址寄存器（-O2 输出里 r11/r13/r15/rdx 全被占），
// 叠加 X / Yv / Xout / Dst / Cb / gCur 等 ⇒ 16 个通用寄存器装不下 ⇒ 溢出到栈。
// 【改动（逐位等价，可证）】把 `gP1[8]/gP2[8][8]/gP3[8][8][8]/gMsk[8][8]` 并成
// `gPB[8+64+512+64]`（gp1@0 / gp2@8 / gp3@72 / 掩码行@584），`gL3P1[16]/gL3P2[16][16]/gL3Inc[16][16]`
// 并成 `gLB[16+256+256]`。因为 **a*8+b == d>>3** 且 **a*64+b*8+c == d**（d<512 恒成立）——
// 三维下标本来就能压成一维下标：
// · 查表 `cnt_q(d) = gPB[d>>6] + gPB[8+((d>>3)&63)] + gPB[72+(d&511)]`（三处都是 base+位移+变址）
// · 插入的 8 车道掩码后缀加，其组基址 = `gPB[72+(d&0x1F8)]` / `gPB[8+((d&0x1C0)>>3)]` / `gPB[0]`
// 同一套做法用于 l3：`l3_q(d) = gLB[d>>4] + gLB[16+(d&255)]`，组基址 `gLB[16+(d&0xF0)]`。
// ★ 32 B 对齐约束原样保持（组基址步长恒为 32 B，`vpaddd/vpaddw (mem),ymm` 的 VEX 要求不破）✓
// 【闸门（判题机同款 gcc9 + `b18z_drv`，双尺度）】n=1e7 ⇒ FNV 6982038405141804923 /
// SUM 14324786532795827997；n=1e8 ⇒ FNV 14022627980954509510 / SUM 8456426522220901567
// —— 与归档基线**逐位相同** ✓
// 【指令数（仪器 = cachegrind，判题机几何 I1/D1=32K,8,64 · LL=6M,12,64）】n=3e6：
// **1,027,128,387 → 985,515,547 = −41,612,840（−4.05%）** ✓（本机墙钟不可用于本题定价 ✓）
// ================#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 << 32) | (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) {
/* [n08e-THP] this engine touches ~1.6 GB of anonymous memory inside the judged
call. With 4 KB pages that is ~400k minor faults, each paying a fault entry
plus a 4 KB clear; with 2 MB pages it is ~800 faults. MADV_HUGEPAGE is only
advice: on a kernel without THP (or with THP=never) it is a no-op, and the
return value is deliberately ignored. No effect on results. */
madvise(p, bytes, MADV_HUGEPAGE);
return p;
}
{
void *q = malloc(bytes); // mmap fails on JudgeDuck -> malloc
if (q) madvise(q, bytes, MADV_HUGEPAGE);
return q;
}
}
// ------------------------------------------------------------------
// 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;
u32 ia;memcpy(&ia,(const unsigned char*)it+1,4);
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-=dup;
if constexpr(!preseed)*xo++=xvv;
u16 high;memcpy(&high,(const unsigned char*)it+6,2);
const u32 d1=(u32)high>>4;
cnt=cnt+cnt_q(d1);
cnt_ins(d1);
const u32 qq=gCur[d1]++;
unsigned char *rank=(unsigned char*)Dst+9u*(size_t)qq;
*(u64*)rank=((u64)cnt<<44)|(((u64)ya<<27)&0x00000fffffffffffull)|(u64)ia;
rank[8]=(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+9*gCur[((u32)v>>16)&2047u]++; /* [n08h-V9] 9-byte record */
*(u64*)q=(v&65535u)|((v>>11)&(M27<<16))|((u64)y<<43);
q[8]=(u8)(y>>21);
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+9u*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+9u*j;u64 v=*(const u64*)p;
unsigned char*q=temp+9*p0[((u32)v)&255u]++;
*(u64*)q=v;q[8]=p[8];
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+9*j;u64 v=*(const u64*)p;
const u32 mid=((u32)v>>8)&255u;
packed[p1[mid]++]=((v>>8)&~255ull)|(v&255u)|((u64)p[8]<<56);
}
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;
if (!preseed) Xout[a] = xvv;
u32 d1 = (ya >> 17) & 511u;
cnt = cnt + cnt_q(d1);
cnt_ins(d1);
u32 qq = gCur[d1]++;
unsigned char *rank=(unsigned char*)Dst+9u*(size_t)qq;
*(u64*)rank=((u64)cnt<<44)|(((u64)ya<<27)&0x00000FFFFFFFFFFFull)|(u64)ia;
rank[8]=(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)));
}
struct L3ByteMasks { u8 data[16][16]; };
static constexpr L3ByteMasks make_l3_byte_masks(){
L3ByteMasks m{};
for(unsigned d=0;d<16;++d)for(unsigned j=0;j<16;++j)m.data[d][j]=(u8)(j>d);
return m;
}
alignas(16) static constexpr L3ByteMasks gL3ByteMasks=make_l3_byte_masks();
alignas(16) static u8 gL3Byte[16+256];
static inline void l3_byte_init(){memset(gL3Byte,0,sizeof(gL3Byte));}
static inline u32 l3_byte_q(u32 d){return (u32)gL3Byte[d>>4]+gL3Byte[16+d];}
static inline void l3_byte_ins(u32 d){
__m128i top=_mm_load_si128((const __m128i*)gL3Byte);
top=_mm_add_epi8(top,_mm_load_si128((const __m128i*)gL3ByteMasks.data[d>>4]));
_mm_store_si128((__m128i*)gL3Byte,top);
u8 *row=gL3Byte+16+(d&240u);
__m128i low=_mm_load_si128((const __m128i*)row);
low=_mm_add_epi8(low,_mm_load_si128((const __m128i*)gL3ByteMasks.data[d&15]));
_mm_store_si128((__m128i*)row,low);
}
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){
const unsigned char *rank=(const unsigned char*)X1+9u*s;
const u64 rec=*(const u64*)rank;
u32 count_bytes;memcpy(&count_bytes,rank+5,4);
PUT_STG(stg,(u32)(rec&M27),count_bytes>>4);
}
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++) {
const unsigned char *rank=(const unsigned char*)X1+9u*t;
u64 rec=*(const u64*)rank;
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 count_bytes;memcpy(&count_bytes,rank+5,4);
u32 cnt = (count_bytes >> 4) + inc;
// 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 <= 255) {
l3_byte_init();
for(size_t a=0;a<nn;++a){
u64 rec=src[a];
u32 d3=(u32)(rec>>27)&255u;
u32 inc=l3_byte_q(d3);
l3_byte_ins(d3);
u32 cnt=((u32)(rec>>37)+inc)&0x7ffffffu;
PUT_STG(stg,(u32)(rec&M27),cnt);
}
} 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 >> 32)] = (u32)v; /* [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;
madvise((void *)out, N * 4, MADV_HUGEPAGE); /* [n08e-THP] see bigalloc */
// 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);
}