/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/126303: retained the accepted36-byte deferred-threshold ten-dimensional engine and all inherited citations.
* Idea: Replay exact ties in128-record batches and prefetch their output counters together with candidate point rows. A batch now fits its threshold records, candidate cache lines and output cache lines within the32KiB L1 budget instead of exceeding it at256 records. Every comparison and output update remains unchanged.
* Purpose: Reduce cache churn and write-miss stalls during deferred exact comparisons.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/126296: retained accepted ten-dimensional deferred tie resolution and all inherited citations.
* saffah_codex_6a_agg3 https://duck.ac/submission/126086 and https://duck.ac/submission/126075: adapted low16/high-bit threshold packing from nine/six-dimensional queues.
* Idea: Pack the first eight thresholds into20 bytes using16-bit low parts and eight4-bit high parts. Together with the two unchanged upper thresholds and query/order indices, each deferred record shrinks48 to36 bytes. BMI2 restores the eight nibbles in one operation. The INT_MAX neutral becomes2^20-1, which remains strictly greater than every valid coordinate (<N<=1000000). Packing is shared by all tied candidates in a signature block.
* Purpose: Reduce queue traffic while preserving exact predicates and the existing overflow fallback.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/126287: retained accepted ten-dimensional compact signature planes and all inherited references.
* saffah_codex_6a_agg3 https://duck.ac/submission/126082 and saffah_cc_v41_agg1 https://duck.ac/submission/123963: adapted deferred ambiguous-signature replay from the eight/seven-dimensional engines.
* Idea: Defer ten-dimensional exact signature ties into48-byte records containing the query, encoded order position and ten exact thresholds. After the dimension group, resolve point IDs in batches of256 and prefetch complete point rows before comparisons. Queue overflow and copied point arenas retain the original inline check.
* Purpose: Hide the dependent order-index and random point-record loads during exact tie resolution, preserving every threshold and output contribution.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/126079: retained the ten-dimensional compact counting and128-point prefix-mask implementation, with all inherited citations.
* saffah_codex_6a_agg3 https://duck.ac/submission/126143 and saffah_cc_v41_agg1 https://duck.ac/submission/122417: adapted storage of only live signature planes from the eight/nine-dimensional variants.
* Idea: Store nine signature planes per omitted coordinate and remap query broadcasts once, so the inner scan has nine fixed offsets and no omitted-dimension tests. Seed the union of saturated differences from the first live plane. Every original comparison remains represented and wide-count fallback is unchanged.
* Purpose: Measure reduced signature scan bookkeeping and storage in the ten-dimensional problem.
*/
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/123954: copied the accepted128-point prefix-mask ten-dimensional engine, with all leading source attributions preserved below.
- saffah_codex_6a_agg3, https://duck.ac/submission/126055 and https://duck.ac/submission/126047: byte histogram validation and combined20-bit-start/count cursor, adapted from our accepted1012a and1016b front ends. No separate software license notice is shown by these public sources.
Idea: Use byte frequencies when their complete sum proves no overflow, then fuse stable scatter, strict-rank output and bucket output with a single32-bit per-value cursor. Preserve the existing two-stage random-destination prefetch by predicting the cursor's exact start plus arrivals. Large tie groups retain the complete original wide-count path. The128-point bitmap algorithm and its four-stream early exit are unchanged.
Purpose: Remove a full million-point random-table pass per dimension and shrink the histogram working set while preserving stable ties and all recorded ranks.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号）提交 #122201 <https://duck.ac/submission/122201>
// （5764.068327 ms = 本账号 1016 现役最好件，全榜最快）= 本文件正文的逐字节来源。
// [2] duck.ac 用户 saffah_codex_6a_agg3 提交 #122014 <https://duck.ac/submission/122014>
// （5771.260688 ms = 实时排行榜 T）= #122201 正文的来源（NF=4 早退 + 桶秩选择排序），
// 其继承的 #121975/#121808/#121426/#120964 引用链原样保留于下。
// [3] 本账号提交 #121749 <https://duck.ac/submission/121749> —— [t16n_ PFXMASK] 掩码维
// 前缀流（把一条位集流降为每字 1 字节的前缀长度行）的来源；本刀把它再压缩一倍。
// 合规：duck.ac 提交正文按站点规则公开可见（题面「你提交的代码将会被公开」）；
// 上述原提交正文均无独立许可证声明。本文件按 RULES §3 逐项署名引用。
// ======================
// ===== 思路 =====
// 【t16t_ 席 · 第 2 发】入场 mine 5764.068327 (#122201) · T 5771.260688 (#122014) ·
// 严支 0.99T+1µs = 5713.549081 ⇒ 缺 50.519246 ms（0.876%）。
// 本发 = **[t16t_ MB128] 池的块编号由 64 点/块改为 128 点/块**，于是掩码维的前缀长度行
// 从「每 64 位 1 字节」变成「每 128 位 1 字节」——**掩码流字节与 s_tm 表双双减半**，
// 而池的位集数组（8 条流）一字未动。
// · 池位号：blk = rank_ds >> 7，块内位序仍按 rank_dmx ⇒ dim g_dm 的条件仍是「该块的前 c 位」，
// c ∈ [0,128] 恰好仍放得进 1 字节（原 c ∈ [0,64]）。
// · 掩码应用：一个字里两块 → 低字掩码 = 前 min(c,64) 位、高字掩码 = 前 max(c-64,0) 位。
// AVX2 上 `~(ones << k)` 对 k≥64 自动给全 1、对 k=0 给 0 ⇒ 只用一个 u8 饱和减法
// (`subs_epu8` 减 (0,64)) 就能同时得到两半，循环只多 ~2 个 uop（该环本为访存受限）。
// · 边缘掩码 s_edg 改为每 (块, 切口) 两个 u64（块内 128 位、切口 0..128），构造扫掠不变；
// ASCR 锚点暂存的拷贝长度 +1 词（切口块有第二个字）——这正是本发第一次闸门 13/62 处
// mismatch 的根因，修好后全绿。
// 闸门：暴力 oracle k=4..10 × mode 0/1/3/4（n=3000, Q=1642）**28/28 OK** + n=10000 mode 0/3 OK；
// 官形状 n=200000 mode 10/13/18 checksum **38652385**、n=400000 checksum **157051574**
// 与基座 #122201 **逐位相同** ✓；`tools/gcc93check.py`（无 LOCAL_TEST 判题形态）COMPILE OK ✓。
// 判负记录（同席第 1 发，勿重走）：给 AND 环/建趟加前瞻（QPF/SDPF/FRPF）= #122413 = +67.5 ms。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3 提交 #122014 <https://duck.ac/submission/122014>
// （5771.260688 ms，当前 T）= 本文件正文的逐字节来源；除本头部注释块外一字未改
// （`diff -u -w` 复核：正文 0 差异）。其继承的 #121975/#121808 引用链原样保留于下。
// [2] 本账号 saffah_cc_v41_agg1 提交 #121902 <https://duck.ac/submission/121902>
// （5967.213975 ms）= #122014 正文的来源（对手用我方前缀掩码引擎 + 早退）。
// 合规：duck.ac 提交正文按站点规则公开可见（题面“你提交的代码将会被公开”）；
// 上述原提交正文均无独立许可证声明。本文件按 RULES §3 逐项署名引用。
// ======================
// ===== 思路 =====
// 【t16p_ 席 · 第 2 发】入场 mine 5967.213975 (#121902) · T 5771.260688 (#122014)。
// 第 1 发 #122126 = #122014 正文逐字节（5771.710135，Accepted）。本发在其上叠**本席新刀
// [CLAMP]**：`and_popcount_pm` 的热环里有 3 条前瞻，其目标地址**越过本查询自己的区域末尾**：
// · 前 4 条流的 `bs[d] + w + 32`（lead=4 迭代）—— 最后 4 个迭代目标的字已 ≥ Wr；
// · `trow + w + 64`（lead=8 迭代）—— 最后 8 个迭代已越过前缀行；
// · 尾 4 条流的同款 lead-32（在生存分支内，按 84.2% 生存率发）。
// 这些行本查询永不读，等于白买 DRAM 行。本刀**只把地址夹回区域**（`(w+32<Wr)?32:0`、
// `(w+64<Wr)?64:0`），不移动任何 load、不改任何取值 ⇒ 逐位等价（本机 n=4e5 两臂 FNV
// 逐位同 = 5620543286167856259）。同进程两臂交替 min-of-3：2860.0 → 2842.5 ms（−0.61%，
// 三轮同向）；按判题机 pass A = 4+3.37+0.26 行/迭代的流量账，夹回的是每查询 ~37 行的白买。
// 闸门：官形状 k=10 n=1e6 Q=1642 mode 10/13/18 checksum 969232827 + out[] FNV 逐位同基座。
// （第 3 发候选 BNE=无分支早退已被本机否决：n=1e6 同进程两臂 min-of-3 11460 vs 10801 ms =
// +6.1%，机理 = 尾流 load 地址被 testz+cmov 拖到测试之后，失去投机发射；未发。）
// 历史（第 1 发）：§2.18.349 `diff -u -w` 我方 #121902 vs #122014 的**全部真改动只有 3 处**：
// (a) **AND 环早退**：先 AND「按桶秩升序的前 4 条流」+ 前缀行，512 位全空则 `continue` 跳过
// 其余 4 条流的载入。本席计数器实测（n=1e6 mode10 Q=1642）：176,891,354 块迭代中跳过
// 27,975,016（15.8%），尾流载入率 84.2% ⇒ 流量 8 → 7.37 行/迭代（NF=3/5 分别 7.61/7.40
// ⇒ 对手的 NF=4 已是最优；纯按生存率算的“NF=2 更省”被实测否决）。
// (b) 每查询对前 4 条流按桶秩（bd[rd]）做 4 趟选择排序（nbx>=8 时），使早退测试打在
// **最选择性的 4 个维**上。
// (c) s_edg 建趟里 `s_npos[sdsr[b0+j]]` 的地址提前 32 点前瞻（= 本账号 [PFXSW] 同族的
// 第二个“零前瞻随机访存”站点，本席复核：该环每点 1 次 4 MB 表随机读）。
// 闸门：官形状 k=10 n=1e6 Q=1642 mode 10/13/18 checksum 969232827 + out[] FNV
// 8996697766657014596 / 2559007080471616052 与基座 #121902 逐位相同 ✓；
// 判题同款编译 gcc9.3 -O2 -static -U_FORTIFY_SOURCE -std=c++17 rc=0 ✓。
// 本发为纯采纳件（本席的自造刀见下一发：BNE 无分支早退，正在定价）。
// ================
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/121975: copied prefix-mask
and prefetched edge-table engine, all inherited notices retained below.
Idea: In ten dimensions, intersect the four smallest-bucket bitsets and the
prefix mask first. If all512 resulting bits are zero, skip the remaining
four stream loads and population count for that block. Exact intersections
are commutative; every surviving block completes all original conditions.
Experiment purpose: Trade per-query pointer ordering and block zero checks
for less traffic from dimensions that cannot change an empty intersection.
*/
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/121902: copied the prefix-mask
engine and its prefetched rank-dimension permutation sweep. All inherited
author, source and license notices remain below; no independent license
notice was declared in the direct source.
Idea: Prefetch the second permutation lookup stream while constructing region
edge masks. Point IDs arrive in rank_ds order, so the address of each later
s_npos read is already known32 positions ahead. Cover across block boundaries
while retaining bounds checks and every edge-mask dependency.
Experiment purpose: Hide random permutation-read latency in edge-table setup.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3 提交 #121808 <https://duck.ac/submission/121808>
// （6023.593438 ms）= 本文件正文逐字节来源；本席相对它仅一处改动：PFXMASK 建表扫描环
// 加 2 条地址前瞻（`[t16o_ PFXSW]`，取值逐位不变）。
// [2] 本账号提交 #121749 <https://duck.ac/submission/121749>（6100.449484 ms）—— #121808
// 正文的来源，其继承的 #121426/#121470 引用链与历史注释原样保留于下方。
// 合规：duck.ac 正文按站点规则公开可见；原正文均无独立许可证声明。
// ======================
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/121749: compressed prefix-length
bitset representation and inherited ten-dimensional engine. All original
attribution and license notices are retained below.
Idea: Expand four prefix lengths directly with AVX2 byte-to-qword widening and
variable shifts of all-ones vectors. Complementing the shifted vector yields
exactly the low t bits, including t=0 and t=64. Replace eight scalar mask-table
loads and vector assembly in the main intersection loop with two expansions.
Experiment purpose: Reduce prefix-mask construction instructions and L1 loads
while preserving every mask, edge condition and count.
*/
// [t16n_] 1016 worker 席 —— pass A 的「掩码维前缀流」：把 AND 链 9 条位集流中的一条降为
// 每字 1 字节的前缀长度行（位集索引/算法形态改动；计数语义与其余相位一字未动）。
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，原提交 **#121426** <https://duck.ac/submission/121426>
// （判题机 6152.580955 ms = 本题 T）—— 本文件正文的逐字节来源（其自带 REFERENCES/思路 区原样保留在下方）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号）提交 **#121470** <https://duck.ac/submission/121470>
// = #121426 的逐字节副本（t16m_ 席落地，6152.591707 ms = 本席入场时的 mine）；其头部
// [t16m_]/[t16j_]/[t16e_]/[t16d_] 段与下方全部历史注释一并保留（仅按判题源码长度上限压掉叙述段）。
// [3] 本仓库 problems/1016/notes.md（t16m_ 席「全相位图 + pass A 带宽墙」与 MF_BUILD_NOPOOL 零页探针）
// 与 .cache/agents.md 全队相位账 —— 本刀方向（能动 66% 的只有减 pass A 字节）来自该账；未取其代码。
// 合规：duck.ac 提交正文按站点规则公开可见（题面「你提交的代码将会被公开」）；上述原提交正文均无独立许可证声明。
// ======================
// [t16m_] 1016 worker 席 —— 采纳当前 T 件（对手 #121426）正文：全值域签名映射 + 第二桶预排序 + MAXQ 上限。
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，原提交 #121426
// <https://duck.ac/submission/121426>（判题机 Time = 6152.580955 ms = 本题实时排行榜第一名 T）
// —— 本文件正文 = 该提交正文的**逐字节副本**（其自带引用/思路注释区原样保留在下方；本席 diff 0 处）。
// 该件自带引用链（其头部原样保留）：其 #121237（full-range cached-query engine = 坐标→255 桶定点映射）、
// 其 #121062 / #121171（second-bucket query locality 预排序）、以及我方 #120964（query/bound 签名 padding）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），原提交 #120964
// <https://duck.ac/submission/120964>（判题机 6196.734906 ms = 本账号 1016 现役最好件）
// —— #121426 正文中继承的 [t16e_ SIGPAD] 段即来自该件；本账号前序件另有 #120444/#120763。
// 合规：duck.ac 提交正文按站点规则公开可见、可直接取用（题面"你提交的代码将会被公开"）；
// 上述原提交正文均无独立许可证声明。本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名作者账号与原提交地址。
// ======================
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/121237: retained accepted full-range cached-query engine.
* saffah_cc_v41_agg1, https://duck.ac/submission/120964: inherited query/bound signature padding.
* All inherited citations retained; no separate license notices were displayed.
* Idea: Apply the signature-minus-one bias once when constructing cached bytes.
* Combine all query/bound bytes with one SIMD minimum and per-dimension masks,
* then broadcast with byte shuffles instead of scalar loads/minima/bias operations.
* Bound-byte loads stay within the same64-byte point record, including lane8.
* Purpose: Official experiment reducing repeated fringe query preparation.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/121171: retained full-range cached-signature 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: Preorder by a second non-fold bucket before the existing stable anchor
* sort, using its ping-pong scratch without a copy-back pass. Preserve exact queries.
* Purpose: Official experiment increasing bitset cache reuse within anchor runs.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/120964: retained cached-query ten-dimensional engine.
* saffah_codex_6a_agg3, https://duck.ac/submission/120939: reused full-range signature quantization.
* All inherited citations retained; no separate license notices were displayed.
* Idea: Map coordinates monotonically across all255 signature buckets using one
* fixed-point multiplier. Transform both candidate planes and cached query/bound
* bytes consistently, retaining exact checks whenever signatures are ambiguous.
* Also cap the inherited grid allocation to its MAXQ table capacity for small N.
* Purpose: Official experiment reducing false ambiguous fringe candidates.
*/
// [t16d_] 1016 worker 席 —— 删除池构建里每 ds 组重建的 s_nid 逆排列（10e6 次随机 RMW）。
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），原提交 #120444
// <https://duck.ac/submission/120444>（判题机 Time = 6240.233886 ms = 本账号 1016 现役最好件）
// —— 本文件正文 = 该提交正文的逐字节副本（其自带引用/思路注释区原样保留在下方）。
// 本席相对它只有"删除 + 改用等价既有数据"的 2 处改动（见下），不含任何新算术。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），原提交 #111659
// <https://duck.ac/submission/111659>（本题 6331.003381 ms）—— #120444 的入场基座，
// 含 MF_QCOEF 20155 / MF_NOFILTB 1 / MF_NOSDMP 1 / MF_NIDPF 12 / MF_PLEN 1 / MF_BUDGETMB 1940
// 等本题旋钮。本文件经 #120444 继承上述全部内容。
// [3] 本席使用本仓库内既有驱动 problems/1016a/work/rc_1016a_jmain.cpp（作者：本账号 a16y_ 席）作为
// "判题机路径"测量口径：它以 -Dmain=sub_main 抑制本题 LOCAL_TEST 自带 main，从而不改 g_Q/g_B，
// 使引擎走真实 auto-Q 预算式 = 判题机实际路径。本文件不含该驱动的任何代码。
// 合规：duck.ac 提交正文按站点规则公开可见、可直接取用（题面"你提交的代码将会被公开"）；
// 上述原提交正文均无独立许可证声明。本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名作者账号与原提交地址。
// ======================
// ================
// [t16e_] 1016 worker 席 —— 把 s_avpool 行载入从 pass B(fringe) 热路径中删除（纯删除随机访存）。
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），原提交 #120763
// <https://duck.ac/submission/120763>（判题机 Time = 6222.129260 ms = 本账号 1016 现役最好件）
// —— 本文件正文 = 该提交正文的逐字节副本（其自带引用/思路注释区原样保留在下方）。
// 本席相对它只有下述"删除 + 改用等价既有数据"的改动，不含任何新算术（无新公式/无新布局）。
// [2] 本席沿用 #120763 继承的全部旋钮与第三方引用（见下方 [t16d_]/[a16x_] 区），未改动它们。
// [a16x_] 1016 worker 席 —— 复刻对手 T (#120337) 的 3 处真改动 + 自带一刀 [FA]
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，原提交 #120337
// <https://duck.ac/submission/120337>（判题机 Time = 6241.633621 ms = 本题实时 T）
// —— 本文件正文 = 该提交的逐字节副本（其自带引用/思路注释区原样保留在下方）。
// 所抄内容（本席逐行 diff 取证：该件 = 本账号 #111659 逐字节 + 恰好 3 处真改动）：
// ① 顶部新增 `#define MF_NOSDMP 1` —— 弃用 s_sdmp 的 400 MB 转置行建表；
// ② 签名面扫描入口去掉 `c0 &&`（并加 `(hi-lo)>=2`），两处精确回退的 `cp` 改 c0 可空三元
// （c0==0 时回退读 s_pm 行）⇒ 无转置行时签名扫描仍然可用；
// ③ 新增独立「点主序 8 位签名字节数组 + 8 记录字节转置写面」块（其自述抄我方 #111877 的
// 预合并与 #114343 的 8 记录转置）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），原提交 #114343
// <https://duck.ac/submission/114343>（1014b，81.107107 ms）—— [FA]「无分支折叠维 argmin」
// 为该件首创；本件把它移植到本题（k=10）的折叠维选择环，是本件相对 #120337 的**唯一**正文改动。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号），原提交 #111659
// <https://duck.ac/submission/111659>（本题 6331.003381 ms = 本账号入场现役最好件）—— 入场基座，
// 含 [x16y2-DEADPF] 死预取删除与 MF_QCOEF 20155 / MF_NOFILTB 1 等本题旋钮；本文件正文经
// #120337 继承上述全部内容。
// 合规：duck.ac 提交正文按站点规则公开可见、可直接取用（题面"你提交的代码将会被公开"）；
// 上述原提交正文均无独立许可证声明。本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名作者账号与原提交地址。
// ======================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/111659:
* copied the accepted min-fold/signature engine and retained inherited references.
* saffah_cc_v41_agg1, https://duck.ac/submission/111877:
* adapted precomputation of compact point-major signature bytes.
* saffah_cc_v41_agg1, https://duck.ac/submission/114343:
* adapted eight-record byte transposition for signature-plane output.
* Public sources displayed no separate software license notices.
* Idea: Extend our successful small-case integration to this million-point task:
* omit full transposed-coordinate rows, gather exact signature survivors through
* point order, precompute signature bytes once, and emit 64-bit plane stores using
* an eight-record transpose. Preserve this engine's exact signature encoding.
* Purpose: Official experiment reducing redundant construction and first-touch work.
*/
#define MF_NOSDMP 1
#ifndef MF_PFXMASK
#define MF_PFXMASK 1
#define MF_PFXSW 24
#endif
#pragma GCC optimize("O3","unroll-loops","rename-registers","no-strict-aliasing","align-functions=64")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#ifndef MF_FRDPF
#define MF_FRDPF 32
#endif
#ifndef MF_CNTPF
#define MF_CNTPF 64
#endif
#ifndef MF_GQPF
#define MF_GQPF 32
#endif
#define MF_QCOEF 20155
#define MF_NOFILTB 1
#include <immintrin.h>
static inline __m256i LDU(const void *p) {
__m256i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline __m256i LDA(const void *p) {
__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)
#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 short u16;
typedef unsigned char u8;
typedef unsigned long long u64;
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#ifndef MAXN
#define MAXN 1000005
#define MF_SPL ((size_t)1000064)
#define MF_PFB1 64
#define MF_PFB2 64
#define MF_PFB3 8
#endif
#define MAXK 10
#define STRIDE 16
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 20155
#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 1940
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
#ifndef MF_PFB
#define MF_PFB 24
#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_PFAOUT
#define MF_PFAOUT 1
#endif
#ifndef MF_PFAOUTD
#define MF_PFAOUTD 24
#endif
#ifndef MF_SCAPF
#define MF_SCAPF 16
#endif
#define Z16S_PFW(p) __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(p)))
#ifndef MF_ORDPF
#define MF_ORDPF 24
#endif
#ifndef MF_NIDPF
#define MF_NIDPF 12
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 4096
#endif
static u32 s_ord[MAXK][MAXN];
static u32 s_frord[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u16 s_bd[(size_t)MAXN * MAXK];
static u16 s_bdT[(size_t)MAXK * MAXN];
static u32 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
alignas(64) static u8 compact_hist[MAXN+8];
alignas(64) static u32 compact_cursor[MAXN+8];
static u16 s_bkt[MAXN];
static u16 s_bka[MAXN];
static u32 s_nid[MAXN];
static u32 s_npos[MAXN] __attribute__((aligned(64)));
static u32 s_blkk[64], s_blkp[64];
static u8 *s_tm = 0; static size_t s_tmcap = 0;
static u64 *s_edg = 0; static size_t s_edgcap = 0;
static u64 s_pref[256] __attribute__((aligned(64)));
static int g_dm = -1;
static u32 g_tstride = 0;
static u32 s_tofs[MAXQ + 1], s_rlen[MAXQ + 1];
static u32 s_cntx[16384];
static int g_mask_on = 0;
static int g_mask_ok = 0;
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, kbadhi, kokhi;
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 u32 *s_avpool = 0;
static size_t s_avcap = 0;
static u8 *s_sig = 0;
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
static u32 g_sigmul;
static inline u32 sigfull(u32 x){return (x*g_sigmul)>>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;
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
return (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
}
static constexpr u32 DEFER_CAP=1u<<19;
alignas(64) static u32 defer_entries[(size_t)DEFER_CAP*9],defer_points[DEFER_CAP];
static u32 defer_count;
static void flush_deferred(u32*out){
u32 n=defer_count;defer_count=0;
for(u32 b=0;b<n;b+=128){
u32 end=b+128<n?b+128:n;
for(u32 j=b;j<end;j++){
const u32*e=defer_entries+(size_t)j*9;
u32 id=s_ord[e[1]>>20][e[1]&0xfffffu];defer_points[j]=id;
__builtin_prefetch(out+e[0],1,3);
_mm_prefetch((const char*)(s_pm+(size_t)id*STRIDE),_MM_HINT_T0);
}
for(u32 j=b;j<end;j++){
const u32*e=defer_entries+(size_t)j*9;
const u32*p=s_pm+(size_t)defer_points[j]*STRIDE;
__m256i lo=_mm256_cvtepu16_epi32(_mm_loadu_si128((const __m128i*)(e+2)));
__m256i hi=_mm256_cvtepu8_epi32(_mm_cvtsi64_si128((long long)_pdep_u64(e[6],0x0f0f0f0f0f0f0f0full)));
__m256i th=_mm256_or_si256(lo,_mm256_slli_epi32(hi,16));
__m256i a=_mm256_cmpgt_epi32(th,_mm256_load_si256((const __m256i*)p));
__m128i h=_mm_cmpgt_epi32(_mm_loadl_epi64((const __m128i*)(e+7)),_mm_loadl_epi64((const __m128i*)(p+8)));
out[e[0]]+=(_mm256_movemask_epi8(a)==-1)&&((_mm_movemask_epi8(h)&255)==255);
}
}
}
__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;
for (; w + 8 <= Wr; w += 8) {
if (KK > 1) _mm_prefetch((const char *)(bs[1] + w + 32), _MM_HINT_T0);
if (KK > 2) _mm_prefetch((const char *)(bs[2] + w + 32), _MM_HINT_T0);
if (KK > 3) _mm_prefetch((const char *)(bs[3] + w + 32), _MM_HINT_T0);
if (KK > 4) _mm_prefetch((const char *)(bs[4] + w + 32), _MM_HINT_T0);
if (KK > 5) _mm_prefetch((const char *)(bs[5] + w + 32), _MM_HINT_T0);
if (KK > 6) _mm_prefetch((const char *)(bs[6] + w + 32), _MM_HINT_T0);
if (KK > 7) _mm_prefetch((const char *)(bs[7] + w + 32), _MM_HINT_T0);
if (KK > 8) _mm_prefetch((const char *)(bs[8] + w + 32), _MM_HINT_T0);
__m256i a = LDU256((bs[0] + w));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < KK; d++) {
a = _mm256_and_si256(a, LDU256((bs[d] + w)));
b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
}
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)));
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)));
if (++cnt == 31) {
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 NS>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_pm(const u64 *const *bs, const u8 *trow, const u64 *edgp, u32 Wf, u32 cut) {
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();
const __m256i ones=_mm256_set1_epi64x(-1);
const __m128i sub64v=_mm_set1_epi16(0x4000);
__m256i tot = z, acc0 = z, acc1 = z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 8 <= Wf; w += 8) {
constexpr int NF=NS<4?NS:4;
const u32 ld32 = (w + 32u < Wf) ? 32u : 0u;
for (int d = 0; d < NF; d++) _mm_prefetch((const char *)(bs[d] + w + ld32), _MM_HINT_T0);
_mm_prefetch((const char *)(trow + (w >> 1) + ((w + 64u < Wf) ? 32u : 0u)), _MM_HINT_T0);
__m256i a = LDU256((bs[0] + w));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < NF; d++) {
a = _mm256_and_si256(a, LDU256((bs[d] + w)));
b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
}
u32 tb = *(const u32 *)(trow + (w >> 1));
__m128i tv=_mm_cvtsi32_si128((int)tb);
__m128i tdup=_mm_unpacklo_epi8(tv,tv);
__m128i tsub=_mm_subs_epu8(tdup,sub64v);
a = _mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(tsub)),ones));
b = _mm256_and_si256(b,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(_mm_srli_si128(tsub,4))),ones));
if constexpr(NS>4){
__m256i live=_mm256_or_si256(a,b);
if(_mm256_testz_si256(live,live))continue;
for(int d=4;d<NS;d++){
_mm_prefetch((const char*)(bs[d]+w+ld32),_MM_HINT_T0);
a=_mm256_and_si256(a,LDU256(bs[d]+w));
b=_mm256_and_si256(b,LDU256(bs[d]+w+4));
}
}
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)));
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)));
if (++cnt == 31) {
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 <= Wf; w += 4) {
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < NS; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
u32 tb = *(const u16 *)(trow + (w >> 1));
__m128i tv=_mm_cvtsi32_si128((int)tb);
__m128i tdup=_mm_unpacklo_epi8(tv,tv);
__m128i tsub=_mm_subs_epu8(tdup,sub64v);
a = _mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(tsub)),ones));
__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 < Wf; w++) {
u64 v = bs[0][w];
for (int d = 1; d < NS; d++) v &= bs[d][w];
u32 c = trow[w >> 1];
u32 k = (w & 1u) ? ((c > 64u) ? c - 64u : 0u) : c;
total += (u32)__builtin_popcountll(v & s_pref[k]);
}
{
u32 c = trow[Wf >> 1];
u64 v0 = bs[0][Wf]; for (int d = 1; d < NS; d++) v0 &= bs[d][Wf];
u64 v1 = bs[0][Wf+1]; for (int d = 1; d < NS; d++) v1 &= bs[d][Wf+1];
total += (u32)__builtin_popcountll(v0 & s_pref[c] & edgp[(size_t)cut * 2u]);
total += (u32)__builtin_popcountll(v1 & ((c > 64u) ? s_pref[c - 64u] : 0ull) & edgp[(size_t)cut * 2u + 1u]);
}
return total;
}
#define T16E_QTCALC() \
const u32 *Avp_ = avpool + (size_t)i * STRIDE; \
__m256i qr_ = _mm256_or_si256(LDA((const __m256i *)rec), khi8v); \
__m256i qT = _mm256_min_epi32(qr_, _mm256_or_si256(LDU256(Avp_), notexlo)); \
__m256i qTh; \
if (k > 8) { \
__m256i qrh_ = _mm256_or_si256(LDU256((rec + 8)), khi8hv); \
__m256i Avmh_ = _mm256_or_si256(LDU256((Avp_ + 8)), notexloh); \
qTh = _mm256_min_epi32(qrh_, Avmh_); \
} else { qTh = all8; }
template<int KK>
__attribute__((target("avx2,popcnt")))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const u32 KS_C = (KK <= 8) ? 8u : (u32)KK;
const int k = KK;
const int K1 = KK - 1;
u32 *qo = s_qord, *qt = s_tmp;
#if MF_ASCR
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
const u32 dmx = (u32)g_dm, emx = (dmx < ds) ? dmx : dmx - 1u;
{const u32 d2=(ds<=1u)?2u:1u;u32 *cv=s_cntv;
for(u32 b=0;b<=Q;b++)cv[b]=0;
for(u32 qi=q0;qi<q1;qi++)cv[s_bdT[(size_t)d2*N+qo[qi]]+1]++;
u32 ac=q0;for(u32 b=0;b<=Q;b++){u32 c=cv[b];cv[b]=ac;ac+=c;}
for(u32 qi=q0;qi<q1;qi++){u32 id=qo[qi];qt[cv[s_bdT[(size_t)d2*N+id]+1]++]=id;}
u32*sw=qo;qo=qt;qt=sw;}
{
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_bdT + (size_t)d1 * N + qo[qi + MF_PFB1]), _MM_HINT_T0);
#endif
{ u32 bb = s_bdT[(size_t)d1 * N + qo[qi]];
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 = qo[qi]; u32 pos = cntv[s_bkt[qi] + 1]++;
qt[pos] = id; s_bka[pos] = s_bkt[qi]; }
{ u32 *swp = qo; qo = qt; qt = swp; }
}
#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_bdT[(size_t)dd * N + qo[qi]] + 1]++;
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) { u32 id = qo[qi]; qt[cntv[s_bdT[(size_t)dd * N + id] + 1]++] = id; }
{ u32 *swp = qo; qo = qt; qt = swp; }
}
for (u32 qi = q0; qi < q1; qi++) s_bka[qi] = s_bdT[(size_t)d1 * N + qo[qi]];
}
#define PF_DIST 4
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = qo[qi + PF_DIST];
const u16 *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + (g_plen_on ? s_pofs[(size_t)e * MAXQ + bdn[d]]
: ((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_bka[qi];
u32 qe = qi + 1;
while (qe < q1 && s_bka[qe] == b1) qe++;
const u64 *bp1 = g_plen_on ? (s_bs + s_pofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
u32 mx = 1;
for (u32 t = qi; t < qe; t++) { u32 rr = s_gt[ds][qo[t]]; if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 2;
if (nwd <= AS_WORDS) {
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = qo[qq];
const u16 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = qo[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
u16 brank[MAXK];
int nbx = 0;
for (int e = 0; e < K1; e++) {
if (g_mask_ok && e == (int)emx) continue;
int rd=(e<ds)?e:e+1; brank[nbx]=bd[rd];
if (e == (int)e1) bp[nbx++] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[nbx++] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
if(g_mask_ok && nbx>=8){
for(int z=0;z<4;z++){
int best=z;
for(int j=z+1;j<nbx;j++)if(brank[j]<brank[best])best=j;
const u64*p=bp[z];bp[z]=bp[best];bp[best]=p;
u16 rr=brank[z];brank[z]=brank[best];brank[best]=rr;
}
}
cnt = g_mask_ok
? and_popcount_pm<KK - 2>(bp, s_tm + s_tofs[bd[dmx]], s_edg + (((size_t)(R >> 7) * 129u) << 1), (R >> 7) << 1, R & 127u)
: and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = qo[qi];
const u16 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = qo[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#if MF_PFAOUT
{ u32 oo = qi + MF_PFAOUTD;
if (oo < q1) _mm_prefetch((const char *)(out + qo[oo]), _MM_HINT_T0); }
#endif
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = qo[qi + PF_DIST];
const u16 *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + (g_plen_on ? s_pofs[(size_t)e * MAXQ + bdn[d]]
: ((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;
u32 *const avpool = s_avpool;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
const u32 *fo = s_frord[d];
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);
const u32 *const odd = s_ord[d];
alignas(16) u8 qmask_bytes[16],dmask_bytes[16];
for(int e=0;e<16;e++){qmask_bytes[e]=(e<d&&e!=(int)ds)?0:255;dmask_bytes[e]=(e==d)?255:0;}
const __m128i qmask=_mm_load_si128((const __m128i*)qmask_bytes),dmask=_mm_load_si128((const __m128i*)dmask_bytes);
for (u32 qi = q0; qi < q1; qi++) {
u32 i = fo[qi] & 0xFFFFFu;
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = fo[qj] & 0xFFFFFu;
_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 *)(out + j), _MM_HINT_T0);
}
}
#endif
const u32 lo = FPOS(d, fo[qi] >> 20);
const u32 hi = s_gt[d][i];
if (lo >= hi) continue;
const u32 *rec = s_pm + (size_t)i * STRIDE;
const u8 *const qpad = ((const u8 *)rec) + 40;
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (size_t)lo * KS_C) : 0;
if (g_use_filt) { T16E_QTCALC();
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_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (g_sig_on && (hi - lo) >= 2) {
__m256i qsb[16];
{
__m128i qr8=_mm_loadu_si128((const __m128i*)qpad);
__m128i qb8=_mm_loadl_epi64((const __m128i*)(qpad+k));
if constexpr(KK>=10)qb8=_mm_insert_epi8(qb8,qpad[k+8],8);
__m128i qs=_mm_or_si128(_mm_min_epu8(qr8,_mm_or_si128(qb8,qmask)),dmask);
__m256i qwide=_mm256_broadcastsi128_si256(qs);
for(int e=0;e<k-1;e++)qsb[e]=_mm256_shuffle_epi8(qwide,_mm256_set1_epi8((char)(e<d?e:e+1)));
}
const u8 *const sgbase = s_sig + (size_t)d * (k-1) * MF_SPL;
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
__m256i ac = _mm256_subs_epu8(LDU256(sgbase + aj), qsb[0]);
for (int e = 1; e < k-1; e++) {
__m256i cv = LDU256(sgbase + (size_t)e * MF_SPL + aj);
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[e]));
}
__m256i lt = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = (u32)_mm256_movemask_epi8(lt);
u32 amb = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)));
aacc += (u32)__builtin_popcount(ltm);
if (amb) { T16E_QTCALC();
__m256i lv=_mm256_and_si256(qT,_mm256_set1_epi32(65535));
__m128i packed_low=_mm_packus_epi32(_mm256_castsi256_si128(lv),_mm256_extracti128_si256(lv,1));
__m256i hv=_mm256_and_si256(_mm256_srli_epi32(qT,16),_mm256_set1_epi32(15));
__m128i hw=_mm_packus_epi32(_mm256_castsi256_si128(hv),_mm256_extracti128_si256(hv,1));
__m128i hb=_mm_packus_epi16(hw,hw);
u32 packed_high=(u32)_pext_u64((u64)_mm_cvtsi128_si64(hb),0x0f0f0f0f0f0f0f0full);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
if(KK==10 && !c0 && defer_count<DEFER_CAP){
u32*e=defer_entries+(size_t)defer_count++*9;
e[0]=i;e[1]=(aj+(u32)t)|((u32)d<<20);
_mm_storeu_si128((__m128i*)(e+2),packed_low);e[6]=packed_high;
_mm_storel_epi64((__m128i*)(e+7),_mm256_castsi256_si128(qTh));
continue;
}
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_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
aacc += okl;
}
}
}
{
u32 r = hi - aj;
if (r) {
u32 base, lk;
if (hi >= 32u) { base = hi - 32u; lk = 32u - r; }
else { base = 0u; lk = lo; }
u32 kmask = ((r >= 32u) ? 0xFFFFFFFFu : ((1u << r) - 1u)) << lk;
__m256i ac = _mm256_subs_epu8(LDU256(sgbase + base), qsb[0]);
for (int e = 1; e < k-1; e++) {
__m256i cv = LDU256(sgbase + (size_t)e * MF_SPL + base);
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[e]));
}
__m256i lt = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = ((u32)_mm256_movemask_epi8(lt)) & kmask;
u32 amb = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(ltm);
if (amb) { T16E_QTCALC();
__m256i lv=_mm256_and_si256(qT,_mm256_set1_epi32(65535));
__m128i packed_low=_mm_packus_epi32(_mm256_castsi256_si128(lv),_mm256_extracti128_si256(lv,1));
__m256i hv=_mm256_and_si256(_mm256_srli_epi32(qT,16),_mm256_set1_epi32(15));
__m128i hw=_mm_packus_epi32(_mm256_castsi256_si128(hv),_mm256_extracti128_si256(hv,1));
__m128i hb=_mm_packus_epi16(hw,hw);
u32 packed_high=(u32)_pext_u64((u64)_mm_cvtsi128_si64(hb),0x0f0f0f0f0f0f0f0full);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
if(KK==10 && !c0 && defer_count<DEFER_CAP){
u32*e=defer_entries+(size_t)defer_count++*9;
e[0]=i;e[1]=(base+(u32)t)|((u32)d<<20);
_mm_storeu_si128((__m128i*)(e+2),packed_low);e[6]=packed_high;
_mm_storel_epi64((__m128i*)(e+7),_mm256_castsi256_si128(qTh));
continue;
}
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_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
aacc += okl;
}
}
}
}
add = aacc;
} else if (c0) {
T16E_QTCALC();
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 {
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_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h0)), kokhi);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8) & (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h1)), kokhi);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8) & (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h2)), kokhi);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8) & (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h3)), kokhi);
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_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
add = a0 + a1 + a2 + a3;
}
} else {
T16E_QTCALC();
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__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_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
add += ok;
}
}
out[i] += add;
}
}
}
template<int K>
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out) {
const int k = 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 * k * ((u64)N + 64);
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
#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;
if(q>=MAXQ)q=MAXQ-1;
g_Q = q;
g_B = (N + q - 1) / q;
}
for (u32 t = 0; t < 256u; t++) s_pref[t] = (t >= 64u) ? ~0ull : ((1ull << t) - 1ull);
g_mask_on = MF_PFXMASK;
const u32 Q = g_Q, B = g_B;
const u64 bmag = ((1ull << 40) + B - 1) / 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], bad8[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;
for (int e = 0; e < 8; e++) bad8[e] = (8 + e < k) ? 0 : -1;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
kbadhi = LDU256(bad8);
kokhi = _mm256_xor_si256(kbadhi, _mm256_set1_epi32(-1));
}
const int K1 = k - 1;
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
u8*hc=compact_hist;u32*ct=compact_cursor;
for(u32 i=0;i<=N;i++)hc[i]=0;
u32 lim=N>MF_ORDPF?N-MF_ORDPF:0;
for(u32 i=0;i<lim;i++){
u32 v=xd[i+MF_ORDPF];if(v>=N)v=0;Z16S_PFW(hc+v);
v=xd[i];hc[v<N?v:0]++;
}
for(u32 i=lim;i<N;i++){u32 v=xd[i];hc[v<N?v:0]++;}
u32 sum=0;for(u32 v=0;v<=N;v++){ct[v]=sum;sum+=hc[v];}
if(sum==N && N<1048576u){
u32*gt=s_gt[d];u32*pp=s_pos[d];
for(u32 i=0;i<N;i++){
if(i+2*MF_SCAPF<N){u32 v=xd[i+2*MF_SCAPF];if(v>=N)v=0;_mm_prefetch((const char*)(ct+v),_MM_HINT_T0);}
if(i+MF_SCAPF<N){u32 v=xd[i+MF_SCAPF];if(v>=N)v=0;u32 w=ct[v];Z16S_PFW(od+(w&0xfffffu)+(w>>20));}
u32 v=xd[i];if(v>=N)v=0;
u32 w=ct[v],st=w&0xfffffu,pos=st+(w>>20);ct[v]=w+(1u<<20);
od[pos]=i;pp[i]=pos;gt[i]=st;
#if MF_FB2
s_posp[(size_t)i*MAXK+d]=pos;
#endif
u32 num=(ct[v+1]&0xfffffu)-1u;if(num>N)num=N;
u32 bb=(u32)(((u64)num*bmag)>>40);
s_bdT[(size_t)d*N+i]=(u16)(bb>Q-1?Q-1:bb);
}
}else{
for (u32 i = 0; i <= N; i++) s_cntv[i] = 0;
#if MF_ORDPF > 0
{ u32 i = 0;
for (; i + MF_ORDPF < N; i++) {
u32 vn = xd[i + MF_ORDPF];
Z16S_PFW(s_cntv + (vn < N ? vn : 0));
u32 v = xd[i]; s_cntv[v < N ? v : 0]++;
}
for (; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
}
#else
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
#endif
{ u32 s = 0; for (u32 i = 0; i <= N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tmp[i] = s; s += c; } }
#if MF_FB2
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i;
s_posp[(size_t)i * MAXK + d] = q;
s_pos[d][i] = q;
}
#else
{
u32 *pp = s_pos[d];
for (u32 i = 0; i < N; i++) {
if (i + 2*MF_SCAPF < N) { u32 v2 = xd[i + 2*MF_SCAPF]; if (v2 >= N) v2 = 0;
_mm_prefetch((const char *)(s_tmp + v2), _MM_HINT_T0); }
if (i + MF_SCAPF < N) { u32 vf = xd[i + MF_SCAPF]; if (vf >= N) vf = 0;
Z16S_PFW(&od[s_tmp[vf]]); }
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i; pp[i] = q;
}
}
#endif
{ u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
{if (i + MF_CNTPF < N) { u32 vp = xd[i + MF_CNTPF]; if (vp >= N) vp = 0;
_mm_prefetch((const char *)(s_cntv + vp), _MM_HINT_T0); }}
u32 v = xd[i]; if (v >= N) v = 0;
gt[i] = s_cntv[v];
u32 num = s_cntv[v + 1] - 1u; if (num > N) num = N;
u32 bb = (u32)(((u64)num * bmag) >> 40);
{ u16 bbx = (bb > Q - 1 ? Q - 1 : (u16)bb);
s_bdT[(size_t)d * N + i] = bbx; }
} }
}
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++) {
u16 *row = s_bd + (size_t)i * MAXK;
for (int d = 0; d < k; d++) row[d] = s_bdT[(size_t)d * N + i];
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
g_sigmul=N>1?(u32)((255ULL<<23)/(N-1)):0;
size_t need = (size_t)k * (k-1) * MF_SPL;
static u8 *sigpool = 0; static size_t sigcap = 0;
if (need > ((size_t)220 << 20) || (size_t)N + 32 > MF_SPL) { 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 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 needav = (size_t)N * STRIDE * 4 + 256;
if (needav > s_avcap) {
s_avpool = (u32 *)pool_alloc(needav);
s_avcap = s_avpool ? needav : 0;
}
}
{
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++) {
u64 bpk = ((u64)s_gt[0][i] << 4);
for (int d = 1; d < k; d++) {
u64 g = ((u64)s_gt[d][i] << 4) | (u64)d;
bpk = (g < bpk) ? g : bpk;
}
u32 best = (u32)(bpk >> 4);
int ds = (int)(bpk & 15u);
s_g[i] = (u32)ds;
const u16 *bdv = s_bd + (size_t)i * MAXK;
u32 *av = s_avpool + (size_t)i * STRIDE;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
for (int d = 0; d < k; d++) {
if (d == ds) { av[d] = 0x7FFFFFFFu; continue; }
u32 b = bdv[d];
av[d] = BVAL(d, b);
if (g_plen_on && best > mr[(size_t)d * MAXQ + b])
mr[(size_t)d * MAXQ + b] = best;
}
}
if (s_sig) {
const u32 sh = g_sigsh;
for (u32 i = 0; i < N; i++) {
u32 *row = s_pm + (size_t)i * STRIDE;
const u32 *av = s_avpool + (size_t)i * STRIDE;
u8 *pad = (u8 *)row + 40;
for (int e = 0; e < k; e++) { u32 v = sigfull(row[e]); pad[e] = (u8)(v ? v-1u : 0u); }
for (int e = 0; e + 1 < k; e++) { u32 v = sigfull(av[e]); pad[k + e] = (u8)(v ? v-1u : 0u); }
}
}
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++) {
{if (i + MF_GQPF < N) _mm_prefetch((const char *)(s_g + s_qord[i + MF_GQPF]), _MM_HINT_T0);}
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; }
for (int d = 0; d < k; d++) {
u32 pos[MAXK];
for (int ds = 0; ds < k; ds++) pos[ds] = gstart[ds];
for (u32 r = 0; r < N; r++) {
{if (r + MF_FRDPF < N) { u32 ip = s_ord[d][r + MF_FRDPF];
_mm_prefetch((const char *)(s_g + ip), _MM_HINT_T0);
_mm_prefetch((const char *)(s_bdT + (size_t)d * N + ip), _MM_HINT_T0); }}
u32 id = s_ord[d][r];
u32 dst = pos[s_g[id]]++;
s_frord[d][dst] = id | ((u32)s_bdT[(size_t)d * N + id] << 20);
}
}
#if MF_NOSDMP
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) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * (k-1) * MF_SPL;
for (int e = 0; e < k; e++) if(e!=d) sb[(size_t)(e-(e>d)) * MF_SPL + j] = (u8)((sigfull(row[e])) | ((sigfull(row[e])) == 0));
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
if (s_sig && !g_sig_merged && k==10) {
static u8 cached_signature[(size_t)MAXN*16] __attribute__((aligned(64)));
for(u32 i=0;i<N;i++)for(int e=0;e<10;e++){u32 v=sigfull(x[e][i]);cached_signature[(size_t)i*16+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_SPL;
u32 j=0;
for(;j+8<=N;j+=8) {
if(j+64<N){const u32*op=od+j+64;
_mm_prefetch((const char *)(cached_signature+(size_t)op[0]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[1]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[2]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[3]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[4]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[5]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[6]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[7]*16),_MM_HINT_T0);}
__m128i r0=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+7]*16));
__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);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
if(d!=0)_mm_storel_epi64((__m128i *)(base+(size_t)(0-(0>d))*MF_SPL+j),v0);
if(d!=1)_mm_storel_epi64((__m128i *)(base+(size_t)(1-(1>d))*MF_SPL+j),_mm_srli_si128(v0,8));
if(d!=2)_mm_storel_epi64((__m128i *)(base+(size_t)(2-(2>d))*MF_SPL+j),v1);
if(d!=3)_mm_storel_epi64((__m128i *)(base+(size_t)(3-(3>d))*MF_SPL+j),_mm_srli_si128(v1,8));
if(d!=4)_mm_storel_epi64((__m128i *)(base+(size_t)(4-(4>d))*MF_SPL+j),v2);
if(d!=5)_mm_storel_epi64((__m128i *)(base+(size_t)(5-(5>d))*MF_SPL+j),_mm_srli_si128(v2,8));
if(d!=6)_mm_storel_epi64((__m128i *)(base+(size_t)(6-(6>d))*MF_SPL+j),v3);
if(d!=7)_mm_storel_epi64((__m128i *)(base+(size_t)(7-(7>d))*MF_SPL+j),_mm_srli_si128(v3,8));
if(d!=8)_mm_storel_epi64((__m128i *)(base+(size_t)(8-(8>d))*MF_SPL+j),v4);
if(d!=9)_mm_storel_epi64((__m128i *)(base+(size_t)(9-(9>d))*MF_SPL+j),_mm_srli_si128(v4,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*16;
for(int e=0;e<k;e++)if(e!=(int)d)base[(size_t)(e-(e>d))*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 * (k-1) * MF_SPL;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) if(e!=d) base[(size_t)(e-(e>d)) * MF_SPL + j] = (u8)((sigfull(row[e])) | ((sigfull(row[e])) == 0));
}
}
}
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
g_mask_ok = 0;
if (g_mask_on) {
{ u32 dd0 = (ds == 0) ? 1u : 0u; g_dm = 0; while ((u32)g_dm == (u32)ds || (u32)g_dm == dd0) g_dm++; }
const u32 dmx = (u32)g_dm;
const u32 *const sdsp = s_pos[ds], *const sdmo = s_pos[dmx], *const sdom = s_ord[dmx], *const sdsr = s_ord[ds];
const u32 *const mrx = s_maxr + ((size_t)ds * MAXK + dmx) * MAXQ;
u32 mblk = 1; { u32 acc = 0;
for (u32 b = 0; b <= Q; b++) {
u32 rr = (mrx[b] >> 7) + 1u; if (rr > Wpool) rr = Wpool; if (rr > mblk) mblk = rr;
u32 rl = (rr + 63u) & ~63u; s_tofs[b] = acc; s_rlen[b] = rl; acc += rl; }
if (mblk > Wpool) mblk = Wpool;
{ size_t tnd = acc;
if (tnd > s_tmcap) { s_tm = (u8 *)pool_alloc(tnd + 4096); s_tmcap = s_tm ? tnd : 0; } } }
{ size_t end = ((size_t)mblk + 1) * 129 * 2;
if (end > s_edgcap) { s_edg = (u64 *)pool_alloc(end * 8 + 4096); s_edgcap = s_edg ? end : 0; } }
if (s_tm && s_edg) {
g_mask_ok = 1; g_tstride = mblk;
{ u32 *const cnt = s_cntx;
for (u32 m = 0; m < mblk; m++) cnt[m] = 0;
u32 j2 = 0;
for (u32 b = 0; b <= Q; b++) {
u32 L = (b == 0u) ? 0u : ((b < Q) ? FPOS(dmx, b) : (u32)N);
while (j2 < L) {
{ u32 pf = j2 + MF_PFXSW; if (pf >= N) pf = j2; u32 pp = sdom[pf];
_mm_prefetch((const char *)(sdsp + pp), _MM_HINT_T0); Z16S_PFW(&s_npos[pp]); }
u32 p = sdom[j2]; j2++;
u32 rk = sdsp[p]; u32 blk = rk >> 7;
if (blk < mblk) { u32 c = cnt[blk]; cnt[blk] = c + 1u; s_npos[p] = (blk << 7) | c; }
else s_npos[p] = rk;
}
{ u8 *row = s_tm + s_tofs[b]; u32 rl = s_rlen[b];
for (u32 m = 0; m < rl; m++) row[m] = (u8)cnt[m]; }
}
}
{ u32 nb = (N + 127u) >> 7; if (mblk > nb) mblk = nb;
for (u32 m = 0; m < mblk; m++) {
u32 b0 = m << 7, nj = N - b0; if (nj > 128u) nj = 128u;
u64 eml = 0, emh = 0; s_edg[(size_t)m * 258] = 0; s_edg[(size_t)m * 258 + 1] = 0;
for (u32 j = 0; j < nj; j++) {
u32 look=b0+j+32;
if(look<N)_mm_prefetch((const char*)(s_npos+sdsr[look]),_MM_HINT_T0);
u32 sg = s_npos[sdsr[b0 + j]] & 127u;
if (sg < 64u) eml |= 1ull << sg; else emh |= 1ull << (sg - 64u);
s_edg[(size_t)m * 258 + (j + 1u) * 2u] = eml; s_edg[(size_t)m * 258 + (j + 1u) * 2u + 1u] = emh;
}
for (u32 j = nj + 1u; j <= 128u; j++) { s_edg[(size_t)m * 258 + j * 2u] = eml; s_edg[(size_t)m * 258 + j * 2u + 1u] = emh; } } }
} else g_mask_on = 0;
}
const u32 *s_nid = g_mask_ok ? s_npos : s_pos[ds];
const u32 emsk = (g_mask_ok ? ((u32)g_dm < (u32)ds ? (u32)g_dm : (u32)g_dm - 1u) : 0xffffffffu);
for (int e = 0; e < K1; e++) {
if ((u32)e == emsk) continue;
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;
#if MF_NIDPF > 0
{ u32 pp = ptr;
for (; pp + MF_NIDPF < pend; pp++) {
_mm_prefetch((const char *)(s_nid + od[pp + MF_NIDPF]), _MM_HINT_T0);
u32 q = s_nid[od[pp]]; cur[q >> 6] |= 1ull << (q & 63);
}
for (; pp < pend; pp++) { u32 q = s_nid[od[pp]]; cur[q >> 6] |= 1ull << (q & 63); }
}
ptr = pend;
#else
while (ptr < pend) { u32 q = s_nid[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
#endif
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
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), LDA((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: defer_count=0;
query_group<10>(N, out, ds, q0, q1, Q, Wpool);
flush_deferred(out); break;
}
}
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf<4>((u32)n, (const unsigned **)x, out); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf<5>((u32)n, (const unsigned **)x, out); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf<6>((u32)n, (const unsigned **)x, out); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf<7>((u32)n, (const unsigned **)x, out); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf<8>((u32)n, (const unsigned **)x, out); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf<9>((u32)n, (const unsigned **)x, out); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf<10>((u32)n, (const unsigned **)x, out); }
#ifdef LOCAL_TEST
static void brute_nd(int n, const u32 **x, u32 *out, int k) {
for (int i = 0; i < n; i++) {
u32 c = 0;
for (int j = 0; j < n; j++) {
int ok = 1;
for (int d = 0; d < k; d++) if (!(x[d][j] < x[d][i])) { ok = 0; break; }
c += ok;
}
out[i] = c;
}
}
int main(int argc, char **argv) {
int k = (argc > 1) ? atoi(argv[1]) : 6;
int n = (argc > 2) ? atoi(argv[2]) : 200;
int seed = (argc > 3) ? atoi(argv[3]) : 1;
int mode = (argc > 4) ? atoi(argv[4]) : 0;
int qq = (argc > 5 && atoi(argv[5]) > 0) ? atoi(argv[5]) : 0;
g_Q = qq ? (u32)qq : (u32)(n / 250 + 1);
g_B = (u32)((n + g_Q - 1) / g_Q);
if (argc > 6) g_nosort = atoi(argv[6]);
if (argc > 7) g_nosdmp = atoi(argv[7]);
if (argc > 8) g_skip_and = atoi(argv[8]);
if (argc > 9) g_skip_fringe = atoi(argv[9]);
if (argc > 10) g_use_filt = atoi(argv[10]);
if (argc > 11) g_lexsort = atoi(argv[11]);
const int dm = (mode >= 10) ? mode - 10 : mode;
srand(seed);
static u32 *xs[10];
u32 *out = (u32 *)malloc(sizeof(u32) * n);
u32 *ref = (u32 *)malloc(sizeof(u32) * n);
for (int d = 0; d < k; d++) {
xs[d] = (u32 *)malloc(sizeof(u32) * n);
for (int i = 0; i < n; i++) {
if (dm == 0) xs[d][i] = (u32)(rand() % n);
else if (dm == 1) xs[d][i] = (u32)(rand() % 3);
else if (dm == 2) xs[d][i] = 0;
else if (dm == 3) xs[d][i] = (u32)(n - 1 - (rand() % n));
else if (dm == 4) xs[d][i] = (u32)(i);
else xs[d][i] = (u32)(rand() % n);
}
}
if (mode >= 10) {
struct timespec t0, t1;
clock_gettime(CLOCK_MONOTONIC, &t0);
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
clock_gettime(CLOCK_MONOTONIC, &t1);
double ms = (t1.tv_sec - t0.tv_sec) * 1e3 + (t1.tv_nsec - t0.tv_nsec) / 1e6;
u64 s = 0; for (int i = 0; i < n; i++) s += out[i];
printf("k=%d n=%d Q=%u TIME %.1f ms checksum %llu\n", k, n, g_Q, ms, s);
return 0;
}
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
brute_nd(n, (const u32 **)xs, ref, k);
int bad = 0;
for (int i = 0; i < n; i++) if (out[i] != ref[i]) { if (bad < 5) printf("MISMATCH i=%d got=%u ref=%u\n", i, out[i], ref[i]); bad++; }
printf("k=%d n=%d mode=%d Q=%u %s (%d mismatches)\n", k, n, mode, g_Q, bad ? "FAIL" : "OK", bad);
return bad != 0;
}
#endif
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 5.663 s | 985 MB + 484 KB | Accepted | Score: 100 | 显示更多 |