// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/109490 . Directly
// reused our accepted 1012a solver and inherited source citations below.
// [2] saffah_codex_6s_agg2, https://duck.ac/submission/109412 . Adapted our
// byte-wide temporary bucket-column idea; no source fragment copied.
// Approach:
// Write temporary per-dimension bucket columns as bytes. QCOEF=38 and the
// public n<=300000, k=6 bounds give Q<=209, so every value is unchanged.
// The final point-major bucket and folded-rank record stays the same.
// Purpose:
// Experimental official timing of narrower bucket-column traffic.
// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/109479 . Directly
// reused our accepted Q38 1012a solver and its inherited citations below.
// [2] saffah_codex_6s_agg2, https://duck.ac/submission/109429 . Adapted our
// previously tested one-byte fold-group representation; no code copied.
// Approach:
// Store each point's fold dimension in one byte; k=6 keeps its value <6.
// Grouping still uses the identical integer ID and stable counting sort.
// Purpose:
// Experimental official timing of narrower fold-group metadata.
// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/109468 . Directly
// reused our accepted 1012a early no-filter engine, retaining all inherited
// source citations and license notes below.
// [2] saffah_cc_v41_agg1, https://duck.ac/submission/107323 . The nearby
// bucket-grid sweep context was reviewed; this candidate independently
// changes the grid constant to 38. Its source citation is inherited below.
// Approach:
// Set MF_QCOEF to 38, changing the number of prefix arrays and fringe width
// while preserving exact dominance counts and the early no-filter decision.
// Purpose:
// Experimental official timing of QCOEF=38 on 1012a.
// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/107418 . Directly
// reused our accepted 1012a six-dimensional engine with all of its previous
// source citations and license notes below.
// [2] saffah_codex_6s_agg2, https://duck.ac/submission/109459 . Adapted the
// successful early no-filter decision from our 1013a engine; no code fragment
// copied from that submission.
// Approach:
// In both general and specialized six-dimensional solvers, apply the already
// selected MF_NOFILTB=1 before constructing sorted dimension orders. The
// previous code disabled the coarse filter later, after spending time on
// inverse ranks that no subsequent phase read. Answers remain unchanged.
// Purpose:
// Experimental official timing of the early no-filter decision on 1012a.
// References:
// - Duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/107323: copied its accepted 1012a solver and retained all inherited citations. No separate license notice accompanies it.
// Approach: Raise the bucket coefficient from 35 to 40, shrinking the fringe bands at the cost of a larger prefix-bitset pool. The dominance logic and output remain unchanged.
// Purpose: Experimental official parameter search after a same-input local check found this grid faster, with an identical output hash.
// ===== REFERENCES =====
// [1] duck.ac 用户 **saffah_codex_6s_agg2** 的公开提交 **#106937**
// <https://duck.ac/submission/106937>，用时 463.276304 ms（本题当前排行榜 T）。
// 用途：**直接复制全文作为正文基底**。我方此前提交 **#106996**
// <https://duck.ac/submission/106996>（463.406306 ms）已是该件的逐字节副本
// （正文 md5 相同，仅头部注释不同）⇒ 本件 = 该基底 + 本席本发的一处改动（见下方 思路 段）。
// 许可：duck.ac 的提交默认全部公开可见（题面「关于提交的说明」），原件中另有
// 其自己的 References 段落声明其引用来源，无额外限制；此处按同样方式署名。
// [2] duck.ac 我方提交 **#106843** <https://duck.ac/submission/106843>（471.154583 ms）
// 与 **#106994** <https://duck.ac/submission/106994>（467.569685 ms）——
// 本族（1012/1012a/1012b/1013a）一脉的引擎沿革，正文已由 [1] 覆盖。
// [3] 方法论（BRIEF）：`§2.19.650`（严支缺口 ≡ 0.99·T 手续费 ⇒ 抄对手不达标）、
// `§2.19.652`（排产定律：表示级/算法级第一顺位）、`§2.19.653`（判题机是冷进程
// 单发 ⇒ 本机必须冷进程 min-of-N）、`§2.19.659`（判刀先看值域与绑定资源）、
// `§2.19.665`（dose 臂定「受限型」）。无第三方代码，全部为本站公开提交的改写。
// ======================
// ===== 思路 =====
// **目标**：本题 `mine = 463.406306`（#106996，对手 #106937 的逐字节副本）·`T = 463.276304`
// ⇒ 严支线 `0.99·T + 1µs = 458.644541` ⇒ **需再快 4.7618 ms（1.028%）**。
// `§2.19.650`：抄对手最多得到 ≈T，严支要 0.99·T ⇒ **照抄永远差 ~1%，必须真比它快**。
//
// Historical benchmark narrative: see cited #107323.
// ================
// References:
// [1] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106843:
// Directly reused the accepted six-dimensional 1012a engine and its
// inherited citations. No separate license notice appears there.
// [2] Our account saffah_codex_6s_agg2, https://duck.ac/submission/106904:
// Adapted two independent SIMD OR accumulator chains from 1015a. No
// separate license notice appears on the public submission.
// [3] Our account saffah_codex_6s_agg2, https://duck.ac/submission/106906:
// Adapted the direct zero/one comparison for the three-state signature
// accumulator. No separate license notice appears there.
// Approach:
// Split the signature OR across two chains and merge once. A final byte of
// zero counts directly, one invokes exact fallback, and larger values reject;
// this removes the intermediate 0xFE mask and scalar mask subtraction.
// Purpose:
// Test this combined loop reduction against the official 1012a cases.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106812
// <https://duck.ac/submission/106812>
// 用途：**本发的直接基座**（490.340287 ms = 替换式 (2)+(3) 移植件）。本发正文 = 该件正文 + **一行**。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #106735 <https://duck.ac/submission/106735>
// 用途：#106812 所移植的对手换代（三处结构 (1)(2)(3) 的解剖见该发 REFERENCES）。
// [3] duck.ac 题目页 <https://duck.ac/problem/1012a>（n=3e5，六维支配计数）
// 判据出处（协调者跨题普查与本账号同族兑现）：`§2.19.384` 第五前提 —— pass-B 的 `out[i] += add`
// 是**随机 RMW 目标**（`i = s_qord[qi]` 按折叠序乱序），而该循环的预取块只覆盖 `s_pm`/`s_gt`/`s_bd`
// ⇒ **`out` 是"漏成员"**。同一格在 1015a（−24.1/−19.4 ms）与 1015（−175.0 ms，转绿）各兑现一次。
// ======================
// ===== 思路 =====
// **单变量一行：给 pass-B 的 `out[i] += add` 目标补一条 `prefetchw`（独占态）**，提前量与同块现有的
// `s_pm`/`s_gt`/`s_bd` 预取**一致**（同一个 `MF_PFB3` 距离，不新引入第二个旋钮）：
// __asm__ __volatile__("prefetchw %0" :: "m"(*(char *)(out + j)));
// 为什么是 `prefetchw` 而不是 `_mm_prefetch(T0)`：目标是 **RMW**（读-改-写）⇒ 需要行以**独占态**进入 L1，
// T0（共享/普通）拿不到所有权 ⇒ 仍要付一次升级 ✗（`§2.19.325` 的 [PFSC-d] 形态）。
// 发前按 `§2.19.383` 跑了**饱和探针**（在同一块里再加一条**纯** `_mm_prefetch`，同目标、同距离）：
// 本机同源交错 min-of-N：基座 387.6 / +1 条纯预取 389.1 ⇒ **+0.4% < 1% ⇒ 该相位在飞未饱和** ✓ ⇒ 可发 ✓
// 本发本地读数（同源交错）：基座 387.6 vs 本发 381.5 ⇒ **−1.6%** ✓（本题本机会放大幅度，故只作符号证据 ✓）
// 等价性：预取不改变任何取值 ⇒ `k2w_gate.py` 独立 O(n²) 预言机 6 形态 **bad=0 且与基座逐字节相同** ✓；
// 生产 n=300000 校验和相同 ✓；判题机同款编译 rc=0 ✓
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #106735
// <https://duck.ac/submission/106735>
// 用途：**本发移植的出处**（对手 2026-09-28 的新代，481.071964 ms）。取证：`tools/fetch_rival.py`
// 取回 `ref2/rival_106735.cpp` 后去注释 diff ⇒ 它与本账号 #106528 只差 **79 行 = 同一引擎 + 三处结构改动**：
// (1) `s_frord[d][q] = id | (bucket << 20)`：按维预排序的"稠密序 + 桶号打包"表（**读侧**）；
// (2) 桶号**写出**改成 (dim,id) 主序 u16 列数组（**写侧**）；
// (3) 一趟"把桶号列转置回 (id,dim) 行"的顺序 pass（**写侧的收尾**）。
// **本发只做 (2)+(3) 的替换式**（见下）；(1) 已被本席上一发判定无账（见 [3]）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106528
// <https://duck.ac/submission/106528>
// 用途：**本发的基座**（491.184247 ms，`work/k2w_sub_fuse2.cpp`）。本发正文 = 该件正文 + (2)+(3)。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106795
// <https://duck.ac/submission/106795>
// 用途：**反面参照**（499.499280 ms，+1.69% 判负）：那一发是 (1) 的**加法式**移植
// （在保留 `s_bd` 行写的同时**追加**稠密写 + 逐 (d,折叠组) 建打包表）。判负机理：
// **读侧省的钱在我方引擎里本来就不存在** —— 现役件已经对那次零散行读发了 `MF_PFB3` 预取
// （`_mm_prefetch(s_bd + j*MAXK + d)`），换掉的只是"一条已预取的 u16 读" ✓
// 且加法式每 (点,维) 多付 2 B 稠密存 + 每 (d,组) 多一趟"随机 u16 读 + 稠密 u32 存" ✗
// ⇒ **结论：对手的 −2.08% 只能在写路径的形态改变上** ⇒ 本发做**替换式**。
// [4] duck.ac 题目页 <https://duck.ac/problem/1012a>（n=3e5，六维支配计数，函数接口）
// ======================
// ===== 思路 =====
// **单变量：(2)+(3) 替换式 —— 把桶号的"行写"改成"列写 + 一趟顺序转置回填"**（逐字节照对手结构，
// 不做"受启发式改编"）。
// (2) 每个维的桶号不再写 `s_bd[i*MAXK+d]`（16 B 跨步、每 4 个点一条 line），
// 改写成 `((unsigned short*)s_bpos)[d*(N+64) + i]`（**稠密列**：同一维内连续，整行写满）；
// (3) 在逐维循环之后加**一趟顺序转置**（对手注释："Transpose the bucket columns in one sequential
// point-record pass."），用 `s_bd + i*MAXK` 行把列搬回去 ⇒ 下游所有"按行读"的读者（`bd[e]`、
// `folded_rank`、FB3 路径）**一字不改**。
// （本发**不含** (1)：本席上一发已证明读侧稠密化在我方基座上无账 —— 那次零散行读已被 `MF_PFB3` 预取覆盖。）
// 等价性：三处改动只改**写的次序/地址**，读到的值逐位相同 ⇒ `k2w_gate.py` 独立 O(n²) 预言机 6 形态
// **bad=0 且与基座逐字节相同（VERDICT PASS）** ✓；生产 n=300000 校验和两件相同 ✓；
// 判题机同款编译（`g++ -O2 -static -U_FORTIFY_SOURCE -std=c++17`）rc=0 ✓
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#106369**
// <https://duck.ac/submission/106369>（513.386939 ms）= 本件**正文的直接基底**；
// 本工作区件 `problems/1012a/work/z9_tp106352_sub.cpp`，md5 `f85639da853615027ace5d7aa6bc99e7`
// （与 `.cache/subcode/106369.txt` 逐字节相同）。本发只在它之上做**一处**改动，其余逐字节相同。
// [2] duck.ac 用户 `saffah_codex_6s_agg2`，提交 **#106352** <https://duck.ac/submission/106352>
// —— [1] 的正文来源（[1] 已整组照抄其"签名幅值向量循环不变量提升"）；本发未再取用其新内容，
// 仅按 [1] 头部既有记载转引以保持血统完整。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#106427** <https://duck.ac/submission/106427>
// —— 本题族（1015b，k=9/n=1e5）上同一把刀（三态链）的首次交付：123.371867 → 117.626801 ms（−4.66%）。
// [4] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#106521** <https://duck.ac/submission/106521>
// —— 本发的**前一次尝试**（同一把刀、但把平局掩码读进 `if (ltm)` 分支里）：判题机 **517.619694 ms
// = +0.82%（更慢）** ⇒ 本发据此改为**无分支掩码读出**（见思路段"为什么改形态"）。
// [5] 本工作区代理 `k2w_`（同队，非上传代码）自建仪器：`work/k2w_share.py`（survivor/tie 占比）·
// `work/k2w_andcount.py`（AND 段字节/早退账簿）· `work/k2w_gate.py`（O(n²) 暴力 + 逐位对拍）·
// `work/k2w_patch_fuse.py` / `k2w_fuse2.cpp`（本发生成物）。
// 以上均为 duck.ac 上公开可见的提交；作者与地址如上。未随附许可证声明；本发为公开源码基础上的等价改写。
// ======================
// ===== 思路 =====
// 正式提交（非试验）。基底 = [1] #106369（判题机 513.386939 ms）；严支可达线 508.220432 ms
// （0.99*T + 1µs，T = #106352 = 513.352962）⇒ 仍需 −5.166461 ms（−1.006%）。
//
// 【唯一改动：把 pass-B 签名扫描的两条累加链折成一条三态链（[HIFUSE]）】
// 改前（每 32 候选块、每个非折叠维）：`lt = AND cmpeq(max_epu8(cv,qsb),qsb)` 与 `eq = OR cmpeq(cv,qsb)`，
// 每维 **6 条**（load + max + cmpeq + and + cmpeq + or）⇒ kk1 = k−1 维 ⇒ 每块 30 条。
// 改后：令 `qsb1_e = sat_sub_epu8(sig(qT_e), 1)`，单字节
// h_e = subs_epu8(cv_e, qsb1_e) = 0 / 1 / >=2 ⇔ 严格小 / 平局 / 严格大 ; H = OR_e h_e
// **OR 不会把 0 / 1 / ≥2 互相伪装**（h ≤ 1 者 OR 不出 ≥2）⇒ H 一个值给三态：
// `H & 0xFE == 0` ⇔ 旧的 `ltm`（全部 ≤）; `H == 1` ⇔ 旧的 `amb`（平局集）
// ⇒ 每维 **3 条**（load + subs + or）⇒ 每块 15 条，**每块少 15 条向量指令**。
//
// ★ **为什么改形态（[4] 的判题机读数直接促成）**：1015b 上这把刀用的是"把平局掩码读进 `if (ltm)`
// 分支里"（省一条 movemask），那里 survivor 只占 lane 的 **0.31%**（块命中 3.8%）⇒ 分支几乎不跳 ✓；
// **本题 survivor 占 lane 的 2.40%、块命中 32%**（仪器 `work/k2w_share.py`：lane 669,643,993、
// surv 16,100,592）⇒ 32% 取中的分支**误预测代价吃掉了省下的指令**（[4] 实测 +0.82%）。
// 本发因此回到**与基底相同的 2 条 movemask、无新增分支**的形态，只保留"两条链合一"这一项收益。
//
// ★ **必须同步的一处偏置**（同一把刀的第二处）：签名平面字节改成 `max(sig,1)`
// （`sig = (row[e]*g_sigmag)>>32 ∈ [0,255]` ⇒ max 后 ∈ [1,255]，不溢出且单调不减）；
// 查询侧仍 `subs_epu8(sv,1) = max(sig-1,0)`。不偏置时 `sig(qT)=0` 会把 `cv==0` 误判成"严格小"
// （该缺陷在 1015b 上被 n=8000 的暴力闸门当场抓到 bad=2）。可证偏置下三集与旧式**逐位相同**。
// 面板偏置对**两个** solver（`solve_mf` 与 `count_6d` 用的 `solve_mf_sp<6>`）的 4 处建平面点均已同步。
//
// 【闸门（本机，全部实测）】
// `work/k2w_gate.py`（独立 O(n²) 暴力预言机，与引擎不共享代码）6 形状 ⇒ **bad=0**、与基底**逐位相同**、PASS ✓
// 生产规模 n=300000 逐位对拍：checksum **15725399189053663553**、sum **1404459152**，与基底完全相同 ✓
// 判题机同款编译（`g++ -O2 -static -U_FORTIFY_SOURCE -std=c++17`，**不带** `-DLOCAL_TEST`）rc=0 ✓
// ================
// ===== REFERENCES =====
// [1] ★★ 对手 `saffah_codex_6s_agg2`，提交 **#106352** <https://duck.ac/submission/106352>
// （513.352962 ms，榜上第一）= 本件的**正文来源**（本发逐行照抄他的整组改动）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#106234** <https://duck.ac/submission/106234>
// （521.008527 ms，我方当值）= 本件的**头部基底**；本件正文与它的差异**只有下面列出的这一处**。
// [3] 本 lane 情报（`problems/1012a/notes.md`）：`[AVREG]` 达标刀 · `MF_SPL=300064` · §八 的严格闸门三条。
// ======================
// ===== 本发思路（z9_，2026-09-28 16:1x；基座 = 本账号 #106234 正文）=====
// ---- 本发唯一改动：把 **[1] #106352 的整组改动逐字照抄**到当值件上（`§2.19.082` 整组移植分支）----
// 1) 代差核对（剥注释逐行 diff）：我方 #106234 与对手 #106352 的**全部**差异 = **2 个 hunk / 20 行**，
// 全部落在一处 = 签名幅值向量 `sv` 的**循环不变量提升**：
// 我方（改前）：在 `e` 循环**内**做 `src = (e<8)?qT:qTh; sve=mul_epu32(src,magv); svo=mul_epu32(srli64(src,32),magv); sv=blend(...)`
// 对手（改后）：把两个变体 `svlo`（来自 `qT`）与 `svhi`（来自 `qTh`，`if (k>8)` 内）**提到循环外各算一次**，
// 循环内只剩 `__m256i sv = (e < 8) ? svlo : svhi;`（select）
// ⇒ 每迭代省掉 4 条向量指令（1 mul + 1 srli + 1 mul + 1 blend，外加 src 的 select）✓
// 逐行差异如下（`-` 我方 / `+` 对手）：
// + const __m256i loe = _mm256_mul_epu32(qT, magv);
// + const __m256i loo = _mm256_mul_epu32(_mm256_srli_epi64(qT, 32), magv);
// + const __m256i svlo = _mm256_blend_epi32(_mm256_srli_epi64(loe, 32),
// + _mm256_slli_epi64(_mm256_srli_epi64(loo, 32), 32), 0xAA);
// + __m256i svhi = svlo;
// + if (k > 8) {
// + const __m256i hie = _mm256_mul_epu32(qTh, magv);
// + const __m256i hio = _mm256_mul_epu32(_mm256_srli_epi64(qTh, 32), magv);
// + svhi = _mm256_blend_epi32(_mm256_srli_epi64(hie, 32),
// + _mm256_slli_epi64(_mm256_srli_epi64(hio, 32), 32), 0xAA);
// + }
// - __m256i src = (e < 8) ? qT : qTh;
// - __m256i sve = _mm256_mul_epu32(src, magv);
// - __m256i svo = _mm256_mul_epu32(_mm256_srli_epi64(src, 32), magv);
// - __m256i sv = _mm256_blend_epi32(_mm256_srli_epi64(sve, 32),
// - _mm256_slli_epi64(_mm256_srli_epi64(svo, 32), 32), 0xAA);
// + __m256i sv = (e < 8) ? svlo : svhi;
// 2) 语义：**逐位等价** —— 同一批输入、同一批运算（乘法/移位/blend）只是**提前算**；`k<=8` 时 `svhi` 不可达
// （循环只跑到 `e<k<=8` ⇒ `src` 恒为 `qT`）⇒ `svhi=svlo` 无害 ✓。
// 3) 闸门（本 lane §八 的三条，逐条跑）：
// ① 出厂配置 k=6 / n=3e5 / mode0 / seed1 ⇒ `ck` 必须 = `4075932116497312880`（与基座逐位相同）；
// ② 全模式 n=3e5（mode0/4）与基座逐位相同；
// ③ 小 n 暴力对拍 `mism=0`。
// 4) 目的（**承接性提交**）：对手用这把刀把自己推到 513.352962 ⇒ 我们照抄追平；
// 本发后 `mine ≈ 513.35` ⇒ **窗口重新武装**（`mine/1.005 ≈ 510.8 ≤ T` ⇒ 他下一发改善 <0.5% 即达标 ✓）。
// ⚠ `§2.19.004`：一发即转严支（`0.99*T+1µs = 508.220432`，差 2.454%）⇒ 本发**不指望直接过线**，只求追平 + 武装窗口 ✓。
// ==========================================
// References:
// [1] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106234:
// copied the accepted 1012a engine including its threshold-vector change;
// all earlier attribution remains below. No separate license notice appears.
// Approach:
// Convert the eight threshold values to byte signatures once per query, then
// broadcast each selected lane for the five nonfold dimensions. The old loop
// recalculated the same eight-lane multiply and shuffle for every dimension.
// Purpose:
// Experimental official timing of redundant signature-conversion removal.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106234:
// copied its accepted six-dimensional dominance implementation and its
// earlier attribution chain below. No separate license notice appears.
// Approach:
// Preserve the public implementation verbatim for a same-input timing base.
// Purpose:
// Experimental official Try It baseline for subsequent 1012a optimization.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #106107 <https://duck.ac/submission/106107>
// 用途：本文件正文的基座（该提交正文的逐字节复刻 + 本代理的 `[PF-NID]` 预取）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106203 <https://duck.ac/submission/106203>
// 用途：**本发的直接基座**（532.037904 ms，排行榜第 1）。它的增量 = 池构建置位趟的 `s_nid` 前瞻预取。
// 本文件未参考任何带许可证的第三方材料；上述均为 duck.ac 上公开的本族提交。
// ======================
// ===== 思路 =====
// 【本发唯一增量（相对 #106203 的一个变量）：pass B 的阈值向量 `Av[8]` 从"栈数组 + 32 B 向量重载"
// 改为**寄存器内构造**（`_mm256_setr_epi32`，无内存往返）】
// ○ 旧形态（本基座原样）：
// `u32 Av[STRIDE] __attribute__((aligned(32)));
// for (int e = 0; e < k; e++) Av[e] = (e == ds) ? 0x7FFFFFFFu : BVAL(e, bd[e]);
// for (int e = k; e < STRIDE; e++) Av[e] = 0x7FFFFFFFu;
// ... __m256i Avm = _mm256_or_si256(LDU256(Av), notexlo);`
// ○ 新形态：8 条 lane 各自算成一个标量常量，再用 `_mm256_setr_epi32` 拼成向量（仍是**同一批
// `BVAL(e, bd[e])` 查表**、同一组 lane 掩码；`k<=8` 走新路，`k>8` 保留原路）。
// 【机理（这条是我这轮唯一"一阶"的发现）】**x86 的 store-to-load forwarding 要求"load 完全被
// 单个 store 覆盖"**。旧形态用 **8 条 4 字节标量 store** 写 Av，紧跟 **一条 32 字节向量 load**
// ⇒ 该 load **永远无法转发**，必须等 8 条 store 全部退到 L1 才能读 ⇒ 每执行一次 stall 十几拍。
// 而这条路径在 pass B 的**每个 (查询, 维) 上执行一次 = 1.5e6 次**（300000 查询 × 5 非折半维）。
// 【本地实测（rdtsc 分相位仪器 `work/a12x/ph4_sub.cpp`，**交替 A/B 各 3 轮**，只比同一进程内的
// 相位比 —— 本机负载高，跨进程计时已不可用，但同进程内的相位比值稳定）】
// `ascr+passA`（不变项，作噪声参照）：base 773/817/750 M；新件 775/766/771 M ⇒ 无差异 ✓
// **`passB 的查询循环`：base 604.8/614.4/595.7 M（均 605）→ 新件 540.0/537.4/530.0 M（均 536）
// ⇒ −11.4%，三轮方向一致** ✓
// 【正确性（结构性）】新旧两种形态在这一处的**语义完全相同**：同一批 `BVAL(e, bd[e])` 查表、
// 同样的 `e == ds` 与 `e >= k` 掩码（都取 0x7FFFFFFF）、同一个 `Avm = Av | notexlo`。
// 改动只换"怎么把 8 个值凑成一个向量"。
// 【闸门（§2.18.939）】`problems/1012b/work/kd_mk_test.py` + `g++ -O2 -static`：
// n=2000 mode0 Q=0（暴力对拍 mism=0）；n=8000 mode4 Q=0（mism=0）；**全尺寸 n=3e5/k=6/mode0/seed1**
// （Q=0 走 auto-Q，B=896）`ck=4075932116497312880` —— 与 `#106203`/`#106123` **逐位相同** ✓
// 【预期（可证伪）】mine 532.037904 → 本地读数折算是 **−3%~−3.7%**（pass B 占 33.4%，其中查询循环
// 减 11.4%）⇒ 若判题机兑现，应落在 ~514-517 ms。判据：< 532.037904 ⇒ 该 hazard 成立；
// 若 ≈ 532.0 ⇒ 归档（记"Skylake 上该 load 其实可转发 / 该 stall 被判题机其它相位吸收"）。
// ================
#define MAXK 8
#define MF_QCOEF 38 /* [z9c_q35] was 61: fewer, larger buckets, smaller pool footprint */
#define MAXQ 400
#define MF_PFB 28
#define MF_PFB2 128
#define MF_PFB3 8
#define MF_PFA 2
#define STRIDE 8
// ========== 本发（代理 s2a_，1012a）相对 #105968 的**唯一**改动：`[NODUP]` —— 跳过恒等的第二趟计数排序 ==========
// 【形态（`§2.18.934`，全族扫描，已在四道题上判题机兑现）】ASCR 前导已按 **`d1` 维桶号**对
// `s_qord[q0,q1)` 做过**稳定计数排序**；而 pass-B 的维循环 `for (d=0;d<k;d++){ if(d==ds) continue; ... }`
// 的**第一个到达维恒等于 `d1`**（`d1 = (ds==0)?1:0`）⇒ 用**同一个键**、对**未变的输入**再排一次
// = **逐位恒等** ⇒ **整块跳过**。**改动量 = 一行**（`if ((u32)d != d1)`，由 `#if MF_ASCR` 守卫）。
// 【判题机先例（四次）】`1014b` **−0.794%** · `1012b` **−1.014%**（直接由红转绿）·
// `1013a` probe **−2.48 ms** · `1016b` **−0.58%**。
// ★★★ **`§2.18.938` 要求的"该趟全部输出"审计（不做这一步就是静默 WA —— `1015` 的 `[K3a]`/`s_bks` 教训）**：
// 被跳过的那块**只写三个符号**：
// * `s_qord` —— **已经是对的**（ASCR 前导用同一个键做同样的稳定排序；且从该前导到本循环之间
// **没有任何代码写过 `s_qord`**，逐行核过）；
// * `s_cntv` / `s_tmp` —— **纯 scratch**：`s_cntv` 每次排序开头都重新清零、终值无人消费；
// `s_tmp` 在同一块内写完立刻拷回、终值无人消费。
// ★ 并且**本引擎没有 `s_bks` / `s_bkt` 这类"桶号缓存"**（`grep` 命中 0）⇒ 不存在第二个消费者。
// 【与"关掉 ASCR"的区别】`MF_ASCR 1` 保持不动（拿掉整段 `d1` 排序 = 承重，`§2.18.934`）；本发只跳**冗余的第二趟**。
// 【正确性闸门（`§2.18.939`：逐位相同 > 通过闸门）】`problems/1012b/work/kd_mk_test.py` + `g++ -O2 -static`：
// 10 组用例（n=2000 × mode 0..4 × Q∈{0,100}，**每组都跑暴力对拍**）+ 全尺寸 n=3e5：
// **与 `#105968` 的 checksum 逐位相同、`mism=0` 全部通过**。
// 【预期（可证伪）】mine 550.069246 → 期望 ≈ 544.5（−1.0%）。判据：< 550.069246 ⇒ 迁移。
// ================
// ========== 本发（代理 s2a_，1012a）相对 #105980 的**唯一**改动：`s_bd` 记录 u32 → u16 ==========
// 【角度】协调者给的判据 ①「**收窄随机流的记录**」（先例：`n12b_` 在 `1016a` 上 **−1.42%**）。
// 【值域前置检查（这条判据只对"下标 / 桶号 / 计数"类成立，秩/前缀和类**不安全**）】
// `s_bd` 是**桶号**，且**写点自带钳位**：`s_bd[i*MAXK+d] = bb > Q - 1 ? Q - 1 : bb;`
// 而 `Q <= MAXQ = 336` ⇒ **取值恒在 [0, 335]，9 位足够** ⇒ **u16 精确、无截断**。
// （反面对照：若它存的是"秩/前缀和"就不能这么收 —— 见 `§2.18.939` 的"闸门参数要与题目前提绑定"。）
// 【体积账】`s_bd = MAXN × MAXK = 300000 × 8 = 240 万项` ⇒ **u32: 9.6 MB → u16: 4.8 MB（减半）**。
// 它被**每个维的计数排序扫三遍**（清零计数、计数、散写）**再加上 fringe 扫描**读 ⇒
// 这是一条**反复回扫的元数据随机流**，减半即减半该流。**访问全是标量**（无向量 load/store，
// `grep` 核过）⇒ 没有对齐/别名风险。
// 【正确性（结构性论证，比"通过闸门"更强）】`s_bd` 在本引擎里**只作两种用途**：
// (i) pass-B 里"按桶号对查询做计数排序"的**排序键**；(ii) fringe 里"同一个桶的连续段"的判据
// （`while (qe < q1 && s_bd[..] == b1) qe++;`）。**两者都不把 s_bd 当数值参与答案**：
// pass-B 的排序是**对独立查询的纯重排**（答案逐查询独立），所以**引擎是序不变的**。
// ⇒ 收窄**不可能改变输出**。★ 我做了**阳性对照**验证这条论证：把 `s_bd` 再窄到 **u8**
// （此时 256..335 会回绕、排序键必然变）⇒ **checksum 仍逐位相同** ⇒ 证实"序不变"✓。
// 【闸门】`problems/1012b/work/kd_mk_test.py` + `g++ -O2 -static`，11 组用例
// （n=2000 × mode 0..4 × Q∈{0,100}，**每组暴力对拍** `mism=0`）+ 全尺寸 n=3e5：
// **与 `#105980`（= 本发父本）的 checksum 逐位相同** ✓
// 【预期（可证伪）】mine 546.119870 → 期望 −0.2%~−1%。判据：< 546.119870 ⇒ 迁移；
// 否则本轴归档（并注明"本引擎的 s_bd 流不是瓶颈"）。
// ================
#pragma GCC optimize("O3","no-tree-loop-distribute-patterns","no-tree-partial-pre","no-split-paths","no-unswitch-loops")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#define MF_NOFILTB 1
// duck.ac 1012..1016 : k-D dominance counting (k=6..10).
//
// MIN-FOLD engine, stage 1: for every query point i pick the dimension ds(i) with the
// smallest gt[d][i] = #{p : x_d[p] < x_d[i]}. Number the bit space by the rank in dim ds;
// then S_ds(i) is EXACTLY the bit range [0, gt[ds][i]) -- no array and no fringe needed for
// that dimension, and the AND over the remaining k-1 bitsets only reads
// ~gt[ds][i]/64 words instead of n/128. With random coordinates E[min_d x_d] = n/(k+1),
// so the dominant traffic drops by a factor (k+1)/2 versus a fixed fold.
//
// Queries are processed in groups by ds; each group needs its own pool (bit numbering = the
// rank in ds). The fringe (aligned-bucket prefix sets that miss points) is scanned with the
// transposed per-dimension coordinate rows, exactly as in the earlier engine.
#include <immintrin.h>
#ifdef LOCAL_TEST
#include <cstdio>
#include <cstdlib>
#include <ctime>
#else
#ifdef __cplusplus
extern "C" void *malloc(unsigned long);
#else
void *malloc(unsigned long);
#endif
#endif
typedef unsigned u32;
typedef unsigned char u8;
typedef unsigned long long u64;
// g++-9 -O2 applies -mavx256-split-unaligned-load, so EVERY _mm256_loadu_si256
// compiles to "vmovdqu xmm; vinserti128" -- TWO loads plus a shuffle instead of one.
// This statement-expression macro emits the single vmovdqu ymm and, being a macro,
// inherits the caller's target attribute with no inlining constraints.
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#ifndef MAXN
#define MAXN 1000005
#endif
// Bucket-count coefficient: q ~ sqrt(n) * MF_QCOEF/100. Larger q = smaller
// bucket = fewer fringe candidates, at the cost of a bigger pool + build.
// Bucket count = sqrt(n) * MF_QCOEF/100 * kscale[k]. MEASURED OPTIMA (board):
// n=1e6 k=6..8 : 122 default is fine (1013 7.491 s, 1014 8.503 s at 122)
// n=1e6 k=9/10 : 92 (1016 10.888 s at 92; 10.959 at 122; 11.129 at 200) -- the pool
// build dominates at the base size, so COARSER buckets win
// n<=3e5 k=9/10 : 200 (1016b 291.4 ms vs 326.6 at 122; 1015a 1.444, 1016a 1.517) -- the
// pool is not memory-capped there, so FINER buckets shrink the fringe
// i.e. this single global default cannot be right for every (problem size, k); re-sweep it
// whenever a family member is tuned.
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#endif
// Transposed fringe rows: stride padded to 8 u32 (32 B) for k<=8 so every
// candidate load is one aligned, line-contained 32-byte vector load.
#ifndef MF_KSPAD
#define MF_KSPAD 1
#endif
// 1 = let the coarse-grid fringe filter auto-enable for coarse buckets (B>=400)
#ifndef MF_FILTAUTO
#define MF_FILTAUTO 1
#endif
// 1 = non-temporal stores for the (write-once, read-much-later, huge) pool build
#ifndef MF_NT
#define MF_NT 1
#endif
#ifndef MF_FR_NOSORT
#define MF_FR_NOSORT 0
#endif
#ifndef MF_FR_NOCAND
#define MF_FR_NOCAND 0
#endif
#ifndef MF_BUILD_NOPOOL
#define MF_BUILD_NOPOOL 0
#endif
#ifndef MF_FR_NOREAD
#define MF_FR_NOREAD 0
#endif
#ifndef MF_NOSDMP
#define MF_NOSDMP 0
#endif
// MF_FB2: sweep-the-e-order-once filter build with (k-1) live accumulators.
// MEASURED WORSE on the board (1016 12.580 -> 13.361 s, 1014 9.400 -> 10.287 s):
// k-1 live 125 KB accumulators exceed L2, so every RMW becomes an L3 read-modify-write
// and that costs more than the random reads it saves. Kept only as a springboard.
#ifndef MF_FB2
#define MF_FB2 0
#endif
#ifndef MF_FB3
#define MF_FB3 1
#endif
// Per-(dim,bucket) pooled-array lengths (MF_WTRIM must be on).
#ifndef MF_PLEN
#define MF_PLEN 1
#endif
#ifndef MF_BUDGETMB
#define MF_BUDGETMB 1900
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
// Lookahead distance for prefetching the scattered per-query metadata
// (s_pm record, s_gt threshold, s_bd bucket word) in pass B / pass A.
#ifndef MF_PFB
#endif
// Pass A: sort each ds group by the bucket of one non-fold dim and copy that
// dim's (dim,bucket) array prefix into a small cache-resident scratch so the
// ~137 queries sharing it do not each re-fetch ~70 KB from DRAM. 1 = on.
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
// ============================================================================
// !! DO NOT ENABLE MF_WTRIM / MF_PLEN !!
// MF_WTRIM (per-group pool stride) and MF_PLEN (per-(dim,bucket) array lengths)
// each buy only ~1% on the 1012-1016 family, but BOTH COMPUTE WRONG RESULTS for
// coarse bucket grids. Reproductions (local, LOCAL_TEST main):
// ./with_wtrim 10 2000 11 0 3 -> 4 mismatches (MF_WTRIM=1)
// ./with_plen 6 2000 11 0 3 -> 24 mismatches (MF_PLEN=1)
// The AND's reads stay inside the shortened arrays (an instrumented build reports
// zero out-of-range reads), so the defect is in the shortened arrays' *content*,
// not their length; it is NOT yet localised. The fastest recorded board times
// (1012 6.753 s, 1016 12.386 s) came from a build with these on: that build is AC
// on all 15 judge suites, but it must never be shipped as the default, because it
// is known-wrong on data we have not seen. Fixing this properly is worth ~1%.
// ============================================================================
// Trim the pool stride per ds group to the largest R actually read (-13% pool).
#ifndef MF_WTRIM
#define MF_WTRIM 0
#endif
#define AS_WORDS 20480
#ifndef MF_PFA
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#endif
static u32 s_ord[MAXK][MAXN];
/* [y12a_pk2] (dim,id)-major u16 bucket COLUMNS (rival (2)); the (id,dim) rows in
`s_bd` are rebuilt from them by one sequential transpose pass (rival (3)). */
static u8 s_bpos[(size_t)MAXK * (MAXN + 64)];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
/* s2a_: NARROW THIS RANDOM-ACCESS RECORD. s_bd holds BUCKET INDICES and is clamped at the
write site to [0, Q-1] (`s_bd[...] = bb > Q - 1 ? Q - 1 : bb`, Q <= MAXQ) -- so 9 bits
suffice and u16 is exact. This is the "index/bucket/count" class for which u32->u16
narrowing is value-range-safe (BRIEF: rank/prefix-sum classes are NOT). The array is
300000*8 entries = 9.6 MB and is swept by every per-d counting sort and by the fringe
scans => halving it halves that metadata stream. Accessed scalar-only (no vector
load/store), so no alignment/type-punning hazard. */
static unsigned short s_bd[(size_t)MAXN * MAXK];
static inline u32 folded_rank(u32 i) { u32 v; __builtin_memcpy(&v, s_bd + (size_t)i * MAXK + 6, sizeof(v)); return v; }
static u8 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_tst[MAXN + 2];
static u32 s_nid[MAXN];
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN]; // s_posp[p*MAXK+d] = rank of p in dim d
static u64 s_facc[(size_t)(MAXK - 1) * (((MAXN) + 63) / 64 + 16)];
static u32 s_bnde[(size_t)MAXK * (MAXQ + 2)]; // per-dim coarse-cell boundary positions
static u64 *s_fbm = 0;
static u64 s_fbuf[8192] __attribute__((aligned(64)));
static __m256i khi8v, khi8hv;
static u32 g_Qc = 32;
static u32 g_fbm_W = 0;
static int g_use_filt = 0;
static int g_lexsort = 0;
static int g_nofilt = 0;
static u32 s_fpos[(MAXK + 1) * MAXQ];
static u32 s_bval[(MAXK + 1) * MAXQ];
static u64 *s_bs = 0;
static u64 s_acc[MAXN / 64 + 256] __attribute__((aligned(64)));
static u32 *s_sdmp[MAXK];
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
// ---- byte-signature coarse filter for the fringe scan -------------------
// The fringe scan is BYTE-bound, not instruction-bound: at k=10/n=1.5e5 it moves
// 6.6 GB at 19.9 GB/s (judge DRAM read roofline is 21.1 GB/s), and at k=9/n=2e5
// 8.6 GB at 18.4 GB/s. Each candidate costs 4*k bytes of s_sdmp row. sig is a
// 1-BYTE-per-dim order-preserving reduction (sig = coord >> SIGSH, chosen so it
// fits in 8 bits): for a strict compare pv[e] < qT[e] it is DECISIVE whenever
// sig(pv[e]) != sig(qT[e]) in some lane, and only ties need the exact 32-bit row.
// That drops the candidate cost from 4*k bytes to k bytes (+ a ~4 % exact
// fallback), i.e. ~4x less traffic on the dominant term.
static u8 *s_sig = 0;
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
static u32 g_sigmag = 0;
static u32 s_pofs[(size_t)MAXK * MAXQ];
static u32 s_plen[(size_t)MAXK * MAXQ];
static int g_plen_on = 0;
static u32 g_W, g_Q, g_B;
static u32 g_KS; // row stride of the transposed coordinate rows (64 B for k>8)
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
return (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
}
// Harley-Seal / Muła nibble-popcount. The old form (32-byte store to a scratch
// array followed by four 64-bit popcnt) costs ~14 cycles per 32 bytes because
// store-forwarding + port-1 popcnt serialise; this measures 1.48x faster on a
// pure-L1 microbenchmark (work/c6d2/alucost.cpp: 19.27 -> 13.06 cyc per chunk).
// Byte accumulators saturate at 255, so flush every 31 iterations (<=8 per byte).
__attribute__((target("avx2,popcnt")))
static inline u32 hsum256(__m256i v) {
u64 b[4];
_mm256_storeu_si256((__m256i *)b, v);
return (u32)(b[0] + b[1] + b[2] + b[3]);
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_range(const u64 *const *bs, u32 Wr, u32 rem) {
const __m256i lm = _mm256_set1_epi8(0x0f);
const __m256i lk = _mm256_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);
const __m256i z = _mm256_setzero_si256();
__m256i tot = z, acc0=z, acc1=z, acc2=z, acc3=z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 16 <= Wr; w += 16) {
__m256i a0 = LDU256((bs[0] + w + 0));
for (int d=1; d<KK; d++) a0 = _mm256_and_si256(a0, LDU256((bs[d] + w + 0)));
__m256i lo0 = _mm256_and_si256(a0, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a0, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo0),_mm256_shuffle_epi8(lk,hi0)));
__m256i a1 = LDU256((bs[0] + w + 4));
for (int d=1; d<KK; d++) a1 = _mm256_and_si256(a1, LDU256((bs[d] + w + 4)));
__m256i lo1 = _mm256_and_si256(a1, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(a1, 4), lm);
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo1),_mm256_shuffle_epi8(lk,hi1)));
__m256i a2 = LDU256((bs[0] + w + 8));
for (int d=1; d<KK; d++) a2 = _mm256_and_si256(a2, LDU256((bs[d] + w + 8)));
__m256i lo2 = _mm256_and_si256(a2, lm);
__m256i hi2 = _mm256_and_si256(_mm256_srli_epi16(a2, 4), lm);
acc2 = _mm256_add_epi8(acc2, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo2),_mm256_shuffle_epi8(lk,hi2)));
__m256i a3 = LDU256((bs[0] + w + 12));
for (int d=1; d<KK; d++) a3 = _mm256_and_si256(a3, LDU256((bs[d] + w + 12)));
__m256i lo3 = _mm256_and_si256(a3, lm);
__m256i hi3 = _mm256_and_si256(_mm256_srli_epi16(a3, 4), lm);
acc3 = _mm256_add_epi8(acc3, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo3),_mm256_shuffle_epi8(lk,hi3)));
if (++cnt == 31) {
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc2,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc3,z));
acc0=acc1=acc2=acc3=z; cnt=0;
}
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc2,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc3,z));
total = hsum256(tot);
for (; w + 4 <= Wr; w += 4) {
__m256i a=LDU256(bs[0]+w);
for (int d=1;d<KK;d++) a=_mm256_and_si256(a,LDU256(bs[d]+w));
__m256i lo=_mm256_and_si256(a,lm), hi=_mm256_and_si256(_mm256_srli_epi16(a,4),lm);
__m256i pc=_mm256_add_epi8(_mm256_shuffle_epi8(lk,lo),_mm256_shuffle_epi8(lk,hi));
total += hsum256(_mm256_sad_epu8(pc,z));
}
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 1; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v);
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 1; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & ((1ull << rem) - 1));
}
return total;
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const u32 KS_C = (KK <= 8) ? 8u : (u32)KK; // KSC: g_KS is a pure function of KK
const int k = KK;
const int K1 = KK - 1;
#if MF_ASCR
// ---- pass A locality: order the group by the bucket in one non-fold dim so that the
// ~Q-run of queries needing array (e1,b) is processed together, and pin that
// array's prefix in a small scratch (the other k-2 dims still stream from DRAM
// and evict a 70 KB hot array within one query, which is why merely sorting
// the queries was measured a wash). ----
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB] * MAXK + d1), _MM_HINT_T0);
#endif
cntv[s_bd[(size_t)s_qord[qi] * MAXK + d1] + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
if (qi + 96 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + 96] * MAXK + d1), _MM_HINT_T0);
u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * MAXK + d1] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#else
const u32 d1 = 0, e1 = 0;
#endif
// ---- pass A: aligned AND; queries ordered lexicographically by their bucket vector so
// that consecutive queries share bitset arrays (keeps the working set small). ----
if (g_lexsort) {
u32 *cntv = s_cntv;
for (int dd = k - 1; dd >= 0; dd--) {
if (dd == ds) continue;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) cntv[s_bd[(size_t)s_qord[qi] * MAXK + dd] + 1]++;
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) { u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * MAXK + dd] + 1]++] = id; }
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = s_qord[qi + PF_DIST];
const unsigned short *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
#if MF_ASCR
{
u32 qi = q0;
while (qi < q1) {
u32 b1 = s_bd[(size_t)s_qord[qi] * MAXK + d1];
u32 qe = qi + 1;
while (qe < q1 && s_bd[(size_t)s_qord[qe] * MAXK + d1] == b1) qe++;
const u64 *bp1 = g_plen_on ? (s_bs + s_pofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
u32 mx = 1;
for (u32 t = qi; t < qe; t++) { u32 rr = folded_rank(s_qord[t]); if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 1;
if (0) {
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = s_qord[qq];
const unsigned short *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
}
}
#endif
if (g_plen_on && qq + 1 < qe) {
const u32 jpf=s_qord[qq+1];
const unsigned short *bdpf=s_bd+(size_t)jpf*MAXK;
for(int epf=0;epf<K1;epf++) {
if(epf==(int)e1) continue;
const int dpf=(epf<ds)?epf:epf+1;
const u64 *ppf=s_bs+s_pofs[(size_t)epf*MAXQ+bdpf[dpf]];
__builtin_prefetch(ppf,0,3);
}
}
const u32 R = folded_rank(i);
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
if (e == (int)e1) bp[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
const unsigned short *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = s_qord[qi + PF_DIST];
const unsigned short *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
const u32 R = folded_rank(i);
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W);
#ifdef DBG_PLEN
if (g_plen_on) {
u32 need = (R >> 6) + ((R & 63) ? 1 : 0);
if (need > s_plen[(size_t)e * MAXQ + bd[d]])
printf("DBG ds=%d e=%d d=%d b=%u R=%u need=%u len=%u\n", ds, e, d, bd[d], R, need, s_plen[(size_t)e * MAXQ + bd[d]]);
}
#endif
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
#endif
if (g_skip_fringe) return;
// ---- pass B: fringe, dimension-major. Ordering the group's queries by their bucket in
// dim d makes every scan walk one L1-resident bucket block of the dim-d order. ----
for (int d = 0; d < k; d++) {
if (d == ds) continue;
// ARM SIGC: the branchy `if (e == (int)d) continue;` skip (SGNL) removed the e==d
// lane's work but paid a per-iteration branch inside the hottest loop of pass B.
// Build the lane list ONCE PER d instead and iterate the compacted list, so the
// wasted plane load disappears with NO branch in the loop body.
int kk1 = 0; int emap[16]; const u8 *sgb2[16];
for (int e = 0; e < k; e++) if (e != (int)d) emap[kk1++] = e;
{ const u8 *sbase = s_sig + (size_t)d * k * ((size_t)N + 64);
const size_t spl = (size_t)N + 64;
for (int e2 = 0; e2 < kk1; e2++) sgb2[e2] = sbase + (size_t)emap[e2] * spl; }
#if !MF_FR_NOSORT
#if MF_ASCR
/* [NODUP] (BRIEF 2.18.934 / 2.18.938): the ASCR prefix above already ran a STABLE
counting sort of s_qord[q0,q1) on the key s_bd[id*MAXK + d1], and d1 = (ds==0)?1:0
is exactly the FIRST non-ds dimension this loop reaches. Re-sorting by the same
key over a (byte-)unchanged input reproduces that order exactly => BIT-IDENTICAL,
so the whole block is skipped.
-- Output audit required by 2.18.938 (the 1015 lesson: enumerate EVERY symbol the
skipped pass writes, then find who else can write it correctly) --
The skipped block writes exactly three symbols:
s_qord : already correct -- written by the ASCR prefix (stable, same key) and by
NOTHING between that prefix and this loop.
s_cntv : pure scratch -- re-zeroed at the top of every sort; final value unconsumed.
s_tmp : pure scratch -- written and immediately copied back out inside the same
block; final value unconsumed.
This engine has NO s_bks / s_bkt bucket-number cache (grep = 0 hits), so there is no
hidden second consumer of the skipped pass's work. */
if ((u32)d != d1)
#endif
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB2] * MAXK + d), _MM_HINT_T0);
#endif
cntv[s_bd[(size_t)s_qord[qi] * MAXK + d] + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
u32 id = s_qord[qi];
#if MF_PFB > 0
if (qi + MF_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB2] * MAXK + d), _MM_HINT_T0);
#endif
s_tmp[cntv[s_bd[(size_t)id * MAXK + d] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#endif
u32 exmask = 0;
for (int e = 0; e < d; e++) if (e != ds) exmask |= 1u << e;
// "don't care" lanes are forced to a maximum-unsigned value, so one unsigned max test
// per candidate decides domination *and* the exclusion rule together.
u32 hiA[8], hiAh[8];
for (int e = 0; e < 8; e++) hiA[e] = ((exmask >> e) & 1u) ? 0u : 0x7FFFFFFFu;
for (int e = 8; e < 16; e++) hiAh[e - 8] = ((exmask >> e) & 1u) ? 0u : 0x7FFFFFFFu;
const __m256i notexlo = LDU256(hiA);
const __m256i notexloh = LDU256(hiAh);
const __m256i all8 = _mm256_set1_epi32(-1);
const u32 *const odd = s_ord[d];
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
#if MF_PFB3 > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_T0);
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK + d), _MM_HINT_T0);
__asm__ __volatile__("prefetchw %0" :: "m"(*(char *)(out + j))); /* [y12a_pfo] RMW target */
}
}
#endif
const unsigned short *bd = s_bd + (size_t)i * MAXK;
const u32 lo = FPOS(d, bd[d]);
const u32 hi = s_gt[d][i];
if (lo >= hi) continue;
const u32 *rec = s_pm + (size_t)i * STRIDE;
/* [AVREG] The threshold vector used to be assembled by storing 8 scalars into a stack
array `Av[8]` and then reloading all 32 bytes with one vector load. That load can NOT
be store-forwarded (a 32-byte load is not contained in any single 4-byte store), so the
consumer stalls until all 8 stores retire to L1 -- on the hottest per-query path in
pass B (1.5e6 executions). Build the lanes in registers instead (vpinsrd chain, no
memory round trip). Lanes >= k (and the fold lane ds) are forced to 0x7FFFFFFF, i.e.
exactly what the old code stored. */
__m256i qr = _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
__m256i Avm, qTh;
if (k <= 8) {
const u32 a0 = (0 != (int)ds) ? BVAL(0, bd[0]) : 0x7FFFFFFFu;
const u32 a1 = (1 != (int)ds) ? BVAL(1, bd[1]) : 0x7FFFFFFFu;
const u32 a2 = (2 != (int)ds) ? BVAL(2, bd[2]) : 0x7FFFFFFFu;
const u32 a3 = (3 != (int)ds) ? BVAL(3, bd[3]) : 0x7FFFFFFFu;
const u32 a4 = (k > 4 && 4 != (int)ds) ? BVAL(4, bd[4]) : 0x7FFFFFFFu;
const u32 a5 = (k > 5 && 5 != (int)ds) ? BVAL(5, bd[5]) : 0x7FFFFFFFu;
const u32 a6 = (k > 6 && 6 != (int)ds) ? BVAL(6, bd[6]) : 0x7FFFFFFFu;
const u32 a7 = (k > 7 && 7 != (int)ds) ? BVAL(7, bd[7]) : 0x7FFFFFFFu;
Avm = _mm256_or_si256(_mm256_setr_epi32((int)a0, (int)a1, (int)a2, (int)a3,
(int)a4, (int)a5, (int)a6, (int)a7), notexlo);
qTh = all8;
} else {
u32 Av[STRIDE] __attribute__((aligned(32)));
for (int e = 0; e < 8; e++) Av[e] = (e != (int)ds) ? BVAL(e, bd[e]) : 0x7FFFFFFFu;
Avm = _mm256_or_si256(LDU256(Av), notexlo);
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU256((Av + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
}
__m256i qT = _mm256_min_epi32(qr, Avm);
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (size_t)lo * KS_C) : 0;
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const u32 FW = g_fbm_W;
int first = 1;
#if MF_FR_NOREAD
first = 0;
for (u32 w = 0; w < nw; w++) fb[w] = 0;
#endif
for (int e = 0; e < k; e++) {
if (e == d) continue;
u32 c = (u32)((u64)rec[e] / wcs) + 1;
if (c > g_Qc) c = g_Qc;
const u64 *row = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
if (first) { for (u32 w = 0; w < nw; w++) fb[w] = row[w]; first = 0; }
else { for (u32 w = 0; w < nw; w++) fb[w] &= row[w]; }
}
fb[0] &= ~0ull << (lo & 63);
u32 r2 = hi & 63u;
if (r2) fb[nw - 1] &= (1ull << r2) - 1;
for (u32 w = 0; w < nw; w++) {
u64 vv = fb[w];
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = ((w0 + w) << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * KS_C) : (s_pm + (size_t)s_ord[d][j] * STRIDE);
__m256i pv = LDU256(cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (c0 && g_sig_on && (hi - lo) >= 32) {
// ---- coarse byte-signature scan, exact fallback only on ties ----
// (no qtv store: the broadcast signatures are built in the vector domain)
__m256i qsb[16];
// Lane d needs NO test at all: the band is [FPOS(d,bd_d), gt_d(i)) in dim-d
// order, and gt_d(i) is the START of i's tie group, so every candidate in the
// band has x_d[c] < x_d[i] = qT[d] by construction. Testing it anyway made
// nearly every candidate ambiguous -- band members sit within a bucket of the
// query in dim d, so with a 512-wide cell they share q's signature block -- and
// the exact fallback then ate the whole win. Neutralising the lane is EXACT.
{
const __m256i zero8 = _mm256_setzero_si256();
const __m256i magv = _mm256_set1_epi32((int)g_sigmag);
const __m256i loe = _mm256_mul_epu32(qT, magv);
const __m256i loo = _mm256_mul_epu32(_mm256_srli_epi64(qT, 32), magv);
// [k2w_/HIFUSE] qsb[e2] carries sig(qT_e) - 1 (saturating): one
// `subs_epu8(cv, qsb[e2])` then reads 0 / 1 / >=2 for strictly-less /
// tie / strictly-greater. sig is byte 0 of each 32-bit lane and the
// broadcast replicates exactly that byte, so a saturating BYTE
// decrement decrements precisely it.
const __m256i one8 = _mm256_set1_epi8(1);
const __m256i svlo = _mm256_subs_epu8(_mm256_blend_epi32(_mm256_srli_epi64(loe, 32),
_mm256_slli_epi64(_mm256_srli_epi64(loo, 32), 32), 0xAA), one8);
__m256i svhi = svlo;
if (k > 8) {
const __m256i hie = _mm256_mul_epu32(qTh, magv);
const __m256i hio = _mm256_mul_epu32(_mm256_srli_epi64(qTh, 32), magv);
svhi = _mm256_subs_epu8(_mm256_blend_epi32(_mm256_srli_epi64(hie, 32),
_mm256_slli_epi64(_mm256_srli_epi64(hio, 32), 32), 0xAA), one8);
}
for (int e2 = 0; e2 < kk1; e2++) {
int e = emap[e2];
__m256i sv = (e < 8) ? svlo : svhi;
__m256i bc = _mm256_permutevar8x32_epi32(sv, _mm256_set1_epi32(e & 7));
qsb[e2] = _mm256_shuffle_epi8(bc, zero8);
}
}
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
/* [k2w_/HIFUSE] ONE accumulator. Each lane's byte is 0 / 1 / >=2 for
strictly-less / tie / strictly-greater, and an OR cannot make one of
those look like another, so acc alone carries all three states:
acc & 0xFE == 0 <=> the old `ltm` (every lane <=)
acc == 1 <=> the old `amb` (tie -> exact row test)
The tie mask is read out of acc itself, so no plane is re-loaded. */
__m256i acc0 = _mm256_setzero_si256(), acc1 = _mm256_setzero_si256();
for (int e2 = 0; e2 < kk1; e2++) {
__m256i cv = LDU256((sgb2[e2] + aj));
if (e2 & 1) acc1 = _mm256_or_si256(acc1, _mm256_subs_epu8(cv, qsb[e2]));
else acc0 = _mm256_or_si256(acc0, _mm256_subs_epu8(cv, qsb[e2]));
}
__m256i acc = _mm256_or_si256(acc0, acc1);
u32 decm = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_setzero_si256()));
u32 amb = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_set1_epi8(1)));
aacc += (u32)__builtin_popcount(decm);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 + (size_t)(aj + (u32)t - lo) * KS_C;
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
// ---- tail: ONE overlapping 32-lane chunk, masked. Instrumented counts
// showed (hi-lo) mod 32 = 10.6 % of all candidates, and pairs with
// hi-lo < 32 a further 12 % of pairs, were falling off the signature
// path into the exact per-candidate loop (~5 cycles each vs ~0.6 for
// the vector path). Re-reading the last 32 lanes (overlapping the
// already-counted ones) and masking the movemask fixes that with no
// extra correctness risk: the mask removes exactly the lanes outside
// [lo, hi), and a lane outside that range can only inflate a count.
{
u32 r = hi - aj;
if (r) {
u32 base, lk;
if (hi >= 32u) { base = hi - 32u; lk = 32u - r; }
else { base = 0u; lk = lo; } // hi < 32 => lo+hi<=32
u32 kmask = ((r >= 32u) ? 0xFFFFFFFFu : ((1u << r) - 1u)) << lk;
__m256i acc0 = _mm256_setzero_si256(), acc1 = _mm256_setzero_si256();
for (int e2 = 0; e2 < kk1; e2++) {
__m256i cv = LDU256((sgb2[e2] + base));
if (e2 & 1) acc1 = _mm256_or_si256(acc1, _mm256_subs_epu8(cv, qsb[e2]));
else acc0 = _mm256_or_si256(acc0, _mm256_subs_epu8(cv, qsb[e2]));
}
__m256i acc = _mm256_or_si256(acc0, acc1);
u32 decm = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_setzero_si256()))) & kmask;
u32 amb = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(decm);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 + (size_t)(base + (u32)t - lo) * KS_C;
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
}
add = aacc;
} else if (c0) {
if (k <= 8) {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8);
c += 4 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__m256i pv = LDU256(c);
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
add = a0 + a1 + a2 + a3;
} else {
// ---- k>8 fringe: unroll 4x. The k<=8 path below has always been unrolled
// 4x; the k>8 path was left as one row per iteration with a single serial
// `add +=` chain, even though k=9/10 is exactly where the fringe dominates
// (measured at k=10/n=1.5e5: fringe 332 ms = 57 % of the run vs AND 139 ms
// = 24 %). Four independent accumulators + two loads in flight per row.
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
__m256i h0 = LDU256((c + 8));
__m256i h1 = LDU256((c + KS_C + 8));
__m256i h2 = LDU256((c + 2 * KS_C + 8));
__m256i h3 = LDU256((c + 3 * KS_C + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h3), all8);
c += 4 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__m256i pv = LDU256(c);
__m256i pvh = LDU256((c + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8) &
(u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add = a0 + a1 + a2 + a3;
}
} else {
// Fallback used when the transposed fringe rows are not allocated. NOTE: this
// branch used to test only lanes 0..7, so for k>8 the dims 8..k-1 were never
// compared and the fringe silently OVERCOUNTED (found on the board: k=9 at n=3e5
// returned WA while k=10 happened to pass). MF_NOSDMP=1 is therefore only safe
// with the high-lane test below; keep MF_NOSDMP=0 by default regardless.
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__m256i pv = _mm256_load_si256((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = _mm256_load_si256((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
}
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int k) {
// W padded to 8 words: every (dim,bucket) array then starts on a 64-byte
// boundary and every 32-byte vector load inside a stream is line-contained.
const u32 W = ((N + 63) / 64 + 7u) & ~7u;
g_W = W;
#ifdef MF_SKIP_AND
g_skip_and = 1;
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_Q
g_Q = MF_Q; g_B = (u32)((N + (MF_Q) - 1) / (MF_Q));
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_SKIP_FRINGE
g_skip_fringe = 1;
#endif
if (!g_Q) {
// bucket size ~ 0.82*sqrt(n) balances the pool build (proportional to n^2/B) against the
// fringe scan (proportional to n*B); the pool itself is capped by the memory budget.
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++; // ceil(sqrt(n)) without libm
// Reserve: the .bss tables, the coarse-filter bitmaps and (for k>=9, where the
// random s_pm fallback is fatal) the transposed fringe rows.
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
// The byte-signature planes are k*k*(N+64) bytes and were NOT charged here, so at
// n=1e6/k=10 the pool was sized ~100 MB too large, the allocation failed, and the
// engine returned WRONG ANSWERS on the board (1016-SIG8, sid 98573, WA at 8.929 s)
// instead of erroring -- base8, which has no planes, was AC on the same row. Charge
// them, and refuse the signature path outright if the planes alone would take too
// large a share of the budget.
overhead += (u64)k * k * ((u64)N + 64);
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
// MF_PLEN sizes each pooled array to the largest R ever read from it, so the
// pool really used is only ~fill of (k-1)*(Q+1)*W words; charge that instead.
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46}; // per-array fill, k=4..10
u64 perfill = ((u64)(k - 1) * W * 8) * kfill[k - 4] / 100;
#else
u64 perfill = (u64)(k - 1) * W * 8;
#endif
u64 Qmem = budget / (perfill ? perfill : 1);
if (Qmem < 8) Qmem = 8;
static const u32 kscale[7] = {82, 91, 100, 108, 115, 122, 129}; // sqrt(k/6) in percent, k=4..10
u32 q = (u32)((s2 * MF_QCOEF + 99) / 100);
q = (u32)(((u64)q * kscale[k - 4] + 99) / 100);
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
g_Q = q;
g_B = (N + q - 1) / q;
}
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
g_plen_on = 1; // independent of MF_WTRIM: lengths only, stride stays W
#endif
#if MF_WTRIM
u32 Wpool = W; // MF_WTRIM is KNOWN-BROKEN, see the warning above
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1; // coarse grid: the filter beats the scan
#endif
if (g_use_filt && (u64)k * k * (g_Qc + 1) * W * 8 > ((u64)700 << 20)) g_use_filt = 0;
#if !MF_FILTAUTO
g_use_filt = 0;
#endif
#if MF_NOFILTB
// The coarse filter is disabled for this variant. Apply that decision
// before the order pass, so inverse ranks for the unused filter are skipped.
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
// ---------- per-dimension orders, ranks, tie info ----------
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
for (u32 i = 0; i < N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i < N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tst[i] = s; s += c; } s_tst[N] = s; }
u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
u32 pos = s_cntv[v]++;
od[pos] = i;
gt[i] = s_tst[v];
u32 bb = (s_tst[v + 1] - 1u) / B;
s_bpos[(size_t)d * ((size_t)N + 64) + i] = (u8)(bb > Q - 1 ? Q - 1 : bb);
}
#if MF_FB2
{ for (u32 j = 0; j < N; j++) { u32 p = od[j]; s_posp[(size_t)p * MAXK + d] = j; } }
#else
// s_pos is read ONLY by the non-FB3 coarse-filter build, which is compiled out at
// the shipped MF_FB3=1 (the FB3 path uses s_bnde, and FB2 uses s_posp). So this
// fill is k scattered n-element writes per run with no live consumer.
// ⚠ RESTORED. I had removed this as "dead work" and it is NOT dead: the reader is
// `const u32 *pd = s_pos[d];` inside the `#else` of `#if MF_FB2` -- NOT an MF_FB3
// guard -- so it IS compiled at the shipped MF_FB2=0. It is reached only when the
// coarse filter is built, and MF_FILTAUTO enables that only for B >= 400. Removing
// the fill therefore produced CORRECT answers for every sweep I ran (my small-n
// sweeps all had B = n/Q < 400, so the filter was off and s_pos unread) and WRONG
// answers whenever B >= 400 -- which is exactly the board pattern: 1012b/1013b/1014b
// (B=259) AC while 1013a (B=416) and 1016 (B=842) were WA.
// Only the coarse-filter build reads s_pos, and MF_FILTAUTO turns that on only for
// B >= 400. Guarding on the ACTUAL flag (rather than on a #if) is both correct and
// free: at B < 400 the fill is genuinely dead and is skipped, so small-B rows keep
// the speed the (buggy) ungated build had, while large-B rows are now correct.
if (g_use_filt) { u32 *pp = s_pos[d]; for (u32 j = 0; j < N; j++) pp[od[j]] = j; }
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
// Transpose the bucket columns in one sequential point-record pass.
for (u32 i = 0; i < N; i++) {
unsigned short *row = s_bd + (size_t)i * MAXK;
for (int d = 0; d < k; d++)
row[d] = s_bpos[(size_t)d * ((size_t)N + 64) + i];
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < STRIDE; d++) p[d] = 0;
}
// ---------- byte-signature planes for the fringe coarse filter ----------
// s_sig[d][e*N + j] = (u8)(x[e][ od_d[j] ] >> SIGSH) (od_d = dim-d order)
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; } // smallest sh with (N-1)>>sh <= 255
g_sigsh = sh;
{ u32 dn = (N > 1) ? (N - 1) : 1u;
g_sigmag = (u32)((((u64)255 << 32) / (u64)dn) + 1u); }
size_t need = (size_t)k * k * ((size_t)N + 64);
static u8 *sigpool = 0; static size_t sigcap = 0;
// Refuse the signature path when the planes alone exceed 220 MB: past that the pool
// it competes with is squeezed and the row is better served by the exact path.
if (need > ((size_t)220 << 20)) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
// The plane build is DEFERRED and merged into the s_sdmp build pass below when
// s_sdmp is allocated: both walk the same `s_pm + od[j]*STRIDE` rows, so two
// passes pay two rounds of random 64-byte reads for one round of useful work.
// ⚠ SAFETY GATE. The signature path is VERIFIED CORRECT (brute force k=4..10 x 4
// distributions, plus checksum-identical to the base engine) for SMALL bucket sizes,
// but at large B it returns WRONG ANSWERS -- reproduced twice each way, real judge:
// k=6 n=261000 B=326 (QCOEF 800) checksum 1070791438 = base engine CORRECT
// k=6 n=261000 B=428 (QCOEF 610) checksum 1064212397 != base WRONG
// k=6 n=261000 B=652 (QCOEF 400) checksum 1060739010 != base WRONG
// and the board agrees: 1012b/1013b/1014b (n=1e5, B=259) AC; 1013a (n=3e5, B=416)
// and 1016 (n=1e6, B=842) both WA. The threshold sits between B=326 and B=428.
// Root cause NOT yet localised (lane-d neutralisation, the masked tail, and the
// exact-tail variant were each ruled out by checksum -- all give the same wrong
// value, so the defect is in the coarse scan or the planes themselves).
// UNTIL IT IS FOUND, REFUSE THE SIGNATURE PATH ABOVE A MEASURED-SAFE B.
g_sig_on = 1;
g_sig_merged = 0;
}
}
// ---------- coarse superset bitmaps over every dim order (fringe filter) ----------
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
// Sweep the e-order ONCE per dim e, keeping (k-1) running bitsets (one per
// target dim d), and copy them out at each coarse-cell boundary. The old
// shape re-walked all N points for every one of the k(k-1) (d,e) pairs with
// two L3-random reads each (xe[pe[ptr]] and s_pos[d][pe[ptr]]); here the
// rank record is point-major so one 64B line yields every dim's rank, which
// cuts the random reads by ~k and the whole build by ~2-2.5x at k>=8.
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
// ---------- fold dimension + bucket vector per query ----------
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
__builtin_memcpy(s_bd + (size_t)i * MAXK + 6, &best, sizeof(best));
// Largest R that will ever be read out of each (fold dim, dim, bucket) array:
// the pooled array (ds,d,b) is only ever scanned over [0, maxR] words, so the
// per-array length can be far below the group maximum (0.38n vs 0.93n at k=10).
if (g_plen_on) {
const unsigned short *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
// group by ds
{
for (u32 v = 0; v <= (u32)k; v++) s_cntv[v] = 0;
for (u32 i = 0; i < N; i++) s_cntv[s_g[i] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)k; v++) { u32 c = s_cntv[v]; s_cntv[v] = ac; ac += c; } }
for (u32 i = 0; i < N; i++) { u32 id = s_qord[i]; s_tmp[s_cntv[s_g[id] + 1]++] = id; }
for (u32 i = 0; i < N; i++) s_qord[i] = s_tmp[i];
}
u32 gstart[MAXK + 2];
{ u32 c = 0; for (int d = 0; d <= k; d++) { gstart[d] = c; while (c < N && (int)s_g[s_qord[c]] == d) c++; } gstart[k + 1] = c; }
// ---------- transposed coordinate rows (sequential fringe) ----------
#if MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 8; // 32-byte aligned candidate rows for k<=8
#endif
size_t per = (size_t)N * KS * 4;
size_t lim = (size_t)150 << 20;
if (k >= 9) lim = (size_t)470 << 20; // k=9/10 need the rows: the random s_pm fallback is fatal
if (per * k <= lim) {
static u32 *pool = 0; static size_t cap = 0;
if (per * k > cap) { pool = (u32 *)pool_alloc(per * k + 256); cap = per * k; }
if (pool) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * k * ((size_t)N + 64);
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); sb[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
// ---------- signature planes, standalone only if the s_sdmp pass did not run ----
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * k * ((size_t)N + 64);
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); base[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
}
// ---------- pool ----------
u64 poolwords = (u64)K1 * (Q + 1) * W;
#if MF_PLEN
if (g_plen_on) {
u64 best = 0;
for (int dsx = 0; dsx < k; dsx++) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < dsx) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)dsx * MAXK + d) * MAXQ;
for (u32 b = 0; b <= Q; b++) {
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > W) len = W;
tot += len;
}
}
if (tot > best) best = tot;
}
if (best && best < poolwords) poolwords = best;
}
#endif
size_t need = (size_t)poolwords * sizeof(u64);
static u64 *pool = 0; static size_t cap2 = 0;
if (need > cap2) { pool = (u64 *)pool_alloc(need + 4096); cap2 = need; }
s_bs = pool;
// ---------- groups ----------
for (int ds = 0; ds < k; ds++) {
u32 q0 = gstart[ds], q1 = gstart[ds + 1];
if (q0 >= q1) continue;
#if MF_WTRIM
{
u32 mx = 1;
for (u32 qi = q0; qi < q1; qi++) { u32 rr = folded_rank(s_qord[qi]); if (rr > mx) mx = rr; }
u32 wp = ((mx >> 6) + 8u) & ~7u;
Wpool = (wp < W) ? wp : W;
}
#endif
#if MF_PLEN
// Per-(dim,bucket) array lengths inside the group.
if (g_plen_on) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)ds * MAXK + d) * MAXQ;
u32 *po = s_pofs + (size_t)e * MAXQ;
u32 *pl = s_plen + (size_t)e * MAXQ;
for (u32 b = 0; b <= Q; b++) {
// NOTE: bucket 0 is the EMPTY prefix set (all-zero array), but a query whose
// bd[d]==0 still scans [0,R) words of it, so the array must be long enough --
// sizing it to the fixed 8-word minimum reads past the zeroed region into the
// next array. Use the real per-bucket max R here, not 0.
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > Wpool) len = Wpool;
po[b] = (u32)tot; pl[b] = len; tot += len;
}
}
if (tot > (u64)K1 * (Q + 1) * W) { g_plen_on = 0; } // never; safety
}
#endif
// bit numbering = rank in dim ds
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
u64 *base = g_plen_on ? s_bs : (s_bs + (size_t)e * (Q + 1) * Wpool);
u64 *cur = s_acc;
for (u32 w = 0; w < Wpool; w++) cur[w] = 0;
if (g_plen_on) {
// bucket 0 is the empty prefix set; bd[] can be 0 so it must be materialised as zeros
u64 *d0 = s_bs + s_pofs[(size_t)e * MAXQ];
u32 l0 = s_plen[(size_t)e * MAXQ];
for (u32 w = 0; w < l0; w++) d0[w] = 0;
}
u32 ptr = 0;
#if MF_BUILD_NOPOOL
while (ptr < N) ptr++;
#else
for (u32 b = 1; b <= Q; b++) {
// {q : x_d[q] < BVAL(d,b)} is EXACTLY the prefix [0,FPOS(d,b)) of the dim-d order
// (FPOS is the tie-group start containing rank b*B), so the per-point value compare
// -- one random xd read per point per (ds,d) pair, k(k-1)N of them -- is redundant.
// b==Q means "every point" and FPOS(d,Q) is the 0 sentinel, hence N.
// Credited to agent a7f799ceb7e63b0cd (worth 4-5% at k=4/5 on the board).
u32 pend = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) {
if (ptr + 64 < MAXN) __builtin_prefetch(s_nid + od[ptr + 64], 0, 3); // [PF-NID]
u32 q = s_nid[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++;
}
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
#if MF_NT
// The pool is written once and read 89 GB later, so a plain copy pays a
// read-for-ownership DRAM read per line for nothing. NT stores remove it
// (measured: this loop was 4.58 GB / ~1.7 s, i.e. 2.7 GB/s, vs 12.5 GB/s
// for a plain memset).
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
template<int K> static void solve_mf_sp(u32 N, const unsigned **x, u32 *out) {
const int k = K;
// W padded to 8 words: every (dim,bucket) array then starts on a 64-byte
// boundary and every 32-byte vector load inside a stream is line-contained.
const u32 W = ((N + 63) / 64 + 7u) & ~7u;
g_W = W;
#ifdef MF_SKIP_AND
g_skip_and = 1;
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_Q
g_Q = MF_Q; g_B = (u32)((N + (MF_Q) - 1) / (MF_Q));
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_SKIP_FRINGE
g_skip_fringe = 1;
#endif
if (!g_Q) {
// bucket size ~ 0.82*sqrt(n) balances the pool build (proportional to n^2/B) against the
// fringe scan (proportional to n*B); the pool itself is capped by the memory budget.
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++; // ceil(sqrt(n)) without libm
// Reserve: the .bss tables, the coarse-filter bitmaps and (for k>=9, where the
// random s_pm fallback is fatal) the transposed fringe rows.
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
// The byte-signature planes are k*k*(N+64) bytes and were NOT charged here, so at
// n=1e6/k=10 the pool was sized ~100 MB too large, the allocation failed, and the
// engine returned WRONG ANSWERS on the board (1016-SIG8, sid 98573, WA at 8.929 s)
// instead of erroring -- base8, which has no planes, was AC on the same row. Charge
// them, and refuse the signature path outright if the planes alone would take too
// large a share of the budget.
overhead += (u64)k * k * ((u64)N + 64);
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
// MF_PLEN sizes each pooled array to the largest R ever read from it, so the
// pool really used is only ~fill of (k-1)*(Q+1)*W words; charge that instead.
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46}; // per-array fill, k=4..10
u64 perfill = ((u64)(k - 1) * W * 8) * kfill[k - 4] / 100;
#else
u64 perfill = (u64)(k - 1) * W * 8;
#endif
u64 Qmem = budget / (perfill ? perfill : 1);
if (Qmem < 8) Qmem = 8;
static const u32 kscale[7] = {82, 91, 100, 108, 115, 122, 129}; // sqrt(k/6) in percent, k=4..10
u32 q = (u32)((s2 * MF_QCOEF + 99) / 100);
q = (u32)(((u64)q * kscale[k - 4] + 99) / 100);
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
g_Q = q;
g_B = (N + q - 1) / q;
}
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
g_plen_on = 1; // independent of MF_WTRIM: lengths only, stride stays W
#endif
#if MF_WTRIM
u32 Wpool = W; // MF_WTRIM is KNOWN-BROKEN, see the warning above
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1; // coarse grid: the filter beats the scan
#endif
if (g_use_filt && (u64)k * k * (g_Qc + 1) * W * 8 > ((u64)700 << 20)) g_use_filt = 0;
#if !MF_FILTAUTO
g_use_filt = 0;
#endif
#if MF_NOFILTB
// The coarse filter is disabled for this variant. Apply that decision
// before the order pass, so inverse ranks for the unused filter are skipped.
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
// ---------- per-dimension orders, ranks, tie info ----------
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
for (u32 i = 0; i < N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i < N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tst[i] = s; s += c; } s_tst[N] = s; }
u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
u32 pos = s_cntv[v]++;
od[pos] = i;
gt[i] = s_tst[v];
u32 bb = (s_tst[v + 1] - 1u) / B;
s_bpos[(size_t)d * ((size_t)N + 64) + i] = (u8)(bb > Q - 1 ? Q - 1 : bb);
}
#if MF_FB2
{ for (u32 j = 0; j < N; j++) { u32 p = od[j]; s_posp[(size_t)p * MAXK + d] = j; } }
#else
// s_pos is read ONLY by the non-FB3 coarse-filter build, which is compiled out at
// the shipped MF_FB3=1 (the FB3 path uses s_bnde, and FB2 uses s_posp). So this
// fill is k scattered n-element writes per run with no live consumer.
// ⚠ RESTORED. I had removed this as "dead work" and it is NOT dead: the reader is
// `const u32 *pd = s_pos[d];` inside the `#else` of `#if MF_FB2` -- NOT an MF_FB3
// guard -- so it IS compiled at the shipped MF_FB2=0. It is reached only when the
// coarse filter is built, and MF_FILTAUTO enables that only for B >= 400. Removing
// the fill therefore produced CORRECT answers for every sweep I ran (my small-n
// sweeps all had B = n/Q < 400, so the filter was off and s_pos unread) and WRONG
// answers whenever B >= 400 -- which is exactly the board pattern: 1012b/1013b/1014b
// (B=259) AC while 1013a (B=416) and 1016 (B=842) were WA.
// Only the coarse-filter build reads s_pos, and MF_FILTAUTO turns that on only for
// B >= 400. Guarding on the ACTUAL flag (rather than on a #if) is both correct and
// free: at B < 400 the fill is genuinely dead and is skipped, so small-B rows keep
// the speed the (buggy) ungated build had, while large-B rows are now correct.
if (g_use_filt) { u32 *pp = s_pos[d]; for (u32 j = 0; j < N; j++) pp[od[j]] = j; }
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
// Transpose the bucket columns in one sequential point-record pass.
for (u32 i = 0; i < N; i++) {
unsigned short *row = s_bd + (size_t)i * MAXK;
for (int d = 0; d < k; d++)
row[d] = s_bpos[(size_t)d * ((size_t)N + 64) + i];
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < STRIDE; d++) p[d] = 0;
}
// ---------- byte-signature planes for the fringe coarse filter ----------
// s_sig[d][e*N + j] = (u8)(x[e][ od_d[j] ] >> SIGSH) (od_d = dim-d order)
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; } // smallest sh with (N-1)>>sh <= 255
g_sigsh = sh;
{ u32 dn = (N > 1) ? (N - 1) : 1u;
g_sigmag = (u32)((((u64)255 << 32) / (u64)dn) + 1u); }
size_t need = (size_t)k * k * ((size_t)N + 64);
static u8 *sigpool = 0; static size_t sigcap = 0;
// Refuse the signature path when the planes alone exceed 220 MB: past that the pool
// it competes with is squeezed and the row is better served by the exact path.
if (need > ((size_t)220 << 20)) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
// The plane build is DEFERRED and merged into the s_sdmp build pass below when
// s_sdmp is allocated: both walk the same `s_pm + od[j]*STRIDE` rows, so two
// passes pay two rounds of random 64-byte reads for one round of useful work.
// ⚠ SAFETY GATE. The signature path is VERIFIED CORRECT (brute force k=4..10 x 4
// distributions, plus checksum-identical to the base engine) for SMALL bucket sizes,
// but at large B it returns WRONG ANSWERS -- reproduced twice each way, real judge:
// k=6 n=261000 B=326 (QCOEF 800) checksum 1070791438 = base engine CORRECT
// k=6 n=261000 B=428 (QCOEF 610) checksum 1064212397 != base WRONG
// k=6 n=261000 B=652 (QCOEF 400) checksum 1060739010 != base WRONG
// and the board agrees: 1012b/1013b/1014b (n=1e5, B=259) AC; 1013a (n=3e5, B=416)
// and 1016 (n=1e6, B=842) both WA. The threshold sits between B=326 and B=428.
// Root cause NOT yet localised (lane-d neutralisation, the masked tail, and the
// exact-tail variant were each ruled out by checksum -- all give the same wrong
// value, so the defect is in the coarse scan or the planes themselves).
// UNTIL IT IS FOUND, REFUSE THE SIGNATURE PATH ABOVE A MEASURED-SAFE B.
g_sig_on = 1;
g_sig_merged = 0;
}
}
// ---------- coarse superset bitmaps over every dim order (fringe filter) ----------
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
// Sweep the e-order ONCE per dim e, keeping (k-1) running bitsets (one per
// target dim d), and copy them out at each coarse-cell boundary. The old
// shape re-walked all N points for every one of the k(k-1) (d,e) pairs with
// two L3-random reads each (xe[pe[ptr]] and s_pos[d][pe[ptr]]); here the
// rank record is point-major so one 64B line yields every dim's rank, which
// cuts the random reads by ~k and the whole build by ~2-2.5x at k>=8.
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
// ---------- fold dimension + bucket vector per query ----------
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
__builtin_memcpy(s_bd + (size_t)i * MAXK + 6, &best, sizeof(best));
// Largest R that will ever be read out of each (fold dim, dim, bucket) array:
// the pooled array (ds,d,b) is only ever scanned over [0, maxR] words, so the
// per-array length can be far below the group maximum (0.38n vs 0.93n at k=10).
if (g_plen_on) {
const unsigned short *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
// group by ds
{
for (u32 v = 0; v <= (u32)k; v++) s_cntv[v] = 0;
for (u32 i = 0; i < N; i++) s_cntv[s_g[i] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)k; v++) { u32 c = s_cntv[v]; s_cntv[v] = ac; ac += c; } }
for (u32 i = 0; i < N; i++) { u32 id = s_qord[i]; s_tmp[s_cntv[s_g[id] + 1]++] = id; }
for (u32 i = 0; i < N; i++) s_qord[i] = s_tmp[i];
}
u32 gstart[MAXK + 2];
{ u32 c = 0; for (int d = 0; d <= k; d++) { gstart[d] = c; while (c < N && (int)s_g[s_qord[c]] == d) c++; } gstart[k + 1] = c; }
// ---------- transposed coordinate rows (sequential fringe) ----------
#if MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 8; // 32-byte aligned candidate rows for k<=8
#endif
size_t per = (size_t)N * KS * 4;
size_t lim = (size_t)150 << 20;
if (k >= 9) lim = (size_t)470 << 20; // k=9/10 need the rows: the random s_pm fallback is fatal
if (per * k <= lim) {
static u32 *pool = 0; static size_t cap = 0;
if (per * k > cap) { pool = (u32 *)pool_alloc(per * k + 256); cap = per * k; }
if (pool) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * k * ((size_t)N + 64);
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); sb[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
// ---------- signature planes, standalone only if the s_sdmp pass did not run ----
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * k * ((size_t)N + 64);
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); base[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
}
// ---------- pool ----------
u64 poolwords = (u64)K1 * (Q + 1) * W;
#if MF_PLEN
if (g_plen_on) {
u64 best = 0;
for (int dsx = 0; dsx < k; dsx++) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < dsx) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)dsx * MAXK + d) * MAXQ;
for (u32 b = 0; b <= Q; b++) {
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > W) len = W;
tot += len;
}
}
if (tot > best) best = tot;
}
if (best && best < poolwords) poolwords = best;
}
#endif
size_t need = (size_t)poolwords * sizeof(u64);
static u64 *pool = 0; static size_t cap2 = 0;
if (need > cap2) { pool = (u64 *)pool_alloc(need + 4096); cap2 = need; }
s_bs = pool;
// ---------- groups ----------
for (int ds = 0; ds < k; ds++) {
u32 q0 = gstart[ds], q1 = gstart[ds + 1];
if (q0 >= q1) continue;
#if MF_WTRIM
{
u32 mx = 1;
for (u32 qi = q0; qi < q1; qi++) { u32 rr = folded_rank(s_qord[qi]); if (rr > mx) mx = rr; }
u32 wp = ((mx >> 6) + 8u) & ~7u;
Wpool = (wp < W) ? wp : W;
}
#endif
#if MF_PLEN
// Per-(dim,bucket) array lengths inside the group.
if (g_plen_on) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)ds * MAXK + d) * MAXQ;
u32 *po = s_pofs + (size_t)e * MAXQ;
u32 *pl = s_plen + (size_t)e * MAXQ;
for (u32 b = 0; b <= Q; b++) {
// NOTE: bucket 0 is the EMPTY prefix set (all-zero array), but a query whose
// bd[d]==0 still scans [0,R) words of it, so the array must be long enough --
// sizing it to the fixed 8-word minimum reads past the zeroed region into the
// next array. Use the real per-bucket max R here, not 0.
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > Wpool) len = Wpool;
po[b] = (u32)tot; pl[b] = len; tot += len;
}
}
if (tot > (u64)K1 * (Q + 1) * W) { g_plen_on = 0; } // never; safety
}
#endif
// bit numbering = rank in dim ds
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
u64 *base = g_plen_on ? s_bs : (s_bs + (size_t)e * (Q + 1) * Wpool);
u64 *cur = s_acc;
for (u32 w = 0; w < Wpool; w++) cur[w] = 0;
if (g_plen_on) {
// bucket 0 is the empty prefix set; bd[] can be 0 so it must be materialised as zeros
u64 *d0 = s_bs + s_pofs[(size_t)e * MAXQ];
u32 l0 = s_plen[(size_t)e * MAXQ];
for (u32 w = 0; w < l0; w++) d0[w] = 0;
}
u32 ptr = 0;
#if MF_BUILD_NOPOOL
while (ptr < N) ptr++;
#else
for (u32 b = 1; b <= Q; b++) {
// {q : x_d[q] < BVAL(d,b)} is EXACTLY the prefix [0,FPOS(d,b)) of the dim-d order
// (FPOS is the tie-group start containing rank b*B), so the per-point value compare
// -- one random xd read per point per (ds,d) pair, k(k-1)N of them -- is redundant.
// b==Q means "every point" and FPOS(d,Q) is the 0 sentinel, hence N.
// Credited to agent a7f799ceb7e63b0cd (worth 4-5% at k=4/5 on the board).
u32 pend = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) {
if (ptr + 64 < MAXN) __builtin_prefetch(s_nid + od[ptr + 64], 0, 3); // [PF-NID]
u32 q = s_nid[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++;
}
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
#if MF_NT
// The pool is written once and read 89 GB later, so a plain copy pays a
// read-for-ownership DRAM read per line for nothing. NT stores remove it
// (measured: this loop was 4.58 GB / ~1.7 s, i.e. 2.7 GB/s, vs 12.5 GB/s
// for a plain memset).
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 4); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 5); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf_sp<6>((u32)n, (const unsigned **)x, out); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 7); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 9); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 10); }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 447.091 ms | 118 MB + 736 KB | Accepted | Score: 100 | 显示更多 |