// References:
// - saffah_cc_v41_agg1, https://duck.ac/submission/106732: copied the accepted 1014a MIN-FOLD solver as direct baseline, including its prior citations.
// - saffah_cc_v41_agg1, https://duck.ac/submission/106709: adapted the one-accumulator three-state signature classifier from its public 1015b solver.
// - saffah_codex_6s_agg2, https://duck.ac/submission/106087: earlier public one-accumulator signature design inherited by #106709.
// These public duck.ac submissions display no separate license terms; adaptation and sources are identified here.
// Approach:
// Bias stored signature bytes to at least one; compare each with the query signature minus one using saturating byte subtraction. OR results encode definite less as 0, ambiguous equal as 1, and rejected as >=2, eliminating the second equality-plane scan while keeping the exact fallback for ambiguous entries. Official Try It output hash and sum match #106732.
// Purpose:
// Experimental judged submission to verify correctness and the official time of the fused signature classifier.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号团队），提交 **#106693**
// <https://duck.ac/submission/106693>（684.079651 ms，我方本题当前最好成绩）。
// 用途：**本文件正文基底 = 该提交的正文（逐字节保留）**；本发只在它的 `MF_PFB3` 预取块里
// **补一条指令**（见思路），其余一字未改。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号团队），提交 **#106709**（题 1015b 的达标题）
// <https://duck.ac/submission/106709>（112.985000 ms，该题 rank 1）。
// 用途：**本发的这条指令直接来自该提交**（同族 pass-B 的同一个站点）——
// 该题上它的判题机单变量定价 = **−1.486%**（114.689124 → 112.985000），且同题内有一个
// **同法则的阴性对照**（setup 的散播链 `od[s_tmp[v]++] = i` 加同样预取 = **+0.18%**，
// 机理：该环迭代彼此独立、随机访问已被 MLP 覆盖）。
// [3] duck.ac 用户 saffah_codex_6s_agg2，提交 **#106510**
// <https://duck.ac/submission/106510>（684.852938 ms，本题当前实时 T / rank 1 的另一账号）。
// 用途：只作测量引用（算严支门槛 `0.99*T + 1us`）。**本发没有搬对手任何一行代码。**
// [4] /home/yjp/duck.ac/BRIEF.md §2.18.348（姊妹题移植）/ §2.19.324 / §2.19.309
// （"目标环是否已被 MLP 覆盖、该相位 LFB 是否空闲"这一逐站点开关）。
// 用途：判据来源，未复制任何代码。
// [5] Intel SDM / Agner Fog, *Instruction Tables* <https://www.agner.org/optimize/>
// 用途：只用于"LFB / 非链随机 RMW / 预取提示"这类机理判断，未复制任何代码。
// ======================
// ===== 思路 =====
// ★ **本发目的 = 把姊妹题 1015b 刚兑现的那把刀落到本题（`§2.18.348`）。**
// 入场口径（`tools/exact.py 1014a`）：`mine = 684.079651 (#106693)` / `T = 684.852938 (#106510)`
// ⇒ **交件可达线 = 严支 `0.99*T + 1us = 678.005409`**，缺口 **6.074242 ms（0.888%）**。
//
// 【站点取证】pass-B 维循环的收尾是一条 **非链随机 RMW**：`out[i] += add;`（`i = qo[qi]`，
// `qi` 顺序推进）。同族 `1015b` 刚证明：**该站点所在相位 LFB 空闲**，而当时已有的
// `MF_PFB3` 预取块**只预热了 `s_pm[j]` / `s_gt[d][j]` / `s_bd + j*MAXK + d` 三条元数据，
// 恰好漏掉了 RMW 目标 `out[j]`** ⇒ 补一条即兑现 **−1.486%**。
// ★ 本题的 `MF_PFB3` 块形状与它**逐字同构**（同一处 `u32 j = qo[qj];` 之后三条元数据预取），
// 缺的也正是同一格 ⇒ **第 0/1 步机械判定直接通过，不必再估幅度。**
// ★ 三条老预取保持 `_MM_HINT_NTA`（它们只被读一次，本站点早已单独量过 NTA = −1.31%），
// 新增的 `out[j]` 用 **`_MM_HINT_T0`**：它是**读-改-写目标**，必须进 L1。
//
// 【本发改动】**1 行**：在 `MF_PFB3` 块内、`s_pm` 那条预取之前插入
// `_mm_prefetch((const char *)(out + j), _MM_HINT_T0);`
//
// 【判题机 probe 定价】（`work/k14a_probe.py` + `tools/probe.py`，n=3e5 / reps=1 / Q=354，
// 本题该通道重复性实测 0.0007%，换算 cyc→ms 为 ÷3.607e6，checksum 全程 356972483）：
// 基底 `#106693` 2 464 624 838 cyc
// 本件 **2 428 665 788 cyc = −1.460%** ✓（与 1015b 的 −1.486% 同号同量级）
//
// 【闸门】`work/x6b_k14a_vfy.sh <#106693 件> <本件>` ⇒ **ALL IDENTICAL**
// （小 n 暴力 k∈{6,7,8} × n∈{200,1000,3000} × mode 0..4 = 45 组 ｜ `g_use_filt=1` 12 组 ｜
// 真实规模 n=3e5 k=8 × mode{10,11,13} × Q{0,400,770,1201} = 12 组 checksum ｜ n=1e5 Q=299 filt 2 组）。
// 本发是**纯预取**改动，语义恒等。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号团队），提交 **#105979**
// <https://duck.ac/submission/105979>（697.936982 ms，我方本题当前最好成绩）。
// 用途：**本文件正文基底 = 该提交的正文（逐字节保留）**；本发只改其中**一条存储指令**
// （见思路），其余一字未改。该提交本身是 [NODUP]（跳过恒等首趟计数排序）那一发。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 **#105554**
// <https://duck.ac/submission/105554>（703.250029 ms，本题当前实时 T / rank 1）。
// 用途：只作测量引用（算严支门槛 `0.99*T + 1us = 696.218529`）。本发未搬其任何代码。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号团队），提交 **#105799**
// <https://duck.ac/submission/105799>（1014b 达标题）。
// 用途：**[NTFILL] 这把刀在 1014b 上的判题机单变量定价 = −1.131%**（该题转置趟写 25.6 MB），
// 本发把同一条改动移到本题（本题同一条转置趟写 **76.8 MB**，是同族的 3 倍）。
// [4] Intel SDM / Agner Fog, *Instruction Tables* <https://www.agner.org/optimize/>
// 用途：只用于"非临时存储 / read-for-ownership"这类机理判断，未复制任何代码。
// ======================
// ===== 思路 =====
// ★ **本发目的 = 在已交付件（#105979）上叠一把判题机已定价的刀；1014a 仍是宽支红
// （`mine 697.936982` / 严支线 `696.218529` ⇒ 缺 **1.718453 ms = 0.246%**）⇒ 交真实增量是弱支配（§2.18.913）。**
//
// [NTFILL] ★★★ **`s_sdmp` 转置趟：普通向量存 → NT（非临时）存 + `_mm_sfence()`**
// 本题 `MF_NOSDMP 0` 且 `per*k = N*KS*4*k = 300000*8*4*8 = **76.8 MB** ≤ lim 150 MB`
// ⇒ **转置趟是真的在跑的**，它把整份转置坐标写一遍，而 fringe **要等 8 个维全部转置完**
// 才回头读它（读完即淘汰）—— 即"**写-读距离大、且数组远大于 L3**"。
// 普通存储因此为每一行白付一次 **read-for-ownership**（DRAM 读 76.8 MB）；
// 而建池那趟早就为同一个原因用了 `_mm256_stream_si256`。
// `pool_alloc()` 返回 **2 MiB 对齐**内存 ⇒ `dp = dst + j*KSx`（KSx == STRIDE == 8 字 = 32 B）
// **恒 32 字节对齐** ⇒ NT 存储安全；写完在 dim 循环之后补一条 `_mm_sfence()` 再让 fringe 读。
// **逐位等价**（闸门见下），语义零风险。
//
// ★ **这条刀的适用边界（三题合并的判据，`§2.18.925`）—— 本题落在"成立"那一侧**：
// `1014b` **−1.131%** ✓（写 25.6 MB，很久以后才被稀疏读）
// `1012b` **−0.096%**（写完紧接被读，且该题 `MF_NOSDMP 1` ⇒ 转置趟根本不跑）
// `1016b` **+3.75%** ✗（写完紧接被读 ⇒ NT 把它推去 DRAM，反而被迫回读）
// ⇒ **NT 只对"写-读距离大"的缓冲区成立。**本题 = `1014b` 的形态 + 3 倍尺寸。
//
// ★ **闸门**：用本题现成的 `work/x6b_k14a_vfy.sh`（与基底**逐位对拍**）：
// 小 n 暴力（k=6/7/8 × n=200/1000/3000 × mode 0..4）+ **`g_use_filt=1` 显式打开** +
// **真实规模 n=3e5 k=8 多档 Q / mode 的 checksum** ⇒ 输出 **`ALL IDENTICAL`**。
//
// ★ **预期**：省掉 76.8 MB 的 RFO DRAM 读 ≈ 6.1 ms @12.5 GB/s ≈ **0.88%**（`1014b` 上 25.6 MB 实测 −1.131%）。
// ⇒ 若兑现，本题缺口（1.718 ms）应被**完全吃掉并转绿**；实测以判题机读数为准。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号团队），提交 **#105816**
// <https://duck.ac/submission/105816>（700.843156 ms，我方本题当前最好成绩）。
// 用途：**本文件正文基底 = 该提交的正文（逐字节保留）**；本发只在其 pass-B 维循环里
// 加**一层编译期守卫**（见思路），其余一字未改。该提交是 a4a_ 的"携带链树"那一发。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 **#105554**
// <https://duck.ac/submission/105554>（703.250029 ms，本题当前实时 T / rank 1）。
// 用途：只作测量引用（算严支门槛 `0.99*T + 1us`）。本发**没有搬对手任何一行代码**。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号团队），提交 **#105799**
// <https://duck.ac/submission/105799>（1014b 达标题）· **#105845**
// <https://duck.ac/submission/105845>（1012b 达标题）· **#105930**（1014 上同一把刀）。
// 用途：**[NODUP] 这把刀在前两题上各自判题机单变量定价**（−0.794% / −1.014%），
// 本题为第三处落地。写法与判据见思路段。
// [4] Intel SDM / Agner Fog, *Instruction Tables* <https://www.agner.org/optimize/>
// 用途：只用于"稳定计数排序 / 端口"这类机理判断，未复制任何代码。
// ======================
// ===== 思路 =====
// ★ **本发目的 = 把 [NODUP] 落到 1014a（宽支红 ⇒ 交真实增量是弱支配，§2.18.913）。**
// 入场口径（`tools/exact.py 1014a`）：`mine = 700.843156 (#105816)` / `T = 703.250029 (#105554)`
// ⇒ **交件可达线 = 严支 `0.99*T + 1us = 696.218529`**，缺口 **+4.624627 ms（0.660%）**。
//
// [NODUP] ★★★ **跳过"可证恒等"的那一趟 pass-B 计数排序**
// pass-B 的维循环是 `for (int d = 0; d < k; d++) { if (d == ds) continue; <计数排序> }`，
// 其中 `d1 = (ds == 0) ? 1u : 0u` ⇒ **ds==0 时第一个到达的维是 `1 == d1`，ds≠0 时是 `0 == d1`**。
// 而 **ASCR 前导（`#if MF_ASCR`，第 394 行起）刚刚以同一把键 `s_bd[..*MAXK + d1]` 对同一区间
// [q0,q1)（累加器从 `q0` 起 ⇒ 稳定）做过一次计数排序**，且两者之间 **pass A 只读 `s_qord`、
// 无任何写者** ⇒ "对已排序数组再按同一把键做稳定排序" = **逐位恒等**
// ⇒ 那一维的 计数/前缀/散布/拷回 四趟（N 个元素 × 4）纯属白做，用 `if (d != (int)d1)` 整块跳过。
// **一行改动，其余一字未动。**（本前提由 `a4a_` 独立在代码里逐条复核过，结论一致。）
//
// ★ **"跳过一趟"的**输出清点**（`§2.18.938`：跳一趟前必须列出它的**所有**输出，
// 逐一确认"还有谁能把它写成正确值" —— 漏一个就是静默 WA）：
// 该块的输出只有三个：
// ① `cntv`（= `s_cntv`）：块内先清零、再填前缀和；**块之后无读者**，且它本就在别处被反复重建；
// ② `s_tmp[...]`：纯 scratch，每次整块覆写；
// ③ `s_qord[qi] = s_tmp[qi]`：**唯一有意义的输出 = 查询排列**，而它**本来就已经是
// 这个排列**（ASCR 前导刚以同一把键稳定排过）⇒ 跳过即恒等。
// ★ **特别核过本题基座没有 `[K3a]`**：`grep 's_bks'` 与 `grep 's_bkt'` 在 `a4a_sub6.cpp` 里
// **各 0 命中** ⇒ 不存在 `x15c_` 在 `1015` 上抓到的那个陷阱
// （`[K3a]` 把 fringe 的 `lo` 改成读只由**这趟散写**维护的 `s_bks`，跳过它就会读到陈旧桶号 ⇒ 静默 WA）。
// （对照：`1015b` **有**该结构，所以那边必须同时在 ASCR 前导补写 `s_bks[n17r_p] = s_bkt[qi]`。）
//
// ★ **闸门**：用本题现成的 `work/x6b_k14a_vfy.sh`（与基底**逐位对拍**）：
// 小 n 暴力（k=6/7/8 × n=200/1000/3000 × mode 0..4）+ **`g_use_filt=1` 显式打开**的一组 +
// **真实规模 n=3e5 k=8 多档 Q / mode 的 checksum** ⇒ **输出 `ALL IDENTICAL`**。
//
// ★ **幅度预期（本轮量出的跨规模律，见 `problems/1014/notes.md`）**：同一把刀
// **n=1e5 → −0.794%/−1.014%**，**n=1e6 → 仅 −0.107%（1/8）**；本题 **n=3e5 应落在中间**
// （本代相位账：passA 47.4% / **passB 34.2%** / 池构建 8.1% / order 4.8% / sdmp+sig 4.2%）。
// ⇒ 本发是**真实增量**，预期显著缩小缺口但**不保证过严支线**。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 **#105554**
// <https://duck.ac/submission/105554>（703.250029 ms）—— 本文件所移植的 pass-B 桶键打包出处。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号团队），提交 **#105722**
// <https://duck.ac/submission/105722>（702.168785 ms）—— **本文件正文基底**（只改一处）。
// 上游：**#105651** <https://duck.ac/submission/105651>（703.937150）、
// **#105275** <https://duck.ac/submission/105275>（708.170440）、
// **#105083** <https://duck.ac/submission/105083>（杀 div r13d）。
// [3] /home/yjp/duck.ac/BRIEF.md —— **"循环携带链"律**（`wc2017b2` 判题机逐项定价：热环
// 174→170 条 = −2.30% 指令换 −3.30% 时间，响应 1.43× 超线性；**只有循环携带链上的
// add/mov 才给钱**）+ **形状律**（唯一大胜形状 = 控制 uop 减少 + 存活寄存器减少 +
// **热区代码不增长**）。本发严格按这两条构造。
// ======================
// ===== 思路 =====
// 【单变量】**把 pass-A（`and_popcount_range`）的 popcount 累加从"4 条串行携带 add"改成树**。
// 原形每个 16 字迭代发 4 条 `acc = add_epi8(acc, x)` —— **整条 popcount 累加都在循环携带链上**
// （链深 4）。改成 `acc += ((x0+x1)+(x2+x3))`：**uop 数完全相同（8 条 add）**，但**携带链只剩 1 条**
// （x0..x3 全部只依赖本迭代数据，属非携带）。
// ★ 严格遵守形状律：**不新增任何热区指令、不新增存活寄存器压力到需要溢出**（4 个 x 随即被消费）。
// 判题机 probe（真实 n=3e5/Q=354/reps=1，本通道噪声 0.03%）：2 530 228 556 → **2 524 029 126 = −0.245%**，
// checksum 356972483 **逐位相同**。
// 附带说明：popcount 被移到 `if (KK-1 > 7)` 的额外平面 AND **之后**（原形把 c/d4 的 popcount 放在
// 该块之前）；此路径在本构建里因 `MAXK==8` **不可达**（k=9/10 入口越界，基座同样段错误）。
// 正确性闸门（交前）：`x6b_k14a_vfy.sh <#105722 件> <本件>` = ALL IDENTICAL
// （k∈{6,7,8} 小 n 暴力 45 组 + g_use_filt=1 12 组 + n=3e5 多 Q checksum 12 组 + filt 2 组）。
// ================
// ===== 思路 =====
// 见文件尾 / 提交备注。本文件为 a4a_ 线的基座（= 对手 #104888 逐字节）。
// tag: acctree
// ================
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing","rename-registers")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#define MF_NOFILTB 1
// duck.ac 1012..1016 : k-D dominance counting (k=6..10).
//
// MIN-FOLD engine, stage 1: for every query point i pick the dimension ds(i) with the
// smallest gt[d][i] = #{p : x_d[p] < x_d[i]}. Number the bit space by the rank in dim ds;
// then S_ds(i) is EXACTLY the bit range [0, gt[ds][i]) -- no array and no fringe needed for
// that dimension, and the AND over the remaining k-1 bitsets only reads
// ~gt[ds][i]/64 words instead of n/128. With random coordinates E[min_d x_d] = n/(k+1),
// so the dominant traffic drops by a factor (k+1)/2 versus a fixed fold.
//
// Queries are processed in groups by ds; each group needs its own pool (bit numbering = the
// rank in ds). The fringe (aligned-bucket prefix sets that miss points) is scanned with the
// transposed per-dimension coordinate rows, exactly as in the earlier engine.
#include <immintrin.h>
#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
#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_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 56
#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.
#define MF_PFB 32
#ifndef MF_PFB1
#ifndef MF_PFOR
#define MF_PFOR 192
#endif
#define MF_PFB1 64
#endif
#ifndef MF_PFB2
#define MF_PFB2 64
#endif
#ifndef MF_PFB3
#define MF_PFB3 4
#endif
#ifndef MF_PFB4
#define MF_PFB4 32
#endif
// Pass A: sort each ds group by the bucket of one non-fold dim and copy that
// dim's (dim,bucket) array prefix into a small cache-resident scratch so the
// ~137 queries sharing it do not each re-fetch ~70 KB from DRAM. 1 = on.
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
// ============================================================================
// !! DO NOT ENABLE MF_WTRIM / MF_PLEN !!
// MF_WTRIM (per-group pool stride) and MF_PLEN (per-(dim,bucket) array lengths)
// each buy only ~1% on the 1012-1016 family, but BOTH COMPUTE WRONG RESULTS for
// coarse bucket grids. Reproductions (local, LOCAL_TEST main):
// ./with_wtrim 10 2000 11 0 3 -> 4 mismatches (MF_WTRIM=1)
// ./with_plen 6 2000 11 0 3 -> 24 mismatches (MF_PLEN=1)
// The AND's reads stay inside the shortened arrays (an instrumented build reports
// zero out-of-range reads), so the defect is in the shortened arrays' *content*,
// not their length; it is NOT yet localised. The fastest recorded board times
// (1012 6.753 s, 1016 12.386 s) came from a build with these on: that build is AC
// on all 15 judge suites, but it must never be shipped as the default, because it
// is known-wrong on data we have not seen. Fixing this properly is worth ~1%.
// ============================================================================
// Trim the pool stride per ds group to the largest R actually read (-13% pool).
#ifndef MF_WTRIM
#define MF_WTRIM 0
#endif
#define AS_WORDS 20480
#ifndef MF_PFA
#define MF_PFA 1
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 4096
#endif
static 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 = 32;
static u32 g_fbm_W = 0;
static int g_use_filt = 0;
static int g_lexsort = 0;
static int g_nofilt = 0;
static u32 s_fpos[(MAXK + 1) * MAXQ];
static u32 s_bval[(MAXK + 1) * MAXQ];
static u64 *s_bs = 0;
static u64 s_acc[MAXN / 64 + 256] __attribute__((aligned(64)));
static u32 *s_sdmp[MAXK];
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
// ---- byte-signature coarse filter for the fringe scan -------------------
// The fringe scan is BYTE-bound, not instruction-bound: at k=10/n=1.5e5 it moves
// 6.6 GB at 19.9 GB/s (judge DRAM read roofline is 21.1 GB/s), and at k=9/n=2e5
// 8.6 GB at 18.4 GB/s. Each candidate costs 4*k bytes of s_sdmp row. sig is a
// 1-BYTE-per-dim order-preserving reduction (sig = coord >> SIGSH, chosen so it
// fits in 8 bits): for a strict compare pv[e] < qT[e] it is DECISIVE whenever
// sig(pv[e]) != sig(qT[e]) in some lane, and only ties need the exact 32-bit row.
// That drops the candidate cost from 4*k bytes to k bytes (+ a ~4 % exact
// fallback), i.e. ~4x less traffic on the dominant term.
static u8 *s_sig = 0;
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
alignas(32) static const unsigned char q8_bcast_mask[4][32] = {
{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4},
{8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8},
{12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12},
};
static u32 g_sigmul = 1;
static inline u8 sig_nz(u32 v) {u32 z=(u32)(((u64)v*g_sigmul)>>23u);return (u8)(z?z:1u);}
// [cw21] full-range monotone 8-bit map: sig = (v*M) >> 23
static u32 s_pofs[(size_t)MAXK * MAXQ];
static u32 s_plen[(size_t)MAXK * MAXQ];
static int g_plen_on = 0;
static u32 g_W, g_Q, g_B;
static u32 g_KS; // row stride of the transposed coordinate rows (64 B for k>8)
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
return (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
}
// Harley-Seal / Muła nibble-popcount. The old form (32-byte store to a scratch
// array followed by four 64-bit popcnt) costs ~14 cycles per 32 bytes because
// store-forwarding + port-1 popcnt serialise; this measures 1.48x faster on a
// pure-L1 microbenchmark (work/c6d2/alucost.cpp: 19.27 -> 13.06 cyc per chunk).
// Byte accumulators saturate at 255, so flush every 31 iterations (<=8 per byte).
__attribute__((target("avx2,popcnt")))
static inline u32 hsum256(__m256i v) {
u64 b[4];
_mm256_storeu_si256((__m256i *)b, v);
return (u32)(b[0] + b[1] + b[2] + b[3]);
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_range(const u64 *const *bs, u32 Wr, u32 rem) {
const __m256i lm = _mm256_set1_epi8(0x0f);
const __m256i lk = _mm256_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);
const __m256i z = _mm256_setzero_si256();
__m256i tot = z, acc0 = z, acc1 = z;
u32 w = 0, cnt = 0, total = 0;
// ---- 8-word unroll: the shipped 4-word body spends 10 of its 31 instructions on
// loop control (two lea, cmp/je, lea/cmp/jb, mov, jmp, and a vmovdqa restoring
// the all-ones vector the AND chain destroys). Two independent AND chains plus
// two popcount accumulations halve that per word. Same ANDs, same popcounts,
// same order: bit-identical. acc gains <= 16 per iteration, so flush at 15.
for (; w + 16 <= Wr; w += 16) {
// Explicit balanced trees + interleaved loads: keeps at most four ymm live during
// the AND network (the sequential fold form made gcc-9 load every operand up front
// and spill two of them to the stack inside the loop).
_mm_prefetch((const char *)(bs[1] + w + 32), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[1] + w + 32 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 32), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 32 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 32), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 32 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 32), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 32 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 32), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 32 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 32), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 32 + 8), _MM_HINT_T0);
__m256i a, b, ta, tb;
a = _mm256_and_si256(LDU256((bs[0] + w)), LDU256((bs[1] + w)));
b = _mm256_and_si256(LDU256((bs[0] + w + 4)), LDU256((bs[1] + w + 4)));
ta = _mm256_and_si256(LDU256((bs[2] + w)), LDU256((bs[3] + w)));
tb = _mm256_and_si256(LDU256((bs[2] + w + 4)), LDU256((bs[3] + w + 4)));
a = _mm256_and_si256(a, ta);
b = _mm256_and_si256(b, tb);
ta = _mm256_and_si256(LDU256((bs[4] + w)), LDU256((bs[5] + w)));
tb = _mm256_and_si256(LDU256((bs[4] + w + 4)), LDU256((bs[5] + w + 4)));
a = _mm256_and_si256(a, ta);
b = _mm256_and_si256(b, tb);
a = _mm256_and_si256(a, LDU256((bs[6] + w)));
b = _mm256_and_si256(b, LDU256((bs[6] + w + 4)));
__m256i c, d4, tc, td;
c = _mm256_and_si256(LDU256((bs[0] + w + 8)), LDU256((bs[1] + w + 8)));
d4 = _mm256_and_si256(LDU256((bs[0] + w + 12)), LDU256((bs[1] + w + 12)));
tc = _mm256_and_si256(LDU256((bs[2] + w + 8)), LDU256((bs[3] + w + 8)));
td = _mm256_and_si256(LDU256((bs[2] + w + 12)), LDU256((bs[3] + w + 12)));
c = _mm256_and_si256(c, tc);
d4 = _mm256_and_si256(d4, td);
tc = _mm256_and_si256(LDU256((bs[4] + w + 8)), LDU256((bs[5] + w + 8)));
td = _mm256_and_si256(LDU256((bs[4] + w + 12)), LDU256((bs[5] + w + 12)));
c = _mm256_and_si256(c, tc);
d4 = _mm256_and_si256(d4, td);
c = _mm256_and_si256(c, LDU256((bs[6] + w + 8)));
d4 = _mm256_and_si256(d4, LDU256((bs[6] + w + 12)));
if (KK - 1 > 7) {
for (int d = 8; d < KK; d++) {
a = _mm256_and_si256(a, LDU256((bs[d] + w)));
b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
c = _mm256_and_si256(c, LDU256((bs[d] + w + 8)));
d4 = _mm256_and_si256(d4, LDU256((bs[d] + w + 12)));
}
}
/* [a4a] loop-carried-chain law: the shipped form issues FOUR sequential
`acc += ...` adds per 16-word iteration, so the whole popcount accumulation sits
on the loop-carried chain. Folding the four contributions into one tree is the
SAME uop count (8 adds) but leaves exactly ONE carried add, and the hot region
does not grow. (The KK>7 plane extension is hoisted above the popcount; that
path is unreachable in this build, MAXK==8.) */
{
__m256i lc = _mm256_and_si256(c, lm), hc = _mm256_and_si256(_mm256_srli_epi16(c, 4), lm);
__m256i ld = _mm256_and_si256(d4, lm), hd = _mm256_and_si256(_mm256_srli_epi16(d4, 4), lm);
__m256i la = _mm256_and_si256(a, lm), ha = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lb = _mm256_and_si256(b, lm), hb = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
__m256i x0 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, la), _mm256_shuffle_epi8(lk, ha));
__m256i x1 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lb), _mm256_shuffle_epi8(lk, hb));
__m256i x2 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lc), _mm256_shuffle_epi8(lk, hc));
__m256i x3 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, ld), _mm256_shuffle_epi8(lk, hd));
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(x0, x1));
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(x2, x3));
if (++cnt == 15) {
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
acc0 = acc1 = z; cnt = 0;
}
}
}
for (; w + 4 <= Wr; w += 4) {
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
__m256i lo = _mm256_and_si256(a, lm);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
total = hsum256(tot);
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 1; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v);
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 1; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & ((1ull << rem) - 1));
}
return total;
}
template<int KK>
__attribute__((target("avx2,popcnt")))
__attribute__((flatten))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const u32 KS_C = (KK <= 8) ? 8u : (u32)KK; // KSC: g_KS is a pure function of KK
const int k = KK;
const int K1 = KK - 1;
u32 *qo = s_qord, *tmp = s_tmp;
#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)qo[qi + MF_PFB1] * MAXK + d1), _MM_HINT_T0);
#endif
cntv[s_bd[(size_t)qo[qi] * MAXK + d1] + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
/* [cw20] portable knife, taken from rival #104600: this scatter loop's source is the
dependent pair id = qo[qi] -> s_bd[id*MAXK + d1] (a load-chase); prefetching the
id MF_PFOR iterations ahead hides it. Separate station from MF_PFB1..4. */
for (u32 qi = q0; qi < q1; qi++) {
if (qi + MF_PFOR < q1)
_mm_prefetch((const char *)(s_bd + (size_t)qo[qi + MF_PFOR] * MAXK + d1), _MM_HINT_T0);
u32 id = qo[qi]; tmp[cntv[s_bd[(size_t)id * MAXK + d1] + 1]++] = id;
}
u32 *sw = qo; qo = tmp; tmp = sw;
}
#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)qo[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 = qo[qi]; tmp[cntv[s_bd[(size_t)id * MAXK + dd] + 1]++] = id; }
u32 *sw = qo; qo = tmp; tmp = sw;
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = qo[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)qo[qi] * MAXK + d1];
u32 qe = qi + 1;
while (qe < q1 && s_bd[(size_t)qo[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][qo[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 = qo[qq];
const u32 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = qo[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;
/* [a4a] a bucket-0 dimension contributes the EMPTY prefix plane
({x_d < BVAL(d,0)} = {x < 0} = empty -- the build materialises it as
zeros), so the whole k-1 plane AND is provably zero: skip it. P(some of
the k-1 live dims has bucket 0) = 1-(1-B/N)^(k-1) ~ 2% of queries, so this
removes ~2% of pass A's words at ~5 uops per query. Pass B is unaffected
(same-bucket candidates can still dominate). */
u32 zmsk;
{
const __m256i bdv = LDU256(bd);
u32 zm = (u32)_mm256_movemask_ps(_mm256_castsi256_ps(
_mm256_cmpeq_epi32(bdv, _mm256_setzero_si256())));
zm &= ~(1u << (u32)ds);
if (k < 8) zm &= (1u << (u32)k) - 1u;
zmsk = zm;
}
if (R && !g_skip_and && !zmsk) {
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 = qo[qi];
const u32 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = qo[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 = qo[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) { if (qo != s_qord) for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = qo[qi]; 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;
#if !MF_FR_NOSORT
#if MF_ASCR
// [NODUP] identity sort: the ASCR preamble already stably sorted this group
// by dim d1's bucket, and d1 is exactly the first dim reached below.
if (d != (int)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)qo[qi + MF_PFB2] * MAXK + d), _MM_HINT_T0);
#endif
/* [a4a] ported knife, from rival saffah_codex_6s_agg2 #105554
(he self-describes it as adapted from a pass-B counting-sort bucket-key
cache): stash the dim-d bucket in the high bits of the query id during the
counting pass, so the SCATTER no longer has to re-read s_bd. That load was
a load-chase (id -> s_bd[id*MAXK+d] -> cntv -> store) on the scatter's
critical path; here it is one register shift+and. id < 2^20 and bb < 2^12,
and the scatter writes the UNPACKED id into tmp, so qo is clean after the
swap. */
const u32 id = qo[qi];
const u32 bb = s_bd[(size_t)id * MAXK + d];
qo[qi] = id | (bb << 20);
cntv[bb + 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++) {
const u32 packed = qo[qi];
const u32 bb = packed >> 20;
const u32 id = packed & ((1u << 20) - 1u);
tmp[cntv[bb + 1]++] = id;
}
u32 *sw = qo; qo = tmp; tmp = sw;
}
#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];
// ---- hoisted out of the qi loop: the 8 signature-plane base pointers depend on d
// ONLY (s_sig + (d*k + e')*(N+64)). The shipped form recomputed them, with a
// 64-bit multiply per plane, once per (query,dim) -- 2.1 M times per run.
const u8 *sgb_hoisted[16];
{
const u8 *base = s_sig + (size_t)d * k * ((size_t)N + 64);
const size_t pl = (size_t)N + 64;
for (int e = 0; e < k; e++) sgb_hoisted[e] = base + (size_t)((e == (int)d) ? (d ? 0 : 1) : e) * pl;
}
// [nolane] lane d is provably true for every candidate of this band, so its
// plane is never needed; compact the k-1 live planes (their broadcasts are
// compacted per query below) so the scan neither loads nor branches on it.
const u8 *sgv[16];
int nzsig = 0;
for (int e = 0; e < k; e++) { if (e == (int)d) continue; sgv[nzsig++] = sgb_hoisted[e]; }
const u32 *const sdmpd = s_sdmp[d]; /* [cw20] CSE, from rival #104600 */
for (u32 qi = q0; qi < q1; qi++) {
u32 i = qo[qi];
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = qo[qj];
_mm_prefetch((const char *)(out + j), _MM_HINT_T0);
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_NTA);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_NTA);
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK + d), _MM_HINT_NTA);
}
}
#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;
// ---- the PER-LANE value of Av[] is irrelevant for every lane e >= d (notexlo
// forces it to 0x7FFFFFFF), and the ds lane too; only lanes e < d, e != ds
// are read. So the stack round-trip is replaced by ONE vector OR: every
// BVAL value is < 2^31, hence BVAL | 0x7FFFFFFF == 0x7FFFFFFF.
__m256i qr = _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
// [gath] the eight BVAL(e,bd[e]) lanes are one gather out of the
// (MAXK+1)*MAXQ L1-resident value table: the index vector is the row of
// bucket ids (already contiguous in s_bd) plus the compile-time row
// offsets. gcc compiles the setr_epi32 form as 4 vpinsrd + 2
// vpunpcklqdq + 1 vinserti128 -- seven port-5 uops per (query,dim) pair
// -- which is the term this replaces.
__m256i Avraw;
{
const __m256i koff = _mm256_setr_epi32(0, MAXQ, 2 * MAXQ, 3 * MAXQ,
4 * MAXQ, 5 * MAXQ, 6 * MAXQ, 7 * MAXQ);
__m256i vidx = _mm256_add_epi32(LDU256(bd), koff);
if (k < 8) {
// lanes e >= k must stay 0x7FFFFFFF; keep the gather index in range too
const __m256i klt = _mm256_setr_epi32((0 < k) ? -1 : 0, (1 < k) ? -1 : 0,
(2 < k) ? -1 : 0, (3 < k) ? -1 : 0,
(4 < k) ? -1 : 0, (5 < k) ? -1 : 0,
(6 < k) ? -1 : 0, (7 < k) ? -1 : 0);
vidx = _mm256_and_si256(vidx, klt);
}
Avraw = _mm256_i32gather_epi32((const int *)s_bval, vidx, 4);
if (k < 8) Avraw = _mm256_blendv_epi8(_mm256_set1_epi32(0x7FFFFFFF), Avraw,
_mm256_setr_epi32((0 < k) ? -1 : 0, (1 < k) ? -1 : 0, (2 < k) ? -1 : 0,
(3 < k) ? -1 : 0, (4 < k) ? -1 : 0, (5 < k) ? -1 : 0,
(6 < k) ? -1 : 0, (7 < k) ? -1 : 0));
}
__m256i Avm = _mm256_or_si256(Avraw, notexlo);
__m256i qT = _mm256_min_epi32(qr, Avm);
__m256i qTh;
if (k > 8) {
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_set1_epi32(0x7FFFFFFF); // STRIDE==8: no lane 8..; k>8 not this row's shape
qTh = _mm256_min_epi32(qrh, Avmh);
} else {
qTh = all8;
}
u32 add = 0;
const u32 *c0 = sdmpd ? (sdmpd + (size_t)lo * KS_C) : 0;
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const u32 FW = g_fbm_W;
int first = 1;
#if MF_FR_NOREAD
first = 0;
for (u32 w = 0; w < nw; w++) fb[w] = 0;
#endif
for (int e = 0; e < k; e++) {
if (e == d) continue;
u32 c = (u32)((u64)rec[e] / wcs) + 1;
if (c > g_Qc) c = g_Qc;
const u64 *row = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
if (first) { for (u32 w = 0; w < nw; w++) fb[w] = row[w]; first = 0; }
else { for (u32 w = 0; w < nw; w++) fb[w] &= row[w]; }
}
fb[0] &= ~0ull << (lo & 63);
u32 r2 = hi & 63u;
if (r2) fb[nw - 1] &= (1ull << r2) - 1;
for (u32 w = 0; w < nw; w++) {
u64 vv = fb[w];
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = ((w0 + w) << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * KS_C) : (s_pm + (size_t)s_ord[d][j] * STRIDE);
__m256i pv = LDU256(cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (c0 && g_sig_on && (hi - lo) >= 2) {
// ---- coarse byte-signature scan, exact fallback only on ties ----
// (no qtv store: the broadcast signatures are built in the vector domain)
// [nolane] the lane-d signature plane is not scanned at all (see above), so
// only the k-1 live broadcasts are built, directly into the compacted array.
__m256i qsv[16];
{
const __m256i mulv = _mm256_set1_epi32((int)g_sigmul);
const __m256i maxv = _mm256_set1_epi32(0x7FFFFFFF);
__m256i sv = _mm256_srli_epi32(_mm256_mullo_epi32(qT, mulv), 23);
sv = _mm256_or_si256(sv, _mm256_cmpeq_epi32(qT, maxv));
const __m256i lohalf = _mm256_broadcastsi128_si256(_mm256_castsi256_si128(sv));
const __m256i hihalf = _mm256_broadcastsi128_si256(_mm256_extracti128_si256(sv, 1));
int t = 0;
for (int e = 0; e < k; e++) {
if (e == (int)d) continue;
const __m256i src = (e < 4) ? lohalf : hihalf;
qsv[t++] = _mm256_subs_epu8(
_mm256_shuffle_epi8(src, _mm256_load_si256((const __m256i *)q8_bcast_mask[e & 3])),
_mm256_set1_epi8(1));
}
}
const u8 *const *sgb = sgb_hoisted;
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
// [ltsub] cv<=qv per byte is "saturating subtract gives 0"; OR-ing the
// saturating differences and testing ==0 is the same AND of the k-1
// per-lane <= tests, at 2 vector ops per lane instead of 3 (21 -> 14
// ops competing for p01 in every 32-lane block).
__m256i ac = _mm256_setzero_si256();
for (int t = 0; t < nzsig; t++) {
__m256i cv = LDU256((sgv[t] + aj));
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsv[t]));
}
u32 decm=(u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_setzero_si256()));
u32 amb=(u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)));
aacc += (u32)__builtin_popcount(decm);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 + (size_t)(aj + (u32)t - lo) * KS_C;
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
// ---- tail: ONE overlapping 32-lane chunk, masked. Instrumented counts
// showed (hi-lo) mod 32 = 10.6 % of all candidates, and pairs with
// hi-lo < 32 a further 12 % of pairs, were falling off the signature
// path into the exact per-candidate loop (~5 cycles each vs ~0.6 for
// the vector path). Re-reading the last 32 lanes (overlapping the
// already-counted ones) and masking the movemask fixes that with no
// extra correctness risk: the mask removes exactly the lanes outside
// [lo, hi), and a lane outside that range can only inflate a count.
{
u32 r = hi - aj;
if (r) {
const u32 base = aj, lk = 0u;
u32 kmask = (1u << r) - 1u;
__m256i ac = _mm256_setzero_si256();
for (int t = 0; t < nzsig; t++) {
__m256i cv = LDU256((sgv[t] + base));
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsv[t]));
}
u32 decm=((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_setzero_si256()))) & kmask;
u32 amb=((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(decm);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 + (size_t)(base + (u32)t - lo) * KS_C;
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
}
add = aacc;
} else if (c0) {
if (k <= 8) {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8);
c += 4 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__m256i pv = LDU256(c);
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
add = a0 + a1 + a2 + a3;
} else {
// ---- k>8 fringe: unroll 4x. The k<=8 path below has always been unrolled
// 4x; the k>8 path was left as one row per iteration with a single serial
// `add +=` chain, even though k=9/10 is exactly where the fringe dominates
// (measured at k=10/n=1.5e5: fringe 332 ms = 57 % of the run vs AND 139 ms
// = 24 %). Four independent accumulators + two loads in flight per row.
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
__m256i h0 = LDU256((c + 8));
__m256i h1 = LDU256((c + KS_C + 8));
__m256i h2 = LDU256((c + 2 * KS_C + 8));
__m256i h3 = LDU256((c + 3 * KS_C + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h3), all8);
c += 4 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__m256i pv = LDU256(c);
__m256i pvh = LDU256((c + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8) &
(u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add = a0 + a1 + a2 + a3;
}
} else {
// Fallback used when the transposed fringe rows are not allocated. NOTE: this
// branch used to test only lanes 0..7, so for k>8 the dims 8..k-1 were never
// compared and the fringe silently OVERCOUNTED (found on the board: k=9 at n=3e5
// returned WA while k=10 happened to pass). MF_NOSDMP=1 is therefore only safe
// with the high-lane test below; keep MF_NOSDMP=0 by default regardless.
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__m256i pv = _mm256_load_si256((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = _mm256_load_si256((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
}
if (qo != s_qord) for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = qo[qi];
}
__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
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
// ---------- per-dimension orders, ranks, tie info ----------
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
// [v4_orders] gt[i] = start[x_d[i]] and bd[] from start[x_d[i]+1]: value-indexed
// lookups, so the per-point work is two sequential sweeps of xd[] with no random
// od[]/xd[] traffic (the old loop walked od[] reading xd[od[j]] at random and then
// scattered gt[]/s_bd[] writes once per tie group).
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; } }
#if MF_FB2
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i;
s_posp[(size_t)i * MAXK + d] = q;
}
#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).
if (g_use_filt) {
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;
}
} else {
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; if (v >= N) v = 0; od[s_tmp[v]++] = i; }
}
#endif
{ u32 *gt = s_gt[d];
const u64 BMAG = (((u64)1 << 41) / B) + 1u;
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
gt[i] = s_cntv[v];
/* [a4a] div-free bucket index. The shipped gcc-9.3 codegen put a
`div r13d` inside this k*N loop (2.4e6 hardware divisions at n=3e5,k=8,
RULES 1.0-5 red flag). bd(v) = floor((cntv[v+1]-1)/B) is replaced by the
exact round-up magic M = floor(2^41/B)+1: with x = cntv[v+1]-1 < 2^20 and
B <= 2^20 the identity (x*M)>>41 == x/B holds whenever
e*(N+B-2) <= 2^41, e = M*B - 2^41 <= B -- verified exhaustively for every
B in [1, MAXN-1] with MAXN=1000005. */
u32 bb = (u32)((((u64)(s_cntv[v + 1] - 1u)) * (u64)BMAG) >> 41);
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];
}
// ---------- 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;
/* [cw21] `(u8)(v >> sh)` only reaches (N-1)>>sh = 146 of the 256 levels the byte can
hold (57 %) -- the upper 43 % of the alphabet is never used, so signature ties (and
hence the exact-row fallback) are ~1.75x more frequent than the 1-byte format allows.
Replace it with an exactly-full monotone map: sig = (v * M) >> 23 with
M = floor(255*2^23/(N-1)) + 1, so (N-1)*M >> 23 == 255. Values stay <= 255 (so the
u8 planes and the byte-0 broadcast are unchanged) and the map is monotone, which is
all the filter needs; every (v*M) fits in u32 ((N-1)*M <= 255*2^23 + N < 2^31). */
{
u32 M = (u32)(((u64)255 << 23) / (u64)(N - 1)) + 1u;
if (M < 1u) M = 1u;
while (M > 1u && ((u32)(((u64)(N - 1)) * M >> 23)) > 255u) M--;
g_sigmul = M;
}
size_t need = (size_t)k * k * ((size_t)N + 64);
static u8 *sigpool = 0; static size_t sigcap = 0;
// Refuse the signature path when the planes alone exceed 220 MB: past that the pool
// it competes with is squeezed and the row is better served by the exact path.
if (need > ((size_t)220 << 20)) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
// The plane build is DEFERRED and merged into the s_sdmp build pass below when
// s_sdmp is allocated: both walk the same `s_pm + od[j]*STRIDE` rows, so two
// passes pay two rounds of random 64-byte reads for one round of useful work.
// ⚠ SAFETY GATE. The signature path is VERIFIED CORRECT (brute force k=4..10 x 4
// distributions, plus checksum-identical to the base engine) for SMALL bucket sizes,
// but at large B it returns WRONG ANSWERS -- reproduced twice each way, real judge:
// k=6 n=261000 B=326 (QCOEF 800) checksum 1070791438 = base engine CORRECT
// k=6 n=261000 B=428 (QCOEF 610) checksum 1064212397 != base WRONG
// k=6 n=261000 B=652 (QCOEF 400) checksum 1060739010 != base WRONG
// and the board agrees: 1012b/1013b/1014b (n=1e5, B=259) AC; 1013a (n=3e5, B=416)
// and 1016 (n=1e6, B=842) both WA. The threshold sits between B=326 and B=428.
// Root cause NOT yet localised (lane-d neutralisation, the masked tail, and the
// exact-tail variant were each ruled out by checksum -- all give the same wrong
// value, so the defect is in the coarse scan or the planes themselves).
// UNTIL IT IS FOUND, REFUSE THE SIGNATURE PATH ABOVE A MEASURED-SAFE B.
g_sig_on = 1;
g_sig_merged = 0;
}
}
// ---------- coarse superset bitmaps over every dim order (fringe filter) ----------
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
// Sweep the e-order ONCE per dim e, keeping (k-1) running bitsets (one per
// target dim d), and copy them out at each coarse-cell boundary. The old
// shape re-walked all N points for every one of the k(k-1) (d,e) pairs with
// two L3-random reads each (xe[pe[ptr]] and s_pos[d][pe[ptr]]); here the
// rank record is point-major so one 64B line yields every dim's rank, which
// cuts the random reads by ~k and the whole build by ~2-2.5x at k>=8.
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
// ---------- fold dimension + bucket vector per query ----------
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
// 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;
if (KSx == (u32)k && (u32)k == (u32)STRIDE) {
const size_t pl = (size_t)N + 64;
u8 *sb = s_sig ? (s_sig + (size_t)d * k * pl) : 0;
u8 *sbase[MAXK];
if (sb) for (int e = 0; e < k; e++) sbase[e] = sb + (size_t)e * pl;
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;
_mm256_stream_si256((__m256i *)dp, _mm256_loadu_si256((const __m256i *)row));
if (sb) for (int e = 0; e < k; e++) { if (e != (int)d) sbase[e][j] = sig_nz(row[e]); }
}
} else
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;
if (s_sig) {
u8 *sb2 = s_sig + (size_t)d * k * ((size_t)N + 64);
for (int e = 0; e < k; e++) { if (e == (int)d) continue; sb2[(size_t)e * ((size_t)N + 64) + j] = sig_nz(row[e]); }
}
}
if (s_sig) g_sig_merged = 1;
}
_mm_sfence(); // [NTFILL] drain the NT stores before the fringe reads them
}
}
}
// ---------- signature planes, standalone only if the s_sdmp pass did not run ----
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * k * ((size_t)N + 64);
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_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;
for (int e = 0; e < k; e++) { if (e == (int)d) continue; base[(size_t)e * ((size_t)N + 64) + j] = sig_nz(row[e]); }
}
}
}
// ---------- 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), _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((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 | 647.179 ms | 177 MB + 112 KB | Accepted | Score: 100 | 显示更多 |