/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140865 supplies streamed byte-prefix snapshots and the current engine; all inherited citations/licenses remain below. No separate license declared.
Idea: during inverse-edge construction, the sequential anchor-order ID stream indexes pfx_pos randomly. Prefetch that position record16 points before its use, crossing64-point block boundaries while checking the global N bound. This keeps the original two-pass construction, all stores and exact inverse permutations; it does not change query-side edge hints.
Purpose: overlap dependent random rank-record loads in the construction phase on the newly faster snapshot implementation. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140688 supplies the current engine, byte prefix snapshots and existing bitmap streaming stores, with all inherited source/license notices below. No separate license was declared.
Idea: the prefix-count snapshot rows are64-byte aligned, every row length is a multiple of64 and the source counter array is also64-byte aligned. Emit their exact copies with two256-bit non-temporal stores per cacheline, then fence once after all rows before any consumer. This keeps live mutable counters resident while avoiding destination read-for-ownership during the one-pass snapshot construction. No geometry, query path, edge representation or byte value changes.
Purpose: test reducing prefix-snapshot write allocation and cache displacement; allocation and fence ordering preserve complete repeated-call semantics. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140582 supplies the current signed-signature engine and all inherited source/license notices below. No separate license was declared.
Idea: extract the equality byte mask once and use that scalar directly for both the zero branch and exact-candidate traversal. On the tail, intersect it with the valid-lane mask before the existing branch. This removes a vector-wide zero test and the later conditional extraction without changing any ambiguity bit, definite count, coordinate or replay order.
Purpose: measure the simpler vector-to-scalar control path in the current seven-dimensional scanner. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140557 supplies the complete seven-coordinate engine and inherited notices below. No new external code is used; no separate license was declared.
* Idea: all eight query-signature coordinates use the same g_sigsh. Replace the lane-variable VPSRLVD and its replicated count vector with the exact uniform-count VPSRLD using one scalar count in XMM. Preserve packing, bias, comparisons, raw coordinates, bounds, source signature construction and every query count.
* Purpose: official experiment shortening per-query signature preparation without changing the representation or output.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140499 supplies the accepted seven-coordinate engine and all inherited notices. Our https://duck.ac/submission/140142 applied the last-coordinate exclusion-bound invariant to a ten-coordinate scanner; this applies that invariant to the separate sequential record producer here. No separate license was declared.
* Idea: a fringe at coordinate d restricts only earlier coordinates e<d, excluding its anchor. Thus bound k-1 is never observed by any fringe. Fill that final bound slot with the same neutral INT_MAX as padded slots instead of reading its bucket index and BVAL entry. Preserve the64-byte record, all raw coordinates, every live bound, query masks and all count logic.
* Purpose: remove one unnecessary bound lookup per point from sequential setup while keeping the existing query engine and memory layout.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140471 supplies the accepted compact edge-prefix engine; all inherited author URLs and license notices follow. No separate license was declared; no new external code is used.
* Idea: Replace each 64-entry u64 cumulative edge table with its exact inverse permutation: one byte per bitmap position stores that point's original offset in its 64-point anchor block. Two signed byte comparisons against rem and two movemasks reconstruct exactly the old cumulative mask; offsets and rem are 0..63, and absent positions use127. Full bitmap words and rank semantics are unchanged.
* Purpose: reduce edge construction stores and its cache footprint eightfold, paying a few SIMD instructions only in the partial-word query. This is an official performance experiment with complete output comparison, including partial blocks and repeated calls.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140465 supplies the accepted byte-prefix construction engine and all inherited notices below. No separate license was declared; no external code is introduced.
* Idea: an edge row is consumed only at rem=R&63, hence only entries0..63 can be read. Store exactly64 entries, omitting the unused64-point endpoint. Limit generated edge blocks and Wpool to the largest rank actually queried by this anchor; obtain that rank from one existing complete maxr bucket row, with the original query scan retained when compact-span metadata is disabled. Keep initialization of ALL byte block counters so repeated calls preserve the full0..64 cardinality invariant even outside the queried edge range.
* Purpose: remove unconsumed edge construction/storage and simplify edge-row addressing. All reachable prefix bits, query ranks, strict comparisons and later fringe data remain identical; no point inside any queried prefix is omitted.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140450 supplies the accepted seven-dimensional prefix-mask engine and all inherited notices below. No separate license was declared; this change introduces no external code.
* Idea: each pfx_counts entry counts the points in one disjoint64-position block. All points are inserted once in sorted order, so every entry is always in0..64, including ties and partial final blocks. Store these exact counters as bytes and copy them directly into the existing byte prefix rows. Keep the original clear extent, mapped positions, row lengths, edge construction, query kernels and fringe code unchanged.
* Purpose: remove three quarters of this construction counter traffic and eliminate the repeated32-bit-to-byte narrowing pass. The bound comes from block cardinality, not an input-distribution assumption.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140257 supplies the accepted bounded-metadata seven-dimensional engine with all inherited notices below. Our https://duck.ac/submission/140337 retains an explicit D0 query-bound omission in the smaller sibling; this candidate extends that proof to the first live fringe dimension when the anchor is dimension0. No separate license was declared.
* Idea: the first fringe dimension other than the anchor has no prior exclusion dimensions. In D0, and in D1 when ds==0, the query threshold therefore equals the raw coordinate vector. Skip the unused bound-vector load, OR and minimum on these two statically selected cases. Later fringe dimensions retain the original exact bound computation.
* Purpose: remove redundant query setup and memory traffic without changing point records, bucket geometry, cached signatures, comparison thresholds, bitmap work or fringe ordering.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/139517 supplies the current million-point seven-dimensional engine and all inherited notices below. Our https://duck.ac/submission/140142 supplies exact capacity-bounded cache-line metadata rows in the ten-dimensional300k sibling. No separate license was declared; inherited credits remain intact.
* Idea: for every1<=N<=1000000, ceil(sqrt(N))<=1000 and the unchanged coefficient130/scale108 imply Q<=1404 before memory capping. Metadata indices0..Q inclusive fit1405 slots, rounded to1408 for64-byte-aligned rows. Replace the physical4096-slot stride only. Runtime Q/B, strict ranks, bitmaps, signatures and all query geometry remain unchanged; the unrelated nine-dimensional override is not executed by this seven-dimensional task.
* Purpose: remove provably unused capacity from initialization and metadata working sets, without reusing the negative joint-original-signature candidate or changing active dimensions.
*/
// References: saffah_v41f1017_agg4 https://duck.ac/submission/134933 supplies the current 7-coordinate engine and inherited notices.
// saffah_codex_6a_agg3 https://duck.ac/submission/139423 provides SIMD packed rank/dimension minima and direct stable fold grouping; https://duck.ac/submission/139506 provides exact sorted bucket intervals and transposed signature production.
// Idea/Purpose: retain each 1024-point maxr phase, but compute eight independent fold minima together and count group membership immediately, then scatter groups directly without a second histogram/copy/rescan.
// Build seven-byte cached signatures eight points at a time with exact shift/XOR and a byte transpose, replacing 56 scalar record stores per batch.
// Original prefix, query, arithmetic and allocation policies are unchanged; scalar tails and all dimensions remain supported.
// ===== 引用区(本发增量;题 1013)=====
// [1] 本账号 saffah_v41f1017_agg4,提交 #134917 <https://duck.ac/submission/134917>(题 1013,3.808 s,Accepted,
// 当前 best):本文件代码正文 = 该提交正文的逐字节副本 + 下述单变量改动(见思路区);#134917 自带的完整
// 引用链/署名/许可说明(含 #134869 的 qg 簿记 slot-1 死存储删除、#134907 的 edges 尾段前瞻、以及其上游全部
// 具名引用)原样保留于下方,未改动其代码正文。
// [2] duck.ac 题目页 <https://duck.ac/problem/1013>(站点 FAQ 明示提交正文公开可见)。
// LICENSE:duck.ac 站内本账号公开提交;本改动为本账号自产,未取自任何外部来源、无第三方代码/资料、无外部仓库/论文。
// ===== 思路区(本发:pass A AND 内核 pfx 前缀行 trow 的核内 +96 字前瞻)=====
// 【单变量】and_popcount_prefix<6> 主环:在 4 条池位图行的 +96 字前瞻旁,为 pfx 前缀行 trow(每查询一条
// 随机行,1 字节/字;行首仅由界面前瞻 2 行、距离 1 查询)增加 `_mm_prefetch(trow + w + 96, T0)`。
// prefetch 无架构状态、不读写内存,输出逐位不变。守卫由环条件 w+8≤Wr−96 保证地址 ≤ trow+Wr−8(行内)。
// 【机理】trow 是 pass A 的第 5 条随机流(1B/字 ≈ Wr 字节/查询,n=1e6 均值 ~1.95KB = 31 行),核内按
// trow[w] / trow[w+4] 逐 8 字推进;行首 2 行由界面前瞻覆盖,其余依赖 HW 流预取器起转(每行一次性
// 随机行,HW 需数个 miss 才锁定)。在核内给 1.5 行(96 字 ≈ 12 次迭代 ≈ 数百拍)前瞻,使行内后续
// 行在 demand 前入 L1/L2。主环为带宽界、uop 余量 3–4×,每 8 字 1 条 prefetch 的代价被隐藏。
// 【判题机定价(experience n=300K 同源 ABBA 交错,2 轮,chk=710654885 逐位同)】
// 底盘(=官方 #134917 正文)441.206878 / 441.085017 ms;本件 438.066752 / 438.051509 ms
// ⇒ −3.140 / −3.034 ms(−0.71% / −0.69%,2/2 同向);同源 in-run rdtsc −0.71%/−0.69%。
// 【值闸门】本地 vendor gcc9 -O2 -static:n=300000 checksum 706398984、n=1000000 checksum 7842486998
// (与底盘逐位同);判题 experience n=300K checksum 710654885(两臂全同)。
// 【预登记】按查询数外推 n=1e6 ≈ −10 ms 量级(同 edges 刀的 300K→1e6 兑现率 ~1.5×,实际以官方 1 ms 分辨率为准)。
// 严阈 3,793,681µs;best 3,808,000µs(#134917)。
// ============================================================
// ===== 引用区(本发增量;题 1013)=====
// [13] 本账号 saffah_v41f1017_agg4,提交 #134818 <https://duck.ac/submission/134818>(题 1013,3.820 s,Accepted)
// —— 本发底盘:下方代码正文 = 该提交正文的逐字节副本 + 下述一处值中性簿记 store 删除;其自带上游引用链保留在下方、未改动。
// [14] duck.ac 题目页 <https://duck.ac/problem/1013>(站点 FAQ 明示提交正文公开可见)。
// LICENSE:以上均为 duck.ac 站内本账号公开提交/实测;本发改动本账号自产,无第三方代码/资料、无外部仓库/论文来源。
// ===== 思路区(本发目的与改动)=====
// 【本发性质】真改动(非同码重发、值逐位精确):qg 簿记 bp2v 六槽指针填充中的死存储删除(store 合并角)。
// 【依据】pfx_live 路径调用 and_popcount_prefix<6>;该内核全部主环/尾环/bsNext 预取只读 bs[0] 与 bs[2..5]
// (源码循环 d=2..KK−1),1 号槽(pfx 维的位图行)已被字节前缀平面替代——setup 环同一约定已见
// `if(pfx_live && e==1)continue;`(本文件 1715 行处)。但组首填充与每查询 next-setup 仍无条件计算并存入
// 1 号槽(s_pofs 随机载 + 地址 lea + 栈存),该槽在内核侧从不被读,属死存储。
// 【改动】两处填充的 e 环首各加 `if (pfx_live && e == 1) continue;`(g_plen_on 的两条直路各一处,共 4 处);
// pfx_live=false 时原式不变(range 内核使用 1 号槽),数组其余槽、调用点、算式全部不变,值逐位同。
// 【值闸门】本机 vendor gcc9 9.3.0 -O2 -std=gnu++17:n=300000 checksum 706398984、n=1000000 checksum 7842486998
// (与底盘 #134818 逐位同)。
// 【判靶】best 3,820,000µs(#134818/#134823 平值);T=#131237 3,832,000µs(早于我方 ⇒ 严支,严阈 3,793,681µs)。
// 本发为簿记桶 store 删除角的官方单发定价件(board 取 min,慢发不伤 best)。
// ============================================================
// ===== 引用区(本发增量)=====
// [11] 本账号 saffah_v41f1017_agg4,提交 #134789 <https://duck.ac/submission/134789>(题 1013,3.822 s,Accepted)
// —— 本发底盘:本文件代码正文 = 该提交正文的逐字节副本 + 下述两处自产重构;该件自带的上游引用链保留在下方,未改动其代码正文。
// 所用内容:整份引擎正文(qg 每查询簿记段)。
// [12] duck.ac 题目页 <https://duck.ac/problem/1013>(站点 FAQ 明示提交正文公开可见)。
// LICENSE:以上均为 duck.ac 站内本账号公开提交;本发两处改动均本账号自产,无第三方代码/资料、无外部仓库/论文来源。
// ===== 思路区 =====
// 【本发性质】真改动(非同码重发):在 #134789 底盘上对 qg 每查询簿记做两处值精确结构重构。
// 定位依据:本账号 2026-10-18 16:0x 判题 experience n=300K 相位仪表 + MF_SKIP_AND 差分:
// qg 簿记 = 21.4M 拍(1.34% of total、71 拍/查询),其中 bp2v 六平面指针 setup 的 g_plen_on 三元
// 与 pfx 偏移查表为结构成本(非预取删除)。
// 【改动 1|g_plen_on 三元外提】next-setup 的 e 环原每查询 6 次运行时三元 `g_plen_on ? s_pofs 表 :
// (e*(Q+1)+bdn)*W`;现每查询一次判定 + 两个 6 迭代直路,乘法臂预计算 q1w=(Q+1)*W 并把 e*q1w 提到
// 循环外。组首 bp2v 初始化同样外提。地址算式 u64 精确等价 ⇒ 输出逐位同。
// 【改动 2|pfx 偏移跨迭代携带】AND 调用处原重新取 bd=s_bd+i*MAXK 与 pfx_offsets[bd[pfx_dim]],
// 而同一值已在前一查询的 next-setup 中为预取行算出 ⇒ 现以 pfx_off_cur 寄存器携带(组首查询就地
// 初始化、每迭代末结转 pfx_off_nxt)。每查询省一条随机 s_bd 行载与一次 pfx_offsets 查表;
// pfx_offsets 在查询环内无写,值恒同。
// 【值闸门】本地 g++:n=300000 checksum 706398984、n=1000000 checksum 7842486998(与底盘逐位同)。
// 【判靶】严阈 3,793,681µs(T=#131237 3,832,000µs,早于我方 ⇒ 严支);best #134789 = 3,822,000µs;
// 本发为真实结构改动的官方单发定价件(board 取 min,慢发不伤 best)。
// ============================================================
// ===== 引用区(本发增量)=====
// [9] 本账号 saffah_v41f1017_agg4,提交 #134780 <https://duck.ac/submission/134780>(题 1013,3.823 s,Accepted)
// —— 本发底盘 = 下方 #134778 逐字节正文 + 其已含的 MF_PFB 6→8(即 #134780 正文);其自带的上游引用链保留在下方,未改动。
// [10] 以下三处改动均为本账号自产、判题 experience 实测的亚阈正项(无第三方来源):
// (a) s_ord 置换散播目的地写意图 hint:本账号 2026-10-18 11:4x 台账,判题 n=300K 两轮交错
// so_base T0@32 448.372/448.392 ms vs so_pw32 448.235/448.160 ms(均值 −0.186 ms,2/2 同向);
// (b) a13_tie_flush 滚动预取:本账号 2026-10-18 11:0x 台账,判题 n=300K 两轮
// D=32 −0.18/−0.11 ms(同向;D=64 −0.14/−0.07);
// (c) s_pm 零环瘦身 d=k..8:本账号 2026-10-18 12:4x 台账,判题 n=300K 三轮 −0.060 ms 均值(3/3 同向);
// 同式已由本账号提交 #134787 官方单发(3.826 s, Accepted)。
// 「1008 律」原始出处:本账号 saffah_v41f1017_agg4,提交 #134685(题 1008)—— 目的地 RFO 预取须用 prefetchw;
// 以上均为 duck.ac 站内本账号公开提交/实测,无第三方许可证声明。
// ===== 思路区(本发目的与改动)=====
// 【本发性质】真改动(非同码重发):在 #134780 底盘上合并三处值中性、判题同向为正的微项。
// 【改动 1】s_ord 置换散播环(od[pp[i]]=i|(s_g[i]<<20))的 pp[i+MF_PFSC2]→od+qp 目的地提示:
// 只读 T0 → prefetchw 写意图(距离 MF_PFSC2=32 不变;目的地为 N 置换、全稠密,每存一次 RFO miss)。
// 【改动 2】a13_tie_flush 滚动预取(D=32 项前瞻):首环对 q+32 的 s_ord[e2[2]]+e2[1] 行发 T0;
// 首环原有 s_pm 行预取收紧为首 32 项;次环对 q+32 的 s_tj 行发 s_pm T0(q+16 out prefetchw 原样保留)。
// 【改动 3】s_pm 填充零环 d=k..STRIDE → d=k..8:p[8..15] 随后被 BVAL/0x7fffffffu 两段全覆盖,
// 行内容逐位同;每行省 8 条冗余零存(判题 n=300K 实测 −0.013%)。
// 【值闸门】本机 n=300000 checksum 706398984(全输出 FNV 12338085432474361989)、
// n=1000000 checksum 7842486998(FNV 17496218108781718547),base 与合并件逐位同;
// 判题 experience n=300000 checksum 710654885(两臂 6 发全同)。
// 【判题机 A/B】n=300K 同源驱动,3 轮交错(正/反/正),每臂 1 发:
// base 447.525218/447.565117/447.457622 ms vs 合并 447.403745/447.327763/447.397057 ms
// ⇒ 3/3 同向,逐轮 −0.121/−0.237/−0.061 ms(Rmed −0.029%)。
// 【判靶】严阈 3,793,681µs(T=#131237 3,832,000µs,早于我方 ⇒ 严支);best #134780 = 3,823,000µs;
// 本发为真实改进的官方单发测量件,兑现率由官方读数标定(board 取 min)。
// ============================================================
// ===== 引用区(本发增量;本文件代码正文 = #134778 的逐字节副本,仅一个宏值不同)=====
// [7] 本账号 saffah_v41f1017_agg4,提交 #134778 <https://duck.ac/submission/134778>(题 1013,3.825 s,Accepted)
// —— 本发底盘(正文逐字节副本;其自带的上游引用链/署名保留在下方,未改动其代码正文)。
// [8] 本账号 saffah_v41f1017_agg4,提交 #134749 / #134734(题 1013,3.826 / 3.828 s)—— 底盘所含的 tie-flush
// 未来项 prefetchw@16 与 pfx_pos 目的地 prefetchw@8 两件,引用链见下方原有区块。
// ===== 思路区(本发目的与改动)=====
// 【本发性质】单变量真改动:fringe 相位的未来项预取距离 MF_PFB 6 → 8(底盘自身常数,无外部来源)。
// 【机理】fringe_dim 的每查询在 `#if MF_PFB>0` 块对 qi+MF_PFB 项的 s_pm 记录、s_gt[d]、out(prefetchw)发三提示;
// 本底盘各查询的 fringe 扫描长度 ≈ 桶内 rank 跨度(1e6: B=713),6→8 项的前瞻更贴合 s_pm 记录的 DRAM 取数延迟,
// 且不明显挤占签名平面流的 L1/L2 窗口。纯提示,语义零改动。
// 【值闸门】判题机 experience n=300000 checksum 710654885(base/f8/f12/f0/a0/s0 六臂逐位同);
// 本地 n=300000 checksum 706398984、n=1000000 checksum 7842486998,与底盘同。
// 【判题机 A/B】n=300K 同源驱动,三轮交错(正/反/正序),每臂 3 发:
// base(#134778 正文) 448.25/448.19/448.14 ms;本发(D=8) 447.57/447.51/447.61 ms ⇒ 三轮同向 −0.64~−0.68ms(−0.15%);
// D=12 对照 447.65/447.74/447.64;删除臂 MF_PFB=0 = 561.6~561.8ms(+25%)⇒ 站点承重;MF_PFA=0 = 457.0ms(+2%)承重;
// pfx 源流预取删除 448.35~448.41(+0.03%,略负)⇒ 现役常数均正确。
// 【判靶】严阈 3,793,681µs(T=#131237 3,832,000µs,早于我方 ⇒ 严支);best #134778 = 3,825,000µs;
// 本发为真实改进的官方单发测量件,兑现率由官方读数标定。
// ============================================================
// ===== 引用区(本发增量;本文件代码正文 = #134749 的逐字节副本,仅下述一处预取块不同)=====
// [5] 本账号 saffah_v41f1017_agg4,提交 #134734 <https://duck.ac/submission/134734>(题 1013,3.828 s,Accepted)
// —— 本发合并的「pfx_pos 目的地 RFO 预取:距离 8 + prefetchw 写意图」即该件对 #134691 的单变量改动,逐字移植。
// [6] 本账号 saffah_v41f1017_agg4,提交 #134749 <https://duck.ac/submission/134749>(题 1013,3.826 s,Accepted)
// —— 本发底盘(正文逐字节副本;该件自带的上游引用链/署名保留在下方,未改动其代码正文)。
// 「1008 律」原始出处:本账号 saffah_v41f1017_agg4,提交 #134685(题 1008)—— 目的地 RFO 预取须用 prefetchw
// (写意图),只读提示无 RFO 语义;距离以 8 为佳。以上均为 duck.ac 站内公开提交,无第三方许可证声明。
// ===== 思路区(本发目的与改动)=====
// 【本发性质】真改动(非同码重发):把两件已各自官方验证的独立改动合并到同一正文。
// 【相对 #134749 的唯一改动】pfx 建表散播环的目的地提示:
// 旧(#134749 正文):if(j+32<limit){ ... _mm_prefetch(pfx_pos+id, T0); } —— 与源预取共用 id 的只读提示
// 新(本发):源预取原样保留;另加 if(j+8<limit){ ... __builtin_prefetch(pfx_pos+idd, 1, 3); }(= #134734 原式)
// 机理:pfx_pos 为 4MB 温表、每点一次随机 4B 存 = 每存一次 RFO miss;prefetchw 在距离 8(≈单点散播延迟)
// 把目的地行以独占态提前放入在飞窗口。
// 【值闸门】判题机 experience n=300000 checksum 710654885(base/合并/D16 三臂逐位同);本地 n=300000 checksum
// 706398984、n=1000000 checksum 7842486998,均与底盘同。
// 【判题机 A/B】n=300K 同源驱动,三轮交错(正/反/正序),每臂 3 发:base 448.839/448.792/448.800 ms;
// 合并臂(D=8) 448.136/448.112/448.197 ms ⇒ 三轮同向,均值 −0.67ms(−0.153%);D=16 对照 448.191/448.205/448.287。
// 【判靶】严阈 3,793,681µs(T=#131237 3,832,000µs,早于我方 ⇒ 严支);best #134749 = 3,826,000µs;
// 本发为真实改进的官方单发测量件,兑现率由官方读数标定。
// ============================================================
// ===== 引用区(REFERENCES)=====
// [1] duck.ac 账号 saffah_cc_v41_agg1,提交 #131237
// <https://duck.ac/submission/131237>(题 1013,3.832 s,Accepted)。
// —— 本文件代码正文(本引用区之后的第一条预处理指令起至文件末尾)= 该提交的代码正文副本;
// 本发相对它的改动见下方【思路区】。该件正文自带的上游引用链/署名在其原提交头部注释中,本发未改动其代码正文。
// [2] 本账号 saffah_v41f1017_agg4,提交 #133053 <https://duck.ac/submission/133053>
// (3.830 s,Accepted)—— 本件直接前驱(底盘 + MF_QCOEF 110→130);本发在其代码正文上仅加下述预取块。
// [3] 本账号 saffah_v41f1017_agg4,提交 #134673 <https://duck.ac/submission/134673>
// (题 1014a,456.396 ms,Accepted)—— 本发所移植的「下一查询行切片预取」思想来源:
// far 相把下一查询的 7 维 × 2×64B 行切片提前一个迭代放进在飞窗口(T0;距离 1 最优,
// 单行负、双行正;官方 −2.38%)。本发在 1013 的 pass-A 前缀(pfx)行上移植该思想,代码为本席自写。
// [4] duck.ac 题目页 <https://duck.ac/problem/1013>(站点 FAQ 明示:提交的代码将会被公开,所有人可见)。
// LICENSE:上述均为 duck.ac 站内公开提交,无第三方许可证声明、非来自任何外部仓库;
// 按「账号 + 原提交 URL + 所用内容」具名引用,不附加额外限制。本发未参考任何论文/外部开源仓库。
// ===== 思路区(本发目的与改动)=====
// 【本发性质】单变量真改动(落在 #133053 代码正文):pass-A 前缀行(pfx)的下一查询切片预取。
// 【机理】and_popcount_prefix 主环已用 bsNext 把下一查询的 4 条位图行按 +96 字距离推进;
// 但每个查询还读第 5 条随机流 trow(pfx 行 = pfx_rows+pfx_offsets[bd[pfx_dim]]),
// 其首行无前瞻(查询首迭代即 demand 读)。本发在构造下一查询 bp2v 的同一处
// (qq+1<qe 且 pfx_live)对下一查询 pfx 行的前两条 64B 行发 _mm_prefetch(T0),
// 把该流首个 demand miss 提前一个查询放进在飞窗口(1014a 距离 1 形态)。语义零改动。
// 【值闸门】判题同款 gcc 9.3.0:本地 n=300000 checksum 706398984、n=1000000 checksum 7842486998,
// 与底盘逐位同;判题机 experience n=300000 checksum 710654885(与底盘逐位同)。
// 【判题机 A/B】experience n=300K 同源交错:底盘 448.447/448.435 ms(漂移 0.012),
// 双行版 447.325/447.330 ms ⇒ −1.12 ms(−0.25%);单行版 447.809 ms(−0.64 ms);
// 4 行版 447.330 ms(与双行同)⇒「双行」即够,剂量方向与 1014a 一致。
// 【1013 第二席追发·相对 #134691 正文的唯一改动】a13_tie_flush 次环(out RMW 环)在项 q 处理前取 q+16 项的 e[0],
// 对 out+idx 发 exclusive prefetchw(destination RFO 前置,距离 16 项);值中性(prefetch 无架构状态。
// 依据:本账号 1008 提交 #134685 的「目的地 RFO 预取律」(写意图必须 prefetchw,只读提示会失去 RFO 语义;
// 前瞻邻点已扫);本相此前 out RFO 64B/项 无前瞻(首环同项预取提前量过大 ⇒ 中性),未来项目的地为未试轴。
// 值闸:本地 n=300000 checksum 706398984;判题 n=300000 checksum 与底盘逐位同(710654885)。
// 目的:判题 n=300K 同源 A/B 定价后空窗官方单发(∝n 兑现)。
// 【判靶】严阈 3,793,681µs(T=3,832,000;我方 best 3,830,000);本发为真实改进的官方单发测量件,
// n=300K→1e6 的兑现率由官方读数标定。
// ============================================================
#define MF_WTRIM 1
#define MF_NOSDMP 1
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing,schedule-insns,no-tree-pre")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#define MF_QCOEF 130
#define MF_NOFILTB 1
#define MF_SPL ((size_t)1000064)
#include <immintrin.h>
#include <type_traits>
#ifdef LOCAL_TEST
#include <cstdio>
#include <cstdlib>
#include <ctime>
#else
#ifdef __cplusplus
extern "C" void *malloc(unsigned long);
#else
void *malloc(unsigned long);
#endif
#endif
typedef unsigned u32;
typedef unsigned char u8;
typedef unsigned short u16;
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
#endif
#define MAXK 8
#define STRIDE 16
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 122
#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 1800
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
#ifndef MF_PFB
#define MF_PFB 8
#define MF_PFSC 32
#define MF_PFSC2 32
#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 2
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 1408
#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)));
#define E6Z_BDR (((size_t)(MAXN) + 31u) & ~(size_t)31u)
static u16 s_bdT[(size_t)MAXK * E6Z_BDR];
static u16 s_bd[(size_t)MAXN * MAXK];
static u8 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_ct[MAXN + 8] __attribute__((aligned(64)));
static u16 s_bkt[MAXN];
static u16 s_bks[MAXN];
static u32 s_nid[MAXN];
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN];
#ifndef MF_INCA
#define MF_INCA 1
#endif
#ifndef MF_INCA_PFD
#define MF_INCA_PFD 0
#endif
static int g_inca = 0;
static u32 *s_ho = 0;
static size_t s_ho_cap = 0;
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 int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
static u32 solve_generation;
#define TIE_CAP 524288
#define TIE_EW 12
static u32 s_tbuf[(size_t)TIE_CAP * TIE_EW];
static u32 s_tj[TIE_CAP];
static u32 s_tn = 0;
static int g_nooff = 0;
static u32 s_pofs[(size_t)MAXK * MAXQ];
static u32 s_plen[(size_t)MAXK * MAXQ];
static int g_plen_on = 0;
static u32 g_W, g_Q, g_B;
static u32 g_KS;
static int g_nosort = 0;
static int g_nosdmp = 0;
static u32 pfx_pos[MAXN+64] __attribute__((aligned(64)));
static u8 pfx_counts[MAXN/64+256] __attribute__((aligned(64)));
static u32 pfx_offsets[MAXQ+1];
static u8 *pfx_rows=0; static size_t pfx_rows_cap=0;
static u8 *pfx_edges=0; static size_t pfx_edges_cap=0;
static bool pfx_live=false; static int pfx_dim=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));
}
__attribute__((target("avx2,popcnt")))
static inline u32 hsum256(__m256i v) {
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i s = _mm_add_epi64(lo, hi);
return (u32)((u64)_mm_cvtsi128_si64(s) + (u64)_mm_extract_epi64(s, 1));
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_range(const u64 *const *bs, u32 Wr, u32 rem, const u64 *const *bsNext) {
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;
u32 wlim = (Wr > 96u) ? (Wr - 96u) : 0u;
for (; w + 8 <= wlim; w += 8) {
if (KK > 1) _mm_prefetch((const char *)(bs[1] + w + 96), _MM_HINT_T0);
if (KK > 2) _mm_prefetch((const char *)(bs[2] + w + 96), _MM_HINT_T0);
if (KK > 3) _mm_prefetch((const char *)(bs[3] + w + 96), _MM_HINT_T0);
if (KK > 4) _mm_prefetch((const char *)(bs[4] + w + 96), _MM_HINT_T0);
if (KK > 5) _mm_prefetch((const char *)(bs[5] + w + 96), _MM_HINT_T0);
if (KK > 6) _mm_prefetch((const char *)(bs[6] + w + 96), _MM_HINT_T0);
if (KK > 7) _mm_prefetch((const char *)(bs[7] + w + 96), _MM_HINT_T0);
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)), z));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < KK; d++) b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)), z));
}
if (bsNext) {
u32 off = 0;
for (; w + 8 <= Wr; w += 8, off += 8) {
if (KK > 1) _mm_prefetch((const char *)(bsNext[1] + off), _MM_HINT_T0);
if (KK > 2) _mm_prefetch((const char *)(bsNext[2] + off), _MM_HINT_T0);
if (KK > 3) _mm_prefetch((const char *)(bsNext[3] + off), _MM_HINT_T0);
if (KK > 4) _mm_prefetch((const char *)(bsNext[4] + off), _MM_HINT_T0);
if (KK > 5) _mm_prefetch((const char *)(bsNext[5] + off), _MM_HINT_T0);
if (KK > 6) _mm_prefetch((const char *)(bsNext[6] + off), _MM_HINT_T0);
if (KK > 7) _mm_prefetch((const char *)(bsNext[7] + off), _MM_HINT_T0);
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)), z));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < KK; d++) b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)), z));
}
}
for (; w + 8 <= Wr; w += 8) {
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)), z));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < KK; d++) b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)), z));
}
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 lo2 = _mm256_and_si256(a, lm);
__m256i hi2 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo2), _mm256_shuffle_epi8(lk, hi2)), z));
}
tot = _mm256_add_epi64(acc0, acc1);
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;
}
static inline __attribute__((always_inline,target("avx2"))) u64 edge_mask(const u8* row,u32 rem) {
const __m256i r=_mm256_set1_epi8((char)rem);
const u32 lo=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(r,_mm256_loadu_si256((const __m256i*)row)));
const u32 hi=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(r,_mm256_loadu_si256((const __m256i*)(row+32))));
return (u64)lo|((u64)hi<<32);
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_prefix(const u64 *const *bs, u32 Wr, u32 rem, const u64 *const *bsNext, const u8* trow, const u8* edges) {
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);
__m256i tot = z, acc0 = z, acc1 = z;
u32 w = 0, cnt = 0, total = 0;
u32 wlim = (Wr > 96u) ? (Wr - 96u) : 0u;
for (; w + 8 <= wlim; w += 8) {
if (KK > 2) _mm_prefetch((const char *)(bs[2] + w + 96), _MM_HINT_T0);
_mm_prefetch((const char *)(trow + w + 96), _MM_HINT_T0);
if (KK > 3) _mm_prefetch((const char *)(bs[3] + w + 96), _MM_HINT_T0);
if (KK > 4) _mm_prefetch((const char *)(bs[4] + w + 96), _MM_HINT_T0);
if (KK > 5) _mm_prefetch((const char *)(bs[5] + w + 96), _MM_HINT_T0);
if (KK > 6) _mm_prefetch((const char *)(bs[6] + w + 96), _MM_HINT_T0);
if (KK > 7) _mm_prefetch((const char *)(bs[7] + w + 96), _MM_HINT_T0);
__m256i a = LDU256((bs[0] + w));
for (int d = 2; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
{ __m128i t=_mm_cvtsi32_si128(*(const u32*)(trow+w)); a=_mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(t)),ones)); }
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)), z));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 2; d < KK; d++) b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
{ __m128i t=_mm_cvtsi32_si128(*(const u32*)(trow+w+4)); b=_mm256_and_si256(b,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(t)),ones)); }
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)), z));
}
if (rem) _mm_prefetch((const char *)edges, _MM_HINT_T0);
if (bsNext) {
u32 off = 0;
for (; w + 8 <= Wr; w += 8, off += 8) {
if (KK > 2) _mm_prefetch((const char *)(bsNext[2] + off), _MM_HINT_T0);
if (KK > 3) _mm_prefetch((const char *)(bsNext[3] + off), _MM_HINT_T0);
if (KK > 4) _mm_prefetch((const char *)(bsNext[4] + off), _MM_HINT_T0);
if (KK > 5) _mm_prefetch((const char *)(bsNext[5] + off), _MM_HINT_T0);
if (KK > 6) _mm_prefetch((const char *)(bsNext[6] + off), _MM_HINT_T0);
if (KK > 7) _mm_prefetch((const char *)(bsNext[7] + off), _MM_HINT_T0);
__m256i a = LDU256((bs[0] + w));
for (int d = 2; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
{ __m128i t=_mm_cvtsi32_si128(*(const u32*)(trow+w)); a=_mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(t)),ones)); }
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)), z));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 2; d < KK; d++) b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
{ __m128i t=_mm_cvtsi32_si128(*(const u32*)(trow+w+4)); b=_mm256_and_si256(b,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(t)),ones)); }
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)), z));
}
}
for (; w + 8 <= Wr; w += 8) {
__m256i a = LDU256((bs[0] + w));
for (int d = 2; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
{ __m128i t=_mm_cvtsi32_si128(*(const u32*)(trow+w)); a=_mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(t)),ones)); }
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)), z));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 2; d < KK; d++) b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
{ __m128i t=_mm_cvtsi32_si128(*(const u32*)(trow+w+4)); b=_mm256_and_si256(b,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(t)),ones)); }
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)), z));
}
for (; w + 4 <= Wr; w += 4) {
__m256i a = LDU256((bs[0] + w));
for (int d = 2; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
{ __m128i t=_mm_cvtsi32_si128(*(const u32*)(trow+w)); a=_mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(t)),ones)); }
__m256i lo2 = _mm256_and_si256(a, lm);
__m256i hi2 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo2), _mm256_shuffle_epi8(lk, hi2)), z));
}
tot = _mm256_add_epi64(acc0, acc1);
total = hsum256(tot);
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 2; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v & (trow[w]>=64?~0ull:((1ull<<trow[w])-1)));
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 2; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & (trow[Wr]>=64?~0ull:((1ull<<trow[Wr])-1)) & edge_mask(edges,rem));
}
return total;
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const int k = KK;
const int K1 = KK - 1;
const u32 ASTRL = STRIDE;
#if MF_ASCR
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
{
const u32* fo=s_frord[d1];
for(u32 qi=q0;qi<q1;qi++)s_qord[qi]=fo[qi]&0xFFFFFu;
}
#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 * E6Z_BDR + s_qord[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 = s_qord[qi]; s_tmp[cntv[s_bdT[(size_t)dd * E6Z_BDR + id] + 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 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 + ((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
{
#if MF_INCA
const u32 *a13_ho = 0;
u32 a13_app = 0;
if (g_inca) {
a13_ho = s_ho;
for (u32 w = 0; w < AS_WORDS; w++) s_ascr[w] = 0;
}
#endif
u32 qi = q0;
while (qi < q1) {
const u32 *fo=s_frord[d1];
u32 b1=fo[qi]>>20;
u32 qe=qi+1;
while(qe<q1 && (fo[qe]>>20)==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);
#if MF_INCA
if (a13_ho) {
u32 pto = (b1 < Q) ? FPOS(d1, b1) : (u32)N;
const u32 *hp = a13_ho;
#if MF_INCA_PFD > 0
while (a13_app < pto) {
if (a13_app + MF_INCA_PFD < pto)
_mm_prefetch((const char *)(s_ascr + (hp[a13_app + MF_INCA_PFD] >> 6)), _MM_HINT_T0);
u32 q = hp[a13_app]; s_ascr[q >> 6] |= 1ull << (q & 63); a13_app++; }
#else
while (a13_app < pto) { u32 q = hp[a13_app]; s_ascr[q >> 6] |= 1ull << (q & 63); a13_app++; }
#endif
bp1 = s_ascr;
} else
#endif
{
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 (nwd <= AS_WORDS) {
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
}
const u64 *bp2v[2][MAXK];
u32 Rcur = s_gt[ds][s_qord[qi]];
u32 pfx_off_cur = pfx_live ? pfx_offsets[s_bd[(size_t)s_qord[qi] * MAXK + pfx_dim]] : 0u;
u32 pfx_off_nxt = 0u;
{
const u16 *bd0 = s_bd + (size_t)s_qord[qi] * MAXK;
if (g_plen_on) {
for (int e = 0; e < K1; e++) {
if (pfx_live && e == 1) continue; /* [ow] prefix kernel skips slot 1 */
if (e == (int)e1) bp2v[qi & 1u][e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp2v[qi & 1u][e] = s_bs + s_pofs[(size_t)e * MAXQ + bd0[d]]; }
}
} else {
const size_t q1w = (size_t)(Q + 1) * (size_t)W;
for (int e = 0; e < K1; e++) {
if (pfx_live && e == 1) continue; /* [ow] prefix kernel skips slot 1 */
if (e == (int)e1) bp2v[qi & 1u][e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp2v[qi & 1u][e] = s_bs + (size_t)e * q1w + (size_t)bd0[d] * (size_t)W; }
}
}
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = s_qord[qq];
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(out + j)));
}
}
#endif
if (qq + 1 < qe) {
const u16 *bdn = s_bd + (size_t)s_qord[qq + 1] * MAXK;
const u64 **bn = bp2v[(qq + 1) & 1u];
if (g_plen_on) {
for (int e = 0; e < K1; e++) {
if (pfx_live && e == 1) continue; /* [ow] prefix kernel skips slot 1 */
if (e == (int)e1) bn[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bn[e] = s_bs + s_pofs[(size_t)e * MAXQ + bdn[d]]; }
}
} else {
const size_t q1w = (size_t)(Q + 1) * (size_t)W;
for (int e = 0; e < K1; e++) {
if (pfx_live && e == 1) continue; /* [ow] prefix kernel skips slot 1 */
if (e == (int)e1) bn[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bn[e] = s_bs + (size_t)e * q1w + (size_t)bdn[d] * (size_t)W; }
}
}
// Next-query prefix-row slice [4]: bsNext above advances the 4 bitmap
// rows of qq+1; the pfx row read at trow+w is a separate random stream,
// so touch its first two 64B lines one query early (distance 1).
if (pfx_live) {
u32 off = pfx_offsets[bdn[pfx_dim]];
const char *pn = (const char *)(pfx_rows + off);
_mm_prefetch(pn, _MM_HINT_T0);
_mm_prefetch(pn + 64, _MM_HINT_T0);
pfx_off_nxt = off;
}
}
const u32 R = Rcur;
if (qq + 1 < qe) Rcur = s_gt[ds][s_qord[qq + 1]];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64*const*next=(qq+1<qe)?bp2v[(qq+1)&1u]:0;
if(pfx_live){
cnt=and_popcount_prefix<KK-1>(bp2v[qq&1u],R>>6,R&63u,next,pfx_rows+pfx_off_cur,pfx_edges+(size_t)(R>>6)*64);
}else cnt=and_popcount_range<KK-1>(bp2v[qq&1u],R>>6,R&63u,next);
}
out[i] = cnt;
if (pfx_live) pfx_off_cur = pfx_off_nxt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
const u16 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
_mm_prefetch((const char *)(out + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = s_qord[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 + ((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, 0);
}
out[i] = cnt;
}
#endif
if (g_skip_fringe) return;
const u32 *const avpool=s_pm+8;
auto fringe_dim = [&](auto dim) {
constexpr int d=decltype(dim)::value;
if(d==ds) return;
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 u8 *const sbase = s_sig + (size_t)d * (k-1) * MF_SPL;
const u32 *const odd = s_ord[d];
for (u32 qi = q0; qi < q1; qi++) {
u32 i = fo[qi] & 0xFFFFFu;
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB;
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);
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(out + j)));
}
}
#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 u32 *Avp = avpool + (size_t)i * ASTRL;
__m256i qr = _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
__m256i qT=qr;
if constexpr(d>1) qT=_mm256_min_epi32(qr,_mm256_or_si256(LDU256(Avp),notexlo));
else if constexpr(d==1) { if(ds!=0) qT=_mm256_min_epi32(qr,_mm256_or_si256(LDU256(Avp),notexlo)); }
__m256i qTh;
if (k > 9) {
__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 if (k > 8) {
qTh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
} else {
qTh = all8;
}
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (size_t)lo * g_KS) : 0;
if (g_use_filt && KK == 9) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
const u64 wcs = (u64)N / g_Qc + 1;
const u32 FW = g_fbm_W;
const u64 *rows[8];
int nr = 0;
for (int e = 0; e < 9; e++) {
if (e == d) continue;
u32 cell = (u32)((u64)rec[e] / wcs) + 1;
if (cell > g_Qc) cell = g_Qc;
rows[nr++] = s_fbm + ((size_t)(d * 9 + e) * (g_Qc + 1) + cell) * FW + w0;
}
const u64 *r0 = rows[0];
const u64 *r1 = rows[1];
const u64 *r2 = rows[2];
const u64 *r3 = rows[3];
const u64 *r4 = rows[4];
const u64 *r5 = rows[5];
const u64 *r6 = rows[6];
const u64 *r7 = rows[7];
const u32 remh = hi & 63u;
{
u64 vv = r0[0] & r1[0] & r2[0] & r3[0] & r4[0] & r5[0] & r6[0] & r7[0];
vv &= ~0ull << (lo & 63);
if (nw == 1u && remh) vv &= (1ull << remh) - 1;
u32 word = w0;
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = (word << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * g_KS) : (s_pm + (size_t)(s_ord[d][j] & 0xFFFFFu) * STRIDE);
__m256i pv = LDU256(cp);
__m256i pvh = LDU256(cp + 8);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
add += ok;
}
}
for (u32 w = 1; w + 1 < nw; w++) {
u64 vv = r0[w] & r1[w] & r2[w] & r3[w] & r4[w] & r5[w] & r6[w] & r7[w];
u32 word = w0 + w;
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = (word << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * g_KS) : (s_pm + (size_t)(s_ord[d][j] & 0xFFFFFu) * STRIDE);
__m256i pv = LDU256(cp);
__m256i pvh = LDU256(cp + 8);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
add += ok;
}
}
if (nw > 1u) {
u32 w = nw - 1u;
u64 vv = r0[w] & r1[w] & r2[w] & r3[w] & r4[w] & r5[w] & r6[w] & r7[w];
if (remh) vv &= (1ull << remh) - 1;
u32 word = w0 + w;
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = (word << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * g_KS) : (s_pm + (size_t)(s_ord[d][j] & 0xFFFFFu) * STRIDE);
__m256i pv = LDU256(cp);
__m256i pvh = LDU256(cp + 8);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
add += ok;
}
}
out[i] += add;
continue;
}
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const u32 FW = g_fbm_W;
int first = 1;
#if MF_FR_NOREAD
first = 0;
for (u32 w = 0; w < nw; w++) fb[w] = 0;
#endif
for (int e = 0; e < k; e++) {
if (e == d) continue;
u32 c = (u32)((u64)rec[e] / wcs) + 1;
if (c > g_Qc) c = g_Qc;
const u64 *row = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
if (first) { for (u32 w = 0; w < nw; w++) fb[w] = row[w]; first = 0; }
else { for (u32 w = 0; w < nw; w++) fb[w] &= row[w]; }
}
fb[0] &= ~0ull << (lo & 63);
u32 r2 = hi & 63u;
if (r2) fb[nw - 1] &= (1ull << r2) - 1;
for (u32 w = 0; w < nw; w++) {
u64 vv = fb[w];
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = ((w0 + w) << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * g_KS) : (s_pm + (size_t)(s_ord[d][j] & 0xFFFFFu) * STRIDE);
__m256i pv = LDU256(cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (g_sig_on && (hi - lo) >= 2) {
__m256i qsb[16];
{
const __m256i zero8 = _mm256_setzero_si256();
const __m128i shv = _mm_cvtsi32_si128((int)g_sigsh);
const __m256i bias8 = _mm256_set1_epi8((char)0x80);
const __m256i sv_lo = _mm256_srl_epi32(qT, shv);
const __m128i u16=_mm_packus_epi32(_mm256_castsi256_si128(sv_lo),_mm256_extracti128_si256(sv_lo,1));
const __m128i u8s=_mm_packus_epi16(u16,u16);
const __m256i packed=_mm256_xor_si256(_mm256_broadcastsi128_si256(u8s),bias8);
for (int e = 0; e < k; e++) {
if(e==(int)d)continue;
qsb[e-(e>(int)d)]=_mm256_shuffle_epi8(packed,_mm256_set1_epi8((char)e));
}
}
u32 aj = lo, aacc = 0, aexact = 0;
if (lo + 32u <= hi) {
const u32 ax_last = hi - 32u;
const int ax_e0 = ((int)d == 0) ? 1 : 0;
do {
__m256i ac=_mm256_subs_epi8(LDU256(sbase+aj),qsb[0]);
for(int e=1;e<k-1;e++){
__m256i cv=LDU256(sbase+(size_t)e*MF_SPL+aj);
ac=_mm256_max_epi8(ac,_mm256_subs_epi8(cv,qsb[e]));
}
u32 decm = (u32)_mm256_movemask_epi8(ac);
aacc += (u32)__builtin_popcount(decm);
{
__m256i eqz = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 amb=(u32)_mm256_movemask_epi8(eqz);
if (amb) {
if (k == 7 and !c0 && !g_nooff && s_tn + 32u <= (u32)TIE_CAP) {
u32 *te = s_tbuf + (size_t)s_tn * TIE_EW;
do {
int t = __builtin_ctz(amb); amb &= amb - 1;
te[0] = i; te[1] = aj + (u32)t; te[2] = (u32)d;
_mm256_storeu_si256((__m256i *)(te + 3), qT);
te += TIE_EW; s_tn++;
} while (amb);
} else {
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 ? (c0 + (size_t)(aj + (u32)t - lo) * g_KS) : (s_pm + (size_t)(odd[aj + (u32)t] & 0xFFFFFu) * STRIDE);
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
}
}
aj += 32u;
} while (aj <= ax_last);
}
{
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;
const int ax_e0t = ((int)d == 0) ? 1 : 0;
__m256i ac=_mm256_subs_epi8(LDU256(sbase+base),qsb[0]);
for(int e=1;e<k-1;e++){
__m256i cv=LDU256(sbase+(size_t)e*MF_SPL+base);
ac=_mm256_max_epi8(ac,_mm256_subs_epi8(cv,qsb[e]));
}
u32 decm = ((u32)_mm256_movemask_epi8(ac)) & kmask;
__m256i ax_eqz = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 amb = ((u32)_mm256_movemask_epi8(ax_eqz)) & kmask;
aacc += (u32)__builtin_popcount(decm);
if (amb) {
if (k == 7 && !c0 && !g_nooff && s_tn + 32u <= (u32)TIE_CAP) {
u32 *te = s_tbuf + (size_t)s_tn * TIE_EW;
do {
int t = __builtin_ctz(amb); amb &= amb - 1;
te[0] = i; te[1] = base + (u32)t; te[2] = (u32)d;
_mm256_storeu_si256((__m256i *)(te + 3), qT);
te += TIE_EW; s_tn++;
} while (amb);
} else {
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 ? (c0 + (size_t)(base + (u32)t - lo) * g_KS) : (s_pm + (size_t)(odd[base + (u32)t] & 0xFFFFFu) * STRIDE);
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
}
}
}
add = aacc;
} else if (c0) {
if (k <= 8) {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + 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] & 0xFFFFFu) * STRIDE;
__m256i pv = _mm256_load_si256((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = _mm256_load_si256((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
};
if constexpr(KK>0) fringe_dim(std::integral_constant<int,0>{});
if constexpr(KK>1) fringe_dim(std::integral_constant<int,1>{});
if constexpr(KK>2) fringe_dim(std::integral_constant<int,2>{});
if constexpr(KK>3) fringe_dim(std::integral_constant<int,3>{});
if constexpr(KK>4) fringe_dim(std::integral_constant<int,4>{});
if constexpr(KK>5) fringe_dim(std::integral_constant<int,5>{});
if constexpr(KK>6) fringe_dim(std::integral_constant<int,6>{});
if constexpr(KK>7) fringe_dim(std::integral_constant<int,7>{});
if constexpr(KK>8) fringe_dim(std::integral_constant<int,8>{});
if constexpr(KK>9) fringe_dim(std::integral_constant<int,9>{});
}
__attribute__((target("avx2,popcnt")))
static void a13_tie_flush(u32 *out) {
u32 n = s_tn;
if (!n) return;
s_tn = 0;
const __m256i all8 = _mm256_set1_epi32(-1);
u32 b0 = 0;
while (b0 < n) {
u32 b1 = b0 + 256u; if (b1 > n) b1 = n;
for (u32 q = b0; q < b1; q++) {
const u32 *e = s_tbuf + (size_t)q * TIE_EW;
if (q + 32 < b1) { const u32 *e2 = s_tbuf + (size_t)(q + 32) * TIE_EW; _mm_prefetch((const char *)(s_ord[e2[2]] + e2[1]), _MM_HINT_T0); }
u32 j = s_ord[e[2]][e[1]] & 0xFFFFFu;
s_tj[q] = j;
if (q < b0 + 32) _mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
}
for (u32 q = b0; q < b1; q++) {
if (q + 32 < b1) _mm_prefetch((const char *)(s_pm + (size_t)s_tj[q + 32] * STRIDE), _MM_HINT_T0);
if (q + 16u < b1) { const u32 *ef = s_tbuf + (size_t)(q + 16u) * TIE_EW; __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(out + ef[0]))); }
const u32 *e = s_tbuf + (size_t)q * TIE_EW;
const u32 *cp = s_pm + (size_t)s_tj[q] * STRIDE;
__m256i pv = LDU256(cp);
__m256i qv = _mm256_loadu_si256((const __m256i *)(e + 3));
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qv, pv), all8);
out[e[0]] += ok;
}
b0 = b1;
}
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int k) {
++solve_generation;
s_tn = 0;
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;
#if !MF_NOFILTB
overhead += (u64)k * k * 4 * N;
#endif
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 (k == 9 && N >= 900000) q = 1400;
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
if (q > (u32)(MAXQ - 1)) q = (u32)(MAXQ - 1);
g_Q = q;
g_B = (N + q - 1) / q;
}
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];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
static u8 hist8[MAXN+1] __attribute__((aligned(64)));
u32 *ct=s_ct;
for(u32 i=0;i<=N;i++)hist8[i]=0;
for(u32 i=0;i<N;i++){u32 v=xd[i];hist8[v<N?v:0]++;}
u32 sum=0;
for(u32 i=0;i<=N;i++){u32 c=hist8[i];ct[i]=sum;sum+=c;}
int pk=(sum==N);
if (!pk) {
for (u32 i = 0; i <= N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i <= N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tmp[i] = s; s += c; } }
{
u32 *pp = s_pos[d];
for (u32 i = 0; i < N; i++) {
#if MF_PFSC2 > 0
if (i + MF_PFSC2 < N) { u32 vp = xd[i + MF_PFSC2]; _mm_prefetch((const char *)(s_tmp + (vp < N ? vp : 0u)), _MM_HINT_T0); }
#endif
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; pp[i] = q;
}
}
{ u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
#if MF_PFSC > 0
if (i + MF_PFSC < N) { u32 vp = xd[i + MF_PFSC]; _mm_prefetch((const char *)(s_cntv + (vp < N ? vp : 0u)), _MM_HINT_T0); }
#endif
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);
u32 bbq = bb > Q - 1 ? Q - 1 : bb;
s_bdT[(size_t)d * E6Z_BDR + i] = bbq;
} }
#if MF_FB2
for (u32 i = 0; i < N; i++) {
#if MF_PFSC2 > 0
if (i + MF_PFSC2 < N) { u32 vp = xd[i + MF_PFSC2]; _mm_prefetch((const char *)(s_tmp + (vp < N ? vp : 0u)), _MM_HINT_T0); }
#endif
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
} else {
{
u32 *pp = s_pos[d];
u32 *gt = s_gt[d];
u32 *ctp = s_ct;
for (u32 i = 0; i < N; i++) {
#if MF_PFSC2 > 0
if (i + MF_PFSC2 < N) { u32 vp = xd[i + MF_PFSC2]; _mm_prefetch((const char *)(ctp + (vp < N ? vp : 0u)), _MM_HINT_T0); }
#endif
u32 v = xd[i]; if (v >= N) v = 0;
u32 w = ctp[v];
u32 st = w & 0xFFFFFu;
pp[i] = st + (w >> 20);
gt[i] = st;
u32 num = (ctp[v + 1] & 0xFFFFFu) - 1u; if (num > N) num = N;
u32 bb = (u32)(((u64)num * bmag) >> 40);
u32 bbq = bb > Q - 1 ? Q - 1 : bb;
s_bdT[(size_t)d * E6Z_BDR + i] = bbq;
ctp[v] = w + (1u << 20);
}
}
}
#endif
}
for (u32 i0 = 0; i0 < N; i0 += 64) {
u32 i1 = i0 + 64; if (i1 > N) i1 = N;
for (u32 i = i0; i < i1; i++) {
u16 *dst = s_bd + (size_t)i * MAXK;
for (int dd = 0; dd < k; dd++) dst[dd] = s_bdT[(size_t)dd * E6Z_BDR + i];
}
}
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
u32 gstart[MAXK+2]={};
{
static u32 y2_dsb[1024], y2_bst[1024];
for (u32 i0 = 0; i0 < N; i0 += 1024) {
u32 i1 = i0 + 1024; if (i1 > N) i1 = N;
for(u32 base=i0;base<i1;base+=8){
alignas(32)u32 folded[8];const u32 end=(base+8<i1)?base+8:i1;
if(end-base==8){
__m256i v=_mm256_slli_epi32(_mm256_loadu_si256((const __m256i*)(s_gt[0]+base)),4);
for(unsigned d=1;d<(unsigned)k;++d){
__m256i w=_mm256_or_si256(_mm256_slli_epi32(_mm256_loadu_si256((const __m256i*)(s_gt[d]+base)),4),_mm256_set1_epi32(d));v=_mm256_min_epu32(v,w);
}_mm256_store_si256((__m256i*)folded,v);
}else for(u32 i=base;i<end;++i){u32 v=s_gt[0][i]<<4;for(unsigned d=1;d<(unsigned)k;++d){u32 w=(s_gt[d][i]<<4)|d;if(w<v)v=w;}folded[i-base]=v;}
for(u32 i=base;i<end;++i){const u32 best=folded[i-base]>>4;const u32 ds=folded[i-base]&15u;++gstart[ds+1];
y2_dsb[i - i0] = (u32)ds;
y2_bst[i - i0] = best;
s_g[i] = (u32)ds;
}
}
if (g_plen_on) {
for (u32 i = i0; i < i1; i++) {
u32 ds = y2_dsb[i - i0];
u32 bb = y2_bst[i - i0];
const u16 *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
{ u32 *mp = mr + (size_t)d * MAXQ + b; u32 mv = *mp; *mp = (bb > mv) ? bb : mv; }
}
}
}
}
}
/* [a13as_ OD-GROUP] od[d][p] now carries the point's query group in its high bits; the bucket it
used to carry is rebuilt from p in K_frord. */
for (int d = 0; d < k; d++) {
u32 *od = s_ord[d];
const u32 *pp = s_pos[d];
for (u32 i = 0; i < N; i++) {
#if MF_PFSC2 > 0
if (i + MF_PFSC2 < N) { u32 qp = pp[i + MF_PFSC2]; __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(od + qp))); }
#endif
od[pp[i]] = i | ((u32)s_g[i] << 20);
}
}
for (int d = 0; d < k; d++) {
/* BVAL/FPOS moved out of the sort loop: it reads od, which loop 2 above fills. */
const unsigned *xd = x[d];
u32 *od = s_ord[d];
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] & 0xFFFFFu];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1] & 0xFFFFFu] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < 8; d++) p[d] = 0;
const u16*bd=s_bd+(size_t)i*MAXK;
for(int d=0;d<k-1;d++)p[8+d]=BVAL(d,bd[d]);
for(int d=k-1;d<8;d++)p[8+d]=0x7fffffffu;
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
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] & 0xFFFFFu] < 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] & 0xFFFFFu] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend] & 0xFFFFFu] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr] & 0xFFFFFu]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
for(int d=1;d<=k+1;++d)gstart[d]+=gstart[d-1];
{u32 cursor[MAXK];for(int d=0;d<k;++d)cursor[d]=gstart[d];for(u32 i=0;i<N;++i)s_qord[cursor[s_g[i]]++]=i;}
for (int d = 0; d < k; d++) {
u32 pos[MAXK];
for (int ds = 0; ds < k; ds++) pos[ds] = gstart[ds];
/* [a13as_ OD-GROUP] the group of the point at d-position r rides in od's high bits; the bucket
is rebuilt here as max{b : FPOS(d,b) <= r}. FPOS(d,b) is the first d-position of bucket b and
the d-order is sorted by value, so bucket(r) is monotone in r and one monotone cursor over the
FPOS table (Q steps per dim, amortised 0.0015 per element) recovers it exactly. */
const u32 *fpo = s_fpos + (size_t)d * MAXQ;
u32 bb = 0, nb = (Q > 1) ? fpo[1] : N;
u32 e_cur = ((Q > 1) ? (fpo[1] - 1u) : (N - 1u));
u32 bl = (u32)(((u64)e_cur * bmag) >> 40);
if (bl > Q - 1) bl = Q - 1;
for (u32 r = 0; r < N; r++) {
if (r >= nb) { /* class boundary: only ~N/B times per dim */
while (bb + 1 < Q && fpo[bb + 1] <= r) bb++;
e_cur = (bb + 1 < Q) ? (fpo[bb + 1] - 1u) : (N - 1u);
bl = (u32)(((u64)e_cur * bmag) >> 40);
if (bl > Q - 1) bl = Q - 1;
nb = (bb + 1 < Q) ? fpo[bb + 1] : N;
}
u32 x = s_ord[d][r];
u32 id = x & 0xFFFFFu;
u32 dst = pos[x >> 20]++;
s_frord[d][dst] = id | (bl << 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] & 0xFFFFFu) * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)(od[j] & 0xFFFFFu) * 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){u32 z_=row[e]>>g_sigsh;sb[(size_t)(e-(e>d))*MF_SPL+j]=(u8)(z_^0x80u);}
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
if (s_sig && !g_sig_merged && k==7) {
static u8 cached_signature[(size_t)MAXN*8] __attribute__((aligned(64)));
{const __m256i gather=_mm256_setr_epi8(0,4,8,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,4,8,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1);
const __m128i sh=_mm_cvtsi32_si128(g_sigsh);u32 i=0;
for(;i+8<=N;i+=8){__m128i c[8];c[7]=_mm_setzero_si128();
for(int e=0;e<7;++e){__m256i v=_mm256_srl_epi32(_mm256_loadu_si256((const __m256i*)(x[e]+i)),sh);v=_mm256_xor_si256(v,_mm256_set1_epi32(128));v=_mm256_shuffle_epi8(v,gather);c[e]=_mm_unpacklo_epi32(_mm256_castsi256_si128(v),_mm256_extracti128_si256(v,1));}
__m128i t0=_mm_unpacklo_epi8(c[0],c[1]),t1=_mm_unpacklo_epi8(c[2],c[3]),t2=_mm_unpacklo_epi8(c[4],c[5]),t3=_mm_unpacklo_epi8(c[6],c[7]);
__m128i u0=_mm_unpacklo_epi16(t0,t1),u1=_mm_unpackhi_epi16(t0,t1),u2=_mm_unpacklo_epi16(t2,t3),u3=_mm_unpackhi_epi16(t2,t3);
_mm_store_si128((__m128i*)(cached_signature+(size_t)i*8),_mm_unpacklo_epi32(u0,u2));
_mm_store_si128((__m128i*)(cached_signature+(size_t)(i+2)*8),_mm_unpackhi_epi32(u0,u2));
_mm_store_si128((__m128i*)(cached_signature+(size_t)(i+4)*8),_mm_unpacklo_epi32(u1,u3));
_mm_store_si128((__m128i*)(cached_signature+(size_t)(i+6)*8),_mm_unpackhi_epi32(u1,u3));
}
for(;i<N;++i)for(int e=0;e<7;e++)cached_signature[(size_t)i*8+e]=(u8)((x[e][i]>>g_sigsh)^128u);
}
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;auto emit=[&](auto tag){constexpr int OMIT=decltype(tag)::value;
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(cached_signature+(size_t)(od[j+64]&0xFFFFFu)*8),_MM_HINT_T0);
__m128i r0=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+0]&0xFFFFFu)*8));
__m128i r1=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+1]&0xFFFFFu)*8));
__m128i r2=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+2]&0xFFFFFu)*8));
__m128i r3=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+3]&0xFFFFFu)*8));
__m128i r4=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+4]&0xFFFFFu)*8));
__m128i r5=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+5]&0xFFFFFu)*8));
__m128i r6=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+6]&0xFFFFFu)*8));
__m128i r7=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)(od[j+7]&0xFFFFFu)*8));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__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 constexpr(OMIT!=0)_mm_storel_epi64((__m128i*)(base+(size_t)(0-(0>OMIT))*MF_SPL+j),v0);
if constexpr(OMIT!=1)_mm_storel_epi64((__m128i*)(base+(size_t)(1-(1>OMIT))*MF_SPL+j),_mm_srli_si128(v0,8));
if constexpr(OMIT!=2)_mm_storel_epi64((__m128i*)(base+(size_t)(2-(2>OMIT))*MF_SPL+j),v1);
if constexpr(OMIT!=3)_mm_storel_epi64((__m128i*)(base+(size_t)(3-(3>OMIT))*MF_SPL+j),_mm_srli_si128(v1,8));
if constexpr(OMIT!=4)_mm_storel_epi64((__m128i*)(base+(size_t)(4-(4>OMIT))*MF_SPL+j),v2);
if constexpr(OMIT!=5)_mm_storel_epi64((__m128i*)(base+(size_t)(5-(5>OMIT))*MF_SPL+j),_mm_srli_si128(v2,8));
if constexpr(OMIT!=6)_mm_storel_epi64((__m128i*)(base+(size_t)(6-(6>OMIT))*MF_SPL+j),v3);
}
};switch(d){case 0:emit(std::integral_constant<int,0>{});break;case 1:emit(std::integral_constant<int,1>{});break;case 2:emit(std::integral_constant<int,2>{});break;case 3:emit(std::integral_constant<int,3>{});break;case 4:emit(std::integral_constant<int,4>{});break;case 5:emit(std::integral_constant<int,5>{});break;case 6:emit(std::integral_constant<int,6>{});break;}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)(od[j]&0xFFFFFu)*8;
for(int e=0;e<k;e++)if(e!=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] & 0xFFFFFu) * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)(od[j] & 0xFFFFFu) * STRIDE;
for(int e=0;e<k;e++)if(e!=d){u32 z_=row[e]>>g_sigsh;base[(size_t)(e-(e>d))*MF_SPL+j]=(u8)(z_^0x80u);}
}
}
}
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;
#if MF_INCA
{
size_t needho = (size_t)N * sizeof(u32) + 4096;
if (needho > s_ho_cap) { s_ho = (u32 *)pool_alloc(needho); s_ho_cap = s_ho ? needho : 0; }
g_inca = (s_ho && (u64)N <= (u64)AS_WORDS * 64ull) ? 1 : 0;
}
#endif
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;
if(g_plen_on){
const u32 md=ds==0?1u:0u;const u32 *mr=s_maxr+((size_t)ds*MAXK+md)*MAXQ;
for(u32 b=0;b<Q;++b)if(mr[b]>mx)mx=mr[b];
}else 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
pfx_live=false;
pfx_dim=0; while(pfx_dim==ds || pfx_dim==((ds==0)?1:0))pfx_dim++;
#ifndef DISABLE_PREFIX_PORT
{
const u32 *mr=s_maxr+((size_t)ds*MAXK+pfx_dim)*MAXQ;
size_t bytes=0;
for(u32 b=0;b<=Q;b++){pfx_offsets[b]=bytes;u32 len=((mr[b]>>6)+64u)&~63u;bytes+=len;}
if(bytes>pfx_rows_cap){pfx_rows=(u8*)pool_alloc(bytes+4096);pfx_rows_cap=pfx_rows?bytes:0;}
size_t edge_words=((size_t)Wpool+1)*64;
if(edge_words>pfx_edges_cap){pfx_edges=(u8*)pool_alloc(edge_words+4096);pfx_edges_cap=pfx_edges?edge_words:0;}
if(pfx_rows && pfx_edges){
pfx_live=true;
for(u32 w=0;w<((N+63u)>>6)+128u;w++)pfx_counts[w]=0;
const u32*sdsp=s_pos[ds],*sdom=s_ord[pfx_dim];u32 j=0;
for(u32 b=0;b<=Q;b++){
u32 limit=b==0?0:(b<Q?FPOS(pfx_dim,b):N);
for(;j<limit;j++){
if(j+32<limit){u32 id=sdom[j+32]&0xfffffu;_mm_prefetch((const char*)(sdsp+id),_MM_HINT_T0);}
/* [PFX-M8] destination RFO (1008 law), distance 8, write intent: merged verbatim from our
official #134734 (standalone -1ms on the #134691 base). Source prefetch above unchanged. */
if(j+8<limit){u32 idd=sdom[j+8]&0xfffffu;__builtin_prefetch((const char*)(pfx_pos+idd),1,3);}
u32 id=sdom[j]&0xfffffu,r=sdsp[id],block=r>>6;
pfx_pos[id]=(block<<6)|pfx_counts[block]++;
}
u32 len=((mr[b]>>6)+64u)&~63u;u8*dst=pfx_rows+pfx_offsets[b];
for(u32 off=0;off<len;off+=64){
_mm256_stream_si256((__m256i*)(dst+off),_mm256_load_si256((const __m256i*)(pfx_counts+off)));
_mm256_stream_si256((__m256i*)(dst+off+32),_mm256_load_si256((const __m256i*)(pfx_counts+off+32)));
}
}
_mm_sfence();
const u32*od=s_ord[ds];
for(u32 block=0;block<(N+63)/64 && block<Wpool+1u;block++){
u32 begin=block*64,remain=N-begin;if(remain>64)remain=64;
u8*dst=pfx_edges+(size_t)block*64;
if(remain<64u)__builtin_memset(dst,127,64);
for(u32 t=0;t<remain;t++){
if(begin+t+16u<N){u32 future=od[begin+t+16u]&0xfffffu;_mm_prefetch((const char*)(pfx_pos+future),_MM_HINT_T0);}
u32 id=od[begin+t]&0xfffffu;dst[pfx_pos[id]&63]=(u8)t;
}
}
}
}
#endif
const u32 *const spp = pfx_live?pfx_pos:s_pos[ds];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
if(pfx_live && e==1)continue;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
#if MF_INCA
if (g_inca && d == ((ds == 0) ? 1 : 0)) {
u32 *dstho = s_ho;
const u32 *hp = spp;
for (u32 j = 0; j < N; j++) dstho[j] = hp[od[j] & 0xFFFFFu];
continue;
}
#endif
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 = spp[od[ptr] & 0xFFFFFu]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
#if MF_NT
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
query_group<7>(N, out, ds, q0, q1, Q, Wpool);
a13_tie_flush(out);
}
}
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 7); }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 3.734 s | 563 MB + 576 KB | Accepted | Score: 100 | 显示更多 |