/* 1011a -- 测测你的五维数点2 (k=5, n=1e5, 函数接口) [a11x-STAIR]
* ============================ REFERENCES ============================
* [R1] 原作者账号 : saffah_cc_v41_260924
* 原提交 URL : https://duck.ac/submission/120327 (Time 51.903594 ms)
* 所用内容 : 本文件正文与该提交**逐字相同**的部分 = 其全部引擎骨架：
* kg:: 下的 rank_compress16 / rank_compress / rank_compress32 /
* g_rk(16 位稠密秩) / g_rk32 / g_rcnt / g_rcnt16 / XDR / XW / hsumw256 /
* VMASKW / W16 门 / CT0 / CT_MAXV / FUSE_OK / ctr_solve /
* base_case_t<ND> / base_case_t16<ND> 的叶数组(SKA/SCA/SCB/SCC)与整段 build、
* kg_mmain（手写 branchless 归并汇编，含 retA/retB 的 ~0ull 哨兵）、
* solve<W> / sort0 / run / count_5d 外壳，以及文件顶部的
* #define（KG_BASE/KG_D01/KG_D2/LEAFNB/MF_* 等）与 #pragma 集合。
* ⇒ 本件 = #120327 原文 + 且仅 + 一处"叶内分组形状"的替换（见 [R2] 与下方思路）。
* [R2] 本队账号 : saffah_cc_v41_agg1
* 原提交 URL : https://duck.ac/submission/119351 (Time 52.648378 ms，本队本题原在册最好件)
* 所用内容 : base_case_t16 的**叶内分组循环体**（"4/3/2 阶梯"：把待查插入数按
* QN[] 的单调性逐段收窄 —— 4 路扫 [0,nA)、3 路扫 [nA,·)、2 路、1 路，
* 段尾用 L16_TAIL1/2/3 掩码宏收口），以及那三个掩码宏本身。
* （该形状溯源自本队 #108154 的 base_case_t；#119351 则是本队 #108154 引擎
* + 对手 #119319 的 16 位叶，本件未沿用其引擎骨架。）
* ============================ 本次思路 ============================
* 【要改的病】#120327 的叶对一组 G 个查询先取 n = min(QN[q..q+G-1]) 跑公共 SIMD 主循环，
* 然后**每个成员各自**从 n 续到自己那一端，各跑一条**独立的 1 路**尾扫描。因为叶的 build 按
* nq 递增写入 ni、而 ni 只增不减 ⇒ QN[] 单调不减 ⇒ 组内 n = QN[q]（第一个），于是那几条尾巴
* 分别是 [nA,nB)、[nA,nC)、[nA,nD) —— 三段**互相重叠**的 1 路扫描，同一批插入被反复载入。
* 【改法】让区间不再重叠、而是像阶梯一样逐段收窄（这正是本队叶的老形状）：
* 4 路扫 [0, nA) → 3 路扫 [j, nB) → 2 路扫 [j, nC) → 1 路扫 [j, nD)，j 逐段承接上一段终点。
* 2 路与 3 路两档共用同一次 16 宽载入，尾巴长度不变但宽度变大 ⇒ 迭代次数约为原来的 43%。
* 【为什么逐位不变】查询 q 被计入的插入前缀集合仍是 [0, QN[q])，每次用的比较仍是同一批
* _mm256_cmpgt_epi16 与同一张 VMASKW 掩码；变的只是"哪几个查询共享同一次 16 宽载入"。
* 子树里的比较总数不变、累加顺序按查询分组不变 ⇒ OUT[] 逐位不变。
* 【闸门（本机 g++ -O2 -std=c++17，与判题机同源）】
* n=3000 暴力：brute_bad=0，FNV=518dd28116438b0a（与 #120327 原值逐位相同）；
* n=100000：FNV=0a33d8eaa7f40d10（本族独立 O(n^2) 暴力验过的生产规模真值）。
* 判题机「立即体验」同样 FNV=0a33d8eaa7f40d10（本件与基座同）。
* 【定价（判题机 rdtsc，best-of-5，同一驱动器）】
* #120327 = 186228364 拍 → 本件 = 181302780 拍 = -2.65%
* （本机同驱动器 best-of-5：134.6M → 126.2M 拍；本机噪声大，只作方向参考。）
*/
#ifdef Q4A_DIAG
#include <stdio.h>
#endif
#define KG_BASE 6400
#ifndef MF_CTRPF
#define MF_CTRPF 0
#endif
#ifndef MF_CTR5
#define MF_CTR5 0
#endif
#ifndef MF_LEAFLP
#define MF_LEAFLP 16
#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 2560
#endif
#ifndef LEAFNB
#define LEAFNB 8192
#endif
#define KG_D01 3456
#define KG_D2 6250 // per-level cutoff: measured -1.25 %
#define KG_MAXK 5
#define KG_ARENA (12u << 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 300020u
static u8 CT0[CT_MAXV];
static u16 CT1[CT_MAXV / 32 + 8];
static u16 CT2[CT_MAXV / 256 + 8];
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);
}
static inline void ctr_ins(u32 v) {
CT0[v]++; CT1[v >> 5]++; CT2[v >> 8]++; CT3[v >> 11]++; CT4[v >> 14]++; CT5[v >> 17]++;
}
/* ARM `CT`: 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 -- and they are
only ever read through ctr_qry, which sums them; restoring them to their initial
zero is a pure RESET, not an inverse walk. Zeroing all five together is one
contiguous memset of ~21 KB per fusion call, against 5*nd scattered decrements.
CT0 stays a per-element decrement because it is the one level that is genuinely
large (293 KB) and densely addressed, so walking `dirt` still beats clearing it.
Semantics are unchanged: ctr_qry reads exactly 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 int RSZN;
static inline void ctr_reset_hi(void) {
/* era-correct extent: n/div + 8, the declaration formula with n in place of CT_MAXV */
memset(CT1, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 5) + 8));
memset(CT2, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 8) + 8));
memset(CT3, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 11) + 8));
memset(CT4, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 14) + 8));
memset(CT5, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 17) + 8));
}
static inline u32 ctr_qry(u32 v) {
u32 k5 = v >> 17, s = 0;
#if MF_CTR5
{
__m256i tv = _mm256_loadu_si256((const __m256i *)CT5);
__m256i tm = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)k5),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
__m128i ta = _mm_add_epi32(_mm256_castsi256_si128(_mm256_and_si256(tv, tm)),
_mm256_extracti128_si256(_mm256_and_si256(tv, tm), 1));
ta = _mm_add_epi32(ta, _mm_shuffle_epi32(ta, _MM_SHUFFLE(1, 0, 3, 2)));
ta = _mm_add_epi32(ta, _mm_shuffle_epi32(ta, _MM_SHUFFLE(2, 3, 0, 1)));
s += (u32)_mm_cvtsi128_si32(ta);
}
#else
for (u32 j = 0; j < k5; j++) s += CT5[j];
#endif
/* ARM `QFUS`: one horizontal reduction for all five levels.
The boarded form ran a separate 3-shuffle reduction inside each of the five
su*_8 helpers; the five results are merely ADDED, so the reduction can be done
once at the end at no semantic cost (all levels are nonneg and their sum is the
count, bounded by n <= CT_MAXV-1, so 32-bit lanes cannot overflow). Work count
is unchanged per level; five serial shuffle chains become one. */
__m128i acc = _mm_setzero_si128();
{
__m256i tv = _mm256_loadu_si256((const __m256i *)(CT4 + ((v >> 17) << 3)));
__m256i tm = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)((v >> 14) & 7)),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
__m256i x = _mm256_and_si256(tv, tm);
acc = _mm_add_epi32(acc, _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1)));
}
{
__m256i tv = _mm256_loadu_si256((const __m256i *)(CT3 + ((v >> 14) << 3)));
__m256i tm = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)((v >> 11) & 7)),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
__m256i x = _mm256_and_si256(tv, tm);
acc = _mm_add_epi32(acc, _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1)));
}
{
__m128i vv = _mm_loadu_si128((const __m128i *)(CT2 + ((v >> 11) << 3)));
__m128i x = _mm_and_si128(vv, _mm_cmpgt_epi16(_mm_set1_epi16((short)((v >> 8) & 7)),
_mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7)));
acc = _mm_add_epi32(acc, _mm_madd_epi16(x, _mm_set1_epi16(1)));
}
{
__m128i vv = _mm_loadu_si128((const __m128i *)(CT1 + ((v >> 8) << 3)));
__m128i x = _mm_and_si128(vv, _mm_cmpgt_epi16(_mm_set1_epi16((short)((v >> 5) & 7)),
_mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7)));
acc = _mm_add_epi32(acc, _mm_madd_epi16(x, _mm_set1_epi16(1)));
}
{
__m256i vv = _mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u)));
__m256i m = _mm256_cmpgt_epi8(_mm256_set1_epi8((char)(v & 31)),
_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 sd = _mm256_sad_epu8(_mm256_and_si256(vv, m), _mm256_setzero_si256());
acc = _mm_add_epi32(acc, _mm_add_epi32(_mm256_castsi256_si128(sd), _mm256_extracti128_si256(sd, 1)));
}
acc = _mm_add_epi32(acc, _mm_shuffle_epi32(acc, _MM_SHUFFLE(1, 0, 3, 2)));
acc = _mm_add_epi32(acc, _mm_shuffle_epi32(acc, _MM_SHUFFLE(2, 3, 0, 1)));
return s + (u32)_mm_cvtsi128_si32(acc);
}
static int K;
#ifndef LFAUTO
#define LFAUTO 0
#endif
#ifndef LF_D
#define LF_D 0
#endif
static int LFSH;
static int LFB;
static int LFD;
static const u32 *XD[KG_MAXK];
static const u16 *XDR[KG_MAXK];
static u32 *OUT;
static u64 *AB, *ABtop;
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 1 - (takeA ^ role); }
static inline u64 mk(u32 key, int q, int id) {
return ((u64)key << 32) | ((u64)(1 - q) << 31) | (u64)id;
}
template<int W> static void solve(u64 *seq, int m, int d, u64 *out, int needOut);
/* asm4: sentinel-terminated merge main loop, 75 statements on ONE line. */
void kg_mmain(const u64*,long,const u64*,long,u64*,u64*,const u16*,int*)
__asm__("kg_mmain");
__asm__(" .text\n .globl kg_mmain\n .type kg_mmain,@function\nkg_mmain:\n pushq %rbx\n pushq %rbp\n pushq %r12\n pushq %r13\n pushq %r14\n pushq %r15\n subq $48, %rsp\n movq 104(%rsp), %r10\n movq 112(%rsp), %rax\n movq %rax, 0(%rsp)\n movq %rdi, 8(%rsp)\n movq %rdx, 16(%rsp)\n movq %r8, 24(%rsp)\n movq %r9, 32(%rsp)\n movq (%rdi), %r11\n movq (%rdx), %rax\n xorl %r15d, %r15d\n movq %rsi, %rcx\n .p2align 4\n.Lkgloop:\n movq 8(%rdi), %r12\n movq 8(%rdx), %r13\n cmpq %rax, %r11\n setb %r15b\n movq %rax, %rbp\n cmovb %r11, %rbp\n cmovb %r12, %r11\n cmovae %r13, %rax\n movl %ebp, %ebx\n andl $0x7FFFFFFF, %ebx\n movl %ebp, (%r9)\n movq %rbp, (%r8)\n movzwl (%r10,%rbx,2), %ebx\n addq $8, %r8\n movl %ebx, 4(%r9)\n shrl $31, %ebp\n xorl %r15d, %ebp\n xorl $1, %ebp\n leaq (%r9,%rbp,8), %r9\n leaq (%rdi,%r15,8), %rdi\n xorl $1, %r15d\n leaq (%rdx,%r15,8), %rdx\n subq $1, %rcx\n jnz .Lkgloop\n.Lkgdone:\n movq 0(%rsp), %rax\n movq %rdi, %rbx\n subq 8(%rsp), %rbx\n sarq $3, %rbx\n movl %ebx, 0(%rax)\n movq %rdx, %rbx\n subq 16(%rsp), %rbx\n sarq $3, %rbx\n movl %ebx, 4(%rax)\n movq %r9, %rbx\n subq 32(%rsp), %rbx\n sarq $3, %rbx\n movl %ebx, 8(%rax)\n movq %r8, %rbx\n subq 24(%rsp), %rbx\n sarq $3, %rbx\n movl %ebx, 12(%rax)\n addq $48, %rsp\n popq %r15\n popq %r14\n popq %r13\n popq %r12\n popq %rbp\n popq %rbx\n ret\n .size kg_mmain, .-kg_mmain\n");
// 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] >> 32);
__builtin_prefetch(&CT0[vf], 1, 0);
__builtin_prefetch(&CT1[vf >> 5], 1, 0);
__builtin_prefetch(&CT2[vf >> 8], 1, 0);
}
#endif
u64 v = cross[i];
u32 val = (u32)(v >> 32);
if ((v >> 31) & 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);
}
/* ===== ARM `DN`: SPECIALISE base_case ON THE REMAINING DIMENSION COUNT =====
`nd = K - (d + 2)` is a RUNTIME value, yet it gates work inside base_case's
INNERMOST loops: three `if (nd > 0/1/2)` tests sit between the AVX2 mask folds
(the loads of X0/X1/X2 and the `_mm256_cmpgt_epi32` chain), so every 8-element
block pays three compare-and-branch on a value that cannot change within the
call. `nd` is in [0, K-2] and K <= KG_MAXK = 5, so there are only FOUR cases;
dispatching once per base_case call and compiling `nd` as a constant removes
every one of those branches and lets the address selection of X0/X1/X2 fold too.
This is a WORK reduction -- instructions are deleted, not rescheduled -- which
is the class whose local sign has held on this row. */
// ===== 16-bit ranked leaf (ported from sibling 1010a lane q4a, md5 3f1021ffc59c) =====
static u16 g_rk[KG_MAXK][300005];
static u16 g_rcnt16[300005];
static u32 g_rk32[KG_MAXK][300005];
static u32 g_rcnt[300005];
static u32 g_mxk, g_mxc0, g_mxc1, g_mxc2, g_mxnd[4];
static int W16 = 0;
static u16 g_dmax;
#ifdef Q4A_DIAG
static void q4adump(void) __attribute__((destructor));
static void q4adump(void){ fprintf(stderr, "Q4A16 W16=%d dmax=%u mxk=%u mxc0=%u mxc1=%u mxc2=%u nd[%u %u %u %u]\n",
W16, (unsigned)g_dmax, g_mxk, g_mxc0, g_mxc1, g_mxc2,
g_mxnd[0], g_mxnd[1], g_mxnd[2], g_mxnd[3]); }
#endif // 1 => the u16 leaf is admissible for this input
static inline u16 XW(u32 x) { return (u16)(x ^ 0x8000u); }
static const u16 VMASKW[17][16] = {
{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFFFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF}};
static inline u32 hsumw256(__m256i v) { // 16 x i16 -> scalar
__m256i w = _mm256_madd_epi16(v, _mm256_set1_epi16(1));
return hsum256(w);
}
// order-preserving rank of each coordinate, per dimension.
// DENSE rank: number of DISTINCT values strictly below v. A count-based prefix (the obvious
// form) is also order-preserving but ranges over [0, n-1], which does NOT fit in 16 bits once
// n > 65536 -- the sibling's first bug, kept here as a warning.
// Dim 0 is skipped: it is consumed only by sort0 (order) and the top-level grouping (equality),
// and the leaf's smallest array index is XD[d+2] >= XD[2] for KT>=4 -- but XD[1] IS compared
// (it is the carried key at d=0), so the sibling's exact loop bound (d >= 1) is kept.
static u32 rank_compress16(int n, int KT) {
u16 *cnt = g_rcnt16;
u32 dmax = 0;
for (int d = 1; d < KT; d++) {
const u32 *xs = XD[d];
memset(cnt, 0, sizeof(u16) * (size_t)(n + 1));
int ok = 1;
for (int i = 0; i < n; i++) { u16 c0 = cnt[xs[i]]; if (c0 == 0xFFFFu) { ok = 0; break; } cnt[xs[i]] = (u16)(c0 + 1); }
if (!ok) return 0;
u32 D = 0, run = 0;
for (int v = 0; v < n; v++) { u32 c = cnt[v]; cnt[v] = (u16)run; if (c) { D++; run++; } }
if (run > 65535u || D > 65535u) return 0;
if (D > dmax) dmax = D;
u16 *dst = g_rk[d];
XDR[d] = dst;
for (int i = 0; i < n; i++) dst[i] = cnt[xs[i]];
}
return dmax ? dmax : 1u;
}
static u32 rank_compress(int n, int KT) {
u32 *cnt = g_rcnt;
u32 dmax = 0;
for (int d = 1; d < KT; d++) {
const u32 *xs = XD[d];
memset(cnt, 0, sizeof(u32) * (size_t)(n + 1));
for (int i = 0; i < n; i++) cnt[xs[i]]++;
u32 D = 0, run = 0;
for (int v = 0; v < n; v++) { u32 c = cnt[v]; cnt[v] = run; if (c) { D++; run++; } }
if (D > dmax) dmax = D;
u16 *dst = g_rk[d];
XDR[d] = dst;
for (int i = 0; i < n; i++) dst[i] = (u16)cnt[xs[i]];
}
return dmax;
}
// the ORIGINAL u32 ranking, redone only when dmax > 65536 (W16 == 0)
static u32 rank_compress32(int n, int KT) {
u32 *cnt = g_rcnt;
u32 dmax = 0;
for (int d = 1; d < KT; d++) {
const u32 *xs = XD[d];
memset(cnt, 0, sizeof(u32) * (size_t)(n + 1));
for (int i = 0; i < n; i++) cnt[xs[i]]++;
u32 D = 0, run = 0;
for (int v = 0; v < n; v++) { u32 c = cnt[v]; cnt[v] = run; if (c) { D++; run++; } }
if (D > dmax) dmax = D;
u32 *dst = g_rk32[d];
for (int i = 0; i < n; i++) dst[i] = cnt[xs[i]];
XD[d] = dst;
}
return dmax;
}
// =========================================================
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];
u32 QK[KG_BASE + 8], QC0[KG_BASE + 8], QC1[KG_BASE + 8], QC2[KG_BASE + 8];
u32 QID[KG_BASE + 8], QN[KG_BASE + 8];
const u32 *X0 = nd > 0 ? XD[d + 2] : 0, *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 = eid(seq[t + MF_LEAFLP]);
if (nd > 0) __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 >> 32);
u32 c0 = 0, c1 = 0, c2 = 0;
if (nd > 0) c0 = X0[id];
if (nd > 1) c1 = X1[id];
if (nd > 2) c2 = X2[id];
int isins = (int)((v >> 31) & 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 ----
int q = 0;
for (; q + 2 <= nq; q += 2) {
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++) {
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) : LFB;
u16 cnt[LEAFNB + 2];
for (int b = 0; b <= NB; b++) cnt[b] = 0;
for (int t = 0; t < m; t++) cnt[(int)(((u32)(seq[t] >> 32)) >> (LFSH + LFD)) + 1]++;
for (int b = 0; b < NB; b++) cnt[b + 1] += cnt[b];
for (int t = 0; t < m; t++) { u64 v = seq[t]; out[cnt[(int)(((u32)(v >> 32)) >> (LFSH + LFD))]++] = v; }
{
// A bucket's span is [cnt[bb-1], cnt[bb]) (cnt is the post-scatter END offset),
// so an unaligned load at b-1 hands all eight `lo`s to one vector compare.
// NOTE bucket 0 is REAL: the scatter writes cnt[bucket]++, so cnt[0] ends at
// hist(0) -- it is not left zero by the prefix. (A first build skipped bucket 0
// and lost 2 counts: SSN2112, caught only by the chk gate.)
int b;
{
const int hi = (int)cnt[0];
for (int t = 1; t < hi; t++) {
u64 v = out[t];
int q = t - 1;
while (q >= 0 && v < out[q]) { out[q + 1] = out[q]; q--; }
out[q + 1] = v;
}
}
for (b = 1; b + 8 <= NB; b += 8) {
__m128i h = _mm_loadu_si128((const __m128i *)(cnt + b));
__m128i p = _mm_loadu_si128((const __m128i *)(cnt + b - 1));
__m128i d = _mm_sub_epi16(h, p);
u32 msk = _pext_u32((u32)_mm_movemask_epi8(_mm_cmpgt_epi16(d, _mm_set1_epi16(1))),
0x5555u);
while (msk) {
int bb = b + __builtin_ctz(msk);
int hi = (int)cnt[bb], lo = (int)cnt[bb - 1];
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
for (int q = t - 1; q >= lo; q--) {
u64 c = out[q];
const int lt = (c < v);
out[q + 1] = lt ? v : c;
v = lt ? c : v;
}
out[lo] = v;
}
msk &= msk - 1;
}
}
for (; b < NB; b++) {
int hi = (int)cnt[b], lo = (int)cnt[b - 1];
if (hi - lo > 1)
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
for (int q = t - 1; q >= lo; q--) {
u64 c = out[q];
const int lt = (c < v);
out[q + 1] = lt ? v : c;
v = lt ? c : v;
}
out[lo] = v;
}
}
}
}
#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
}
#define LDAW(p) _mm256_load_si256((const __m256i *)(p))
#define LDUW(k) _mm256_loadu_si256((const __m256i *)(VMASKW[k]))
#define L16_TAIL1(ACC, KVQ, C0Q, NX, JS) do { \
int jj = (JS); \
for (; jj + 16 <= (NX); jj += 16) { \
__m256i msk = _mm256_cmpgt_epi16(KVQ, LDAW(IK + jj)); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C0Q, LDAW(IC0 + jj))); \
ACC = _mm256_sub_epi16(ACC, msk); \
} \
if (jj < (NX)) { \
const int nw_ = (NX); \
__m256i msk = _mm256_cmpgt_epi16(KVQ, LDAW(IK + jj)); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C0Q, LDAW(IC0 + jj))); \
msk = _mm256_and_si256(msk, LDUW(nw_ - jj)); \
ACC = _mm256_sub_epi16(ACC, msk); \
} \
} while (0)
#define L16_TAIL2(ACC, KVQ, C0Q, C1Q, NX, JS) do { \
int jj = (JS); \
for (; jj + 16 <= (NX); jj += 16) { \
__m256i msk = _mm256_cmpgt_epi16(KVQ, LDAW(IK + jj)); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C0Q, LDAW(IC0 + jj))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C1Q, LDAW(IC1 + jj))); \
ACC = _mm256_sub_epi16(ACC, msk); \
} \
if (jj < (NX)) { \
const int nw_ = (NX); \
__m256i msk = _mm256_cmpgt_epi16(KVQ, LDAW(IK + jj)); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C0Q, LDAW(IC0 + jj))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C1Q, LDAW(IC1 + jj))); \
msk = _mm256_and_si256(msk, LDUW(nw_ - jj)); \
ACC = _mm256_sub_epi16(ACC, msk); \
} \
} while (0)
#define L16_TAIL3(ACC, KVQ, C0Q, C1Q, C2Q, NX, JS) do { \
int jj = (JS); \
for (; jj + 16 <= (NX); jj += 16) { \
__m256i msk = _mm256_cmpgt_epi16(KVQ, LDAW(IK + jj)); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C0Q, LDAW(IC0 + jj))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C1Q, LDAW(IC1 + jj))); \
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C2Q, LDAW(IC2 + jj))); \
ACC = _mm256_sub_epi16(ACC, msk); \
} \
if (jj < (NX)) { \
const int nw_ = (NX); \
__m256i msk = _mm256_cmpgt_epi16(KVQ, LDAW(IK + jj)); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C0Q, LDAW(IC0 + jj))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C1Q, LDAW(IC1 + jj))); \
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(C2Q, LDAW(IC2 + jj))); \
msk = _mm256_and_si256(msk, LDUW(nw_ - jj)); \
ACC = _mm256_sub_epi16(ACC, msk); \
} \
} while (0)
template<int ND>
static void base_case_t16(u64 *seq, int m, int d, u64 *out, int needOut) {
const int nd = ND;
__attribute__((aligned(32)))
u16 SKA[2 * (KG_BASE + 16)], SCA[2 * (KG_BASE + 16)],
SCB[2 * (KG_BASE + 16)], SCC[2 * (KG_BASE + 16)];
u32 QID[KG_BASE + 8], QN[KG_BASE + 8];
const u16 *X0 = nd > 0 ? XDR[d + 2] : 0, *X1 = nd > 1 ? XDR[d + 3] : 0, *X2 = nd > 2 ? XDR[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 = eid(seq[t + MF_LEAFLP]);
if (nd > 0) __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 >> 32);
u32 c0 = 0, c1 = 0, c2 = 0;
if (nd > 0) c0 = X0[id];
if (nd > 1) c1 = X1[id];
if (nd > 2) c2 = X2[id];
int isins = (int)((v >> 31) & 1u);
#ifdef Q4A_DIAG
if (nd < 4 && ++g_mxnd[nd] == 0) {}
if (k > g_mxk) g_mxk = k; if (c0 > g_mxc0) g_mxc0 = c0; if (c1 > g_mxc1) g_mxc1 = c1; if (c2 > g_mxc2) g_mxc2 = c2;
#endif
{ const int sx = isins ? ni : (KG_BASE + 16) + nq;
SKA[sx] = XW(k); SCA[sx] = XW(c0);
if (nd > 1) SCB[sx] = XW(c1);
if (nd > 2) SCC[sx] = XW(c2); }
QID[nq] = (u32)id; QN[nq] = (u32)ni; ni += isins; nq += 1 - isins;
}
// [a11x-STAIR] aliases: their merged SKA/SCA/SCB/SCC -> our IK/IC0/IC1/IC2 + query window
const u16 *const IK = SKA, *const IC0 = SCA, *const IC1 = SCB, *const IC2 = SCC;
const u16 *const QK = SKA + (KG_BASE + 16), *const QC0 = SCA + (KG_BASE + 16),
*const QC1 = SCB + (KG_BASE + 16), *const QC2 = SCC + (KG_BASE + 16);
int q = 0;
#if MF_LEAFG >= 4
/* ---- 4-wide staircase leaf, nd <= 1 (A = nd+1 <= 2 insert arrays) ---- */
if (nd <= 1) {
for (; q + 4 <= nq; q += 4) {
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 0]])));
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 1]])));
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 2]])));
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 3]])));
const int nA4 = (int)QN[q], nB4 = (int)QN[q + 1],
nC4 = (int)QN[q + 2], nD4 = (int)QN[q + 3];
__m256i kA4 = _mm256_set1_epi16((short)QK[q]), dA4 = _mm256_set1_epi16((short)QC0[q]),
aA4 = _mm256_setzero_si256();
__m256i kB4 = _mm256_set1_epi16((short)QK[q + 1]), dB4 = _mm256_set1_epi16((short)QC0[q + 1]),
aB4 = _mm256_setzero_si256();
__m256i kC4 = _mm256_set1_epi16((short)QK[q + 2]), dC4 = _mm256_set1_epi16((short)QC0[q + 2]),
aC4 = _mm256_setzero_si256();
__m256i kD4 = _mm256_set1_epi16((short)QK[q + 3]), dD4 = _mm256_set1_epi16((short)QC0[q + 3]),
aD4 = _mm256_setzero_si256();
int j4 = 0;
for (; j4 + 16 <= nA4; j4 += 16) { /* A, B, C, D */
__m256i IKv = LDAW(IK + j4);
__m256i mA = _mm256_cmpgt_epi16(kA4, IKv), mB = _mm256_cmpgt_epi16(kB4, IKv),
mC = _mm256_cmpgt_epi16(kC4, IKv), mD = _mm256_cmpgt_epi16(kD4, IKv);
if (nd > 0) {
__m256i C0v = LDAW(IC0 + j4);
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(dA4, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(dB4, C0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(dC4, C0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi16(dD4, C0v));
}
aA4 = _mm256_sub_epi16(aA4, mA); aB4 = _mm256_sub_epi16(aB4, mB);
aC4 = _mm256_sub_epi16(aC4, mC); aD4 = _mm256_sub_epi16(aD4, mD);
}
const int pA4 = j4;
for (; j4 + 16 <= nB4; j4 += 16) { /* B, C, D */
__m256i IKv = LDAW(IK + j4);
__m256i mB = _mm256_cmpgt_epi16(kB4, IKv), mC = _mm256_cmpgt_epi16(kC4, IKv),
mD = _mm256_cmpgt_epi16(kD4, IKv);
if (nd > 0) {
__m256i C0v = LDAW(IC0 + j4);
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(dB4, C0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(dC4, C0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi16(dD4, C0v));
}
aB4 = _mm256_sub_epi16(aB4, mB); aC4 = _mm256_sub_epi16(aC4, mC);
aD4 = _mm256_sub_epi16(aD4, mD);
}
const int pB4 = j4;
for (; j4 + 16 <= nC4; j4 += 16) { /* C, D */
__m256i IKv = LDAW(IK + j4);
__m256i mC = _mm256_cmpgt_epi16(kC4, IKv), mD = _mm256_cmpgt_epi16(kD4, IKv);
if (nd > 0) {
__m256i C0v = LDAW(IC0 + j4);
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(dC4, C0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi16(dD4, C0v));
}
aC4 = _mm256_sub_epi16(aC4, mC); aD4 = _mm256_sub_epi16(aD4, mD);
}
const int pC4 = j4;
for (; j4 + 16 <= nD4; j4 += 16) { /* D only */
__m256i mD = _mm256_cmpgt_epi16(kD4, LDAW(IK + j4));
if (nd > 0) mD = _mm256_and_si256(mD, _mm256_cmpgt_epi16(dD4, LDAW(IC0 + j4)));
aD4 = _mm256_sub_epi16(aD4, mD);
}
const int pD4 = j4;
L16_TAIL1(aA4, kA4, dA4, nA4, pA4);
L16_TAIL1(aB4, kB4, dB4, nB4, pB4);
L16_TAIL1(aC4, kC4, dC4, nC4, pC4);
L16_TAIL1(aD4, kD4, dD4, nD4, pD4);
OUT[(u32)QID[q]] += hsumw256(aA4);
OUT[(u32)QID[q + 1]] += hsumw256(aB4);
OUT[(u32)QID[q + 2]] += hsumw256(aC4);
OUT[(u32)QID[q + 3]] += hsumw256(aD4);
}
}
#endif
#if MF_LEAFG >= 3
/* ---- 3-wide staircase leaf, nd <= 2 (A = nd+1 <= 3 insert arrays) ---- */
if (nd <= 2) {
for (; q + 3 <= nq; q += 3) {
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 0]])));
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 1]])));
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 2]])));
const int nA3 = (int)QN[q], nB3 = (int)QN[q + 1], nC3 = (int)QN[q + 2];
__m256i kA3 = _mm256_set1_epi16((short)QK[q]), dA3 = _mm256_set1_epi16((short)QC0[q]),
eA3 = _mm256_set1_epi16((short)QC1[q]), aA3 = _mm256_setzero_si256();
__m256i kB3 = _mm256_set1_epi16((short)QK[q + 1]), dB3 = _mm256_set1_epi16((short)QC0[q + 1]),
eB3 = _mm256_set1_epi16((short)QC1[q + 1]), aB3 = _mm256_setzero_si256();
__m256i kC3 = _mm256_set1_epi16((short)QK[q + 2]), dC3 = _mm256_set1_epi16((short)QC0[q + 2]),
eC3 = _mm256_set1_epi16((short)QC1[q + 2]), aC3 = _mm256_setzero_si256();
int j3 = 0;
for (; j3 + 16 <= nA3; j3 += 16) { /* A, B, C */
__m256i IKv = LDAW(IK + j3);
__m256i mA = _mm256_cmpgt_epi16(kA3, IKv), mB = _mm256_cmpgt_epi16(kB3, IKv),
mC = _mm256_cmpgt_epi16(kC3, IKv);
if (nd > 0) {
__m256i C0v = LDAW(IC0 + j3);
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(dA3, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(dB3, C0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(dC3, C0v));
}
if (nd > 1) {
__m256i C1v = LDAW(IC1 + j3);
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(eA3, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(eB3, C1v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(eC3, C1v));
}
aA3 = _mm256_sub_epi16(aA3, mA); aB3 = _mm256_sub_epi16(aB3, mB);
aC3 = _mm256_sub_epi16(aC3, mC);
}
const int pA3 = j3;
for (; j3 + 16 <= nB3; j3 += 16) { /* B, C */
__m256i IKv = LDAW(IK + j3);
__m256i mB = _mm256_cmpgt_epi16(kB3, IKv), mC = _mm256_cmpgt_epi16(kC3, IKv);
if (nd > 0) {
__m256i C0v = LDAW(IC0 + j3);
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(dB3, C0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(dC3, C0v));
}
if (nd > 1) {
__m256i C1v = LDAW(IC1 + j3);
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(eB3, C1v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(eC3, C1v));
}
aB3 = _mm256_sub_epi16(aB3, mB); aC3 = _mm256_sub_epi16(aC3, mC);
}
const int pB3 = j3;
for (; j3 + 16 <= nC3; j3 += 16) { /* C only */
__m256i mC = _mm256_cmpgt_epi16(kC3, LDAW(IK + j3));
if (nd > 0) mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(dC3, LDAW(IC0 + j3)));
if (nd > 1) mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(eC3, LDAW(IC1 + j3)));
aC3 = _mm256_sub_epi16(aC3, mC);
}
const int pC3 = j3;
L16_TAIL2(aA3, kA3, dA3, eA3, nA3, pA3);
L16_TAIL2(aB3, kB3, dB3, eB3, nB3, pB3);
L16_TAIL2(aC3, kC3, dC3, eC3, nC3, pC3);
OUT[(u32)QID[q]] += hsumw256(aA3);
OUT[(u32)QID[q + 1]] += hsumw256(aB3);
OUT[(u32)QID[q + 2]] += hsumw256(aC3);
}
}
#endif
#if MF_LEAFG >= 2
/* ---- paired-query leaf: one insert-block load serves two queries ---- */
for (; q + 2 <= nq; q += 2) {
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 0]])));
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&OUT[(u32)QID[q + 1]])));
const int nA = (int)QN[q], nB = (int)QN[q + 1];
const int nn = nA < nB ? nA : nB;
__m256i kvA = _mm256_set1_epi16((short)QK[q]), kvB = _mm256_set1_epi16((short)QK[q + 1]);
__m256i v0A = _mm256_set1_epi16((short)QC0[q]), v0B = _mm256_set1_epi16((short)QC0[q + 1]);
__m256i v1A = _mm256_set1_epi16((short)QC1[q]), v1B = _mm256_set1_epi16((short)QC1[q + 1]);
__m256i v2A = _mm256_set1_epi16((short)QC2[q]), v2B = _mm256_set1_epi16((short)QC2[q + 1]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256();
int j = 0;
if (nd >= 3) {
for (; j + 16 <= nn; j += 16) {
__m256i IKv = LDAW(IK + j);
__m256i C0v = LDAW(IC0 + j);
__m256i C1v = LDAW(IC1 + j);
__m256i C2v = LDAW(IC2 + j);
__m256i mA = _mm256_and_si256(_mm256_cmpgt_epi16(kvA, IKv), _mm256_cmpgt_epi16(v0A, C0v));
__m256i mB = _mm256_and_si256(_mm256_cmpgt_epi16(kvB, IKv), _mm256_cmpgt_epi16(v0B, C0v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(v1B, C1v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(v2A, C2v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(v2B, C2v));
accA = _mm256_sub_epi16(accA, mA);
accB = _mm256_sub_epi16(accB, mB);
}
} else {
for (; j + 16 <= nn; j += 16) {
__m256i IKv = LDAW(IK + j);
__m256i mA = _mm256_cmpgt_epi16(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi16(kvB, IKv);
if (nd > 0) {
__m256i C0v = LDAW(IC0 + j);
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(v0A, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(v0B, C0v));
}
if (nd > 1) {
__m256i C1v = LDAW(IC1 + j);
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(v1B, C1v));
}
accA = _mm256_sub_epi16(accA, mA);
accB = _mm256_sub_epi16(accB, mB);
}
}
L16_TAIL3(accA, kvA, v0A, v1A, v2A, nA, j);
OUT[(u32)QID[q]] += hsumw256(accA);
L16_TAIL3(accB, kvB, v0B, v1B, v2B, nB, j);
OUT[(u32)QID[q + 1]] += hsumw256(accB);
}
#endif
for (; q < nq; q++) {
const int n = (int)QN[q];
__m256i kv = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q]);
__m256i v0 = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q]);
__m256i v1 = _mm256_set1_epi16((short)SCB[(KG_BASE + 16) + q]);
__m256i v2 = _mm256_set1_epi16((short)SCC[(KG_BASE + 16) + q]);
__m256i acc = _mm256_setzero_si256();
int j = 0;
switch (nd) {
case 0:
for (; j + 16 <= n; j += 16)
acc = _mm256_sub_epi16(acc, _mm256_cmpgt_epi16(kv, _mm256_load_si256((const __m256i *)(SKA + j))));
break;
case 1:
for (; j + 16 <= n; j += 16) {
__m256i msk = _mm256_cmpgt_epi16(kv, _mm256_load_si256((const __m256i *)(SKA + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(SCA + j))));
acc = _mm256_sub_epi16(acc, msk);
}
break;
case 2:
for (; j + 16 <= n; j += 16) {
__m256i msk = _mm256_cmpgt_epi16(kv, _mm256_load_si256((const __m256i *)(SKA + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(SCA + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1, _mm256_load_si256((const __m256i *)(SCB + j))));
acc = _mm256_sub_epi16(acc, msk);
}
break;
default:
for (; j + 16 <= n; j += 16) {
__m256i msk = _mm256_cmpgt_epi16(kv, _mm256_load_si256((const __m256i *)(SKA + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(SCA + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1, _mm256_load_si256((const __m256i *)(SCB + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2, _mm256_load_si256((const __m256i *)(SCC + j))));
acc = _mm256_sub_epi16(acc, msk);
}
break;
}
if (j < n) { // partial block, masked by a table lookup
__m256i msk = _mm256_cmpgt_epi16(kv, _mm256_load_si256((const __m256i *)(SKA + j)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(SCA + j))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1, _mm256_load_si256((const __m256i *)(SCB + j))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2, _mm256_load_si256((const __m256i *)(SCC + j))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[n - j]));
acc = _mm256_sub_epi16(acc, msk);
}
u32 c = hsumw256(acc);
OUT[QID[q]] += c;
}
if (!needOut) return;
#if LEAFNB
{
const int NB = (LFAUTO && m <= LEAFNB) ? (m < 8 ? 8 : m) : LFB;
u16 cnt[LEAFNB + 2];
for (int b = 0; b <= NB; b++) cnt[b] = 0;
for (int t = 0; t < m; t++) cnt[(int)(((u32)(seq[t] >> 32)) >> (LFSH + LFD)) + 1]++;
for (int b = 0; b < NB; b++) cnt[b + 1] += cnt[b];
for (int t = 0; t < m; t++) { u64 v = seq[t]; out[cnt[(int)(((u32)(v >> 32)) >> (LFSH + LFD))]++] = v; }
{
// A bucket's span is [cnt[bb-1], cnt[bb]) (cnt is the post-scatter END offset),
// so an unaligned load at b-1 hands all eight `lo`s to one vector compare.
// NOTE bucket 0 is REAL: the scatter writes cnt[bucket]++, so cnt[0] ends at
// hist(0) -- it is not left zero by the prefix. (A first build skipped bucket 0
// and lost 2 counts: SSN2112, caught only by the chk gate.)
int b;
{
const int hi = (int)cnt[0];
for (int t = 1; t < hi; t++) {
u64 v = out[t];
int q = t - 1;
while (q >= 0 && v < out[q]) { out[q + 1] = out[q]; q--; }
out[q + 1] = v;
}
}
for (b = 1; b + 8 <= NB; b += 8) {
__m128i h = _mm_loadu_si128((const __m128i *)(cnt + b));
__m128i p = _mm_loadu_si128((const __m128i *)(cnt + b - 1));
__m128i d = _mm_sub_epi16(h, p);
u32 msk = _pext_u32((u32)_mm_movemask_epi8(_mm_cmpgt_epi16(d, _mm_set1_epi16(1))),
0x5555u);
while (msk) {
int bb = b + __builtin_ctz(msk);
int hi = (int)cnt[bb], lo = (int)cnt[bb - 1];
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
for (int q = t - 1; q >= lo; q--) {
u64 c = out[q];
const int lt = (c < v);
out[q + 1] = lt ? v : c;
v = lt ? c : v;
}
out[lo] = v;
}
msk &= msk - 1;
}
}
for (; b < NB; b++) {
int hi = (int)cnt[b], lo = (int)cnt[b - 1];
if (hi - lo > 1)
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
for (int q = t - 1; q >= lo; q--) {
u64 c = out[q];
const int lt = (c < v);
out[q + 1] = lt ? v : c;
v = lt ? c : v;
}
out[lo] = v;
}
}
}
}
#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) {
if (W16) {
switch (K - (d + 2)) {
case 0: base_case_t16<0>(seq, m, d, out, needOut); break;
case 1: base_case_t16<1>(seq, m, d, out, needOut); break;
case 2: base_case_t16<2>(seq, m, d, out, needOut); break;
default: base_case_t16<3>(seq, m, d, out, needOut); break;
}
return;
}
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;
}
}
template<int W>
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 >> 31) & 1) cnt++; // insert (bit31 = 1)
else OUT[eid(v)] += (u32)cnt; // query
}
return;
}
if (m <= (d == 2 ? KG_D2 : KG_D01)) { 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 + 1; /* +1: the sentinel slot (asm4) */
u64 *retB = ABtop; ABtop += r + 1;
solve<W>(seq, h, d, retA, 1);
solve<W>(seq + h, r, d, retB, 1);
u64 *cross = ABtop; ABtop += m + 1;
const int rekey = (d + 2 <= K - 1);
const u16 *XN16 = W ? (rekey ? XDR[d + 2] : XDR[0]) : 0;
const u32 *XN32 = W ? 0 : (rekey ? XD[d + 2] : XD[0]);
int ia = 0, ib = 0, cx = 0, t = 0;
retA[h] = ~0ull; retB[r] = ~0ull;
if (W && rekey) {
int mres[4];
kg_mmain(retA, (long)m, retB, 0, out, cross, XN16, mres);
ia = mres[0]; ib = mres[1]; cx = mres[2]; t = mres[3];
} else
while (ia < h && ib < r) {
#if KG_PF
if (!(ia & (KG_PF - 1)) && ia + KG_PF < h) __builtin_prefetch(&XN[eid(retA[ia + KG_PF])], 0, 1);
if (!(ib & (KG_PF - 1)) && ib + KG_PF < r) __builtin_prefetch(&XN[eid(retB[ib + KG_PF])], 0, 1);
#endif
u64 ea = retA[ia], eb = retB[ib];
int takeA = (ea < eb);
ia += takeA; ib += 1 - takeA;
u64 v = takeA ? ea : eb;
if (needOut) out[t] = v;
t++;
int role = (int)((v >> 31) & 1);
#if KG_COND
if (keep_test(takeA, role)) {
cross[cx++] = rekey ? (((u64)(W ? (u32)XN16[eid(v)] : XN32[eid(v)]) << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
}
#else
u64 c = rekey ? (((u64)(W ? (u32)XN16[eid(v)] : XN32[eid(v)]) << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
cross[cx] = c;
cx += keep_test(takeA, role);
#endif
}
while (ia < h) {
u64 v = retA[ia++];
if (needOut) out[t] = v;
t++;
if ((v >> 31) & 1) cross[cx++] = rekey ? (((u64)(W ? (u32)XN16[eid(v)] : XN32[eid(v)]) << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
}
while (ib < r) {
u64 v = retB[ib++];
if (needOut) out[t] = v;
t++;
if (!((v >> 31) & 1)) cross[cx++] = rekey ? (((u64)(W ? (u32)XN16[eid(v)] : XN32[eid(v)]) << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
}
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<W>(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;
RSZN = n;
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 = LF_D;
}
{
u32 dm = rank_compress16(n, k);
if (!dm) dm = rank_compress(n, k);
g_dmax = (u16)(dm > 65535u ? 65535u : dm);
W16 = (dm <= 65536u);
if (!W16) { rank_compress32(n, k); }
#ifdef Q4A_DIAG
g_dmax = (u16)(dm > 65535u ? 65535u : dm);
#endif
LFB = (int)((dm ? dm - 1u : 0u) >> (LFSH + LFD)) + 1;
if (LFB > LEAFNB + 1) LFB = LEAFNB + 1;
if (LFB < 1) LFB = 1;
}
sort0(n);
const u32 *X0 = XD[0]; const u16 *X1r = XDR[1]; const u32 *X1u = 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++;
if (W16) {
for (int t = g; t < g1; t++) ev[m++] = mk((u32)X1r[g_ord[t]], 1, (int)g_ord[t]);
for (int t = g; t < g1; t++) ev[m++] = mk((u32)X1r[g_ord[t]], 0, (int)g_ord[t]);
} else {
for (int t = g; t < g1; t++) ev[m++] = mk(X1u[g_ord[t]], 1, (int)g_ord[t]);
for (int t = g; t < g1; t++) ev[m++] = mk(X1u[g_ord[t]], 0, (int)g_ord[t]);
}
g = g1;
}
ABtop = AB + m;
u64 *res = ABtop; ABtop += m;
if (W16) solve<1>(ev, m, 0, res, 0); else solve<0>(ev, m, 0, res, 0);
}
} // namespace kg
static u64 g_arena[KG_ARENA];
static u32 g_ord_st[300005], g_x_st[300005];
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 | 50.885 ms | 11 MB + 972 KB | Accepted | Score: 100 | 显示更多 |