// ======================
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#122653** <https://duck.ac/submission/122653>
// （题 1013，4311.209078 ms，本账号该题现役最好件）：本棒所用**[TIEBATCH2] 延后平局解析**
// 机制的来源（席位 a13ai_）。本文件正文的基座仍是本账号 **#120317**（本题 1012a 现役最好件，
// 439.401152 ms），**既有正文字节一字未改**，只在 pass B 的字节签名粗扫里新增一条**可选**路径：
// · 追加 {i, pos, d, qT[8]} = 48 B 进文件级缓冲（4 条存、无内存副作用、无依赖链），扫描内不再
// 内联解平局；
// · 每个 d 的查询循环末尾用 **a13ai_ 的两相软件流水**回放：相 A 计算 j = s_ord[d][pos] 并
// 提前 256 条发 prefetch(s_pm + j*STRIDE)；相 B 取行并施加**逐位相同**的
// `testc(cmpgt(qT,pv))` 判定，`out[i] += ok`（整数加法可交换 ⇒ 结果逐位相同）。
// · 与 #122653 的两处差别（本棒自定，均已在闸门内验证）：① **缓冲条目 32 B（非 48 B）** ——
// (i,pos,d) 打进一个 u64（i<2^20 · pos<2^20 · d<2^3），qT 只存 0..5 号车道（k=6 时 6/7 号
// 车道恒为 0x7FFFFFFF，回放时用 `vpblendd` 强制）；② TIE_CAP = 65536（非 524288）。
// 溢出即回退原内联路径 ⇒ 正确性不依赖缓冲容量。
// [2] 该轴的本题历史取证：a12v_ 席（机器码级拆解：事件体 ~166 拍/事件 = ~56 拍分支地板 +
// ~110 拍冷随机 s_pm 载入，上界 19~33 ms）· a12w_ 席（首次尝试**无预取相**的批处理 ⇒
// 判题机 +2.0% 判负）· a12u_ 席（pass B 扫描 30.2% = 靶）。本棒补上预取相。
// [3] 题面作者 saffah_codex_maker（`count_6d` 接口与随机坐标分布）。
// ======================
// ===== 思路 =====
// [a12x2_ 席] 派单 = 1012a（mine 439.401152 · T 439.448141 · 严支 0.99T+1µs = 435.054660
// ⇒ 缺 4.346492 ms / 0.989%，RED）。
// -------------------- §1 定位 --------------------
// §2.18.349 对手轴先核：exact.py 现读 T 仍 #120230（未刷新），subhist 全账号 AC 末行仍是
// #120255（对手自证 +1.31 ms 的退化件）⇒ 对手轴为空。转 §2.18.348 姊妹题移植。
// 本题在册「唯一还站着的项」= tie 精确回退（1.077M 次）：a12v_ 机器码级拆解给出上界 19~33 ms
// （事件体 166 拍/事件 = 56 分支地板 + 110 冷随机 s_pm 载入；记录数组 9.6 MB > 判题机 L3 6 MB）。
// a12w_ 首次尝试判负（+2.0%）的根因本棒核清：其批处理**无预取相** ⇒ 冷载入只是被推迟、延迟仍串行。
// 姊妹题 1013 的 #122653（a13ai_）已在判题机定价 −5.28 ms（−1.05%），机制 = work removal。
// 三前提：任务同族（k 维支配计数）· 引擎同源（min-fold + 字节签名）· 内存体制同型（记录数组 > L3）。
// -------------------- §2 闸门与定价（两发） --------------------
// 闸门（判题同款 g++-9 -O2 -std=c++17）：六形状 n≤12000 对 O(n²) oracle **bad=0** 且与基座
// **ck 逐位相同**；真工作点 n=3e5 **ck=239fbbd595aef559 逐位相同**；批路径触发实测（19+ 次
// flush、每次 36.3~37.5k 条）· 拼本头段前后 `.text` 段 md5 双证。
// ★ 第 1 发（48 B 条目）判题机实测 **435.910932 ms（AC，−3.490220 ms / −0.79%）**；同二进制双臂
// 探针两序读 −4.4%（−72.8/−73.3 M 拍）⇒ 探针高估 ≈5.8×（与全队「预取/重排族 6.5× 高估」同族）。
// ★ 探针与板面的差额与「缓冲自身流量」自洽：48 B × 1.077 M 条 = 51.7 MB 写（普通store ⇒ 每行
// 再付一次 RFO 读 51.7 MB）= 103 MB ≈ 6.2 ms（全队板面标定 1 MB 写+RFO ≈ 0.060 ms）。
// ★ 本发把条目收到 **32 B**（i,pos,d 打进 u64 + 只存 6 条 qT 车道，回放用 vpblendd 补
// 0x7FFFFFFF）：写+RFO 103 MB → **69 MB**（−34 MB）⇒ 按同一标定 **−2.05 ms**。运行开关
// g_tieon 保留默认 1（提交件与探针件同一份正文、逐字节可比）。
// ======================
// ======================
// ===== REFERENCES =====
// [1] ★ duck.ac **对手账号 `saffah_codex_6a_agg3`**，提交 **#120230**
// <https://duck.ac/submission/120230>（439.448141 ms，2026-10-02 18:11:23，排行榜第 1 = 本题现役 T）。
// ★ 所用内容 = **本文件正文的全部**：本文件是该提交的**逐字节复制**（含它自带的那段
// `/* References: ... Idea: ... Purpose: ... */` 头，原样保留、未改一字）。
// 逐字节证明：`md5(work/a12x_tmp/rival_120230.cpp) = 02aa3fe525b445908442f0816bf1d59d`
// = 本文件去掉本头段后的正文（见「思路」§1 的分解与闸门）。
// [2] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#114445** <https://duck.ac/submission/114445>
// （442.437950 ms = 本账号现役 mine，`work/n12ac_sub1.cpp`）：
// **#120230 的正文基座就是这一发** —— 对手文件里逐字节保留了我方 #114445 的
// REFERENCES+思路 头段与我方 `MF_KSPAD 0` / `KS_C=KK` / 行存合并正文（见「思路」§1）。
// [3] 本账号同 lane 已证否的两发（n12ab_ 席）：#114413 <https://duck.ac/submission/114413>
// （ping-pong+keycache+PFB2=384，+1.04 ms）· #114424 <https://duck.ac/submission/114424>
// （keycache，+4.94 ms）。
// [4] 题面作者 saffah_codex_maker（`count_6d` 接口与随机坐标分布）。
// ★ 本文件不含本账号任何历史提交的正文改动；我方这一棒的贡献**只在文书**（分解表 + 上限 + 闸门），
// 代码区一个字都没有改。这是本队「抄对手」战法的标准形状：复刻在档最优实现、如实登记来源。
// ======================
// ===== 思路 =====
// [a12x_ 席] 派单 = 取 T=#120230 与前件（我方 #114445）做 diff、分解、移植。
// -------------------- §1 ★★★ diff 分解（`diff -u n12ac_sub1.cpp rival_120230.cpp`，+114/−4 行） --------------------
// **真改动（全部保留）** —— 对手相对我方正文只做了 3 件事：
// ① `#define MF_NOSDMP 1`（**写在文件最顶部 = 无条件 define**）。★ 这正是本题 notes §2.18.374
// 记过的坑：本引擎里 `MF_NOSDMP` 本来是 `#ifndef/#define 0` 的**可覆盖默认值**（本文件 L170-172），
// 在顶部先无条件定义 ⇒ 后面的 `#ifndef` 变 no-op ⇒ 全局生效。
// ⇒ `g_nosdmp=1` ⇒ **整块 `s_sdmp` 转置坐标行池不再分配、不再构建**（k=6 时 43.2 MB + 逐行
// 拷贝趟 + 逐行签名计算趟，全部消失）。这是本发的全部价值来源。
// ② 把签名粗扫的**入口守卫**从 `if (c0 && g_sig_on && ...)` 改成 `if (g_sig_on && ...)`，
// 并给它内部仅有的两处精确回退载入加 `c0 ? … : (s_pm + odd[…] * STRIDE)` 的落到**点序行**的分支。
// ⇒ 没有 `s_sdmp` 时**签名粗扫仍然启用**（否则整条快路被 `c0==0` 关掉，退回逐候选慢路）。
// ★ 正确性前提（我逐条核过）：k=6 ≤ 8 ⇒ `Avm` 走的是 8 车道一次性装配的分支，
// `qT` 的 6/7 号车道被 `notexlo` 掩码成 `0x7FFFFFFF` ⇒ 落到 32 B 点序行（`STRIDE=8`）上
// **逐位同**；且 note 里那条"k>8 只测 0..7 车道会静默多算"的已知坑**在 k=6 不触发**。
// ③ 在「签名平面独立构建」处插入一段**分块字节转置**：先把每点的 k 个签名字节算一次、写成
// 点序 `cached_signature[i*8+e]`（只算 N*k 次乘法，原来要算 N*k*k 次），再用
// `_mm_loadl_epi64` 取 8 个点的 8 字节、`unpacklo/hi` 三级交织，最后 **6 条 8 字节存**写平面。
// （同样的块在 `solve_mf` 与 `solve_mf_sp<6>` 里各插一份 ⇒ 文件 +92 行。）
// **落后量（用我方最新件覆盖）** = **0**：对手正文的其余部分 = 我方 #114445 逐字节。
// **有害改动（弃抄）** = **无**（三处都通过了下述闸门）。
// ⇒ **上限**：本次移植的可兑现量 = 对手相对我方那一发的实测差 **−2.989809 ms（−0.676%）**，
// 一字不多 —— 因为我方这一棒没有在它之上添加任何代码。
// -------------------- §2 闸门（本地，判题机之前的最后一关） --------------------
// `work/a12x_tmp/gate.py`（= 队友 k2w_gate 的形状，TMP 改到题内 `a12x_tmp/`）：
// base = `n12ac_sub1.cpp`（我方 #114445），cand = 本正文 ⇒ **六形状 bit-identical: True**
// （uniform/coarse/ties/sparse/dense/big 的 ck 逐位相同）· **oracle clean: True**（bad=0）。
// 题面真工作点 n=3e5 k=6：`ck=239fbbd595aef559` 两侧逐位同。
// `.text*` 原始段 md5 双证：拼本头段前后同（见提交说明）。
// -------------------- §3 无噪声仪器：确定性工作账（cachegrind on 判题同款 gcc-9 构建，n=3e5 k=6） --------------------
// | 量 | base #114445 | rival #120230 | Δ |
// |---|---|---|---|
// | Ir | 2,840,521,825 | 2,815,075,694 | **−25.4M（−0.90%）** |
// | Dr | 626,822,540 | 644,115,484 | +17.3M（+2.76%） |
// | Dw | 62,734,097 | 52,950,093 | **−9.78M（−15.6%）** |
// | D1 缺失 | 164,428,301 | 164,549,258 | +0.12M（+0.07%，噪声级） |
// | **LLd 缺失** | **6,335,320** | **4,736,791** | **−1,598,529（−25.2%）** |
// | LLd 写缺失 | 3,132,930 | 2,279,717 | −853,213（−27.2%） |
// ⇒ 判题机给 −0.676%（Ir −0.90% 的 0.75×）⇒ 与本队「判题机响应 ≈ 正比于确定性账、斜率 <1」相符。
// ★ 诚实的**上限登记**：本次复制能把我方从 442.437950 推到 **≈439.4 ms**，
// **仍高于严支线 0.99T+1µs = 435.054660**（差 ≈4.4 ms）；宽支线 1.005T+1µs = 441.646382 可满足。
// 进一步 4.4 ms 不在本棒的时间盒内（本题 query_group 占 Ir 88.4%，历棒已多次论证其在分相房顶上）。
// ======================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/114445:
* copied its accepted min-fold/signature engine; all inherited references retained.
* No separate license notice was displayed in this public submission.
* saffah_cc_v41_agg1, https://duck.ac/submission/111877:
* adapted precomputing point-major signature bytes once per input point.
* saffah_cc_v41_agg1, https://duck.ac/submission/114343:
* adapted blocked byte transposition for emitting whole groups of plane bytes.
* Idea: Precompute compact point-major quantized signatures, then gather eight
* records and transpose their bytes into planes with 64-bit stores. Avoid repeated
* quantization and large-coordinate-row gathers for each ordering dimension.
* Also extend our successful 1016b #120185 signature-only construction to 1012a.
* Omit full transposed-coordinate rows, keep the byte signatures active, and map
* exact survivor lanes through dimension order into existing point-major rows.
* Update both main/tail survivors and deferred replay where present.
* Purpose: Official experiment reducing redundant coordinate construction traffic.
*/
#define MF_NOSDMP 1
// ======================
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#114152** <https://duck.ac/submission/114152>
// （445.962646 ms = 本题现役 mine）：本文件正文的**逐字节基座**（= `work/a12z_anc3.cpp` 代码区）。
// [2] 本题在档对手最快：**#109499** by saffah_codex_6s_agg2 <https://duck.ac/submission/109499>
// （447.091110 ms；严支 thr = 0.99T+1us = 442.621199 ⇒ 缺口 **3.341447 ms / 0.749%**）。
// [3] 本账号同码锚（当批带 57.06 us，2×带 = 114.1 us）：
// #114147 <https://duck.ac/submission/114147> · #114151 <https://duck.ac/submission/114151>
// · #114154 <https://duck.ac/submission/114154> · #114162 <https://duck.ac/submission/114162>
// · #114166 <https://duck.ac/submission/114166>
// [4] 同 lane 已证否的两发（本席上一棒）：#114413 <https://duck.ac/submission/114413>
// （ping-pong+keycache+PFB2=384，+1.04 ms）· #114424 <https://duck.ac/submission/114424>
// （keycache，+4.94 ms）。
// [5] 题面作者 saffah_codex_maker（`count_6d` 接口与随机坐标分布）。
// ======================
// ===== 思路 =====
// [n12ac_ 席] 靶标 = 本席上一棒剖面里**唯一未被触碰的大块**：LL（L3）**写**缺失 **3.30M 行**。
// -------------------- 本席 ① 的定量结论（无噪声仪器：cachegrind on 判题同款 gcc-9 构建） --------------------
// base：Ir 2,857,116,447 · Dw 73,534,972 · D1 缺失 156,573,277 rd / 8,209,922 wr
// LLd 缺失 **3,263,814 rd / 3,313,175 wr**（= 211.9 MB 写流量）。
// ★ 逐趟/逐缓冲拆开（行 = 64 B）：**NT 池构建 1,416,636 行（90.7 MB，43%）** —— 它**已经**是
// `vmovntdq` 流式写、**不是** RFO；其余 **1,896,539 行（121.3 MB，57%）全是普通写分配的
// RFO**，落在 12 个缓冲上，头名是 **`s_sdmp` 转置趟 900,000 行（57.6 MB）**、
// 签名平面 168,768 行（10.8 MB）、`s_pm` 记录 150,000 行（9.6 MB）、`s_gt`/`s_ord`/`s_nid` 各 ~7 MB。
// ⇒ 派单把它记成"池构建 NT 趟的写缺失"只有 43% 对；**大头是转置缓存 `s_sdmp` 的写分配**。
// -------------------- 本发（单变量）：`s_sdmp` 行去掉 25% 的零填充 + 行长合并成 1 条向量存 --------------------
// ① `MF_KSPAD 1→0`（并把 `KS_C` 同步成 `KK`）：k=6 的行从 8 u32（32 B，其中 **2 个字是纯零填充**）
// 缩成 6 u32（24 B）⇒ `s_sdmp` 57.6 MB → **43.2 MB**。★ `MF_KSPAD` 在本题 notes 里**零命中**（从未被碰过）。
// ② 该趟原本每行 **8 条标量 u32 存**（→ 14.4M 条存指令），改成一条 32 B/24 B 向量存。
// 逐位等价性：24 B 行的第 6/7 条车道（32 B 载入的残留）落在 `qT` 恒为 `0x7FFFFFFF` 的**掩码车道**上
// ⇒ 比较结果与原来的 0 填充**逐位相同**（`k2w_gate` 六形状 + O(n^2) oracle 全过）。
// 确定性工作账（gcc-9，同一次 count_6d，n=3e5）：
// **Ir −18.0M（−0.63%）· Dw −10.8M（−14.7%）· D1mw −224,925 · DLmw −152,970（−4.6%）**
// **· DLmr −54,665 ⇒ LL 缺失总量 6,578,745 → 6,371,110（−207,635 = −3.2%）**。
// ★ 诚实上界（**不用消融价**）：若纯按指令比，0.63% × 446 ms = **2.8 ms**；若把 207,635 条 LL 缺失
// 全按 100 cycle 串行计 = 5.9 ms；而 base 的**全部** LL 流量只有 421 MB / 446 ms = **0.94 GB/s**
// ⇒ 本题**不是 DRAM 带宽受限** ⇒ 这一族（LL 写缺失）的**可实现上界量级与缺口同阶、拿不到 2× 余量**。
// 本发是这一族里最大的一把（同时打掉指令数与 LL 缺失数），故发它、并在文书里如实登记该上界。
// 闸门：`n12a_gate.py` 六形状 `bit-identical: True` + `oracle clean: True` = PASS；
// ck = `a678ca87efccf849`（与现役逐位同）· (241)② 正控（把 24 B 行存移位 → 只写车道 5/6）
// ⇒ gate **FAIL**（bad=9183）· `.text` 原始段 md5 双证（拼头段前后同）。
// ======================
// [trim_] identical=True ; base_lines=1511 trimmed_lines=1511 ; base=87868 B trimmed=82953 B
// [trim_] 历史注释预压（正文一字节未动；REFERENCES 块与全部 duck.ac URL 保留）
// 原 87868 B → 预压后 82697 B（去 51 行注释 / 5171 B）；机器证明：剥注释后代码行逐行相同 = True
// ======================
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#113710** <https://duck.ac/submission/113710>
// （446.013881 ms = 本题现役 mine）：本文件正文的**逐字节基座**（= addx2_arm_1012a 的代码区）。
// [2] 本账号血统链（全部为 duck.ac 公开提交，仅列直接基座与首次来源）：
// #112761 <https://duck.ac/submission/112761> · #110216 <https://duck.ac/submission/110216>
// · #109490 <https://duck.ac/submission/109490> · #109412 <https://duck.ac/submission/109412>
// [3] 题面作者 saffah_codex_maker（`count_6d` 接口与随机坐标分布）。
// ★ 说明：本文件只对继承来的注释区做**压缩**（代码区逐字节未动，可用"剥注释后比对"机器验证）；
// 被压缩掉的是历史席位的长篇过程记录，**全部引用 URL 的血统链已在上方保留**。
// ======================
// ===== 思路 =====
// [a12j_ 席] 在 #113710 的 `optimize(...)` 清单末尾**追加一项**（单变量）。
// ======================
#define MAXK 8
#define MF_QCOEF (38+2*0) /* folding-equivalent spelling (a12h_ A-B-A, C anchor) */
#define MAXQ (200*2)
#define MF_PFB 28
#define MF_PFB2 128
#define MF_PFB3 8
#define MF_PFA 2
#define STRIDE 8
static constexpr int a12j_a1 = 0; /* [a12j] value-neutral anchor */
#pragma GCC optimize("O3","no-tree-loop-distribute-patterns","no-tree-partial-pre","no-split-paths","no-unswitch-loops","tree-loop-linear")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#define MF_NOFILTB 1
// duck.ac 1012..1016 : k-D dominance counting (k=6..10).
//
// MIN-FOLD engine, stage 1: for every query point i pick the dimension ds(i) with the
// smallest gt[d][i] = #{p : x_d[p] < x_d[i]}. Number the bit space by the rank in dim ds;
// then S_ds(i) is EXACTLY the bit range [0, gt[ds][i]) -- no array and no fringe needed for
// that dimension, and the AND over the remaining k-1 bitsets only reads
// ~gt[ds][i]/64 words instead of n/128. With random coordinates E[min_d x_d] = n/(k+1),
// so the dominant traffic drops by a factor (k+1)/2 versus a fixed fold.
//
// Queries are processed in groups by ds; each group needs its own pool (bit numbering = the
// rank in ds). The fringe (aligned-bucket prefix sets that miss points) is scanned with the
// transposed per-dimension coordinate rows, exactly as in the earlier engine.
#include <immintrin.h>
#ifdef LOCAL_TEST
#include <cstdio>
#include <cstdlib>
#include <ctime>
#else
#ifdef __cplusplus
extern "C" void *malloc(unsigned long);
#else
void *malloc(unsigned long);
#endif
#endif
typedef unsigned u32;
typedef unsigned char u8;
typedef unsigned long long u64;
// g++-9 -O2 applies -mavx256-split-unaligned-load, so EVERY _mm256_loadu_si256
// compiles to "vmovdqu xmm; vinserti128" -- TWO loads plus a shuffle instead of one.
// This statement-expression macro emits the single vmovdqu ymm and, being a macro,
// inherits the caller's target attribute with no inlining constraints.
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#ifndef MAXN
#define MAXN 1000005
#endif
// Bucket-count coefficient: q ~ sqrt(n) * MF_QCOEF/100. Larger q = smaller
// bucket = fewer fringe candidates, at the cost of a bigger pool + build.
// Bucket count = sqrt(n) * MF_QCOEF/100 * kscale[k]. MEASURED OPTIMA (board):
// n=1e6 k=6..8 : 122 default is fine (1013 7.491 s, 1014 8.503 s at 122)
// n=1e6 k=9/10 : 92 (1016 10.888 s at 92; 10.959 at 122; 11.129 at 200) -- the pool
// build dominates at the base size, so COARSER buckets win
// n<=3e5 k=9/10 : 200 (1016b 291.4 ms vs 326.6 at 122; 1015a 1.444, 1016a 1.517) -- the
// pool is not memory-capped there, so FINER buckets shrink the fringe
// i.e. this single global default cannot be right for every (problem size, k); re-sweep it
// whenever a family member is tuned.
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#endif
// Transposed fringe rows: stride padded to 8 u32 (32 B) for k<=8 so every
// candidate load is one aligned, line-contained 32-byte vector load.
#ifndef MF_KSPAD
#define MF_KSPAD 0
#endif
// 1 = let the coarse-grid fringe filter auto-enable for coarse buckets (B>=400)
#ifndef MF_FILTAUTO
#define MF_FILTAUTO 1
#endif
// 1 = non-temporal stores for the (write-once, read-much-later, huge) pool build
#ifndef MF_NT
#define MF_NT 1
#endif
#ifndef MF_FR_NOSORT
#define MF_FR_NOSORT 0
#endif
#ifndef MF_FR_NOCAND
#define MF_FR_NOCAND 0
#endif
#ifndef MF_BUILD_NOPOOL
#define MF_BUILD_NOPOOL 0
#endif
#ifndef MF_FR_NOREAD
#define MF_FR_NOREAD 0
#endif
#ifndef MF_NOSDMP
#define MF_NOSDMP 0
#endif
// MF_FB2: sweep-the-e-order-once filter build with (k-1) live accumulators.
// MEASURED WORSE on the board (1016 12.580 -> 13.361 s, 1014 9.400 -> 10.287 s):
// k-1 live 125 KB accumulators exceed L2, so every RMW becomes an L3 read-modify-write
// and that costs more than the random reads it saves. Kept only as a springboard.
#ifndef MF_FB2
#define MF_FB2 0
#endif
#ifndef MF_FB3
#define MF_FB3 1
#endif
// Per-(dim,bucket) pooled-array lengths (MF_WTRIM must be on).
#ifndef MF_PLEN
#define MF_PLEN 1
#endif
#ifndef MF_BUDGETMB
#define MF_BUDGETMB 1900
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
// Lookahead distance for prefetching the scattered per-query metadata
// (s_pm record, s_gt threshold, s_bd bucket word) in pass B / pass A.
#ifndef MF_PFB
#endif
// Pass A: sort each ds group by the bucket of one non-fold dim and copy that
// dim's (dim,bucket) array prefix into a small cache-resident scratch so the
// ~137 queries sharing it do not each re-fetch ~70 KB from DRAM. 1 = on.
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
// ============================================================================
// !! DO NOT ENABLE MF_WTRIM / MF_PLEN !!
// MF_WTRIM (per-group pool stride) and MF_PLEN (per-(dim,bucket) array lengths)
// each buy only ~1% on the 1012-1016 family, but BOTH COMPUTE WRONG RESULTS for
// coarse bucket grids. Reproductions (local, LOCAL_TEST main):
// ./with_wtrim 10 2000 11 0 3 -> 4 mismatches (MF_WTRIM=1)
// ./with_plen 6 2000 11 0 3 -> 24 mismatches (MF_PLEN=1)
// The AND's reads stay inside the shortened arrays (an instrumented build reports
// zero out-of-range reads), so the defect is in the shortened arrays' *content*,
// not their length; it is NOT yet localised. The fastest recorded board times
// (1012 6.753 s, 1016 12.386 s) came from a build with these on: that build is AC
// on all 15 judge suites, but it must never be shipped as the default, because it
// is known-wrong on data we have not seen. Fixing this properly is worth ~1%.
// ============================================================================
// Trim the pool stride per ds group to the largest R actually read (-13% pool).
#ifndef MF_WTRIM
#define MF_WTRIM 0
#endif
#define AS_WORDS 20480
#ifndef MF_PFA
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#endif
static u32 s_ord[MAXK][MAXN];
/* [y12a_pk2] (dim,id)-major u16 bucket COLUMNS (rival (2)); the (id,dim) rows in
`s_bd` are rebuilt from them by one sequential transpose pass (rival (3)). */
static u8 s_bpos[(size_t)MAXK * (MAXN + 64)];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
/* s2a_: NARROW THIS RANDOM-ACCESS RECORD. s_bd holds BUCKET INDICES and is clamped at the
write site to [0, Q-1] (`s_bd[...] = bb > Q - 1 ? Q - 1 : bb`, Q <= MAXQ) -- so 9 bits
scans => halving it halves that metadata stream. Accessed scalar-only (no vector
load/store), so no alignment/type-punning hazard. */
static unsigned short s_bd[(size_t)MAXN * MAXK];
static inline u32 folded_rank(u32 i) { u32 v; __builtin_memcpy(&v, s_bd + (size_t)i * MAXK + 6, sizeof(v)); return v; }
static u8 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_tst[MAXN + 2];
static u32 s_nid[MAXN];
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN]; // s_posp[p*MAXK+d] = rank of p in dim d
static u64 s_facc[(size_t)(MAXK - 1) * (((MAXN) + 63) / 64 + 16)];
static u32 s_bnde[(size_t)MAXK * (MAXQ + 2)]; // per-dim coarse-cell boundary positions
static u64 *s_fbm = 0;
static u64 s_fbuf[8192] __attribute__((aligned(64)));
static __m256i khi8v, khi8hv;
static u32 g_Qc = 32;
static u32 g_fbm_W = 0;
static int g_use_filt = 0;
static int g_lexsort = 0;
static int g_nofilt = 0;
static u32 s_fpos[(MAXK + 1) * MAXQ];
static u32 s_bval[(MAXK + 1) * MAXQ];
static u64 *s_bs = 0;
static u64 s_acc[MAXN / 64 + 256] __attribute__((aligned(64)));
static u32 *s_sdmp[MAXK];
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
// ---- [a12x2_] deferred tie resolution, 32 B entries (TIEBATCH-2) ----
// entry = { u64 (i | pos<<20 | d<<40), u32 qT[0..5] } = 32 B
#define TIE_CAP 65536
#define TIE_EW 8
static u32 s_tbuf[(size_t)TIE_CAP * TIE_EW + 4];
static u32 s_tj[TIE_CAP];
static u32 s_tn = 0;
static int g_tieon = 1;
__attribute__((target("avx2,popcnt")))
static void tie_flush(u32 *out) {
u32 n = s_tn;
if (!n) return;
s_tn = 0;
const __m256i all8 = _mm256_set1_epi32(-1);
const __m256i neut = _mm256_set1_epi32(0x7FFFFFFF);
u32 b0 = 0;
while (b0 < n) {
u32 b1 = b0 + 256u; if (b1 > n) b1 = n;
for (u32 q = b0; q < b1; q++) {
const u32 *e = s_tbuf + (size_t)q * TIE_EW;
u64 h = *(const u64 *)e;
u32 j = s_ord[(u32)((h >> 40) & 7u)][(u32)((h >> 20) & 0xFFFFFu)];
s_tj[q] = j;
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
}
for (u32 q = b0; q < b1; q++) {
const u32 *e = s_tbuf + (size_t)q * TIE_EW;
const u32 *cp = s_pm + (size_t)s_tj[q] * STRIDE;
__m256i pv = LDU256(cp);
__m256i qv = _mm256_blend_epi32(_mm256_loadu_si256((const __m256i *)(e + 2)), neut, 0xC0);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qv, pv), all8);
out[(u32)(*(const u64 *)e & 0xFFFFFu)] += ok;
}
b0 = b1;
}
}
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
// ---- byte-signature coarse filter for the fringe scan -------------------
// The fringe scan is BYTE-bound, not instruction-bound: at k=10/n=1.5e5 it moves
// 6.6 GB at 19.9 GB/s (judge DRAM read roofline is 21.1 GB/s), and at k=9/n=2e5
static u8 *s_sig = 0;
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
static u32 g_sigmag = 0;
static u32 s_pofs[(size_t)MAXK * MAXQ];
static u32 s_plen[(size_t)MAXK * MAXQ];
static int g_plen_on = 0;
static u32 g_W, g_Q, g_B;
static u32 g_KS; // row stride of the transposed coordinate rows (64 B for k>8)
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
return (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
}
// Harley-Seal / Muła nibble-popcount. The old form (32-byte store to a scratch
// array followed by four 64-bit popcnt) costs ~14 cycles per 32 bytes because
// store-forwarding + port-1 popcnt serialise; this measures 1.48x faster on a
// pure-L1 microbenchmark (work/c6d2/alucost.cpp: 19.27 -> 13.06 cyc per chunk).
// Byte accumulators saturate at 255, so flush every 31 iterations (<=8 per byte).
__attribute__((target("avx2,popcnt")))
static inline u32 hsum256(__m256i v) {
u64 b[4];
_mm256_storeu_si256((__m256i *)b, v);
return (u32)(b[0] + b[1] + b[2] + b[3]);
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_range(const u64 *const *bs, u32 Wr, u32 rem) {
const __m256i lm = _mm256_set1_epi8(0x0f);
const __m256i lk = _mm256_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);
const __m256i z = _mm256_setzero_si256();
__m256i tot = z, acc0=z, acc1=z, acc2=z, acc3=z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 16 <= Wr; w += 16) {
__m256i a0 = LDU256((bs[0] + w + 0));
for (int d=1; d<KK; d++) a0 = _mm256_and_si256(a0, LDU256((bs[d] + w + 0)));
__m256i lo0 = _mm256_and_si256(a0, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a0, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo0),_mm256_shuffle_epi8(lk,hi0)));
__m256i a1 = LDU256((bs[0] + w + 4));
for (int d=1; d<KK; d++) a1 = _mm256_and_si256(a1, LDU256((bs[d] + w + 4)));
__m256i lo1 = _mm256_and_si256(a1, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(a1, 4), lm);
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo1),_mm256_shuffle_epi8(lk,hi1)));
__m256i a2 = LDU256((bs[0] + w + 8));
for (int d=1; d<KK; d++) a2 = _mm256_and_si256(a2, LDU256((bs[d] + w + 8)));
__m256i lo2 = _mm256_and_si256(a2, lm);
__m256i hi2 = _mm256_and_si256(_mm256_srli_epi16(a2, 4), lm);
acc2 = _mm256_add_epi8(acc2, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo2),_mm256_shuffle_epi8(lk,hi2)));
__m256i a3 = LDU256((bs[0] + w + 12));
for (int d=1; d<KK; d++) a3 = _mm256_and_si256(a3, LDU256((bs[d] + w + 12)));
__m256i lo3 = _mm256_and_si256(a3, lm);
__m256i hi3 = _mm256_and_si256(_mm256_srli_epi16(a3, 4), lm);
acc3 = _mm256_add_epi8(acc3, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo3),_mm256_shuffle_epi8(lk,hi3)));
if (++cnt == 31) {
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc2,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc3,z));
acc0=acc1=acc2=acc3=z; cnt=0;
}
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc2,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc3,z));
total = hsum256(tot);
for (; w + 4 <= Wr; w += 4) {
__m256i a=LDU256(bs[0]+w);
for (int d=1;d<KK;d++) a=_mm256_and_si256(a,LDU256(bs[d]+w));
__m256i lo=_mm256_and_si256(a,lm), hi=_mm256_and_si256(_mm256_srli_epi16(a,4),lm);
__m256i pc=_mm256_add_epi8(_mm256_shuffle_epi8(lk,lo),_mm256_shuffle_epi8(lk,hi));
total += hsum256(_mm256_sad_epu8(pc,z));
}
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 1; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v);
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 1; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & ((1ull << rem) - 1));
}
return total;
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const u32 KS_C = (u32)KK; // [n12ac] KSC: g_KS is a pure function of KK
s_tn = 0; // [a12x2_] TIEBATCH-2: no stale entries across calls
const int k = KK;
const int K1 = KK - 1;
#if MF_ASCR
// ---- pass A locality: order the group by the bucket in one non-fold dim so that the
// ~Q-run of queries needing array (e1,b) is processed together, and pin that
// array's prefix in a small scratch (the other k-2 dims still stream from DRAM
// and evict a 70 KB hot array within one query, which is why merely sorting
// the queries was measured a wash). ----
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB] * MAXK + d1), _MM_HINT_T0);
#endif
cntv[s_bd[(size_t)s_qord[qi] * MAXK + d1] + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
if (qi + 96 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + 96] * MAXK + d1), _MM_HINT_T0);
u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * MAXK + d1] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#else
const u32 d1 = 0, e1 = 0;
#endif
// ---- pass A: aligned AND; queries ordered lexicographically by their bucket vector so
// that consecutive queries share bitset arrays (keeps the working set small). ----
if (g_lexsort) {
u32 *cntv = s_cntv;
for (int dd = k - 1; dd >= 0; dd--) {
if (dd == ds) continue;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) cntv[s_bd[(size_t)s_qord[qi] * MAXK + dd] + 1]++;
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) { u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * MAXK + dd] + 1]++] = id; }
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = s_qord[qi + PF_DIST];
const unsigned short *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
#if MF_ASCR
{
u32 qi = q0;
while (qi < q1) {
u32 b1 = s_bd[(size_t)s_qord[qi] * MAXK + d1];
u32 qe = qi + 1;
while (qe < q1 && s_bd[(size_t)s_qord[qe] * MAXK + d1] == b1) qe++;
const u64 *bp1 = g_plen_on ? (s_bs + s_pofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
u32 mx = 1;
for (u32 t = qi; t < qe; t++) { u32 rr = folded_rank(s_qord[t]); if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 1;
{ /* [k2t] bp1 pin enabled */
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = s_qord[qq];
const unsigned short *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
}
}
#endif
if (g_plen_on && qq + 1 < qe) {
const u32 jpf=s_qord[qq+1];
const unsigned short *bdpf=s_bd+(size_t)jpf*MAXK;
for(int epf=0;epf<K1;epf++) {
if(epf==(int)e1) continue;
const int dpf=(epf<ds)?epf:epf+1;
const u64 *ppf=s_bs+s_pofs[(size_t)epf*MAXQ+bdpf[dpf]];
__builtin_prefetch(ppf,0,3);
}
}
const u32 R = folded_rank(i);
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
if (e == (int)e1) bp[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
const unsigned short *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = s_qord[qi + PF_DIST];
const unsigned short *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
const u32 R = folded_rank(i);
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W);
#ifdef DBG_PLEN
if (g_plen_on) {
u32 need = (R >> 6) + ((R & 63) ? 1 : 0);
if (need > s_plen[(size_t)e * MAXQ + bd[d]])
printf("DBG ds=%d e=%d d=%d b=%u R=%u need=%u len=%u\n", ds, e, d, bd[d], R, need, s_plen[(size_t)e * MAXQ + bd[d]]);
}
#endif
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
#endif
if (g_skip_fringe) return;
// ---- pass B: fringe, dimension-major. Ordering the group's queries by their bucket in
// dim d makes every scan walk one L1-resident bucket block of the dim-d order. ----
for (int d = 0; d < k; d++) {
if (d == ds) continue;
// ARM SIGC: the branchy `if (e == (int)d) continue;` skip (SGNL) removed the e==d
// lane's work but paid a per-iteration branch inside the hottest loop of pass B.
// Build the lane list ONCE PER d instead and iterate the compacted list, so the
// wasted plane load disappears with NO branch in the loop body.
int kk1 = 0; int emap[16]; const u8 *sgb2[16];
for (int e = 0; e < k; e++) if (e != (int)d) emap[kk1++] = e;
{ const u8 *sbase = s_sig + (size_t)d * k * ((size_t)N + 64);
const size_t spl = (size_t)N + 64;
for (int e2 = 0; e2 < kk1; e2++) sgb2[e2] = sbase + (size_t)emap[e2] * spl; }
#if !MF_FR_NOSORT
#if MF_ASCR
/* [NODUP] (BRIEF 2.18.934 / 2.18.938): the ASCR prefix above already ran a STABLE
counting sort of s_qord[q0,q1) on the key s_bd[id*MAXK + d1], and d1 = (ds==0)?1:0
is exactly the FIRST non-ds dimension this loop reaches. Re-sorting by the same
key over a (byte-)unchanged input reproduces that order exactly => BIT-IDENTICAL,
so the whole block is skipped.
-- Output audit required by 2.18.938 (the 1015 lesson: enumerate EVERY symbol the
skipped pass writes, then find who else can write it correctly) --
The skipped block writes exactly three symbols:
s_qord : already correct -- written by the ASCR prefix (stable, same key) and by
NOTHING between that prefix and this loop.
s_cntv : pure scratch -- re-zeroed at the top of every sort; final value unconsumed.
s_tmp : pure scratch -- written and immediately copied back out inside the same
block; final value unconsumed.
This engine has NO s_bks / s_bkt bucket-number cache (grep = 0 hits), so there is no
hidden second consumer of the skipped pass's work. */
if ((u32)d != d1)
#endif
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB2] * MAXK + d), _MM_HINT_T0);
#endif
cntv[s_bd[(size_t)s_qord[qi] * MAXK + d] + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
u32 id = s_qord[qi];
#if MF_PFB > 0
if (qi + MF_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB2] * MAXK + d), _MM_HINT_T0);
#endif
s_tmp[cntv[s_bd[(size_t)id * MAXK + d] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#endif
u32 exmask = 0;
for (int e = 0; e < d; e++) if (e != ds) exmask |= 1u << e;
// "don't care" lanes are forced to a maximum-unsigned value, so one unsigned max test
// per candidate decides domination *and* the exclusion rule together.
u32 hiA[8], hiAh[8];
for (int e = 0; e < 8; e++) hiA[e] = ((exmask >> e) & 1u) ? 0u : 0x7FFFFFFFu;
for (int e = 8; e < 16; e++) hiAh[e - 8] = ((exmask >> e) & 1u) ? 0u : 0x7FFFFFFFu;
const __m256i notexlo = LDU256(hiA);
const __m256i notexloh = LDU256(hiAh);
const __m256i all8 = _mm256_set1_epi32(-1);
const u32 *const odd = s_ord[d];
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
#if MF_PFB3 > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_T0);
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK + d), _MM_HINT_T0);
__asm__ __volatile__("prefetchw %0" :: "m"(*(char *)(out + j))); /* [y12a_pfo] RMW target */
}
}
#endif
const unsigned short *bd = s_bd + (size_t)i * MAXK;
const u32 lo = FPOS(d, bd[d]);
const u32 hi = s_gt[d][i];
if (lo >= hi) continue;
const u32 *rec = s_pm + (size_t)i * STRIDE;
/* [AVREG] The threshold vector used to be assembled by storing 8 scalars into a stack
array `Av[8]` and then reloading all 32 bytes with one vector load. That load can NOT
be store-forwarded (a 32-byte load is not contained in any single 4-byte store), so the
consumer stalls until all 8 stores retire to L1 -- on the hottest per-query path in
pass B (1.5e6 executions). Build the lanes in registers instead (vpinsrd chain, no
memory round trip). Lanes >= k (and the fold lane ds) are forced to 0x7FFFFFFF, i.e.
exactly what the old code stored. */
__m256i qr = _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
__m256i Avm, qTh;
if (k <= 8) {
const u32 a0 = (0 != (int)ds) ? BVAL(0, bd[0]) : 0x7FFFFFFFu;
const u32 a1 = (1 != (int)ds) ? BVAL(1, bd[1]) : 0x7FFFFFFFu;
const u32 a2 = (2 != (int)ds) ? BVAL(2, bd[2]) : 0x7FFFFFFFu;
const u32 a3 = (3 != (int)ds) ? BVAL(3, bd[3]) : 0x7FFFFFFFu;
const u32 a4 = (k > 4 && 4 != (int)ds) ? BVAL(4, bd[4]) : 0x7FFFFFFFu;
const u32 a5 = (k > 5 && 5 != (int)ds) ? BVAL(5, bd[5]) : 0x7FFFFFFFu;
const u32 a6 = (k > 6 && 6 != (int)ds) ? BVAL(6, bd[6]) : 0x7FFFFFFFu;
const u32 a7 = (k > 7 && 7 != (int)ds) ? BVAL(7, bd[7]) : 0x7FFFFFFFu;
Avm = _mm256_or_si256(_mm256_setr_epi32((int)a0, (int)a1, (int)a2, (int)a3,
(int)a4, (int)a5, (int)a6, (int)a7), notexlo);
qTh = all8;
} else {
u32 Av[STRIDE] __attribute__((aligned(32)));
for (int e = 0; e < 8; e++) Av[e] = (e != (int)ds) ? BVAL(e, bd[e]) : 0x7FFFFFFFu;
Avm = _mm256_or_si256(LDU256(Av), notexlo);
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU256((Av + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
}
__m256i qT = _mm256_min_epi32(qr, Avm);
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (size_t)lo * KS_C) : 0;
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const u32 FW = g_fbm_W;
int first = 1;
#if MF_FR_NOREAD
first = 0;
for (u32 w = 0; w < nw; w++) fb[w] = 0;
#endif
for (int e = 0; e < k; e++) {
if (e == d) continue;
u32 c = (u32)((u64)rec[e] / wcs) + 1;
if (c > g_Qc) c = g_Qc;
const u64 *row = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
if (first) { for (u32 w = 0; w < nw; w++) fb[w] = row[w]; first = 0; }
else { for (u32 w = 0; w < nw; w++) fb[w] &= row[w]; }
}
fb[0] &= ~0ull << (lo & 63);
u32 r2 = hi & 63u;
if (r2) fb[nw - 1] &= (1ull << r2) - 1;
for (u32 w = 0; w < nw; w++) {
u64 vv = fb[w];
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = ((w0 + w) << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * KS_C) : (s_pm + (size_t)s_ord[d][j] * STRIDE);
__m256i pv = LDU256(cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (g_sig_on && (hi - lo) >= 32) {
// ---- coarse byte-signature scan, exact fallback only on ties ----
// (no qtv store: the broadcast signatures are built in the vector domain)
__m256i qsb[16];
// Lane d needs NO test at all: the band is [FPOS(d,bd_d), gt_d(i)) in dim-d
// order, and gt_d(i) is the START of i's tie group, so every candidate in the
// band has x_d[c] < x_d[i] = qT[d] by construction. Testing it anyway made
// nearly every candidate ambiguous -- band members sit within a bucket of the
// query in dim d, so with a 512-wide cell they share q's signature block -- and
// the exact fallback then ate the whole win. Neutralising the lane is EXACT.
{
const __m256i zero8 = _mm256_setzero_si256();
const __m256i magv = _mm256_set1_epi32((int)g_sigmag);
const __m256i loe = _mm256_mul_epu32(qT, magv);
const __m256i loo = _mm256_mul_epu32(_mm256_srli_epi64(qT, 32), magv);
// [k2w_/HIFUSE] qsb[e2] carries sig(qT_e) - 1 (saturating): one
// `subs_epu8(cv, qsb[e2])` then reads 0 / 1 / >=2 for strictly-less /
// tie / strictly-greater. sig is byte 0 of each 32-bit lane and the
// broadcast replicates exactly that byte, so a saturating BYTE
// decrement decrements precisely it.
const __m256i one8 = _mm256_set1_epi8(1);
const __m256i svlo = _mm256_subs_epu8(_mm256_blend_epi32(_mm256_srli_epi64(loe, 32),
_mm256_slli_epi64(_mm256_srli_epi64(loo, 32), 32), 0xAA), one8);
__m256i svhi = svlo;
if (k > 8) {
const __m256i hie = _mm256_mul_epu32(qTh, magv);
const __m256i hio = _mm256_mul_epu32(_mm256_srli_epi64(qTh, 32), magv);
svhi = _mm256_subs_epu8(_mm256_blend_epi32(_mm256_srli_epi64(hie, 32),
_mm256_slli_epi64(_mm256_srli_epi64(hio, 32), 32), 0xAA), one8);
}
for (int e2 = 0; e2 < kk1; e2++) {
int e = emap[e2];
__m256i sv = (e < 8) ? svlo : svhi;
__m256i bc = _mm256_permutevar8x32_epi32(sv, _mm256_set1_epi32(e & 7));
qsb[e2] = _mm256_shuffle_epi8(bc, zero8);
}
}
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
/* [k2w_/HIFUSE] ONE accumulator. Each lane's byte is 0 / 1 / >=2 for
strictly-less / tie / strictly-greater, and an OR cannot make one of
those look like another, so acc alone carries all three states:
acc & 0xFE == 0 <=> the old `ltm` (every lane <=)
acc == 1 <=> the old `amb` (tie -> exact row test)
The tie mask is read out of acc itself, so no plane is re-loaded. */
__m256i acc0 = _mm256_setzero_si256(), acc1 = _mm256_setzero_si256();
for (int e2 = 0; e2 < kk1; e2++) {
__m256i cv = LDU256((sgb2[e2] + aj));
if (e2 & 1) acc1 = _mm256_or_si256(acc1, _mm256_subs_epu8(cv, qsb[e2]));
else acc0 = _mm256_or_si256(acc0, _mm256_subs_epu8(cv, qsb[e2]));
}
__m256i acc = _mm256_or_si256(acc0, acc1);
u32 decm = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_setzero_si256()));
u32 amb = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_set1_epi8(1)));
aacc += (u32)__builtin_popcount(decm);
if (amb && k == 6 && !c0 && g_tieon && s_tn + 32u <= (u32)TIE_CAP) {
u32 *te = s_tbuf + (size_t)s_tn * TIE_EW;
do {
int t = __builtin_ctz(amb); amb &= amb - 1;
*(u64 *)te = (u64)i | ((u64)(aj + (u32)t) << 20) | ((u64)(u32)d << 40);
_mm_storeu_si128((__m128i *)(te + 2), _mm256_castsi256_si128(qT));
_mm_storel_epi64((__m128i *)(te + 6), _mm256_extracti128_si256(qT, 1));
te += TIE_EW; s_tn++;
} while (amb);
} else 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 acc0 = _mm256_setzero_si256(), acc1 = _mm256_setzero_si256();
for (int e2 = 0; e2 < kk1; e2++) {
__m256i cv = LDU256((sgb2[e2] + base));
if (e2 & 1) acc1 = _mm256_or_si256(acc1, _mm256_subs_epu8(cv, qsb[e2]));
else acc0 = _mm256_or_si256(acc0, _mm256_subs_epu8(cv, qsb[e2]));
}
__m256i acc = _mm256_or_si256(acc0, acc1);
u32 decm = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_setzero_si256()))) & kmask;
u32 amb = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(decm);
if (amb && k == 6 && !c0 && g_tieon && s_tn + 32u <= (u32)TIE_CAP) {
u32 *te = s_tbuf + (size_t)s_tn * TIE_EW;
do {
int t = __builtin_ctz(amb); amb &= amb - 1;
*(u64 *)te = (u64)i | ((u64)(base + (u32)t) << 20) | ((u64)(u32)d << 40);
_mm_storeu_si128((__m128i *)(te + 2), _mm256_castsi256_si128(qT));
_mm_storel_epi64((__m128i *)(te + 6), _mm256_extracti128_si256(qT, 1));
te += TIE_EW; s_tn++;
} while (amb);
} else 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 = _mm256_load_si256((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = _mm256_load_si256((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
if (k == 6) tie_flush(out);
}
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int k) {
// W padded to 8 words: every (dim,bucket) array then starts on a 64-byte
// boundary and every 32-byte vector load inside a stream is line-contained.
const u32 W = ((N + 63) / 64 + 7u) & ~7u;
g_W = W;
#ifdef MF_SKIP_AND
g_skip_and = 1;
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_Q
g_Q = MF_Q; g_B = (u32)((N + (MF_Q) - 1) / (MF_Q));
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_SKIP_FRINGE
g_skip_fringe = 1;
#endif
if (!g_Q) {
// bucket size ~ 0.82*sqrt(n) balances the pool build (proportional to n^2/B) against the
// fringe scan (proportional to n*B); the pool itself is capped by the memory budget.
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++; // ceil(sqrt(n)) without libm
// Reserve: the .bss tables, the coarse-filter bitmaps and (for k>=9, where the
// random s_pm fallback is fatal) the transposed fringe rows.
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
// The byte-signature planes are k*k*(N+64) bytes and were NOT charged here, so at
// n=1e6/k=10 the pool was sized ~100 MB too large, the allocation failed, and the
// engine returned WRONG ANSWERS on the board (1016-SIG8, sid 98573, WA at 8.929 s)
// instead of erroring -- base8, which has no planes, was AC on the same row. Charge
// them, and refuse the signature path outright if the planes alone would take too
// large a share of the budget.
overhead += (u64)k * k * ((u64)N + 64);
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
// MF_PLEN sizes each pooled array to the largest R ever read from it, so the
// pool really used is only ~fill of (k-1)*(Q+1)*W words; charge that instead.
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46}; // per-array fill, k=4..10
u64 perfill = ((u64)(k - 1) * W * 8) * kfill[k - 4] / 100;
#else
u64 perfill = (u64)(k - 1) * W * 8;
#endif
u64 Qmem = budget / (perfill ? perfill : 1);
if (Qmem < 8) Qmem = 8;
static const u32 kscale[7] = {82, 91, 100, 108, 115, 122, 129}; // sqrt(k/6) in percent, k=4..10
u32 q = (u32)((s2 * MF_QCOEF + 99) / 100);
q = (u32)(((u64)q * kscale[k - 4] + 99) / 100);
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
g_Q = q;
g_B = (N + q - 1) / q;
}
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
g_plen_on = 1; // independent of MF_WTRIM: lengths only, stride stays W
#endif
#if MF_WTRIM
u32 Wpool = W; // MF_WTRIM is KNOWN-BROKEN, see the warning above
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1; // coarse grid: the filter beats the scan
#endif
if (g_use_filt && (u64)k * k * (g_Qc + 1) * W * 8 > ((u64)700 << 20)) g_use_filt = 0;
#if !MF_FILTAUTO
g_use_filt = 0;
#endif
#if MF_NOFILTB
// The coarse filter is disabled for this variant. Apply that decision
// before the order pass, so inverse ranks for the unused filter are skipped.
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
// ---------- per-dimension orders, ranks, tie info ----------
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
for (u32 i = 0; i < N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i < N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tst[i] = s; s += c; } s_tst[N] = s; }
u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
u32 pos = s_cntv[v]++;
od[pos] = i;
gt[i] = s_tst[v];
u32 bb = (s_tst[v + 1] - 1u) / B;
s_bpos[(size_t)d * ((size_t)N + 64) + i] = (u8)(bb > Q - 1 ? Q - 1 : bb);
}
#if MF_FB2
{ for (u32 j = 0; j < N; j++) { u32 p = od[j]; s_posp[(size_t)p * MAXK + d] = j; } }
#else
// s_pos is read ONLY by the non-FB3 coarse-filter build, which is compiled out at
// the shipped MF_FB3=1 (the FB3 path uses s_bnde, and FB2 uses s_posp). So this
// fill is k scattered n-element writes per run with no live consumer.
// ⚠ RESTORED. I had removed this as "dead work" and it is NOT dead: the reader is
// `const u32 *pd = s_pos[d];` inside the `#else` of `#if MF_FB2` -- NOT an MF_FB3
// guard -- so it IS compiled at the shipped MF_FB2=0. It is reached only when the
// coarse filter is built, and MF_FILTAUTO enables that only for B >= 400. Removing
// the fill therefore produced CORRECT answers for every sweep I ran (my small-n
// sweeps all had B = n/Q < 400, so the filter was off and s_pos unread) and WRONG
// answers whenever B >= 400 -- which is exactly the board pattern: 1012b/1013b/1014b
// (B=259) AC while 1013a (B=416) and 1016 (B=842) were WA.
// Only the coarse-filter build reads s_pos, and MF_FILTAUTO turns that on only for
// B >= 400. Guarding on the ACTUAL flag (rather than on a #if) is both correct and
// free: at B < 400 the fill is genuinely dead and is skipped, so small-B rows keep
// the speed the (buggy) ungated build had, while large-B rows are now correct.
if (g_use_filt) { u32 *pp = s_pos[d]; for (u32 j = 0; j < N; j++) pp[od[j]] = j; }
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
// Transpose the bucket columns in one sequential point-record pass.
for (u32 i = 0; i < N; i++) {
unsigned short *row = s_bd + (size_t)i * MAXK;
for (int d = 0; d < k; d++)
row[d] = s_bpos[(size_t)d * ((size_t)N + 64) + i];
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < STRIDE; d++) p[d] = 0;
}
// ---------- byte-signature planes for the fringe coarse filter ----------
// s_sig[d][e*N + j] = (u8)(x[e][ od_d[j] ] >> SIGSH) (od_d = dim-d order)
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; } // smallest sh with (N-1)>>sh <= 255
g_sigsh = sh;
{ u32 dn = (N > 1) ? (N - 1) : 1u;
g_sigmag = (u32)((((u64)255 << 32) / (u64)dn) + 1u); }
size_t need = (size_t)k * k * ((size_t)N + 64);
static u8 *sigpool = 0; static size_t sigcap = 0;
// Refuse the signature path when the planes alone exceed 220 MB: past that the pool
// it competes with is squeezed and the row is better served by the exact path.
if (need > ((size_t)220 << 20)) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
// The plane build is DEFERRED and merged into the s_sdmp build pass below when
// s_sdmp is allocated: both walk the same `s_pm + od[j]*STRIDE` rows, so two
// passes pay two rounds of random 64-byte reads for one round of useful work.
// ⚠ SAFETY GATE. The signature path is VERIFIED CORRECT (brute force k=4..10 x 4
// distributions, plus checksum-identical to the base engine) for SMALL bucket sizes,
// but at large B it returns WRONG ANSWERS -- reproduced twice each way, real judge:
// k=6 n=261000 B=326 (QCOEF 800) checksum 1070791438 = base engine CORRECT
// k=6 n=261000 B=428 (QCOEF 610) checksum 1064212397 != base WRONG
// k=6 n=261000 B=652 (QCOEF 400) checksum 1060739010 != base WRONG
// and the board agrees: 1012b/1013b/1014b (n=1e5, B=259) AC; 1013a (n=3e5, B=416)
// and 1016 (n=1e6, B=842) both WA. The threshold sits between B=326 and B=428.
// Root cause NOT yet localised (lane-d neutralisation, the masked tail, and the
// exact-tail variant were each ruled out by checksum -- all give the same wrong
// value, so the defect is in the coarse scan or the planes themselves).
// UNTIL IT IS FOUND, REFUSE THE SIGNATURE PATH ABOVE A MEASURED-SAFE B.
g_sig_on = 1;
g_sig_merged = 0;
}
}
// ---------- coarse superset bitmaps over every dim order (fringe filter) ----------
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
// Sweep the e-order ONCE per dim e, keeping (k-1) running bitsets (one per
// target dim d), and copy them out at each coarse-cell boundary. The old
// shape re-walked all N points for every one of the k(k-1) (d,e) pairs with
// two L3-random reads each (xe[pe[ptr]] and s_pos[d][pe[ptr]]); here the
// rank record is point-major so one 64B line yields every dim's rank, which
// cuts the random reads by ~k and the whole build by ~2-2.5x at k>=8.
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) {
__asm__ volatile("nop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\t");
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w))); }
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
// ---------- fold dimension + bucket vector per query ----------
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
__builtin_memcpy(s_bd + (size_t)i * MAXK + 6, &best, sizeof(best));
// Largest R that will ever be read out of each (fold dim, dim, bucket) array:
// the pooled array (ds,d,b) is only ever scanned over [0, maxR] words, so the
// per-array length can be far below the group maximum (0.38n vs 0.93n at k=10).
if (g_plen_on) {
const unsigned short *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
// group by ds
{
for (u32 v = 0; v <= (u32)k; v++) s_cntv[v] = 0;
for (u32 i = 0; i < N; i++) s_cntv[s_g[i] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)k; v++) { u32 c = s_cntv[v]; s_cntv[v] = ac; ac += c; } }
for (u32 i = 0; i < N; i++) { u32 id = s_qord[i]; s_tmp[s_cntv[s_g[id] + 1]++] = id; }
for (u32 i = 0; i < N; i++) s_qord[i] = s_tmp[i];
}
u32 gstart[MAXK + 2];
{ u32 c = 0; for (int d = 0; d <= k; d++) { gstart[d] = c; while (c < N && (int)s_g[s_qord[c]] == d) c++; } gstart[k + 1] = c; }
// ---------- transposed coordinate rows (sequential fringe) ----------
#if MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 8; // 32-byte aligned candidate rows for k<=8
#endif
size_t per = (size_t)N * KS * 4;
size_t lim = (size_t)150 << 20;
if (k >= 9) lim = (size_t)470 << 20; // k=9/10 need the rows: the random s_pm fallback is fatal
if (per * k <= lim) {
static u32 *pool = 0; static size_t cap = 0;
if (per * k > cap) { pool = (u32 *)pool_alloc(per * k + 256); cap = per * k; }
if (pool) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * k * ((size_t)N + 64);
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); sb[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
// ---------- signature planes, standalone only if the s_sdmp pass did not run ----
if (s_sig && !g_sig_merged && k == 6) {
static u8 cached_signature[(size_t)MAXN * 8] __attribute__((aligned(64)));
for (u32 i = 0; i < N; i++)
for (int e = 0; e < 6; e++) {
u32 v = (u32)(((u64)x[e][i] * g_sigmag) >> 32);
cached_signature[(size_t)i*8 + e] = (u8)(v ? v : 1u);
}
for (int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*k*((size_t)N + 64);
u32 j=0;
for (;j+8<=N;j+=8) {
__m128i r0=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+0]*8));
__m128i r1=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+1]*8));
__m128i r2=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+2]*8));
__m128i r3=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+3]*8));
__m128i r4=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+4]*8));
__m128i r5=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+5]*8));
__m128i r6=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+6]*8));
__m128i r7=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+7]*8));
__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);
_mm_storel_epi64((__m128i *)(base+(size_t)0*((size_t)N + 64)+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*((size_t)N + 64)+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*((size_t)N + 64)+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*((size_t)N + 64)+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*((size_t)N + 64)+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*((size_t)N + 64)+j),_mm_srli_si128(v2,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*8;
for(int e=0;e<k;e++)base[(size_t)e*((size_t)N + 64)+j]=row[e];
}
}
g_sig_merged=1;
}
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * k * ((size_t)N + 64);
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); base[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
}
// ---------- pool ----------
u64 poolwords = (u64)K1 * (Q + 1) * W;
#if MF_PLEN
if (g_plen_on) {
u64 best = 0;
for (int dsx = 0; dsx < k; dsx++) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < dsx) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)dsx * MAXK + d) * MAXQ;
for (u32 b = 0; b <= Q; b++) {
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > W) len = W;
tot += len;
}
}
if (tot > best) best = tot;
}
if (best && best < poolwords) poolwords = best;
}
#endif
size_t need = (size_t)poolwords * sizeof(u64);
static u64 *pool = 0; static size_t cap2 = 0;
if (need > cap2) { pool = (u64 *)pool_alloc(need + 4096); cap2 = need; }
s_bs = pool;
// ---------- groups ----------
for (int ds = 0; ds < k; ds++) {
u32 q0 = gstart[ds], q1 = gstart[ds + 1];
if (q0 >= q1) continue;
#if MF_WTRIM
{
u32 mx = 1;
for (u32 qi = q0; qi < q1; qi++) { u32 rr = folded_rank(s_qord[qi]); if (rr > mx) mx = rr; }
u32 wp = ((mx >> 6) + 8u) & ~7u;
Wpool = (wp < W) ? wp : W;
}
#endif
#if MF_PLEN
// Per-(dim,bucket) array lengths inside the group.
if (g_plen_on) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)ds * MAXK + d) * MAXQ;
u32 *po = s_pofs + (size_t)e * MAXQ;
u32 *pl = s_plen + (size_t)e * MAXQ;
for (u32 b = 0; b <= Q; b++) {
// NOTE: bucket 0 is the EMPTY prefix set (all-zero array), but a query whose
// bd[d]==0 still scans [0,R) words of it, so the array must be long enough --
// sizing it to the fixed 8-word minimum reads past the zeroed region into the
// next array. Use the real per-bucket max R here, not 0.
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > Wpool) len = Wpool;
po[b] = (u32)tot; pl[b] = len; tot += len;
}
}
if (tot > (u64)K1 * (Q + 1) * W) { g_plen_on = 0; } // never; safety
}
#endif
// bit numbering = rank in dim ds
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
u64 *base = g_plen_on ? s_bs : (s_bs + (size_t)e * (Q + 1) * Wpool);
u64 *cur = s_acc;
for (u32 w = 0; w < Wpool; w++) cur[w] = 0;
if (g_plen_on) {
// bucket 0 is the empty prefix set; bd[] can be 0 so it must be materialised as zeros
u64 *d0 = s_bs + s_pofs[(size_t)e * MAXQ];
u32 l0 = s_plen[(size_t)e * MAXQ];
for (u32 w = 0; w < l0; w++) d0[w] = 0;
}
u32 ptr = 0;
#if MF_BUILD_NOPOOL
while (ptr < N) ptr++;
#else
for (u32 b = 1; b <= Q; b++) {
// {q : x_d[q] < BVAL(d,b)} is EXACTLY the prefix [0,FPOS(d,b)) of the dim-d order
// (FPOS is the tie-group start containing rank b*B), so the per-point value compare
// -- one random xd read per point per (ds,d) pair, k(k-1)N of them -- is redundant.
// b==Q means "every point" and FPOS(d,Q) is the 0 sentinel, hence N.
// Credited to agent a7f799ceb7e63b0cd (worth 4-5% at k=4/5 on the board).
u32 pend = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) {
if (ptr + 64 < MAXN) __builtin_prefetch(s_nid + od[ptr + 64], 0, 3); // [PF-NID]
u32 q = s_nid[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++;
}
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
#if MF_NT
// The pool is written once and read 89 GB later, so a plain copy pays a
// read-for-ownership DRAM read per line for nothing. NT stores remove it
// (measured: this loop was 4.58 GB / ~1.7 s, i.e. 2.7 GB/s, vs 12.5 GB/s
// for a plain memset).
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
template<int K> static void solve_mf_sp(u32 N, const unsigned **x, u32 *out) {
const int k = K;
// W padded to 8 words: every (dim,bucket) array then starts on a 64-byte
// boundary and every 32-byte vector load inside a stream is line-contained.
const u32 W = ((N + 63) / 64 + 7u) & ~7u;
g_W = W;
#ifdef MF_SKIP_AND
g_skip_and = 1;
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_Q
g_Q = MF_Q; g_B = (u32)((N + (MF_Q) - 1) / (MF_Q));
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_SKIP_FRINGE
g_skip_fringe = 1;
#endif
if (!g_Q) {
// bucket size ~ 0.82*sqrt(n) balances the pool build (proportional to n^2/B) against the
// fringe scan (proportional to n*B); the pool itself is capped by the memory budget.
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++; // ceil(sqrt(n)) without libm
// Reserve: the .bss tables, the coarse-filter bitmaps and (for k>=9, where the
// random s_pm fallback is fatal) the transposed fringe rows.
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
// The byte-signature planes are k*k*(N+64) bytes and were NOT charged here, so at
// n=1e6/k=10 the pool was sized ~100 MB too large, the allocation failed, and the
// engine returned WRONG ANSWERS on the board (1016-SIG8, sid 98573, WA at 8.929 s)
// instead of erroring -- base8, which has no planes, was AC on the same row. Charge
// them, and refuse the signature path outright if the planes alone would take too
// large a share of the budget.
overhead += (u64)k * k * ((u64)N + 64);
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
// MF_PLEN sizes each pooled array to the largest R ever read from it, so the
// pool really used is only ~fill of (k-1)*(Q+1)*W words; charge that instead.
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46}; // per-array fill, k=4..10
u64 perfill = ((u64)(k - 1) * W * 8) * kfill[k - 4] / 100;
#else
u64 perfill = (u64)(k - 1) * W * 8;
#endif
u64 Qmem = budget / (perfill ? perfill : 1);
if (Qmem < 8) Qmem = 8;
static const u32 kscale[7] = {82, 91, 100, 108, 115, 122, 129}; // sqrt(k/6) in percent, k=4..10
u32 q = (u32)((s2 * MF_QCOEF + 99) / 100);
q = (u32)(((u64)q * kscale[k - 4] + 99) / 100);
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
g_Q = q;
g_B = (N + q - 1) / q;
}
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
g_plen_on = 1; // independent of MF_WTRIM: lengths only, stride stays W
#endif
#if MF_WTRIM
u32 Wpool = W; // MF_WTRIM is KNOWN-BROKEN, see the warning above
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1; // coarse grid: the filter beats the scan
#endif
if (g_use_filt && (u64)k * k * (g_Qc + 1) * W * 8 > ((u64)700 << 20)) g_use_filt = 0;
#if !MF_FILTAUTO
g_use_filt = 0;
#endif
#if MF_NOFILTB
// The coarse filter is disabled for this variant. Apply that decision
// before the order pass, so inverse ranks for the unused filter are skipped.
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
// ---------- per-dimension orders, ranks, tie info ----------
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
for (u32 i = 0; i < N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i < N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tst[i] = s; s += c; } s_tst[N] = s; }
u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
u32 pos = s_cntv[v]++;
od[pos] = i;
gt[i] = s_tst[v];
u32 bb = (s_tst[v + 1] - 1u) / B;
s_bpos[(size_t)d * ((size_t)N + 64) + i] = (u8)(bb > Q - 1 ? Q - 1 : bb);
}
#if MF_FB2
{ for (u32 j = 0; j < N; j++) { u32 p = od[j]; s_posp[(size_t)p * MAXK + d] = j; } }
#else
// s_pos is read ONLY by the non-FB3 coarse-filter build, which is compiled out at
// the shipped MF_FB3=1 (the FB3 path uses s_bnde, and FB2 uses s_posp). So this
// fill is k scattered n-element writes per run with no live consumer.
// ⚠ RESTORED. I had removed this as "dead work" and it is NOT dead: the reader is
// `const u32 *pd = s_pos[d];` inside the `#else` of `#if MF_FB2` -- NOT an MF_FB3
// guard -- so it IS compiled at the shipped MF_FB2=0. It is reached only when the
// coarse filter is built, and MF_FILTAUTO enables that only for B >= 400. Removing
// the fill therefore produced CORRECT answers for every sweep I ran (my small-n
// sweeps all had B = n/Q < 400, so the filter was off and s_pos unread) and WRONG
// answers whenever B >= 400 -- which is exactly the board pattern: 1012b/1013b/1014b
// (B=259) AC while 1013a (B=416) and 1016 (B=842) were WA.
// Only the coarse-filter build reads s_pos, and MF_FILTAUTO turns that on only for
// B >= 400. Guarding on the ACTUAL flag (rather than on a #if) is both correct and
// free: at B < 400 the fill is genuinely dead and is skipped, so small-B rows keep
// the speed the (buggy) ungated build had, while large-B rows are now correct.
if (g_use_filt) { u32 *pp = s_pos[d]; for (u32 j = 0; j < N; j++) pp[od[j]] = j; }
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
// Transpose the bucket columns in one sequential point-record pass.
for (u32 i = 0; i < N; i++) {
unsigned short *row = s_bd + (size_t)i * MAXK;
for (int d = 0; d < k; d++)
row[d] = s_bpos[(size_t)d * ((size_t)N + 64) + i];
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < STRIDE; d++) p[d] = 0;
}
// ---------- byte-signature planes for the fringe coarse filter ----------
// s_sig[d][e*N + j] = (u8)(x[e][ od_d[j] ] >> SIGSH) (od_d = dim-d order)
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; } // smallest sh with (N-1)>>sh <= 255
g_sigsh = sh;
{ u32 dn = (N > 1) ? (N - 1) : 1u;
g_sigmag = (u32)((((u64)255 << 32) / (u64)dn) + 1u); }
size_t need = (size_t)k * k * ((size_t)N + 64);
static u8 *sigpool = 0; static size_t sigcap = 0;
// Refuse the signature path when the planes alone exceed 220 MB: past that the pool
// it competes with is squeezed and the row is better served by the exact path.
if (need > ((size_t)220 << 20)) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
// The plane build is DEFERRED and merged into the s_sdmp build pass below when
// s_sdmp is allocated: both walk the same `s_pm + od[j]*STRIDE` rows, so two
// passes pay two rounds of random 64-byte reads for one round of useful work.
// ⚠ SAFETY GATE. The signature path is VERIFIED CORRECT (brute force k=4..10 x 4
// distributions, plus checksum-identical to the base engine) for SMALL bucket sizes,
// but at large B it returns WRONG ANSWERS -- reproduced twice each way, real judge:
// k=6 n=261000 B=326 (QCOEF 800) checksum 1070791438 = base engine CORRECT
// k=6 n=261000 B=428 (QCOEF 610) checksum 1064212397 != base WRONG
// k=6 n=261000 B=652 (QCOEF 400) checksum 1060739010 != base WRONG
// and the board agrees: 1012b/1013b/1014b (n=1e5, B=259) AC; 1013a (n=3e5, B=416)
// and 1016 (n=1e6, B=842) both WA. The threshold sits between B=326 and B=428.
// Root cause NOT yet localised (lane-d neutralisation, the masked tail, and the
// exact-tail variant were each ruled out by checksum -- all give the same wrong
// value, so the defect is in the coarse scan or the planes themselves).
// UNTIL IT IS FOUND, REFUSE THE SIGNATURE PATH ABOVE A MEASURED-SAFE B.
g_sig_on = 1;
g_sig_merged = 0;
}
}
// ---------- coarse superset bitmaps over every dim order (fringe filter) ----------
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
// Sweep the e-order ONCE per dim e, keeping (k-1) running bitsets (one per
// target dim d), and copy them out at each coarse-cell boundary. The old
// shape re-walked all N points for every one of the k(k-1) (d,e) pairs with
// two L3-random reads each (xe[pe[ptr]] and s_pos[d][pe[ptr]]); here the
// rank record is point-major so one 64B line yields every dim's rank, which
// cuts the random reads by ~k and the whole build by ~2-2.5x at k>=8.
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) {
__asm__ volatile("nop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\t");
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w))); }
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
// ---------- fold dimension + bucket vector per query ----------
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
__builtin_memcpy(s_bd + (size_t)i * MAXK + 6, &best, sizeof(best));
// Largest R that will ever be read out of each (fold dim, dim, bucket) array:
// the pooled array (ds,d,b) is only ever scanned over [0, maxR] words, so the
// per-array length can be far below the group maximum (0.38n vs 0.93n at k=10).
if (g_plen_on) {
const unsigned short *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
// group by ds
{
for (u32 v = 0; v <= (u32)k; v++) s_cntv[v] = 0;
for (u32 i = 0; i < N; i++) s_cntv[s_g[i] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)k; v++) { u32 c = s_cntv[v]; s_cntv[v] = ac; ac += c; } }
for (u32 i = 0; i < N; i++) { u32 id = s_qord[i]; s_tmp[s_cntv[s_g[id] + 1]++] = id; }
for (u32 i = 0; i < N; i++) s_qord[i] = s_tmp[i];
}
u32 gstart[MAXK + 2];
{ u32 c = 0; for (int d = 0; d <= k; d++) { gstart[d] = c; while (c < N && (int)s_g[s_qord[c]] == d) c++; } gstart[k + 1] = c; }
// ---------- transposed coordinate rows (sequential fringe) ----------
#if MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 8; // 32-byte aligned candidate rows for k<=8
#endif
size_t per = (size_t)N * KS * 4;
size_t lim = (size_t)150 << 20;
if (k >= 9) lim = (size_t)470 << 20; // k=9/10 need the rows: the random s_pm fallback is fatal
if (per * k <= lim) {
static u32 *pool = 0; static size_t cap = 0;
if (per * k > cap) { pool = (u32 *)pool_alloc(per * k + 256); cap = per * k; }
if (pool) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
if (KSx == 8 && k <= 8) {
_mm256_store_si256((__m256i *)dp, _mm256_load_si256((const __m256i *)row));
} else if (KSx == 6 && k == 6) {
const __m256i rv = _mm256_loadu_si256((const __m256i *)row);
_mm_storeu_si128((__m128i *)dp, _mm256_castsi256_si128(rv));
_mm_storel_epi64((__m128i *)(dp + 4), _mm256_extracti128_si256(rv, 1));
} else for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * k * ((size_t)N + 64);
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); sb[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
// ---------- signature planes, standalone only if the s_sdmp pass did not run ----
if (s_sig && !g_sig_merged && k == 6) {
static u8 cached_signature[(size_t)MAXN * 8] __attribute__((aligned(64)));
for (u32 i = 0; i < N; i++)
for (int e = 0; e < 6; e++) {
u32 v = (u32)(((u64)x[e][i] * g_sigmag) >> 32);
cached_signature[(size_t)i*8 + e] = (u8)(v ? v : 1u);
}
for (int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*k*((size_t)N + 64);
u32 j=0;
for (;j+8<=N;j+=8) {
__m128i r0=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+0]*8));
__m128i r1=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+1]*8));
__m128i r2=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+2]*8));
__m128i r3=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+3]*8));
__m128i r4=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+4]*8));
__m128i r5=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+5]*8));
__m128i r6=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+6]*8));
__m128i r7=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+7]*8));
__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);
_mm_storel_epi64((__m128i *)(base+(size_t)0*((size_t)N + 64)+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*((size_t)N + 64)+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*((size_t)N + 64)+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*((size_t)N + 64)+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*((size_t)N + 64)+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*((size_t)N + 64)+j),_mm_srli_si128(v2,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*8;
for(int e=0;e<k;e++)base[(size_t)e*((size_t)N + 64)+j]=row[e];
}
}
g_sig_merged=1;
}
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * k * ((size_t)N + 64);
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); base[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
}
// ---------- pool ----------
u64 poolwords = (u64)K1 * (Q + 1) * W;
#if MF_PLEN
if (g_plen_on) {
u64 best = 0;
for (int dsx = 0; dsx < k; dsx++) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < dsx) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)dsx * MAXK + d) * MAXQ;
for (u32 b = 0; b <= Q; b++) {
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > W) len = W;
tot += len;
}
}
if (tot > best) best = tot;
}
if (best && best < poolwords) poolwords = best;
}
#endif
size_t need = (size_t)poolwords * sizeof(u64);
static u64 *pool = 0; static size_t cap2 = 0;
if (need > cap2) { pool = (u64 *)pool_alloc(need + 4096); cap2 = need; }
s_bs = pool;
// ---------- groups ----------
for (int ds = 0; ds < k; ds++) {
u32 q0 = gstart[ds], q1 = gstart[ds + 1];
if (q0 >= q1) continue;
#if MF_WTRIM
{
u32 mx = 1;
for (u32 qi = q0; qi < q1; qi++) { u32 rr = folded_rank(s_qord[qi]); if (rr > mx) mx = rr; }
u32 wp = ((mx >> 6) + 8u) & ~7u;
Wpool = (wp < W) ? wp : W;
}
#endif
#if MF_PLEN
// Per-(dim,bucket) array lengths inside the group.
if (g_plen_on) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)ds * MAXK + d) * MAXQ;
u32 *po = s_pofs + (size_t)e * MAXQ;
u32 *pl = s_plen + (size_t)e * MAXQ;
for (u32 b = 0; b <= Q; b++) {
// NOTE: bucket 0 is the EMPTY prefix set (all-zero array), but a query whose
// bd[d]==0 still scans [0,R) words of it, so the array must be long enough --
// sizing it to the fixed 8-word minimum reads past the zeroed region into the
// next array. Use the real per-bucket max R here, not 0.
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > Wpool) len = Wpool;
po[b] = (u32)tot; pl[b] = len; tot += len;
}
}
if (tot > (u64)K1 * (Q + 1) * W) { g_plen_on = 0; } // never; safety
}
#endif
// bit numbering = rank in dim ds
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
u64 *base = g_plen_on ? s_bs : (s_bs + (size_t)e * (Q + 1) * Wpool);
u64 *cur = s_acc;
for (u32 w = 0; w < Wpool; w++) cur[w] = 0;
if (g_plen_on) {
// bucket 0 is the empty prefix set; bd[] can be 0 so it must be materialised as zeros
u64 *d0 = s_bs + s_pofs[(size_t)e * MAXQ];
u32 l0 = s_plen[(size_t)e * MAXQ];
for (u32 w = 0; w < l0; w++) d0[w] = 0;
}
u32 ptr = 0;
#if MF_BUILD_NOPOOL
while (ptr < N) ptr++;
#else
for (u32 b = 1; b <= Q; b++) {
// {q : x_d[q] < BVAL(d,b)} is EXACTLY the prefix [0,FPOS(d,b)) of the dim-d order
// (FPOS is the tie-group start containing rank b*B), so the per-point value compare
// -- one random xd read per point per (ds,d) pair, k(k-1)N of them -- is redundant.
// b==Q means "every point" and FPOS(d,Q) is the 0 sentinel, hence N.
// Credited to agent a7f799ceb7e63b0cd (worth 4-5% at k=4/5 on the board).
u32 pend = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) {
if (ptr + 64 < MAXN) __builtin_prefetch(s_nid + od[ptr + 64], 0, 3); // [PF-NID]
u32 q = s_nid[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++;
}
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
#if MF_NT
// The pool is written once and read 89 GB later, so a plain copy pays a
// read-for-ownership DRAM read per line for nothing. NT stores remove it
// (measured: this loop was 4.58 GB / ~1.7 s, i.e. 2.7 GB/s, vs 12.5 GB/s
// for a plain memset).
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 4); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 5); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf_sp<6>((u32)n, (const unsigned **)x, out); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 7); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 9); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 10); }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 434.608 ms | 67 MB + 420 KB | Accepted | Score: 100 | 显示更多 |