// ARMVAR KG_BASE=2048 MF_CTRPF=32
#define KG_BASE 1536
#define KGB0 1536
#define KGB1 1536
#define KGB2 1536
#define MF_CTRPF 0
#ifndef MF_CTRPF
#define MF_CTRPF 0
#endif
#ifndef MF_CTR5
#define MF_CTR5 1
#endif
#ifndef MF_LEAFLP
#define MF_LEAFLP 32
#endif
#ifndef MF_LEAFG
#define MF_LEAFG 4
#endif
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
// duck.ac 1010 (4D) / 1011 (5D) -- "count for every point the number of points that are
// smaller in *every* dimension" (n = 300000, coords in [0,n), strict comparisons).
//
// Algorithm: D&C over dimensions on an event sequence, keys carried inline.
// * 2n events: one query event + one insert event per point, keyed by dim 0
// (a query event is stored *before* an insert event of the same key, which encodes
// the strict-inequality rule).
// * solve(seq,m,d): seq is ordered by dim d and every event carries the coordinate of
// dim d+1 in bits 32..63 of its u64. Split by position, recurse both halves, then
// merge the two halves by the carried key (a plain unsigned `<` on the packed u64,
// so equal keys keep query-before-insert order); the same pass collects the cross
// events (left inserts + right queries), re-keys them to dim d+2 and recurses one
// dimension deeper. Ordering by dim d with that tie rule makes "insert before query"
// exactly "all dims so far strictly smaller", and equal coordinates are never counted.
// * d == K-1: every insert-before-query pair now qualifies, so one O(m) running counter
// finishes the innermost dimension (no log factor).
// * m <= 32: direct O(m^2) counting plus an insertion sort.
//
// Local timings (i5-6500 @ 3.2GHz, best of 3, n = 300000 uniform random):
// 4D 1.57 s 5D 3.85 s
// (Board #1: 508.596 ms / 745.358 ms -- this is still 3x / 5x behind.)
// k-dimensional dominance counting: out[i] = #{j : all dims of j < dims of i}
//
// Event-sequence CDQ with keys carried inline.
// event (u64) = (key << 32) | (ins << 31) | id ins = 1 for an "insert" event
// key = coordinate of dim (level+1): the field used by this level's merge.
// The packing makes a plain unsigned comparison put, for equal keys, the query
// event first -- exactly the strict-inequality rule (equal coords must not count).
//
// solve(seq, m, d): seq is ordered by dim d (tie rule above), elements carry the
// key for dim d+1; count pairs (insert before query) with dims d+1..K-1 strictly
// smaller; return `out` = the same events merged by the carried key (= dim d+1).
// Cross = A-inserts + B-queries re-keyed to dim d+2, recursed at level d+1.
// d == K-1 : every (insert before query) pair qualifies -> running counter, O(m).
#include <cstring>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC target("avx2,popcnt,bmi2")
#include <cstdint>
typedef unsigned u32;
typedef uint64_t u64;
typedef unsigned char u8;
typedef unsigned short u16;
#ifndef KG_COND
#define KG_COND 0
#endif
#ifndef KG_PF
#define KG_PF 0
#endif
#ifndef KG_BASE
#define KG_BASE 4096
#endif
#ifndef LEAFNB
#define LEAFNB 8192
#endif
#define KG_MAXK 5
#define KG_ARENA (32u << 20)
namespace kg {
// ------------ fused innermost-dimension counter (prefix count over [0,n)) -------------
// AVX2. Levels: u8 leaf /32 values, u16 per 32 and per 256, u32 per 2048/16384/131072.
// Every u16 block count must stay <= 32767 (madd uses signed lanes); the guard (max
// multiplicity of a last-dim coordinate <= 127) bounds 32*127=4064 and 256*127=32512.
#define CT_MAXV 1048580u
// q9: FIVE levels covering 19 bits (n <= 524288). Two u8 levels replace four
// u16/u32 ones: the u8 levels are the only ones whose SIMD reduction costs ~9
// instructions for FIVE index bits (su8_32 = one 256-bit load + sad), against
// ~10-11 for THREE bits at u16/u32. Index widths: 5/5/3/3/3.
#define CT_LIM (1u << 19)
/* r2i_rv ARENA */
#define A_CT0 0u
#define A_CT1 (A_CT0 + (CT_LIM + 64u))
#define A_MT8 (A_CT1 + (CT_LIM / 32u + 64u))
#define A_MT16A (A_MT8 + 33u * 32u)
#define A_MT16B (A_MT16A + 33u * 32u)
#define A_CMT3 (A_MT16B + 33u * 32u)
#define A_CMT2 (A_CMT3 + 32u * 32u * 2u)
#define A_CT2 (A_CMT2 + 16u * 16u * 2u)
#define A_CT3 (A_CT2 + (2u * (CT_LIM / 1024u + 16u)))
#define A_CT2P (A_CT3 + (2u * (CT_LIM / 16384u + 32u)))
#define A_P3 (A_CT2P + (2u * (CT_LIM / 1024u + 32u)))
#define A_GSC (A_P3 + 64u)
#define A_END (A_GSC + 4u * (1u << 14))
static u8 ARENA[A_END] __attribute__((aligned(64)));
#define CT0 ((u8 *)(ARENA + A_CT0))
#define CT1 ((u8 *)(ARENA + A_CT1))
#define MT8 ((u8 (*)[32])(ARENA + A_MT8))
#define MT16A ((u16 (*)[16])(ARENA + A_MT16A))
#define MT16B ((u16 (*)[16])(ARENA + A_MT16B))
#define CMT3 ((u16 (*)[32])(ARENA + A_CMT3))
#define CMT2 ((u16 (*)[16])(ARENA + A_CMT2))
#define CT2 ((u16 *)(ARENA + A_CT2))
#define CT3 ((u16 *)(ARENA + A_CT3))
#define CT2P ((u16 *)(ARENA + A_CT2P))
#define P3 ((u16 *)(ARENA + A_P3))
#define GSC ((u32 *)(ARENA + A_GSC))
static int FUSE_OK;
static void mt_init(void){ for(int c=0;c<=32;c++){ for(int j=0;j<32;j++) MT8[c][j]=(u8)((j<c)?0xFF:0);
for(int j=0;j<16;j++){ MT16A[c][j]=(u16)((j<c)?0xFFFF:0); MT16B[c][j]=(u16)((16+j<c)?0xFFFF:0); } }
for(int k=0;k<32;k++) for(int j=0;j<32;j++) CMT3[k][j]=(u16)((j>k)?1:0);
for(int k=0;k<16;k++) for(int j=0;j<16;j++) CMT2[k][j]=(u16)((j>k)?1:0); }
static inline u32 su8_32(const u8 *p, int cnt) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
__m256i m = _mm256_cmpgt_epi8(_mm256_set1_epi8((char)cnt),
_mm256_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31));
__m256i x = _mm256_and_si256(v, m);
__m256i s = _mm256_sad_epu8(x, _mm256_setzero_si256());
__m128i t = _mm_add_epi64(_mm256_castsi256_si128(s), _mm256_extracti128_si256(s, 1));
return (u32)((unsigned long long)_mm_cvtsi128_si64(t) + (unsigned long long)_mm_extract_epi64(t, 1));
}
static inline u32 su16_8(const u16 *p, int cnt) {
__m128i v = _mm_loadu_si128((const __m128i *)p);
__m128i x = _mm_and_si128(v, _mm_cmpgt_epi16(_mm_set1_epi16((short)cnt),
_mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7)));
__m128i s = _mm_madd_epi16(x, _mm_set1_epi16(1));
s = _mm_add_epi32(s, _mm_shuffle_epi32(s, _MM_SHUFFLE(1, 0, 3, 2)));
s = _mm_add_epi32(s, _mm_shuffle_epi32(s, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(s);
}
static inline u32 su32_8(const u32 *p, int cnt) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
__m256i x = _mm256_and_si256(v, _mm256_cmpgt_epi32(_mm256_set1_epi32(cnt),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7)));
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
static inline void ctr_ins(u32 v) {
CT0[v]++; CT1[v >> 5]++; CT2[v >> 10]++; CT3[v >> 14]++;
}
/* r2i_rv fast insert: CT2 -> CT2P (prefix suffix-update, one 32 B RMW). */
static inline void ctr_ins_fast(u32 v) {
CT0[v]++; CT1[v >> 5]++;
const u32 w = v >> 10;
__m256i t = _mm256_load_si256((const __m256i *)(CT2P + (w & ~15u)));
t = _mm256_add_epi16(t, _mm256_load_si256((const __m256i *)CMT2[w & 15u]));
_mm256_store_si256((__m256i *)(CT2P + (w & ~15u)), t);
}
static inline void ctr_clr(u32 v) { CT0[v]--; }
static int RSZN;
static inline void ctr_reset_hi(void) {
/* q9: extents = the highest index the query path can READ, +1 */
memset(CT1, 0, sizeof(u8) * ((size_t)((unsigned)RSZN >> 5) + 64));
memset(CT2, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 10) + 16));
memset(CT3, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 14) + 32));
memset(CT2P, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 10) + 32));
memset(P3, 0, 64);
}
#define QIDX32 _mm256_setr_epi32(0,1,2,3,4,5,6,7)
#define QIDX16 _mm256_setr_epi16(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
#define QIDX16A _mm256_setr_epi16(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
#define QIDX16B _mm256_setr_epi16(16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31)
#define QIDX8 _mm256_setr_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15, \
16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31)
/* v1: ONE __m256i accumulator across all five levels, ONE horizontal reduce.
Lane assignment is irrelevant because the final reduce sums every lane.
The u16 levels load 16 entries; cnt <= 7 masks the upper 8 off (the array is
padded so the extra 16 bytes are always inside the object). */
static inline u32 ctr_qry(u32 v) {
const __m256i one16 = _mm256_set1_epi16(1);
const __m256i one8 = _mm256_set1_epi8(1); (void)one8;
__m256i acc = _mm256_setzero_si256();
{
/* L3: CT3, u16, 32 entries = 64 B window, base identically 0 (v < 2^19).
mask table index (c3<<5) = ((v>>9) & 0x3e0). */
const u32 t3 = (v >> 9) & 0x3e0u;
__m256i mk = _mm256_load_si256((const __m256i *)MT16A[t3 >> 5]);
__m256i x0 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)CT3), mk);
__m256i x1 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT3 + 16)),
_mm256_load_si256((const __m256i *)MT16B[t3 >> 5]));
acc = _mm256_add_epi32(acc, _mm256_madd_epi16(x0, one16));
acc = _mm256_add_epi32(acc, _mm256_madd_epi16(x1, one16));
}
{
const u32 w2 = v >> 10;
__m256i m = _mm256_load_si256((const __m256i *)MT16A[w2 & 15u]);
__m256i x = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT2 + (w2 & ~15u))), m);
acc = _mm256_add_epi32(acc, _mm256_madd_epi16(x, one16));
}
{
/* r2h_rv: the two u8 levels go through vpsadbw. sad is 5 cyc from the AND
(vs maddubs 5 + madd 5 = 10), and it zero-extends to u64 lanes, so the
two u64 halves of each lane are zero and a single vpaddd folds them into
the u32 accumulator exactly. */
const u32 w1 = v >> 5;
__m256i m = _mm256_load_si256((const __m256i *)MT8[w1 & 31u]);
__m256i x1 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT1 + (w1 & ~31u))), m);
__m256i m0 = _mm256_load_si256((const __m256i *)MT8[v & 31u]);
__m256i x0 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u))), m0);
__m256i s = _mm256_sad_epu8(_mm256_add_epi8(x1, x0), _mm256_setzero_si256());
acc = _mm256_add_epi32(acc, s);
}
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(acc), _mm256_extracti128_si256(acc, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
/* r2i_rv fast query. c2 = CT2P[a>>10], c3 = P3[a>>14] are read by the caller (the
P3 read is preceded by materialising the two P3 registers to memory). The two u8
levels stay a masked 32 B window summed with vpsadbw; the two prefix levels enter as
ONE zero-extended scalar. Every contribution therefore sits in a u64 lane whose
upper 32 bits are 0, so the whole horizontal reduce is vpaddq + movq/pextrq/add. */
static inline u32 ctr_qry_fast(u32 v, u32 c2, u32 c3) {
const u32 w1 = v >> 5;
__m256i m = _mm256_load_si256((const __m256i *)MT8[w1 & 31u]);
__m256i x1 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT1 + (w1 & ~31u))), m);
__m256i m0 = _mm256_load_si256((const __m256i *)MT8[v & 31u]);
__m256i x0 = _mm256_and_si256(
_mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u))), m0);
__m256i acc = _mm256_sad_epu8(_mm256_add_epi8(x1, x0), _mm256_setzero_si256());
__m128i t = _mm_add_epi64(_mm256_castsi256_si128(acc), _mm256_extracti128_si256(acc, 1));
unsigned long long s64 = (unsigned long long)_mm_cvtsi128_si64(t)
+ (unsigned long long)_mm_extract_epi64(t, 1);
return (u32)((u32)s64 + c2 + c3);
}
static int K;
#ifndef LFAUTO
#define LFAUTO 0
#endif
static int G_LIM;
static int LFSH;
static int LFD;
static const u32 *XD[KG_MAXK];
static u32 *OUT;
static u64 *AB, *ABtop;
// ===== v2 PACKED EVENT =====================================================
// v = [KA:20][TY:1][KB:20][ID:20] (bits 43..62, 42, 22..41, 0..19)
// KA = the compare key for THIS level (dim D+1), KB = the key for the NEXT
// level (dim D+2), carried so the cross needs no re-key gather.
// A plain u64 compare gives (KA, TY, KB, ID) -- KA primary, TY breaks equal
// KA with query(TY=0) before insert(TY=1) = the strict-inequality rule.
// Requires id < 2^20 and every coordinate < 2^20 (true: n = 1e6 < 2^20).
#define E_ID(v) ((u32)((v) & 0xFFFFFu))
#define E_KB(v) ((u32)(((v) >> 22) & 0xFFFFFu))
#define E_TY(v) ((u32)(((v) >> 42) & 1u))
#define E_KA(v) ((u32)((v) >> 43))
static inline u64 mk2(u32 ka, u32 ty, u32 kb, u32 id) {
return ((u64)ka << 43) | ((u64)ty << 42) | ((u64)kb << 22) | (u64)id;
}
static inline int eid(u64 v) { return (int)(v & 0x7FFFFFFFu); }
static inline int lfbucket(u64 v) { return (int)(((u32)(v >> 32)) >> LFSH); }
static inline int eq_(u64 v) { return (int)((v >> 31) & 1); }
static inline int keep_test(int takeA, int role) { return (takeA ^ role) ^ 1; }
static inline u64 mk(u32 key, int q, int id) {
return ((u64)key << 32) | ((u64)(1 - q) << 31) | (u64)id;
}
// resolve a fusion-level cross: events carry the last dim (K-1) coordinate in bits 32..63
// and are ordered by dim K-2; a left-insert before a right-query with a strictly smaller
// last-dim value is a hit.
/* ---------------- r2i_rv: the original (fallback) engine, byte-for-byte ----------- */
static void ctr_solve_slow(u64 *cross, int cx, u32 *dirt) {
u32 *dirtp = dirt;
const u64 *ce = cross + cx;
const u64 *cp = cross;
for (; cp + 4 <= ce; cp += 4) {
{ u64 v=cp[0]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
{ u64 v=cp[1]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
{ u64 v=cp[2]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
{ u64 v=cp[3]; u32 a=E_KB(v); if (E_TY(v)) { ctr_ins(a); *dirtp++ = a; } else OUT[E_ID(v)] += ctr_qry(a); }
}
for (; cp < ce; cp++) {
u64 v = *cp;
u32 val = E_KB(v);
if (E_TY(v)) { ctr_ins(val); *dirtp++ = val; }
else OUT[E_ID(v)] += ctr_qry(val);
}
int nd = (int)(dirtp - dirt);
for (int i = 0; i < nd; i++) {
#if MF_CTRPF > 0
if (i + MF_CTRPF < nd) {
u32 vf = dirt[i + MF_CTRPF];
__builtin_prefetch(&CT0[vf], 1, 1);
__builtin_prefetch(&CT1[vf >> 5], 1, 1);
__builtin_prefetch(&CT2[vf >> 8], 1, 1);
}
#endif
ctr_clr(dirt[i]);
}
ctr_reset_hi();
}
/* ---------------- r2i_rv: the prefix-top engine ----------------------------------- */
#define P3_STEP(a) do { const u32 k_ = (a) >> 14; \
p3v0 = _mm256_add_epi16(p3v0, _mm256_load_si256((const __m256i *)CMT3[k_])); \
p3v1 = _mm256_add_epi16(p3v1, _mm256_load_si256((const __m256i *)(CMT3[k_] + 16))); \
} while (0)
#define P3_MAT() do { \
_mm256_store_si256((__m256i *)P3, p3v0); \
_mm256_store_si256((__m256i *)(P3 + 16), p3v1); } while (0)
#define EV_F(v) do { const u64 v_ = (v); const u32 a_ = E_KB(v_); \
if (E_TY(v_)) { ctr_ins_fast(a_); P3_STEP(a_); *dirtp++ = a_; } \
else { P3_MAT(); OUT[E_ID(v_)] += ctr_qry_fast(a_, CT2P[a_ >> 10], P3[a_ >> 14]); } } while (0)
static u64 CROSSB[65536 + 8];
static void ctr_solve_split2(u64 *ca, int nca, u64 *cb, int ncb, u32 *dirt) {
u32 *dirtp = dirt;
__m256i p3v0 = _mm256_setzero_si256(), p3v1 = _mm256_setzero_si256();
const u64 *ce = ca + nca, *cp = ca;
for (; cp + 4 <= ce; cp += 4) {
EV_F(cp[0]); EV_F(cp[1]); EV_F(cp[2]); EV_F(cp[3]);
}
for (; cp < ce; cp++) EV_F(*cp);
ce = cb + ncb; cp = cb;
for (; cp + 4 <= ce; cp += 4) {
EV_F(cp[0]); EV_F(cp[1]); EV_F(cp[2]); EV_F(cp[3]);
}
for (; cp < ce; cp++) EV_F(*cp);
int nd = (int)(dirtp - dirt);
for (int i = 0; i < nd; i++) ctr_clr(dirt[i]);
ctr_reset_hi();
(void)p3v0; (void)p3v1;
}
static void ctr_solve_fast(u64 *cross, int cx, u32 *dirt) {
u32 *dirtp = dirt;
const u64 *ce = cross + cx;
const u64 *cp = cross;
__m256i p3v0 = _mm256_setzero_si256(), p3v1 = _mm256_setzero_si256();
for (; cp + 4 <= ce; cp += 4) {
EV_F(cp[0]); EV_F(cp[1]); EV_F(cp[2]); EV_F(cp[3]);
}
for (; cp < ce; cp++) EV_F(*cp);
int nd = (int)(dirtp - dirt);
for (int i = 0; i < nd; i++) ctr_clr(dirt[i]);
ctr_reset_hi();
(void)p3v0; (void)p3v1;
}
/* P3's entries are u16 prefixes of the WHOLE call's insert count, so the fast path is
used exactly when no call can push them past 65535: inserts <= cx. */
static void ctr_solve(u64 *cross, int cx, u32 *dirt) {
if (cx <= 65535) ctr_solve_fast(cross, cx, dirt);
else ctr_solve_slow(cross, cx, dirt);
}
/* lane r2k_rv: 8-wide inclusive prefix sum over u16, replacing
`for (b = 0; b < n; b++) c[b + 1] += c[b];` (an inclusive prefix over n+1 entries).
The scalar form is one loop-carried load->add->store per entry; this one is one
paddw + one pshufb per 8 entries. Counts here are bounded by m (<= KGB) or by
cx (<= 65535), so the u16 arithmetic cannot wrap. */
static inline void pfx16(u16 *c, int n) {
const __m128i lastmask = _mm_set1_epi16(0x0F0E);
__m128i carry = _mm_setzero_si128();
int b = 0;
for (; b + 8 <= n; b += 8) {
__m128i v = _mm_loadu_si128((const __m128i *)(c + b));
__m128i x = _mm_add_epi16(v, _mm_slli_si128(v, 2));
x = _mm_add_epi16(x, _mm_slli_si128(x, 4));
x = _mm_add_epi16(x, _mm_slli_si128(x, 8));
x = _mm_add_epi16(x, carry);
_mm_storeu_si128((__m128i *)(c + b), x);
carry = _mm_shuffle_epi8(x, lastmask);
}
if (b < n) {
u16 cc = (u16)_mm_extract_epi16(carry, 0);
for (; b < n; b++) { cc = (u16)(cc + c[b]); c[b] = cc; }
}
}
// small base: compact insert list + SIMD comparison, then insertion sort by the key
static const u32 VMASK[8][8] = {
{0,0,0,0,0,0,0,0},
{~0u,0,0,0,0,0,0,0},
{~0u,~0u,0,0,0,0,0,0},
{~0u,~0u,~0u,0,0,0,0,0},
{~0u,~0u,~0u,~0u,0,0,0,0},
{~0u,~0u,~0u,~0u,~0u,0,0,0},
{~0u,~0u,~0u,~0u,~0u,~0u,0,0},
{~0u,~0u,~0u,~0u,~0u,~0u,~0u,0}};
static inline u32 hsum256(__m256i v) {
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(v), _mm256_extracti128_si256(v, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
#define L4Q_TAIL(acc, kv, v0, v1, QNx, j, qid) do { \
int jj = (j); __m256i a_ = (acc); \
const int nx_ = (int)(QNx); \
for (; jj + 8 <= nx_; jj += 8) { \
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + jj))); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + jj)))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + jj)))); \
a_ = _mm256_sub_epi32(a_, msk); \
} \
if (jj < nx_) { \
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + jj))); \
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + jj)))); \
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + jj)))); \
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[nx_ - jj])); \
a_ = _mm256_sub_epi32(a_, msk); \
} \
OUT[(qid)] += hsum256(a_); \
} while (0)
template<int ND>
static void base_case_t(u64 *seq, int m, int d, u64 *out, int needOut) {
const int nd = ND;
__attribute__((aligned(32)))
u32 IK[KG_BASE + 8], IC0[KG_BASE + 8], IC1[KG_BASE + 8], IC2[KG_BASE + 8];
// r2g_rv: the query-side u32 pairs are PACKED into u64 arrays. The partition loop is
// STORE-UOP bound (1 store/cycle), so 4 u32 stores -> 2 u64 stores per event. The layout is
// chosen so the sweep's reads are still single 4-byte loads (vpbroadcastd from memory).
u64 QKC[KG_BASE + 8]; // hi32 = k (dim D+1) lo32 = c0 (dim D+2)
u64 QIN[KG_BASE + 8]; // hi32 = ni (insert prefix) lo32 = id
u64 QC12[KG_BASE + 8]; // hi32 = c1 (dim D+3), lo32 = c2 (dim D+4)
const u32 *X1 = nd > 1 ? XD[d + 3] : 0, *X2 = nd > 2 ? XD[d + 4] : 0;
int ni = 0, nq = 0;
for (int t = 0; t < m; t++) {
#if MF_LEAFLP > 0
if (t + MF_LEAFLP < m) {
int idf = (int)E_ID(seq[t + MF_LEAFLP]);
if (nd > 1) __builtin_prefetch(&X1[idf], 0, 3);
if (nd > 2) __builtin_prefetch(&X2[idf], 0, 3);
}
#endif
u64 v = seq[t];
int id = (int)E_ID(v);
u32 k = E_KA(v);
u32 c0 = 0, c1 = 0, c2 = 0;
if (nd > 0) c0 = E_KB(v); // dim d+2 is IN the event -- no gather
if (nd > 1) c1 = X1[id];
if (nd > 2) c2 = X2[id];
int isins = (int)E_TY(v);
IK[ni] = k; IC0[ni] = c0; IC1[ni] = c1; IC2[ni] = c2; ni += isins;
QKC[nq] = ((u64)k << 32) | (u64)c0; QC12[nq] = ((u64)c1 << 32) | (u64)c2;
QIN[nq] = ((u64)(u32)ni << 32) | (u64)(u32)id; nq += 1 - isins;
}
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
int q = 0;
#if MF_LEAFG >= 4
// ================= FOUR-QUERY MERGED LEAF (lane m188) =================
// Gated on `nd <= 1` BY REGISTER ARITHMETIC, not by taste: at nd=2 the live set is
// 3 + 4*4 = 19 ymm against 16. See rig/mk4q.py's header for the full table.
if (nd <= 1) {
for (; q + 4 <= nq; q += 4) {
const int nA = (int)(u32)(QIN[q] >> 32), nB = (int)(u32)(QIN[q + 1] >> 32),
nC = (int)(u32)(QIN[q + 2] >> 32), nD = (int)(u32)(QIN[q + 3] >> 32);
__m256i kvA = _mm256_set1_epi32((int)(QKC[q] >> 32)), kvB = _mm256_set1_epi32((int)(QKC[q + 1] >> 32)),
kvC = _mm256_set1_epi32((int)(QKC[q + 2] >> 32)), kvD = _mm256_set1_epi32((int)(QKC[q + 3] >> 32));
__m256i v0A = _mm256_set1_epi32((int)QKC[q]), v0B = _mm256_set1_epi32((int)QKC[q + 1]),
v0C = _mm256_set1_epi32((int)QKC[q + 2]),v0D = _mm256_set1_epi32((int)QKC[q + 3]);
#if (1) /* nd==2 never reaches here, so v1* is only emitted when nd can be 1 -- see below */
__m256i v1A = _mm256_set1_epi32((int)(QC12[q] >> 32)), v1B = _mm256_set1_epi32((int)(QC12[q + 1] >> 32)),
v1C = _mm256_set1_epi32((int)(QC12[q + 2] >> 32)),v1D = _mm256_set1_epi32((int)(QC12[q + 3] >> 32));
#endif
__builtin_prefetch(&OUT[(u32)QIN[q + 8]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 10]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 11]], 0, 1);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(),
accC = _mm256_setzero_si256(), accD = _mm256_setzero_si256();
int j = 0, jA, jB, jC, jD;
for (; j + 8 <= nA; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mA = _mm256_cmpgt_epi32(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v0A, v0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
jA = j;
for (; j + 8 <= nB; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
jB = j;
for (; j + 8 <= nC; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
jC = j;
for (; j + 8 <= nD; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mD = _mm256_cmpgt_epi32(kvD, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v0D, v0v));
accD = _mm256_sub_epi32(accD, mD);
}
jD = j;
L4Q_TAIL(accA, kvA, v0A, v1A, nA, jA, (u32)QIN[q]);
L4Q_TAIL(accB, kvB, v0B, v1B, nB, jB, (u32)QIN[q + 1]);
L4Q_TAIL(accC, kvC, v0C, v1C, nC, jC, (u32)QIN[q + 2]);
L4Q_TAIL(accD, kvD, v0D, v1D, nD, jD, (u32)QIN[q + 3]);
}
}
#endif
#if MF_LEAFG >= 4
if (nd == 2) { // q5b/tri: three queries share one insert-block load (nd==3 needs a 4th dim: IC2)
for (; q + 3 <= nq; q += 3) {
const int nA = (int)(u32)(QIN[q] >> 32), nB = (int)(u32)(QIN[q + 1] >> 32), nC = (int)(u32)(QIN[q + 2] >> 32);
__m256i kvA = _mm256_set1_epi32((int)(QKC[q] >> 32)), kvB = _mm256_set1_epi32((int)(QKC[q + 1] >> 32)),
kvC = _mm256_set1_epi32((int)(QKC[q + 2] >> 32));
__m256i v0A = _mm256_set1_epi32((int)QKC[q]), v0B = _mm256_set1_epi32((int)QKC[q + 1]),
v0C = _mm256_set1_epi32((int)QKC[q + 2]);
__m256i v1A = _mm256_set1_epi32((int)(QC12[q] >> 32)), v1B = _mm256_set1_epi32((int)(QC12[q + 1] >> 32)),
v1C = _mm256_set1_epi32((int)(QC12[q + 2] >> 32));
__builtin_prefetch(&OUT[(u32)QIN[q + 8]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 10]], 0, 1);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256(),
accC = _mm256_setzero_si256();
int j = 0, jA, jB, jC;
for (; j + 8 <= nA; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mA = _mm256_cmpgt_epi32(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v0A, v0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
__m256i v1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, v1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, v1v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v1C, v1v));
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
}
jA = j;
for (; j + 8 <= nB; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, v0v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
__m256i v1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, v1v));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v1C, v1v));
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
}
jB = j;
for (; j + 8 <= nC; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mC = _mm256_cmpgt_epi32(kvC, IKv);
__m256i v0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v0C, v0v));
__m256i v1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mC = _mm256_and_si256(mC, _mm256_cmpgt_epi32(v1C, v1v));
accC = _mm256_sub_epi32(accC, mC);
}
jC = j;
L4Q_TAIL(accA, kvA, v0A, v1A, nA, jA, (u32)QIN[q]);
L4Q_TAIL(accB, kvB, v0B, v1B, nB, jB, (u32)QIN[q + 1]);
L4Q_TAIL(accC, kvC, v0C, v1C, nC, jC, (u32)QIN[q + 2]);
}
}
#endif
for (; q + 2 <= nq; q += 2) {
const int n = ((int)(u32)(QIN[q] >> 32) < (int)(u32)(QIN[q + 1] >> 32)) ? (int)(u32)(QIN[q] >> 32) : (int)(u32)(QIN[q + 1] >> 32);
const int nA = (int)(u32)(QIN[q] >> 32), nB = (int)(u32)(QIN[q + 1] >> 32);
__m256i kvA = _mm256_set1_epi32((int)(QKC[q] >> 32)), kvB = _mm256_set1_epi32((int)(QKC[q + 1] >> 32));
__m256i v0A = _mm256_set1_epi32((int)QKC[q]), v0B = _mm256_set1_epi32((int)QKC[q + 1]);
__m256i v1A = _mm256_set1_epi32((int)(QC12[q] >> 32)), v1B = _mm256_set1_epi32((int)(QC12[q + 1] >> 32));
__m256i v2A = _mm256_set1_epi32((int)QC12[q]), v2B = _mm256_set1_epi32((int)QC12[q + 1]);
__builtin_prefetch(&OUT[(u32)QIN[q + 8]], 0, 1); __builtin_prefetch(&OUT[(u32)QIN[q + 9]], 0, 1);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256();
int j = 0;
if (nd >= 3) {
for (; j + 8 <= n; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i C0v = _mm256_load_si256((const __m256i *)(IC0 + j));
__m256i C1v = _mm256_load_si256((const __m256i *)(IC1 + j));
__m256i C2v = _mm256_load_si256((const __m256i *)(IC2 + j));
__m256i mA = _mm256_and_si256(_mm256_cmpgt_epi32(kvA, IKv), _mm256_cmpgt_epi32(v0A, C0v));
__m256i mB = _mm256_and_si256(_mm256_cmpgt_epi32(kvB, IKv), _mm256_cmpgt_epi32(v0B, C0v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, C1v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v2A, C2v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v2B, C2v));
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
}
} else {
for (; j + 8 <= n; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mA = _mm256_cmpgt_epi32(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
if (nd > 0) {
__m256i C0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v0A, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, C0v));
}
if (nd > 1) {
__m256i C1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, C1v));
}
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
}
}
// per-query tails: each member finishes [j, QN[.]) on its own
{
int jj = j;
__m256i acc = accA;
for (; jj + 8 <= nA; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nA) {
int n = nA;
__m256i msk = _mm256_cmpgt_epi32(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[(u32)QIN[q]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = accB;
for (; jj + 8 <= nB; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nB) {
int n = nB;
__m256i msk = _mm256_cmpgt_epi32(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[(u32)QIN[q + 1]] += hsum256(acc);
}
}
#endif
for (; q < nq; q++) {
const int n = (int)(u32)(QIN[q] >> 32);
__m256i kv = _mm256_set1_epi32((int)(QKC[q] >> 32));
__m256i v0 = _mm256_set1_epi32((int)QKC[q]);
__m256i v1 = _mm256_set1_epi32((int)(QC12[q] >> 32));
__m256i v2 = _mm256_set1_epi32((int)QC12[q]);
__m256i acc = _mm256_setzero_si256();
int j = 0;
switch (nd) {
case 0:
for (; j + 8 <= n; j += 8)
acc = _mm256_sub_epi32(acc, _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j))));
break;
case 1:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
case 2:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
default:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC2 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
}
if (j < n) { // partial block, masked by a table lookup
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC2 + j))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - j]));
acc = _mm256_sub_epi32(acc, msk);
}
u32 c = hsum256(acc);
OUT[(u32)QIN[q]] += c;
}
if (!needOut) return;
#if LEAFNB
{
const int NB = (LFAUTO && m <= LEAFNB) ? (m < 8 ? 8 : m) : LEAFNB;
const int LB = (G_LIM < NB) ? G_LIM : NB; // used-bucket bound
(void)LB; (void)NB;
{ // q5b: stable 2-pass LSD radix on the 21-bit (KA, TY) key (bits 42..62)
u64 *tmp = ABtop; ABtop += m;
u16 c0_[1024 + 1], c1_[2048 + 1];
for (int b = 0; b <= 1024; b++) c0_[b] = 0;
for (int t = 0; t < m; t++) c0_[((u32)(seq[t] >> 42) & 1023u) + 1]++;
pfx16(c0_, 1024 + 1);
for (int t = 0; t < m; t++) { u64 v = seq[t]; tmp[c0_[(u32)(v >> 42) & 1023u]++] = v; }
for (int b = 0; b <= 2048; b++) c1_[b] = 0;
for (int t = 0; t < m; t++) c1_[((u32)(tmp[t] >> 52) & 2047u) + 1]++;
pfx16(c1_, 2048 + 1);
for (int t = 0; t < m; t++) { u64 v = tmp[t]; out[c1_[(u32)(v >> 52) & 2047u]++] = v; }
ABtop -= m;
}
}
#else
for (int t = 1; t < m; t++) {
u64 v = seq[t];
int q = t - 1;
while (q >= 0 && v < seq[q]) { seq[q + 1] = seq[q]; q--; }
seq[q + 1] = v;
}
memcpy(out, seq, sizeof(u64) * (size_t)m);
#endif
}
static void base_case(u64 *seq, int m, int d, u64 *out, int needOut) {
switch (K - (d + 2)) {
case 0: base_case_t<0>(seq, m, d, out, needOut); break;
case 1: base_case_t<1>(seq, m, d, out, needOut); break;
case 2: base_case_t<2>(seq, m, d, out, needOut); break;
default: base_case_t<3>(seq, m, d, out, needOut); break;
}
}
// ---- ND specialisation dispatcher (the fleet's `DN` recipe) ----
static void base_case(u64 *seq, int m, int d, u64 *out, int needOut);
// ================= depth-templated solve: D and KT are compile-time =================
template <bool B, bool NO = true> struct TB {};
template <int V> struct NCL { static const int v = (V < 0) ? 0 : ((V > 3) ? 3 : V); };
template <int KT, int D, bool NO> static inline void solve_t(u64 *seq, int m, u64 *out);
// ---- terminal level (D == KT-1): every insert-before-query pair qualifies ----
template <int KT, int D, bool NO>
static void solve_b(u64 *seq, int m, u64 *out, TB<true, NO>) {
int cnt = 0;
for (int t = 0; t < m; t++) {
u64 v = seq[t];
if (E_TY(v)) cnt++; // insert
else OUT[E_ID(v)] += (u32)cnt; // query
}
}
// ---- row 1011 lane r2: the merge split into TWO independent merge-path chains --------
// One step of the reference two-pointer merge, verbatim from the bound bytes, so every
// chain emits exactly the reference stream for its own t-range.
#define KG_MSTEP(IA, IB, CX, CRS, OIDX) do { \
if (GAT) { __builtin_prefetch(&XN[E_ID(retA[(IA)+2])], 0, 3); __builtin_prefetch(&XN[E_ID(retB[(IB)+2])], 0, 3); } \
u64 ea = retA[(IA)], eb = retB[(IB)]; \
int takeA = (ea < eb); \
(IA) += takeA; (IB) += takeA ^ 1; \
u64 v = takeA ? ea : eb; \
if (NO) out[(OIDX)] = v; \
int role = (int)E_TY(v); \
u64 c = GAT ? ((((v << 21) & (0xFFFFFull << 43)) | (v & (1ull << 42)) | (v & 0xFFFFFull)) | ((u64)XN[E_ID(v)] << 22)) \
: (PROM ? (((v << 21) & (0xFFFFFull << 43)) | (v & (1ull << 42)) | (v & 0xFFFFFull)) : v); \
(CRS)[(CX)] = (KT - 3 == D) ? v : c; \
(CX) += keep_test(takeA, role); \
} while (0)
// (I,J) = the state the merge passes through when it has consumed every element whose
// key (v>>42) is <= k, for the largest such k with I+J <= target. SOUND on inputs that
// are only key-sorted (see the header). I = #A with key <= k, J = #B with key <= k.
static inline void kg_kbound(const u64 *retA, int h, const u64 *retB, int r, int d,
int *pI, int *pJ) {
int lo = (d - r > 0) ? (d - r) : 0;
int hi = (d < h) ? d : h;
while (lo < hi) {
int mid = (lo + hi + 1) >> 1;
if (mid == 0 || (retA[mid - 1] >> 42) < (retB[d - mid] >> 42)) lo = mid; else hi = mid - 1;
}
if (lo == 0) { *pI = 0; *pJ = 0; return; }
u64 kk = retA[lo - 1] >> 42;
if (d > lo) { u64 kb = retB[d - lo - 1] >> 42; if (kb > kk) kk = kb; }
{ int a = lo; while (a < h && (retA[a] >> 42) <= kk) a++; *pI = a; }
{ int a = (d - lo > 0) ? (d - lo) : 0; while (a < r && (retB[a] >> 42) <= kk) a++; *pJ = a; }
}
// ---- general level (D < KT-1) ----
template <int KT, int D, bool NO>
static void solve_b(u64 *seq, int m, u64 *out, TB<false, NO>) {
const int KGB = (D == 0) ? KGB0 : ((D == 1) ? KGB1 : KGB2);
if (m <= KGB) {
base_case_t<NCL< KT - (D + 2) >::v>(seq, m, D, out, NO ? 1 : 0);
return;
}
int h = m >> 1, r = m - h;
u64 *save = ABtop;
u64 *retA = ABtop; ABtop += h + 1 + 64;
u64 *retB = ABtop; ABtop += r + 1 + 64;
solve_t<KT, D, true>(seq, h, retA);
solve_t<KT, D, true>(seq + h, r, retB);
retA[h] = ~0ull; retB[r] = ~0ull;
u64 *cross = ABtop; ABtop += m + 1;
// v2: PROMOTE moves the carried KB up into the KA slot (a new level's key);
// GATHER fetches the next KB (dim D+3) -- the ONLY remaining per-level gather,
// and it exists only at D = 0 for KT = 4.
const int PROM = (D + 2 <= KT - 1);
const int GAT = (D + 3 <= KT - 1);
const u32 *XN = XD[(D + 3 <= KG_MAXK - 1) ? (D + 3) : 0];
long ia = 0, ib = 0, cx = 0, cx1 = 0;
/* lane r2k_rv: chain 2 writes its OWN segment and the counter body consumes the
concatenation directly, so the joining memmove is never issued. Both conditions
are the ones the call site ALREADY evaluates, and m <= 65535 bounds
cx + cx1 <= m <= 65535 -- inside the u16 fast engine's range. */
const int SPLIT = (FUSE_OK && D == KT - 3 && m <= 65535);
{
// ---- 2 INDEPENDENT MERGE-PATH CHAINS (key-boundary split) -------------
int dd[3], jj[3];
dd[0] = 0; jj[0] = 0;
{ int xi, xj; kg_kbound(retA, h, retB, r, (int)(((long)m * 1) / 2), &xi, &xj);
if (xi + xj < dd[0]) { xi = dd[0] - jj[0]; xj = jj[0]; }
dd[1] = xi + xj; jj[1] = xj; }
dd[2] = m; jj[2] = r;
long ia1 = dd[1] - jj[1], ib1 = jj[1]; cx1 = 0;
u64 *cb = SPLIT ? CROSSB : (cross + dd[1]);
int nmin = dd[1] - dd[0];
if (dd[2] - dd[1] < nmin) nmin = dd[2] - dd[1];
int u = 0;
for (; u < nmin; u++) {
KG_MSTEP(ia, ib, cx, cross, u);
KG_MSTEP(ia1, ib1, cx1, cb, dd[1] + u);
}
for (; u < dd[1]; u++) KG_MSTEP(ia, ib, cx, cross, u);
for (u = nmin; u < dd[2] - dd[1]; u++) KG_MSTEP(ia1, ib1, cx1, cb, dd[1] + u);
long cp = cx;
if (!SPLIT && cx1) memmove(cross + cp, cross + dd[1], sizeof(u64) * (size_t)cx1);
cp += cx1;
cx = cp;
}
if (FUSE_OK && D == KT - 3) {
u32 *dirt = (u32 *)ABtop; ABtop += (cx + 3) / 2 + 1;
if (SPLIT) ctr_solve_split2(cross, (int)(cx - cx1), CROSSB, (int)cx1, dirt);
else ctr_solve(cross, cx, dirt);
} else {
u64 *cres = ABtop; ABtop += cx;
solve_t<KT, D + 1, false>(cross, cx, cres);
}
ABtop = save;
}
template <int KT, int D, bool NO>
static inline void solve_t(u64 *seq, int m, u64 *out) {
if (m <= 1) return;
solve_b<KT, D, NO>(seq, m, out, TB<(D >= KT - 1), NO>());
}
// ====================================================================================
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 && (u32)n <= CT_LIM) {
const u32 *XL = XD[KT - 1];
memset(CT0, 0, (size_t)n);
int ok = 1;
for (int i = 0; i < n; i++) if (++CT0[XL[i]] > 127) { ok = 0; break; }
memset(CT0, 0, (size_t)n);
/* q9: the two u8 levels must not wrap. The counter is reset per ctr_solve
call and a call's inserts are a SUBSET of all 2n events, so the GLOBAL
bucket counts bound them. CT1 u8 needs <= 255, CT2 u16 needs <= 65535. */
if (ok) {
memset(GSC, 0, sizeof(u32) * ((size_t)(((unsigned)n) >> 5) + 64));
for (int i = 0; i < n && ok; i++) {
u32 b = XL[i] >> 5;
if (GSC[b] >= 255u) ok = 0; else GSC[b]++;
}
memset(GSC, 0, sizeof(u32) * ((size_t)(((unsigned)n) >> 10) + 16));
for (int i = 0; i < n && ok; i++) {
u32 b = XL[i] >> 10;
if (GSC[b] >= 65535u) ok = 0; else GSC[b]++;
}
memset(GSC, 0, sizeof(u32) * ((size_t)(((unsigned)n) >> 14) + 32));
for (int i = 0; i < n && ok; i++) {
u32 b = XL[i] >> 14;
if (GSC[b] >= 65535u) ok = 0; else GSC[b]++;
}
}
FUSE_OK = ok;
}
{
int bit = 0;
while ((n >> bit) >= LEAFNB && bit < 30) bit++;
LFSH = bit; LFD = 0;
}
G_LIM = (int)(((u32)(n - 1)) >> (LFSH + LFD)) + 2;
sort0(n);
const u32 *X0 = XD[0], *X1 = XD[1];
u64 *ev = AB; // top-level event array: 2n entries
int m = 0, g = 0;
while (g < n) {
int g1 = g;
u32 v = X0[g_ord[g]];
while (g1 < n && X0[g_ord[g1]] == v) g1++;
const u32 *X2 = XD[2];
for (int t = g; t < g1; t++) { int q = (int)g_ord[t]; ev[m++] = mk2(X1[q], 0, X2[q], (u32)q); }
for (int t = g; t < g1; t++) { int q = (int)g_ord[t]; ev[m++] = mk2(X1[q], 1, X2[q], (u32)q); }
g = g1;
}
ABtop = AB + m;
u64 *res = ABtop; ABtop += m;
solve_t<KT, 0, false>(ev, m, res);
}
} // namespace kg
static u64 g_arena[KG_ARENA];
static u32 g_ord_st[1048576], g_x_st[1048576];
static int g_ready = 0;
static void setup() {
if (g_ready) return;
g_ready = 1;
kg::mt_init();
kg::AB = g_arena;
kg::g_ord = g_ord_st;
kg::g_x = g_x_st;
}
void count_4d(int n, const unsigned *x[4], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run_t<4>(n, (const unsigned **)x, out);
}
void count_5d(int n, const unsigned *x[5], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run_t<5>(n, (const unsigned **)x, out);
}
#pragma GCC pop_options
#include <cstring>
#include <cstdio>
#include <cstdlib>
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 276.672 ms | 28 MB + 88 KB | Accepted | Score: 100 | 显示更多 |