// ARMVAR KG_BASE=2048 MF_CTRPF=32
#define KG_BASE 2816
#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 0
#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 4096
#endif
#define KG_MAXK 5
#ifndef E6_FUSE
#define E6_FUSE 1
#endif
#ifndef E6_KB2
#define E6_KB2 0
#endif
#define KG_ARENA (32u << 20)
#ifndef E6_KB1
#define E6_KB1 0
#endif
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
/* ---- lane_e6v: ONE contiguous block, fixed displacements off a single base ---- */
#define CT_SZ1 ((CT_MAXV / 32) + 8)
#define CT_SZ2 ((CT_MAXV / 256) + 8)
#define CT_SZ3 ((CT_MAXV / 2048) + 8)
#define CT_SZ4 ((CT_MAXV / 16384) + 8)
#define CT_SZ5 ((CT_MAXV / 131072) + 8)
#define CT_AL 64
#define CT_ALN(x) (((x) + (CT_AL - 1)) & ~((size_t)(CT_AL - 1)))
#define CT_O0 ((size_t)0)
#define CT_O1 (CT_ALN(CT_O0 + (size_t)CT_MAXV))
#define CT_O2 (CT_ALN(CT_O1 + (size_t)2 * CT_SZ1))
#define CT_O3 (CT_ALN(CT_O2 + (size_t)2 * CT_SZ2))
#define CT_O4 (CT_ALN(CT_O3 + (size_t)4 * CT_SZ3))
#define CT_O5 (CT_ALN(CT_O4 + (size_t)4 * CT_SZ4))
#define CT_TOT (CT_ALN(CT_O5 + (size_t)4 * CT_SZ5) + 64u)
typedef u8 ct_a8 __attribute__((__may_alias__));
typedef u16 ct_a16 __attribute__((__may_alias__));
typedef u32 ct_a32 __attribute__((__may_alias__));
static ct_a8 CTBLK[CT_TOT] __attribute__((aligned(64)));
#define CT0 (CTBLK + CT_O0)
#define CT1 ((ct_a16 *)(void *)(CTBLK + CT_O1))
#define CT2 ((ct_a16 *)(void *)(CTBLK + CT_O2))
#define CT3 ((ct_a32 *)(void *)(CTBLK + CT_O3))
#define CT4 ((ct_a32 *)(void *)(CTBLK + CT_O4))
#define CT5 ((ct_a32 *)(void *)(CTBLK + CT_O5))
static int FUSE_OK;
static const u8 MT8[32][32] __attribute__((aligned(32))) = {
{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,0},
{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0},
};
static const u16 MT16[8][8] __attribute__((aligned(16))) = {
{0,0,0,0,0,0,0,0},
{0xFFFF,0,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0,0},
{0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0xFFFF,0},
};
static const u32 MT32[8][8] __attribute__((aligned(32))) = {
{0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFFFFFu,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFFFFFu,0xFFFFFFFFu,0u,0u,0u,0u,0u,0u},
{0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0u,0u,0u,0u,0u},
{0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0u,0u,0u,0u},
{0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0u,0u,0u},
{0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0u,0u},
{0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0xFFFFFFFFu,0u},
};
static inline u32 su8_32(const u8 *p, int cnt) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
__m256i m = _mm256_load_si256((const __m256i *)MT8[cnt & 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_load_si128((const __m128i *)MT16[cnt & 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_load_si256((const __m256i *)MT32[cnt & 7]));
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
static inline void ctr_ins(u32 v) {
CT0[v]++; CT1[v >> 5]++; CT2[v >> 8]++; CT3[v >> 11]++; CT4[v >> 14]++; CT5[v >> 17]++;
}
static inline void ctr_clr(u32 v) { CT0[v]--; }
static int RSZN;
static inline void ctr_reset_hi(void) {
/* era-correct extent: n/div + 8, the declaration formula with n in place of CT_MAXV */
memset(CT1, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 5) + 8));
memset(CT2, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 8) + 8));
memset(CT3, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 11) + 8));
memset(CT4, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 14) + 8));
memset(CT5, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 17) + 8));
}
static const int IDXV[8] __attribute__((aligned(32))) = {0, 2, 4, 6, 7, 7, 7, 7};
static const signed char LANE8[32] __attribute__((aligned(32))) = {
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};
static const short LANE16[8] __attribute__((aligned(16))) = {0,1,2,3,4,5,6,7};
static const int LANE32[8] __attribute__((aligned(32))) = {0,1,2,3,4,5,6,7};
#define MK32(k) _mm256_cmpgt_epi32(_mm256_set1_epi32((int)((k) & 7)), _mm256_load_si256((const __m256i *)LANE32))
#define MK16(k) _mm_cmpgt_epi16(_mm_set1_epi16((short)((k) & 7)), _mm_load_si128((const __m128i *)LANE16))
#define MK8(k) _mm256_cmpgt_epi8(_mm256_set1_epi8((signed char)((k) & 31)), _mm256_load_si256((const __m256i *)LANE8))
static inline u32 ctr_qry(u32 v) {
u32 k5 = v >> 17;
__m256i K = _mm256_setzero_si256();
__m128i L = _mm_set1_epi16(1);
__m256i a5 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)CT5),
_mm256_load_si256((const __m256i *)MT32[k5 & 7]));
__m256i a4 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT4 + ((v >> 17) << 3))),
_mm256_load_si256((const __m256i *)MT32[(v >> 14) & 7]));
__m256i a3 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT3 + ((v >> 14) << 3))),
_mm256_load_si256((const __m256i *)MT32[(v >> 11) & 7]));
__m128i x2 = _mm_and_si128(_mm_loadu_si128((const __m128i *)(CT2 + ((v >> 11) << 3))),
_mm_load_si128((const __m128i *)MT16[(v >> 8) & 7]));
__m128i x1 = _mm_and_si128(_mm_loadu_si128((const __m128i *)(CT1 + ((v >> 8) << 3))),
_mm_load_si128((const __m128i *)MT16[(v >> 5) & 7]));
__m256i x0 = _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u))), MK8(v));
__m256i acc = _mm256_add_epi32(a5, a4);
acc = _mm256_add_epi32(acc, a3);
acc = _mm256_add_epi32(acc, _mm256_inserti128_si256(K, _mm_madd_epi16(x2, L), 0));
acc = _mm256_add_epi32(acc, _mm256_inserti128_si256(K, _mm_madd_epi16(x1, L), 0));
acc = _mm256_add_epi32(acc, _mm256_permutevar8x32_epi32(
_mm256_sad_epu8(x0, K), _mm256_load_si256((const __m256i *)IDXV)));
__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);
}
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;
// ================= lane_red SIDE ARRAY (mkside.py) =================
// g_s3 is parallel to the arena AB: g_s3[(p - AB)] is p's side value (KT==4: dim3).
static u32 *g_s3;
#define S3Q(p) (g_s3 + (size_t)((p) - AB))
// ===================================================================
// ===== 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 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) {
int nd = 0;
u32 *dp = dirt;
for (int i = 0; i < cx; i++) {
#if MF_CTRPF > 0
if (i + MF_CTRPF < cx) {
u32 vf = E_KA(cross[i + MF_CTRPF]);
__builtin_prefetch(&CT0[vf], 1, 0);
__builtin_prefetch(&CT1[vf >> 5], 1, 0);
__builtin_prefetch(&CT2[vf >> 8], 1, 0);
}
#endif
u64 v = cross[i];
u32 val = E_KB(v); // lane_g7: the fuse level now emits the VERBATIM event, so
// the dim K-1 counter key rides in KB (see solve_b).
if (v & (1ull << 42)) { ctr_ins(val); *dp++ = val; }
else OUT[E_ID(v)] += ctr_qry(val);
}
nd = (int)(dp - 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();
}
// 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)
#ifndef MF_S3
#define MF_S3 1
#endif
template<int ND, bool TS3 = false>
static void base_case_t(u64 *seq, int m, int d, u64 *out, int needOut) {
const int nd = ND;
__attribute__((aligned(32)))
u32 IK[KG_BASE + 8], IC0[KG_BASE + 8], IC1[KG_BASE + 8], IC2[KG_BASE + 8];
u32 QK[KG_BASE + 8], QC0[KG_BASE + 8], QC1[KG_BASE + 8], QC2[KG_BASE + 8];
u32 QID[KG_BASE + 8], QN[KG_BASE + 8];
const u32 *X1 = nd > 1 ? XD[d + 3] : 0, *X2 = nd > 2 ? XD[d + 4] : 0;
const u32 *S3S = TS3 ? S3Q(seq) : 0; // lane_red: dim3 arrives sequentially
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]);
/* lane_e6u PF: when TS3 supplies c1, X1 is NEVER read in this leaf -- the shipped
prefetch fetches a 4 MB array that no load consumes. Gate it on the same
condition that gates the load. */
if (nd > 1 && !TS3) __builtin_prefetch(&X1[idf], 0, 1);
if (nd > 2) __builtin_prefetch(&X2[idf], 0, 1);
}
#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 (TS3) c1 = S3S[t];
else 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;
QK[nq] = k; QC0[nq] = c0; QC1[nq] = c1; QC2[nq] = c2; QID[nq] = (u32)id; QN[nq] = (u32)ni; nq += 1 - isins;
}
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
int q = 0;
#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)QN[q], nB = (int)QN[q + 1],
nC = (int)QN[q + 2], nD = (int)QN[q + 3];
int n = nA; if (nB < n) n = 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]);
#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)QC1[q]), v1B = _mm256_set1_epi32((int)QC1[q + 1]),
v1C = _mm256_set1_epi32((int)QC1[q + 2]),v1D = _mm256_set1_epi32((int)QC1[q + 3]);
#endif
__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));
}
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));
mD = _mm256_and_si256(mD, _mm256_cmpgt_epi32(v1D, C1v));
}
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
accC = _mm256_sub_epi32(accC, mC);
accD = _mm256_sub_epi32(accD, mD);
}
L4Q_TAIL(accA, kvA, v0A, v1A, QN[q], j, QID[q]);
L4Q_TAIL(accB, kvB, v0B, v1B, QN[q + 1], j, QID[q + 1]);
L4Q_TAIL(accC, kvC, v0C, v1C, QN[q + 2], j, QID[q + 2]);
L4Q_TAIL(accD, kvD, v0D, v1D, QN[q + 3], j, QID[q + 3]);
}
}
#endif
#if MF_LEAFG >= 3
// ================= THREE-QUERY MERGED LEAF (lane e6w) =================
// Gated on `nd == 2` BY REGISTER ARITHMETIC: (nd+1) + k*((nd+1)+1) = 3 + 3*4 = 15 <= 16.
// This is EXACTLY the population m188's `nd <= 1` gate leaves out, and at KT=4 it is
// the D==0 leaf -- the largest one. k=4 needs 19 here; k=3 fits.
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; if (nB < n) n = 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);
__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));
__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);
}
L4Q_TAIL(accA, kvA, v0A, v1A, QN[q], j, QID[q]);
L4Q_TAIL(accB, kvB, v0B, v1B, QN[q + 1], j, QID[q + 1]);
L4Q_TAIL(accC, kvC, v0C, v1C, QN[q + 2], j, QID[q + 2]);
}
}
#endif
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)(E_KA(seq[t]) >> (LFSH + LFD)) + 1]++;
for (int b = 0; b + 1 < LB; b++) cnt[b + 1] += cnt[b];
u32 *S3O = TS3 ? S3Q(out) : 0;
for (int t = 0; t < m; t++) {
u64 v = seq[t];
int k = cnt[(int)(E_KA(v) >> (LFSH + LFD))]++;
out[k] = v;
if (TS3) S3O[k] = S3S[t];
}
int lo = 0;
for (int b = 0; b + 1 < LB; b++) {
int hi = (int)cnt[b];
if (hi - lo > 1)
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
u32 sv = TS3 ? S3O[t] : 0;
int q = t - 1;
while (q >= lo && v < out[q]) { out[q + 1] = out[q]; if (TS3) S3O[q + 1] = S3O[q]; q--; }
out[q + 1] = v;
if (TS3) S3O[q + 1] = sv;
}
lo = hi;
}
}
#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> 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 (E_TY(v)) cnt++; // insert
else OUT[E_ID(v)] += (u32)cnt; // query
}
}
// ================= lane_e6w: THE REPACK (see mkrepack.py header) =================
// At (KT==4, D==0) the side array holds dim2 and the event holds [dim1][TY][dim3][ID],
// so the cross stamp is one shift and one mask instead of an eleven-instruction promotion.
#define S3_STAMP(v) ((((v) << 21) & (0xFFFFFull << 43)) | ((v) & (1ull << 42)) | ((v) & 0xFFFFFull))
#define S3X(v, sv) ((KT == 4 && D == 0) \
? ((((u64)(u32)(sv)) << 43) | ((v) & 0x7FFFFCFFFFFull)) \
: (GAT ? (S3_STAMP(v) | ((u64)XN[E_ID(v)] << 22)) \
: (PROM ? S3_STAMP(v) : (v))))
// ==================================================================================
// ---- 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 <= ((E6_KB1 && D >= 1) ? E6_KB1 : ((E6_KB2 && D >= 2) ? E6_KB2 : KG_BASE))) {
#if MF_S3
if (KT == 4 && D == 0) base_case_t<NCL< KT - (D + 2) >::v, true>(seq, m, D, out, needOut);
else base_case_t<NCL< KT - (D + 2) >::v, false>(seq, m, D, out, needOut);
#else
base_case_t<NCL< KT - (D + 2) >::v>(seq, m, D, out, needOut);
#endif
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;
// 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 int FUS = (D == KT - 3);
const u32 *XN = XD[(D + 3 <= KG_MAXK - 1) ? (D + 3) : 0];
int ia = 0, ib = 0, cx = 0, t = 0;
while (ia < h && ib < r) {
#if KG_PF
if (!(ia & (KG_PF - 1)) && ia + KG_PF < h) __builtin_prefetch(&XN[E_ID(retA[ia + KG_PF])], 0, 1);
if (!(ib & (KG_PF - 1)) && ib + KG_PF < r) __builtin_prefetch(&XN[E_ID(retB[ib + KG_PF])], 0, 1);
#endif
#if MF_S3
// lane_red: the side load's address depends on ia/ib ALONE, so it issues with the element
// load instead of after the compare -- that is the whole point (see header).
const u32 sa = (KT == 4 && D == 0) ? S3Q(retA)[ia] : 0u;
const u32 sb = (KT == 4 && D == 0) ? S3Q(retB)[ib] : 0u;
#endif
u64 ea = retA[ia], eb = retB[ib];
int takeA = (ea < eb);
ia += takeA; ib += 1 - takeA;
u64 v = takeA ? ea : eb;
#if MF_S3
const u32 sv = takeA ? sa : sb;
if (KT == 4 && D == 0) { if (needOut) S3Q(out)[t] = sv; }
#endif
if (needOut) out[t] = v;
t++;
int role = (int)E_TY(v);
#if KG_COND
if (keep_test(takeA, role)) {
cross[cx] = S3X(v, sv);
cx++;
}
#else
u64 c = FUS ? v : S3X(v, sv);
cross[cx] = c;
cx += keep_test(takeA, role);
#endif
}
while (ia < h) {
int _i = ia++;
u64 v = retA[_i];
#if MF_S3
u32 _s = (KT == 4 && D == 0) ? S3Q(retA)[_i] : 0u;
if (KT == 4 && D == 0) { if (needOut) S3Q(out)[t] = _s; }
#endif
if (needOut) out[t] = v;
t++;
if (E_TY(v)) cross[cx++] = FUS ? v : S3X(v, _s);
}
while (ib < r) {
int _i = ib++;
u64 v = retB[_i];
#if MF_S3
if (KT == 4 && D == 0) { if (needOut) S3Q(out)[t] = S3Q(retB)[_i]; }
#endif
if (needOut) out[t] = v;
t++;
if (!E_TY(v)) cross[cx++] = FUS ? v : S3X(v, (KT == 4 && D == 0) ? S3Q(retB)[_i] : 0u);
}
if (FUSE_OK && D == KT - 3) {
u32 *dirt = (u32 *)ABtop; ABtop += (cx + 3) / 2 + 1;
ctr_solve(cross, cx, dirt);
} else {
/* lane_g7: when the fuse does NOT fire (FUSE_OK == 0) the D+1 level still
wants the promoted key, so restore the shipped form here. Compile-time
dead for every D != KT-3, and never executed when FUSE_OK == 1. */
if (D == KT - 3) for (int _z = 0; _z < cx; _z++) cross[_z] = S3_STAMP(cross[_z]);
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;
/* ================= lane_e6u S0: 2-pass LSD radix sort0 =================
key = ((u64)dim0 << 20) | index -- sorting it ascending sorts by dim0 and, for
equal dim0, by index, i.e. EXACTLY the stable counting sort the shipped code does
with `g_ord[cnt[X0[i]]++] = i` over ascending i.
Two passes of 10 bits; both tables are 4 KB (L1) and both arrays are streamed. */
static u64 *g_s0a, *g_s0b;
static u64 *g_pk; // lane_e6u PK: (dim1<<40)|(dim2<<20)|dim3, by point id
static u32 *g_x0s; // lane_e6u PK: dim0 of the point at each g_ord position
static void sort0(int n) {
const u32 *X0 = XD[0];
u64 *A = g_s0a, *B = g_s0b;
u32 cnt[1024];
for (int i = 0; i < 1024; i++) cnt[i] = 0;
for (int i = 0; i < n; i++) cnt[X0[i] & 1023u]++;
{ u32 s = 0; for (int i = 0; i < 1024; i++) { u32 c = cnt[i]; cnt[i] = s; s += c; } }
for (int i = 0; i < n; i++) { u32 k = X0[i]; A[cnt[k & 1023u]++] = ((u64)k << 20) | (u32)i; }
for (int i = 0; i < 1024; i++) cnt[i] = 0;
for (int i = 0; i < n; i++) cnt[(u32)(A[i] >> 30)]++;
{ u32 s = 0; for (int i = 0; i < 1024; i++) { u32 c = cnt[i]; cnt[i] = s; s += c; } }
for (int i = 0; i < n; i++) { u64 v = A[i]; B[cnt[(u32)(v >> 30)]++] = v; }
u32 *o = g_ord, *xs = g_x0s;
for (int i = 0; i < n; i++) { u64 v = B[i]; o[i] = (u32)(v & 0xFFFFFu); xs[i] = (u32)(v >> 20); }
}
/* ======================================================================= */
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 & E6_FUSE;
}
{
int bit = 0;
while ((n >> bit) >= LEAFNB && bit < 30) bit++;
LFSH = bit; LFD = 0;
}
G_LIM = (int)(((u32)(n - 1)) >> (LFSH + LFD)) + 2;
{
// lane_e6u PK: one sequential pass over the three arrays the event build needs.
const u32 *Y1 = XD[1], *Y2 = XD[2], *Y3 = XD[3];
u64 *pk = g_pk;
for (int i = 0; i < n; i++)
pk[i] = ((u64)Y1[i] << 40) | ((u64)Y2[i] << 20) | (u64)Y3[i];
}
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 = g_x0s[g];
while (g1 < n && g_x0s[g1] == v) g1++;
const u32 *X2 = XD[2];
const u32 *X3 = XD[3];
#if MF_S3
if (KT == 4) {
u32 *es = S3Q(ev);
const int gk = g1 - g;
for (int t = g; t < g1; t++) {
int q = (int)g_ord[t];
u64 p = g_pk[q]; // ONE gather, not three
u64 e = ((p & (0xFFFFFull << 40)) << 3) | (u64)q | ((p & 0xFFFFFull) << 22);
es[m + (t - g)] = (u32)((p >> 20) & 0xFFFFFull);
es[m + gk + (t - g)] = (u32)((p >> 20) & 0xFFFFFull);
ev[m + (t - g)] = e;
ev[m + gk + (t - g)] = e | (1ull << 42);
}
m += 2 * gk;
} else
#endif
{
/* lane_e6u EV: the shipped pair gathers X1[q]/X2[q] ONCE PER LOOP and does the
same index load twice. One pass, two events, ONE gather each -- and the
shipped event order (all queries, then all inserts, groups unchanged) is kept. */
const int gk = g1 - g;
for (int t = g; t < g1; t++) {
int q = (int)g_ord[t];
u64 p = g_pk[q]; // ONE gather, not two
u64 e = ((p & (0xFFFFFull << 40)) << 3) | (u64)q | (((p >> 20) & 0xFFFFFull) << 22);
ev[m + (t - g)] = e;
ev[m + gk + (t - g)] = e | (1ull << 42);
}
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_s3_st[KG_ARENA];
static u32 g_ord_st[1048576], g_x_st[1048576];
static u64 g_s0a_st[1048576], g_s0b_st[1048576];
static u64 g_pk_st[1048576];
static u32 g_x0s_st[1048576];
static int g_ready = 0;
static void setup() {
if (g_ready) return;
g_ready = 1;
kg::AB = g_arena;
kg::g_s3 = g_s3_st;
kg::g_ord = g_ord_st;
kg::g_x = g_x_st;
kg::g_s0a = g_s0a_st;
kg::g_s0b = g_s0b_st;
kg::g_pk = g_pk_st;
kg::g_x0s = g_x0s_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 | 816.011 ms | 127 MB + 960 KB | Accepted | Score: 100 | 显示更多 |