/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/126046: retained accepted five-dimensional engine including its corrected dense fallback and inherited citations.
* saffah_codex_6a_agg3, https://duck.ac/submission/121747: adapted joint horizontal reduction of multiple leaf-query accumulators.
* Idea: Keep masked tails in their vector accumulators, then reduce each group
* of three or four queries together. Three horizontal additions and one half
* reduction share the work previously repeated for every query separately.
* Purpose: Official experiment reducing SIMD leaf reduction overhead.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/125871: retained compact output arena and exact merge/leaf engine with all inherited citations.
* saffah_codex_6a_agg3, https://duck.ac/submission/125836: reused byte-histogram overflow detection from one-dimensional counting.
* Idea: Use one-byte frequencies for the initial coordinate ordering, reducing
* its random histogram footprint fourfold. The sum of wrapped frequencies can
* equal N only if no byte overflowed; otherwise rebuild the original full-width
* histogram. Prefix starts and stable point placement are unchanged.
* Purpose: Official experiment reducing the initial counting sort's memory traffic.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/123464: retained accepted dual merge chains, flag-preserving selection and dead-store elimination.
* saffah_cc_v41_260924, https://duck.ac/submission/123434: inherited event CDQ engine; all original citations retained below.
* Idea: Remove arena reservations for result streams whose template NO is false.
* The accepted dead-store specialization never writes those streams and its leaf
* returns before sorting when needOut is false. Pass nullptr there instead of
* leaving unused holes between live recursive scratch buffers. Keep every live
* input, child output and cross stream disjoint under the existing stack discipline.
* Also fix the inherited dense-input fallback: preserve the unpromoted event
* only when counter fusion is actually active; recursive fallback needs its
* next coordinate promoted. Also route ND=0 leaves through the exact paired
* loop: the inherited four-query loop unconditionally checked an absent
* second coordinate (all zeros), incorrectly rejecting every full vector. The reference and prior candidate both failed
* seed82731, N300000, five coordinate streams mod1000 (point0:21681 vs18861).
* Purpose: Official experiment reducing arena working-set extent and address work,
* with a general-correctness repair for the inherited non-fused path.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_260924，提交 #123434 <https://duck.ac/submission/123434>
// 用途：**本文件正文 = 该提交的完整源码为基座**（本题当前排行榜 rank-1，判题机 285.196 ms）。
// 本发在它之上**只改两处机制**（见「思路」），算法本体（事件序列 CDQ / 双键内联打包 /
// 双链归并 / 融合前缀计数器 / 叶四查询比较环 / LSD 基数重排）**一律逐字未动**。
// [2] 本账号 saffah_cc_v41_agg1，提交 #123339 <https://duck.ac/submission/123339>
// 用途：本发第一把刀 `[k11c3-SELADV]` 的刀法出处（该发以同一把刀把归并主环从
// `cmp/setb/sbb/and/add` 链改成"一次比较供三处消费者"，判题机 292.669255 → 287.896806 ms）。
// [3] 本工作区 `problems/1011/notes.md` 【k11c2_ 席】节与 `BRIEF.md` §2.19.010
// 用途：[2] 的判题机兑现率实录（探针 −1.681% → 板面 −1.631%）与"同进程双臂探针"
// 的仪器规程（本发的判题机定价即按该规程执行）。
// 许可：以上均为 duck.ac 上的公开提交与本人工作区材料，未见任何附带软件许可证声明；
// 本文件按 RULES §3 逐项署名"原作者账号 + 原始提交地址 + 所用内容"。
// ======================
// ===== 思路 =====
// 本发 = `#123434` 正文 + **两处机制改动**（同轴：归并主环的两笔纯开销），零算法改动。
//
// 【判型】判题机同款 gcc 9.3 `-O2 -std=c++17` 的 `.s`（godbolt g93 实编）显示：`count_5d` 里
// 展开后的 **双链归并主环 = 86 条/2 事件**，其中「选小 + 游标前进」一段是
// `cmp / setb / adc`（IA）+ `cmp / setnb / add`（IB）+ `cmp / cmovbe`（取 v）
// —— **同一次 `cmpq` 被重发三遍**：`xor/add` 改 CF，gcc 只好在第二个、第三个消费者前各补一条。
// 本机 rdtsc 相位账（同件同 n=300000）：叶比较环 36.5% · **归并 30.2%** · 融合计数器 15.1% ·
// 叶基数重排 10.2% · 分割 7.8% ⇒ 归并主环是第二大单块，且其绑定资源是前端条数。
//
// 【刀 1 · `[k11c3-SELADV]` = 一次比较供三处消费者（护 CF）】
// 把「选小 + 两条游标前进」交内联汇编，令 CF 只被 `sbb/cmov` 读、而这三者都不写 CF：
// movq (pa),ea / movq (pb),eb / cmpq eb,ea / sbbq tt,tt / cmovb ea,eb
// / leaq 8(pa),t1 / cmovb t1,pa / leaq 8(pb),t1 / cmovae t1,pb
// `sbbq tt,tt` 给出 tt = −takeA 且 **CF_out ≡ CF_in** ⇒ 三处消费者共用**同一次**比较。
// 逐步等价：`v = min(ea,eb)`；`keep = 1 − (takeA ^ TY(v)) = 1 − ((u32)tt ^ (v>>42)) & 1` ✓。
// 事件次序 / `out` 流 / `cross` 流 / `cx` 逐位不变。
// ⇒ 双链归并主环 **86 → 72 条/2 事件（−16%）**，`cmpq` 6 → 2。
//
// 【刀 2 · `[k11c3-DEADST]` = 归并的 `out` 存储与 `NO` 解耦（删死存）】
// 本引擎 `out` 只有两种读者：① 半区递归 `solve_t<KT,D,true>(seq,h,retA/retB)` —— 返回后
// 被归并当输入读；② 叶 `base_case_t(..., needOut=NO)` 的基数重排输出。而 `NO=false` 的
// "交叉脊"（`run_t` 的 `res`、各级 `cres = ABtop; ABtop += cx;`）**写出去没有任何读者**
// （`run_t` 在 `solve_t<KT,0,false>(ev,m,res)` 之后直接返回）。故新增第 4 个模板参数
// `bool ST = NO`，把 `if (NO) out[(OIDX)] = v;` 改为 `if (ST) out[(OIDX)] = v;`，
// 沿脊的调用点保持实参 `NO=false`（叶的 `needOut` 行为逐字不变）⇒ 删掉 ~1.05M 笔死存、
// 约 7 MB 无效写带宽，语义与数据流**逐位不变**。
//
// 【判题机定价（`tools/probe.py`，同进程双臂 rdtscp，n=300000，1 rep/臂，本文件的两个臂）】
// A（`#123434` 逐字去注释副本）= 1,025,736,488 ticks
// B（本发 = A + 上述两刀） = 991,920,722 ticks
// ⇒ **Δ = −3.297%**（臂序修正 ≈ −3.15%），两臂 CK/SUM 逐位相同（13176558850641716771 / 2821571298）
// ⇒ 预期板面 ≈ 285.196 × 0.967 ≈ **275.8 ms**（阈值 = 0.99·285.195812 + 1 µs = 282.344854）。
//
// 【闸门】本机 n=300 / 1000 / 3000 / 30000 / 300000 五档 CK+SUM 与基座**逐位相同** ✓ ·
// 判题同款 gcc 9.3 `-O2 -std=c++17` 实编 exit=0 零错零警 ✓ · 件内无 `main` ✓
// （本发为**正式性能提交**：目的是把本题判题机用时从 287.897 ms 压到阈值 282.345 ms 以下）。
// ======================
// ARMVAR KG_BASE=2048 MF_CTRPF=32
#define KG_BASE 1536
#define KGB0 1536
#define KGB1 1536
#define KGB2 1536
#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);
}
#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 const u32 *XD[KG_MAXK];
static u32 *OUT;
static u64 *AB, *ABtop;
// ===== 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_hi();
(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);
}
// 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)
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;
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, 1);
if (nd > 2) __builtin_prefetch(&X2[idf], 0, 1);
}
#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;
QIN[nq] = ((u64)(u32)ni << 32) | (u64)(u32)id; nq += 1 - isins;
}
#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]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 10]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 11]], 0, 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;
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]);
__m256i h0=_mm256_hadd_epi32(accA,accB),h1=_mm256_hadd_epi32(accC,accD);
__m256i h=_mm256_hadd_epi32(h0,h1);
__m128i sums=_mm_add_epi32(_mm256_castsi256_si128(h),_mm256_extracti128_si256(h,1));
OUT[(u32)QIN[q]]+=(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]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 10]], 0, 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;
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]);
__m256i h0=_mm256_hadd_epi32(accA,accB),h1=_mm256_hadd_epi32(accC,accC);
__m256i h=_mm256_hadd_epi32(h0,h1);
__m128i sums=_mm_add_epi32(_mm256_castsi256_si128(h),_mm256_extracti128_si256(h,1));
OUT[(u32)QIN[q]]+=(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]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 0, 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;
u16 c0_[1024 + 1], c1_[2048 + 1];
for (int b = 0; b <= 1024; b++) c0_[b] = 0;
for (int t = 0; t < m; t++) c0_[((u32)(seq[t] >> 42) & 1023u) + 1]++;
for (int b = 0; b < 1024; b++) c0_[b + 1] += c0_[b];
for (int t = 0; t < m; t++) { u64 v = seq[t]; tmp[c0_[(u32)(v >> 42) & 1023u]++] = v; }
for (int b = 0; b <= 2048; b++) c1_[b] = 0;
for (int t = 0; t < m; t++) c1_[((u32)(tmp[t] >> 52) & 2047u) + 1]++;
for (int b = 0; b < 2048; b++) c1_[b + 1] += c1_[b];
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, 1); __builtin_prefetch(&XN[E_ID((PB)[24])], 0, 1); } \
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;
{
// ---- 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;
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, cross + dd[1], 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, cross + dd[1], dd[1] + u);
long cp = cx;
if (cx1) memmove(cross + cp, cross + dd[1], sizeof(u64) * (size_t)cx1);
cp += cx1;
cx = cp;
}
if (FUSE_OK && D == KT - 3) {
u32 *dirt = (u32 *)ABtop; ABtop += (cx + 3) / 2 + 1;
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];
memset(CT0, 0, (size_t)n);
int ok = 1;
for (int i = 0; i < n; i++) if (++CT0[XL[i]] > 127) { ok = 0; break; }
memset(CT0, 0, (size_t)n);
/* q9: the two u8 levels must not wrap. The counter is reset per ctr_solve
call and a call's inserts are a SUBSET of all 2n events, so the GLOBAL
bucket counts bound them. CT1 u8 needs <= 255, CT2 u16 needs <= 65535. */
if (ok) {
memset(GSC, 0, sizeof(u32) * ((size_t)(((unsigned)n) >> 5) + 64));
for (int i = 0; i < n && ok; i++) {
u32 b = XL[i] >> 5;
if (GSC[b] >= 255u) ok = 0; else GSC[b]++;
}
memset(GSC, 0, sizeof(u32) * ((size_t)(((unsigned)n) >> 10) + 16));
for (int i = 0; i < n && ok; i++) {
u32 b = XL[i] >> 10;
if (GSC[b] >= 65535u) ok = 0; else GSC[b]++;
}
memset(GSC, 0, sizeof(u32) * ((size_t)(((unsigned)n) >> 14) + 32));
for (int i = 0; i < n && ok; i++) {
u32 b = XL[i] >> 14;
if (GSC[b] >= 65535u) ok = 0; else GSC[b]++;
}
}
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;
sort0(n);
const u32 *X0 = XD[0], *X1 = XD[1];
u64 *ev = AB; // top-level event array: 2n entries
int m = 0, g = 0;
while (g < n) {
int g1 = g;
u32 v = X0[g_ord[g]];
while (g1 < n && X0[g_ord[g1]] == v) g1++;
const u32 *X2 = XD[2];
for (int t = g; t < g1; t++) { int q = (int)g_ord[t]; ev[m++] = mk2(X1[q], 0, X2[q], (u32)q); }
for (int t = g; t < g1; t++) { int q = (int)g_ord[t]; ev[m++] = mk2(X1[q], 1, X2[q], (u32)q); }
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 | 275.067 ms | 24 MB + 920 KB | Accepted | Score: 100 | 显示更多 |