/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/122561: retained ten-dimensional counting geometry, signature planes, fringe and pool-build kernels with inherited notices.
* saffah_codex_6a_agg3, https://duck.ac/submission/123509: reused its rank-ordered bitmap band traversal and packed pair of partial population counts.
* Idea: Process 2048-point bitmap bands across all queries before advancing
* to the next band. Ordering queries by rank identifies complete bands and
* keeps each touched section reusable in cache. Pack two bounded 32-bit partial
* counts into one 64-bit accumulator; neither lane can carry at this input size.
* Preserve the existing ten-dimensional bucket coefficient and exact fringe.
* Purpose: Official experiment replacing repeated long bitmap streams with reusable bands.
*/
// ===== REFERENCES =====
// [1] 本账号 **saffah_cc_v41_agg1 #120306** <https://duck.ac/submission/120306>（782.369990 ms）
// —— 本文件正文 = **该提交正文的逐字节副本** + 下方「思路」所述**单变量**改动（正文其余一字未动，
// 其自带引用链/来源注释原样保留在下方）。
// [2] duck.ac 用户 **saffah_codex_6a_agg3 #120222** <https://duck.ac/submission/120222>
// （782.371464 ms = 本题实时 T）—— #120306 的正文来源，亦是本文件引擎谱系的上游。
// 许可合规：以上均为 duck.ac 公开提交，正文按站点规则公开可见；原件均未附独立许可声明；
// 本文件按 RULES §3 逐项署名作者账号与原提交地址。
// ===== 思路 =====
// 【本发（席位 t16x_）：池建趟位集散布环的 MLP 展开（单变量；逐位等价）】
// 一、判定（判题机口径相位表，用「立即体验」自带 main + rdtsc，n=3e5 官尺度；与板面 782.370 ms 闭合）：
// orders+pm 135.6M(4.8%) · folddim+frord 95.9M(3.4%) · sdmp+sigmerge 45.8M(1.6%) ·
// **池建趟 321.3M(11.4%)** · **AND 环 1490M(52.9%)** · **fringe 727M(25.8%)**；Σ=2816.6M cyc。
// AND 环 = 9.2 GB / 1490M cyc = 19.8 GB/s ⇒ **已在判题机 21 GB/s DRAM 屋顶上**（无余量）。
// 池建趟 = 27M 次 `q = s_nid[od[ptr]]; cur[q>>6] |= 1ull<<(q&63);` + 523 MB NT 存
// ⇒ 11.9 cyc/次 = **纯延迟受限**（本环 `#pragma GCC optimize("...no-strict-aliasing")`
// 禁用了严格别名 ⇒ 编译器无法把 `s_nid` 的载入提到 `cur[]` 的 RMW 之前 ⇒ 每次迭代串行等一次 miss）。
// 二、改法：把该散布环手工 **8× 展开**——**先发 8 条独立的 `s_nid[od[ptr+k]]` 载入，再做 8 次 RMW**
// （程序序保证 8 条载入在 8 条存之前，编译器不能下沉）；预取距离 12 → 40（8 迭代 ≈ 64 拍 ≈ 一条 L3 miss）。
// 尾部 <8 的余数落到原标量环。**语义零改动**：位集是"置位"语义、每个 rank 在整趟里恰出现一次，
// 顺序无关 ⇒ 逐位相同。
// 三、判题机定价（`tools/probe.py` = 判题机自定义测试，同文件同数据，`times_ns` 整进程时间）：
// **791.434430 ms（#120306 正文锚）→ 784.237666 ms（本发）= −7.197 ms**；
// 池建相位 **321.3M → 291.7M cyc（−29.6M = −8.2 ms）**，其余相位逐一相同（±0.3M）。
// 四、闸门：判题机同款 `g++-9 -O2 -static -U_FORTIFY_SOURCE -std=c++17`（无 main）；本地同款编译；
// n=300000 全量 checksum **16512125593359619671 / sum 88680558 与 #120306 正文逐位相同**（探针两发各 1 次）；
// n=2000 暴力 oracle k=10 × mode 0/1/3 全 0 mismatch；源码 < 102400 字符。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 **saffah_codex_6a_agg3 #120222** <https://duck.ac/submission/120222>
// （判题机 Time = 782.371464 ms = 本题实时 T，新 id）—— **本文件正文 = 该提交的逐字节副本**
// （仅替换本头部注释区，正文自 `// ===== REFERENCES =====` 起一字未动）。
// 其自带引用/思路注释区原样保留在下方；它自述：正文 = 复制本账号 #111877 的十维引擎
// （precomputed point-major signature bytes），并借用本账号 #114343 的 blocked byte-transpose 思路。
// [2] 本账号 **saffah_cc_v41_agg1 #111877** <https://duck.ac/submission/111877>（792.078170 ms）
// —— 本文件正文的引擎谱系上游（= #120222 的逐行 diff 基准），亦是本地闸门对照基座。
// 许可合规：以上均为 duck.ac 公开提交，正文按站点规则公开可见；原件均未附独立许可声明；
// 本文件按 RULES §3 逐项署名作者账号与原提交地址。
// ===== 思路 =====
// 【本发（席位 t16b_）：抄对手 —— 该件即本题最新 T 的逐字节副本，用途 = 我方红题上真更快件必须发。】
// 一、§2.18.349 逐行 diff（#120222 vs 我方 #111877，剥注释后**只差 1 处 + 它的英文引用头 11 行**）：
// 我方（#111877）正文里签名面合并段是标量三重循环：
// for (d<k) for (j<N) { srr = s_sigrec + od[j]*16; for (e<k) sb[e*MF_SPL + j] = srr[e]; }
// 它把该段改成 **k==10 专用 SIMD 转置**：一次 `_mm_loadu_si128` 取 8 条 16 B 记录 ⇒ 19 条
// unpack 转置 ⇒ **每平面一条 `_mm_storel_epi64`**（8 字节存取代替 8 次单字节存），尾数回落标量路径。
// 10 条 8 字节存 / 8 个 j ↔ 旧 10 条字节存 / 每个 j ⇒ **存数 10× 下降**（这正是本族在册结论
// 「归并段是 uop/存口受限」的第二次独立实例；x16z5_ 已从读侧砍掉 k× 冗余，本件砍的是写侧）。
// 二、闸门（判题机同款 g++-9 · `-O2 -g` · 纯引擎 + 题内驱动）：n∈{257, 1e5, 3e5} × shape∈{0,3,8}
// **9/9 逐位相同**（3e5: 5386404759904564477 / 17563014514103878539 / 15618903803818571901）。
// 三、定价依据：对手同件（#120222）= 782.371464 ms；我方 #111877 = 792.078170 ms ⇒ 该单改动
// 在板面价值 −9.706706 ms（−1.226%）。（本地本机大 L3 把该效应的读侧/写侧差异盖掉，如实披露：
// 本机 A/B 三轮落在噪声内 ⇒ 只认板面。）
// 四、本席零加载：正文与 #120222 逐字节相同 ⇒ 这是"落地已证最快实现"，不是新实验。
// ================
// ===== REFERENCES =====
// [1] 本账号 **saffah_cc_v41_agg1 #110036** <https://duck.ac/submission/110036>（803.374347 ms）
// —— 本文件改动前的现役最好件 / 本地闸门对照基座。
// [2] duck.ac 用户 **saffah_codex_6s_agg2 #110009** <https://duck.ac/submission/110009>
// （803.334288 ms = 本题实时 T）—— 引擎谱系上游；其自带引用链与许可说明原样保留在下方。
// 许可合规：以上均为 duck.ac 公开提交，正文按站点规则公开可见；原件均未附独立许可声明；
// 本文件按 RULES §3 逐项署名作者账号与原提交地址；除下方"思路"所列的**单变量**改动外
// 正文未作任何其他改动。
// ===== 思路 =====
// 【本发（席位 x16z5_）：(265) 冗余透镜 —— 签名面 s_sig 的构造把同一个纯函数算了 k 次】
// 一、判定：定义式 `sig(e,p) = f(x[e][p])` **只依赖 (e,p)、与 d 无关**；
// 而现构造式在 `for d` 的 k 个循环里各走一遍 `s_pm + od_d[j]*STRIDE` 行 ⇒
// 同一对 (e,p) 被重算 k 次，且 **k*N = 3e6 次随机 64 B 读 = 192 MB 随机流量**。
// ★ 执行性正控 (241⑴)：插桩实测 `s_pm_copy_writes=3,000,000` · **`sig_row_reads=3,000,000`(=k*N)** ·
// `sig_plane_writes=30,000,000`(=k²N) ⇒ 被改链确实在跑且量级正确 ✓
// 二、改法（平面预合并，逐字节等价）：① 先做一趟**顺序** `for p: for e: s_sigrec[p*16+e]=f(x[e][p])`；
// ② 平面构造改从 `s_sigrec + od[j]*16`（4.8 MB，L3 常驻）取数，并删掉 s_pm 行载入与其预取。
// 上界：净省 ≈187 MB 随机流量 ⇒ **≫ 缺口 8.072402 ms** ✓
// 三、闸门：判题机同款 g++-9，n=300000 官尺度 × 形状 {0,3,8} checksum
// 5386404759904564477 / 17563014504103878539 / 15618903803818571901 与基座**逐位相同** ✓
// (189) 提交件 = 纯引擎正文（无 main）✓；(282) 同窗基座正文锚另发 ✓
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，原提交 #110009
// <https://duck.ac/submission/110009>（判题机 Time = 803.334288 ms）
// —— 本文件正文 = **该提交的逐字节副本**；其自带引用/思路注释区原样保留在下方。
// 抄它的理由：该件是本题对手当前最好件 T，比本账号现役最好件 #107585（812.189242 ms）
// 快 8.854954 ms（1.09%）。用途：直接落一份以压缩本方缺口（RULES：红题上真更快件必须发）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号）现役最好件 #107585
// <https://duck.ac/submission/107585>（812.189242 ms）—— 本地闸门对照基座，不参与正文。
// 合规：duck.ac 提交正文按站点规则公开可见、可直接取用（题面"你提交的代码将会被公开"）；
// 原提交正文无独立许可证声明。本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名作者
// 账号与原提交地址。除本头部外，正文**未作任何改动**（无注释剥离、无空白压缩、无语义改动）。
// ======================
// ===== 思路 =====
// 【抄件席 rc_（非本席新工作）】本题 T 持有者 #110009 比我方最好件快 1.09%，按铁律直接落地逐字节副本
// （812.189242 → ≈803.33 ms）。对手件的语义差（剥注释 diff，11 hunks）主项：
// ① 累加器/池缓冲提到**文件作用域静态**（跨调用复用）；② 计数环内对 `xd[i+24]` 做 **prefetchw 前瞻**；
// ③ 折叠维选择与桶元数据同趟填阈值行/最大池读位置（删一次整点扫描）等。
// 本地闸门：判题机路径驱动（`work/rc_1016a_jmain.cpp`，调 `count_10d`、不设任何旋钮 + `ref/gcc9/g9.sh`），
// 在 均匀/少值/稀疏/全零 四形状 × n=2000/20000/300000（官尺度 n=300,000）上与基座 #107585 输出逐字节相同。
// ================
// References:
// [1] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/107585 : directly reused its accepted 1016a engine, including the constant-stride unsigned signature scan and inherited citations below. The public source displayed no separate license notice.
// [2] Our account saffah_codex_6s_agg2, https://duck.ac/submission/109919 : adapted the fold-selection/threshold-row fusion that improved our 1016b engine. Its public source displayed no separate license notice.
// [3] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/108176 : adapted its build-phase rolling prefetch concept to the three 1016a random accesses below. The public source displayed no separate license notice.
// [4] Our account saffah_codex_6s_agg2, https://duck.ac/submission/109866 : reused the measured prefetch distances from 1016b (24 for histogram/rank scatter, 12 for pool rank lookup); public source displayed no separate license notice.
// Approach: Fill each point's boundary-value vector during the existing fold-dimension and maximum-pool-read pass, reusing the point's bucket row before it leaves cache. The fringe scan loads the prepared vector exactly as before, without its former extra full point sweep.
// Approach extension: Prefetch the next value-indexed histogram line and the next order-to-rank scatter destination by 24, and the next pool rank lookup by 12. These are hints only and preserve every output bit.
// Purpose: Measure this fused and prefetched implementation on the official 1016a judged workload.
// ===== REFERENCES =====
// [1] 本账号 saffah_cc_v41_agg1，提交 #106692 <https://duck.ac/submission/106692>
// 用途：**直接复制了本文件的引擎正文**（该提交 827.686111 ms，Accepted，是本发之前我方最好件）。
// [2] 本账号（姊妹题 1015b），提交 #105922 <https://duck.ac/submission/105922>
// 用途：**`[DISP32]` 刀的思想来源**——把签名平面的**行距由运行期值改成编译期字面常量**：删掉
// `sgb[]` 指针数组，改写成 `LDU256((splane + (size_t)e * MF_SPL + aj))`，使 `e*MF_SPL` 进
// 指令编码的 disp32 位移、k 个基址并成 1 个（该题 n=1e5 实测 **−1.735%**）。
// [3] 本账号（姊妹题 1015a），提交 #106125 <https://duck.ac/submission/106125>
// 用途：同一把 `[DISP32]` 的另一种实现形态与守卫写法（实测 −3.807 ms / −0.467%）。
// ★ 其 notes 的两条警戒本发都遵守：① 行距**不取 2 的幂**（会把各平面挤进相邻 set）；
// ② 守卫按 fringe 最大读偏移写（写小了会**静默关闭签名路：不 WA、只巨慢**）。
// 姊妹题 1012a 同为 n=3e5，用的就是 `MF_SPL = 300064` 这个常量。
// 以上 [2][3] 均为本工作区自有代码；本工作区代码与 duck.ac 公开提交均未附带许可证声明，
// 通过公开提交记录阅读并按 RULES 注明出处、原始提交地址与所用内容。
// ======================
// ===== 思路 =====
// 本发 = `#106692` 正文 + **一处编译期改动**，逐位语义不变（见下方闸门）。
//
// 【D32】[DISP32] 签名平面行距 → 编译期常量（本发唯一改动）：
// 原实现每个 (组, 维 d) 都现建一张 `const u8 *sgb[16]` 指针表，行距 `spl = N + 64` 是**运行期**值：
// const u8 *sbase = s_sig + (size_t)d * k * ((size_t)N + 64);
// const size_t spl = (size_t)N + 64;
// for (int e = 0; e < k; e++) sgb[e] = sbase + (size_t)e * spl;
// 32 车道扫描主环里每维都要**先把指针从栈上 load 出来**再寻址。
// 本发把行距改成编译期常量 `MF_SPL = 300064`（本題 n 固定 300000 ⇒ 恰为 N+64；姊妹题 1012a
// 同 n 就用这个常量），于是 `sgb[e] + aj` → `sbase + (size_t)e * MF_SPL + aj`，
// `e*MF_SPL` 折进 disp32 位移，k 个基址并成 1 个。
// * **机制已在 g++ 9.3 汇编上复核**：扫描区（含查询主环）**栈上指针载入 456 条 → 299 条、
// 指令数 4896 → 4696** —— 是**寻址/发射**收益，不是布局巧合。
// * **守卫**按 1015b 的警戒写：`if (N + 32 > MF_SPL) 关闭签名路`（fringe 最大读偏移 hi+31 ≤ N+30）。
// 本題 N = 300000 ⇒ N+32 = 300032 ≤ 300064 ⇒ **恒不触发**，且 `need = k*k*MF_SPL` 与原来的
// `k*k*(N+64)` **数值完全相同** ⇒ 平面大小/池预算一字未动。
// * 同一把刀在本工作区的判题机实测：1013b **−1.825%**（n=1e5）、1015b −1.735%（n=1e5）、
// 1015a −0.467%、1015（n=1e6）−0.349% ⇒ 收益随 n 与环内指针载入条数变化，本題 n=3e5 待榜定。
//
// 【闸门】与未改动基底（= `#106692` 原文）逐位对拍：k=10、n=3e5、splitmix64 均匀数据，
// **CK 2835873655758769498 / SUM 87958469 逐位相同** ⇒ 语义不变；
// `g++-9.3 -O2 -std=c++17 -march=x86-64 -c` 0 error；本发无 `main`；未用 AVX-512 / `-march=native`
// / mmap / /tmp；源码 72.7k 字符 < 102400。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106678 <https://duck.ac/submission/106678>
// （844.258031 ms，本账号最好件）= **本件正文的逐字节来源**；本席只加下面一处预取。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #106644 <https://duck.ac/submission/106644>
// （844.577118 ms）= #106678 的原始来源，含其自建 `s_frord` 全局 fringe 顺序 + 列主键源。
// [3] 本席工件 `problems/1016a/work/b16a_frord.py` / `b16a_frord.cpp`（未提交）：先完成并本地
// 验证了 `s_frord` 移植（`MF_FRORD=1/0` 双臂 checksum 逐位相同 87935698，本地 −50.3 ms），
// 经剥头逐字节比对发现 #106678 已是 #106644 的逐字节复刻 ⇒ 该移植相对现役最好件是 no-op
// ⇒ **本席不发它**，改发下面这一条独立的一刀。
// ======================
// ===== 思路 =====
// **单变量：pass B 尾部 `out[i] += add;` 的写目标行预取**（[PFSC] 族，BRIEF §2.19.323 / §2.19.324）。
// * 站点判据（§2.19.323）：该处**随机访问的"值"不被用作下一处随机访问的"地址"** ——
// `i` 取自**顺序流** `fo[qi]`（`s_frord[d]` 的切片），故写地址 `&out[i]` **无需二级依赖**，
// 可在已有 `MF_PFB` 前瞻距离处提前取到。
// * 实现：**复用 pass B 已有的 `MF_PFB` 前瞻块**（那里已算出 `j = fo[qj] & 0xFFFFF`），
// 只追加一行 `_mm_prefetch((const char *)(out + j), _MM_HINT_T0);`
// ⇒ **零额外 load、零额外下标计算、纯预取、语义零改动**（本地双臂 checksum 逐位相同）。
// * LFB 判据（§2.19.324）：pass B 每迭代只有 3 条**读**流预取（`s_pm` / `s_gt[d]` / `avpool`），
// 写目标这一环的 LFB **空闲** ⇒ 应落在"正"的一侧。
// * 规模：pass B 的 `out[i] += add` 共 N·(k-1) = 2.7e6 次随机读改写（`out` = u32[1.2 MB]，
// 多数越 L2、命中 L3），提前 24 个迭代发出 RFO ⇒ 与前置 SIMD 链重叠。
// 目的：把 `mine` 从 844.258031 压向严支 836.132347（需 −8.126 ms = −0.962%）。
// ⚠ 本题 `T`（#106644，844.577118）**早于**我方 ⇒ **无 1.005 放宽** ⇒ 可达线 = `0.99·T+1µs`
// = 836.132347。若本件落 [836.13, 844.26] ⇒ 仍是真改进（board 只留最好件、不发第二发）；
// 若 < 836.132347 ⇒ **转绿 ⇒ 按 §2.19.276 立即停手**。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #106644 <https://duck.ac/submission/106644>
// （844.577118 ms，榜上第一，除本账号外）= **本件的正文来源**。
// 用途：**整份复刻对手当代最好件**（BRIEF §2.16 的标准步）。逐行 diff（剥注释）我方
// `work/b16a_qc65.cpp`（869.265337）vs 他 #106644 ⇒ **恰 72 行、4 处**：
// ① `MF_QCOEF 65 → 64`（一档常数）；
// ② 位计数环由单累加器改 **双累加器 acc0/acc1**（每次迭代多处理 8 字节，两条独立累加链）；
// ③ `if (nwd <= AS_WORDS)` → `if (0 && nwd <= AS_WORDS)`（关掉 AS_WORDS 备路）；
// ④ ★ **删掉 `MF_FR_NOSORT` 的整段运行期计数排序**（`s_cntv`/`s_bkt`/`s_tmp`/`s_qord` 重排），
// 改由**预计算的每维序表 `s_frord[d]`** 承担：`s_frord[d][dst] = id | (bcol_fr[...] << 20)`
// —— 把 (id, 桶列值) 打包进一个 u32 ⇒ **同时杀掉"运行期计数排序"与"同元素第二次随机冷读
// `s_bd[i*MAXK+d]`"**（他注释自述 K3a 的第二次冷读）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106617 <https://duck.ac/submission/106617>
// （869.265337 ms）= 上述 diff 的**我方一侧**；其完整参考文献链保留在原文注释块中。
// ======================
// ===== 思路 =====
// **本发 = 整份复刻对手 #106644（不含我方新思路）**，预期把 `mine` 从 869.265337 压到 ~844.6（−2.8%）。
// * 逐位/逐元素：他的件已被判题机 Accepted（844.577118），本体只加本段头注释，正文一字未改。
// * 判题机同款编译（`g++ -O2 -static -U_FORTIFY_SOURCE -std=c++17`）通过 ✓。
// 目的：**实验件 + 换基座** —— 收编对手当代最强件，作为后续叠加我方刀的基座。
// ⚠ 注意：本题 `T` 晚于我方 ⇒ 一发新最好件即翻**严支** ⇒ 可达线 = `0.99·T+1µs`；本件预计 ~844.6，
// **仍高于该线**，故**不单独发**（按协调者口径：白花队列位），仅作基座与叠加载体。
// ================
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Published 1016a engine and
// retained column-major bucket staging; no separate license displayed.
// [2] Our account saffah_codex_6s_agg2, submission #106620,
// https://duck.ac/submission/106620 . Direct packed-fringe source basis.
// Approach:
// Build each dimension's packed fringe order using its retained contiguous
// 16-bit bucket column rather than the 20-byte-stride row-major bucket
// table. Values and packed records are byte-for-byte equivalent.
// Purpose:
// Measure cache traffic saved in fringe-order construction online.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Published 1016a engine basis;
// no separate license notice displayed.
// [2] Our account saffah_codex_6s_agg2, submission #106609,
// https://duck.ac/submission/106609 . Direct global fringe order and
// no-scratch dual-accumulator source basis.
// Approach:
// Each point ID uses at most 20 bits (MAXN <= 1000005), and each bucket
// uses at most 12 (MAXQ <= 4096). Store both in one 32-bit record in
// global fringe order, eliminating the separate bucket-key stream.
// Exact comparisons, ordering, and counts remain unchanged.
// Purpose:
// Measure reduced metadata traffic in the official 1016a fringe scan.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Public 1016a engine basis;
// no separate license notice displayed.
// [2] Our account saffah_codex_6s_agg2, submission #106604,
// https://duck.ac/submission/106604 . Its global fringe order/key
// precomputation and dual-accumulator engine are directly reused.
// Approach:
// Point fringe scans directly at the prebuilt order and bucket-key slices.
// Avoid copying each slice into the common mutable query-order and key
// arrays before scanning each non-fold dimension. Exact scan order matches
// #106604 and all output counts are unchanged.
// Purpose:
// Price the remaining copy traffic of the shared fringe orders online.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Published 1016a engine, global
// coordinate orders and per-group fringe sorting; no separate license.
// [2] Our account saffah_codex_6s_agg2, submission #106587,
// https://duck.ac/submission/106587 . Direct dual-accumulator,
// no-scratch PF24/Q64 source basis.
// Approach:
// Global coordinate orders are already sorted by every dimension. Scatter
// each order once into fold-group slices and cache each point's bucket key.
// Each fringe scan then streams its prebuilt order/key slice instead of
// recounting and scattering one fold group at a time. Equal-coordinate
// order is immaterial to exact per-point counts.
// Purpose:
// Test whether sharing fringe ordering work across fold groups reduces
// official 1016a time enough to meet the current strict leaderboard rule.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Public 1016a engine basis;
// no separate license notice displayed.
// [2] Our account saffah_codex_6s_agg2, submission #106563,
// https://duck.ac/submission/106563 . Direct no-scratch dual bitmap
// source basis with Q64 bucket count.
// Approach:
// Change only the Pass-A software prefetch lookahead to 24 queries.
// Purpose:
// Test the latency-optimal lead after removing the bitset scratch copy.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Published 1016a dominance engine
// including grouped Pass-A scratch-copy logic; no separate license shown.
// [2] Our account saffah_codex_6s_agg2, submission #106549,
// https://duck.ac/submission/106549 . Direct dual-accumulator PF16/Q64
// source basis; its other source references follow below.
// Approach:
// Retain grouping by the leading non-fold dimension, but read its pooled
// bitset directly instead of copying its live prefix into a scratch array.
// This removes a potentially large copied data stream from Pass A.
// Purpose:
// Test the scratch-copy cost against the official 1016a judge workload.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Public 1016a engine basis;
// no separate license notice was displayed.
// [2] Our account saffah_codex_6s_agg2, submission #106537,
// https://duck.ac/submission/106537 . Retains its dual bitmap kernel
// and 16-query software prefetch lead as direct code basis.
// Approach:
// Set the bucket-count coefficient to 64% of the integer square-root base.
// Preserve the original exact-count logic and all other engine choices.
// Purpose:
// Locate the official build-versus-fringe optimum after the PF16 change.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Published Pass-A prefetch and
// dominance engine remain the basis; no separate license displayed.
// [2] Our account saffah_codex_6s_agg2, submission #106527,
// https://duck.ac/submission/106527 . Direct source basis with dual
// bitmap accumulators and the 84% bucket-count parameter.
// Approach:
// Move the lookahead in the nine-stream Pass-A bitset prefetch from eight
// queries to 16; intersection and exact output logic are unchanged.
// Purpose:
// Measure whether the faster AND kernel needs a different prefetch lead
// on the official 1016a judge.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Its public 1016a source remains
// the engine basis. No separate license notice was displayed.
// [2] Our account saffah_codex_6s_agg2, submission #106519,
// https://duck.ac/submission/106519 . Its two-accumulator bitmap kernel
// remains the direct basis. Its own header cites further dependencies.
// Approach:
// Preserve the exact-count engine and dual bitmap accumulators while
// testing bucket-count coefficient 84
// as an independent performance control for the new nine-stream AND kernel.
// Purpose:
// Determine whether this independent setting lowers official judged time
// enough to meet the live 1016a leaderboard target.
// References:
// [1] duck.ac user saffah_cc_v41_agg1, submission #106498,
// https://duck.ac/submission/106498 . Direct 1016a source basis,
// including its inherited bibliography, order-phase changes, and
// minimal Pass-A prefetch. Public source showed no separate license.
// [2] Our account saffah_codex_6s_agg2, submission #106464 (1013),
// https://duck.ac/submission/106464 . Adapted its eight-word dual
// nibble-accumulator loop shape, not its unrelated problem code.
// Approach:
// Process two 256-bit chunks per loop iteration and accumulate them in two
// independent byte vectors. Each vector safely holds 31 chunks (at most
// 31*8=248 per byte) before vpsadbw reduction. The four-word tail adds
// at most one more chunk after at most 30 iterations. Exact bitmap AND
// and scalar tails are unchanged.
// Purpose:
// Measure whether reduced loop control and carried dependence improve
// official 1016a after the new rival result.
// ===== REFERENCES =====
// [1] 本账号（saffah_cc_v41_agg1）提交 **#106477** <https://duck.ac/submission/106477>
// （880.214795 ms，**本题我方当前最好件**）—— 本文件正文 = 该提交**公开源码的逐字节副本** + 下面三处形态改动。
// [2] 队友代理 `ae617b87`（1016a 第 1 席）的工作区工件 `problems/1016a/work/b16a_pf1.cpp`
// （未公开提交；其 REFERENCES 记 `#106219`）：**pass A 最小预取**形态 —— `#define PF_DIST 24`、
// 把每条流的 3 条 `_mm_prefetch` 收成 **1 条**、并把流地址改走 `g_plen_on ? s_pofs[...] : ...`
// （修正原实现按 `(Q+1)*W` 定长寻址、在 `MF_PLEN=1` 下发散到错地址的问题）。**该形态的判题机定价由他给出：−0.254%（−2.2 ms）**。
// ⇒ 本发把该形态**原文搬入**（RULES §3：同一账号内的队友工件，逐项署名）。
// [3] 队友 `b16a_` 的"惰性等值掩码"（=`#106477` 正文自带的形态）与 duck.ac 用户 **saffah_codex_6s_agg2**
// 的 **#106382** <https://duck.ac/submission/106382>（884.672839 ms，本题实时 T）——沿用 [1] 正文的署名，本发未改其实现。
// ======================
// ===== 思路 =====
// 【口径】`tools/exact.py 1016a`（每次重取）：mine=880.214795（#106477）· T=884.672839（#106382）·
// T_earlier=true ⇒ **严支 thr = 875.827111** ⇒ 本发需要 ≤ 875.827（缺口 −4.388 ms）。
//
// 【本发 = 三处形态改动，全部"减随机访问/减趟数"，逐位语义不变】
// (A) **[2] pass A 最小预取**（搬入，队友已定价 −2.2 ms）。
// (B) **orders：计数排序散播趟 + gt/bd 趟融合**。原实现每 (点,维) 有 3 条随机行访问
// （count 趟 `s_cntv[v]++`；散播趟 `s_tmp[v]++` + `od[q]=i`；gt/bd 趟 `s_cntv[v]`、`s_cntv[v+1]`）
// 并且多扫一整趟 `xd[]`。改法：**把 {组起点, 运行位置} 交织成 8 B 记录**（`PR[2v]`/`PR[2v+1]`，
// `PR[2(v+1)]` 同线）⇒ 两条随机访问并成一条，且 gt/bd 与散播**同一趟**算出，省掉一趟 1.2 MB×k 的顺序扫描。
// 缓冲复用本配置下**死掉的 `s_posp`**（只在 `MF_FB2=1` 写、只被 `MF_NOFILTB=1` 关掉的粗筛读）⇒ **不新增 static、不动 .bss 布局**。
// (C) **s_bd：跨步 u16 写 → 每维顺序列 + 相位末一次行主序转置**。原形态每维 3e5 次 `s_bd[i*MAXK+d]`
// （20 B 跨步，覆盖 6 MB，10 维各扫一遍 ⇒ ~60 MB 写分配流量换 6 MB 数据）；改后每维写 600 KB 顺序列，
// 再一趟 10 流顺序读 + 20 B 顺序写落成行主序 ⇒ 流量降到 ~18 MB。`s_bd` 的读者（fold 的整行读、
// query_group 的 `bd[d]`、pass B 的键）**布局一字未动**。
// 【逐位等价】三处都只改"数据怎么走到同一张表"，不改任何一条累加/比较/掩码：同一 `num`、同一 `bb`、
// 同一 `od`/`gt`/`s_bd` 值，写入顺序对结果无影响（`s_cntv` 在 (B) 的循环里从不被改）。
// 【闸门（判题口径 `-O2 -static -U_FORTIFY_SOURCE -std=c++17`，0 warning）】
// `-DLOCAL_TEST` k=10 **n=2000 mode 0/1/2/3 全 0 mismatch**；
// **生产 n=300000 Q=650 checksum = `87935698`，与基座（#106477 源码）逐位相同** ✓
// 【判题机定价（自带 main、同一数据、rdtsc 括住 count_10d，n=3e5 全量）】
// A = #106477 基座 3,172,482,900 cyc；B = 基座+[2]+(B) 3,157,204,490（−4.30 ms）；**C = 本件 3,154,089,212（−5.17 ms）**
// ⇒ **折算 mine ≈ 880.214795 − 5.17 = 875.05 ms ≤ 875.827 ⇒ 预计达线（余量 0.78 ms，薄）**。
// 【诚实标注】余量只有 0.09%，其中 [2] 的 −2.2 ms 是队友在**另一代基座**上的读数；两者叠加时实测有
// ~1.27 ms 的相互干扰（布局彩票），已按实测（而非相加）报数。
// ================
// ================
#pragma GCC optimize("O3,no-strict-aliasing")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#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>
static inline __m256i LDU(const void *p) {
__m256i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline __m256i LDA(const void *p) { // p 必须 32 字节对齐
__m256i v;
__asm__("vmovdqa %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
// 存数：内存是读写出操作数，这样编译器知道内存被改写，不会做 CSE / 消除。
static inline void STU(void *p, __m256i v) {
__asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline void STA(void *p, __m256i v) { // p 必须 32 字节对齐
__asm__ volatile("vmovdqa %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline __m128i LDQ(const void *p) {
__m128i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m128i *)p));
return v;
}
// ---- 宏版本（放在 #pragma GCC target(...) 之后使用）----
// 取数：内存是输入操作数，语义正确。
// 存数：`"+m"` 读写出 + volatile，语义正确。
#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)
#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 short u16;
typedef unsigned char u8;
typedef unsigned long long u64;
// g++-9 -O2 applies -mavx256-split-unaligned-load, so EVERY LDU
// 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
// [DISP32] compile-time signature-plane row pitch. This problem's n is fixed at
// 300000, so MF_SPL = 300064 = N+64 exactly; the guard below refuses the
// signature path (correct, merely slow) if N ever grew past it.
#define MF_SPL ((size_t)300064)
#endif
#define MAXK 10
#define STRIDE 16
// 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
#define MF_QCOEF 64
#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
#define MF_PFB 24
#endif
// [b16a PFSOUT] pass B 尾部 `out[i] += add` 的写目标行预取（i 取自顺序流 ⇒ 地址无二级依赖）。
// 1 = 开（默认），0 = 关（同源 A/B 对照臂）。纯预取，语义零改动。
#ifndef MF_PFSOUT
#define MF_PFSOUT 1
#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
#define MF_PFA 1
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 4096
#endif
static u8 s_sigrec[(size_t)MAXN * 16];
static u32 s_bandrec[(size_t)MAXN*(MAXK+1)],s_bandacc[MAXN];
static u64 s_bandvec[MAXN];
static u32 s_ord[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u16 s_bd[(size_t)MAXN * MAXK];
static u32 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_frord[MAXK][MAXN];
static u16 s_bkt[MAXN];
static u32 s_bks[MAXN]; // K3a: 排序后顺序存放的桶键（与 s_qord 同序）
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];
static u32 *s_avpool = 0;
static size_t s_avcap = 0;
// ---- 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; // [z9] full-range scale, BRIEF 2.18.821
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;
u32 w = 0, cnt = 0, total = 0;
for (; w + 8 <= Wr; w += 8) {
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
__m256i la = _mm256_and_si256(a, lm);
__m256i ha = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, la), _mm256_shuffle_epi8(lk, ha)));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < KK; d++) b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
__m256i lb = _mm256_and_si256(b, lm);
__m256i hb = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lb), _mm256_shuffle_epi8(lk, hb)));
if (++cnt == 31) {
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
acc0 = acc1 = z; cnt = 0;
}
}
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);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
total = hsum256(tot);
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;
const size_t RS=KK+1;
if (!g_skip_and) {
const u32 nq = q1 - q0;
u32 *rec = s_bandrec;
{
u32 *cv = s_cntv;
for (u32 v = 0; v <= (u32)N; v++) cv[v] = 0;
for (u32 t = 0; t < nq; t++) cv[s_gt[ds][s_qord[q0 + t]] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)N; v++) { u32 c = cv[v]; cv[v] = ac; ac += c; } }
for (u32 t = 0; t < nq; t++) {
u32 i = s_qord[q0 + t];
const u16 *bd = s_bd + (size_t)i * MAXK;
u32 R = s_gt[ds][i];
u32 *rp = rec + (size_t)(cv[R + 1]++) * RS;
rp[0] = R;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
rp[1 + e] = g_plen_on ? s_pofs[(size_t)e * MAXQ + bd[d]]
: (u32)(((size_t)e * (Q + 1) + bd[d]) * W);
}
rp[KK] = i;
}
}
{
u32 *acc = s_bandacc;
for (u32 t = 0; t < nq; t++){ acc[t] = 0; s_bandvec[t]=0; }
if (nq) {
const u32 maxR = rec[(size_t)(nq - 1) * RS];
const u32 maxband = (maxR + 2047u) / 2048u;
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 zz = _mm256_setzero_si256();
u32 ps = 0;
for (u32 bk = 0; bk <= maxband; bk++) {
const u32 lim = bk * 2048u;
while (ps < nq && rec[(size_t)ps * RS] <= lim) ps++;
if (ps >= nq) break;
u32 psp = ps;
while (psp < nq && rec[(size_t)psp * RS] < lim + 2048u) psp++;
const size_t off8 = (size_t)bk * 32u;
for (u32 t = psp; t < nq; t++) {
if(t+4<nq){const u32* next=rec+(size_t)(t+4)*RS;
for(int e=0;e<K1;e++){const u64* pp=s_bs+(size_t)next[1+e]+off8;_mm_prefetch((const char*)pp,_MM_HINT_T0);_mm_prefetch((const char*)(pp+8),_MM_HINT_T0);}
}
const u32 *rp = rec + (size_t)t * RS;
__m256i bytes=zz;
for (u32 lane=0;lane<32u;lane+=4) {
const u64 *p0=s_bs+(size_t)rp[1]+off8+lane;
__m256i a0=LDU256(p0);
for(int e=1;e<K1;e++)a0=_mm256_and_si256(a0,LDU256(s_bs+(size_t)rp[1+e]+off8+lane));
__m256i low=_mm256_and_si256(a0,lm),high=_mm256_and_si256(_mm256_srli_epi16(a0,4),lm);
bytes=_mm256_add_epi8(bytes,_mm256_add_epi8(_mm256_shuffle_epi8(lk,low),_mm256_shuffle_epi8(lk,high)));
}
__m256i sum64=_mm256_sad_epu8(bytes,zz);
__m128i lanes=_mm_add_epi32(_mm256_castsi256_si128(sum64),_mm256_extracti128_si256(sum64,1));
lanes=_mm_shuffle_epi32(lanes,0x88);
// Both packed 32-bit totals are bounded by N, so a 64-bit add cannot carry between lanes.
s_bandvec[t] += (u64)_mm_cvtsi128_si64(lanes);
}
for (u32 t = ps; t < psp; t++) {
const u32 *rp = rec + (size_t)t * RS;
const u32 R = rp[0];
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) bp[e] = s_bs + (size_t)rp[1 + e] + off8;
acc[t] += and_popcount_range<KK - 1>(bp, (R >> 6) - (u32)off8, R & 63u);
}
}
for (u32 t = 0; t < nq; t++){__m128i lanes=_mm_loadl_epi64((const __m128i*)(s_bandvec+t));out[rec[(size_t)t * RS + KK]] = acc[t]+(u32)_mm_cvtsi128_si32(lanes)+(u32)_mm_extract_epi32(lanes,1);}
}
}
}
if (g_skip_fringe) return;
u32 *const avpool = s_avpool;
// ---- 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;
const u32 *fo = s_frord[d];
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);
// [f5/setup-hoist] the k byte-plane base pointers depend only on d, so build them
// once per dimension instead of once per (query,dimension).
const u8 *const sbase = s_sig + (size_t)d * k * MF_SPL;
const u32 *const odd = s_ord[d];
for (u32 qi = q0; qi < q1; qi++) {
u32 i = fo[qi] & 0xFFFFFu;
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB;
if (qj < q1) {
u32 j = fo[qj] & 0xFFFFFu;
_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 *)(avpool + (size_t)j * STRIDE), _MM_HINT_T0);
#if MF_PFSOUT
_mm_prefetch((const char *)(out + j), _MM_HINT_T0); // [b16a PFSOUT] 写目标行（尾 store）
#endif
}
}
#endif
const u32 lo = FPOS(d, (fo[qi] >> 20)); // K3a: 原来是 s_bd[i*MAXK+d]（同元素第二次随机冷读）
const u32 hi = s_gt[d][i];
if (lo >= hi) continue;
const u32 *rec = s_pm + (size_t)i * STRIDE;
const u32 *Avp = avpool + (size_t)i * STRIDE;
__m256i qr = _mm256_or_si256(LDA((const __m256i *)rec), khi8v);
__m256i Avm = _mm256_or_si256(LDU256(Avp), notexlo);
__m256i qT = _mm256_min_epi32(qr, Avm);
__m256i qTh;
if (k > 8) {
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU256((Avp + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
} else {
qTh = all8;
}
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 (g_sig_on && (hi - lo) >= 2) {
// ---- 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 capv = _mm256_set1_epi32(255);
// [f5/setup-hoist] sv depends on the half (e<8) only, not on e itself.
// [z9] 2.18.821 full-range map: high 32 bits of v*MAG, even lanes via mul_epu32
// and odd lanes via a 32-bit shift, blended, then clamped at 255 (the sentinel
// lanes used to land on 0xFF through the old shift form -- keep that exactly).
const __m256i sv_lo = _mm256_min_epu32(_mm256_blend_epi32(
_mm256_srli_epi64(_mm256_mul_epu32(qT, magv), 32),
_mm256_slli_epi64(_mm256_srli_epi64(_mm256_mul_epu32(_mm256_srli_epi64(qT, 32), magv), 32), 32), 0xAA), capv);
const __m256i sv_hi = _mm256_min_epu32(_mm256_blend_epi32(
_mm256_srli_epi64(_mm256_mul_epu32(qTh, magv), 32),
_mm256_slli_epi64(_mm256_srli_epi64(_mm256_mul_epu32(_mm256_srli_epi64(qTh, 32), magv), 32), 32), 0xAA), capv);
for (int e = 0; e < k; e++) {
if (e == (int)d) { qsb[e] = _mm256_set1_epi8((char)0xFF); continue; }
__m256i src = (e < 8) ? sv_lo : sv_hi;
__m256i bc = _mm256_permutevar8x32_epi32(src, _mm256_set1_epi32(e & 7));
qsb[e] = _mm256_shuffle_epi8(bc, zero8);
}
}
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
__m256i ac = _mm256_setzero_si256();
{ /* SGNL: the e==d lane is a PROVABLE NO-OP -- qsb[d]=0xFF,
so max_epu8(cv,0xFF)=0xFF and cmpeq(0xFF,0xFF)=-1 (lt unchanged),
and cv = coord>>sh <= (N-1)>>sh <= 255 < 0xFF so the eq bit never
fires. Dropping it shortens the serial lt/eq chain by one link. */
for (int e = 0; e < k; e++) {
if (e == (int)d) continue;
__m256i cv = LDU256((sbase + (size_t)e * MF_SPL + aj));
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[e]));
}
}
__m256i lt = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = (u32)_mm256_movemask_epi8(lt);
u32 amb = 0;
if (ltm) { // [a15y] lazy equality mask: rare, planes are L1-hot
__m256i eq = _mm256_setzero_si256();
for (int e = 0; e < k; e++) {
if (e == (int)d) continue;
__m256i cv = LDU256((sbase + (size_t)e * MF_SPL + aj));
eq = _mm256_or_si256(eq, _mm256_cmpeq_epi8(cv, qsb[e]));
}
amb = ltm & (u32)_mm256_movemask_epi8(eq);
}
aacc += (u32)__builtin_popcount(ltm & ~amb);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 ? (c0 + (size_t)(aj + (u32)t - lo) * KS_C)
: (s_pm + (size_t)odd[aj + (u32)t] * STRIDE);
__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 ac = _mm256_setzero_si256();
{ /* SGNL tail: same provable no-op */
for (int e = 0; e < k; e++) {
if (e == (int)d) continue;
__m256i cv = LDU256((sbase + (size_t)e * MF_SPL + base));
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[e]));
}
}
__m256i lt = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = ((u32)_mm256_movemask_epi8(lt)) & kmask;
u32 amb = 0;
if (ltm) { // [a15y] lazy equality mask (tail chunk)
__m256i eq = _mm256_setzero_si256();
for (int e = 0; e < k; e++) {
if (e == (int)d) continue;
__m256i cv = LDU256((sbase + (size_t)e * MF_SPL + base));
eq = _mm256_or_si256(eq, _mm256_cmpeq_epi8(cv, qsb[e]));
}
amb = ltm & ((u32)_mm256_movemask_epi8(eq));
}
aacc += (u32)__builtin_popcount(ltm & ~amb);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 ? (c0 + (size_t)(base + (u32)t - lo) * KS_C)
: (s_pm + (size_t)odd[base + (u32)t] * STRIDE);
__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 = LDA((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDA((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
}
}
template<int K>
__attribute__((target("avx2")))
static void solve_mf(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)MF_SPL;
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;
const u64 bmag = ((1ull << 40) + B - 1) / 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
{
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];
// [v4_orders] gt[i] = start[x_d[i]] and bd[] from start[x_d[i]+1]: value-indexed
// lookups, so the per-point work is two sequential sweeps of xd[] with no random
// od[]/xd[] traffic (the old loop walked od[] reading xd[od[j]] at random and then
// scattered gt[]/s_bd[] writes once per tie group).
// [b16y_ orders-merge] pass3 (counting-sort scatter) and pass4 (gt/bd) fused,
// with {start, running-position} interleaved so one line serves both random
// accesses. Bit-identical: same ops, same order, s_cntv never mutated in-loop.
u32 *PR = s_posp; // PR[2v]=start (was s_cntv[v]), PR[2v+1]=run (was s_tmp[v])
for (u32 i = 0; i <= N; i++) PR[2u * i] = 0;
for (u32 i = 0; i < N; i++) {
if (i + 24 < N) { u32 vf = xd[i + 24]; if (vf >= N) vf = 0;
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(PR + 2u * vf))); }
u32 v = xd[i]; PR[2u * (v < N ? v : 0)]++;
}
{ u32 s = 0; for (u32 i = 0; i <= N; i++) { u32 c = PR[2u * i]; PR[2u * i] = s; PR[2u * i + 1] = s; s += c; } }
{
u32 *gt = s_gt[d];
u32 *pp = s_pos[d];
u16 *col = (u16 *)(s_posp + 2u * (size_t)MAXN) + (size_t)d * ((size_t)N + 64);
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
u32 c0 = PR[2u * v];
u32 q = PR[2u * v + 1]++;
u32 num = PR[2u * (v + 1)] - 1u; if (num > N) num = N;
od[q] = i;
gt[i] = c0;
u32 bb = (u32)(((u64)num * bmag) >> 40);
col[i] = (u16)(bb > Q - 1 ? Q - 1 : bb);
if (g_use_filt) pp[i] = q;
} }
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;
}
// [b16y] s_bd: one sequential row-major pass from the per-dim columns written above
// (the old form stored u16 at i*MAXK+d, i.e. 3e5 stores over a 6 MB span per dim and
// 10 passes over it => ~60 MB of write-allocate traffic for 6 MB of data).
{
const u16 *bcol = (const u16 *)(s_posp + 2u * (size_t)MAXN);
const size_t cstride = (size_t)N + 64;
for (u32 i = 0; i < N; i++) {
u16 *row = s_bd + (size_t)i * MAXK;
for (int dd = 0; dd < k; dd++) row[dd] = bcol[(size_t)dd * cstride + 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];
}
// ---------- 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; // [z9] exact monotone full-range map
g_sigmag = (u32)((((u64)255 << 32) / (u64)dn) + 1u); }
size_t need = (size_t)k * k * MF_SPL;
static u8 *sigpool = 0; static size_t sigcap = 0;
if ((size_t)N + 32 > MF_SPL) { sigpool = 0; sigcap = 0; } // fringe max offset hi+31 <= N+30
// 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), LDA((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), LDA((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 needav = (size_t)N * STRIDE * 4 + 256;
if (needav > s_avcap) {
s_avpool = (u32 *)pool_alloc(needav);
s_avcap = s_avpool ? needav : 0;
}
}
{
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;
// 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).
const u16 *bdv = s_bd + (size_t)i * MAXK;
u32 *av = s_avpool + (size_t)i * STRIDE;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
for (int d = 0; d < k; d++) {
if (d == ds) { av[d] = 0x7FFFFFFFu; continue; }
u32 b = bdv[d];
av[d] = BVAL(d, b);
if (g_plen_on && best > mr[(size_t)d * MAXQ + b])
mr[(size_t)d * MAXQ + b] = best;
}
}
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; }
// Sort each dimension's query IDs by its existing coordinate order once.
// The same ID is written once per dimension into its selected fold group.
const u16 *bcol_fr = (const u16 *)(s_posp + 2u * (size_t)MAXN);
for (int d = 0; d < k; d++) {
u32 pos[MAXK];
for (int ds = 0; ds < k; ds++) pos[ds] = gstart[ds];
for (u32 r = 0; r < N; r++) {
u32 id = s_ord[d][r];
u32 dst = pos[s_g[id]]++;
s_frord[d][dst] = id | ((u32)bcol_fr[(size_t)d * ((size_t)N + 64) + id] << 20);
}
}
// ---------- 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 (s_sig) {
/* [x16z5] (265) redundancy lens: sig(e,p) = f(x[e][p]) is INDEPENDENT of d, but the
plane build below recomputes it k times per point while walking s_pm rows
(k*N random 64B reads = 192 MB). Compute it ONCE per (e,p) here, sequentially
from x[][], then let the plane build scatter it out of this small table. */
u8 *sr = s_sigrec;
for (u32 p = 0; p < N; p++)
for (int e = 0; e < k; e++)
sr[(size_t)p * 16 + e] = (u8)((u32)(((u64)x[e][p] * g_sigmag) >> 32));
}
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] = 0; // [b16a] transposed rows no longer materialised (97.7% never read)
const u32 *od = s_ord[d];
if (s_sig) {
u8 *sb = s_sig + (size_t)d * k * MF_SPL;
u32 j = 0;
if (k == 10) {
for (; j + 8 <= N; j += 8) {
__m128i r0 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+0] * 16));
__m128i r1 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+1] * 16));
__m128i r2 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+2] * 16));
__m128i r3 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+3] * 16));
__m128i r4 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+4] * 16));
__m128i r5 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+5] * 16));
__m128i r6 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+6] * 16));
__m128i r7 = _mm_loadu_si128((const __m128i *)(s_sigrec + (size_t)od[j+7] * 16));
__m128i t0 = _mm_unpacklo_epi8(r0,r1);
__m128i t1 = _mm_unpacklo_epi8(r2,r3);
__m128i t2 = _mm_unpacklo_epi8(r4,r5);
__m128i t3 = _mm_unpacklo_epi8(r6,r7);
__m128i u0 = _mm_unpacklo_epi16(t0,t1);
__m128i u1 = _mm_unpackhi_epi16(t0,t1);
__m128i u2 = _mm_unpacklo_epi16(t2,t3);
__m128i u3 = _mm_unpackhi_epi16(t2,t3);
__m128i v0 = _mm_unpacklo_epi32(u0,u2);
__m128i v1 = _mm_unpackhi_epi32(u0,u2);
__m128i v2 = _mm_unpacklo_epi32(u1,u3);
__m128i v3 = _mm_unpackhi_epi32(u1,u3);
__m128i h0 = _mm_unpackhi_epi8(r0,r1);
__m128i h1 = _mm_unpackhi_epi8(r2,r3);
__m128i h2 = _mm_unpackhi_epi8(r4,r5);
__m128i h3 = _mm_unpackhi_epi8(r6,r7);
__m128i hu0 = _mm_unpacklo_epi16(h0,h1);
__m128i hu1 = _mm_unpacklo_epi16(h2,h3);
__m128i v4 = _mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(sb + (size_t)0 * MF_SPL + j),v0);
_mm_storel_epi64((__m128i *)(sb + (size_t)1 * MF_SPL + j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(sb + (size_t)2 * MF_SPL + j),v1);
_mm_storel_epi64((__m128i *)(sb + (size_t)3 * MF_SPL + j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(sb + (size_t)4 * MF_SPL + j),v2);
_mm_storel_epi64((__m128i *)(sb + (size_t)5 * MF_SPL + j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(sb + (size_t)6 * MF_SPL + j),v3);
_mm_storel_epi64((__m128i *)(sb + (size_t)7 * MF_SPL + j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(sb + (size_t)8 * MF_SPL + j),v4);
_mm_storel_epi64((__m128i *)(sb + (size_t)9 * MF_SPL + j),_mm_srli_si128(v4,8));
}
}
for (; j < N; j++) {
const u8 *srr = s_sigrec + (size_t)od[j] * 16;
for (int e = 0; e < k; e++) sb[(size_t)e * MF_SPL + j] = srr[e];
}
}
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 * MF_SPL;
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++) base[(size_t)e * MF_SPL + j] = (u8)((u32)(((u64)row[e] * g_sigmag) >> 32));
}
}
}
// ---------- 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 = s_gt[ds][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++) {
if (j + 24 < N) __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(s_nid + od[j + 24])));
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 + 8 <= pend) {
if (ptr + 40 < pend) __asm__ __volatile__("prefetcht0 %0" :: "m"(*(const char *)(s_nid + od[ptr + 40])));
u32 q0 = s_nid[od[ptr+0]], q1 = s_nid[od[ptr+1]], q2 = s_nid[od[ptr+2]], q3 = s_nid[od[ptr+3]];
u32 q4 = s_nid[od[ptr+4]], q5 = s_nid[od[ptr+5]], q6 = s_nid[od[ptr+6]], q7 = s_nid[od[ptr+7]];
cur[q0>>6] |= 1ull<<(q0&63); cur[q1>>6] |= 1ull<<(q1&63); cur[q2>>6] |= 1ull<<(q2&63); cur[q3>>6] |= 1ull<<(q3&63);
cur[q4>>6] |= 1ull<<(q4&63); cur[q5>>6] |= 1ull<<(q5&63); cur[q6>>6] |= 1ull<<(q6&63); cur[q7>>6] |= 1ull<<(q7&63);
ptr += 8;
}
while (ptr < pend) {
if (ptr + 12 < pend) __asm__ __volatile__("prefetcht0 %0" :: "m"(*(const char *)(s_nid + od[ptr + 12])));
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), LDA((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<4>((u32)n, (const unsigned **)x, out); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf<5>((u32)n, (const unsigned **)x, out); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf<6>((u32)n, (const unsigned **)x, out); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf<7>((u32)n, (const unsigned **)x, out); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf<8>((u32)n, (const unsigned **)x, out); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf<9>((u32)n, (const unsigned **)x, out); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf<10>((u32)n, (const unsigned **)x, out); }
#ifdef LOCAL_TEST
static void brute_nd(int n, const u32 **x, u32 *out, int k) {
for (int i = 0; i < n; i++) {
u32 c = 0;
for (int j = 0; j < n; j++) {
int ok = 1;
for (int d = 0; d < k; d++) if (!(x[d][j] < x[d][i])) { ok = 0; break; }
c += ok;
}
out[i] = c;
}
}
int main(int argc, char **argv) {
int k = (argc > 1) ? atoi(argv[1]) : 6;
int n = (argc > 2) ? atoi(argv[2]) : 200;
int seed = (argc > 3) ? atoi(argv[3]) : 1;
int mode = (argc > 4) ? atoi(argv[4]) : 0;
int qq = (argc > 5 && atoi(argv[5]) > 0) ? atoi(argv[5]) : 0;
g_Q = qq ? (u32)qq : (u32)(n / 250 + 1);
g_B = (u32)((n + g_Q - 1) / g_Q);
if (argc > 6) g_nosort = atoi(argv[6]);
if (argc > 7) g_nosdmp = atoi(argv[7]);
if (argc > 8) g_skip_and = atoi(argv[8]);
if (argc > 9) g_skip_fringe = atoi(argv[9]);
if (argc > 10) g_use_filt = atoi(argv[10]);
if (argc > 11) g_lexsort = atoi(argv[11]);
const int dm = (mode >= 10) ? mode - 10 : mode;
srand(seed);
static u32 *xs[10];
u32 *out = (u32 *)malloc(sizeof(u32) * n);
u32 *ref = (u32 *)malloc(sizeof(u32) * n);
for (int d = 0; d < k; d++) {
xs[d] = (u32 *)malloc(sizeof(u32) * n);
for (int i = 0; i < n; i++) {
if (dm == 0) xs[d][i] = (u32)(rand() % n);
else if (dm == 1) xs[d][i] = (u32)(rand() % 3);
else if (dm == 2) xs[d][i] = 0;
else if (dm == 3) xs[d][i] = (u32)(n - 1 - (rand() % n));
else if (dm == 4) xs[d][i] = (u32)(i);
else xs[d][i] = (u32)(rand() % n);
}
}
if (mode >= 10) {
struct timespec t0, t1;
clock_gettime(CLOCK_MONOTONIC, &t0);
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
clock_gettime(CLOCK_MONOTONIC, &t1);
double ms = (t1.tv_sec - t0.tv_sec) * 1e3 + (t1.tv_nsec - t0.tv_nsec) / 1e6;
u64 s = 0; for (int i = 0; i < n; i++) s += out[i];
printf("k=%d n=%d Q=%u TIME %.1f ms checksum %llu\n", k, n, g_Q, ms, s);
return 0;
}
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
brute_nd(n, (const u32 **)xs, ref, k);
int bad = 0;
for (int i = 0; i < n; i++) if (out[i] != ref[i]) { if (bad < 5) printf("MISMATCH i=%d got=%u ref=%u\n", i, out[i], ref[i]); bad++; }
printf("k=%d n=%d mode=%d Q=%u %s (%d mismatches)\n", k, n, mode, g_Q, bad ? "FAIL" : "OK", bad);
return bad != 0;
}
#endif
// [a16z] 最小驱动：n / mod 参数，clock_gettime 计时 + checksum
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cstdint>
#include <ctime>
void count_10d(int n, const unsigned *x[10], unsigned *out);
static unsigned gs = 12345u;
static inline unsigned nxr(void){ gs^=gs<<13; gs^=gs>>17; gs^=gs<<5; return gs; }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 682.089 ms | 186 MB + 748 KB | Accepted | Score: 100 | 显示更多 |