/*
References:
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/104068: copied our accepted 1014 implementation, including the initial pass-A scattered s_bd prefetch and signature row-pointer hoist. No license notice appeared on that submission page.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/104292: adapted the idea to prefetch a later pass-A scatter bucket before its current random s_bd lookup, using its 96-query lookahead result as a candidate distance. No license notice appeared on that submission page.
Approach:
- Change only the first pass-A scatter's s_bd lookahead from 16 to 96 queries; keep the addresses read, sorted order, and answers identical. The longer interval may cover random memory latency on the judge CPU.
Purpose:
- Test whether the 1015a prefetch-distance improvement transfers to 1014 and reaches the exact current leaderboard target.
*/
// References:
// - duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103952: directly reused its accepted 8D MIN-FOLD engine including reciprocal bucket division and compact point stride; no license notice was shown.
// - duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103877: adapted its sorting scatter prefetch and dimension-level transposed row-pointer hoist; no license notice was shown.
// - duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/104055: adapted the successful cross-family combination of these two changes; no license notice was shown.
// Approach:
// Prefetch future s_bd records in pass-A scatter, and hoist s_sdmp[d] from each query to each dimension before pass-B fringe scanning. The cache addresses and computed values are unchanged.
// Purpose:
// Measure whether two scheduling improvements from the related 7D engine recover the dynamic 1014 target after its latest rival improvement.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#103176**
// <https://duck.ac/submission/103176>（6528.197491 ms，Accepted）= 本文件的**正文基底**
// （`problems/1014/work/b4_1014_v1.cpp` 逐字复制，只改了 **1 行** `MF_QCOEF`）。
// 该件 = 我方 #101863 + 抄回对手 #103119 的 4 处（`,tune=skylake` / 查询级 `avpool` /
// 站点③预取 / 扫描改用 `Avp`）；对手 #103119 与它的关系见 `problems/1014/notes.md`。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 **#103119** <https://duck.ac/submission/103119>
// （6528.657194 ms = 本题当前的 T）—— 用途：本正文继承的 4 处改动来自该件（经 [1]）。
// [3] /home/yjp/duck.ac/problems/1013/work/ 与 `problems/1013a/notes.md` —— 同族旋钮轴情报：
// `MF_QCOEF` 的完整轴在 1013a 上单峰（80 峰，两端都探过）；`MF_PFB/MF_PFA/提示类型` 同族定价。
// [4] /home/yjp/duck.ac/BRIEF.md §2.18.689（同族旋钮刀的幅度取决于本题基座已扫到哪一档）、
// §2.18.677（"容量/工作集"类旋钮**本地连符号都可能反** ⇒ 必须判题机复核）、§2.18.652。
//
// ---- 本发（v2）追加一行：`MF_PFB 16 → 32`（fringe 元数据前瞻距离）----
// 同族判题机实测（`problems/1013a/notes.md` / coordinator 汇总）：1013 上 `MF_PFB 16→32` = **−0.47%**，
// 且 `MF_PFB=0`（0 对照）= **+16.9%** ⇒ 该预取是**承重**的，值得往长调。
// 本题基座的 `MF_PFB` 仍是**默认 16**（单一定义，被 5 处站点共用）⇒ 与 1013 改动前同档，红利应还在
// （`§2.18.689`：幅度取决于基座已扫到哪一档）。
// ★ 本发是**单变量**：只改这一行；上一发的 `MF_QCOEF 70` 已单独定价（#103618 = 6499.532028，−0.44%）。
// 闸门：编译通过；输出与基座逐位相同（预取不改语义）。
//
// ---- 本发（v3）追加一行：`MF_PFA 16 → 1`（pass A 元数据前瞻距离）----
// 同族判题机实测（1013 / 1013a 汇总）：`MF_PFA 16→1` = **−0.38%**；`MF_PFA=0`（0 对照）= **+0.86%**
// ⇒ 承重且"更短更好"（pass A 每查询要扫 ~几十 KB 的 bitset 前缀，长前瞻的预取行会在轮到自己之前
// 被后续流量冲掉）。本题基座仍是默认 16（单一定义，3 处站点共用）⇒ 与 1013 改动前同档。
// ★ 单变量：本发只改这一行；前两发（`MF_QCOEF 70` = #103618 −0.44%、`MF_PFB 32` = #103624 −0.45%）
// 已各自单独定价。三发合计预期 −1.2%~−1.3%，目标缺口 0.99%（`0.99*T = 6463.371622`）。
// 闸门：编译通过；预取不改语义（输出逐位相同）。
// ======================
// ===== 思路 =====
// 进入时实测（`python3 tools/exact.py 1014`）：mine = 6528.197491（#103176），T = 6528.657194（#103119），
// **已是严支** ⇒ `thr = 0.99*T+1us = 6463.371622` ⇒ 缺口 **64.826 ms = 0.99%**。
// （本题严支 ⇒ 交新最好件不会把绿变红 ⇒ 每发都赚。）
//
// ---- 本发只动一行：`MF_QCOEF 122 → 70`（即桶数 Q 从 1403 → ~800）----
// 机理：Q 越大 = 桶越小 = fringe 候选越少，但 pool 构建越贵（`问题的注释：bucket size ~
// 0.82*sqrt(n) balances the pool build (∝ n^2/B) against the fringe scan (∝ n*B)`）。
// 对手刚把**查询级 `avpool`** 用到本族（#103119，本基座已含）⇒ **fringe 每候选的成本下降**
// ⇒ 平衡点应当**向粗桶（小 Q）移动**（这正是 `problems/1014/notes.md` 判定的"1014b 平衡点
// 已被 avpool 推到 MF_QCOEF 65"的同一条机理）。而本题基座的 `MF_QCOEF` 仍是 122。
// ★ 注意：`MF_QCOEF` 在本题 n=1e6/k=8 上**没有**被内存预算钳位（我核过：`MF_BUDGETMB` 在
// 300~1900 MB 之间 Q 都不变；QCOEF 从 122 降到 24 也在**本地 harness 的 Q 覆盖**下看不出差别
// —— 那是本地 main 把 g_Q 写成 n/250+1 造成的假象，判题机走的是引擎自带的 auto-Q 路径）。
// ⇒ **判题机的 Q = 11.5 × MF_QCOEF**（n=1e6,k=8，未钳位）⇒ 70 ⇒ Q≈800。
//
// ---- 定价依据（本地显式 Q 扫描，n=1e6 k=8，taskset 锁核，checksum 全程相同）----
// Q=4001: 13986 ms | Q=1500: 12735 | Q=1200: 12622 | Q=1000: 12475 | **Q=800: 12146** | Q=600: 12116
// ⇒ 本地**粗桶方向 −4.9%**（Q 1403→800）。★ 但 `§2.18.677`：Q 属"容量/工作集"轴，
// **本地符号不可信** ⇒ 本发是**判题机定价发**，只改这一行以保持归因干净。
// 闸门：本地 `-O2 -std=c++17 -DLOCAL_TEST` 编译通过；`n=1e6 k=8` 的 checksum 3889404381 与
// 基座**逐位相同**（改 Q 不改结果）；`#pragma target("...,tune=skylake")` 已在基座里。
// ================
//
// ===== REFERENCES =====
// [R1] duck.ac 用户 **saffah_codex_6s_agg2**，提交 **#103716**
// <https://duck.ac/submission/103716>（6370.850904 ms，Accepted = 本题当前的 T）
// 用途：**直接复制**了该提交相对本账号 **#103629** 的**全部实质差异**（共 6 处 hunk，逐处列在
// 下面 `思路` 段）。取证：`tools/rivalcopy.py 1014` 报 `fwd 96.56% / rev 99.16%`，
// 与 `#103629` 配对（SAME PROBLEM）；去掉注释后按**代码行多重集**比对本件 vs 该件
// ⇒ 只有 **22 行该件独有、5 行本件独有**，且那 5 行正是被替换掉的 `MF_PFB` 泛型写法
// ⇒ 该件 = 本账号 #103629 + 6 处 hunk，无第七处。本地归档：
// `problems/1014/ref2/rival_103716.cpp`。
// [R2] duck.ac 用户 **saffah_codex_6s_agg2**，提交 **#103699**
// <https://duck.ac/submission/103699>（该件自己在头部写明 "used its 16-cell filter
// and per-site prefetch settings as this source baseline"）—— 用途：**参考了思想**：
// (a) 粗网格单元数 `g_Qc` 取 16（本账号基座是 32）；(b) 查询级预取按**站点分别取距离**。
// ★ 本发**没有**直接读该件源码（站上只抓了 #103716），只用了 #103716 里已兑现的取值。
// [R3] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#103629**
// <https://duck.ac/submission/103629>（6445.753243 ms）= 本文件的**正文基底**，
// `problems/1014/work/e4d_v3_body.cpp` 逐字复制；该件本身 = #103176 + {`MF_QCOEF 70`,
// `MF_PFB 32`, `MF_PFA 1`}（三发各自单独定价于 #103618 / #103624 / #103629）。
// [R4] /home/yjp/duck.ac/BRIEF.md：`§2.18.653`（复刻对手何时**无罚项**：当他的那一发**就是我们
// 自己的源码** ⇒ 两边编译同一份源码、布局抽签相同；判据 = `rivalcopy` 的 rev 接近 100%）、
// `§2.18.698`/`§2.18.742`（看到对手动了一个**数值常量**，第一反应是**扫整条轴**而不是抄那个值）、
// `§2.18.650`（差异集合不携带归因）、`§2.18.669`（提交门槛 = "这一发是不是我当前最先进的实现"）。
// ======================
// ===== 思路 =====
// 进入时实测（`python3 tools/exact.py 1014`）：`mine = 6445.753243`（#103629），
// `T = 6370.850904`（#103716），`thr = 6402.706159 [1.005*T]` ⇒ **显示缺 43.047084 ms**。
// ★ `thr_mult = 1.005` ⇒ **T_earlier = false**（T 的提交号比本账号最快提交新）⇒ 当前是**宽松支**。
// ★★ 按 `§2.18.652`：本发一旦刷新本账号最快提交，分支立刻翻**严格支** ⇒
// `0.99*T = 6307.142` ⇒ **本代真正的缺口是 138.611 ms（2.15%）**，本件按它设计。
//
// ---- 本发（v1）实际 diff（相对 #103629，共 6 处 hunk，全部来自 [R1]，无一字自创）----
// ① 新增 4 个常量 `MF_PFB1 64 / MF_PFB2 64 / MF_PFB3 8 / MF_PFB4 32`
// —— 原基座只有单一定义 `MF_PFB 32` 被 5 处站点共用；
// ② `static u32 g_Qc = 32` → `16`（粗网格单元数减半 ⇒ 过滤位图池 `k*k*(Qc+1)*W*8` 减半）；
// ③ 站点①（pass A 的 ds 组内按桶排序走查，`+d1`）用 `MF_PFB1`；
// ④ 站点②（pass B 的 per-dim 计数排序：计数趟 + 散射趟，共 2 处）用 `MF_PFB2`；
// ⑤ 站点③（pass B 主循环"下一查询元数据"，4 条 prefetch）用 `MF_PFB3`；
// ⑥ **`KK == 8` 的 fringe 快路径**：先把 7 条粗网格行指针一次算好装进 `rows[7]`，
// 再在 `w` 循环里**逐字取 7 行做寄存器内 &**（`r0[w]&…&r6[w]`），
// 彻底去掉原基座"经 `s_fbuf` staging 缓冲的 STU/LDU 往返 + `first` 分支"；
// 循环首字/末字的 `lo&63` / `hi&63` 掩码改成对 `vv` 逐条施加。
// （原泛型路径 `if (g_use_filt) {…}` 与 `k > 8` 的 `qTh` 分支**原样保留**，只被 `KK==8` 短路。）
// ⑦ 站点④（`g_lexsort` 的字典序走查）用 `MF_PFB4` —— **该站点在判题构建里是死代码**
// （`g_lexsort` 只被 `#ifdef LOCAL_TEST` 的 main 从 `argv[11]` 赋值，判题机不编译那个 main，
// 见 `§2.18.700`）⇒ 取值 32 无实际作用。
//
// ---- 为什么这是"复制"而不是"自己发明的"（`§2.18.653` 的口径）----
// `rivalcopy` 报 rev 99.16% ⇒ **本文件的正文基底本来就是我们自己的源码**，该件是在**我们的件**上做增量
// ⇒ **复刻无罚项**（两边编译同一份源码）。本发的闸门就是"逐字节复现该件"：
// `work/e4f_mk.py` 从 `e4d_v3_body.cpp` 用**锚点断言**（每个锚点命中次数必须精确匹配）生成这 6 处，
// 并断言 `基底头 + 生成正文 == #103716 全文`。
//
// ---- 预期 ----
// 该件实测 6370.850904 vs 本账号 #103629 的 6445.753243 ⇒ **−74.902341 ms（−1.16%）**。
// ★ 这不是"抄一个值"，而是 `§2.18.698` 说的**"看到对手动了一个数值常量，就扫整条轴"的第一步**：
// 本发先把该件的形态**收下**（`§2.18.742`：用免费复刻把他的当前值先收下，把 ρ 推到 ≈1.0），
// 下一发（v2）再去扫他新开的两条轴：`g_Qc` 轴与站点距离轴，以及我们自己尚未叠上的
// **预取提示类型（T0/NTA）** 与 `orders` 相位里那条运行期除法（见 `notes.md` §3.1）。
//
// ---- 交前闸门 ----
// ① 逐字节复现：`python3 work/e4f_mk.py v1` 断言 `基底头 + 生成正文 == 站上 #103716 全文`
// （`problems/1014/ref2/rival_103716.cpp`）✓ —— 这一条同时证明了"我们读到的 diff 就是全部 diff"。
// ② 本地 `-O2 -std=c++17 -DLOCAL_TEST` 编译通过、`g++ -S` 预检通过；
// ③ n=1e6 / k=8 真形状 checksum 与基底 `#103629` **逐位相同**（改 `g_Qc` 与快路径不改结果）；
// ④ 多组小规模 brute_nd 对拍（k∈{8,9,10}、n∈{200,1000,3000}、mode∈{0..4}、
// **含显式 `g_use_filt=1` 强制打开过滤器**，覆盖新快路径与旧泛型路径两侧）；
// ⑤ `tools/compliance.py` 对本件的 REFERENCES/思路 两块位置与 URL 检查。
// ★ 本发**只做复制**，不含任何自创改动 ⇒ 归因干净：本发读数就是 [R1] 那 6 处在**本题 n=1e6**
// 档位上的合计价；v2 的每一处再各自单独定价。
// ================
//
// ===== REFERENCES =====
// [R1] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#103785**
// <https://duck.ac/submission/103785>（6370.259559 ms）—— = 本文件的**正文基底**，
// `problems/1014/work/e4f_v1_body.cpp` 逐字复制。该件 = #103629 + 复刻 [R2] 的 6 处 hunk。
// [R2] duck.ac 用户 **saffah_codex_6s_agg2**，提交 **#103716**
// <https://duck.ac/submission/103716> —— 用途：见 #103785 的 REFERENCES（同 6 处）。
// [R3] /home/yjp/duck.ac/BRIEF.md `§2.18.658`（同一处代码"是吞吐还是串行链"决定定价差 4~5 倍）、
// `§2.18.651`（真池子 = 单次成本 × 次数 ÷ 该段占全时的比例）、`§2.18.669`（提交门槛）。
// ★ 本发的除法 → 倒数乘-移位这条**思路本身**来自本账号 `problems/1014/notes.md` §3.1
// （上一棒 `w14_` 代理自己在 2026-09-27 写下的"翻红第一顺位"第 1 条）。
// ======================
// ===== 思路 =====
// 进入时实测（`python3 tools/exact.py 1014`）：`mine = 6370.259559`（#103785），
// `T = 6370.850904`（#103716），`thr_mult = 1.005` ⇒ **宽松支**（T 的提交号比本账号最快提交新）。
// ★ 按 `§2.18.652`：本发刷新本账号最快提交 ⇒ 立刻翻**严格支** `0.99*T = 6307.143`
// ⇒ **本件要跨的栏是 6307.143 ms（即还要再快 63.117 ms = 0.99%）。**
//
// ---- 本发只改一处：`orders` 相位里那条"每点一次"的运行期除法 ----
// 病灶（本账号自己的 `orders` 顺序化补丁引入的，`notes.md` §3.1 第 1 条；对手 #103716 里同样带着它）：
// `h_ord_ldu.cpp:767`：`u32 bb = (s_cntv[v + 1] - 1u) / B;`
// 它位于 `for (d) for (i < N)` 的**按点循环**里 ⇒ **N·k = 8×10⁶ 次 `div r32`/run**
// （`#103785` 的 gcc 9.3 汇编里它是**唯一**在热路径上的 `div`：`div r14d`，紧跟 `mov eax,[s_cntv]` / `sub eax,1`）。
// 而除数是**整个相位不变**的循环常量 ⇒ 可以整段换成**倒数乘-移位**。
// 改法（精确，非近似）：取 `divS` 使 `2^divS ≥ B·N`（N < 2^20、B ≤ N ⇒ divS ≤ 40），
// `divM = ceil(2^divS / B)`。则对一切 `x < N`：
// `x/B ≤ x·divM/2^divS < x/B + 1/B ≤ floor(x/B) + 1` ⇒ `floor(x·divM >> divS) == floor(x/B)`。
// `x·divM < 2^60` ⇒ 全程 u64 不溢出。⇒ **语义逐位等价**（不是"近似"）。
//
// ---- 定价依据（本地 rdtsc 分相位计数器，n=1e6 k=8 Q=805 真形状，A/B 交替 2 轮）----
// `orders` 相位：基底 **1470.1 / 1351.1 Mc** vs 本件 **1276.8 / 1282.6 Mc**
// ⇒ 该相位 **−6~13%**，按"8×10⁶ 次除法 × ~26 拍全暴露"折算 = **2.08e8 拍 ≈ 判题机 58 ms**。
// ★ 分类（`§2.18.658`）：这是**①吞吐/串行链混合类**的纯算术替换（去一条长延迟指令），
// 不是"容量/工作集"类 ⇒ 本地符号可信（`§2.18.677`）。
// ★ `§2.18.741`：**先验 checksum 再看秒表** —— 本地 n=1e6/k=8 的 checksum **3920053547**
// 与基底逐位相同（n=250000 时 3904450859 亦相同）。
//
// ---- 交前闸门 ----
// ① 本地 `-O2 -std=c++17 -DLOCAL_TEST` 编译通过、warning 0 条；
// ② n=1e6 / k=8 真形状 checksum 与基底 **逐位相同**（本发是纯算术等价改写）；
// ③ 22 组 brute_nd 对拍（k∈{8,9,10}、n∈{128..3000}、mode∈{0..4}、`g_use_filt` 强制 0/1 两侧）；
// ④ `tools/compliance.py` 检查 REFERENCES/思路 两块位置与 URL。
// ★ 本发**单变量**（只动这一处），归因干净。
// ================
// ===== REFERENCES =====
// [A1] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#103830**
// <https://duck.ac/submission/103830>（6362.789607 ms）= 本文件的**正文基底**，
// `problems/1014/work/e4f_v2_body.cpp` 逐字复制。
// 该件 = #103629 + 复刻 [A2] 的 6 处 hunk + `orders` 相位除法改倒数乘-移位。
// [A2] duck.ac 用户 **saffah_codex_6s_agg2**，提交 **#103716**
// <https://duck.ac/submission/103716> —— 用途：见 #103785 的 REFERENCES（6 处 hunk）。
// [A3] duck.ac 用户 **saffah_codex_6s_agg2**，提交 **#103699**
// <https://duck.ac/submission/103699> —— 用途：**参考了思想**：该件自述用 16 格粗网格
// 过滤器 + 逐站点预取距离；本发在其取值上继续扫同一条轴。
// [A4] /home/yjp/duck.ac/BRIEF.md：`§2.18.742`（对手在一格一格爬轴 ⇒ 我们**一次扫完整条轴、
// 两端都取点**、直接跳到峰）、`§2.18.698`（看到对手动了一个**数值常量**，第一反应是扫整条轴）、
// `§2.18.684`（"只改常量/未系统扫过"的候选应当**当天扫掉**）、`§2.18.655`（"往某方向推"的
// 建议要先问那个方向是从最优点往哪边走）、`§2.18.669`（提交门槛 = 这一发是不是我当前最先进的实现）。
// ★ `g_Qc` 在全族（1012/1013/1014/1015/1016）都被 `tools/deadswitch.py` 标为
// "**never assigned anywhere; value is the initialiser**" ⇒ 它是**编译期常量**而不是运行期旋钮
// ⇒ **要扫它必须改源码字面量**，因此这条轴在本族**从来没有被真正扫过**。
// ======================
//
// ===== 思路 =====
// 进入时实测（`python3 tools/exact.py 1014`）：`mine = 6359.841840`（#103863），`T = 6370.850904`（#103716），
// **严格支** `thr = 0.99*T+1us = 6307.143` ⇒ 缺口 **52.699 ms**。
//
// ---- 本发只改两行：`#define MAXK 10 → 8`、`#define STRIDE 16 → 8` ----
// 题面（`problems/1014/statement.txt`）与接口都写明**本题是八维且仅八维**：
// `void count_8d(int n, const unsigned *x[8], unsigned *out)`，`n` 恒等于 1,000,000，
// 每维坐标在 `[0, n)`。⇒ **k 恒为 8，`MAXK`/`STRIDE` 取 8 是安全且紧的。**
// * `MAXK` 是 `s_bd[]`/`s_ord[]`/`s_gt[]`/`s_pos[]`/`s_posp[]`/`s_bnde[]`/`s_maxr[]`/
// `s_pofs[]`/`s_plen[]`/`s_fpos[]`/`s_bval[]` 的行距 ⇒ **10 → 8 直接把这些表的足迹砍 20%**；
// * `STRIDE` 是 `s_pm[]`（每点 16 个 u32 的记录）与 `avpool[]` 的行距 ⇒ **16 → 8 把 `s_pm`
// 从 64 MB 砍到 32 MB**，而 k=8 只用到前 8 个 lane（`for (int d = k; d < STRIDE; d++) p[d]=0;`
// 的填充循环随之消失）。
//
// ---- 定价依据（**判题机真数，同一引擎血脉**）----
// `problems/1014/notes.md` §3.2 记录了对手对本族做过的**纯旋钮消融**（判题机真数，
// 且他当时的 `MF_PFA/PFB` 档与我们的基座完全相同）：
// **`MAXK 10→8 + STRIDE 16→8` = −44 ms**（`#101160 → #101692`），
// 且 `mem 891732 → 852080 KB` **逐 KB 对得上**（`s_pm −31250` + `s_bd −7812` + `s_maxr −576` KB）
// ⇒ 这是一条**有判题机读数、且内存账自洽**的改动。
// ★ 同族 `1013` 上的同两项（`MAXK/STRIDE`）= **−62 ms**。
// ★ `§2.18.689`：同族旋钮刀的幅度取决于本题基座已扫过哪一档 ⇒ **本题基座的 `MAXK/STRIDE`
// 仍是默认 10/16（从未扫过）** ⇒ 红利应当完整在。
// ★ `§2.18.684`：这是"只改常量、`.text` 语义不变、零正确性风险"的那一类 —— 本族 notes
// 点名已久却从未做过 ⇒ 本发把它当天做掉。
//
// ---- 适用边界（**必须留档**，`§2.18.685`：判死要带重试条件，同理"适用范围"也要带）----
// ★★ **`MAXK=8` / `STRIDE=8` 只对 `k ≤ 8` 正确**：`k=9` 时 `query_group` 会读
// `LDU((const __m256i *)(rec + 8))`（`k > 8` 分支）与 `bp[]` 的第 9 项 ⇒ 越界读。
// **本题 k 恒为 8（题面+接口），所以成立**；本文件里的 `count_9d/count_10d` 入口
// 在本题不会被调用，但**若要移植到 k>8 的姊妹题（1015/1016 家族）必须先逐处核 `k > 8` 分支**。
//
// ---- 交前闸门 ----
// ① 本地 `-O2 -std=c++17 -DLOCAL_TEST` 编译通过；
// ② **k=8 全档** `brute_nd` 对拍（n∈{128,200,512,1000,3000}、mode∈{0..4}、`g_use_filt` 强制 0/1 两侧）
// 与 n=1e6/k=8 真形状 checksum **逐位相同**；
// ③ ★ 已知**预期差异**：本改动**只对 k=8 有效**，所以 k=9/k=10 的本地对拍**会失败**（这是设计边界，
// 不是缺陷）—— 闸门脚本本发只跑 k≤8，并在 notes 里写明；
// ④ `tools/compliance.py` 检查 REFERENCES/思路 两块位置与 URL。
// ================
// ===== REFERENCES =====
// [1] duck.ac user saffah_cc_v41_260924, rank-1 submission of this problem
// (see problems/1014/ref/01_*.cpp, e.g. 1012 #89064, 1014 #89074, 1015 #89078,
// 1016 #89165). Use: the whole file is a derivative of that account's MIN-FOLD
// engine for k-D dominance counting -- per-query fold dimension ds = argmin_d
// gt[d][i], bit space numbered by the rank in dim ds, per-(fold,dim,bucket) bitset
// pool, coarse-grid (s_fbm) filtered fringe scan, MF_ASCR scratch reuse, nibble
// popcount. Structure copied verbatim; the changes below are additive.
// [2] This workspace, agent-written file problems/1012/work/v4_orders.cpp
// (no duck.ac submission): the "start[v]" rewrite of the per-dimension orders
// phase -- gt[i] = start[x_d[i]] and bd follows from start[x_d[i]+1] -- replacing
// the random xd[od[j]] walk / random gt[], bd[] scatter. Copied verbatim.
// [3] This workspace, shared header tools/fastload.h (no duck.ac submission) and
// tools/gcc93check.py: under the judge's gcc 9.3 -O2 the pointers given to
// _mm256_load[u]_si256 cannot be proved 32-byte aligned, so each 32-byte load is
// emitted as vmovdqu xmm + vinserti128 (one extra load uop plus a port-5 uop).
// The inline-asm forms force a single vmovdqu ymm. Load version copied verbatim;
// the store version is used with the corrected "+m" volatile operand form.
// [4] /home/yjp/duck.ac/BRIEF.md section 2.5: the judge compiles with -O2 and no
// -march, so source-level #pragma GCC optimize/target is allowed; measured +3~6%
// on this family. Multi-threading is unavailable.
// [5] Judge Duck FAQ <https://duck.ac/faq> (flags, machine: i3-8100, 6 MB L3, 2x4 GB
// DDR4-2400, memory accounted as pages written).
// [6] duck.ac user saffah_cc_v41_agg1, submission #97898 (problem 1012)
// <https://duck.ac/submission/97898>
// local archive: problems/1012/work/agg.cpp == problems/1013/work/w2_fringe.cpp
// Use: copied verbatim the [W2] rewrite of the fringe coarse-cell index -- the
// per-(query,dim,other-dim) 64-bit division `(u64)rec[e] / (N/Qc+1) + 1` replaced by
// the multiply-shift `((u64)rec[e] * Qc) >> g_cshift`, with the boundary ladder
// `s_bnde` moved onto the same shift (`CBOUND(c) = (c << g_cshift) / Qc`) and
// g_cshift = ceil(log2 N) fixed before any boundary is computed. That submission
// measured it at -13.5% on the fringe phase on the real judge, and the same
// rewrite measured -202.4 ms on problem 1013 (#101853).
// ======================
// [1] duck.ac 用户 saffah_codex_6s_agg2, 提交 #103119 <https://duck.ac/submission/103119>
// （本题实时 T = 6528.657194 ms）。用途：**直接复制**了该提交相对本账号 #101863 的**全部实质差异**
// （`tools/rivalcopy.py 1014` 报 fwd 97.65% / rev 98.98%，SAME PROBLEM；
// `tools/fnmd5.py` 报 **per-function SAME=16 DIFF=1（只有 query_group）、全量 22 行全在函数体内、
// 无 MISSED**）——共 4 处：
// ① `#pragma GCC target(...)` 末尾加 `,tune=skylake`；
// ② 新增**查询级 `avpool` 预计算**（`Av[]` 从每 (query,dim) 对重算一次 → 每 query 一次）；
// ③ 站点③预取块加一行 `avpool + j*STRIDE`；
// ④ 扫描里 `u32 Av[STRIDE]` 三条构建语句 → `const u32 *Avp = avpool + i*STRIDE` + 两处 LDU 改址。
// ★ 与 #103058（1016b）的 4 处**逐处同型** ⇒ "一招通吃一个族"（本族 1012/1013/1014/1015/1016 都被同一刀翻红）。
// ===== 思路 =====
// 本发（b4_1014_v1）：把对手 #103119 的 4 处差异**原样搬回**（见 REFERENCES [1]），其余一字未改。
// 依据：对手 T=6528.657194 vs 我方 #101863=6634.406733 ⇒ **−105.750 ms（−1.59%）**。
// ★ 同一刀在 1016b 上已由我实测 = −4.497 ms（闸门 105/105 逐位）⇒ 机理已验证，本发是**跨题复用**。
// 正确性：`Av` 数据逐元素相同（同 `BVAL(e,bdi[e])`、同 0x7FFFFFFF 掩码），只是**从栈上重算改为常驻缓冲复用**
// ⇒ 逐位等价；闸门见下。
// The algorithm is untouched; this submission stacks only changes that were each
// measured separately on the real judge, on top of the rank-1 source:
// (a) file-scope #pragma GCC optimize("O3,unroll-loops,no-strict-aliasing") and
// #pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt"), because the judge passes
// only -O2 and no -march (ref [4]).
// (b) the sequential orders phase (ref [2]): the previous code, for every point and
// dimension, walked od[] and read xd[od[j]] at random, then wrote gt[p] and
// s_bd[p*MAXK+d] at random. The value-indexed formulation computes the same
// numbers with two sequential sweeps of xd[] and lookups in a (n+1)-entry table.
// (c) 32-byte loads through the inline-asm helpers (ref [3]) so that gcc 9.3 emits a
// single vmovdqu ymm instead of vmovdqu xmm + vinserti128.
// (d) the [W2] fringe coarse-cell index by multiply-shift (ref [6]). The fringe's
// filter row lookup did one 64-bit division per (query, other-dim) -- k-1 = 7 per
// query at k=8 -- and a 64-bit `div` is a multi-uop, non-pipelined instruction that
// sits directly on the "value -> cell -> row address -> load" dependency chain.
// Both the cell index and the boundary ladder now share one shift, so the filter
// stays the same tight superset. This is an instruction-selection change only;
// the cell width moves by <5% for n=1e6 (31251 -> 32768), so the candidate count
// is essentially unchanged.
// Verified by checksum equality against the unmodified engine on local random data
// (several k and several value distributions, including heavy ties) before submitting.
// ================
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#define MF_QCOEF 70
// 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 long long u64;
#ifndef MAXN
#define MAXN 1000005
#endif
#define MAXK 8
#define STRIDE 8
// 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_QCOEF
#define MF_QCOEF 70
#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 32
#endif
#define MF_PFB1 64
#define MF_PFB2 64
#define MF_PFB3 8
#define MF_PFB4 32
// 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 u32 s_ord[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u32 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_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 = 16;
static int g_cshift = 20; // cell(x) = (x*g_Qc)>>g_cshift ; boundary(c) = (c<<g_cshift)/g_Qc
#define CBOUND(c) ((u64)(c) << g_cshift) / g_Qc
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_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];
STU((__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, acc = z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 4 <= Wr; w += 4) {
__m256i a = LDU((const __m256i *)(bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU((const __m256i *)(bs[d] + w)));
__m256i lo = _mm256_and_si256(a, lm);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc = _mm256_add_epi8(acc, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
if (++cnt == 31) { tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z)); acc = z; cnt = 0; }
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, 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 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_PFB1 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB1] * 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 u32 *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 = s_gt[ds][s_qord[t]]; if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 1;
if (nwd <= AS_WORDS) {
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 u32 *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);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][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 u32 *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 u32 *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 = s_gt[ds][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;
// Reuse each query's bucket-boundary vector across all non-fold dimensions.
static u32 *avpool = 0; static size_t avcap = 0;
{
size_t needav = (size_t)N * STRIDE * 4 + 256;
if (needav > avcap) { avpool = (u32 *)pool_alloc(needav); avcap = avpool ? needav : 0; }
if (avpool) for (u32 qi = q0; qi < q1; qi++) {
u32 ii = s_qord[qi];
const u32 *bdi = s_bd + (size_t)ii * MAXK;
u32 *av = avpool + (size_t)ii * STRIDE;
for (int e = 0; e < k; e++) av[e] = (e == (int)ds) ? 0x7FFFFFFFu : BVAL(e, bdi[e]);
for (int e = k; e < STRIDE; e++) av[e] = 0x7FFFFFFFu;
}
}
// ---- 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;
#if !MF_FR_NOSORT
{
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 = LDU((const __m256i *)hiA);
const __m256i notexloh = LDU((const __m256i *)hiAh);
const __m256i all8 = _mm256_set1_epi32(-1);
const u32 *const odd = s_ord[d];
const u32 *const sdmpd = s_sdmp[d];
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
#if MF_PFB > 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);
_mm_prefetch((const char *)(avpool + (size_t)j * STRIDE), _MM_HINT_T0);
}
}
#endif
const u32 *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;
const u32 *Avp = avpool + (size_t)i * STRIDE;
__m256i qr = _mm256_or_si256(LDA((const __m256i *)rec), khi8v);
__m256i Avm = _mm256_or_si256(LDU((const __m256i *)Avp), notexlo);
__m256i qT = _mm256_min_epi32(qr, Avm);
__m256i qTh;
if (k > 8) {
__m256i qrh = _mm256_or_si256(LDU((const __m256i *)(rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU((const __m256i *)(Avp + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
} else {
qTh = all8;
}
u32 add = 0;
const u32 *c0 = sdmpd ? (sdmpd + (size_t)lo * g_KS) : 0;
if (g_use_filt && KK == 8) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
const u32 csh = (u32)g_cshift, cQc = (u32)g_Qc;
const u32 FW = g_fbm_W;
const u64 *rows[7];
int nr = 0;
for (int e = 0; e < 8; e++) {
if (e == d) continue;
u32 cell = (u32)(((u64)rec[e] * cQc) >> csh) + 1;
if (cell > cQc) cell = cQc;
rows[nr++] = s_fbm + ((size_t)(d * 8 + e) * (g_Qc + 1) + cell) * FW + w0;
}
const u64 *r0=rows[0], *r1=rows[1], *r2=rows[2], *r3=rows[3];
const u64 *r4=rows[4], *r5=rows[5], *r6=rows[6];
for (u32 w = 0; w < nw; w++) {
u64 vv = r0[w] & r1[w] & r2[w] & r3[w] & r4[w] & r5[w] & r6[w];
if (w == 0) vv &= ~0ull << (lo & 63);
u32 rem = hi & 63u;
if (rem && w + 1 == nw) vv &= (1ull << rem) - 1;
#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) * g_KS) : (s_pm + (size_t)s_ord[d][j] * STRIDE);
__m256i pv = LDU((const __m256i *)cp);
add += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
}
out[i] += add;
continue;
}
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u32 csh = (u32)g_cshift, cQc = (u32)g_Qc;
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] * cQc) >> csh) + 1;
if (c > cQc) c = cQc;
const u64 *row = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
if (first) {
u32 w = 0;
for (; w + 4 <= nw; w += 4) STU(fb + w, LDU(row + w));
for (; w < nw; w++) fb[w] = row[w];
first = 0;
} else {
u32 w = 0;
for (; w + 4 <= nw; w += 4) STU(fb + w, _mm256_and_si256(LDU(fb + w), LDU(row + w)));
for (; 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) * g_KS) : (s_pm + (size_t)s_ord[d][j] * STRIDE);
__m256i pv = LDU((const __m256i *)cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU((const __m256i *)(cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
continue;
}
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 = LDU((const __m256i *)(c));
__m256i p1 = LDU((const __m256i *)(c + g_KS));
__m256i p2 = LDU((const __m256i *)(c + 2 * g_KS));
__m256i p3 = LDU((const __m256i *)(c + 3 * g_KS));
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 * g_KS;
}
for (; j < hi; j++, c += g_KS) {
__m256i pv = LDU((const __m256i *)c);
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
add = a0 + a1 + a2 + a3;
} else {
const u32 *c = c0;
for (u32 j = lo; j < hi; j++, c += g_KS) {
__m256i pv = LDU((const __m256i *)c);
__m256i pvh = LDU((const __m256i *)(c + 8));
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8) &
(u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
add += ok;
}
}
} else {
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__m256i pv = LDA((const __m256i *)c);
add += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
}
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;
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;
// [e4f] exact reciprocal for the orders-phase floor(x/B): with 2^divS >= B*N and
// divM = ceil(2^divS/B), floor(x*divM/2^divS) == floor(x/B) for every x < N
// (proof: x*divM/2^divS < x/B + 1/B <= floor(x/B) + 1). N < 2^20 and B <= N, so
// divS <= 40 and x*divM < 2^60 fits u64.
u64 divM; int divS;
{ u64 t = (u64)B * (u64)N; divS = 0;
while (divS < 64 && ((u64)1 << divS) < t) divS++;
divM = ((u64)1 << divS) / B + 1; }
#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
// [W2] the coarse-cell ladder must be fixed up-front: s_bnde (boundaries in
// value space) and the per-query cell index must use the SAME shift, otherwise
// the filter stops being a tight superset (measured: 7.5x more candidates).
if (g_use_filt) { int sh = 0; while ((sh < 32) && (((u64)1 << sh) < (u64)N)) sh++;
g_cshift = sh; }
{
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 = LDU((const __m256i *)lo8);
khi8hv = LDU((const __m256i *)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];
// Counting sort + rank record in sequential passes: start[v] = #{q : x_d[q] < v}
// is exactly gt[] of every point with x_d == v, and end[v]-1 = start[v+1]-1 is the
// last rank of that tie group, so bd[] follows from start[v+1] too. Both are value
// lookups, so the per-point work is a sequential sweep of xd[] instead of the
// previous random xd[od[j]] read plus random gt[]/bd[] writes.
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_tmp[i] = s; s += c; } }
{ u32 *pp = s_pos[d];
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i; pp[i] = q; } }
{ u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
gt[i] = s_cntv[v];
u32 bb = (u32)(((u64)(s_cntv[v + 1] - 1u) * divM) >> divS);
s_bd[(size_t)i * MAXK + d] = bb > Q - 1 ? Q - 1 : bb;
} }
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;
}
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;
}
// ---------- 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;
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 = CBOUND(c);
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 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).
if (g_plen_on) {
const u32 *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_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB4] * 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;
}
}
}
}
}
// ---------- 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++) 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) { 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((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((u32)n, (const unsigned **)x, out, 6); }
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); }
#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
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 6.296 s | 583 MB + 544 KB | Accepted | Score: 100 | 显示更多 |