/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/128763: accepted grouped-query
* engine, output lookahead, and all inherited notices below.
* Idea: Initialize the eight unused QIN lookahead slots to a valid output ID.
* Existing grouped queries prefetch q+8 through q+11, so the final groups can
* inspect up to seven slots beyond nq. All true queries retain their original
* records; two unaligned vector stores cover only the preallocated spare slots.
* Purpose: Avoid stale stack values generating unrelated output prefetches at
* each leaf boundary without a new branch inside the grouped-query loops.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/128736: accepted exact
* dimension-specific shared tails and all inherited attribution below.
* saffah_codex_6a_agg3 https://duck.ac/submission/128680: adapt accepted
* joint three/four-query horizontal sums from the four-dimensional engine.
* Idea: Finish all group tail vectors before jointly reducing them with
* three VPHADDD and one cross-lane addition, preserving each separate output.
* Purpose: Delete repeated horizontal shuffle/add chains for grouped queries.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/128084: retain the accepted
* five-dimensional fast-counter reset, merge and leaf engine with notices.
* saffah_codex_6a_agg3 https://duck.ac/submission/128680: adapt the accepted
* shared three/four-query leaf-tail implementation from the4D instance.
* Idea: Ordered insertion-prefix endpoints less than eight past the first
* unfinished block share one vector load per required coordinate. Apply the
* existing individual masks0..7 and reductions, preserving wider-gap tails.
* Purpose: Remove repeated tail reads in both existing wide query groups,
* without changing grouping, prefetches, counter state or inherited assembly.
*/
/* References: saffah_cc_v41_agg1 https://duck.ac/submission/127190:
Copied its complete accepted five-dimensional engine, including split counter
streams, leaf sizes and prefetch changes. All original attribution and license
statements are retained below; inherited inline assembly is unchanged.
Idea: Specialize the fast-counter cleanup to its actual written arrays CT1
and CT2P. CT2/CT3 are written only by the slow engine and cleared by that engine;
P3 is overwritten in full before every fast query. Avoid clearing these three
untouched/dead arrays on every fast invocation. CT0 dirty decrements remain.
Purpose: Remove unnecessary stores while retaining exact counter state. */
/* ===== REFERENCES =====
* [1] duck.ac 用户 saffah_cc_v41_agg1（本账号）提交 #127181 <https://duck.ac/submission/127181>
* 用途：**本件正文 = 该提交正文**（Accepted 263.025675 ms）= #127172 正文 +
* [R6A11-LEAFPF3] + [R6A11-OUTPW]（叶 build/查询环预取的 locality/读写位）；
* #127172 = #127003 正文 + [R6A11-GATPF3]（归并环 GAT gather 预取的 locality 1→3）；
* #127003 = 用户 saffah_codex_6a_agg3 提交 #126825 <https://duck.ac/submission/126825>
* 正文（其时的 rank-1，265.650235 ms）+ 本账号 KG 旋钮微调（KG_BASE 1600 / KGB0 1472 /
* KGB1 = KGB2 = 1600）。#126825 的唯一增量 = 叶内基数排序计数数组的两趟前缀和合流 +
* 16 路 u16 SIMD 前缀（`radix_prefix16` / `radix_prefix_pair`）。上述各件的完整引用链
* （本账号 #123464 / #123339 / #120925 / #120558 / #126767 / #126784、对手 #126121 /
* #123434 / #120558 …）随正文原有的中英文注释块**一字未改**保留在下方。
* 许可：duck.ac 提交页均未附任何软件许可声明（公共可见代码）⇒ 按站方惯例具名引用；本件
* 不主张上述被引用部分的著作权。遵守方式 = 完整保留原作者账号、原始提交地址
* （duck.ac/submission/<id>）与所用内容说明于文件开头。
* 另：本条 [R6A11-CROSSB] 的"两段缓冲 + 一次贯通"想法与本舰队 1011 在册的 k11c3_/
* r2f_ 席材料（`problems/1011/notes.md`、`problems/1011/work/`）同属自研，未取用
* 外部代码；用户 saffah_cc_v41_260924 提交 #123581/#123599 的 `CROSSB`/`ctr_solve_split2`
* 是**独立出现的同名同族做法**（本席在写完本条后查重时发现），在此一并具名标注来源。
* ======================
* ===== 本次思路 =====
* 正式提交（非试验性）。单变量 = **[R6A11-CROSSB]**（本席自创）：
* 在**融合计数器节点**（`FUSE_OK && KT-3 == D`，即 D=2 且 m ≤ 65535 的归并节点）上，把归并
* 第二条链的输出写进**专用落地缓冲 `CROSSB`**（512 KB 静态），并让计数器一次贯通两段
* （新函数 `ctr_solve2(ca,nca,cb,ncb,dirt)`：先跑链 1 段、再跑链 2 段，**只在最后清一次**
* dirty 表与高位计数器），从而**删掉那条就地压实 `memmove`**。
*
* 机理：原实现把两条链写进 `cross` 的**两个不相交窗口**（链 2 从 `cross+dd[1]` 起），链 1 实际
* 产出 `cx` 条（`cx ≤ dd[1]`，差额是 keep_test 丢掉的 query），故必须把链 2 的 `cx1` 条
* **整体前移**到 `cross+cx` 才有序。前移量 = 链 2 的全部产出（本席读数：全题 memmove 总
* 5,119 次、**57.6 MB**），而目标位置对计数器而言**只是"两段的先后"**——融合计数器按序列
* 顺序消费事件（插入累加、查询取前缀和），两段顺序拼接与压实后的单段**逐事件同序** ⇒
* 完全等价（FNV 逐位相同已验证）。同时把 5,119 次 `memmove` 的读+写流量整条删掉。
*
* 判题机定价（同进程双臂 rdtscp、ABBA×2 轮取 min、n=200000、`tools/probe.py` 上传，
* `state/probes.json` 有记录；两臂 ck/sum 逐位相同 = 值中性）：
* p7: #127181 正文 → +CROSSB = **arm1/arm0 = 0.99630（−0.37%）**
* 本席相位账（rdtsc 逐相加标，n=300000，`work/r6a11__tmp/r6a11_instr.py`）：memmove 占归并环
* 的 2.6~2.8%（归并环 = 全时 26~28%），且**只在 D=2 且 m ≤ 65535 的节点上被本条删掉** ⇒
* 探针读数与该占比自洽。
*
* 闸门（判题同款 gcc9 `-O2 -std=c++17 -static -U_FORTIFY_SOURCE`）：① 编译 exit=0 零错；
* ② n=300000 uniform 5D FNV `274e5232a0c6a620` / sum 2818095627 与 #127003/#127181 逐位相同；
* ③ n∈{1,2,3,4,5,17,100,1000,4000,50000} × {uniform, 三值} 共 20 档对独立 O(n²·5) 暴力
* oracle **bad=0**；④ 件内无 main ✓。
* ================================
*/
/* ===== REFERENCES =====
* duck.ac 用户 saffah_codex_6a_agg3，提交 #126825 <https://duck.ac/submission/126825>
* 用途：**本件正文 = 该提交正文**（该题排行榜 rank-1，t = 265.650235 ms）。逐字沿用其
* 唯一的增量 = 叶内基数排序计数数组的两趟**前缀和合流 + 16 路 u16 SIMD 前缀**
* （`radix_prefix16` / `radix_prefix_pair`，替换原本两段串行的标量 store→load 链），
* 以及它自身继承的全部引用链（其正文下方原有的中文注释头与 `References` 块一字未改，
* 随正文保留）。许可证：duck.ac 提交页未附任何软件许可声明（公共可见代码），故按站方
* 惯例沿用并在此具名引用；本件不主张该部分的著作权。
* duck.ac 用户 saffah_cc_v41_agg1（本账号）提交 #126784 <https://duck.ac/submission/126784>
* （`r2g_` 席）：#126825 的底盘即本账号 #126784 + 上述 1 处改动（本席以非注释行 diff
* 证实：两件正文仅 2 个 hunk，其余 841 行逐字相同）⇒ 本件同时沿用 #126784 的全部内容
* （4-查询合并叶、r2g_ 的 u64 对打包、u8 直方图 sort0、双归并链 KG_MSTEP 内联汇编等），
* 其引用链（本账号 #123464 / #123339 / #120925 / #120558 …、对手 #120558 …）随正文
* 注释块原样保留。
* ======================
* ===== 本次思路 =====
* 正式提交（非试验性）。单变量 = **[R2F-KGB] 叶阈值上移**（本席自创，判题机探针定价）：
* `KGB0/KGB1/KGB2` 1536 → 2048（连带叶数组上界 `KG_BASE` 1536 → 2048）。
*
* 机理：对手 #126825 刚把叶内前缀和/清零这一段**每次调用的定常开销**砍掉（−1.75% 探针），
* 叶的边际成本结构随之改变 —— 叶尺寸的**最优点因此上移**：取更大的叶 ⇒ 叶调用次数与
* 递归层数同时减少（bg `KG_MSTEP`/`ctr_solve` 的每层固定成本），而叶内 O(m²) 比较工作量
* 只按 m 线性增长（Σ nq·ni/8 = M·m/32）。旧最优 1536 是在标量前缀和时代扫出来的。
*
* 判题机定价（`tools/probe.py`，同引擎 + 运行期旋钮 g_kgb，4 臂随机旋转 min-of-R，
* 同设置双臂 = 噪声底，n=100000）：
* 1536（原值）1.0000 · **1792 = 0.9736（−2.64%）** · **2048 = 0.9734（−2.66%）**
* 2560 相对 2048 = 1.0002 · 3072 相对 2048 = 1.0022 ⇒ **峰在 2048~2560，取 2048**
* 噪声底 = 0.9999~1.0001（0.01%）；四臂输出 FNV 逐位相同（值中性）✓
* 同轮另测 [R2F-ZERO]（叶计数数组改 256 位清零）= **1.0001 无效应** ⇒ 不采纳。
*
* 闸门（本席独立跑）：判题同款 gcc9 `-O2 -std=c++17 -static -U_FORTIFY_SOURCE` exit=0 零错；
* n=300000 uniform 5D **FNV `82b08bc287e1217a` / sum 2814957600 与在册 #126784 逐位相同**；
* n∈{1,2,3,17,100,1000,4000} × {uniform, 三值} 对独立 O(n²·5) 暴力 bad=0；件内无 main ✓。
* ================================
*/
/* References: saffah_cc_v41_agg1 https://duck.ac/submission/126784:
copied its accepted five-dimensional event engine, partition-fused radix
histograms and bounded upper bucket range. All inherited citations and license
notices, including our saffah_codex_6a_agg3 https://duck.ac/submission/126121,
remain below.
Idea: Replace the two scalar leaf-radix prefix recurrences by exact sixteen-lane
unsigned-short scans. Three intra-half shift/add stages, one lower-half carry
and one previous-vector carry produce the same inclusive bucket offsets. Zip
the independent arrays, then finish each bounded tail without overreading.
Every leaf has at most1536 events, so all partial sums fit sixteen bits.
Purpose: Reduce serialized bucket-prefix work without changing event ordering,
histogram counts, scatter addresses, comparisons or fallback behavior. */
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #126121 <https://duck.ac/submission/126121>
// 用途：**本文件正文 = 该提交的完整源码为基座**（本题现役排行榜 rank-1，判题机 272.455 ms）。
// 本发在它之上只改 `base_case_t` 叶基数的**记账与桶区间**（见「思路」），
// 算法与数据流（事件序列 CDQ / 双链归并 / 融合前缀计数器 / 叶四查询比较环 / LSD 基数重排）逐字未动。
// [2] 本账号 saffah_cc_v41_agg1，提交 #123464 <https://duck.ac/submission/123464> 与 #126767
// <https://duck.ac/submission/126767>
// 用途：继承链标注。#126121 的第 3 个引用块自述取自本号 #123464（两把刀 `[k11c3-SELADV]`
// 护 CF 的选小/双游标前进、`[k11c3-DEADST]` 模板第 4 参 `ST` 解耦死存）；#126767 是本号
// 对 #126121 的**逐字节取回件**（它把口径钉到严支 0.99·T）。本发接在其后。
// [3] duck.ac 用户 saffah_cc_v41_260924，提交 #123434 <https://duck.ac/submission/123434>
// 用途：事件序列 CDQ 引擎原始出处，经 [2] 隔两级继承进入本文件，本号只作来源标注。
// [4] 本工作区 `problems/1011/notes.md`（【k11c3_ 席】/【k11b_ 席】节）与 `BRIEF.md` §2.19.010
// 用途：判题机同进程双臂 rdtscp 探针的仪器规程与兑现率实录；本发的全部定价按该规程执行。
// 许可：以上均为 duck.ac 上的公开提交与本工作区材料，均未附带任何软件许可证声明；
// 本文件按 RULES §3 逐项署名「原作者账号 + 原始提交地址（duck.ac/submission/<id>）+ 所用内容」，
// 遵守方式 = 完整保留原作者账号、提交 URL 与用途说明于文件开头。
// ======================
// ===== 思路 =====
// 本发 = `#126121` 正文 + **一把刀、三处改动**，全部在 `base_case_t` 的**叶 LSD 基数重排**上；
// 比较/计数结果逐位不变（闸门：n=400/3000/30000/300000 四档 CK+SUM 与基座逐位相同，n=400 另过
// O(n²) 朴素 oracle）。
//
// 【判型】判题机同款 gcc 9.3 `-O2 -std=c++17` 建件 + 本机 rdtscp 独占相位计时（n=300000）实测相位账：
// 叶查询 sweep 34.6% · 归并 28.7% · 融合计数器 ctr_solve 15.9% · **叶 LSD 基数重排 9.94%** ·
// 叶 partition 8.02% · 事件构造 1.51% · sort0 0.99%。
// 把 9.94% 的基数重排拆开：**两趟计数 2.54%**（r1_p1count 1.22 + r4_p2count 1.32）+ 两趟散射 5.24% +
// 桶清零/前缀和 1.92%。⇒ 两趟计数是**纯重复扫描**：partition 循环已经逐事件读过 `seq[t]`，
// 手里就有 `k = E_KA(v)` 与 `TY` 位，而 LSD 两趟要的桶号完全由 (k, TY) 决定。
//
// 【刀 `[r2g-RADIX-FUSE]` 三处】
// ① **两趟计数折进 partition**。两趟计数的桶号可解析地用 partition 已有的量表示：
// pass1 桶号 = `((v>>42) & 1023) + 1` = `(((k<<1)|TY) & 1023) + 1`（v>>42 = (KA<<1)|TY）；
// pass2 桶号 = `((v>>52) & 2047) + 1` = `(k>>9) + 1`（v>>52 = KA>>9，KB/ID 都在低位）。
// ★ 关键：**pass2 的直方图与顺序无关**（多重集相同）⇒ 两趟计数都可在 partition 里一次完成，
// 免掉两趟独立的 `for (t < m)` 扫描（各含一次顺序读 seq/tmp + 一次随机自增）。
// 位置：`if (needOut)` 守卫（脊上 `needOut=0` 的叶不做重排，逐字沿用原语义）。
// ② **pass2 的桶区间按 n 收窄**。键 KA 是 dim-(D+1) 坐标，恒 `< n`（run_t 的顶层事件与各级
// `((v<<21)&(0xFFFFF<<43))` 提升出的新 KA 都是坐标）⇒ 桶号 `k>>9 ≤ ((n-1)>>9)`。
// 新增运行期 `NB1 = ((n-1)>>9)+2`（`if (NB1>2048) NB1=2048`），清零与前缀和只跑 `[0,NB1]`；
// n=300000 ⇒ NB1=588，2049 → 588（省 1460 个桶的清零 + 前缀和）。
// **安全性**：scatter 只读 `c1_[k>>9] ≤ c1_[NB1-2]`，其写游标上界 = `c1_[NB1-1] ≤ m`，
// 全部落在被清零/前缀和覆盖的区间内，其余桶**永不读**（不是"未清零就危险"，是"根本读不到"）。
// ③ **两趟前缀和合到一个循环体**（zipped）。①完成后两次计数在散射前都已完备 ⇒ 两条前缀和
// 互相独立 ⇒ 同一个循环体里交替推进，吃满乱序 ILP（原本是两段串行的 store→load 转发链）。
//
// 【判题机定价（`tools/probe.py`，同进程双臂 rdtscp，n=300000，1 rep/臂，本文件的两个臂）】
// A（`#126121` 逐字副本）= 977,736,290 ticks
// B（本发 = A + 上述三处） = 960,642,006 ticks
// ⇒ **Δ = −17,094,284 ticks = −1.748%**；两臂 CK/SUM 逐位相同（13176558850641716771 / 2821571298）。
// 折判题机 TSC（A 臂 977.736 M ticks ↔ 榜面 272.455 ms ⇒ 0.27866 ms/Mtick）⇒ 预期榜面 ≈ **267.7 ms**，
// 严支可达线 = 0.99·272.455033+1µs = **269.731483 ms** ⇒ 预期落在线内 ≈ 2.0 ms。
//
// 【闸门】判题同款 gcc 9.3 `-O2 -std=c++17 -static -U_FORTIFY_SOURCE` 实编 exit=0 零错零警 ✓ ·
// 件内无 `main` ✓ · n=400 过 O(n²·5) 朴素 oracle（bad=0）✓ · 四档 CK+SUM 与基座逐位相同 ✓ ·
// `tools/undefsw.py` 无未定义宏名 ✓ · 源 47 062 字符 < 102 400 ✓ · 零 `/tmp` ✓。
// （本发为正式性能提交：目的是把本题判题机用时从 272.495 ms 压到严支 269.731 ms 以下。）
// ======================
// ARMVAR KG_BASE=2048 MF_CTRPF=32
#define KG_BASE 1600
#define KGB0 1472
#define KGB1 1600
#define KGB2 1600
#define MF_CTRPF 0
#ifndef MF_CTRPF
#define MF_CTRPF 0
#endif
#ifndef MF_CTR5
#define MF_CTR5 1
#endif
#ifndef MF_LEAFLP
#define MF_LEAFLP 32
#endif
#ifndef MF_LEAFG
#define MF_LEAFG 4
#endif
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
// duck.ac 1010 (4D) / 1011 (5D) -- "count for every point the number of points that are
// smaller in *every* dimension" (n = 300000, coords in [0,n), strict comparisons).
//
// Algorithm: D&C over dimensions on an event sequence, keys carried inline.
// * 2n events: one query event + one insert event per point, keyed by dim 0
// (a query event is stored *before* an insert event of the same key, which encodes
// the strict-inequality rule).
// * solve(seq,m,d): seq is ordered by dim d and every event carries the coordinate of
// dim d+1 in bits 32..63 of its u64. Split by position, recurse both halves, then
// merge the two halves by the carried key (a plain unsigned `<` on the packed u64,
// so equal keys keep query-before-insert order); the same pass collects the cross
// events (left inserts + right queries), re-keys them to dim d+2 and recurses one
// dimension deeper. Ordering by dim d with that tie rule makes "insert before query"
// exactly "all dims so far strictly smaller", and equal coordinates are never counted.
// * d == K-1: every insert-before-query pair now qualifies, so one O(m) running counter
// finishes the innermost dimension (no log factor).
// * m <= 32: direct O(m^2) counting plus an insertion sort.
//
// Local timings (i5-6500 @ 3.2GHz, best of 3, n = 300000 uniform random):
// 4D 1.57 s 5D 3.85 s
// (Board #1: 508.596 ms / 745.358 ms -- this is still 3x / 5x behind.)
// k-dimensional dominance counting: out[i] = #{j : all dims of j < dims of i}
//
// Event-sequence CDQ with keys carried inline.
// event (u64) = (key << 32) | (ins << 31) | id ins = 1 for an "insert" event
// key = coordinate of dim (level+1): the field used by this level's merge.
// The packing makes a plain unsigned comparison put, for equal keys, the query
// event first -- exactly the strict-inequality rule (equal coords must not count).
//
// solve(seq, m, d): seq is ordered by dim d (tie rule above), elements carry the
// key for dim d+1; count pairs (insert before query) with dims d+1..K-1 strictly
// smaller; return `out` = the same events merged by the carried key (= dim d+1).
// Cross = A-inserts + B-queries re-keyed to dim d+2, recursed at level d+1.
// d == K-1 : every (insert before query) pair qualifies -> running counter, O(m).
#include <cstring>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC target("avx2,popcnt,bmi2")
#include <cstdint>
typedef unsigned u32;
typedef uint64_t u64;
typedef unsigned char u8;
typedef unsigned short u16;
#ifndef KG_COND
#define KG_COND 0
#endif
#ifndef KG_PF
#define KG_PF 0
#endif
#ifndef KG_BASE
#define KG_BASE 4096
#endif
#ifndef LEAFNB
#define LEAFNB 8192
#endif
#define KG_MAXK 5
#define KG_ARENA (32u << 20)
namespace kg {
// ------------ fused innermost-dimension counter (prefix count over [0,n)) -------------
// AVX2. Levels: u8 leaf /32 values, u16 per 32 and per 256, u32 per 2048/16384/131072.
// Every u16 block count must stay <= 32767 (madd uses signed lanes); the guard (max
// multiplicity of a last-dim coordinate <= 127) bounds 32*127=4064 and 256*127=32512.
#define CT_MAXV 1048580u
// q9: FIVE levels covering 19 bits (n <= 524288). Two u8 levels replace four
// u16/u32 ones: the u8 levels are the only ones whose SIMD reduction costs ~9
// instructions for FIVE index bits (su8_32 = one 256-bit load + sad), against
// ~10-11 for THREE bits at u16/u32. Index widths: 5/5/3/3/3.
#define CT_LIM (1u << 19)
/* r2i_rv ARENA */
#define A_CT0 0u
#define A_CT1 (A_CT0 + (CT_LIM + 64u))
#define A_MT8 (A_CT1 + (CT_LIM / 32u + 64u))
#define A_MT16A (A_MT8 + 33u * 32u)
#define A_MT16B (A_MT16A + 33u * 32u)
#define A_CMT3 (A_MT16B + 33u * 32u)
#define A_CMT2 (A_CMT3 + 32u * 32u * 2u)
#define A_CT2 (A_CMT2 + 16u * 16u * 2u)
#define A_CT3 (A_CT2 + (2u * (CT_LIM / 1024u + 16u)))
#define A_CT2P (A_CT3 + (2u * (CT_LIM / 16384u + 32u)))
#define A_P3 (A_CT2P + (2u * (CT_LIM / 1024u + 32u)))
#define A_GSC (A_P3 + 64u)
#define A_END (A_GSC + 4u * (1u << 14))
static u8 ARENA[A_END] __attribute__((aligned(64)));
#define CT0 ((u8 *)(ARENA + A_CT0))
#define CT1 ((u8 *)(ARENA + A_CT1))
#define MT8 ((u8 (*)[32])(ARENA + A_MT8))
#define MT16A ((u16 (*)[16])(ARENA + A_MT16A))
#define MT16B ((u16 (*)[16])(ARENA + A_MT16B))
#define CMT3 ((u16 (*)[32])(ARENA + A_CMT3))
#define CMT2 ((u16 (*)[16])(ARENA + A_CMT2))
#define CT2 ((u16 *)(ARENA + A_CT2))
#define CT3 ((u16 *)(ARENA + A_CT3))
#define CT2P ((u16 *)(ARENA + A_CT2P))
#define P3 ((u16 *)(ARENA + A_P3))
#define GSC ((u32 *)(ARENA + A_GSC))
static int FUSE_OK;
static void mt_init(void){ for(int c=0;c<=32;c++){ for(int j=0;j<32;j++) MT8[c][j]=(u8)((j<c)?0xFF:0);
for(int j=0;j<16;j++){ MT16A[c][j]=(u16)((j<c)?0xFFFF:0); MT16B[c][j]=(u16)((16+j<c)?0xFFFF:0); } }
for(int k=0;k<32;k++) for(int j=0;j<32;j++) CMT3[k][j]=(u16)((j>k)?1:0);
for(int k=0;k<16;k++) for(int j=0;j<16;j++) CMT2[k][j]=(u16)((j>k)?1:0); }
static inline u32 su8_32(const u8 *p, int cnt) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
__m256i m = _mm256_cmpgt_epi8(_mm256_set1_epi8((char)cnt),
_mm256_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31));
__m256i x = _mm256_and_si256(v, m);
__m256i s = _mm256_sad_epu8(x, _mm256_setzero_si256());
__m128i t = _mm_add_epi64(_mm256_castsi256_si128(s), _mm256_extracti128_si256(s, 1));
return (u32)((unsigned long long)_mm_cvtsi128_si64(t) + (unsigned long long)_mm_extract_epi64(t, 1));
}
static inline u32 su16_8(const u16 *p, int cnt) {
__m128i v = _mm_loadu_si128((const __m128i *)p);
__m128i x = _mm_and_si128(v, _mm_cmpgt_epi16(_mm_set1_epi16((short)cnt),
_mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7)));
__m128i s = _mm_madd_epi16(x, _mm_set1_epi16(1));
s = _mm_add_epi32(s, _mm_shuffle_epi32(s, _MM_SHUFFLE(1, 0, 3, 2)));
s = _mm_add_epi32(s, _mm_shuffle_epi32(s, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(s);
}
static inline u32 su32_8(const u32 *p, int cnt) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
__m256i x = _mm256_and_si256(v, _mm256_cmpgt_epi32(_mm256_set1_epi32(cnt),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7)));
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
static inline void ctr_ins(u32 v) {
CT0[v]++; CT1[v >> 5]++; CT2[v >> 10]++; CT3[v >> 14]++;
}
/* r2i_rv fast insert: CT2 -> CT2P (prefix suffix-update, one 32 B RMW). */
static inline void ctr_ins_fast(u32 v) {
CT0[v]++; CT1[v >> 5]++;
const u32 w = v >> 10;
__m256i t = _mm256_load_si256((const __m256i *)(CT2P + (w & ~15u)));
t = _mm256_add_epi16(t, _mm256_load_si256((const __m256i *)CMT2[w & 15u]));
_mm256_store_si256((__m256i *)(CT2P + (w & ~15u)), t);
}
static inline void ctr_clr(u32 v) { CT0[v]--; }
static int RSZN;
static inline void ctr_reset_hi(void) {
/* q9: extents = the highest index the query path can READ, +1 */
memset(CT1, 0, sizeof(u8) * ((size_t)((unsigned)RSZN >> 5) + 64));
memset(CT2, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 10) + 16));
memset(CT3, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 14) + 32));
memset(CT2P, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 10) + 32));
memset(P3, 0, 64);
}
static inline void ctr_reset_fast(void) {
// Fast inserts touch CT1 and CT2P; CT2/CT3 remain zero after their slow reset.
memset(CT1, 0, sizeof(u8) * ((size_t)((unsigned)RSZN >> 5) + 64));
memset(CT2P, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 10) + 32));
// P3 is always fully materialized before a fast query, never read otherwise.
}
#define QIDX32 _mm256_setr_epi32(0,1,2,3,4,5,6,7)
#define QIDX16 _mm256_setr_epi16(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
#define QIDX16A _mm256_setr_epi16(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
#define QIDX16B _mm256_setr_epi16(16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31)
#define QIDX8 _mm256_setr_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15, \
16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31)
/* v1: ONE __m256i accumulator across all five levels, ONE horizontal reduce.
Lane assignment is irrelevant because the final reduce sums every lane.
The u16 levels load 16 entries; cnt <= 7 masks the upper 8 off (the array is
padded so the extra 16 bytes are always inside the object). */
static inline u32 ctr_qry(u32 v) {
const __m256i one16 = _mm256_set1_epi16(1);
const __m256i one8 = _mm256_set1_epi8(1); (void)one8;
__m256i acc = _mm256_setzero_si256();
{
/* L3: CT3, u16, 32 entries = 64 B window, base identically 0 (v < 2^19).
mask table index (c3<<5) = ((v>>9) & 0x3e0). */
const u32 t3 = (v >> 9) & 0x3e0u;
__m256i mk = _mm256_load_si256((const __m256i *)MT16A[t3 >> 5]);
__m256i x0 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)CT3), mk);
__m256i x1 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT3 + 16)),
_mm256_load_si256((const __m256i *)MT16B[t3 >> 5]));
acc = _mm256_add_epi32(acc, _mm256_madd_epi16(x0, one16));
acc = _mm256_add_epi32(acc, _mm256_madd_epi16(x1, one16));
}
{
const u32 w2 = v >> 10;
__m256i m = _mm256_load_si256((const __m256i *)MT16A[w2 & 15u]);
__m256i x = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT2 + (w2 & ~15u))), m);
acc = _mm256_add_epi32(acc, _mm256_madd_epi16(x, one16));
}
{
/* r2h_rv: the two u8 levels go through vpsadbw. sad is 5 cyc from the AND
(vs maddubs 5 + madd 5 = 10), and it zero-extends to u64 lanes, so the
two u64 halves of each lane are zero and a single vpaddd folds them into
the u32 accumulator exactly. */
const u32 w1 = v >> 5;
__m256i m = _mm256_load_si256((const __m256i *)MT8[w1 & 31u]);
__m256i x1 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT1 + (w1 & ~31u))), m);
__m256i m0 = _mm256_load_si256((const __m256i *)MT8[v & 31u]);
__m256i x0 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u))), m0);
__m256i s = _mm256_sad_epu8(_mm256_add_epi8(x1, x0), _mm256_setzero_si256());
acc = _mm256_add_epi32(acc, s);
}
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(acc), _mm256_extracti128_si256(acc, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
/* r2i_rv fast query. c2 = CT2P[a>>10], c3 = P3[a>>14] are read by the caller (the
P3 read is preceded by materialising the two P3 registers to memory). The two u8
levels stay a masked 32 B window summed with vpsadbw; the two prefix levels enter as
ONE zero-extended scalar. Every contribution therefore sits in a u64 lane whose
upper 32 bits are 0, so the whole horizontal reduce is vpaddq + movq/pextrq/add. */
static inline u32 ctr_qry_fast(u32 v, u32 c2, u32 c3) {
const u32 w1 = v >> 5;
__m256i m = _mm256_load_si256((const __m256i *)MT8[w1 & 31u]);
__m256i x1 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT1 + (w1 & ~31u))), m);
__m256i m0 = _mm256_load_si256((const __m256i *)MT8[v & 31u]);
__m256i x0 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u))), m0);
__m256i acc = _mm256_sad_epu8(_mm256_add_epi8(x1, x0), _mm256_setzero_si256());
__m128i t = _mm_add_epi64(_mm256_castsi256_si128(acc), _mm256_extracti128_si256(acc, 1));
unsigned long long s64 = (unsigned long long)_mm_cvtsi128_si64(t)
+ (unsigned long long)_mm_extract_epi64(t, 1);
return (u32)((u32)s64 + c2 + c3);
}
static int K;
#ifndef LFAUTO
#define LFAUTO 0
#endif
static int G_LIM;
static int LFSH;
static int LFD;
static int NB1;
static const u32 *XD[KG_MAXK];
static u32 *OUT;
static u64 *AB, *ABtop;
static u64 CROSSB[65536 + 8]; // [R6A11-CROSSB] chain-2 landing buffer
// ===== v2 PACKED EVENT =====================================================
// v = [KA:20][TY:1][KB:20][ID:20] (bits 43..62, 42, 22..41, 0..19)
// KA = the compare key for THIS level (dim D+1), KB = the key for the NEXT
// level (dim D+2), carried so the cross needs no re-key gather.
// A plain u64 compare gives (KA, TY, KB, ID) -- KA primary, TY breaks equal
// KA with query(TY=0) before insert(TY=1) = the strict-inequality rule.
// Requires id < 2^20 and every coordinate < 2^20 (true: n = 1e6 < 2^20).
#define E_ID(v) ((u32)((v) & 0xFFFFFu))
#define E_KB(v) ((u32)(((v) >> 22) & 0xFFFFFu))
#define E_TY(v) ((u32)(((v) >> 42) & 1u))
#define E_KA(v) ((u32)((v) >> 43))
static inline u64 mk2(u32 ka, u32 ty, u32 kb, u32 id) {
return ((u64)ka << 43) | ((u64)ty << 42) | ((u64)kb << 22) | (u64)id;
}
static inline int eid(u64 v) { return (int)(v & 0x7FFFFFFFu); }
static inline int lfbucket(u64 v) { return (int)(((u32)(v >> 32)) >> LFSH); }
static inline int eq_(u64 v) { return (int)((v >> 31) & 1); }
static inline int keep_test(int takeA, int role) { return (takeA ^ role) ^ 1; }
static inline u64 mk(u32 key, int q, int id) {
return ((u64)key << 32) | ((u64)(1 - q) << 31) | (u64)id;
}
// resolve a fusion-level cross: events carry the last dim (K-1) coordinate in bits 32..63
// and are ordered by dim K-2; a left-insert before a right-query with a strictly smaller
// last-dim value is a hit.
/* ---------------- r2i_rv: the original (fallback) engine, byte-for-byte ----------- */
static void ctr_solve_slow(u64 *cross, int cx, u32 *dirt) {
u32 *dirtp = dirt;
const u64 *ce = cross + cx;
const u64 *cp = cross;
for (; cp + 4 <= ce; cp += 4) {
{ u64 v=cp[0]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
{ u64 v=cp[1]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
{ u64 v=cp[2]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
{ u64 v=cp[3]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
}
for (; cp < ce; cp++) {
u64 v = *cp;
u32 val = E_KB(v);
if (E_TY(v)) { ctr_ins(val); *dirtp++ = val; }
else OUT[E_ID(v)] += ctr_qry(val);
}
int nd = (int)(dirtp - dirt);
for (int i = 0; i < nd; i++) {
#if MF_CTRPF > 0
if (i + MF_CTRPF < nd) {
u32 vf = dirt[i + MF_CTRPF];
__builtin_prefetch(&CT0[vf], 1, 1);
__builtin_prefetch(&CT1[vf >> 5], 1, 1);
__builtin_prefetch(&CT2[vf >> 8], 1, 1);
}
#endif
ctr_clr(dirt[i]);
}
ctr_reset_hi();
}
/* ---------------- r2i_rv: the prefix-top engine ----------------------------------- */
#define P3_STEP(a) do { const u32 k_ = (a) >> 14; \
p3v0 = _mm256_add_epi16(p3v0, _mm256_load_si256((const __m256i *)CMT3[k_])); \
p3v1 = _mm256_add_epi16(p3v1, _mm256_load_si256((const __m256i *)(CMT3[k_] + 16))); \
} while (0)
#define P3_MAT() do { \
_mm256_store_si256((__m256i *)P3, p3v0); \
_mm256_store_si256((__m256i *)(P3 + 16), p3v1); } while (0)
#define EV_F(v) do { const u64 v_ = (v); const u32 a_ = E_KB(v_); \
if (E_TY(v_)) { ctr_ins_fast(a_); P3_STEP(a_); *dirtp++ = a_; } \
else { P3_MAT(); OUT[E_ID(v_)] += ctr_qry_fast(a_, CT2P[a_ >> 10], P3[a_ >> 14]); } } while (0)
static void ctr_solve_fast(u64 *cross, int cx, u32 *dirt) {
u32 *dirtp = dirt;
const u64 *ce = cross + cx;
const u64 *cp = cross;
__m256i p3v0 = _mm256_setzero_si256(), p3v1 = _mm256_setzero_si256();
for (; cp + 4 <= ce; cp += 4) {
EV_F(cp[0]); EV_F(cp[1]); EV_F(cp[2]); EV_F(cp[3]);
}
for (; cp < ce; cp++) EV_F(*cp);
int nd = (int)(dirtp - dirt);
for (int i = 0; i < nd; i++) ctr_clr(dirt[i]);
ctr_reset_fast();
(void)p3v0; (void)p3v1;
}
/* P3's entries are u16 prefixes of the WHOLE call's insert count, so the fast path is
used exactly when no call can push them past 65535: inserts <= cx. */
static void ctr_solve(u64 *cross, int cx, u32 *dirt) {
if (cx <= 65535) ctr_solve_fast(cross, cx, dirt);
else ctr_solve_slow(cross, cx, dirt);
}
// [R6A11-CROSSB] same as ctr_solve_fast, but the merged sequence arrives as two contiguous
// pieces (chain-1 prefix, then chain-2) so the in-place compaction memmove is unnecessary.
static void ctr_solve2(u64 *ca, int nca, u64 *cb, int ncb, u32 *dirt) {
u32 *dirtp = dirt;
__m256i p3v0 = _mm256_setzero_si256(), p3v1 = _mm256_setzero_si256();
{ const u64 *ce = ca + nca, *cp = ca;
for (; cp + 4 <= ce; cp += 4) { EV_F(cp[0]); EV_F(cp[1]); EV_F(cp[2]); EV_F(cp[3]); }
for (; cp < ce; cp++) EV_F(*cp); }
{ const u64 *ce = cb + ncb, *cp = cb;
for (; cp + 4 <= ce; cp += 4) { EV_F(cp[0]); EV_F(cp[1]); EV_F(cp[2]); EV_F(cp[3]); }
for (; cp < ce; cp++) EV_F(*cp); }
int nd = (int)(dirtp - dirt);
for (int i = 0; i < nd; i++) ctr_clr(dirt[i]);
ctr_reset_fast();
(void)p3v0; (void)p3v1;
}
// small base: compact insert list + SIMD comparison, then insertion sort by the key
static const u32 VMASK[8][8] = {
{0,0,0,0,0,0,0,0},
{~0u,0,0,0,0,0,0,0},
{~0u,~0u,0,0,0,0,0,0},
{~0u,~0u,~0u,0,0,0,0,0},
{~0u,~0u,~0u,~0u,0,0,0,0},
{~0u,~0u,~0u,~0u,~0u,0,0,0},
{~0u,~0u,~0u,~0u,~0u,~0u,0,0},
{~0u,~0u,~0u,~0u,~0u,~0u,~0u,0}};
static inline u32 hsum256(__m256i v) {
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(v), _mm256_extracti128_si256(v, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
#define L4Q_TAIL(acc, kv, v0, v1, QNx, j, qid) do { \
int jj = (j); __m256i a_ = (acc); \
const int nx_ = (int)(QNx); \
for (; jj + 8 <= nx_; jj += 8) { \
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + jj))); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + jj)))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + jj)))); \
a_ = _mm256_sub_epi32(a_, msk); \
} \
if (jj < nx_) { \
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + jj))); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + jj)))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + jj)))); \
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nx_ - jj])); \
a_ = _mm256_sub_epi32(a_, msk); \
} \
(acc) = a_; \
} while (0)
static inline void radix_prefix16(u16 *p, __m256i &carry) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
v = _mm256_add_epi16(v, _mm256_slli_si256(v, 2));
v = _mm256_add_epi16(v, _mm256_slli_si256(v, 4));
v = _mm256_add_epi16(v, _mm256_slli_si256(v, 8));
const __m256i last = _mm256_set1_epi16(0x0f0e);
__m256i lower = _mm256_permute2x128_si256(v, v, 0x08);
lower = _mm256_shuffle_epi8(lower, last);
v = _mm256_add_epi16(_mm256_add_epi16(v, lower), carry);
_mm256_storeu_si256((__m256i *)p, v);
carry = _mm256_shuffle_epi8(_mm256_permute2x128_si256(v, v, 0x11), last);
}
static inline void radix_prefix_pair(u16 *a, int na, u16 *b, int nb) {
__m256i ca = _mm256_setzero_si256(), cb = ca;
int i = 1, both = na < nb ? na : nb;
for (; i + 15 <= both; i += 16) {
radix_prefix16(a+i, ca); radix_prefix16(b+i, cb);
}
int ia=i, ib=i;
for (; ia+15<=na; ia+=16) radix_prefix16(a+ia,ca);
for (; ib+15<=nb; ib+=16) radix_prefix16(b+ib,cb);
for (; ia<=na; ++ia) a[ia] = (u16)(a[ia]+a[ia-1]);
for (; ib<=nb; ++ib) b[ib] = (u16)(b[ib]+b[ib-1]);
}
template<int ND>
static void base_case_t(u64 *seq, int m, int d, u64 *out, int needOut) {
const int nd = ND;
__attribute__((aligned(32)))
u32 IK[KG_BASE + 8], IC0[KG_BASE + 8], IC1[KG_BASE + 8], IC2[KG_BASE + 8];
// r2g_rv: the query-side u32 pairs are PACKED into u64 arrays. The partition loop is
// STORE-UOP bound (1 store/cycle), so 4 u32 stores -> 2 u64 stores per event. The layout is
// chosen so the sweep's reads are still single 4-byte loads (vpbroadcastd from memory).
u64 QKC[KG_BASE + 8]; // hi32 = k (dim D+1) lo32 = c0 (dim D+2)
u64 QIN[KG_BASE + 8]; // hi32 = ni (insert prefix) lo32 = id
u64 QC12[KG_BASE + 8]; // hi32 = c1 (dim D+3), lo32 = c2 (dim D+4)
const u32 *X1 = nd > 1 ? XD[d + 3] : 0, *X2 = nd > 2 ? XD[d + 4] : 0;
u16 c0_[1024 + 1], c1_[2048 + 1];
if (needOut) { for (int b = 0; b <= 1024; b++) c0_[b] = 0; for (int b = 0; b <= NB1; b++) c1_[b] = 0; }
int ni = 0, nq = 0;
for (int t = 0; t < m; t++) {
#if MF_LEAFLP > 0
if (t + MF_LEAFLP < m) {
int idf = (int)E_ID(seq[t + MF_LEAFLP]);
if (nd > 1) __builtin_prefetch(&X1[idf], 0, 3);
if (nd > 2) __builtin_prefetch(&X2[idf], 0, 3);
}
#endif
u64 v = seq[t];
int id = (int)E_ID(v);
u32 k = E_KA(v);
u32 c0 = 0, c1 = 0, c2 = 0;
if (nd > 0) c0 = E_KB(v); // dim d+2 is IN the event -- no gather
if (nd > 1) c1 = X1[id];
if (nd > 2) c2 = X2[id];
int isins = (int)E_TY(v);
IK[ni] = k; IC0[ni] = c0; IC1[ni] = c1; IC2[ni] = c2; ni += isins;
QKC[nq] = ((u64)k << 32) | (u64)c0; QC12[nq] = ((u64)c1 << 32) | (u64)c2;
if (needOut) { c0_[(((k << 1) | (u32)isins) & 1023u) + 1]++; c1_[(k >> 9) + 1]++; }
QIN[nq] = ((u64)(u32)ni << 32) | (u64)(u32)id; nq += 1 - isins;
}
_mm256_storeu_si256((__m256i*)(QIN+nq),_mm256_setzero_si256());
_mm256_storeu_si256((__m256i*)(QIN+nq+4),_mm256_setzero_si256());
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
int q = 0;
#if MF_LEAFG >= 4
// ================= FOUR-QUERY MERGED LEAF (lane m188) =================
// Gated on `nd <= 1` BY REGISTER ARITHMETIC, not by taste: at nd=2 the live set is
// 3 + 4*4 = 19 ymm against 16. See rig/mk4q.py's header for the full table.
if (nd == 1) {
for (; q + 4 <= nq; q += 4) {
const int nA = (int)(u32)(QIN[q] >> 32), nB = (int)(u32)(QIN[q + 1] >> 32),
nC = (int)(u32)(QIN[q + 2] >> 32), nD = (int)(u32)(QIN[q + 3] >> 32);
__m256i kvA = _mm256_set1_epi32((int)(QKC[q] >> 32)), kvB = _mm256_set1_epi32((int)(QKC[q + 1] >> 32)),
kvC = _mm256_set1_epi32((int)(QKC[q + 2] >> 32)), kvD = _mm256_set1_epi32((int)(QKC[q + 3] >> 32));
__m256i v0A = _mm256_set1_epi32((int)QKC[q]), v0B = _mm256_set1_epi32((int)QKC[q + 1]),
v0C = _mm256_set1_epi32((int)QKC[q + 2]),v0D = _mm256_set1_epi32((int)QKC[q + 3]);
#if (1) /* nd==2 never reaches here, so v1* is only emitted when nd can be 1 -- see below */
__m256i v1A = _mm256_set1_epi32((int)(QC12[q] >> 32)), v1B = _mm256_set1_epi32((int)(QC12[q + 1] >> 32)),
v1C = _mm256_set1_epi32((int)(QC12[q + 2] >> 32)),v1D = _mm256_set1_epi32((int)(QC12[q + 3] >> 32));
#endif
__builtin_prefetch(&OUT[(u32)QIN[q + 8]], 1, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 1, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 10]], 1, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 11]], 1, 1);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(),
accC = _mm256_setzero_si256(), accD = _mm256_setzero_si256();
int j = 0, jA, jB, jC, jD;
for (; j + 8 <= nA; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mA = _mm256_cmpgt_epi32(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v0A, v0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
jA = j;
for (; j + 8 <= nB; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
jB = j;
for (; j + 8 <= nC; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
jC = j;
for (; j + 8 <= nD; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accD = _mm256_sub_epi32(accD, mD);
}
jD = j;
if(nD-jA<8){
const __m256i ik=_mm256_load_si256((const __m256i*)(IK+jA));
__m256i mA=_mm256_cmpgt_epi32(kvA,ik);
__m256i mB=_mm256_cmpgt_epi32(kvB,ik);
__m256i mC=_mm256_cmpgt_epi32(kvC,ik);
__m256i mD=_mm256_cmpgt_epi32(kvD,ik);
const __m256i p0=_mm256_load_si256((const __m256i*)(IC0+jA));
mA=_mm256_and_si256(mA,_mm256_cmpgt_epi32(v0A,p0));
mB=_mm256_and_si256(mB,_mm256_cmpgt_epi32(v0B,p0));
mC=_mm256_and_si256(mC,_mm256_cmpgt_epi32(v0C,p0));
mD=_mm256_and_si256(mD,_mm256_cmpgt_epi32(v0D,p0));
accA=_mm256_sub_epi32(accA,_mm256_and_si256(mA,_mm256_loadu_si256((const __m256i*)VMASK[nA-jA])));
accB=_mm256_sub_epi32(accB,_mm256_and_si256(mB,_mm256_loadu_si256((const __m256i*)VMASK[nB-jA])));
accC=_mm256_sub_epi32(accC,_mm256_and_si256(mC,_mm256_loadu_si256((const __m256i*)VMASK[nC-jA])));
accD=_mm256_sub_epi32(accD,_mm256_and_si256(mD,_mm256_loadu_si256((const __m256i*)VMASK[nD-jA])));
}else{
L4Q_TAIL(accA, kvA, v0A, v1A, nA, jA, (u32)QIN[q]);
L4Q_TAIL(accB, kvB, v0B, v1B, nB, jB, (u32)QIN[q + 1]);
L4Q_TAIL(accC, kvC, v0C, v1C, nC, jC, (u32)QIN[q + 2]);
L4Q_TAIL(accD, kvD, v0D, v1D, nD, jD, (u32)QIN[q + 3]);
}
const __m256i h0=_mm256_hadd_epi32(accA,accB);
const __m256i h1=_mm256_hadd_epi32(accC,accD);
const __m256i h=_mm256_hadd_epi32(h0,h1);
const __m128i sums=_mm_add_epi32(_mm256_castsi256_si128(h),_mm256_extracti128_si256(h,1));
OUT[(u32)QIN[q+0]]+=(u32)_mm_cvtsi128_si32(sums);
OUT[(u32)QIN[q+1]]+=(u32)_mm_extract_epi32(sums,1);
OUT[(u32)QIN[q+2]]+=(u32)_mm_extract_epi32(sums,2);
OUT[(u32)QIN[q+3]]+=(u32)_mm_extract_epi32(sums,3);
}
}
#endif
#if MF_LEAFG >= 4
if (nd == 2) { // q5b/tri: three queries share one insert-block load (nd==3 needs a 4th dim: IC2)
for (; q + 3 <= nq; q += 3) {
const int nA = (int)(u32)(QIN[q] >> 32), nB = (int)(u32)(QIN[q + 1] >> 32), nC = (int)(u32)(QIN[q + 2] >> 32);
__m256i kvA = _mm256_set1_epi32((int)(QKC[q] >> 32)), kvB = _mm256_set1_epi32((int)(QKC[q + 1] >> 32)),
kvC = _mm256_set1_epi32((int)(QKC[q + 2] >> 32));
__m256i v0A = _mm256_set1_epi32((int)QKC[q]), v0B = _mm256_set1_epi32((int)QKC[q + 1]),
v0C = _mm256_set1_epi32((int)QKC[q + 2]);
__m256i v1A = _mm256_set1_epi32((int)(QC12[q] >> 32)), v1B = _mm256_set1_epi32((int)(QC12[q + 1] >> 32)),
v1C = _mm256_set1_epi32((int)(QC12[q + 2] >> 32));
__builtin_prefetch(&OUT[(u32)QIN[q + 8]], 1, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 1, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 10]], 1, 1);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(),
accC = _mm256_setzero_si256();
int j = 0, jA, jB, jC;
for (; j + 8 <= nA; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mA = _mm256_cmpgt_epi32(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v0A, v0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
__m256i v1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, v1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, v1v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v1C, v1v));
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
}
jA = j;
for (; j + 8 <= nB; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
__m256i v1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, v1v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v1C, v1v));
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
}
jB = j;
for (; j + 8 <= nC; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
__m256i v1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v1C, v1v));
accC = _mm256_sub_epi32(accC, mC);
}
jC = j;
if(nC-jA<8){
const __m256i ik=_mm256_load_si256((const __m256i*)(IK+jA));
__m256i mA=_mm256_cmpgt_epi32(kvA,ik);
__m256i mB=_mm256_cmpgt_epi32(kvB,ik);
__m256i mC=_mm256_cmpgt_epi32(kvC,ik);
const __m256i p0=_mm256_load_si256((const __m256i*)(IC0+jA));
mA=_mm256_and_si256(mA,_mm256_cmpgt_epi32(v0A,p0));
mB=_mm256_and_si256(mB,_mm256_cmpgt_epi32(v0B,p0));
mC=_mm256_and_si256(mC,_mm256_cmpgt_epi32(v0C,p0));
const __m256i p1=_mm256_load_si256((const __m256i*)(IC1+jA));
mA=_mm256_and_si256(mA,_mm256_cmpgt_epi32(v1A,p1));
mB=_mm256_and_si256(mB,_mm256_cmpgt_epi32(v1B,p1));
mC=_mm256_and_si256(mC,_mm256_cmpgt_epi32(v1C,p1));
accA=_mm256_sub_epi32(accA,_mm256_and_si256(mA,_mm256_loadu_si256((const __m256i*)VMASK[nA-jA])));
accB=_mm256_sub_epi32(accB,_mm256_and_si256(mB,_mm256_loadu_si256((const __m256i*)VMASK[nB-jA])));
accC=_mm256_sub_epi32(accC,_mm256_and_si256(mC,_mm256_loadu_si256((const __m256i*)VMASK[nC-jA])));
}else{
L4Q_TAIL(accA, kvA, v0A, v1A, nA, jA, (u32)QIN[q]);
L4Q_TAIL(accB, kvB, v0B, v1B, nB, jB, (u32)QIN[q + 1]);
L4Q_TAIL(accC, kvC, v0C, v1C, nC, jC, (u32)QIN[q + 2]);
}
const __m256i h0=_mm256_hadd_epi32(accA,accB);
const __m256i h1=_mm256_hadd_epi32(accC,accC);
const __m256i h=_mm256_hadd_epi32(h0,h1);
const __m128i sums=_mm_add_epi32(_mm256_castsi256_si128(h),_mm256_extracti128_si256(h,1));
OUT[(u32)QIN[q+0]]+=(u32)_mm_cvtsi128_si32(sums);
OUT[(u32)QIN[q+1]]+=(u32)_mm_extract_epi32(sums,1);
OUT[(u32)QIN[q+2]]+=(u32)_mm_extract_epi32(sums,2);
}
}
#endif
for (; q + 2 <= nq; q += 2) {
const int n = ((int)(u32)(QIN[q] >> 32) < (int)(u32)(QIN[q + 1] >> 32)) ? (int)(u32)(QIN[q] >> 32) : (int)(u32)(QIN[q + 1] >> 32);
const int nA = (int)(u32)(QIN[q] >> 32), nB = (int)(u32)(QIN[q + 1] >> 32);
__m256i kvA = _mm256_set1_epi32((int)(QKC[q] >> 32)), kvB = _mm256_set1_epi32((int)(QKC[q + 1] >> 32));
__m256i v0A = _mm256_set1_epi32((int)QKC[q]), v0B = _mm256_set1_epi32((int)QKC[q + 1]);
__m256i v1A = _mm256_set1_epi32((int)(QC12[q] >> 32)), v1B = _mm256_set1_epi32((int)(QC12[q + 1] >> 32));
__m256i v2A = _mm256_set1_epi32((int)QC12[q]), v2B = _mm256_set1_epi32((int)QC12[q + 1]);
__builtin_prefetch(&OUT[(u32)QIN[q + 8]], 1, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 1, 1);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256();
int j = 0;
if (nd >= 3) {
for (; j + 8 <= n; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i C0v = _mm256_load_si256((const __m256i *)(IC0 + j));
__m256i C1v = _mm256_load_si256((const __m256i *)(IC1 + j));
__m256i C2v = _mm256_load_si256((const __m256i *)(IC2 + j));
__m256i mA = _mm256_and_si256(_mm256_cmpgt_epi32(kvA, IKv), _mm256_cmpgt_epi32(v0A, C0v));
__m256i mB = _mm256_and_si256(_mm256_cmpgt_epi32(kvB, IKv), _mm256_cmpgt_epi32(v0B, C0v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, C1v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v2A, C2v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v2B, C2v));
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
}
} else {
for (; j + 8 <= n; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mA = _mm256_cmpgt_epi32(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
if (nd > 0) {
__m256i C0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v0A, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, C0v));
}
if (nd > 1) {
__m256i C1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, C1v));
}
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
}
}
// per-query tails: each member finishes [j, QN[.]) on its own
{
int jj = j;
__m256i acc = accA;
for (; jj + 8 <= nA; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nA) {
int n = nA;
__m256i msk = _mm256_cmpgt_epi32(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[(u32)QIN[q]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = accB;
for (; jj + 8 <= nB; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nB) {
int n = nB;
__m256i msk = _mm256_cmpgt_epi32(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[(u32)QIN[q + 1]] += hsum256(acc);
}
}
#endif
for (; q < nq; q++) {
const int n = (int)(u32)(QIN[q] >> 32);
__m256i kv = _mm256_set1_epi32((int)(QKC[q] >> 32));
__m256i v0 = _mm256_set1_epi32((int)QKC[q]);
__m256i v1 = _mm256_set1_epi32((int)(QC12[q] >> 32));
__m256i v2 = _mm256_set1_epi32((int)QC12[q]);
__m256i acc = _mm256_setzero_si256();
int j = 0;
switch (nd) {
case 0:
for (; j + 8 <= n; j += 8)
acc = _mm256_sub_epi32(acc, _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j))));
break;
case 1:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
case 2:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
default:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC2 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
}
if (j < n) { // partial block, masked by a table lookup
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC2 + j))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - j]));
acc = _mm256_sub_epi32(acc, msk);
}
u32 c = hsum256(acc);
OUT[(u32)QIN[q]] += c;
}
if (!needOut) return;
#if LEAFNB
{
const int NB = (LFAUTO && m <= LEAFNB) ? (m < 8 ? 8 : m) : LEAFNB;
const int LB = (G_LIM < NB) ? G_LIM : NB; // used-bucket bound
(void)LB; (void)NB;
{ // q5b: stable 2-pass LSD radix on the 21-bit (KA, TY) key (bits 42..62)
u64 *tmp = ABtop; ABtop += m;
radix_prefix_pair(c0_, 1024, c1_, NB1);
for (int t = 0; t < m; t++) { u64 v = seq[t]; tmp[c0_[(u32)(v >> 42) & 1023u]++] = v; }
for (int t = 0; t < m; t++) { u64 v = tmp[t]; out[c1_[(u32)(v >> 52) & 2047u]++] = v; }
ABtop -= m;
}
}
#else
for (int t = 1; t < m; t++) {
u64 v = seq[t];
int q = t - 1;
while (q >= 0 && v < seq[q]) { seq[q + 1] = seq[q]; q--; }
seq[q + 1] = v;
}
memcpy(out, seq, sizeof(u64) * (size_t)m);
#endif
}
static void base_case(u64 *seq, int m, int d, u64 *out, int needOut) {
switch (K - (d + 2)) {
case 0: base_case_t<0>(seq, m, d, out, needOut); break;
case 1: base_case_t<1>(seq, m, d, out, needOut); break;
case 2: base_case_t<2>(seq, m, d, out, needOut); break;
default: base_case_t<3>(seq, m, d, out, needOut); break;
}
}
// ---- ND specialisation dispatcher (the fleet's `DN` recipe) ----
static void base_case(u64 *seq, int m, int d, u64 *out, int needOut);
// ================= depth-templated solve: D and KT are compile-time =================
template <bool B, bool NO = true> struct TB {};
template <int V> struct NCL { static const int v = (V < 0) ? 0 : ((V > 3) ? 3 : V); };
template <int KT, int D, bool NO, bool ST = NO> static inline void solve_t(u64 *seq, int m, u64 *out);
// ---- terminal level (D == KT-1): every insert-before-query pair qualifies ----
template <int KT, int D, bool NO, bool ST = NO>
static void solve_b(u64 *seq, int m, u64 *out, TB<true, NO>) {
int cnt = 0;
for (int t = 0; t < m; t++) {
u64 v = seq[t];
if (E_TY(v)) cnt++; // insert
else OUT[E_ID(v)] += (u32)cnt; // query
}
}
// ---- row 1011 lane r2: the merge split into TWO independent merge-path chains --------
// One step of the reference two-pointer merge, verbatim from the bound bytes, so every
// chain emits exactly the reference stream for its own t-range.
/* [k11c3-SELADV] one compare feeds all three consumers (v, pa, pb): CF is read
only by sbb/cmov, none of which writes CF, so gcc9 no longer re-emits the same
`cmpq` three times (it did: setb+rcl/adc, setnb+add, cmovbe each demanded their
own compare because xor/add clobber CF). Same events, same order, same cross
stream; only the codegen changes. */
#define KG_MSTEP(PA, PB, CX, CRS, OIDX) do { \
if (GAT) { __builtin_prefetch(&XN[E_ID((PA)[24])], 0, 3); __builtin_prefetch(&XN[E_ID((PB)[24])], 0, 3); } \
u64 ea_, eb_, tt_, t1_; \
__asm__ ("movq (%[pa]), %[ea]\n\t" \
"movq (%[pb]), %[eb]\n\t" \
"cmpq %[eb], %[ea]\n\t" \
"sbbq %[tt], %[tt]\n\t" \
"cmovb %[ea], %[eb]\n\t" \
"leaq 8(%[pa]), %[t1]\n\t" \
"cmovb %[t1], %[pa]\n\t" \
"leaq 8(%[pb]), %[t1]\n\t" \
"cmovae %[t1], %[pb]\n\t" \
: [ea] "=&r"(ea_), [eb] "=&r"(eb_), [tt] "=&r"(tt_), [t1] "=&r"(t1_), \
[pa] "+r"(PA), [pb] "+r"(PB) \
: : "cc"); \
u64 v = eb_; \
if (ST) out[(OIDX)] = v; \
u64 c = GAT ? ((((v << 21) & (0xFFFFFull << 43)) | (v & (1ull << 42)) | (v & 0xFFFFFull)) | ((u64)XN[E_ID(v)] << 22)) \
: (PROM ? (((v << 21) & (0xFFFFFull << 43)) | (v & (1ull << 42)) | (v & 0xFFFFFull)) : v); \
(CRS)[(CX)] = (FUSE_OK && KT - 3 == D) ? v : c; \
(CX) += (long)(1u - (((u32)(v >> 42) ^ (u32)tt_) & 1u)); \
} while (0)
// (I,J) = the state the merge passes through when it has consumed every element whose
// key (v>>42) is <= k, for the largest such k with I+J <= target. SOUND on inputs that
// are only key-sorted (see the header). I = #A with key <= k, J = #B with key <= k.
static inline void kg_kbound(const u64 *retA, int h, const u64 *retB, int r, int d,
int *pI, int *pJ) {
int lo = (d - r > 0) ? (d - r) : 0;
int hi = (d < h) ? d : h;
while (lo < hi) {
int mid = (lo + hi + 1) >> 1;
if (mid == 0 || (retA[mid - 1] >> 42) < (retB[d - mid] >> 42)) lo = mid; else hi = mid - 1;
}
if (lo == 0) { *pI = 0; *pJ = 0; return; }
u64 kk = retA[lo - 1] >> 42;
if (d > lo) { u64 kb = retB[d - lo - 1] >> 42; if (kb > kk) kk = kb; }
{ int a = lo; while (a < h && (retA[a] >> 42) <= kk) a++; *pI = a; }
{ int a = (d - lo > 0) ? (d - lo) : 0; while (a < r && (retB[a] >> 42) <= kk) a++; *pJ = a; }
}
// ---- general level (D < KT-1) ----
template <int KT, int D, bool NO, bool ST = NO>
static void solve_b(u64 *seq, int m, u64 *out, TB<false, NO>) {
const int KGB = (D == 0) ? KGB0 : ((D == 1) ? KGB1 : KGB2);
if (m <= KGB) {
base_case_t<NCL< KT - (D + 2) >::v>(seq, m, D, out, NO ? 1 : 0);
return;
}
int h = m >> 1, r = m - h;
u64 *save = ABtop;
u64 *retA = ABtop; ABtop += h + 1 + 64;
u64 *retB = ABtop; ABtop += r + 1 + 64;
solve_t<KT, D, true>(seq, h, retA);
solve_t<KT, D, true>(seq + h, r, retB);
retA[h] = ~0ull; retB[r] = ~0ull;
u64 *cross = ABtop; ABtop += m + 1;
// v2: PROMOTE moves the carried KB up into the KA slot (a new level's key);
// GATHER fetches the next KB (dim D+3) -- the ONLY remaining per-level gather,
// and it exists only at D = 0 for KT = 4.
const int PROM = (D + 2 <= KT - 1);
const int GAT = (D + 3 <= KT - 1);
const u32 *XN = XD[(D + 3 <= KG_MAXK - 1) ? (D + 3) : 0];
const u64 *pa_ = retA, *pb_ = retB;
long cx = 0, cx1_out = 0;
int SPLIT2 = 0;
{
// ---- 2 INDEPENDENT MERGE-PATH CHAINS (key-boundary split) -------------
int dd[3], jj[3];
dd[0] = 0; jj[0] = 0;
{ int xi, xj; kg_kbound(retA, h, retB, r, (int)(((long)m * 1) / 2), &xi, &xj);
if (xi + xj < dd[0]) { xi = dd[0] - jj[0]; xj = jj[0]; }
dd[1] = xi + xj; jj[1] = xj; }
dd[2] = m; jj[2] = r;
const u64 *pa1_ = retA + (dd[1] - jj[1]), *pb1_ = retB + jj[1];
long cx1 = 0;
SPLIT2 = (FUSE_OK && KT - 3 == D && m <= 65535);
u64 *cb_ = SPLIT2 ? CROSSB : (cross + dd[1]);
int nmin = dd[1] - dd[0];
if (dd[2] - dd[1] < nmin) nmin = dd[2] - dd[1];
int u = 0;
for (; u < nmin; u++) {
KG_MSTEP(pa_, pb_, cx, cross, u);
KG_MSTEP(pa1_, pb1_, cx1, cb_, dd[1] + u);
}
for (; u < dd[1]; u++) KG_MSTEP(pa_, pb_, cx, cross, u);
for (u = nmin; u < dd[2] - dd[1]; u++) KG_MSTEP(pa1_, pb1_, cx1, cb_, dd[1] + u);
long cp = cx;
if (!SPLIT2 && cx1) memmove(cross + cp, cross + dd[1], sizeof(u64) * (size_t)cx1);
cp += cx1;
cx = cp;
cx1_out = cx1;
}
if (FUSE_OK && D == KT - 3) {
u32 *dirt = (u32 *)ABtop; ABtop += (cx + 3) / 2 + 1;
if (SPLIT2) ctr_solve2(cross, (int)(cx - cx1_out), CROSSB, (int)cx1_out, dirt);
else ctr_solve(cross, cx, dirt);
} else {
solve_t<KT, D + 1, false>(cross, cx, nullptr);
}
ABtop = save;
}
template <int KT, int D, bool NO, bool ST>
static inline void solve_t(u64 *seq, int m, u64 *out) {
if (m <= 1) return;
solve_b<KT, D, NO, ST>(seq, m, out, TB<(D >= KT - 1), NO>());
}
// ====================================================================================
static u32 *g_ord, *g_x;
static unsigned char initial_hist8[1048576];
static void sort0(int n) {
const u32 *X0=XD[0];u32 *cnt=g_x;
memset(initial_hist8,0,(size_t)n+1);
for(int i=0;i<n;i++)++initial_hist8[X0[i]];
u32 sum=0;
for(int i=0;i<=n;i++){u32 c=initial_hist8[i];cnt[i]=sum;sum+=c;}
if(sum!=(u32)n){
memset(cnt,0,sizeof(u32)*(size_t)(n+1));
for(int i=0;i<n;i++)cnt[X0[i]+1]++;
for(int i=0;i<n;i++)cnt[i+1]+=cnt[i];
}
for(int i=0;i<n;i++)g_ord[cnt[X0[i]]++]=(u32)i;
}
template <int KT>
static void run_t(int n, const u32 **x, u32 *out) {
LFSH = 0; LFD = 0;
for (int d = 0; d < KT; d++) XD[d] = x[d];
RSZN = n;
OUT = out;
ABtop = AB;
if (n <= 1) return;
FUSE_OK = 0;
if (KT >= 4 && n > 2 && (u32)n < CT_MAXV && (u32)n <= CT_LIM) {
const u32 *XL = XD[KT - 1];
const unsigned padded=((unsigned)n+31u)&~31u;
memset(CT0,0,padded);
int ok=1;
for(int i=0;i<n;i++)if(++CT0[XL[i]]>127){ok=0;break;}
if(ok){
unsigned group16384=0;
const __m256i zero=_mm256_setzero_si256();
for(unsigned base=0;base<padded;base+=32){
__m256i sums=_mm256_sad_epu8(_mm256_loadu_si256((const __m256i*)(CT0+base)),zero);
__m128i ss=_mm_add_epi64(_mm256_castsi256_si128(sums),_mm256_extracti128_si256(sums,1));
ss=_mm_add_epi64(ss,_mm_srli_si128(ss,8));
unsigned c=(unsigned)_mm_cvtsi128_si32(ss);
if(c>255){ok=0;break;}
group16384+=c;
if((base&16383u)==16352u){if(group16384>65535){ok=0;break;}group16384=0;}
}
if(group16384>65535)ok=0;
}
memset(CT0,0,padded);
FUSE_OK = ok;
}
{
int bit = 0;
while ((n >> bit) >= LEAFNB && bit < 30) bit++;
LFSH = bit; LFD = 0;
}
G_LIM = (int)(((u32)(n - 1)) >> (LFSH + LFD)) + 2;
NB1 = (int)((((u32)(n - 1)) >> 9)) + 2;
if (NB1 > 2048) NB1 = 2048;
sort0(n);
const u32 *X0 = XD[0], *X1 = XD[1];
u64 *ev = AB; // top-level event array: 2n entries
int m=0,g=0;
const u32 *X2=XD[2];
for(int value=0;value<=n;value++){
const int g1=(int)g_x[value],count=g1-g;
if(!count)continue;
for(int t=g;t<g1;t++){
if(t+16<n){u32 idf=g_ord[t+16];__builtin_prefetch(X1+idf,0,1);__builtin_prefetch(X2+idf,0,1);}
u32 q=g_ord[t],a=X1[q],b=X2[q];
ev[m+t-g]=mk2(a,0,b,q);
ev[m+count+t-g]=mk2(a,1,b,q);
}
m+=2*count;g=g1;
}
ABtop = AB + m;
solve_t<KT, 0, false>(ev, m, nullptr);
}
} // namespace kg
static u64 g_arena[KG_ARENA];
static u32 g_ord_st[1048576], g_x_st[1048576];
static int g_ready = 0;
static void setup() {
if (g_ready) return;
g_ready = 1;
kg::mt_init();
kg::AB = g_arena;
kg::g_ord = g_ord_st;
kg::g_x = g_x_st;
}
void count_4d(int n, const unsigned *x[4], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run_t<4>(n, (const unsigned **)x, out);
}
void count_5d(int n, const unsigned *x[5], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run_t<5>(n, (const unsigned **)x, out);
}
#pragma GCC pop_options
#include <cstring>
#include <cstdio>
#include <cstdlib>
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 259.897 ms | 24 MB + 960 KB | Accepted | Score: 100 | 显示更多 |