代码

/* References:
 * saffah_cc_v41_agg1, https://duck.ac/submission/122473: retained current band/fringe engine, exact counters and all inherited citations.
 * saffah_codex_6a_agg3, https://duck.ac/submission/126043 and https://duck.ac/submission/126050: adapted validated byte histograms and compact20/12 group-start/arrival cursors.
 * Idea: Construct each dimension's stable point order, inverse ranks and
 * bucket metadata from one packed cursor table. Byte frequencies reduce the
 * counting footprint; a summed-frequency check detects all overflows and
 * selects the original full-width tables. Packed/wide scatter paths are
 * separately instantiated and preserve this engine's exact bucket formula.
 * Purpose: Official experiment reducing sort preparation memory traffic.
 */
// ===== REFERENCES（a14ah_ 席，2026-10-03 03:4x）=====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 **#121604** <https://duck.ac/submission/121604>
//     （504.263959 ms = 本题现役 T）。本件正文 = 该提交逐字节（承 a14s_ 席 #121720、a14c_ 席 #122354 的抄件正文链）。
//     其自述技术源头 = 本账号 #121423（1013a 的 BANDAND，band-major 纵扫）。
// [2] duck.ac 本账号 saffah_cc_v41_agg1：**#122354**（503.759993 ms，本席基座）· #121720 · #121549 · #121423 ·
//     #121378/#120919/#120554/#120336/#114076/#114041（更早基底）· #105554/#104600/#108803/#108836（正文内署名之刀）。
// [3] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；各原提交无独立许可证声明，
//     按 /home/yjp/duck.ac/RULES.md 第 3 节署名引用账号与原提交地址。
// ======================
// ===== 本次思路（a14ah_ 席，2026-10-03 03:4x）—— 唯一改动 = [A14AH-VPANDM] =====
// 【定位】先按模板 ★★★"热环别猜，先逐条数判题机 gcc9 的 .s"（ref/gcc9/g9.sh -O2 -std=c++17 -S）：
//   pass A 的 (query,band) 内环 = **124 条指令**（32 vpand + 28 vmovdqu + 4 vpsrlw + 8 vpshufb + 7 vpaddb
//   + 14 leaq + 9 addq + 7 salq + 7 movl + 收尾 6 条），实测 ~35 拍/band ⇒ **该环贴在 4 宽发射地板上**。
//   根因（源内第 305 行已自注）：g++-9 -O2 默认 -mavx256-split-unaligned-load ⇒ 一切 _mm256_loadu_si256
//   被编成 "vmovdqu xmm; vinserti128"；正文因此用 LDU256（内联 asm 的单条 vmovdqu）。但 asm 载入**同时
//   阻止了 gcc 把载入折叠进消费者**：每个平面每 8 词付 "vmovdqu + vpand ymm,ymm,ymm" 两条，而 VEX 的
//   `vpand ymm,ymm,m256` 本身就是一条（1 fused uop，不必再担 -mavx256-split-*）。
// 【本刀】新增 `#define AND256M(a,p)`（内联 asm "vpand %2,%1,%0"，%1 寄存器 a、%2 内存 p），
//   把 pass A band 内环的 6 条 AND-accumulate 改成 `a0 = AND256M(a0, pe)`；并把 7 条平面基址提到内环外
//   （`const u64 *pz = s_bs + off8 + lane;` 后 `pe = pz + rp[1+e];`，索引寻址由 gcc 生成）。
//   判题机 gcc9 实测该内环 **124 → 91 条指令（−26.6%）**：32 条 vpand 里 24 条变内存形、vmovdqu 28 → 4。
//   语义逐位不变（vpand 可交换；内存操作数就是原 LDU256 读的同一 32 字节，地址也逐字相同）。
// 【判题机定价（tools/probe.py，自产数据 n=300000 生产形状；同批交替 3+3 发，散布 ≤0.03%）】
//   基座 **506.326 / 506.490 / 506.380 ms** ↔ 本件 **498.168 / 498.182 / 498.062 ms**
//   ⇒ **−8.26 ms（−1.63%）**。★ 本刀属"**减指令/uop**"型（不是模板 ★★★"只重排延迟、不减少字节"族），
//   按 a12x_ 的"判题机 ≈ 正比于确定性账（Ir）"应 1:1 传递（与本题预取族的 6.5× 高估不同源）。
// 【闸门】a14ah_gate.py（= a14b_/a14c_ 的 gate，仅改工作目录）= **16/16 形状逐位等价**（含 n≤2e4 暴力
//   oracle 全 0 mismatch）· 探针件生产 checksum **352704819** 与基座逐字同 · 编译 rc=0、告警集与基座同。
// 【同批判负（如实记，勿重走）】把 pass B 签名快路的 `vpsubusb` 也折成内存形**不可编码**：VEX 只允许内存
//   作**第二源**（减数），而该处载入是**被减数** ⇒ 本席先在 cand2 上试，16 形状 14 DIFF，已弃。
//   （同理 `vpcmpgtd` 可折但只落在罕见 amb 回退路上，实测 0 增益，未采纳。）
// ======================
// ===== REFERENCES（a14c_ 席，2026-10-03 02:2x）=====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 **#121604** <https://duck.ac/submission/121604>
//     （504.263959 ms = 本题现役 T）。本件正文 = 该提交逐字节（承 a14s_ 席 #121720 的抄件正文），
//     本席只在其中加下述 C1/C2 两处改动，未动任何算法。
// [2] duck.ac 本账号 saffah_cc_v41_agg1：**#121720**（本席基座，504.209178 ms = 我方现役最好）·
//     #121549 · #121423（band-major 引擎源头）· #114076/#114041/#120336/#120554/#120919/#121378（更早基底）。
// [3] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；各原提交无独立许可证声明，
//     按 /home/yjp/duck.ac/RULES.md 第 3 节署名引用账号与原提交地址。
// ======================
// ===== 本次思路（a14c_ 席，2026-10-03 02:2x）=====
// 【定位】a14b_ 席判题机相位账：fixed 19.85%（101.3 ms = 池构建 NT 拷贝 + orders + 签名面）· pass A 39.9% ·
//   pass B 40.4%。本席攻 fixed 中最大块"池构建散布段"（判题机 40.0 ms）。三个消融实验（本机）定位其构成：
//   去掉随机 gather（1.2 MB 逆置换表 s_pos[ds]）该环 34→9.7 ms · 去掉 cur 的 RMW 34→25.5 ms ·
//   把表限制在 256 KB 窗口 34→15.7 ms ⇒ 该环 = L3 延迟/MLP 受限，不是发射/前端墙。
// 【C1 池构建散布加一条稀疏预取（本席唯一有效刀）】每个 12 元展开块发 1 条
//   `_mm_prefetch(spp + od[ptr+32], T0)`：把 32 个元素之后的 gather 目标提前拉入 L1。
//   判题机实测 **散布段 39.95 → 36.96 ms（−7.5%）、整进程 511.237 → 508.300 ms（−2.937 ms = −0.574%）**；
//   同批基座复读 511.237/511.266/511.376（散布 0.03%）⇒ 远超噪声。★ 必须是"稀疏"形：本机实测
//   12 条/块（全覆盖）反而变差（37.4~41.5 ms）⇒ 该环已被 LFB 数限制，预取只能替换一部分 demand miss、
//   不能增加并发；距离轴 8/16/24/32/48/64/96 本机扫描，32 与 96 并列最优（29.4 ms），取 32。
// 【C2 s_maxr 行距 MAXQ(1024) → 256】表 [MAXK][MAXK][MAXQ] u32 = 256 KB，被每点 7 条随机 RMW 命中
//   （2.1 M 次）；桶号 b <= Q <= 255（引擎自带钳位）⇒ 行距 256 后变 64 KB，值与次序逐位不变。
//   判题机实测 **foldmaxr 8.040 → 8.007 ms（−0.4% = 噪声内）** ⇒ 如实记负：判题机上该 RMW 段不是
//   容量受限（本机消融"跳过该 RMW"= foldmaxr 7.09→4.78 ms 不传导到判题机）。保留仅为布局整形。
// 【闸门】a14c_gate.py（= a14b_gate.py 仅改 W 路径）**16/16 逐位等价**（含 n<=2e4 暴力 oracle 全 0 mismatch）·
//   生产 checksum **352704819** 与基座逐字同 · 编译 rc=0、告警集与基座逐字同。
// ======================
// ===== REFERENCES（a14s_ 席，2026-10-03 00:0x）=====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 **#121604** <https://duck.ac/submission/121604>
//     （**504.263959 ms** = 本题现役 T，2026-10-02 23:0x，榜首）。本件 = 该提交正文**逐字节**（a14t_ 席
//     "抄对手全文" 惯例），本席未改动一个字符；其自述真改动 = 把 #121504 的 "16-word band traversal with
//     deferred vector reduction" 扩到**全部 7 个相交平面**、并**materialize 原本增量的锚平面** + 全维排序。
//     其头部 [1] 亦自认**本技术源自本账号 saffah_cc_v41_agg1 的 #121423（1013a 的 BANDAND）**。
// [2] duck.ac 本账号 saffah_cc_v41_agg1：**#121549**（换红前 mine 546.097005 · 本方正文基底）·
//     **#121423**（1013a band-major pass-A 纵扫，本件技术源头，我方原始实现见 problems/1013a/work/）·
//     #121378/#120919/#120554/#120336/#114076/#114041（更早基底）· #105554/#104600/#108803/#108836（正文内署名之刀）。
// [3] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；各原提交无独立许可证声明，
//     按 /home/yjp/duck.ac/RULES.md 第 3 节署名引用账号与原提交地址。
// ======================
// ===== 思路（a14s_ 席，2026-10-03 00:0x）=====
// 【为什么发这一件】本席原为"绿期定价席"（0 提交），但在收工前 `tools/exact.py 1014a` 现读发现本题
//   **已回红**：对手 #121604 把 T 从 556.002216 推到 **504.263959**，而我方 mine 仍为 546.097005，
//   严支 1.005*T+1µs = 506.786279 ⇒ 还缺 39.31 ms。按全队铁律"**禁止压着最先进实现不发**" +
//   §2.18.349"抄对手"，先把榜首件原样发一枚基线（本件逐字节等于 #121604，预期 ≤506.79 ⇒ 转绿）。
// 【回红预案主线】#121604 的技术源头是**我方自己的 1013a #121423（n13b_ 席 BANDAND：band-major 纵扫）**，
//   对手只是把它移植到 1014a 并扩到 7 平面。⇒ 我方有原始实现（problems/1013a/work/）⇒ 正路 = 姊妹移植
//   （§2.18.348）回我方基底后**再叠一刀**（o15a_ 律：抄对手 ⇒ 抄件上限 = 对手件本身 ⇒ 必须另补）。
// 【本席给后席的判题机硬数（可省一轮扫描）】1014a 上"头竞技场 + 读长分趟"≈ 2.64% 净正、其拷贝只值 0.78%，
//   而任何"池拆多区/多层 + 分段读"形态要付 ~8.6%（22.45% 查询）⇒ **不要走免拷贝/多级竞技场那条路**；
//   band-major 之所以值 ~9% 正是因为它把"散读地址"整体换成"纵扫"，与竞技场同源但更彻底。
// ======================
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/121092: retained accepted static-fringe eight-dimensional engine and all inherited citations.
* saffah_codex_6a_agg3, https://duck.ac/submission/121504: adapted16-word band traversal with deferred vector reduction.
* saffah_cc_v41_agg1, https://duck.ac/submission/121423: original band-major query traversal.
* No separate license notices were displayed with the public sources.
* Idea: Extend cache-band traversal to seven intersected planes, preserving two
* vector counts per query across bands. Materialize the formerly incremental
* anchor plane and sort every fringe dimension since pass A now uses rank order.
* Purpose: Official experiment reducing repeated bitset traffic in eight dimensions.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/121092: retained combined static-fringe/incremental-anchor engine.
* saffah_codex_6a_agg3, https://duck.ac/submission/121062: adapted second-bucket query locality.
* All inherited citations retained; no separate license notices were displayed.
* Idea: Within each existing anchor run, order queries by a second plane's bucket
* using one extra stable count-sort pass, improving neighboring-query bitset reuse.
* Purpose: Official experiment reducing random bitset reads in the eight-dimensional engine.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/120953:
* retained its accepted grid32 engine and inherited citations.
* saffah_codex_6a_agg3, https://duck.ac/submission/120922:
* reused the exact shared incremental-anchor construction.
* saffah_codex_6a_agg3, https://duck.ac/submission/121022:
* reused compile-time fringe dimensions (that submitted variant used grid40).
* Public sources displayed no separate software license notices.
* Idea: Combine the independently accepted anchor-traffic and dimension-constant
* improvements at grid coefficient32. Preserve exact queries and final query order.
* Purpose: Official experiment combining complementary construction/scan savings.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/120953:
* retained its accepted eight-dimensional grid40 engine and inherited credits.
* saffah_codex_6a_agg3, https://duck.ac/submission/120949:
* adapted compile-time fringe dimension specialization from the large task.
* Public sources displayed no separate software license notices.
* Idea: Specialize each outer fringe dimension so signature plane selection and
* exclusion masks can be compiled once instead of recomputed in the scan.
* Preserve the final query-order restoration after all dimensions finish.
* Purpose: Official experiment reducing invariant scan instructions.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/120919:
* retained the accepted engine and inherited citations.
* Public source displayed no separate software license notice.
* Idea: Change grid coefficient from36 to40 after the new cache/codegen
* improvements, balancing dense bitset reads against signature fringe lengths.
* All ranks, boundary strictness and exact ambiguity checks remain in place.
* Purpose: Official experiment locating the new geometry crossover.
*/
// ===== REFERENCES（a14y_ 席，2026-10-02）=====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 **#120507** <https://duck.ac/submission/120507>
//     （**566.842498 ms** = 本题现役 T）。本件 = 该提交正文 + 本席注释；其相对我方 #120336 的**真改动 4 处**全数采用
//     （后两处是 #120507 在 #120475 之上新增的）：(d) 热环 `qr` 在 `KK==8` 时去掉 `| khi8v`（k==8 时 khi8v 全 0 ⇒ 可证无用）；
//     (e) 建表用**一条 8-lane gather + 一条 32 B store** 取代 8 次标量 BVAL 存取（d>=k 的车道随后被 0x7FFFFFFF 覆盖，
//     索引上限 255+7*256=2047 < 8*256 ⇒ 无越界）。逐位等价由 15 组对拍实测确认。
// [1b] duck.ac 用户 saffah_codex_6a_agg3，提交 **#120475** <https://duck.ac/submission/120475>
//     （**567.910149 ms** = 本题现役 T，2026-10-02 18:55:08 提交）。本件 = 该提交正文 + 本席的一处注释/闸门补写；
//     **采用该提交的内容**（其头部自述 "Idea: Cache the eight bucket-bound coordinates in the same 64-byte record
//     as the original eight coordinates ... replaces its eight-lane gather with the adjacent vector load"）：
//       (a) `#define STRIDE 8 -> 16`（每点记录 32 B -> 64 B = 恰好一条 cache line）；
//       (b) 建表循环写 `p[8+d] = BVAL(d, dst[d])`（d<k；d>=k 写 0x7FFFFFFF，与旧 blendv 同值）；
//       (c) 热环把 8-lane `_mm256_i32gather_epi32(s_bval, vidx, 4)` 换成同一条线内的 `_mm256_load_si256(rec+8)`。
//     **逐位等价论证**：被删 gather 的 lane e 取值 `BVAL(e, bd_e(i))`，而 `bd` 就是 `s_bd + i*MAXK`、
//     即建表时写下的 `dst[e]`（同源同值），车道 e>=k 两式都取 0x7FFFFFFF ⇒ 语义逐位不变。
//     （本席 15 组对拍实测 stdout 全部相同，见下。）
// [2] duck.ac 用户 saffah_codex_6a_agg3，提交 **#120231** <https://duck.ac/submission/120231>
//     （579.986914 ms；本件与 [1] 的共同正文基底，其头部自述基底 = 本账号 #114076）。
// [3] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#120336** <https://duck.ac/submission/120336>
//     （**573.085029 ms** = 本席开工时的 mine；[1] 的正文基底即此件**逐字节**，本件即在该件正文上叠 [1] 的两处真改动）。
// [4] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；原提交无独立许可证声明，按 RULES §3 署名。
// ======================
// ===== 本次思路（a14y_ 席）=====
// 【分解：把 T 件 #120475 拆成"真改动 / 落后量 / 有害"】逐行 diff(#120475, #120231) 与 diff(#120475, 我方 #120336)：
//   · 真改动 2 处 = 上文 (a)(b)(c)（`STRIDE 8->16` 与记录内 Avraw 缓存，后者是其自述"参考我方 #114343 的 Avraw 预计算"的移植）；
//   · 落后量 = **0**（它的基底就是我方最新件）；有害改动 = **0**；
//   · 判题机读数：它 567.910149 vs 我方 #120336 573.085029 ⇒ 该真改动值 **−5.175 ms（−0.903%）**。
// 【本发】= 我方 #120336 正文 + 该两处真改动（不加任何未定价项）⇒ 目标 = 追平 T（567.91），再由图外刀法拉开。
// 【本席另两条实测，写给后席（装置与读数都在 work/a14y_tmp/）】
//   (a) 本机相位账 + 消融（同一形状 k=8/n=3e5/mode=10；`g_ph4`/`g_noA`/`g_noB`）：**pass A
//       (`and_popcount_range`：池位图 AND+popcount) ≈ 45%（消融 1050→573 ms）· pass B（fringe 签名扫描）≈ 30%
//       （→746 ms）· pool 构建 ≈ 11% · orders ≈ 8%** ⇒ ★ 本题**第一大项是 pass A**。
//   (b) ★ 与在档旧读数（`a14z_cand.md`：and_popcount_range 只占整题 0.94%）矛盾的根因 = **规模**：那次量的是小 n 形状
//       （整题 146 ms 级，池 ≈6 MB ⇒ L3 驻留）；n=3e5 时池 = K1*(Q+1)*W 词 ≈ **60 MB**，而 pass A 每 query 读
//       `K1 * (R>>6)` 词 ≈ 29 KB ⇒ 它从"算"变成"**访存流**"（本机 ≈19 GB/s 有效）。⇒ 后续刀应打在
//       **每 query 的 K1×R/64 词读取**或**池的 60 MB 体积**上，而不是 and_popcount_range 的指令数。
//   (c) `MAXQ` 轴已到底：`solve_mf` 的 `if (g_Q > 255u) g_Q = 255u;` 是无条件钳位（先于所有 MAXQ 使用），
//       `BVAL(d,Q)` 在 b=Q=255 有哨兵写 ⇒ MAXQ=256 是**下界**，128/64 会越界（不可作刀）。
// 闸门：判题同款 `ref/gcc9/g9.sh -O2 -std=c++17 -static -U_FORTIFY_SOURCE`；本件与 #120336 原版**双载 -DLOCAL_TEST
//   15 组逐位对拍**（n=3e5 × mode{10,11,13} × Q{自主,200,228,255} + n=2000 × mode{0..4} + k=7/10）
//   stdout **全部相同**，**生产 checksum 351867771 不变**；(241) 正控走输出命中。
// ======================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/114076:
* 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 1014a.
* 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，提交 **#114041** <https://duck.ac/submission/114041>
//     (**590.143241 ms**) —— 本文件正文 = 该提交正文**逐字节** + 下方"思路"的唯一改动 (该件含 `no-tree-reassoc`,
//     在本题实测 -723.851 us)。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 **#109857** <https://duck.ac/submission/109857>
//     (**595.814870 ms** = 本题 T) => 严支 `0.99*T+1us = 589.857721 ms`。
// [3] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；原提交无独立许可证声明，按 RULES §3 署名。
// ======================
// ===== 思路 =====
// 【本发只有一个变量：在含 `no-tree-reassoc` 的 optimize 串尾**再叠一项** `live-range-shrinkage`】
//   = 叠件 (协调者派单③: 与已晋级的 `no-tree-reassoc` 叠成一件)。
//   选它的理由：它属于"纯扰动池" —— callgrind 实测 Ir 变化 -0.11% (不改工作量)，而 `.text` 真改机器码。
//   ★ (936) 第 (3) 条已核对：**Ir 可加 ≠ 时间可加**，故本件单独真发取判题机读数，不作纸面外推。
// 闸门：判题同款 `ref/gcc9/g9.sh -O2 -std=c++17` (judge_flags.txt 缺 => 未新建) rc=0；
//   双载 -DLOCAL_TEST stdout 逐位同 (5 形状含 mode>=10)；(241) 正控走输出。
// ======================
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#112420** <https://duck.ac/submission/112420>
//     (**590.848619 ms** = 本题现役 mine) —— 本文件正文 = 该提交正文**逐字节**。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 **#109857** <https://duck.ac/submission/109857>
//     (**595.814870 ms** = 本题 T) => 严支 `0.99*T+1us = 589.857721 ms`。
// [3] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；原提交无独立许可证声明，按 RULES §3 署名。
// ======================
// ===== 思路 =====
// 【本发只有一个变量：optimize 串尾追加 `no-tree-reassoc`】
//   依据 = 1014 上 `no-tree-fre` 换来 -1.443201 ms 且**在热环里一条指令都没删** (只改归约链结合次序 + 寄存器映射)
//   => "扰动彩票" (改 SSA/IRA/调度而不改工作量) 是真实收益来源。`no-tree-reassoc` 直击该机制本体。
//   ★ 同批实测 (本席 a14y_ 第一轮, 同窗 4 秒内): 以本批自身基线臂为参照, **`no-tree-reassoc` = -768 us**
//     (同批 `web` +206 us / `no-tree-pre` +116 us / `tracer`族无 / `no-tree-fre` +5575 us / `rename-registers` +2591 us)
//   => 它是本批唯一的 <= -0.3 ms 项, 按协调者派单③ 的晋级规则立刻单独复测定价。
//   ★ 本机 Ir 中性: callgrind 1,210,524 -> 1,211,498 (+0.08%), 即"纯扰动、不改工作量" ✓
// 闸门：判题同款 `ref/gcc9/g9.sh -O2 -std=c++17` rc=0；双载 -DLOCAL_TEST stdout 与基座逐位同 (5 形状含 mode>=10)；
//   (241) 正控走输出。**不重开预取轴 / s_sdmp 消费轴 / 页故障轴**。
// ======================
// [trim_] identical=True ; base_lines=1502 trimmed_lines=1502 ; base=103293 B trimmed=87206 B
// [trim_] 历史注释预压（正文一字节未动；REFERENCES 块与全部 duck.ac URL 保留）
//   原 103293 B → 预压后 86947 B（去 202 行注释 / 16346 B）；机器证明：剥注释后代码行逐行相同 = True
// ======================
//     The copied engine inherited our signature transpose.
// Approach:
// Set bucket-count coefficient to 36 in the cited engine to trade pool build
// traffic against fringe candidate scan traffic.
// Purpose:
// Measure the new bucket point on duck.ac official 1014a test workload.
#define MF_QCOEF 32
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #108836 <https://duck.ac/submission/108836>
//     （612.869165 ms）与 #108803 <https://duck.ac/submission/108803>：**本文件正文基底**
//     —— 直接复制其整份公开正文（他自述基底为我方 #108160；继承引用链见其头部原样保留）。
//     合规：其正文**无任何许可证声明**；duck.ac 提交正文按站点规则公开可见、可直接取用；
//     本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名引用账号与原提交地址。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号团队）提交 #108160
//     <https://duck.ac/submission/108160>：经 [1] 继承的正文基底。
// ======================
// ===== 思路 =====
// register pair, skipping the excluded dimension before the final permute.
// This avoids the intermediate column array and its stores and reloads.
// Purpose:
// Measure whether compiler option set plain speeds our cited base.
// ===== REFERENCES（a14w_ 席，2026-10-02）=====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#120554** <https://duck.ac/submission/120554>
//     （**566.851569 ms** = 本席开工时的 mine）。**本文件正文 = 该提交正文逐字节**，唯一改动见下。
// [2] duck.ac 用户 saffah_codex_6a_agg3，提交 **#120507** <https://duck.ac/submission/120507>
//     （**566.842498 ms** = 本题现役 T）/#120475 <https://duck.ac/submission/120475>（567.910149 ms）
//     /#120231 <https://duck.ac/submission/120231>（579.986914 ms）：#120554 正文中的对手改动的来源。
// [3] duck.ac 用户 saffah_codex_6s_agg2，提交 #105554 <https://duck.ac/submission/105554> 与
//     #104600 <https://duck.ac/submission/104600>：pass-A/pass-B counting-sort bucket-key cache 与
//     scatter load-chase 预取两把刀的来源（正文内已逐处署名）。
// [4] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；原提交无独立许可证声明，
//     按 /home/yjp/duck.ac/RULES.md 第 3 节署名引用账号与原提交地址。本件继承引用链见正文头部。
// ======================
// ===== 本次思路（a14w_ 席，2026-10-02）=====
// 【本发只有一个变量】在最后一处 `#pragma GCC optimize(...)` 串尾**只追加一项 `no-tree-ch`**
//   （= 关掉 loop header copying / 循环头复制趟）。正文其余部分与我的 #120554 **逐字节相同**。
// 【为什么是它】本席先做了判题机一侧的分相计时（自建 MF_PHASE 探针装置）：
//   真规模 n=3e5/k=8 下 **pass A `and_popcount_range` = 942 Mcyc（46.0%）· pass B 签名扫描 = 723 Mcyc
//   （35.4%）· 池构建 = 167.6 M（8.2%）· orders = 122 M（6.0%）**，与在档相位账一致。
//   随后逐项证否了三条大轴 —— (i) **Q/池体积两向**（Q 228→152 判题机实测：paand −18.9% 但 pb +33.8%，
//   净 **+2.6%**；Q→255 在档 +0.162%）；(ii) **一切预取形态**（T0 内环预取 12→6→0 分别 **+0.93%/+1.63%**；
//   改 T1 提示 **+5.6%**；twosort 复用 **−0.04%** ⇒ 全部判负，与在档"本引擎预取必负"一致）；
//   (iii) **一次改写闸门**：pass A 实测 **8.75 GB / 942 Mcyc = 137 M 条 L1-miss line，12 个 LFB ⇒
//   等效延迟 82 cyc**，与 L3 延迟同量级 ⇒ **pass A 是 L1-miss/LFB 吞吐墙（不是 DRAM 带宽墙）**，
//   而 7 个平面 × R/64 词的读取量与重用结构均已最优（`R = min_d rank` 已最优、7 平面是结构、
//   逐词 AND 后非零率 90% ⇒ 无剪枝）⇒ **pass A 无剩余空间**。
//   ⇒ 于是回到"扰动彩票"轴（不改工作量、只改 SSA/IRA/调度/代码布局），用**判题机侧 rdtsc 分相探针**定价
//   （同文件重复性实测 **0.016%**，两次独立读数 **−1.066% / −1.065%**）。本批 13 个真改机器码项里
//   `no-tree-ch` 是唯一 ≤ −0.3 ms 者（同批 no-tree-dominator-opts +0.5% · no-tree-ter +0.8% ·
//   no-tree-dce +0.6% · no-tree-vrp +1.0% · no-schedule-insns2 +2.6% · no-ira-share-spill-slots +0.14% ·
//   no-peephole2 +0.09% 全为负）。
// 【分解出的机理】**pass B 签名扫描桶 723.0 → 701.8 M（−2.9%）**、pass A 基本不动（941.4 → 940.1 M）
//   ⇒ `no-tree-ch` 消掉的正是签名扫描那 7 条平面流内环的**循环头复制**（同一循环体被复制成"首轮 + 稳态"
//   两份，I-cache/分支预测双双受压）。
// 【闸门】判题同款 `ref/gcc9/g9.sh -O2 -std=c++17 -static -U_FORTIFY_SOURCE` rc=0（零告警）；
//   与 #120554 原版**双载 -DLOCAL_TEST 15 组逐位对拍**（n=3e5 × mode{10,11,13} × Q{自主,200,228,255}
//   + n=2000 × mode{0..4} + k=7/10）stdout **全部逐位相同**，生产 checksum 不变；(241) 正控走输出命中。
// ======================
// ===== 追加（a14w_ 席第三发，2026-10-02）=====
// 【相对本席第一发 #120788（561.367301 ms，mine）的唯一变量】同在 optimize 串尾再追加两项
//   `unroll-all-loops`、`no-sched-critical-path-heuristic`；正文其余与 #120788 逐字节相同。
// 【依据】#120788 后缺口只剩 0.192228 ms（0.034%），故用判题机侧 rdtsc 分相探针在 no-tree-ch 基座上
//   继续扫剩余真改机器码项（同文件重复性 0.016%）。本组合三次独立读数 TOTAL 2022.5~2022.8 →
//   **2012.53 / 2012.55 / 2012.82 Mcyc（−0.49%）**；其两个分量 `unroll-all-loops` 单独
//   2017.85/2017.51（−0.24%）、`no-sched-critical-path-heuristic` 单独 2018.39/2018.96/2018.78（−0.19%）。
//   同批其它全为负或持平：no-gcse(−0.18% 但与上二者相叠反而 2020.24) · no-tree-ter −0.05% ·
//   no-early-inlining −0.05% · no-caller-saves/no-if-conversion/no-optimize-sibling-calls 0 ·
//   no-tree-loop-im/no-tree-pta ±0 · no-tree-ccp +0.10% · schedule-insns +0.35% · align-loops=32 +0.55% ·
//   ira-region=one +0.55% · no-guess-branch-probability +0.78% · twosort +0.05%。
// 【闸门】判题同款 flag rc=0；与 #120554 原版**双载 -DLOCAL_TEST 15 组逐位对拍** stdout 全同。
// ======================
#pragma GCC optimize("O3,no-strict-aliasing","no-tree-reassoc","live-range-shrinkage","no-tree-forwprop","no-tree-ch","unroll-all-loops","no-sched-critical-path-heuristic")
#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>
#include <type_traits>
#include <type_traits>
#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; })
// [y14q16] u16 bucket row -> 8 x u32 zero-extended lanes.  Produces byte-for-byte the same
// eight values the old 256-bit row load produced (every stored bucket is clamped to <= Q-1),
// from one 128-bit load; asm form mirrors LDU256 so it needs no target attribute.
#define AND256M(a,p) ({ __m256i _v; __asm__("vpand %2, %1, %0" : "=x"(_v) : "x"(a), "m"(*(const __m256i *)(p))); _v; })
#define LDU8U16(p) ({ __m256i _v; \
__asm__("vpmovzxbd %1, %0" : "=x"(_v) : "m"(*(const u64 *)(p))); _v; })
#ifndef MAXN
#define MAXN 1000005
#endif
#define MAXK 8
#define STRIDE 16
// Bucket-count coefficient: q ~ sqrt(n) * MF_QCOEF/100.  Larger q = smaller
// bucket = fewer fringe candidates, at the cost of a bigger pool + build.
// Bucket count = sqrt(n) * MF_QCOEF/100 * kscale[k].  MEASURED OPTIMA (board):
//   n=1e6  k=6..8 : 122 default is fine (1013 7.491 s, 1014 8.503 s at 122)
//   n=1e6  k=9/10 : 92  (1016 10.888 s at 92; 10.959 at 122; 11.129 at 200) -- the pool
//                        build dominates at the base size, so COARSER buckets win
//   n<=3e5 k=9/10 : 200 (1016b 291.4 ms vs 326.6 at 122; 1015a 1.444, 1016a 1.517) -- the
//                        pool is not memory-capped there, so FINER buckets shrink the fringe
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 42
#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
#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 128
#endif
#ifndef MF_PFB2
#define MF_PFB2 128
#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 256
#endif
#define MF_SIGPL ((((size_t)MAXN + 63u) & ~(size_t)63u) + 0x840u)
static u32 s_ord[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u8 s_bd[(size_t)MAXN * MAXK];   // [y14q16] u16: the only write site clamps to <= Q-1 (Q=354 => <=353) => lossless; table 32 MB -> 16 MB
static u8 s_bdv[(size_t)MAXN * MAXK];   // [x14z/d4b-split] dim-major staging plane
static u32 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_bandrec[(size_t)MAXN*(MAXK+1)],s_bandacc[MAXN];
static __m128i s_bandvec[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 u32 *inc_anchor=0;static size_t inc_cap=0;static bool inc_on=false;
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
#define MAXRS 256u   // [a14c] s_maxr row stride (b <= Q <= 255)
static u32 s_maxr[(size_t)MAXK * MAXK * MAXRS];
// ---- 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)]
extern "C" int madvise(void *, unsigned long, int);
#define MF_MADV_HUGEPAGE 14
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
void *p = (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
/* [a14e-THP] the big pools are freshly mmapped and written once front-to-back, so every
4 KB page pays a minor fault + kernel zeroing (~14.5k faults before sdmp, ~19k for the
sdmp pool alone; 34 k total ~= 3.6% of the run).  MADV_HUGEPAGE gives 2 MB pages: 512x
fewer traps, same zeroing bytes.  Data layout and results are untouched. */
madvise(p, bytes, MF_MADV_HUGEPAGE);
return p;
}
__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) {
{ __asm__ __volatile__("nop;nop;nop;nop;nop;nop" :: "m"(*(const char(*)[32])bs[0]) : "ymm15"); }
// 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 + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[1] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 40 + 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 (!g_skip_and) {
const u32 nq = q1 - q0;
u32 *rec = s_bandrec;
{
u32 *cv = s_cntv;
for (u32 v = 0; v <= (u32)N; v++) cv[v] = 0;
for (u32 t = 0; t < nq; t++) cv[s_gt[ds][s_qord[q0 + t]] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)N; v++) { u32 c = cv[v]; cv[v] = ac; ac += c; } }
for (u32 t = 0; t < nq; t++) {
u32 i = s_qord[q0 + t];
const u8 *bd = s_bd + (size_t)i * MAXK;
u32 R = s_gt[ds][i];
u32 *rp = rec + (size_t)(cv[R + 1]++) * (KK+1);
rp[0] = R;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
rp[1 + e] = g_plen_on ? s_pofs[(size_t)e * MAXQ + bd[d]]
: (u32)(((size_t)e * (Q + 1) + bd[d]) * W);
}
rp[KK] = i;
}
}
{
u32 *acc = s_bandacc;
for (u32 t = 0; t < nq; t++){ acc[t] = 0; s_bandvec[t]=_mm_setzero_si128(); }
if (nq) {
const u32 maxR = rec[(size_t)(nq - 1) * (KK+1)];
const u32 maxband = (maxR + 1023u) >> 10;
const __m256i lm = _mm256_set1_epi8(0x0f);
const __m256i lk = _mm256_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);
const __m256i zz = _mm256_setzero_si256();
u32 ps = 0;
for (u32 bk = 0; bk <= maxband; bk++) {
const u32 lim = bk << 10;
while (ps < nq && rec[(size_t)ps * (KK+1)] <= lim) ps++;
if (ps >= nq) break;
u32 psp = ps;
while (psp < nq && rec[(size_t)psp * (KK+1)] < lim + 1024u) psp++;
const size_t off8 = (size_t)bk << 4;
for (u32 t = psp; t < nq; t++) {
const u32 *rp = rec + (size_t)t * (KK+1);
__m256i bytes=zz;
for (u32 lane=0;lane<16;lane+=8) {
const u64 *pz = s_bs + off8 + lane;
__m256i a0 = LDU256(pz + rp[1]), a1 = LDU256(pz + rp[1] + 4);
for (int e = 1; e < K1; e++) {
const u64 *pe = pz + rp[1 + e];
a0 = AND256M(a0, pe);
a1 = AND256M(a1, pe + 4);
}
__m256i lo0 = _mm256_and_si256(a0, lm), hi0 = _mm256_and_si256(_mm256_srli_epi16(a0, 4), lm);
__m256i lo1 = _mm256_and_si256(a1, lm), hi1 = _mm256_and_si256(_mm256_srli_epi16(a1, 4), lm);
__m256i s0 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0));
__m256i s1 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1));
bytes = _mm256_add_epi8(bytes,_mm256_add_epi8(s0,s1));
}
__m256i sum64=_mm256_sad_epu8(bytes,zz);
__m128i lanes=_mm_add_epi32(_mm256_castsi256_si128(sum64),_mm256_extracti128_si256(sum64,1));
s_bandvec[t]=_mm_add_epi32(s_bandvec[t],lanes);
}
for (u32 t = ps; t < psp; t++) {
const u32 *rp = rec + (size_t)t * (KK+1);
const u32 R = rp[0];
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) bp[e] = s_bs + (size_t)rp[1 + e] + off8;
acc[t] += and_popcount_range<KK - 1>(bp, (R >> 6) - (u32)off8, R & 63u);
}
}
for (u32 t = 0; t < nq; t++){__m128i lanes=s_bandvec[t];out[rec[(size_t)t * (KK+1) + KK]] = acc[t]+(u32)_mm_cvtsi128_si32(lanes)+(u32)_mm_extract_epi32(lanes,2);}
}
}
}
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. ----
auto fringe_dim=[&](auto dim){ constexpr int d=decltype(dim)::value;
if (d == ds) return;
#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 (true)
#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.
// [x14z DISP32] no plane pointer table: dim d's k-1 planes live at t*MF_SIGPL,
// t = 0..k-2 in scan order (e = t<d ? t : t+1), so every plane load in the ring is
// `disp32(%pbase_sig,%rsi)` off ONE base register instead of a stack pointer reload.
const u8 *const pbase_sig = s_sig + (size_t)d * (k - 1) * MF_SIGPL;
// [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 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];
__asm__ __volatile__("prefetchw %0" ::"m"(*(const char *)(&out[j])));   // [y14q] write-intent (RMW target; was _mm_prefetch(out+j, T0))
_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);
}
}
#endif
const u8 *bd = s_bd + (size_t)i * MAXK;   // [y14q16] u16 row
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.
// [KK8] for KK==8 every lane of khi8v is 0 (lo8[e] = (e<k) ? 0 : 0x7FFFFFFF), so the
// OR is a provable no-op there; keep it for KK<8 where lanes >= k need the 0x7FFFFFFF.
__m256i qr = (KK == 8) ? _mm256_load_si256((const __m256i *)rec)
: _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.
// [AVC] Avraw is materialised in the record itself (rec[8..15]): the value
// BVAL(e, bd_e(i)) depends ONLY on (i,e) and is unchanged during the scan, so
// the eight-lane gather is a redundant recomputation of a stored constant.
__m256i Avraw = _mm256_load_si256((const __m256i *)(rec + 8));
__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 (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);
if (k > 8) sv = _mm256_or_si256(sv, _mm256_cmpeq_epi32(qT, maxv));
if (k == 8) sv = _mm256_subs_epu8(sv, _mm256_set1_epi8(1));
const __m256i lohalf = _mm256_broadcastsi128_si256(_mm256_castsi256_si128(sv));
const __m256i hihalf = _mm256_broadcastsi128_si256(_mm256_extracti128_si256(sv, 1));
if (k == 8) {
switch (d) {
case 0:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 1:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 2:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 3:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 4:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 5:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 6:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 7:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
break;
}
} else {
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));
}
}
}
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 ac0 = _mm256_setzero_si256(), ac1 = ac0;
for (int t = 0; t < KK - 1; t++) {
__m256i cv = LDU256((const u8 *)(pbase_sig + (size_t)t * MF_SIGPL + aj));
__m256i df = _mm256_subs_epu8(cv, qsv[t]);
if (t & 1) ac1 = _mm256_or_si256(ac1, df);
else ac0 = _mm256_or_si256(ac0, df);
}
__m256i ac = _mm256_or_si256(ac0, ac1);
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 ? (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) {
const u32 base = aj, lk = 0u;
u32 kmask = (1u << r) - 1u;
__m256i ac0 = _mm256_setzero_si256(), ac1 = ac0;
for (int t = 0; t < KK - 1; t++) {
__m256i cv = LDU256((const u8 *)(pbase_sig + (size_t)t * MF_SIGPL + base));
__m256i df = _mm256_subs_epu8(cv, qsv[t]);
if (t & 1) ac1 = _mm256_or_si256(ac1, df);
else ac0 = _mm256_or_si256(ac0, df);
}
__m256i ac = _mm256_or_si256(ac0, ac1);
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 ? (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 constexpr(KK>0) fringe_dim(std::integral_constant<int,0>{});
if constexpr(KK>1) fringe_dim(std::integral_constant<int,1>{});
if constexpr(KK>2) fringe_dim(std::integral_constant<int,2>{});
if constexpr(KK>3) fringe_dim(std::integral_constant<int,3>{});
if constexpr(KK>4) fringe_dim(std::integral_constant<int,4>{});
if constexpr(KK>5) fringe_dim(std::integral_constant<int,5>{});
if constexpr(KK>6) fringe_dim(std::integral_constant<int,6>{});
if constexpr(KK>7) fringe_dim(std::integral_constant<int,7>{});
if constexpr(KK>8) fringe_dim(std::integral_constant<int,8>{});
if constexpr(KK>9) fringe_dim(std::integral_constant<int,9>{});
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 - 1) * (u64)MF_SIGPL;   // [x14z DISP32] k*(k-1) planes
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;
}
if (g_Q > 255u) { g_Q = 255u; g_B = (N + g_Q - 1) / g_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
// Apply the existing no-filter policy before the per-dimension order pass.
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).
static u8 hist8[MAXN+1] __attribute__((aligned(64)));
for(u32 i=0;i<=N;i++)hist8[i]=0;
for(u32 i=0;i<N;i++){u32 v=xd[i];hist8[v<N?v:0]++;}
u32 sum=0;
for(u32 i=0;i<=N;i++){u32 c=hist8[i];s_tmp[i]=sum;sum+=c;}
bool compact=(sum==N);
if(!compact){
  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]++;}
  sum=0;for(u32 i=0;i<=N;i++){u32 c=s_cntv[i];s_cntv[i]=sum;s_tmp[i]=sum;sum+=c;}
}
u32 *gt=s_gt[d],*pp=s_pos[d];
const u64 BMAG=(((u64)1<<41)/B)+1u;
auto scatter=[&](auto tag){
  for(u32 i=0;i<N;i++){
    u32 v=xd[i];if(v>=N)v=0;
    u32 st,q,num;
    if constexpr(decltype(tag)::value){
      u32 w=s_tmp[v];st=w&0xfffffu;q=st+(w>>20);s_tmp[v]=w+(1u<<20);
      num=(s_tmp[v+1]&0xfffffu)-1u;
    }else{st=s_cntv[v];q=s_tmp[v]++;num=s_cntv[v+1]-1u;}
    od[q]=i;pp[i]=q;gt[i]=st;
#if MF_FB2
    s_posp[(size_t)i*MAXK+d]=q;
#endif
    u32 bb=(u32)(((u64)num*BMAG)>>41);
    s_bdv[(size_t)d*N+i]=bb>Q-1?Q-1:bb;
  }
};
if(compact)scatter(std::true_type{});else scatter(std::false_type{});
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;
}
// [x14z/d4b-split] one dense transpose: dim-major -> record layout.  The per-dim
// pass wrote only 2 B per 16-B record (8/64 B of a line) across 8 passes, so it
// touched 8*N/4 lines with RFO; the staging plane is dense and this single dense
// transpose replaces that line traffic 1:1.
// [a14c FUSE2] the two dense staging sweeps below write disjoint destinations from
// disjoint sources, so they share one pass over i.
for (u32 i = 0; i < N; i++) {
u8 *dst = s_bd + (size_t)i * MAXK;
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) { dst[d] = s_bdv[(size_t)d * N + i]; p[d] = x[d][i]; }
// [AVC] same record, adjacent 32 B: Avraw[e] = BVAL(e, bd_e(i)); lanes e >= k stay
// at the 0x7FFFFFFF the masked gather+blend produced, so every consumer is bit-equal.
{
// [AVB1] one 8-lane gather + one 32 B store replaces eight scalar BVAL loads/stores.
// Index = (bucket id) + d*MAXQ, dst[d] <= 255 so idx <= 2047 < (MAXK+1)*MAXQ: always in
// bounds even for the stale dst[d], d >= k, which the mask below overwrites anyway.
__m256i idx = _mm256_add_epi32(_mm256_cvtepu8_epi32(_mm_loadl_epi64((const __m128i *)dst)),
_mm256_setr_epi32(0, MAXQ, 2*MAXQ, 3*MAXQ,
4*MAXQ, 5*MAXQ, 6*MAXQ, 7*MAXQ));
_mm256_store_si256((__m256i *)(p + 8),
_mm256_i32gather_epi32((const int *)s_bval, idx, 4));
if (k < 8) for (int d = k; d < 8; d++) p[8 + d] = 0x7fffffffu;
}
}
// ---------- 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 - 1) * (size_t)MF_SIGPL;   // [x14z DISP32]
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;
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 * MAXRS;
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 u8 *bdv = s_bd + (size_t)i * MAXK;   // [y14q16] u16 row
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXRS;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXRS + b]) mr[(size_t)d * MAXRS + 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) {
u8 *sb = s_sig ? (s_sig + (size_t)d * (k - 1) * MF_SIGPL) : 0;
u32 j = 0;
if (sb && k == 8) {
const __m256i sigmul = _mm256_set1_epi32((int)g_sigmul);
const __m256i one = _mm256_set1_epi32(1);
for (; j + 8 <= N; j += 8) {
// Signature row 0 of the eight-point block.
#if MF_PFB > 0
if (j + 0 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 0 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row0 = s_pm + (size_t)od[j + 0] * STRIDE;
__m256i rv0 = _mm256_loadu_si256((const __m256i *)row0);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 0) * KSx), rv0);
__m256i r0 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv0, sigmul), 23), one);
// Signature row 1 of the eight-point block.
#if MF_PFB > 0
if (j + 1 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 1 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row1 = s_pm + (size_t)od[j + 1] * STRIDE;
__m256i rv1 = _mm256_loadu_si256((const __m256i *)row1);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 1) * KSx), rv1);
__m256i r1 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv1, sigmul), 23), one);
// Signature row 2 of the eight-point block.
#if MF_PFB > 0
if (j + 2 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 2 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row2 = s_pm + (size_t)od[j + 2] * STRIDE;
__m256i rv2 = _mm256_loadu_si256((const __m256i *)row2);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 2) * KSx), rv2);
__m256i r2 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv2, sigmul), 23), one);
// Signature row 3 of the eight-point block.
#if MF_PFB > 0
if (j + 3 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 3 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row3 = s_pm + (size_t)od[j + 3] * STRIDE;
__m256i rv3 = _mm256_loadu_si256((const __m256i *)row3);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 3) * KSx), rv3);
__m256i r3 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv3, sigmul), 23), one);
// Signature row 4 of the eight-point block.
#if MF_PFB > 0
if (j + 4 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 4 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row4 = s_pm + (size_t)od[j + 4] * STRIDE;
__m256i rv4 = _mm256_loadu_si256((const __m256i *)row4);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 4) * KSx), rv4);
__m256i r4 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv4, sigmul), 23), one);
// Signature row 5 of the eight-point block.
#if MF_PFB > 0
if (j + 5 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 5 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row5 = s_pm + (size_t)od[j + 5] * STRIDE;
__m256i rv5 = _mm256_loadu_si256((const __m256i *)row5);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 5) * KSx), rv5);
__m256i r5 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv5, sigmul), 23), one);
// Signature row 6 of the eight-point block.
#if MF_PFB > 0
if (j + 6 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 6 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row6 = s_pm + (size_t)od[j + 6] * STRIDE;
__m256i rv6 = _mm256_loadu_si256((const __m256i *)row6);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 6) * KSx), rv6);
__m256i r6 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv6, sigmul), 23), one);
// Signature row 7 of the eight-point block.
#if MF_PFB > 0
if (j + 7 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 7 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row7 = s_pm + (size_t)od[j + 7] * STRIDE;
__m256i rv7 = _mm256_loadu_si256((const __m256i *)row7);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 7) * KSx), rv7);
__m256i r7 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv7, sigmul), 23), one);
__m256i t0 = _mm256_unpacklo_epi32(r0, r1);
__m256i t1 = _mm256_unpackhi_epi32(r0, r1);
__m256i t2 = _mm256_unpacklo_epi32(r2, r3);
__m256i t3 = _mm256_unpackhi_epi32(r2, r3);
__m256i t4 = _mm256_unpacklo_epi32(r4, r5);
__m256i t5 = _mm256_unpackhi_epi32(r4, r5);
__m256i t6 = _mm256_unpacklo_epi32(r6, r7);
__m256i t7 = _mm256_unpackhi_epi32(r6, r7);
__m256i u0 = _mm256_unpacklo_epi64(t0, t2);
__m256i u1 = _mm256_unpackhi_epi64(t0, t2);
__m256i u2 = _mm256_unpacklo_epi64(t1, t3);
__m256i u3 = _mm256_unpackhi_epi64(t1, t3);
__m256i u4 = _mm256_unpacklo_epi64(t4, t6);
__m256i u5 = _mm256_unpackhi_epi64(t4, t6);
__m256i u6 = _mm256_unpacklo_epi64(t5, t7);
__m256i u7 = _mm256_unpackhi_epi64(t5, t7);
if ((int)d != 0) {
__m256i v = _mm256_permute2x128_si256(u0, u4, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((0 < (int)d) ? 0 : 0 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 1) {
__m256i v = _mm256_permute2x128_si256(u1, u5, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((1 < (int)d) ? 1 : 1 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 2) {
__m256i v = _mm256_permute2x128_si256(u2, u6, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((2 < (int)d) ? 2 : 2 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 3) {
__m256i v = _mm256_permute2x128_si256(u3, u7, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((3 < (int)d) ? 3 : 3 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 4) {
__m256i v = _mm256_permute2x128_si256(u0, u4, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((4 < (int)d) ? 4 : 4 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 5) {
__m256i v = _mm256_permute2x128_si256(u1, u5, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((5 < (int)d) ? 5 : 5 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 6) {
__m256i v = _mm256_permute2x128_si256(u2, u6, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((6 < (int)d) ? 6 : 6 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 7) {
__m256i v = _mm256_permute2x128_si256(u3, u7, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((7 < (int)d) ? 7 : 7 - 1) * MF_SIGPL + j,
&bytes, 8);
}
}
}
for (; 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;
__m256i rv = _mm256_loadu_si256((const __m256i *)row);
_mm256_stream_si256((__m256i *)dp, rv);
if (sb) {
/* [y14x] store-side signature, vectorised.  sig_nz(v) = max((v*M)>>23, 1)
and (N-1)*M <= 255*2^23+N < 2^31 (the engine's own invariant), so the
32-bit product is exact: one vpmulld + vpsrld + vpmaxud replaces k
scalar 64-bit multiplies per (dim, point).  The bytes are read back
from a store-forwarded 32-byte slot. */
__m256i sg = _mm256_srli_epi32(
_mm256_mullo_epi32(rv, _mm256_set1_epi32((int)g_sigmul)), 23);
sg = _mm256_max_epu32(sg, _mm256_set1_epi32(1));
u32 tmpk[MAXK];
_mm256_storeu_si256((__m256i *)tmpk, sg);
for (int e = 0; e < k; e++) { if (e != (int)d) sb[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + j] = (u8)tmpk[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 - 1) * MF_SIGPL;
for (int e = 0; e < k; e++) { if (e == (int)d) continue; sb2[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + 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 && k == 8) {
static u8 cached_signature[(size_t)MAXN * 8] __attribute__((aligned(64)));
for (u32 i = 0; i < N; i++)
for (int e = 0; e < 8; e++) {
u32 v = (u32)sig_nz(x[e][i]);
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-1)*MF_SIGPL;
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);
if(d!=0) _mm_storel_epi64((__m128i *)(base+(size_t)(0 < d ? 0 : 0-1)*MF_SIGPL+j),v0);
if(d!=1) _mm_storel_epi64((__m128i *)(base+(size_t)(1 < d ? 1 : 1-1)*MF_SIGPL+j),_mm_srli_si128(v0,8));
if(d!=2) _mm_storel_epi64((__m128i *)(base+(size_t)(2 < d ? 2 : 2-1)*MF_SIGPL+j),v1);
if(d!=3) _mm_storel_epi64((__m128i *)(base+(size_t)(3 < d ? 3 : 3-1)*MF_SIGPL+j),_mm_srli_si128(v1,8));
if(d!=4) _mm_storel_epi64((__m128i *)(base+(size_t)(4 < d ? 4 : 4-1)*MF_SIGPL+j),v2);
if(d!=5) _mm_storel_epi64((__m128i *)(base+(size_t)(5 < d ? 5 : 5-1)*MF_SIGPL+j),_mm_srli_si128(v2,8));
if(d!=6) _mm_storel_epi64((__m128i *)(base+(size_t)(6 < d ? 6 : 6-1)*MF_SIGPL+j),v3);
if(d!=7) _mm_storel_epi64((__m128i *)(base+(size_t)(7 < d ? 7 : 7-1)*MF_SIGPL+j),_mm_srli_si128(v3,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*8;
for(int e=0;e<k;e++)if(e!=d)base[(size_t)(e<d?e:e-1)*MF_SIGPL+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 - 1) * MF_SIGPL;
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 < (int)d) ? e : e - 1) * MF_SIGPL + 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) * MAXRS;
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;
{size_t bytes=(size_t)N*sizeof(u32)+4096;if(bytes>inc_cap){inc_anchor=(u32*)pool_alloc(bytes);inc_cap=inc_anchor?bytes:0;}inc_on=inc_anchor&&(u64)N<=(u64)AS_WORDS*64;}
// ---------- 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) * MAXRS;
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
// [a14c NIDFUSE] bit numbering = rank in dim ds.  s_pos[ds] is exactly that
// inverse permutation, filled once per dim in the order phase; the N-random-store
// rebuild that used to happen here is redundant (probe: 22.34 Mcyc judge-side).
const u32 *const spp = s_pos[ds];
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 + 12 <= pend) {
if (ptr+32 < pend) _mm_prefetch((const char *)(spp + od[ptr+32]), _MM_HINT_T0);
u32 q = spp[od[ptr+0]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+1]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+2]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+3]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+4]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+5]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+6]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+7]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+8]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+9]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+10]];   cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+11]];   cur[q >> 6] |= 1ull << (q & 63);
ptr += 12; }
while (ptr < pend) { u32 q = spp[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
}
query_group<8>(N, out, ds, q0, q1, Q, Wpool);
}
}
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
#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

评测结果