#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 negq %rcx\n leaq (%r8,%rsi,8), %r8\n .p2align 4\n.Lkgloop:\n cmpq %rax, %r11\n setb %r15b\n movq %rax, %rbp\n cmovb %r11, %rbp\n cmovb 8(%rdi), %r11\n cmovae 8(%rdx), %rax\n movl %ebp, %ebx\n andl $0x7FFFFFFF, %ebx\n movl %ebp, (%r9)\n movq %rbp, (%r8,%rcx,8)\n movzwl (%r10,%rbx,2), %ebx\n movl %ebx, 4(%r9)\n shrl $31, %ebp\n leaq (%rdi,%r15,8), %rdi\n xorl $1, %r15d\n leaq (%rdx,%r15,8), %rdx\n xorl %r15d, %ebp\n leaq (%r9,%rbp,8), %r9\n addq $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
}
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)];
u64 QIQN[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); }
QIQN[nq] = (u64)(u32)id | ((u64)(u32)ni << 32); ni += isins; nq += 1 - isins;
}
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
int q = 0;
#if MF_LEAFG >= 3
#if MF_LEAFG >= 4
// ---- l1010a Q4: QUAD-query leaf, nd <= 1 ONLY (one insert-block load, four queries) ----
// Register bound from the row's own corrected inequality, G(nd+2)+4: nd=1 -> 4*3+4 = 16,
// which is EXACTLY the ymm file; nd=2 -> 20 and does NOT fit, so nd=2 falls through to the
// triple rung below (measured 27.3 % of this row's qloop, so the rung is genuinely split).
// nd=0 never occurs for K=4 (measured: PHQ[0] = 0), so `nd <= 1` is exactly nd == 1 here.
// uops per (queries x 8 inserts): triple-nd1 14/24 = 0.5833, quad 18/32 = 0.5625 => -3.57 %.
if (nd <= 1) {
for (; q + 4 <= nq; q += 4) {
const int nA = (int)(QIQN[q] >> 32), nB = (int)(QIQN[q + 1] >> 32), nC = (int)(QIQN[q + 2] >> 32), nD = (int)(QIQN[q + 3] >> 32);
int n = nA < nB ? nA : nB; if (nC < n) n = nC; if (nD < n) n = nD;
__m256i kvA = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q]), kvB = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q + 1]),
kvC = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q + 2]), kvD = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q + 3]);
__m256i v0A = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q]), v0B = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q + 1]),
v0C = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q + 2]), v0D = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q + 3]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(),
accC = _mm256_setzero_si256(), accD = _mm256_setzero_si256();
__builtin_prefetch(&OUT[(u32)QIQN[q]], 1, 1);
__builtin_prefetch(&OUT[(u32)QIQN[q + 1]], 1, 1);
__builtin_prefetch(&OUT[(u32)QIQN[q + 2]], 1, 1);
__builtin_prefetch(&OUT[(u32)QIQN[q + 3]], 1, 1);
int j = 0;
for (; j + 16 <= n; j += 16) {
__m256i IKv = _mm256_load_si256((const __m256i *)(SKA + j));
__m256i mA = _mm256_cmpgt_epi16(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi16(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi16(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi16(kvD, IKv);
if (nd > 0) {
__m256i C0v = _mm256_load_si256((const __m256i *)(SCA + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(v0A, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(v0B, C0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(v0C, C0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi16(v0D, C0v));
}
accA = _mm256_sub_epi16(accA, mA);
accB = _mm256_sub_epi16(accB, mB);
accC = _mm256_sub_epi16(accC, mC);
accD = _mm256_sub_epi16(accD, mD);
}
// per-query tails: each member finishes [j, QN[.]) on its own
{
int jj = j;
__m256i acc = accA;
for (; jj + 16 <= nA; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(SCA + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nA) {
int nn = nA;
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(SCA + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q]] += hsumw256(acc);
}
{
int jj = j;
__m256i acc = accB;
for (; jj + 16 <= nB; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(SCA + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nB) {
int nn = nB;
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(SCA + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q + 1]] += hsumw256(acc);
}
{
int jj = j;
__m256i acc = accC;
for (; jj + 16 <= nC; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvC, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(SCA + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nC) {
int nn = nC;
__m256i msk = _mm256_cmpgt_epi16(kvC, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(SCA + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q + 2]] += hsumw256(acc);
}
{
int jj = j;
__m256i acc = accD;
for (; jj + 16 <= nD; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvD, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0D, _mm256_load_si256((const __m256i *)(SCA + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nD) {
int nn = nD;
__m256i msk = _mm256_cmpgt_epi16(kvD, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0D, _mm256_load_si256((const __m256i *)(SCA + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q + 3]] += hsumw256(acc);
}
}
}
#endif
if (nd <= 2) {
for (; q + 3 <= nq; q += 3) {
const int nA = (int)(QIQN[q] >> 32), nB = (int)(QIQN[q + 1] >> 32), nC = (int)(QIQN[q + 2] >> 32);
int n = nA < nB ? nA : nB; if (nC < n) n = nC;
__m256i kvA = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q]), kvB = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q + 1]), kvC = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q + 2]);
__m256i v0A = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q]), v0B = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q + 1]), v0C = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q + 2]);
__m256i v1A = _mm256_set1_epi16((short)SCB[(KG_BASE + 16) + q]), v1B = _mm256_set1_epi16((short)SCB[(KG_BASE + 16) + q + 1]), v1C = _mm256_set1_epi16((short)SCB[(KG_BASE + 16) + q + 2]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(), accC = _mm256_setzero_si256();
__builtin_prefetch(&OUT[(u32)QIQN[q]], 1, 1);
__builtin_prefetch(&OUT[(u32)QIQN[q + 1]], 1, 1);
__builtin_prefetch(&OUT[(u32)QIQN[q + 2]], 1, 1);
int j = 0;
for (; j + 16 <= n; j += 16) {
__m256i IKv = _mm256_load_si256((const __m256i *)(SKA + j));
__m256i mA = _mm256_cmpgt_epi16(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi16(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi16(kvC, IKv);
if (nd > 0) {
__m256i C0v = _mm256_load_si256((const __m256i *)(SCA + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(v0A, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(v0B, C0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(v0C, C0v));
}
if (nd > 1) {
__m256i C1v = _mm256_load_si256((const __m256i *)(SCB + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi16(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi16(v1B, C1v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi16(v1C, C1v));
}
accA = _mm256_sub_epi16(accA, mA);
accB = _mm256_sub_epi16(accB, mB);
accC = _mm256_sub_epi16(accC, mC);
}
// per-query tails: each member finishes [j, QN[.]) on its own
{
int jj = j;
__m256i acc = accA;
for (; jj + 16 <= nA; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(SCB + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nA) {
int nn = nA;
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(SCB + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q]] += hsumw256(acc);
}
{
int jj = j;
__m256i acc = accB;
for (; jj + 16 <= nB; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(SCB + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nB) {
int nn = nB;
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(SCB + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q + 1]] += hsumw256(acc);
}
{
int jj = j;
__m256i acc = accC;
for (; jj + 16 <= nC; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvC, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1C, _mm256_load_si256((const __m256i *)(SCB + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nC) {
int nn = nC;
__m256i msk = _mm256_cmpgt_epi16(kvC, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1C, _mm256_load_si256((const __m256i *)(SCB + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q + 2]] += hsumw256(acc);
}
}
}
#endif
for (; q + 2 <= nq; q += 2) {
const int n = ((int)(QIQN[q] >> 32) < (int)(QIQN[q + 1] >> 32)) ? (int)(QIQN[q] >> 32) : (int)(QIQN[q + 1] >> 32);
const int nA = (int)(QIQN[q] >> 32), nB = (int)(QIQN[q + 1] >> 32);
__m256i kvA = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q]), kvB = _mm256_set1_epi16((short)SKA[(KG_BASE + 16) + q + 1]);
__builtin_prefetch(&OUT[(u32)QIQN[q]], 1, 1);
__builtin_prefetch(&OUT[(u32)QIQN[q + 1]], 1, 1);
__m256i v0A = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q]), v0B = _mm256_set1_epi16((short)SCA[(KG_BASE + 16) + q + 1]);
__m256i v1A = _mm256_set1_epi16((short)SCB[(KG_BASE + 16) + q]), v1B = _mm256_set1_epi16((short)SCB[(KG_BASE + 16) + q + 1]);
__m256i v2A = _mm256_set1_epi16((short)SCC[(KG_BASE + 16) + q]), v2B = _mm256_set1_epi16((short)SCC[(KG_BASE + 16) + q + 1]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256();
int j = 0;
if (nd >= 3) {
for (; j + 16 <= n; j += 16) {
__m256i IKv = _mm256_load_si256((const __m256i *)(SKA + j));
__m256i C0v = _mm256_load_si256((const __m256i *)(SCA + j));
__m256i C1v = _mm256_load_si256((const __m256i *)(SCB + j));
__m256i C2v = _mm256_load_si256((const __m256i *)(SCC + 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 <= n; j += 16) {
__m256i IKv = _mm256_load_si256((const __m256i *)(SKA + j));
__m256i mA = _mm256_cmpgt_epi16(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi16(kvB, IKv);
if (nd > 0) {
__m256i C0v = _mm256_load_si256((const __m256i *)(SCA + 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 = _mm256_load_si256((const __m256i *)(SCB + 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);
}
}
// per-query tails: each member finishes [j, QN[.]) on its own
{
int jj = j;
__m256i acc = accA;
for (; jj + 16 <= nA; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(SCB + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2A, _mm256_load_si256((const __m256i *)(SCC + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nA) {
int n = nA;
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(SCB + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2A, _mm256_load_si256((const __m256i *)(SCC + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[n - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q]] += hsumw256(acc);
}
{
int jj = j;
__m256i acc = accB;
for (; jj + 16 <= nB; jj += 16) {
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(SCB + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2B, _mm256_load_si256((const __m256i *)(SCC + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nB) {
int n = nB;
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(SKA + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(SCA + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(SCB + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2B, _mm256_load_si256((const __m256i *)(SCC + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[n - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[(u32)QIQN[q + 1]] += hsumw256(acc);
}
}
#endif
for (; q < nq; q++) {
const int n = (int)(QIQN[q] >> 32);
__builtin_prefetch(&OUT[(u32)QIQN[q]], 1, 1);
__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[(u32)QIQN[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
// ONE pass over the dim-0 order: each point's rank is gathered ONCE and both of its
// events are written. Within a dim-0 group the QUERIES occupy [b, b+G) and the
// INSERTs [b+G, b+2G): the group's query-before-insert tie order is what makes an
// equal-dim-0 pair unseparable by the merge's left/right split.
int g = 0;
const int m = 2 * n;
while (g < n) {
int g1 = g;
u32 v = X0[g_ord[g]];
while (g1 < n && X0[g_ord[g1]] == v) g1++;
const int G = g1 - g;
u64 *b = ev + 2 * (long)g;
if (W16) {
for (int t = g; t < g1; t++) { const int id = (int)g_ord[t]; const u32 kk = (u32)X1r[id];
b[t - g] = mk(kk, 1, id); b[G + (t - g)] = mk(kk, 0, id); }
} else {
for (int t = g; t < g1; t++) { const int id = (int)g_ord[t]; const u32 kk = X1u[id];
b[t - g] = mk(kk, 1, id); b[G + (t - g)] = mk(kk, 0, id); }
}
g = g1;
}
ABtop = AB + 2 * n;
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 | 49.789 ms | 11 MB + 972 KB | Accepted | Score: 100 | 显示更多 |