// References: saffah_cc_v41_agg1, https://duck.ac/submission/106300 ; copied its full accepted 1011e6 engine, including suffix counter, thresholds, and inherited references listed below. No license notice appeared with that submission. The original engine itself incorporates ideas from saffah_codex_6s_agg2 submissions cited in its inherited reference block.
// Approach: In depth-two leaves (one remaining coordinate), process four adjacent queries together; they load each insert key and coordinate once per vector block and apply four independent comparisons. Each query then handles its own short tail and output update. Other depths keep the accepted paired-query kernel.
// Purpose: Experimental judged submission to test whether the official n=1000000 workload gains enough from four-query SIMD load sharing; official Try It at n=400000 improved 537.107 ms to 521.569 ms with identical output hash.
// ===== REFERENCES =====
// [1] 本账号（saffah_cc_v41_agg1）当前最好件 = 同 lane `sm7_` 的 `work/sm7_kg3_2048.cpp`
// （提交 #106204 = 2027.205180 ms，= 用户 saffah_codex_6s_agg2 #106156 的逐字节复刻 + `KG3 2560→2048`）
// —— **本文件的正文基座**（本发只在它的 `ctr_ins`/`ctr_qry` 两处做单变量替换）。
// [2] ★ **本刀的来源与本账号的已兑现读数**（同族 `1010` / `1010a` 两题，同一份计数器结构）：
// duck.ac 用户 saffah_codex_6s_agg2 提交 #103094 <https://duck.ac/submission/103094>
// —— "把顶层 CT5 改成**后缀数组**（`CT5[0..7] += CT5INC[v>>17]`，查询退化成一条标量载入）"
// 这一恒等变形的出处；本账号把**同一条恒等变形推广到 CT1..CT4**（每题只改计数器两处）：
// * 本账号 #103203 <https://duck.ac/submission/103203>（题 1010a）—— 判题机 **−3.20%**（181.443 → 175.637）；
// * 本账号 #106210 <https://duck.ac/submission/106210>（题 1010）—— 该件的直接后继成分；
// * 本账号 #106218 <https://duck.ac/submission/106218>（题 1010a）。
// 本文件把同一变形移植到本题（`KGB`/`CT0..CT5` 阶梯与 4D 完全同构）。
// **语义恒等**：`CT_k[p] := 本组内块号 < p 的插入数` ⇒ 查询 = 一条标量载入；
// 插入 = 该组内一组 8 车道的宽加（行 r 的掩码在"严格高于 r"的车道上放 1）。
// 各层贡献互不重叠、合起来恰好覆盖 `[0, v)`，与原"前缀归约"逐位等价。
// [3] 同 lane `sm7_` 的构建/闸门工具：`problems/1010e6/work/ke6_build.py`（`--brute` / `--probe`）。
// [4] 本文件未取用其它他人代码。
// ======================
// ===== 思路 =====
// 现状：`mine = 2027.205180`（#106204）· `T = 2027.488114`（对手 #106156）· 严支 `thr = 2007.214233`
// ⇒ 缺口 **19.990947 ms（0.986%）**，红态（`mine < T` 但严支要求 ≤0.99·T）。
//
// ★ 唯一改动 = **[SFX] 计数器后缀化**（只碰 `ctr_ins` / `ctr_qry` 两个函数 + CT2 的声明）：
// 原来的每一层都是「块内计数 + 查询时做一次掩码前缀归约」：
// `ctr_ins`: `CT1[v>>5]++; CT2[v>>8]++; ...`（6 条标量 RMW）
// `ctr_qry`: `su32_8(CT4+…)+su32_8(CT3+…)+su16_8(CT2+…)+su16_8(CT1+…)+su8_32(CT0+…)`（4 条向量归约 + 1 条）
// 改成**后缀数组**后：查询 = **5 条标量载入**（CT5/CT4/CT3/CT2/CT1）；插入 = 5 条 16/32 B 宽加。
// 唯一约束（与 1010 相同）：组内后缀值上界 = 8 × 本层块跨度 × 单值重数 ⇒
// `CT1`(跨度32)=32512 仍 u16 ✓；**`CT2`(跨度256)=260096 ⇒ 必须 u16→u32**；CT3/CT4/CT5 本就是 u32 ✓。
//
// 闸门：
// 1) `ke6_build.py --brute --kbrute 5`（独立 O(n²k) 暴力，3 种分布）⇒ **mism=0 / BRUTE OK** ✓；
// 2) **判题机 probe（n=4e5 k=5，checksum 恒为 4981929175 = 本 lane 的历史常数）**：
// 基座 **1967.71 Mcyc** → 本发 **1912.46 Mcyc** ⇒ **−2.81%**，**两臂 checksum 逐位相同** ✓。
// 未改动：引擎其余部分、`KGB_(d)`、`LEAFNB 8192`、`MF_*` 全部常量、叶/归并/事件构造。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #106156 <https://duck.ac/submission/106156>（2027.488114 ms，当代 T）。
// 用途：**本文件正文 = 该提交正文逐字节复制**（RULES §2.16 第一步）。它自述其父本 = **本账号 #102291**
// （= 我方最好件 2048.201396 ms；含我方 #100459 的"链式内联 XE"）⇒ 本次复刻 = 追平/反超对手当代件。
// [2] 本账号 #102291 <https://duck.ac/submission/102291>（`work/e11_sub_1011e6.cpp`）—— 对手文件里已含我方 XE 刀 ✓。
// [3] duck.ac 用户 saffah_cc_v41_260924，提交 #102263 <https://duck.ac/submission/102263> —— 更早的 CDQ 引擎父本。
// ======================
// ===== 思路 =====
// 入场：`mine = 2048.201396`（#102291）/ `T = 2027.488114`（#106156）/ 宽支 `1.005*T = 2037.626555` ⇒ 红 10.574841 ms；
// 严支（交新最好后）0.99*T+1µs ⇒ 缺 ~40.99 ms。
// ★ 本发 = **纯复刻**（协调者 §2.18.993 的"按是否结构件排序地抄"）。逐行 diff `work/e11_rival_102263.cpp → ref2/rival_106156.cpp`
// 得到对手相对其*旧*父本的全部增量（其中 `[XE]` 那几段是**我方 #102291 的刀**，他照抄了 ✓），**他自己的两件**是：
// ① **事件构造改成计数排序的组尾 + 两种角色一次载入双写**：`for(v=0..n) { g1=g_x[v]; ... ev[m+t-g]=mk(key,1,id,nk); ev[m+cnt+t-g]=mk(key,0,id,nk); }`
// —— 原来的 `while(g<n){ v=X0[g_ord[g]]; ... }` 每事件要**随机 gather 一次首维** ⇒ 本件把它换成**顺序扫 `g_x`（计数排序组尾）** ✓（对手自述 "avoid random first-coordinate reads"）
// 且把两次 `mk()` 的**公共载入合并** ✓。按本 lane 的兑现律（"省 4 MB 级随机 gather"类，probe 低估 5.7×），这是**结构件** ⇒ 优先 ✓。
// ② **`KG2 2560 → 2048`**（二级叶阈值，参数件 ✓）。
// ★ 为何整份复刻：我方 2048.201 > T 2027.488（慢 1.01%）⇒ 按 §2.18.566② "mine > T ⇒ 整份复刻有效" ⇒ 一次把这 1.01% 压掉 ✓。
// ★ 正确性闸门：① 与 `work/e11_sub_1011e6.cpp`（#102291）**逐位对拍**（多 mode × 多 n）；② gcc93 编译；③ 配置回显 + 断言输出非空。
// 对手件在判题机上 Accepted ✓；其自述的实现坑（`cnt[]` 静态化）已在 [1] 里修好 ✓。
// ================
#ifndef MF_CTRPF
#define MF_CTRPF 32
#endif
#ifndef MF_CTR5
#define MF_CTR5 1
#endif
#ifndef OUTPF
#define OUTPF 4
#endif
#ifndef MF_LEAFLP
#define MF_LEAFLP 64
#endif
#ifndef MF_LEAFG
#define MF_LEAFG 2
#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>
static volatile unsigned long long gCV_ = 0;
#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 2560
#endif
/* ---- DEPTH-DEPENDENT LEAF THRESHOLD (Lane 1011e6, 2026-09-27) ----------------
The leaf's cost is ~ (number of comparisons + 1) * m per event, because it is
an O(m^2/8) all-pairs scan whose inner block does 2*(1+nd) ALU uops per 8 pairs.
nd = K-(d+2) is a COMPILE-TIME constant per depth, so the optimal leaf size m*(d)
is NOT the same at every depth. MEASURED at n=3e5,k=5 (judge, custom_test):
d=0 : 256 leaf calls, 600 000 events, nd=3 (4 comparisons)
d=1 : 1024 leaf calls, 2 400 008 events, nd=2 (3)
d=2 : 1792 leaf calls, 4 199 134 events, nd=1 (2)
=> 58 % of all leaf events are at d=2, where the leaf is cheap per event, and
only 8 % are at d=0 where it is dearest. The single-valued KG_BASE sweep that
the fleet ran could not see this: it moves all three thresholds together, and a
sum of curves whose optima differ looks FLAT (measured: total varies +-0.16 %
over 1280..4096 at n=3e5 while base_case itself varies 4.6 %).
The arrays are sized by the LARGEST threshold. */
#define KG0 1920
#define KG1 2560
#define KG2 2048
#define KG3 2048
static int KGB[8] = { KG0, KG1, KG2, KG3, KG3, KG3, KG3, KG3 };
#define KG_BMAX (KG0>KG1 ? (KG0>KG2? (KG0>KG3?KG0:KG3) : (KG2>KG3?KG2:KG3)) : (KG1>KG2? (KG1>KG3?KG1:KG3) : (KG2>KG3?KG2:KG3)))
#define KGB_(dd) (KGB[(dd) < 0 ? 0 : ((dd) > 7 ? 7 : (dd))])
#ifndef LEAFNB
#define LEAFNB 8192
#endif
#define KG_MAXK 5
#define KG_ARENA (16000000)
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
static u8 CT0[CT_MAXV];
static u16 CT1[CT_MAXV / 32 + 8];
static u32 CT2[CT_MAXV / 256 + 8]; /* [SFX] u16->u32：组内后缀上界 260096 */
static u32 CT3[CT_MAXV / 2048 + 8];
static u32 CT4[CT_MAXV / 16384 + 8];
static u32 CT5[CT_MAXV / 131072 + 8];
static int FUSE_OK;
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);
}
/* [SFX] 后缀计数器：行 r 在"严格高于 r"的车道上放 1 ⇒ 插入落在块 r 时，
把它所在 8 块组里所有需要开始统计它的后缀位上各加一。查询退化成一条标量载入。 */
static const u16 INC16[8][8] __attribute__((aligned(32))) = {
{0,1,1,1,1,1,1,1},{0,0,1,1,1,1,1,1},{0,0,0,1,1,1,1,1},{0,0,0,0,1,1,1,1},
{0,0,0,0,0,1,1,1},{0,0,0,0,0,0,1,1},{0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0}};
static const u32 INC32[8][8] __attribute__((aligned(32))) = {
{0,1,1,1,1,1,1,1},{0,0,1,1,1,1,1,1},{0,0,0,1,1,1,1,1},{0,0,0,0,1,1,1,1},
{0,0,0,0,0,1,1,1},{0,0,0,0,0,0,1,1},{0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0}};
#define SFX16(p, r) do { __m128i _a = _mm_load_si128((const __m128i *)(p)); \
__m128i _b = _mm_load_si128((const __m128i *)(INC16 + (r))); \
_mm_store_si128((__m128i *)(p), _mm_add_epi16(_a, _b)); } while (0)
#define SFX32(p, r) do { __m256i _a = _mm256_load_si256((const __m256i *)(p)); \
__m256i _b = _mm256_load_si256((const __m256i *)(INC32 + (r))); \
_mm256_store_si256((__m256i *)(p), _mm256_add_epi32(_a, _b)); } while (0)
static inline void ctr_ins(u32 v) {
CT0[v]++;
SFX16(CT1 + ((v >> 8) << 3), (v >> 5) & 7);
SFX32(CT2 + ((v >> 11) << 3), (v >> 8) & 7);
SFX32(CT3 + ((v >> 14) << 3), (v >> 11) & 7);
SFX32(CT4 + ((v >> 17) << 3), (v >> 14) & 7);
SFX32(CT5, (v >> 17) & 7); /* CT5 活跃项 <= 8 ⇒ 整组一个 SFX32 */
}
/* ARM `CT` (ported from the sister row 1011a): RESET THE HIERARCHY INSTEAD OF DECREMENTING IT.
ctr_clr mirrored all six of ctr_ins's read-modify-writes, so every element cost SIX random
RMWs to insert and SIX more to clear. Five of those six arrays are SMALL -- CT1 18 KB,
CT2 2.4 KB, CT3 592 B, CT4 76 B, CT5 16 B at n=1e5 (62 KB / 8 KB / 2 KB / 244 B / 60 B at
n=1e6) -- and they are only ever read through ctr_qry, which sums them, so restoring them to
their initial zero is a pure RESET, not an inverse walk. CT0 stays a per-element decrement
because it is the one level that is genuinely large and densely addressed, so walking `dirt`
still beats clearing it. Semantics are unchanged: ctr_qry reads the same six levels and the
structure is left in its initial state for the next call, which is what the decrements were for. */
static inline void ctr_clr(u32 v) { CT0[v]--; }
static inline void ctr_reset_hi(void) {
memset(CT1, 0, sizeof(CT1)); memset(CT2, 0, sizeof(CT2));
memset(CT3, 0, sizeof(CT3)); memset(CT4, 0, sizeof(CT4));
memset(CT5, 0, sizeof(CT5));
}
static inline u32 ctr_qry(u32 v) {
u32 k5 = v >> 17, s = 0;
s += CT5[v >> 17]; /* [SFX] 标量后缀载入替代向量归约 */
s += CT4[v >> 14];
s += CT3[v >> 11];
s += CT2[v >> 8];
s += CT1[v >> 5];
s += su8_32(CT0 + (v & ~31u), (int)(v & 31));
return s;
}
static int K;
#ifndef LFAUTO
#define LFAUTO 0
#endif
static int LFSH;
static int LFD;
static const u32 *XD[KG_MAXK];
static u32 *OUT;
static u64 *AB, *ABtop;
static inline int eid(u64 v) { return (int)(v & 0x1FFFFFu); }
static inline int lfbucket(u64 v) { return (int)(((u32)(v >> 43)) >> LFSH); }
static inline int eq_(u64 v) { return (int)((v >> 42) & 1); }
static inline int keep_test(int takeA, int role) { return 1 - (takeA ^ role); }
static inline u32 embnk(u64 v) { return (u32)((v >> 21) & 0x1FFFFFu); }
static inline u64 mk(u32 key, int q, int id, u32 nk) {
return ((u64)key << 43) | ((u64)(1 - q) << 42) | ((u64)nk << 21) | (u64)id;
}
static void solve(u64 *seq, int m, int d, u64 *out, int needOut);
// 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.
static void ctr_solve(u64 *cross, int cx, u32 *dirt) {
int nd = 0;
for (int i = 0; i < cx; i++) {
#if MF_CTRPF > 0
if (i + MF_CTRPF < cx) {
u32 vf = (u32)(cross[i + MF_CTRPF] >> 43);
__builtin_prefetch(&CT0[vf], 1, 0);
__builtin_prefetch(&CT1[vf >> 5], 1, 0);
__builtin_prefetch(&CT2[vf >> 8], 1, 0);
__builtin_prefetch(&OUT[cross[i + MF_CTRPF] & 0x1FFFFFu], 1, 0);
}
#endif
u64 v = cross[i];
u32 val = (u32)(v >> 43);
if ((v >> 42) & 1) { ctr_ins(val); dirt[nd++] = val; }
else OUT[eid(v)] += ctr_qry(val);
}
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();
}
// 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);
}
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_BMAX + 8], IC0[KG_BMAX + 8], IC1[KG_BMAX + 8], IC2[KG_BMAX + 8];
u32 QK[KG_BMAX + 8], QC0[KG_BMAX + 8], QC1[KG_BMAX + 8], QC2[KG_BMAX + 8];
u32 QID[KG_BMAX + 8], QN[KG_BMAX + 8];
const u32 *X0 = nd > 0 ? XD[d + 2] : 0, *X1 = nd > 1 ? XD[d + 3] : 0, *X2 = nd > 2 ? XD[d + 4] : 0;
/* [XE] the X0 dimension (dim d+2) is already inside the event: at d == 0 by mk()'s nk,
at d >= 1 by the level-(d-1) cross's embed (valid exactly while d + 2 <= K - 1) */
const int useEmb = (nd > 0) && (d == 0 || d <= K - 3);
int ni = 0, nq = 0;
for (int t = 0; t < m; t++) {
#if MF_LEAFLP > 0
if (t + MF_LEAFLP < m) {
int idf = eid(seq[t + MF_LEAFLP]);
if (nd > 0 && !useEmb) __builtin_prefetch(&X0[idf], 0, 1);
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 = eid(v);
u32 k = (u32)(v >> 43);
u32 c0 = 0, c1 = 0, c2 = 0;
if (nd > 0) c0 = useEmb ? embnk(v) : X0[id];
if (nd > 1) c1 = X1[id];
if (nd > 2) c2 = X2[id];
int isins = (int)((v >> 42) & 1u);
IK[ni] = k; IC0[ni] = c0; IC1[ni] = c1; IC2[ni] = c2; ni += isins;
QK[nq] = k; QC0[nq] = c0; QC1[nq] = c1; QC2[nq] = c2; QID[nq] = (u32)id; QN[nq] = (u32)ni; nq += 1 - isins;
}
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
/* ARM L_outpf (lane l1011e6_ctr): the leaf's counting loop ends in a RANDOM
read-modify-write, `OUT[QID[q]] += c`. A random LOAD is absorbed by the OOO window (on
this row the X-gather id-packing is an exact null, +0.014 %), but a random RMW is NOT:
the store must own the line before it can retire. `L_outhot` (the same loop, destination
forced onto two cache lines) prices the whole component at -3.86 % on the judge at
n=3e5. This arm SW-prefetches the destination 2*OUTPF queries ahead -- the mechanism the
construction lane banked at -5.90 % on this row for the merge's RMWs. A prefetch has no
semantics, so the checksum must be unchanged. */
int q = 0;
if (nd == 1) {
for (; q + 4 <= nq; q += 4) {
if (q + 4 * OUTPF < nq) {
__builtin_prefetch(&OUT[QID[q + 4 * OUTPF + 0]], 1, 0);
__builtin_prefetch(&OUT[QID[q + 4 * OUTPF + 1]], 1, 0);
__builtin_prefetch(&OUT[QID[q + 4 * OUTPF + 2]], 1, 0);
__builtin_prefetch(&OUT[QID[q + 4 * OUTPF + 3]], 1, 0);
}
int n = (int)QN[q];
if ((int)QN[q + 1] < n) n = (int)QN[q + 1];
if ((int)QN[q + 2] < n) n = (int)QN[q + 2];
if ((int)QN[q + 3] < n) n = (int)QN[q + 3];
const int n0 = (int)QN[q + 0];
__m256i kv0 = _mm256_set1_epi32((int)QK[q + 0]);
__m256i v0 = _mm256_set1_epi32((int)QC0[q + 0]);
__m256i ac0 = _mm256_setzero_si256();
const int n1 = (int)QN[q + 1];
__m256i kv1 = _mm256_set1_epi32((int)QK[q + 1]);
__m256i v1 = _mm256_set1_epi32((int)QC0[q + 1]);
__m256i ac1 = _mm256_setzero_si256();
const int n2 = (int)QN[q + 2];
__m256i kv2 = _mm256_set1_epi32((int)QK[q + 2]);
__m256i v2 = _mm256_set1_epi32((int)QC0[q + 2]);
__m256i ac2 = _mm256_setzero_si256();
const int n3 = (int)QN[q + 3];
__m256i kv3 = _mm256_set1_epi32((int)QK[q + 3]);
__m256i v3 = _mm256_set1_epi32((int)QC0[q + 3]);
__m256i ac3 = _mm256_setzero_si256();
int j = 0;
for (; j + 8 <= n; j += 8) {
__m256i ik = _mm256_load_si256((const __m256i *)(IK + j));
__m256i ic = _mm256_load_si256((const __m256i *)(IC0 + j));
__m256i m0 = _mm256_and_si256(_mm256_cmpgt_epi32(kv0, ik), _mm256_cmpgt_epi32(v0, ic));
ac0 = _mm256_sub_epi32(ac0, m0);
__m256i m1 = _mm256_and_si256(_mm256_cmpgt_epi32(kv1, ik), _mm256_cmpgt_epi32(v1, ic));
ac1 = _mm256_sub_epi32(ac1, m1);
__m256i m2 = _mm256_and_si256(_mm256_cmpgt_epi32(kv2, ik), _mm256_cmpgt_epi32(v2, ic));
ac2 = _mm256_sub_epi32(ac2, m2);
__m256i m3 = _mm256_and_si256(_mm256_cmpgt_epi32(kv3, ik), _mm256_cmpgt_epi32(v3, ic));
ac3 = _mm256_sub_epi32(ac3, m3);
}
{
int jj = j;
__m256i acc = ac0;
for (; jj + 8 <= n0; jj += 8) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv0, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi32(acc, mask);
}
if (jj < n0) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv0, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + jj))));
mask = _mm256_and_si256(mask, _mm256_loadu_si256((const __m256i *)VMASK[n0 - jj]));
acc = _mm256_sub_epi32(acc, mask);
}
OUT[QID[q + 0]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = ac1;
for (; jj + 8 <= n1; jj += 8) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv1, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi32(acc, mask);
}
if (jj < n1) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv1, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC0 + jj))));
mask = _mm256_and_si256(mask, _mm256_loadu_si256((const __m256i *)VMASK[n1 - jj]));
acc = _mm256_sub_epi32(acc, mask);
}
OUT[QID[q + 1]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = ac2;
for (; jj + 8 <= n2; jj += 8) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv2, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi32(acc, mask);
}
if (jj < n2) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv2, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC0 + jj))));
mask = _mm256_and_si256(mask, _mm256_loadu_si256((const __m256i *)VMASK[n2 - jj]));
acc = _mm256_sub_epi32(acc, mask);
}
OUT[QID[q + 2]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = ac3;
for (; jj + 8 <= n3; jj += 8) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv3, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v3, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi32(acc, mask);
}
if (jj < n3) {
__m256i mask = _mm256_and_si256(_mm256_cmpgt_epi32(kv3, _mm256_load_si256((const __m256i *)(IK + jj))), _mm256_cmpgt_epi32(v3, _mm256_load_si256((const __m256i *)(IC0 + jj))));
mask = _mm256_and_si256(mask, _mm256_loadu_si256((const __m256i *)VMASK[n3 - jj]));
acc = _mm256_sub_epi32(acc, mask);
}
OUT[QID[q + 3]] += hsum256(acc);
}
}
}
for (; q + 2 <= nq; q += 2) {
if (q + 2 * OUTPF < nq) {
__builtin_prefetch(&OUT[QID[q + 2 * OUTPF]], 1, 0);
__builtin_prefetch(&OUT[QID[q + 2 * OUTPF + 1]], 1, 0);
}
const int n = ((int)QN[q] < (int)QN[q + 1]) ? (int)QN[q] : (int)QN[q + 1];
const int nA = (int)QN[q], nB = (int)QN[q + 1];
__m256i kvA = _mm256_set1_epi32((int)QK[q]), kvB = _mm256_set1_epi32((int)QK[q + 1]);
__m256i v0A = _mm256_set1_epi32((int)QC0[q]), v0B = _mm256_set1_epi32((int)QC0[q + 1]);
__m256i v1A = _mm256_set1_epi32((int)QC1[q]), v1B = _mm256_set1_epi32((int)QC1[q + 1]);
__m256i v2A = _mm256_set1_epi32((int)QC2[q]), v2B = _mm256_set1_epi32((int)QC2[q + 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[QID[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[QID[q + 1]] += hsum256(acc);
}
}
#endif
for (; q < nq; q++) {
if (q + OUTPF < nq) __builtin_prefetch(&OUT[QID[q + OUTPF]], 1, 0);
const int n = (int)QN[q];
__m256i kv = _mm256_set1_epi32((int)QK[q]);
__m256i v0 = _mm256_set1_epi32((int)QC0[q]);
__m256i v1 = _mm256_set1_epi32((int)QC1[q]);
__m256i v2 = _mm256_set1_epi32((int)QC2[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[QID[q]] += c;
}
if (!needOut) return;
#if LEAFNB
{
const int NB = (LFAUTO && m <= LEAFNB) ? (m < 8 ? 8 : m) : LEAFNB;
static u32 cnt[LEAFNB + 2]; /* FIX: was a stack local relying on the previous call's bitmap clear; first-ever use of a bucket at a stack slot read uninitialised bytes (verified: segfault at n=5e5 in the rig). static = .bss = zero-initialised, so the "cnt[] is all-zero at entry" invariant is true by induction. Zero cost: the clear-via-bitmap is unchanged. */
/* ---- SPARSE REBUILD ----------------------------------------------------------
The shipped counting sort walks all NB = 8192 buckets twice (an NB-iteration
prefix and an NB-iteration bucket walk) for a leaf that holds at most
KGB_(d) <= 3840 events, so ~60% of those iterations look at EMPTY buckets.
A 1 KB occupancy bitmap, set during the count pass, lets the prefix, the walk
and the clear each visit only the buckets that actually hold something:
O(NB) -> O(#non-empty + NB/64).
EQUIVALENCE IS A PROOF: cnt[] is only ever read for a bucket we set a bit for,
which is exactly a bucket the scatter addresses; and `lo` for a non-empty
bucket b' is the end of the previous NON-EMPTY bucket, which equals the end
of b'-1 because every bucket in between is empty. Same ranges, same order,
same insertion sorts, bit-identical out[].
cnt[] is left zeroed via the same bitmap walk, so the next call starts clean
(base_case is a leaf: it never re-enters). */
static u64 lbm[(LEAFNB >> 6) + 4];
const int NW = (NB >> 6) + 2;
memset(lbm, 0, sizeof(u64) * (size_t)NW);
for (int t = 0; t < m; t++) {
const int b = (int)(((u32)(seq[t] >> 43)) >> (LFSH + LFD));
lbm[b >> 6] |= 1ull << (b & 63);
cnt[b]++;
}
u32 acc = 0;
for (int w = 0; w < NW; w++) {
u64 x = lbm[w];
while (x) {
const int b = (w << 6) + __builtin_ctzll(x);
const u32 c = cnt[b];
cnt[b] = acc; acc += c;
x &= x - 1;
}
}
for (int t = 0; t < m; t++) { u64 v = seq[t]; out[cnt[(int)(((u32)(v >> 43)) >> (LFSH + LFD))]++] = v; }
int lo = 0;
for (int w = 0; w < NW; w++) {
u64 x = lbm[w];
while (x) {
const int b = (w << 6) + __builtin_ctzll(x);
const int hi = (int)cnt[b];
if (hi - lo > 1)
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
int q = t - 1;
while (q >= lo && v < out[q]) { out[q + 1] = out[q]; q--; }
out[q + 1] = v;
}
lo = hi;
cnt[b] = 0; /* clear via the bitmap: no 32 KB memset */
x &= x - 1;
}
}
}
#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;
}
}
/* [XE] cross construction: the rekey key comes from the event's OWN embedded field
(bits 21..41 == XD[d+2][id]) instead of gathering XN[id], and the cross EMBEDS
XE[id] = XD[d+3][id] into bits 21..41 for the next level. XE == 0 at the levels that
do not rekey, so the null test folds away (XEMB is loop-invariant). */
#define MKE_CROSS(V) (((RK) ? (((((V) >> 21) & 0x1FFFFFull) << 43) | ((V) & MASK)) \
: ((V) & MASK)) \
| (XEMB ? ((u64)XEMB[(u32)(V) & 0x1FFFFFu] << 21) : 0ull))
template<int NO, int RK, int DD>
static inline void merge_loops(const u64 *retA, int h, const u64 *retB, int r,
u64 *out, u64 *cross, const u32 *XN, const u32 *XE, int &cx) {
/* Tight pointer-based merge. The stock body spends 34 uops per merged event:
`t` and the `out`/`cross` base pointers get SPILLED every iteration
(`movl %r9d,36(%rsp)` / `movq 8(%rsp),%rcx`), the 43-bit mask constant is
re-materialised with `movabsq` inside the loop, and the same `cmpq` is executed
twice (once for setb, once for cmovbe). Indices are also sign-extended twice
per event. Measured: the merge is 21.3 % of the row (idempotent duplication
probe, judge-side). Here: pointers only, no `t` (t == ia+ib), the mask lives
in a register, and `keep` is one xor. */
const u64 MASK = 4398048608255ull;
const u32 *const XEMB = XE; /* [XE] loop-invariant (may be null) */
const u64 *pa = retA, * const pa_end = retA + h;
const u64 *pb = retB, * const pb_end = retB + r;
u64 *po = out;
u64 *pc = cross;
while (pa < pa_end && pb < pb_end) {
const u64 ea = *pa, eb = *pb;
const u32 takeA = (u32)(ea < eb);
pa += takeA; pb += 1 - takeA;
const u64 v = takeA ? ea : eb;
if (NO) *po = v;
++po;
u64 c = MKE_CROSS(v);
*pc = c;
/* keep = 1 - (takeA ^ role) = 1 - (takeA ^ ((v>>42)&1)) -- as one xor pair */
pc += (size_t)(1 - ((takeA ^ (u32)(v >> 42)) & 1u));
}
while (pa < pa_end) {
const u64 v = *pa++;
if (NO) *po = v;
++po;
if ((v >> 42) & 1) {
*pc++ = MKE_CROSS(v);
}
}
while (pb < pb_end) {
const u64 v = *pb++;
if (NO) *po = v;
++po;
if (!((v >> 42) & 1)) {
*pc++ = MKE_CROSS(v);
}
}
cx = (int)(pc - cross);
}
// ---- ND specialisation dispatcher (the fleet's `DN` recipe) ----
static void base_case(u64 *seq, int m, int d, u64 *out, int needOut);
static void solve(u64 *seq, int m, int d, u64 *out, int needOut) {
if (m <= 1) return;
if (d >= K - 1) {
int cnt = 0;
for (int t = 0; t < m; t++) {
u64 v = seq[t];
if ((v >> 42) & 1) cnt++; // insert (bit31 = 1)
else OUT[eid(v)] += (u32)cnt; // query
}
return;
}
if (m <= KGB_(d)) { base_case(seq, m, d, out, needOut); return; }
(void)0;
int h = m >> 1, r = m - h;
u64 *save = ABtop;
u64 *retA = ABtop; ABtop += h;
u64 *retB = ABtop; ABtop += r;
solve(seq, h, d, retA, 1);
solve(seq + h, r, d, retB, 1);
u64 *cross = ABtop; ABtop += m + 1;
const int rekey = (d + 2 <= K - 1);
const u32 *XN = XD[rekey ? d + 2 : 0];
/* [XE] the level-(d+1) rekey/c0 needs XD[d+3]; embed it in this level's cross entries */
const u32 *XE = (d + 3 <= K - 1) ? XD[d + 3] : 0;
int cx = 0; /* ia/ib/t live inside merge_loops */
if (needOut) {
switch (d) {
case 0: merge_loops<1, 1, 0>(retA, h, retB, r, out, cross, XN, XE, cx); break;
case 1: merge_loops<1, 1, 1>(retA, h, retB, r, out, cross, XN, XE, cx); break;
case 2: merge_loops<1, 1, 2>(retA, h, retB, r, out, cross, XN, XE, cx); break;
case 3: merge_loops<1, 0, 3>(retA, h, retB, r, out, cross, XN, XE, cx); break;
default: merge_loops<1, 0, 4>(retA, h, retB, r, out, cross, XN, XE, cx); break;
}
} else {
switch (d) {
case 0: merge_loops<0, 1, 0>(retA, h, retB, r, out, cross, XN, XE, cx); break;
case 1: merge_loops<0, 1, 1>(retA, h, retB, r, out, cross, XN, XE, cx); break;
case 2: merge_loops<0, 1, 2>(retA, h, retB, r, out, cross, XN, XE, cx); break;
case 3: merge_loops<0, 0, 3>(retA, h, retB, r, out, cross, XN, XE, cx); break;
default: merge_loops<0, 0, 4>(retA, h, retB, r, out, cross, XN, XE, cx); break;
}
}
if (FUSE_OK && d == K - 3) {
u32 *dirt = (u32 *)ABtop; ABtop += (cx + 3) / 2 + 1;
ctr_solve(cross, cx, dirt);
} else {
u64 *cres = ABtop; ABtop += cx;
solve(cross, cx, d + 1, cres, 0);
}
ABtop = save;
}
static u32 *g_ord, *g_x;
static void sort0(int n) {
const u32 *X0 = XD[0];
u32 *cnt = g_x;
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;
}
void run(int n, const u32 **x, int k, u32 *out) {
LFSH = 0; LFD = 0;
K = k;
for (int d = 0; d < k; d++) XD[d] = x[d];
OUT = out;
ABtop = AB;
if (n <= 1) return;
FUSE_OK = 0;
if (K >= 4 && n > 2 && (u32)n < CT_MAXV) {
const u32 *XL = XD[K - 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);
FUSE_OK = ok;
}
{
int bit = 0;
while ((n >> bit) >= LEAFNB && bit < 30) bit++;
LFSH = bit; LFD = 0;
}
sort0(n);
const u32 *X0 = XD[0], *X1 = XD[1];
const u32 *nk2 = (K >= 3) ? XD[2] : X0; /* dim 2's coordinate, zero-only when K < 3 */
u64 *ev = AB; // top-level event array: 2n entries
int m = 0, g = 0;
for (int v = 0; v < n; ++v) {
int g1 = (int)g_x[v];
if (g1 == g) continue;
int cnt = g1 - g;
for (int t = g; t < g1; ++t) {
int id = (int)g_ord[t];
u32 key = X1[id], nk = nk2[id];
ev[m + t - g] = mk(key, 1, id, nk);
ev[m + cnt + t - g] = mk(key, 0, id, nk);
}
m += cnt * 2;
g = g1;
}
ABtop = AB + m;
u64 *res = ABtop; ABtop += m;
solve(ev, m, 0, res, 0);
}
} // 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::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(n, (const unsigned **)x, 4, out);
}
void count_5d(int n, const unsigned *x[5], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run(n, (const unsigned **)x, 5, out);
}
#pragma GCC pop_options
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 1.931 s | 93 MB + 220 KB | Accepted | Score: 100 | 显示更多 |