/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/126014: retained validated page-aligned signature planes and seven-dimensional engine, with all inherited citations.
* saffah_cc_v41_agg1, https://duck.ac/submission/123963: adapted compact group-start/arrival fields for counting-sort placement.
* Idea: Keep the existing byte histogram. When its overflow check passes,
* each group has at most255 arrivals, so pack its immutable20-bit start and
* running12-bit arrival count into one32-bit cursor. The unfiltered placement
* pass then uses one random table instead of two, and skips a redundant prefix
* copy. Overflow and filtered cases retain their original exact paths.
* Purpose: Official experiment shrinking random sort-table traffic.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/122903: retained seven-dimensional compact signature planes, blocked transpose, counters and inherited citations.
* saffah_cc_v41_agg1, https://duck.ac/submission/122788: transfer page-aligned signature-plane spacing from the accepted eight-dimensional engine.
* Idea: Pad each compact byte signature plane to a complete 4096-byte page.
* Every logical index and fringe bound is preserved; the changed inter-plane
* addresses alter simultaneous signature stream bank/cache mapping.
* Purpose: Official experiment measuring layout on this six-plane fringe kernel.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#120245**
// <https://duck.ac/submission/120245>（判题机 Time = 69.889315 ms = 本题实时 T）
// —— 所用内容：**本文件正文的基座（逐字节）**。该件 = 本账号 #111503（`b13z_sp`）的副本
// + 它新增的「8 路 SSE 字节转置签名平面构建块」（`if (s_sig && !g_sig_merged && k==7)`）。
// 其自带的引用/思路注释区**原样保留在下方**（含 #111503 → #110262 → #109794 全链）。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#111503**
// <https://duck.ac/submission/111503>（70.664007 ms）= 本账号现役最好件 = #120245 的基座。
// [3] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#120339 / #111492 / #111478**
// <https://duck.ac/submission/120339> 等 —— 本席回看本问题历史件时的对照。
// 合规：duck.ac 提交正文按站点规则公开可见、可直接取用；原提交正文无独立许可证声明。
// 本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名作者账号与原提交地址。
// 除下方「思路」里列出的**三个整数常量**外，正文与 [1] **逐字节相同**（无注释剥离、无空白压缩）。
// ======================
// ===== 思路（[z2-pf]，本席 b13z2_ 的唯一改动）=====
// 【本发三个变量，全部是**前瞻距离常量**，不改变任何语义（CK 逐位相同）】
// `MF_PFB3` 8 → **12**：pass B 主循环里给「第 qi+D 个 query」预取的 4 条前瞻
// （`s_pm` 行 / `s_gt[d]` / `s_bd` 行 / `prefetchw out[j]`）的距离。
// `MF_PFB1` 64 → **16**、`MF_PFB2` 64 → **16**：pass B 每维计数排序两趟里的
// `s_bd` 前瞻距离（循环长度 ≈ N/7 ≈ 14,286 ⇒ 64 已超出整段）。
//
// 【定价通道（本席新交付，全队可用）= 判题机**多臂单进程**探针】
// `tools/probe.py C++17 <file>` 一次进程内跑 N 个「运行期选臂」的完整 solve，
// 每臂用 `__rdtscp` 包夹、值全部写 stderr ⇒ **一次探针 = N 个判题机原值 CYC**，
// 不占发次、单次进程 ~0.5–0.7 s（限 1000 ms）。
// 本席实测**同臂重复的噪声底 = 0.06%**（同批 ARM=0 两次：294,555,620 / 294,550,206），
// 跨批漂移 ~0.15% ⇒ **≥0.2% 的效应单批可判**，比"一发一件"省 ~8× 发次。
// （★ 口径：判题机侧 `CLOCK_MONOTONIC` 不可用（台阶 8000 ms）⇒ 必须用 rdtscp；
// 本席实测 250 M 拍量级上 rdtscp 重复性 0.06%，可用；相位必须靠**运行期选臂**差分。）
//
// 【判题机原值读数（同一批、臂序交错、CK 全部 `13586311538739599108` 逐位相同）】
// pf3=8(基座) 249,610,406 · **pf3=12 249,002,738（−0.24%）** · pf3=10 249,098,170(−0.21%)
// pf3=14 249,007,818(−0.24%) · pf3=16 247,569,806(−0.21%,另一批)
// pf1=16 + pf3=12 **248,918,914（−0.28%，本发取值）** · pf1=128 + pf3=12 249,073,472(−0.22%)
// 基座重复臂 249,877,300 / 249,601,936 ⇒ 同批噪声底 ≤0.11%
// ⇒ **本发 = 三项前瞻距离重标定，合计 −0.25~0.28%（判题机口径）**。
// ★ pf3=2/4/24/32 全部更差（+4.4%/+0.75%/+0.13%/+0.22%）⇒ 12 是新的局部最优。
//
// 【本轮一并**判负**的轴（全部判题机原值，留档给下一棒，勿重走）】
// · pass A「下一个 query 的 6 条池行」软件前瞻（1013a 姊妹题的成名刀）：判题机 **+1.4~+2.5%**
// （1/2/3 线、距离 2/4/8 全扫）⇒ **本题 pass A 是带宽墙、不是延迟墙，前瞻只会加压力**。
// · 「把 s_pm 的 32B 行载入提前一轮」（[z2-qr]，运行期可切的模板臂）：**+0.10%** ⇒ 无效。
// · 「用两条 8B 压缩签名行替换 32B 随机载入 + VBv」（[z2-crow]）：逐位等价（CK 同），
// 但判题机 **+16.7%**（294.5 vs 252.4 Mcyc）⇒ 该重构在本循环里撑大了环体，**净负**。
// （★ 但分量定价证明上限是真的：只删 `s_pm` 载入 **−3.64%**、只常量替 `Avm` **−1.28%**、
// 只删 qsb 的 6 条 shuffle **−0.93%** ⇒ **该池 ≥缺口，缺的是"不撑大环体"的形态**。）
// · 「去掉 pass B 四条前瞻中的任一条」= +1.9%~+12.1% ⇒ 四条全是承重件。
// · 判题机相位表（消融，基座 249.9 M）：pass A **32.0%** · pass B **37.4%**
// （主扫描 14.6% / 体后段 26.7%）· 计数排序**值 −11.4%（删掉更慢）** · 前瞻块值 **19.5%**。
// ======================
// ===== 思路（[b13z4-tf]，本席 b13z4_ 的唯一改动；正文其余部分与 [1] 逐字相同）=====
// 【本发两处改动，都不改任何计算结果：24 形状逐位闸门（n=100000/31000/20000/12000/
// 8000/5000/3000/1024/900/200/64/33 × span 100000/99991/512/37/7/3，含 9 个小规模
// 暴力 oracle 对照）与基座 CK 全同、oracle_bad=0】
//
// 【改动 1 = `[z4-tf]` 掩码尾块提前（本发大头）】
// 把 pass B 每-(i,d) 的「掩码尾块」（处理 r = (hi-lo) mod 32 的余数车道）从**主扫描循环
// 之后**移到**之前**：`aj` 预先置为 `ajEnd = lo + 32*floor((hi-lo)/32)`，尾块文本一字未
// 改（它只读 `aj`），随后主循环从 `lo` 正常跑 `ajEnd`。语义零变化（`aacc` 是求和、
// `amb_n` 是只追加表，两者都与顺序无关）。
// ★ 依据（本席判题机多臂探针实测，同批同臂重复 0.06%、跨批 0.02%）：
// · 相位定价：`掩码尾块 off` = **−3.63%**（9.5 M 拍 / 578,599 次 = **16.4 拍/次**），
// 而 `pass A AND 链 off` = −30.6%；**16.4 拍 ≈ 该块完整依赖链长度**
// （载入→vpsubsb→vpmaxsb×3→vpcmpeqb→vpmovmskb→popcnt→add ≈ 17 拍）
// ⇒ 尾块**延迟完全暴露**：它排在 (i,d) 最后，后面没有可重叠的活。
// · 把整块搬到主循环之前后，其延迟落在随后的 5~6 次主循环迭代（≈50 拍）之下。
// ★ 实测（判题机原值拍，臂序交错，base 双锚 244.50/244.54 M）：基座 244.50 M
// → **尾块提前 242.30 M（−0.88%）** → 再叠改动 2 **241.36 M（−1.30%）**。
//
// 【改动 2 = `[z4-tb]` 尾块 `base = aj` + 32 项 kmask 表】
// 尾块窗口取 `[aj, aj+32)`（原为 `[hi-32, hi)`），掩码 = `kmaskT[r]`（32 项表，1 次载入）
// 取代 `if (hi>=32u){base=hi-32;lk=32-r;} else {base=0;lk=lo;}` + 变量移位。
// ★ 安全性：`hi <= N <= MAXN` ⇒ `aj+31 <= MAXN+30 < MF_SPL = MAXN+64`，平面行内 64 B
// 余量内，无越界；窗口外车道的 **载入字节只落在被 kmask 屏蔽的位置**，不进 ltm/amb。
// 实测该项单独（旧结构下）−0.148%，叠在改动 1 上再拿 **−0.42%**。
//
// 【本轮一并判负的轴（判题机原值，留档勿重走）】
// · `[z4-fold]` 把 6 条平面载入折进 `vpsubsb`（把比较改成 `min_sl subs(qsb_sl, plane_sl)`，
// 使载入成为**被减数**、可用 `vpsubsb ymm,ymm,m256` 内存形，6 条载入变 1 uop/条）：
// 逐位等价、汇编确认已折叠，但判题机 **+2.03%**（244.50→249.47 M）⇒ **本环不是发射墙**：
// 折叠把 6 条独立载入绑进 max/min 树、**毁掉 MLP**（与 a14ah_「折叠只减 fused-domain
// uop、load 端口压力一字未变」互补：本环代价在**载入并行度**上）。
// ★ 旁证：把 `LDU256` 的**内联汇编**换成 `_mm256_loadu_si256` 后 gcc 因
// `-mavx256-split-unaligned-load` 把每条 32 B 载入拆成 `vmovdqu xmm`+`vinserti128`
// （1 条变 2 条）⇒ **原作的内联汇编宏正是在防这个**，别动它。
// · `[z4-unroll2]` 主扫描 2× 展开（每次迭代两个独立 32 车道窗口，12 条载入在飞）：
// **+4.49%**（244.50→255.44 M）⇒ 印证「**环体字节数不增**」：本环**正好驻留前端**，
// 撑大环体即付费（与上面「不是发射墙」并不矛盾：本环同时踩前端驻留与载入并行两条线）。
// · `[z4-pa-kz]` pass A 32 B 块全零短路（实测 10.12% 的块全零、26.74% 的 u64 全零；
// 逐位等价）：**+0.52%**；加「最稀疏行前置」（按桶号升序排 6 行）**+1.50%** ⇒
// 跳掉的载入省不过分支成本，pass A 是**带宽墙**（30.6%，940 MB，与上一棒一致）。
//
// 【仪器】判题机 `tools/probe.py` + 自建「运行期选臂」小驱动（臂 = z4 位掩码，
// 同一进程内多次完整 solve、`__rdtscp` 包夹回传原值）；同批同臂重复 ≤0.06%，
// 本席跨批锚差 ≤0.02% ⇒ ≥0.1% 的效应可判。
// ======================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/111503:
* copied the accepted signature/min-fold engine and preserved its original credits.
* saffah_cc_v41_agg1, https://duck.ac/submission/114343:
* adapted the eight-point blocked byte-transpose idea for signature-plane output.
* Both public sources displayed no separate software license notices.
* Idea: Replace scalar signature-plane scattering with an eight-record SSE byte
* transpose and 64-bit stores, retaining the existing compact cached signatures
* and excluding the ordering dimension exactly as in the original emitter.
* Purpose: Official experiment reducing plane-build store instructions and indexing.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，原提交 #110262
// <https://duck.ac/submission/110262>（判题机 Time = 71.676060 ms = 本题实时 T）
// —— 本文件正文 = **该提交的逐字节副本**；其自带引用/思路注释区原样保留在下方。
// 抄它的理由：**发前活取**（`§2.19.1348`）发现 T 已从派单给的 #110257（71.877608）推进到
// 本件（71.676060）⇒ 以活取件为准；比本席刚落地的 #110263（71.894579）快 0.218519 ms。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号）现役最好件 #110263
// <https://duck.ac/submission/110263>（71.894579 ms）= 上一代对手件 #110257 的副本 —— 闸门对照基座。
// [3] duck.ac 用户 saffah_cc_v41_agg1 #109794 <https://duck.ac/submission/109794> —— 本系上游，原样保留。
// 合规：duck.ac 提交正文按站点规则公开可见、可直接取用；原提交正文无独立许可证声明。
// 本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名作者账号与原提交地址。
// 本席（b13z_）在正文上做了**两处逐位等价的改动**（`[b13z-unroll-v2]` + `[b13z-sp]`，见下"思路"）：
// ★ [b13z-sp] = **pass B 的软件流水**：把第 qi+1 个 query 的 setup（`s_qord` → `s_bd[i*SBK+d]`
// → `FPOS` → `lo`，以及 `s_gt[d][i]`、`s_pm[i*SMP]` 的地址）提到**第 qi 个 query 的扫描之前**。
// 依据：pass B 的"每-(i,d) 固定部"实测 **77 cyc**（其中前端 uop 只有 ~20 cyc 当量）⇒ 约
// **57 cyc/(i,d) 是暴露的载入延迟**（`mask 128` 去 4 条前瞻 ⇒ +68 Mcyc 反证同一批载入是延迟源）；
// 而同一 (i,d) 的扫描只有 41 cyc 且是**前端受限** ⇒ 两者正好互补，可重叠。
// 只携带 5 个标量（V/i/lo/hi/rec），**不携带 qsb**（16×32 B ymm）⇒ ymm 压力不变；
// `V==0` ⇔ 原 `lo >= hi` 的 `continue`（逐位等价 ✓）。
// 除该处外，正文未作任何改动（无注释剥离、无空白压缩、无语义改动）。
// 【本席的刀】`solve_mf` 的签名平面转置里两个**运行期行程**的循环
// （`for(e<d) base[e*MF_SPL+j]=row[e]` 与 `for(e>d) base[(e-1)*MF_SPL+j]=row[e]`，
// 实测 22.0 M Ir / 4.2 M 次存储 = 5.2 条/存储，全部是循环簿记）展开成 7 个 case。
// 依据：`MF_SPL = (size_t)(MAXN + 64)` 是**编译期常量** ⇒ 六个目标是 `bp = base + j`
// 的常量位移（`mov disp(%rbp)`），只有源码下标 `e + (e>=d)` 是运行期的。
// switch 放在 **j 环之外**（放环内会被 -O2 编成比较链——反汇编 `jmp *` 计数 = 0 ⇒ 前端每轮多 ~7 分支）。
// 闸门：判题机同款 gcc9，CK = 13586311538739599108 = 基座（逐位相同）✓
// Ir（cachegrind n=20000 同形）：60,475,873 → ~58.8 M（−2.7%）
// ======================
// ===== 思路 =====
// 【抄件席 rc_（非本席新工作）】与 #110257 的差值（剥注释 diff **仅 1 hunk**）= **`#define MF_Q 256 → 254`**
// （桶数单变量；判题机上值约 −0.2 ms）⇒ 本发即该单变量的活取落地。
// 本地闸门：本题现成装置 14+ 形状差异测试（只取**题面域内** span<=n，含 tie-heavy 与 tiny n），
// **CK 非空且逐位相同**；正控由「域外形状会报 FAIL」提供（见 notes 记录）。
// ================
// References:
// - duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/109794 :
// Directly reuse the accepted compact-signature-plane 1013b engine and its
// inherited citations; the public source showed no separate license notice.
// - Our own experimental source, work/1013b/hist8_109794.cpp : reuse its
// 8-bit coordinate histogram with exact overflow fallback, officially
// accepted as duck.ac submission https://duck.ac/submission/110252 .
// ===== REFERENCES =====
// [1] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **duck.ac #108259**
// <https://duck.ac/submission/108259>（73.277790 ms，本题我方最好件）
// —— 用途：**本文件正文的基座**（正文逐字节取自该件，本发只在上面做下面这一处改动）。
// 该件本身的来源链见它自己的引用块（#107977 → … → #108259，全部为本账号/本团队自产）。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **duck.ac #102187**
// <https://duck.ac/submission/102187>（99.689594 ms）——`[ltsub]` 刀的出处（历史件）。
// [3] `/home/yjp/duck.ac/BRIEF.md` `§2.19.889`（消融表 × 指令账判型）· `§2.19.963`
// （判型必须用剂量臂，不能用占比相除）· `§2.19.1008/1009`（读写遍历序不同的数组先量写侧）。
// **未复制任何外部来源的代码**；本发是自产刀。
// [4] 本 work 区量具（团队内部工件）：`problems/1013b/work/b13y_instr.py`（逐相位 rdtsc 插桩，
// 判题机等价 flag）· `work/b13r_slot.py`（本发正文生成器，锚点断言）·
// `work/b13r_difftest.py`（差分闸门）· `work/b13r_slot_ab.txt`（配对 A/B 原始读数）。
// ======================
// ===== 思路 =====
// 【本发只有一个变量：pass-B 签名平面的**槽位布局**，[CSLOT]】
// 基座把签名平面存成 `s_sig[d][e*MF_SPL + j]`，`e = 0..k-1`，即**每个折叠维 d 都带一列 e==d 的平面**。
// 那一列的字节**从不被使用**：候选带 `[FPOS(d,bd_d), gt_d(i))` 按构造保证 `x_d[c] < x_d[i]`，
// 所以扫描时把第 d 条车道中性化（基座用哨兵 0x7F 把它并进 max 链）。代价：**每条 (query,dim)
// 都要多读 1/7 的平面字节**，多跑 1 条 `vpsubsb` + 1 条 `vpmaxsb`，并且 8 条车道的地址各需一个寄存器
// （反汇编里表现为每轮 4 条 `movq …(%rsp)` 的地址重载）。
//
// 【改法（逐位等价，算术一字未动）】
// * 平面改存 `s_sig[d][slot*MF_SPL + j]`，slot = 0..k-2，slot ↔ 真实车道 e = slot + (slot >= d)；
// * 两处建面循环把车道 e 写到槽 `(e < d ? e : e-1)`（循环按 e<d / e>d 拆成两段，各自步长常量）；
// * 主扫描块与带掩码的尾块只读 slot 0..k-2，`qsb[slot]` 仍取 sv0 的**第 e 个字节**
// （`t13_mskb[e & 3]` 与 ylo/yhi/zlo/zhi 的选择都按 e）⇒ 测试语义不变；
// * 车道-d 中性化哨兵与 k<=8 尾块里的 `d != 0 / d != 1` 守卫**随之删除**（它们唯一的作用就是中性化车道 d）；
// * 平面存储 4.9 MB → 4.2 MB（预算/分配口径同步）。
// 收益按"字节/访存"计：pass B 的平面读量少 1/7，热环少 3 条 uop/轮。**不涉及任何算法/桶/阈值。**
//
// 【正确性闸门（三条全过）】
// ① 生产形 n=1e5 的 FNV CK **与基座逐位相同** `17109020819928187212`；
// ② 差分闸门 `work/b13r_difftest.py`：**14 个形状**（n=1e5/2e4/5e3/1024/900/64/37/33/31 ×
// 坐标跨度 100000/99991/512/100/37/7/5/3 ⇒ 含**重并列**与极小 n）**两侧 CK 逐位全同，DIFFS 0**；
// ③ 两版都用判题端同款 flag（`-O2 -static -U_FORTIFY_SOURCE -std=c++17 -march=x86-64`）编过，
// 0 error（warning 与基座同源）。
//
// 【定价（本地配对 A/B）】判题机等价 flag、同机、`work/b13y_instr.py` 逐相位计数，
// 16 轮 **base/slot 交替**配对：**PB 13/16 轮为负**，均值 **−3.78 Mcyc（PB −4.9%）**，
// 整轮 TOT 均值 −2.62 Mcyc（−1.1%）。★ 如实标注：本机对 <3% 的效应只有方向意义
// （`§2.19.1242`：本机反号 5 例），**判题机读数是唯一终审** —— 本发的价格由这一次提交给出。
//
// 【注】本件按 `§2.19.096` **剥掉正文全部注释**（基座 103.4 KB 已贴近判题端长度上限），
// 再附本引用/思路块 ⇒ 45.9 KB + 2.5 KB。计时行与 checksum 行分开（`§2.19.092`）。
// ======================
// Approach:
// Keep the 8-bit coordinate histogram, compile seven separate query loops
// for fixed grouped dimensions, and fold the fixed seven-dimensional
// signature shift to nine bits. All three changes preserve exact results.
// Purpose:
// Test whether Q=254 bucket geometry combined with the histogram,
// signature shift and dimension specialization reaches the live goal.
#define MF_QCOEF 74
#define MF_Q 254
#define MF_NOSDMPB 1
#pragma GCC optimize("O3,no-strict-aliasing,unroll-loops,rename-registers,align-loops=32")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#include <immintrin.h>
#include <cstring>
static inline __m256i LDU(const void *p) {
__m256i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline __m256i LDA(const void *p) {
__m256i v;
__asm__("vmovdqa %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline void STU(void *p, __m256i v) {
__asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline void STA(void *p, __m256i v) {
__asm__ volatile("vmovdqa %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline __m128i LDQ(const void *p) {
__m128i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m128i *)p));
return v;
}
#define LDU_M(p) ({ __m256i _v; __asm__("vmovdqu %1, %0" : "=x"(_v) : "m"(*(const __m256i *)(p))); _v; })
#define LDA_M(p) ({ __m256i _v; __asm__("vmovdqa %1, %0" : "=x"(_v) : "m"(*(const __m256i *)(p))); _v; })
#define STU_M(p, v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i *)(p)) : "x"(_vv)); } while (0)
#define STA_M(p, v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqa %1, %0" : "+m"(*(__m256i *)(p)) : "x"(_vv)); } while (0)
extern "C" void *malloc(unsigned long);
typedef unsigned u32;
typedef unsigned short u16;
typedef unsigned char u8;
typedef unsigned long long u64;
#define AMBD_CAP 4096
static u32 s_abi[AMBD_CAP], s_abj[AMBD_CAP];
static __m256i s_abqT[AMBD_CAP], s_abqTh[AMBD_CAP];
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#ifndef MAXN
#define MAXN 100016
#define MF_SPL ((size_t)((MAXN + 64 + 4095) & ~4095))
#endif
#define MAXK 8
#define STRIDE 8
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 90
#endif
#ifndef MF_KSPAD
#define MF_KSPAD 1
#endif
#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
#ifndef MF_FB2
#define MF_FB2 0
#endif
#ifndef MF_FB3
#define MF_FB3 1
#endif
#ifndef MF_PLEN
#define MF_PLEN 1
#endif
#ifndef MF_BUDGETMB
#define MF_BUDGETMB 1900
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
#ifndef MF_PFB
#define MF_PFB 16
#define MF_PFB1 16
#define MF_PFB2 16
#define MF_PFB3 12
#define MF_PFB4 16
#define MF_PFB5 16
#define MF_PFB5 16
#endif
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
#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 512
#endif
static u32 s_ord[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
#define BDS 7
static u8 s_bd[(size_t)MAXN * BDS];
static const u32 g_SBK = BDS;
static u32 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u8 b13cnt8[MAXN];
static u16 s_bkt[MAXN];
static u32 s_nid[MAXN];
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN];
static u64 s_facc[(size_t)(MAXK - 1) * (((MAXN) + 63) / 64 + 16)];
static u32 s_bnde[(size_t)MAXK * (MAXQ + 2)];
static u64 *s_fbm = 0;
static u64 s_fbuf[8192] __attribute__((aligned(64)));
static __m256i khi8v, khi8hv;
static u32 g_Qc = 32;
static u32 g_fbm_W = 0;
static int g_use_filt = 0;
static int g_lexsort = 0;
static int g_nofilt = 0;
static u32 s_fpos[(MAXK + 1) * MAXQ];
static u32 s_bval[(MAXK + 1) * MAXQ];
static u64 *s_bs = 0;
static u64 s_acc[MAXN / 64 + 256] __attribute__((aligned(64)));
static u32 *s_sdmp[MAXK];
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
static u8 *s_sig = 0;
static u8 b13_sigrow[(size_t)MF_SPL * 7] __attribute__((aligned(64)));
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
static const u32 g_kmaskT[32] = {0u,1u,3u,7u,15u,31u,63u,127u,255u,511u,1023u,2047u,
4095u,8191u,16383u,32767u,65535u,131071u,262143u,524287u,1048575u,2097151u,
4194303u,8388607u,16777215u,33554431u,67108863u,134217727u,268435455u,
536870911u,1073741823u,2147483647u};
alignas(32) static const unsigned char t13_mskb[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}
};
#define T13SIGV(x) _mm256_srlv_epi32((x), shv)
#define T13THR(x) ((x) >= 1u)
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 const u32 g_KS = 8;
static unsigned long long g_bmag = 0;
static const u32 g_SMP = 8;
__attribute__((aligned(64))) static __m256i g_blm[MAXK][MAXK];
static int g_blm_ready = 0;
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
void *p = (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
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, acc = z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 4 <= Wr; w += 4) {
__m256i a0 = LDU256((bs[0] + w));
__m256i a1 = LDU256((bs[KK > 1 ? 1 : 0] + w));
for (int d = 2; d < KK; d += 2) a0 = _mm256_and_si256(a0, LDU256((bs[d] + w)));
for (int d = 3; d < KK; d += 2) a1 = _mm256_and_si256(a1, LDU256((bs[d] + w)));
__m256i a = (KK > 1) ? _mm256_and_si256(a0, a1) : a0;
__m256i lo = _mm256_and_si256(a, lm);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc = _mm256_add_epi8(acc, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
if (++cnt == 31) { tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z)); acc = z; cnt = 0; }
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z));
total = hsum256(tot);
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 1; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v);
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 1; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & ((1ull << rem) - 1));
}
return total;
}
template<int KK, int DS = -1>
__attribute__((target("avx2,popcnt")))
__attribute__((flatten)) static void query_group(u32 N, u32 *out, u32 ds_runtime, u32 q0, u32 q1, u32 Q, u32 W) {
const u32 ds = DS >= 0 ? (u32)DS : ds_runtime;
const int k = KK;
const int K1 = KK - 1;
#if MF_ASCR
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB1 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB1] * g_SBK + d1), _MM_HINT_T0);
#endif
{ u32 _bb = s_bd[(size_t)s_qord[qi] * g_SBK + d1];
s_bkt[qi] = _bb; 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++) {
u32 id = s_qord[qi];
s_tmp[cntv[s_bkt[qi] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#else
const u32 d1 = 0, e1 = 0;
#endif
if (g_lexsort) {
u32 *cntv = s_cntv;
for (int dd = k - 1; dd >= 0; dd--) {
if (dd == ds) continue;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) cntv[s_bd[(size_t)s_qord[qi] * g_SBK + dd] + 1]++;
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
u32 id = s_qord[qi];
if (qi + 16 < q1) __builtin_prefetch((const char *)(s_bd + (size_t)s_qord[qi + 16] * g_SBK + dd), 0, 3);
s_tmp[cntv[s_bd[(size_t)id * g_SBK + dd] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = s_qord[qi + PF_DIST];
const u8 *bdn = s_bd + (size_t)j * g_SBK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
#if MF_ASCR
{
u32 qi = q0;
while (qi < q1) {
u32 b1 = s_bd[(size_t)s_qord[qi] * g_SBK + d1];
u32 qe = qi + 1;
while (qe < q1 && s_bd[(size_t)s_qord[qe] * g_SBK + d1] == b1) qe++;
const u64 *bp1 = g_plen_on ? (s_bs + s_pofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
u32 mx = 1;
for (u32 t = qi; t < qe; t++) { u32 rr = s_gt[ds][s_qord[t]]; if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 1;
if (0) {
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = s_qord[qq];
const u8 *bd = s_bd + (size_t)i * g_SBK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * g_SBK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
if (e == (int)e1) bp[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
const u8 *bd = s_bd + (size_t)i * g_SBK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * g_SBK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = s_qord[qi + PF_DIST];
const u8 *bdn = s_bd + (size_t)j * g_SBK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W);
#ifdef DBG_PLEN
if (g_plen_on) {
u32 need = (R >> 6) + ((R & 63) ? 1 : 0);
if (need > s_plen[(size_t)e * MAXQ + bd[d]])
printf("DBG ds=%d e=%d d=%d b=%u R=%u need=%u len=%u\n", ds, e, d, bd[d], R, need, s_plen[(size_t)e * MAXQ + bd[d]]);
}
#endif
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
#endif
if (g_skip_fringe) return;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
#if !MF_FR_NOSORT
if (!(MF_ASCR && !g_lexsort && d == (int)d1)) {
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB2] * g_SBK + d), _MM_HINT_T0);
#endif
u32 bb = s_bd[(size_t)s_qord[qi] * g_SBK + d];
s_bkt[qi] = bb;
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++) {
u32 id = s_qord[qi];
s_tmp[cntv[s_bkt[qi] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#endif
u32 exmask = 0;
for (int e = 0; e < d; e++) if (e != ds) exmask |= 1u << e;
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);
int amb_n = 0;
const u32 *const odd = s_ord[d];
const u32 *const sdmpd = s_sdmp[d];
const u8 *const sbase = s_sig + (size_t)d * (size_t)(k - 1) * MF_SPL;
/* [b13z-sp] 软件流水：把第 (qi+1) 个 query 的 **setup**（含 4 条依赖载入
`s_qord -> s_bd[i*SBK+d] -> FPOS` 与 `s_gt[d][i]`、`s_pm[i*SMP]` 的地址）提到
**第 qi 个 query 的扫描之前**执行。扫描本身与下一轮的 setup 完全独立 ⇒
OoO 可在扫描执行期间把 setup 的载入跑完（扫描是前端受限、setup 是延迟受限，
两者正好互补）。只携带 5 个标量（V/i/lo/hi/rec），**不携带 qsb**（16×32 B 的
ymm 数组），故 ymm 压力不变；qr/qT/qTh 仍在 `if (cV)` 内由 rec 现算。
语义逐位不变：`V==0` ⇔ 原 `lo >= hi` 的 `continue`（跳过本轮全部剩余工作）✓ */
#define B13Z_SETUP(P, V, II, LL, HH, RR) \
do { \
V = 0u; \
if ((P) < q1) { \
u32 _i = s_qord[(P)]; \
const u8 *_b = s_bd + (size_t)_i * g_SBK; \
u32 _l = FPOS(d, _b[d]); \
u32 _h = s_gt[d][_i]; \
if (_l < _h) { V = 1u; II = _i; LL = _l; HH = _h; \
RR = s_pm + (size_t)_i * g_SMP; } \
} \
} while (0)
u32 cV, cI, cLo, cHi; const u32 *cRec;
B13Z_SETUP(q0, cV, cI, cLo, cHi, cRec);
for (u32 qi = q0; qi < q1; qi++) {
u32 nV, nI, nLo, nHi; const u32 *nRec;
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_pm + (size_t)(j) * g_SMP), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_T0);
_mm_prefetch((const char *)(s_bd + (size_t)j * g_SBK + d), _MM_HINT_T0);
__asm__ __volatile__("prefetchw %0" :: "m"(out[j]));
}
}
#endif
/* [b13z-sp] 下一轮的 setup 先发（其载入与下面的扫描重叠） */
B13Z_SETUP(qi + 1, nV, nI, nLo, nHi, nRec);
if (cV) {
const u32 i = cI, lo = cLo, hi = cHi;
const u32 *rec = cRec;
const u8 *bd = s_bd + (size_t)i * g_SBK; /* [b13z-sp] k>8 路径与死码仍引用 bd */
__m256i qr = _mm256_or_si256(LDA((const __m256i *)rec), khi8v);
__m256i qT, qTh;
if (k <= 8) {
__m256i VBv = _mm256_setr_epi32(
(int)((0 < d && 0 != (int)ds) ? BVAL(0, bd[0]) : 0x7FFFFFFFu),
(int)((1 < d && 1 != (int)ds) ? BVAL(1, bd[1]) : 0x7FFFFFFFu),
(int)((2 < d && 2 != (int)ds) ? BVAL(2, bd[2]) : 0x7FFFFFFFu),
(int)((3 < d && 3 != (int)ds) ? BVAL(3, bd[3]) : 0x7FFFFFFFu),
(int)((4 < d && 4 != (int)ds) ? BVAL(4, bd[4]) : 0x7FFFFFFFu),
(int)((5 < d && 5 != (int)ds) ? BVAL(5, bd[5]) : 0x7FFFFFFFu),
(int)((6 < d && 6 != (int)ds) ? BVAL(6, bd[6]) : 0x7FFFFFFFu),
(int)((7 < d && 7 != (int)ds) ? BVAL(7, bd[7]) : 0x7FFFFFFFu));
__m256i Avm = VBv;
qT = _mm256_min_epi32(qr, Avm);
qTh = all8;
} else {
u32 Av[STRIDE] __attribute__((aligned(32)));
for (int e = 0; e < k; e++) Av[e] = (e == ds) ? 0x7FFFFFFFu : BVAL(e, bd[e]);
for (int e = k; e < STRIDE; e++) Av[e] = 0x7FFFFFFFu;
__m256i Avm = _mm256_or_si256(LDU256(Av), notexlo);
qT = _mm256_min_epi32(qr, Avm);
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU256((Av + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
}
u32 add = 0;
const u32 *c0 = sdmpd ? (sdmpd + (size_t)lo * g_KS) : 0;
if (0) {
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) * g_KS) : (s_pm + (size_t)(s_ord[d][j]) * g_SMP);
__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 MF_NOSDMPB
if (1) {
#else
if (c0 && g_sig_on && T13THR(hi - lo)) {
#endif
__m256i qsb[16];
{
const __m256i shv = _mm256_set1_epi32((int)g_sigsh); (void)shv;
__m256i sv0 = _mm256_xor_si256((k == 7) ? _mm256_srli_epi32(qT, 9) : T13SIGV(qT), _mm256_set1_epi8((char)0x80));
const __m256i ylo = _mm256_broadcastsi128_si256(_mm256_castsi256_si128(sv0));
const __m256i yhi = _mm256_broadcastsi128_si256(_mm256_extracti128_si256(sv0, 1));
__m256i zlo = ylo, zhi = yhi;
if (k > 8) {
__m256i sv1 = _mm256_xor_si256(T13SIGV(qTh), _mm256_set1_epi8((char)0x80));
zlo = _mm256_broadcastsi128_si256(_mm256_castsi256_si128(sv1));
zhi = _mm256_broadcastsi128_si256(_mm256_extracti128_si256(sv1, 1));
}
for (int sl = 0; sl < k - 1; sl++) {
const int e = (sl < (int)d) ? sl : sl + 1;
__m256i src = (e < 8) ? ((e < 4) ? ylo : yhi) : ((e < 12) ? zlo : zhi);
qsb[sl] = _mm256_shuffle_epi8(src, _mm256_load_si256((const __m256i *)t13_mskb[e & 3]));
}
}
u32 aj = lo, aacc = 0, aexact = 0;
u32 ajEnd_ = aj + (((hi - aj) >> 5) << 5);
{ u32 ajSav_ = aj; aj = ajEnd_;
{
u32 r = hi - aj;
if (r) {
u32 base = aj, lk = 0u;
u32 kmask = g_kmaskT[r];
__m256i ac0 = _mm256_set1_epi8((char)0x80), ac1 = _mm256_set1_epi8((char)0x80);
ac0 = _mm256_subs_epi8(LDU256(sbase + base), qsb[0]);
ac1 = _mm256_subs_epi8(LDU256(sbase + MF_SPL + base), qsb[1]);
for (int sl = 2; sl < k - 1; sl += 2) {
{
__m256i cv = LDU256((sbase + (size_t)sl * MF_SPL + base));
ac0 = _mm256_max_epi8(ac0, _mm256_subs_epi8(cv, qsb[sl]));
}
if (sl + 1 < k - 1) {
__m256i cv = LDU256((sbase + (size_t)(sl + 1) * MF_SPL + base));
ac1 = _mm256_max_epi8(ac1, _mm256_subs_epi8(cv, qsb[sl + 1]));
}
}
__m256i ac = _mm256_max_epi8(ac0, ac1);
__m256i eqv = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = ((u32)_mm256_movemask_epi8(ac)) & kmask;
u32 amb = ((u32)_mm256_movemask_epi8(eqv)) & kmask;
aacc += (u32)__builtin_popcount(ltm);
if (amb) {
int a2 = amb;
while (a2) {
int t = __builtin_ctz(a2); a2 &= a2 - 1;
if (amb_n < AMBD_CAP) {
s_abi[amb_n] = i; s_abj[amb_n] = base + (u32)t;
_mm256_store_si256(&s_abqT[amb_n], qT);
if (k > 8) _mm256_store_si256(&s_abqTh[amb_n], qTh);
amb_n++;
} else {
const u32 *cp = c0 ? (c0 + (size_t)(base + (u32)t - lo) * g_KS)
: (s_pm + (size_t)(odd[base + (u32)t]) * g_SMP);
__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;
}
}
}
}
}
aj = ajSav_; }
for (; aj + 32 <= hi; aj += 32) {
__m256i ac0 = _mm256_set1_epi8((char)0x80), ac1 = _mm256_set1_epi8((char)0x80);
ac0 = _mm256_subs_epi8(LDU256(sbase + aj), qsb[0]);
ac1 = _mm256_subs_epi8(LDU256(sbase + MF_SPL + aj), qsb[1]);
for (int sl = 2; sl < k - 1; sl += 2) {
{
__m256i cv = LDU256((sbase + (size_t)sl * MF_SPL + aj));
ac0 = _mm256_max_epi8(ac0, _mm256_subs_epi8(cv, qsb[sl]));
}
if (sl + 1 < k - 1) {
__m256i cv = LDU256((sbase + (size_t)(sl + 1) * MF_SPL + aj));
ac1 = _mm256_max_epi8(ac1, _mm256_subs_epi8(cv, qsb[sl + 1]));
}
}
__m256i ac = _mm256_max_epi8(ac0, ac1);
__m256i eqv = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = (u32)_mm256_movemask_epi8(ac);
u32 amb = (u32)_mm256_movemask_epi8(eqv);
aacc += (u32)__builtin_popcount(ltm);
if (amb) {
int a2 = amb;
while (a2) {
int t = __builtin_ctz(a2); a2 &= a2 - 1;
if (amb_n < AMBD_CAP) {
s_abi[amb_n] = i; s_abj[amb_n] = aj + (u32)t;
_mm256_store_si256(&s_abqT[amb_n], qT);
if (k > 8) _mm256_store_si256(&s_abqTh[amb_n], qTh);
amb_n++;
} else {
const u32 *cp = c0 ? (c0 + (size_t)(aj + (u32)t - lo) * g_KS)
: (s_pm + (size_t)(odd[aj + (u32)t]) * g_SMP);
__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 + g_KS));
__m256i p2 = LDU256((c + 2 * g_KS));
__m256i p3 = LDU256((c + 3 * g_KS));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8);
c += 4 * g_KS;
}
for (; j < hi; j++, c += g_KS) {
__m256i pv = LDU256(c);
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
add = a0 + a1 + a2 + a3;
} else {
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 + g_KS));
__m256i p2 = LDU256((c + 2 * g_KS));
__m256i p3 = LDU256((c + 3 * g_KS));
__m256i h0 = LDU256((c + 8));
__m256i h1 = LDU256((c + g_KS + 8));
__m256i h2 = LDU256((c + 2 * g_KS + 8));
__m256i h3 = LDU256((c + 3 * g_KS + 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 * g_KS;
}
for (; j < hi; j++, c += g_KS) {
__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 {
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)(odd[j]) * g_SMP;
__m256i pv = LDA((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDA((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
cV = nV; cI = nI; cLo = nLo; cHi = nHi; cRec = nRec;
}
#undef B13Z_SETUP
if (amb_n) {
for (int z = 0; z + 12 < amb_n; z++)
__builtin_prefetch(s_pm + (size_t)odd[s_abj[z + 12]] * g_SMP, 0, 3);
for (int z = 0; z < amb_n; z++) {
const u32 *cp = s_pm + (size_t)odd[s_abj[z]] * g_SMP;
__m256i pv = LDU256(cp);
__m256i qv = _mm256_load_si256((const __m256i *)&s_abqT[z]);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qv, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
__m256i qvh = _mm256_load_si256((const __m256i *)&s_abqTh[z]);
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qvh, pvh), all8);
}
out[s_abi[z]] += ok;
}
amb_n = 0;
}
}
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int k) {
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) {
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++;
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
overhead += (u64)k * (u64)(k - 1) * (u64)MF_SPL;
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46};
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};
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;
if (g_Q > (u32)MAXQ) g_Q = (u32)MAXQ;
g_B = (N + g_Q - 1) / g_Q;
}
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
g_plen_on = 1;
#endif
#if MF_WTRIM
u32 Wpool = W;
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1;
#endif
if (g_use_filt && (u64)k * k * (g_Qc + 1) * W * 8 > ((u64)700 << 20)) g_use_filt = 0;
#if !MF_FILTAUTO
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
g_bmag = ((1ull << 40) + g_B - 1) / g_B;
const int K1 = k - 1;
if (!g_blm_ready) {
for (int a = 0; a < MAXK; a++)
for (int b = 0; b < MAXK; b++) {
u32 m[8];
for (int e = 0; e < 8; e++) m[e] = ((e < b && e != a) ? 0xFFFFFFFFu : 0u);
g_blm[a][b] = LDU256(m);
}
g_blm_ready = 1;
}
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
for (u32 i = 0; i <= N; i++) b13cnt8[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; b13cnt8[v < N ? v : 0]++; }
bool compact_counts=true;
{ u32 s = 0; for (u32 i = 0; i <= N; i++) { u32 c = b13cnt8[i]; s_tmp[i] = s; s += c; }
if(g_use_filt || MF_FB2)memcpy(s_cntv,s_tmp,(size_t)(N+1)*sizeof(u32));
if (s != N) {
compact_counts=false;
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]++; }
s = 0; for (u32 i = 0; i <= N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tmp[i] = s; s += c; }
}
}
#if MF_FB2
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i;
s_posp[(size_t)i * MAXK + d] = q;
}
#else
if (g_use_filt) {
u32 *pp = s_pos[d];
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i; pp[i] = q; }
} else if(compact_counts){
u32 *gt=s_gt[d];
for(u32 i=0;i<N;i++){
u32 v=xd[i];if(v>=N)v=0;
u32 w=s_tmp[v],st=w&0xfffffu;
gt[i]=st;
u32 n0=(s_tmp[v+1]&0xfffffu)-1u;if(n0>N)n0=N;
u32 bb=(u32)(((u64)n0*g_bmag)>>40);
s_bd[(size_t)i*g_SBK+d]=bb>Q-1?Q-1:bb;
od[st+(w>>20)]=i;
s_tmp[v]=w+(1u<<20);
}
} else {
u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
gt[i] = s_cntv[v];
u32 n0 = s_cntv[v + 1] - 1u; if (n0 > N) n0 = (u32)N;
u32 bb = (u32)(((u64)n0 * g_bmag) >> 40);
s_bd[(size_t)i * g_SBK + d] = bb > Q - 1 ? Q - 1 : bb;
od[s_tmp[v]++] = i;
}
}
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * g_SMP;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (u32 d = (u32)k; d < g_SMP; d++) p[d] = 0;
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
size_t need = (size_t)k * (size_t)(k - 1) * MF_SPL;
if ((size_t)N + 64 > MF_SPL) { s_sig = 0; g_sig_on = 0; }
static u8 *sigpool = 0; static size_t sigcap = 0;
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;
}
}
if (s_sig && k == 7 && (size_t)N <= MF_SPL) {
for (u32 i = 0; i < N; ++i) {
const u32 *row = s_pm + (size_t)i * g_SMP;
u8 *dst = b13_sigrow + (size_t)i * 7;
for (int e = 0; e < 7; ++e) dst[e] = (u8)((u8)(row[e] >> 9) ^ (u8)0x80);
}
}
#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
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), LDA((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), LDA((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
{
const u8 *bdv = s_bd + (size_t)i * g_SBK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
{
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; }
#if MF_NOSDMP
g_nosdmp = 1;
#endif
#if MF_NOSDMPB
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 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;
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) {
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)(od[j + MF_PFB4]) * g_SMP), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)(od[j]) * g_SMP;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * (size_t)(k - 1) * MF_SPL;
for (int e = 0; e < (int)d; e++) sb[(size_t)e * MF_SPL + j] = (u8)((u8)(row[e] >> g_sigsh) ^ (u8)0x80);
for (int e = (int)d + 1; e < k; e++) sb[(size_t)(e - 1) * MF_SPL + j] = (u8)((u8)(row[e] >> g_sigsh) ^ (u8)0x80);
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
if (s_sig && !g_sig_merged && k==7) {
for(int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*(k-1)*MF_SPL;
u32 j=0;
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(b13_sigrow+(size_t)od[j+64]*7),_MM_HINT_T0);
__m128i r0=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+0]*7));
__m128i r1=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+1]*7));
__m128i r2=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+2]*7));
__m128i r3=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+3]*7));
__m128i r4=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+4]*7));
__m128i r5=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+5]*7));
__m128i r6=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+6]*7));
__m128i r7=_mm_loadl_epi64((const __m128i *)(b13_sigrow+(size_t)od[j+7]*7));
__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_SPL+j),v0);
if(d!=1)_mm_storel_epi64((__m128i *)(base+(size_t)(1<d?1:1-1)*MF_SPL+j),_mm_srli_si128(v0,8));
if(d!=2)_mm_storel_epi64((__m128i *)(base+(size_t)(2<d?2:2-1)*MF_SPL+j),v1);
if(d!=3)_mm_storel_epi64((__m128i *)(base+(size_t)(3<d?3:3-1)*MF_SPL+j),_mm_srli_si128(v1,8));
if(d!=4)_mm_storel_epi64((__m128i *)(base+(size_t)(4<d?4:4-1)*MF_SPL+j),v2);
if(d!=5)_mm_storel_epi64((__m128i *)(base+(size_t)(5<d?5:5-1)*MF_SPL+j),_mm_srli_si128(v2,8));
if(d!=6)_mm_storel_epi64((__m128i *)(base+(size_t)(6<d?6:6-1)*MF_SPL+j),v3);
}
for(;j<N;j++) {
const u8 *row=b13_sigrow+(size_t)od[j]*7;
for(int e=0;e<k;e++)if(e!=d)base[(size_t)(e<d?e:e-1)*MF_SPL+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 * (size_t)(k - 1) * MF_SPL;
if (k == 7 && (size_t)N <= MF_SPL) {
/* [b13z-unroll-v2] (见首部"思路")：MF_SPL 是编译期常量 ⇒ 六个目标是 `bp = base + j`
的常量位移；把两个运行期行程的循环展开成 7 个 case，且 switch 提到 j 环之外
（放环内会被 -O2 编成比较链，`jmp *` 计数 = 0）。语义逐位不变：
槽 s 取 row[s]（s<d）或 row[s+1]（s>=d）⇔ "row 去掉第 d 个元素"。 */
#define B13Z_ROW_BODY(K0, K1, K2, K3, K4, K5) \
for (u32 j = 0; j < N; j++) { \
if (j + MF_PFB5 < N) \
_mm_prefetch((const char *)(b13_sigrow + (size_t)od[j + MF_PFB5] * 7), _MM_HINT_T0); \
const u8 *row = b13_sigrow + (size_t)od[j] * 7; \
u8 *bp = base + j; \
bp[(size_t)0 * MF_SPL] = row[K0]; bp[(size_t)1 * MF_SPL] = row[K1]; \
bp[(size_t)2 * MF_SPL] = row[K2]; bp[(size_t)3 * MF_SPL] = row[K3]; \
bp[(size_t)4 * MF_SPL] = row[K4]; bp[(size_t)5 * MF_SPL] = row[K5]; \
}
switch (d) {
case 0: B13Z_ROW_BODY(1, 2, 3, 4, 5, 6) break;
case 1: B13Z_ROW_BODY(0, 2, 3, 4, 5, 6) break;
case 2: B13Z_ROW_BODY(0, 1, 3, 4, 5, 6) break;
case 3: B13Z_ROW_BODY(0, 1, 2, 4, 5, 6) break;
case 4: B13Z_ROW_BODY(0, 1, 2, 3, 5, 6) break;
case 5: B13Z_ROW_BODY(0, 1, 2, 3, 4, 6) break;
default: B13Z_ROW_BODY(0, 1, 2, 3, 4, 5) break;
}
#undef B13Z_ROW_BODY
} else
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB5 < N) _mm_prefetch((const char *)(s_pm + (size_t)(od[j + MF_PFB5]) * g_SMP), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)(od[j]) * g_SMP;
for (int e = 0; e < (int)d; e++) base[(size_t)e * MF_SPL + j] = (u8)((u8)(row[e] >> g_sigsh) ^ (u8)0x80);
for (int e = (int)d + 1; e < k; e++) base[(size_t)(e - 1) * MF_SPL + j] = (u8)((u8)(row[e] >> g_sigsh) ^ (u8)0x80);
}
}
}
u64 poolwords = (u64)K1 * (Q + 1) * W;
#if MF_PLEN
if (g_plen_on) {
u64 best = 0;
for (int dsx = 0; dsx < k; dsx++) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < dsx) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)dsx * MAXK + d) * MAXQ;
for (u32 b = 0; b <= Q; b++) {
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > W) len = W;
tot += len;
}
}
if (tot > best) best = tot;
}
if (best && best < poolwords) poolwords = best;
}
#endif
size_t need = (size_t)poolwords * sizeof(u64);
static u64 *pool = 0; static size_t cap2 = 0;
if (need > cap2) { pool = (u64 *)pool_alloc(need + 4096); cap2 = need; }
s_bs = pool;
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
if (g_plen_on) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)ds * MAXK + d) * MAXQ;
u32 *po = s_pofs + (size_t)e * MAXQ;
u32 *pl = s_plen + (size_t)e * MAXQ;
for (u32 b = 0; b <= Q; b++) {
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; }
}
#endif
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
u64 *base = g_plen_on ? s_bs : (s_bs + (size_t)e * (Q + 1) * Wpool);
u64 *cur = s_acc;
for (u32 w = 0; w < Wpool; w++) cur[w] = 0;
if (g_plen_on) {
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++) {
u32 pend = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) { u32 q = s_nid[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = Wpool;
{
u32 mr_ = s_maxr[((size_t)ds * MAXK + (size_t)d) * MAXQ + b];
u32 nl = ((mr_ >> 6) + 1u + 7u) & ~7u;
if (nl < wlen) wlen = nl;
}
#if MF_NT
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), LDA((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7:
switch (ds) {
case 0: query_group<7,0>(N, out, ds, q0, q1, Q, Wpool); break;
case 1: query_group<7,1>(N, out, ds, q0, q1, Q, Wpool); break;
case 2: query_group<7,2>(N, out, ds, q0, q1, Q, Wpool); break;
case 3: query_group<7,3>(N, out, ds, q0, q1, Q, Wpool); break;
case 4: query_group<7,4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<7,5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<7,6>(N, out, ds, q0, q1, Q, Wpool); break;
}
break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 4); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 5); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 6); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 7); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 9); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 10); }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 67.261 ms | 23 MB + 352 KB | Accepted | Score: 100 | 显示更多 |