#define KG_MU4 1
#define MF_LEAFLP 128
#define LEAFNB 10240
#define KG_BASE 3840
#include <cstdio>
#ifndef KG_BASE
#define KG_BASE 2560
#endif
#ifndef MF_LEAFLP
#define MF_LEAFLP 32
#endif
#ifndef MF_LEAFG
#define MF_LEAFG 4
#endif
#pragma GCC optimize("O3","rename-registers","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 3072
#endif
#ifndef LEAFNB
#define LEAFNB 4096
#endif
#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 const u32 MV32[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 const u16 MV16[8][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},
};
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 int CT5DEAD;
static inline void ctr_ins(u32 v) {
/* q4a A3: CT5[v>>17] is read only by ctr_qry's `for (j < v>>17) s += CT5[j]` and by
the CT4 base pointer `(v>>17)<<3`. Both are empty/zero while every coordinate is
below 131072, i.e. while the DENSE rank max is below 131072 -- which rank_compress
makes checkable instead of assumed. */
CT0[v]++; CT1[v >> 5]++; CT2[v >> 8]++; CT3[v >> 11]++; CT4[v >> 14]++;
if (!CT5DEAD) CT5[v >> 17]++;
}
static inline void ctr_clr(u32 v) { CT0[v]--; }
static int RSZN;
static u32 CTZ; // a1010a_rv: counter index range = dense rank max (rank values are < CTZ <= 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;
for (u32 j = 0; j < k5; j++) s += CT5[j];
/* lane a1010a_rv R1: the five partial sums were reduced to scalars SEPARATELY (five
horizontal reductions); four of them live in the SAME lane geometry, so they are
accumulated as vectors and reduced ONCE. Same terms, same values, u32 lanes so the
u16 levels cannot overflow. */
const __m256i ix32 = _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7);
const __m128i ix16 = _mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7);
(void)ix32;
__m256i a = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT4 + ((v >> 17) << 3))),
_mm256_loadu_si256((const __m256i *)MV32[(v >> 14) & 7]));
a = _mm256_add_epi32(a, _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT3 + ((v >> 14) << 3))),
_mm256_loadu_si256((const __m256i *)MV32[(v >> 11) & 7])));
(void)ix16;
__m256i b2 = _mm256_cvtepu16_epi32(_mm_and_si128(_mm_loadu_si128((const __m128i *)(CT2 + ((v >> 11) << 3))),
_mm_loadu_si128((const __m128i *)MV16[(v >> 8) & 7])));
__m256i b1 = _mm256_cvtepu16_epi32(_mm_and_si128(_mm_loadu_si128((const __m128i *)(CT1 + ((v >> 8) << 3))),
_mm_loadu_si128((const __m128i *)MV16[(v >> 5) & 7])));
a = _mm256_add_epi32(a, _mm256_add_epi32(b2, b1));
__m256i z = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u))),
_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 a64 = _mm256_add_epi64(_mm256_and_si256(a, _mm256_set1_epi64x(0xFFFFFFFFll)),
_mm256_srli_epi64(a, 32));
__m256i t = _mm256_add_epi64(a64, _mm256_sad_epu8(z, _mm256_setzero_si256()));
__m128i t2 = _mm_add_epi64(_mm256_castsi256_si128(t), _mm256_extracti128_si256(t, 1));
return (u32)((unsigned long long)_mm_cvtsi128_si64(t2)
+ (unsigned long long)_mm_extract_epi64(t2, 1)) + s;
}
static int K;
#ifndef LFAUTO
#define LFAUTO 0
#endif
static int G_LIM;
static int LFSH;
static int LFD;
static const u32 *XD[KG_MAXK];
static u32 *OUT;
static u64 *AB, *ABtop;
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;
}
// 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) {
/* lane a1010a_rv R2: the clear side. Every insert raised CT0[val] with val = the DENSE
rank of dim KT-1, so val < CTZ (rank max) <= RSZN always holds -- the same fact
rank_compress makes checkable. Undoing 5.3e5 inserts by a random RMW per insert is
therefore replaced by ONE contiguous memset of the used extent. */
(void)dirt;
for (int i = 0; i < cx; i++) {
u64 v = cross[i];
u32 val = (u32)(v >> 32);
if ((v >> 31) & 1) ctr_ins(val);
else OUT[eid(v)] += ctr_qry(val);
}
ctr_reset_hi();
memset(CT0, 0, (size_t)(CTZ ? CTZ : (u32)RSZN));
}
// 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);
}
// ===== q4a 16-bit ranked leaf (lane q4a, 2026-10-01) =====
// rank compression: one order-preserving pass per dimension.
static u32 g_rk[KG_MAXK][300005];
static u32 g_rcnt[300005];
static u32 g_mxk, g_mxc0, g_mxc1, g_mxc2;
static int W16 = 0;
#ifdef Q4A_DIAG
static void q4adump(void) __attribute__((destructor));
static void q4adump(void){ fprintf(stderr, "Q4A16 W16=%d mxk=%u mxc0=%u mxc1=%u mxc2=%u\n", W16, g_mxk, g_mxc0, g_mxc1, g_mxc2); }
#endif // 1 => the u16 leaf is admissible for this input
static u16 g_dmax;
static inline u16 XW(u32 x) { return (u16)(x ^ 0x8000u); }
static const u16 VMASKW[16][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}};
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; g_x_st is reused as the counter
static u32 rank_compress(int n, int KT) {
u32 *cnt = g_rcnt;
u32 dmax = 0;
/* q4a A4: dim 0 is consumed only by sort0 (order) and by the top-level grouping
(equality), and the leaf's smallest array index is XD[d+2] >= XD[2]. So dim 0
needs no rank at all -- one of the four compression passes is pure overhead. */
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]]++;
// 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 -- that was this arm's first bug.
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_rk[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; (void)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;
}
int q = 0;
#if MF_LEAFG >= 3
// ---- triple-query leaf: one insert-block load serves three queries ----
// Register accounting, from THIS ROW'S OWN CORRECTED BOUND (README, "THE REGISTER INEQUALITY,
// CORRECTED"): 3 broadcast constants per query x 3 queries + 3 accumulators + 4 data vectors
// = 3*nd + 10 -> nd=2: 16, nd=1: 13, nd=0: 10 => ALL FIT (<= 16 ymm).
// nd = K - (d+2) and at this leaf (K=4) nd in {0,1,2}; the nd>=3 instantiation would need 19
// and is unreachable here, so it keeps the paired path (guarded below, compile-time per ND).
#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)QN[q], nB = (int)QN[q + 1], nC = (int)QN[q + 2], nD = (int)QN[q + 3];
int n = nA < nB ? nA : nB; if (nC < n) n = nC; if (nD < n) n = nD;
__m256i kvA = _mm256_set1_epi32((int)QK[q]), kvB = _mm256_set1_epi32((int)QK[q + 1]),
kvC = _mm256_set1_epi32((int)QK[q + 2]), kvD = _mm256_set1_epi32((int)QK[q + 3]);
__m256i v0A = _mm256_set1_epi32((int)QC0[q]), v0B = _mm256_set1_epi32((int)QC0[q + 1]),
v0C = _mm256_set1_epi32((int)QC0[q + 2]), v0D = _mm256_set1_epi32((int)QC0[q + 3]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(),
accC = _mm256_setzero_si256(), accD = _mm256_setzero_si256();
int j = 0;
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);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi32(kvD, 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));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, C0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, C0v));
}
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
// 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))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nA) {
int nn = 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))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nn - 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))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nB) {
int nn = 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))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nn - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[QID[q + 1]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = accC;
for (; jj + 8 <= nC; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvC, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nC) {
int nn = nC;
__m256i msk = _mm256_cmpgt_epi32(kvC, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nn - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[QID[q + 2]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = accD;
for (; jj + 8 <= nD; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvD, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0D, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nD) {
int nn = nD;
__m256i msk = _mm256_cmpgt_epi32(kvD, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0D, _mm256_load_si256((const __m256i *)(IC0 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nn - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[QID[q + 3]] += hsum256(acc);
}
}
}
#endif
if (nd <= 2) {
for (; q + 3 <= nq; q += 3) {
const int nA = (int)QN[q], nB = (int)QN[q + 1], nC = (int)QN[q + 2];
int n = nA < nB ? nA : nB; if (nC < n) n = nC;
__m256i kvA = _mm256_set1_epi32((int)QK[q]), kvB = _mm256_set1_epi32((int)QK[q + 1]), kvC = _mm256_set1_epi32((int)QK[q + 2]);
__m256i v0A = _mm256_set1_epi32((int)QC0[q]), v0B = _mm256_set1_epi32((int)QC0[q + 1]), v0C = _mm256_set1_epi32((int)QC0[q + 2]);
__m256i v1A = _mm256_set1_epi32((int)QC1[q]), v1B = _mm256_set1_epi32((int)QC1[q + 1]), v1C = _mm256_set1_epi32((int)QC1[q + 2]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(), accC = _mm256_setzero_si256();
int j = 0;
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);
__m256i mC = _mm256_cmpgt_epi32(kvC, 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));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, 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));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v1C, C1v));
}
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
}
// 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))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nA) {
int nn = 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))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nn - 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))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nB) {
int nn = 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))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nn - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[QID[q + 1]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = accC;
for (; jj + 8 <= nC; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvC, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1C, _mm256_load_si256((const __m256i *)(IC1 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nC) {
int nn = nC;
__m256i msk = _mm256_cmpgt_epi32(kvC, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1C, _mm256_load_si256((const __m256i *)(IC1 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nn - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[QID[q + 2]] += hsum256(acc);
}
}
}
#endif
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
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) : LEAFNB;
const int LB = (G_LIM < NB) ? G_LIM : NB; // used-bucket bound
u32 cnt[LEAFNB + 2];
for (int b = 0; b < LB; 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 + 1 < LB; 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; }
// ---- l1010a S1c: single-pass insertion sort over the bucket-sorted prefix ----
{
int lim = (LB >= 2) ? (int)cnt[LB - 2] : 0;
for (int t = 1; t < lim; t++) {
u64 v = out[t];
if (v < out[t - 1]) {
int q = t - 1;
while (q >= 0 && v < out[q]) { out[q + 1] = out[q]; q--; }
out[q + 1] = 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; (void)nd;
__attribute__((aligned(32)))
u16 IK[KG_BASE + 16], IC0[KG_BASE + 16], IC1[KG_BASE + 16], IC2[KG_BASE + 16];
u16 QK[KG_BASE + 16], QC0[KG_BASE + 16], QC1[KG_BASE + 16], QC2[KG_BASE + 16];
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);
#ifdef Q4A_DIAG
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
IK[ni] = XW(k); IC0[ni] = XW(c0); IC1[ni] = XW(c1); IC2[ni] = XW(c2); ni += isins;
QK[nq] = XW(k); QC0[nq] = XW(c0); QC1[nq] = XW(c1); QC2[nq] = XW(c2); QID[nq] = (u32)id; QN[nq] = (u32)ni; nq += 1 - isins;
}
int q = 0;
#if MF_LEAFG >= 3
// ---- triple-query leaf: one insert-block load serves three queries ----
// Register accounting, from THIS ROW'S OWN CORRECTED BOUND (README, "THE REGISTER INEQUALITY,
// CORRECTED"): 3 broadcast constants per query x 3 queries + 3 accumulators + 4 data vectors
// = 3*nd + 10 -> nd=2: 16, nd=1: 13, nd=0: 10 => ALL FIT (<= 16 ymm).
// nd = K - (d+2) and at this leaf (K=4) nd in {0,1,2}; the nd>=3 instantiation would need 19
// and is unreachable here, so it keeps the paired path (guarded below, compile-time per ND).
#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)QN[q], nB = (int)QN[q + 1], nC = (int)QN[q + 2], nD = (int)QN[q + 3];
int n = nA < nB ? nA : nB; if (nC < n) n = nC; if (nD < n) n = nD;
__m256i kvA = _mm256_set1_epi16((short)QK[q]), kvB = _mm256_set1_epi16((short)QK[q + 1]),
kvC = _mm256_set1_epi16((short)QK[q + 2]), kvD = _mm256_set1_epi16((short)QK[q + 3]);
__m256i v0A = _mm256_set1_epi16((short)QC0[q]), v0B = _mm256_set1_epi16((short)QC0[q + 1]),
v0C = _mm256_set1_epi16((short)QC0[q + 2]), v0D = _mm256_set1_epi16((short)QC0[q + 3]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(),
accC = _mm256_setzero_si256(), accD = _mm256_setzero_si256();
int j = 0;
for (; j + 16 <= n; j += 16) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + 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 *)(IC0 + 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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nA) {
int nn = nA;
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nB) {
int nn = nB;
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nC) {
int nn = nC;
__m256i msk = _mm256_cmpgt_epi16(kvC, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0D, _mm256_load_si256((const __m256i *)(IC0 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nD) {
int nn = nD;
__m256i msk = _mm256_cmpgt_epi16(kvD, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0D, _mm256_load_si256((const __m256i *)(IC0 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[q + 3]] += hsumw256(acc);
}
}
}
#endif
if (nd <= 2) {
for (; q + 3 <= nq; q += 3) {
const int nA = (int)QN[q], nB = (int)QN[q + 1], nC = (int)QN[q + 2];
int n = nA < nB ? nA : nB; if (nC < n) n = nC;
__m256i kvA = _mm256_set1_epi16((short)QK[q]), kvB = _mm256_set1_epi16((short)QK[q + 1]), kvC = _mm256_set1_epi16((short)QK[q + 2]);
__m256i v0A = _mm256_set1_epi16((short)QC0[q]), v0B = _mm256_set1_epi16((short)QC0[q + 1]), v0C = _mm256_set1_epi16((short)QC0[q + 2]);
__m256i v1A = _mm256_set1_epi16((short)QC1[q]), v1B = _mm256_set1_epi16((short)QC1[q + 1]), v1C = _mm256_set1_epi16((short)QC1[q + 2]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(), accC = _mm256_setzero_si256();
int j = 0;
for (; j + 16 <= n; j += 16) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + 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 *)(IC0 + 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 *)(IC1 + 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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nA) {
int nn = nA;
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nB) {
int nn = nB;
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1C, _mm256_load_si256((const __m256i *)(IC1 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nC) {
int nn = nC;
__m256i msk = _mm256_cmpgt_epi16(kvC, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0C, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1C, _mm256_load_si256((const __m256i *)(IC1 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[nn - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[q + 2]] += hsumw256(acc);
}
}
}
#endif
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
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_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 <= n; j += 16) {
__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_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 *)(IK + 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 *)(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 = _mm256_load_si256((const __m256i *)(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);
}
}
// 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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nA) {
int n = nA;
__m256i msk = _mm256_cmpgt_epi16(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[n - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[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 *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi16(acc, msk);
}
if (jj < nB) {
int n = nB;
__m256i msk = _mm256_cmpgt_epi16(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASKW[n - jj]));
acc = _mm256_sub_epi16(acc, msk);
}
OUT[QID[q + 1]] += hsumw256(acc);
}
}
#endif
for (; q < nq; q++) {
const int n = (int)QN[q];
__m256i kv = _mm256_set1_epi16((short)QK[q]);
__m256i v0 = _mm256_set1_epi16((short)QC0[q]);
__m256i v1 = _mm256_set1_epi16((short)QC1[q]);
__m256i v2 = _mm256_set1_epi16((short)QC2[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 *)(IK + j))));
break;
case 1:
for (; j + 16 <= n; j += 16) {
__m256i msk = _mm256_cmpgt_epi16(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(IC0 + 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 *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1, _mm256_load_si256((const __m256i *)(IC1 + 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 *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2, _mm256_load_si256((const __m256i *)(IC2 + 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 *)(IK + j)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi16(v2, _mm256_load_si256((const __m256i *)(IC2 + 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) : LEAFNB;
const int LB = (G_LIM < NB) ? G_LIM : NB; // used-bucket bound
u32 cnt[LEAFNB + 2];
for (int b = 0; b < LB; 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 + 1 < LB; 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; }
// ---- l1010a S1c: single-pass insertion sort over the bucket-sorted prefix ----
{
int lim = (LB >= 2) ? (int)cnt[LB - 2] : 0;
for (int t = 1; t < lim; t++) {
u64 v = out[t];
if (v < out[t - 1]) {
int q = t - 1;
while (q >= 0 && v < out[q]) { out[q + 1] = out[q]; q--; }
out[q + 1] = 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
}
// ================= depth-templated solve: D and KT are compile-time =================
template <bool B> struct TB {};
template <int V> struct NCL { static const int v = (V < 0) ? 0 : ((V > 3) ? 3 : V); };
template <int KT, int D> static inline void solve_t(u64 *seq, int m, u64 *out, int needOut);
// ---- terminal level (D == KT-1): every insert-before-query pair qualifies ----
template <int KT, int D>
static void solve_b(u64 *seq, int m, u64 *out, int needOut, TB<true>) {
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
}
}
// ---- general level (D < KT-1) ----
template <int KT, int D>
static void solve_b(u64 *seq, int m, u64 *out, int needOut, TB<false>) {
if (m <= KG_BASE) {
if (W16) base_case_t16<NCL< KT - (D + 2) >::v>(seq, m, D, out, needOut);
else base_case_t<NCL< KT - (D + 2) >::v>(seq, m, D, out, needOut);
return;
}
int h = m >> 1, r = m - h;
u64 *save = ABtop;
u64 *retA = ABtop; ABtop += h;
u64 *retB = ABtop; ABtop += r;
solve_t<KT, D>(seq, h, retA, 1);
solve_t<KT, D>(seq + h, r, retB, 1);
u64 *cross = ABtop; ABtop += m + 1;
const int rekey = (D + 2 <= KT - 1);
const u32 *XN = XD[rekey ? D + 2 : 0];
size_t ia = 0, ib = 0, cx = 0, t = 0;
#ifdef KG_MU4
/* MU4: the balanced merge loop, unrolled four times with the two bound tests HOISTED to
once per four elements. Legal because ia+4 <= h and ib+4 <= r at the top of the block,
and each of the four steps advances exactly one side by one -- so neither side can be
exhausted (and no store can go out of range) inside the block. The scalar loop below
still finishes the remainder, byte-for-byte as before. */
while (ia + 4 <= (size_t)h && ib + 4 <= (size_t)r) {
#define KG_STEP \
{ u64 ea = retA[ia], eb = retB[ib]; \
int takeA = (ea < eb); \
ia += (size_t)takeA; ib += (size_t)(ea >= eb); \
u64 v = takeA ? ea : eb; \
if (needOut) out[t] = v; \
t++; \
u64 c = rekey ? (((u64)XN[eid(v)] << 32) | (v & 0xFFFFFFFFull)) \
: (v & 0xFFFFFFFFull); \
cross[cx] = c; \
cx += (size_t)(takeA == (int)((v >> 31) & 1u)); \
}
KG_STEP KG_STEP KG_STEP KG_STEP
#undef KG_STEP
}
#endif
while (ia < (size_t)h && ib < (size_t)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 += (size_t)takeA; ib += (size_t)(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)XN[eid(v)] << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
}
#else
u64 c = rekey ? (((u64)XN[eid(v)] << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
cross[cx] = c;
cx += keep_test(takeA, role);
#endif
}
while (ia < (size_t)h) {
u64 v = retA[ia++];
if (needOut) out[t] = v;
t++;
if ((v >> 31) & 1) cross[cx++] = rekey ? (((u64)XN[eid(v)] << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
}
while (ib < (size_t)r) {
u64 v = retB[ib++];
if (needOut) out[t] = v;
t++;
if (!((v >> 31) & 1)) cross[cx++] = rekey ? (((u64)XN[eid(v)] << 32) | (v & 0xFFFFFFFFull)) : (v & 0xFFFFFFFFull);
}
if (FUSE_OK && D == KT - 3) {
u32 *dirt = (u32 *)ABtop; ABtop += (cx + 3) / 2 + 1;
ctr_solve(cross, cx, dirt);
} else {
u64 *cres = ABtop; ABtop += cx;
solve_t<KT, D + 1>(cross, cx, cres, 0);
}
ABtop = save;
}
template <int KT, int D>
static inline void solve_t(u64 *seq, int m, u64 *out, int needOut) {
if (m <= 1) return;
solve_b<KT, D>(seq, m, out, needOut, TB<(D >= KT - 1)>());
}
// ====================================================================================
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;
}
template <int KT>
static void run_t(int n, const u32 **x, u32 *out) {
LFSH = 0; LFD = 0;
for (int d = 0; d < KT; d++) XD[d] = x[d];
RSZN = n;
OUT = out;
ABtop = AB;
if (n <= 1) return;
FUSE_OK = 0;
if (KT >= 4 && n > 2 && (u32)n < CT_MAXV) {
const u32 *XL = XD[KT - 1];
memset(CT0, 0, (size_t)n);
int ok = 1;
for (int i = 0; i < n; i++) if (++CT0[XL[i]] > 127) { ok = 0; break; }
memset(CT0, 0, (size_t)n);
FUSE_OK = ok;
}
{
int bit = 0;
while ((n >> bit) >= LEAFNB && bit < 30) bit++;
LFSH = bit; LFD = 0;
}
G_LIM = (int)(((u32)(n - 1)) >> (LFSH + LFD)) + 2;
{
u32 dm = rank_compress(n, KT);
g_dmax = (u16)(dm > 65535u ? 65535u : dm);
W16 = (dm <= 65536u);
/* q4a A5: the bucket key is a compressed rank, so the used-bucket bound is
(D-1)>>LFSH + 2, not (n-1)>>LFSH + 2. With D=63188 against n=100000 the shipped
form walks 3127 buckets where 1977 are reachable -- twice per needOut leaf
(zero + prefix). */
if (dm) G_LIM = (int)(((u32)(dm - 1)) >> (LFSH + LFD)) + 2;
CT5DEAD = (dm <= 131072u);
CTZ = dm;
#ifdef Q4A_DIAG
g_dmax = (u16)dm;
#endif
}
sort0(n);
const u32 *X0 = XD[0], *X1 = XD[1];
u64 *ev = AB; // top-level event array: 2n entries
int m = 0, g = 0;
while (g < n) {
int g1 = g;
u32 v = X0[g_ord[g]];
while (g1 < n && X0[g_ord[g1]] == v) g1++;
/* q4a A6: the two loops gathered the SAME X1[g_ord[t]] and the same g_ord[t] once
each; one pass writes both events from a single gather, preserving the shipped
queries-then-inserts order exactly. */
{
const int gk = g1 - g;
for (int t = g; t < g1; t++) {
u32 id = g_ord[t];
u32 k = X1[id];
ev[m + (t - g)] = mk(k, 1, (int)id);
ev[m + gk + (t - g)] = mk(k, 0, (int)id);
}
m += 2 * gk;
}
g = g1;
}
ABtop = AB + m;
u64 *res = ABtop; ABtop += m;
solve_t<KT, 0>(ev, m, 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_t<4>(n, (const unsigned **)x, out);
}
void count_5d(int n, const unsigned *x[5], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run_t<5>(n, (const unsigned **)x, out);
}
#pragma GCC pop_options
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 34.549 ms | 9 MB + 720 KB | Accepted | Score: 100 | 显示更多 |