#ifndef MF_CTRPF
#define MF_CTRPF 0
#endif
#ifndef MF_CTR5
#define MF_CTR5 0
#endif
#ifndef MF_LEAFLP
#define MF_LEAFLP 16
#endif
#ifndef MF_LEAFG
#define MF_LEAFG 2
#endif
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
// duck.ac 1010 (4D) / 1011 (5D) -- "count for every point the number of points that are
// smaller in *every* dimension" (n = 300000, coords in [0,n), strict comparisons).
//
// Algorithm: D&C over dimensions on an event sequence, keys carried inline.
// * 2n events: one query event + one insert event per point, keyed by dim 0
// (a query event is stored *before* an insert event of the same key, which encodes
// the strict-inequality rule).
// * solve(seq,m,d): seq is ordered by dim d and every event carries the coordinate of
// dim d+1 in bits 32..63 of its u64. Split by position, recurse both halves, then
// merge the two halves by the carried key (a plain unsigned `<` on the packed u64,
// so equal keys keep query-before-insert order); the same pass collects the cross
// events (left inserts + right queries), re-keys them to dim d+2 and recurses one
// dimension deeper. Ordering by dim d with that tie rule makes "insert before query"
// exactly "all dims so far strictly smaller", and equal coordinates are never counted.
// * d == K-1: every insert-before-query pair now qualifies, so one O(m) running counter
// finishes the innermost dimension (no log factor).
// * m <= 32: direct O(m^2) counting plus an insertion sort.
//
// Local timings (i5-6500 @ 3.2GHz, best of 3, n = 300000 uniform random):
// 4D 1.57 s 5D 3.85 s
// (Board #1: 508.596 ms / 745.358 ms -- this is still 3x / 5x behind.)
// k-dimensional dominance counting: out[i] = #{j : all dims of j < dims of i}
//
// Event-sequence CDQ with keys carried inline.
// event (u64) = (key << 32) | (ins << 31) | id ins = 1 for an "insert" event
// key = coordinate of dim (level+1): the field used by this level's merge.
// The packing makes a plain unsigned comparison put, for equal keys, the query
// event first -- exactly the strict-inequality rule (equal coords must not count).
//
// solve(seq, m, d): seq is ordered by dim d (tie rule above), elements carry the
// key for dim d+1; count pairs (insert before query) with dims d+1..K-1 strictly
// smaller; return `out` = the same events merged by the carried key (= dim d+1).
// Cross = A-inserts + B-queries re-keyed to dim d+2, recursed at level d+1.
// d == K-1 : every (insert before query) pair qualifies -> running counter, O(m).
#include <cstring>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC target("avx2,popcnt,bmi2")
#include <cstdint>
typedef unsigned u32;
typedef uint64_t u64;
typedef unsigned char u8;
typedef unsigned short u16;
#ifndef KG_COND
#define KG_COND 0
#endif
#ifndef KG_PF
#define KG_PF 0
#endif
#ifndef KG_BASE
#define KG_BASE 2560
#endif
#ifndef LEAFNB
#define LEAFNB 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 inline u32 su8_32(const u8 *p, int cnt) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
__m256i m = _mm256_cmpgt_epi8(_mm256_set1_epi8((char)cnt),
_mm256_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31));
__m256i x = _mm256_and_si256(v, m);
__m256i s = _mm256_sad_epu8(x, _mm256_setzero_si256());
__m128i t = _mm_add_epi64(_mm256_castsi256_si128(s), _mm256_extracti128_si256(s, 1));
return (u32)((unsigned long long)_mm_cvtsi128_si64(t) + (unsigned long long)_mm_extract_epi64(t, 1));
}
static inline u32 su16_8(const u16 *p, int cnt) {
__m128i v = _mm_loadu_si128((const __m128i *)p);
__m128i x = _mm_and_si128(v, _mm_cmpgt_epi16(_mm_set1_epi16((short)cnt),
_mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7)));
__m128i s = _mm_madd_epi16(x, _mm_set1_epi16(1));
s = _mm_add_epi32(s, _mm_shuffle_epi32(s, _MM_SHUFFLE(1, 0, 3, 2)));
s = _mm_add_epi32(s, _mm_shuffle_epi32(s, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(s);
}
static inline u32 su32_8(const u32 *p, int cnt) {
__m256i v = _mm256_loadu_si256((const __m256i *)p);
__m256i x = _mm256_and_si256(v, _mm256_cmpgt_epi32(_mm256_set1_epi32(cnt),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7)));
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
static inline void ctr_ins(u32 v) {
CT0[v]++; CT1[v >> 5]++; CT2[v >> 8]++; CT3[v >> 11]++; CT4[v >> 14]++; CT5[v >> 17]++;
}
/* ARM `CT`: RESET THE HIERARCHY INSTEAD OF DECREMENTING IT.
ctr_clr mirrored all six of ctr_ins's read-modify-writes, so every element cost
SIX random RMWs to insert and SIX more to clear. Five of those six arrays are
SMALL -- CT1 18 KB, CT2 2.4 KB, CT3 592 B, CT4 76 B, CT5 16 B -- and they are
only ever read through ctr_qry, which sums them; restoring them to their initial
zero is a pure RESET, not an inverse walk. Zeroing all five together is one
contiguous memset of ~21 KB per fusion call, against 5*nd scattered decrements.
CT0 stays a per-element decrement because it is the one level that is genuinely
large (293 KB) and densely addressed, so walking `dirt` still beats clearing it.
Semantics are unchanged: ctr_qry reads exactly the same six levels and the
structure is left in its initial state for the next call, which is what the
decrements were for. */
static inline void ctr_clr(u32 v) { CT0[v]--; }
static int RSZN;
static inline void ctr_reset_hi(void) {
/* era-correct extent: n/div + 8, the declaration formula with n in place of CT_MAXV */
memset(CT1, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 5) + 8));
memset(CT2, 0, sizeof(u16) * ((size_t)((unsigned)RSZN >> 8) + 8));
memset(CT3, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 11) + 8));
memset(CT4, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 14) + 8));
memset(CT5, 0, sizeof(u32) * ((size_t)((unsigned)RSZN >> 17) + 8));
}
static inline u32 ctr_qry(u32 v) {
u32 k5 = v >> 17, s = 0;
#if MF_CTR5
{
__m256i tv = _mm256_loadu_si256((const __m256i *)CT5);
__m256i tm = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)k5),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
__m128i ta = _mm_add_epi32(_mm256_castsi256_si128(_mm256_and_si256(tv, tm)),
_mm256_extracti128_si256(_mm256_and_si256(tv, tm), 1));
ta = _mm_add_epi32(ta, _mm_shuffle_epi32(ta, _MM_SHUFFLE(1, 0, 3, 2)));
ta = _mm_add_epi32(ta, _mm_shuffle_epi32(ta, _MM_SHUFFLE(2, 3, 0, 1)));
s += (u32)_mm_cvtsi128_si32(ta);
}
#else
for (u32 j = 0; j < k5; j++) s += CT5[j];
#endif
/* ARM `QFUS`: one horizontal reduction for all five levels.
The boarded form ran a separate 3-shuffle reduction inside each of the five
su*_8 helpers; the five results are merely ADDED, so the reduction can be done
once at the end at no semantic cost (all levels are nonneg and their sum is the
count, bounded by n <= CT_MAXV-1, so 32-bit lanes cannot overflow). Work count
is unchanged per level; five serial shuffle chains become one. */
__m128i acc = _mm_setzero_si128();
{
__m256i tv = _mm256_loadu_si256((const __m256i *)(CT4 + ((v >> 17) << 3)));
__m256i tm = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)((v >> 14) & 7)),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
__m256i x = _mm256_and_si256(tv, tm);
acc = _mm_add_epi32(acc, _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1)));
}
{
__m256i tv = _mm256_loadu_si256((const __m256i *)(CT3 + ((v >> 14) << 3)));
__m256i tm = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)((v >> 11) & 7)),
_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
__m256i x = _mm256_and_si256(tv, tm);
acc = _mm_add_epi32(acc, _mm_add_epi32(_mm256_castsi256_si128(x), _mm256_extracti128_si256(x, 1)));
}
{
__m128i vv = _mm_loadu_si128((const __m128i *)(CT2 + ((v >> 11) << 3)));
__m128i x = _mm_and_si128(vv, _mm_cmpgt_epi16(_mm_set1_epi16((short)((v >> 8) & 7)),
_mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7)));
acc = _mm_add_epi32(acc, _mm_madd_epi16(x, _mm_set1_epi16(1)));
}
{
__m128i vv = _mm_loadu_si128((const __m128i *)(CT1 + ((v >> 8) << 3)));
__m128i x = _mm_and_si128(vv, _mm_cmpgt_epi16(_mm_set1_epi16((short)((v >> 5) & 7)),
_mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7)));
acc = _mm_add_epi32(acc, _mm_madd_epi16(x, _mm_set1_epi16(1)));
}
{
__m256i vv = _mm256_loadu_si256((const __m256i *)(CT0 + (v & ~31u)));
__m256i m = _mm256_cmpgt_epi8(_mm256_set1_epi8((char)(v & 31)),
_mm256_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31));
__m256i sd = _mm256_sad_epu8(_mm256_and_si256(vv, m), _mm256_setzero_si256());
acc = _mm_add_epi32(acc, _mm_add_epi32(_mm256_castsi256_si128(sd), _mm256_extracti128_si256(sd, 1)));
}
acc = _mm_add_epi32(acc, _mm_shuffle_epi32(acc, _MM_SHUFFLE(1, 0, 3, 2)));
acc = _mm_add_epi32(acc, _mm_shuffle_epi32(acc, _MM_SHUFFLE(2, 3, 0, 1)));
return s + (u32)_mm_cvtsi128_si32(acc);
}
static int K;
#ifndef LFAUTO
#define LFAUTO 0
#endif
#ifndef LF_D
#define LF_D 0
#endif
static int LFSH;
static int LFB;
static int LFD;
static const u32 *XD[KG_MAXK];
static 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;
}
static void solve(u64 *seq, int m, int d, u64 *out, int needOut);
// resolve a fusion-level cross: events carry the last dim (K-1) coordinate in bits 32..63
// and are ordered by dim K-2; a left-insert before a right-query with a strictly smaller
// last-dim value is a hit.
static void ctr_solve(u64 *cross, int cx, u32 *dirt) {
int nd = 0;
for (int i = 0; i < cx; i++) {
#if MF_CTRPF > 0
if (i + MF_CTRPF < cx) {
u32 vf = (u32)(cross[i + MF_CTRPF] >> 32);
__builtin_prefetch(&CT0[vf], 1, 0);
__builtin_prefetch(&CT1[vf >> 5], 1, 0);
__builtin_prefetch(&CT2[vf >> 8], 1, 0);
}
#endif
u64 v = cross[i];
u32 val = (u32)(v >> 32);
if ((v >> 31) & 1) { ctr_ins(val); dirt[nd++] = val; }
else OUT[eid(v)] += ctr_qry(val);
}
for (int i = 0; i < nd; i++) {
#if MF_CTRPF > 0
if (i + MF_CTRPF < nd) {
u32 vf = dirt[i + MF_CTRPF];
__builtin_prefetch(&CT0[vf], 1, 1);
__builtin_prefetch(&CT1[vf >> 5], 1, 1);
__builtin_prefetch(&CT2[vf >> 8], 1, 1);
}
#endif
ctr_clr(dirt[i]);
}
ctr_reset_hi();
}
// small base: compact insert list + SIMD comparison, then insertion sort by the key
static const u32 VMASK[8][8] = {
{0,0,0,0,0,0,0,0},
{~0u,0,0,0,0,0,0,0},
{~0u,~0u,0,0,0,0,0,0},
{~0u,~0u,~0u,0,0,0,0,0},
{~0u,~0u,~0u,~0u,0,0,0,0},
{~0u,~0u,~0u,~0u,~0u,0,0,0},
{~0u,~0u,~0u,~0u,~0u,~0u,0,0},
{~0u,~0u,~0u,~0u,~0u,~0u,~0u,0}};
static inline u32 hsum256(__m256i v) {
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(v), _mm256_extracti128_si256(v, 1));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
return (u32)_mm_cvtsi128_si32(t);
}
/* ===== ARM `DN`: SPECIALISE base_case ON THE REMAINING DIMENSION COUNT =====
`nd = K - (d + 2)` is a RUNTIME value, yet it gates work inside base_case's
INNERMOST loops: three `if (nd > 0/1/2)` tests sit between the AVX2 mask folds
(the loads of X0/X1/X2 and the `_mm256_cmpgt_epi32` chain), so every 8-element
block pays three compare-and-branch on a value that cannot change within the
call. `nd` is in [0, K-2] and K <= KG_MAXK = 5, so there are only FOUR cases;
dispatching once per base_case call and compiling `nd` as a constant removes
every one of those branches and lets the address selection of X0/X1/X2 fold too.
This is a WORK reduction -- instructions are deleted, not rescheduled -- which
is the class whose local sign has held on this row. */
template<int ND>
static void base_case_t(u64 *seq, int m, int d, u64 *out, int needOut) {
const int nd = ND;
__attribute__((aligned(32)))
u32 IK[KG_BASE + 8], IC0[KG_BASE + 8], IC1[KG_BASE + 8], IC2[KG_BASE + 8];
u32 QK[KG_BASE + 8], QC0[KG_BASE + 8], QC1[KG_BASE + 8], QC2[KG_BASE + 8];
u32 QID[KG_BASE + 8], QN[KG_BASE + 8];
const u32 *X0 = nd > 0 ? XD[d + 2] : 0, *X1 = nd > 1 ? XD[d + 3] : 0, *X2 = nd > 2 ? XD[d + 4] : 0;
int ni = 0, nq = 0;
for (int t = 0; t < m; t++) {
#if MF_LEAFLP > 0
if (t + MF_LEAFLP < m) {
int idf = eid(seq[t + MF_LEAFLP]);
if (nd > 0) __builtin_prefetch(&X0[idf], 0, 1);
if (nd > 1) __builtin_prefetch(&X1[idf], 0, 1);
if (nd > 2) __builtin_prefetch(&X2[idf], 0, 1);
}
#endif
u64 v = seq[t];
int id = eid(v);
u32 k = (u32)(v >> 32);
u32 c0 = 0, c1 = 0, c2 = 0;
if (nd > 0) c0 = X0[id];
if (nd > 1) c1 = X1[id];
if (nd > 2) c2 = X2[id];
int isins = (int)((v >> 31) & 1u);
IK[ni] = k; IC0[ni] = c0; IC1[ni] = c1; IC2[ni] = c2; ni += isins;
QK[nq] = k; QC0[nq] = c0; QC1[nq] = c1; QC2[nq] = c2; QID[nq] = (u32)id; QN[nq] = (u32)ni; nq += 1 - isins;
}
#if MF_LEAFG >= 2
// ---- paired-query leaf: one insert-block load serves two queries ----
int q = 0;
for (; q + 2 <= nq; q += 2) {
const int n = ((int)QN[q] < (int)QN[q + 1]) ? (int)QN[q] : (int)QN[q + 1];
const int nA = (int)QN[q], nB = (int)QN[q + 1];
__m256i kvA = _mm256_set1_epi32((int)QK[q]), kvB = _mm256_set1_epi32((int)QK[q + 1]);
__m256i v0A = _mm256_set1_epi32((int)QC0[q]), v0B = _mm256_set1_epi32((int)QC0[q + 1]);
__m256i v1A = _mm256_set1_epi32((int)QC1[q]), v1B = _mm256_set1_epi32((int)QC1[q + 1]);
__m256i v2A = _mm256_set1_epi32((int)QC2[q]), v2B = _mm256_set1_epi32((int)QC2[q + 1]);
__m256i accA = _mm256_setzero_si256(), accB = _mm256_setzero_si256();
int j = 0;
if (nd >= 3) {
for (; j + 8 <= n; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i C0v = _mm256_load_si256((const __m256i *)(IC0 + j));
__m256i C1v = _mm256_load_si256((const __m256i *)(IC1 + j));
__m256i C2v = _mm256_load_si256((const __m256i *)(IC2 + j));
__m256i mA = _mm256_and_si256(_mm256_cmpgt_epi32(kvA, IKv), _mm256_cmpgt_epi32(v0A, C0v));
__m256i mB = _mm256_and_si256(_mm256_cmpgt_epi32(kvB, IKv), _mm256_cmpgt_epi32(v0B, C0v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, C1v));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v2A, C2v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v2B, C2v));
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
}
} else {
for (; j + 8 <= n; j += 8) {
__m256i IKv = _mm256_load_si256((const __m256i *)(IK + j));
__m256i mA = _mm256_cmpgt_epi32(kvA, IKv);
__m256i mB = _mm256_cmpgt_epi32(kvB, IKv);
if (nd > 0) {
__m256i C0v = _mm256_load_si256((const __m256i *)(IC0 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v0A, C0v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v0B, C0v));
}
if (nd > 1) {
__m256i C1v = _mm256_load_si256((const __m256i *)(IC1 + j));
mA = _mm256_and_si256(mA, _mm256_cmpgt_epi32(v1A, C1v));
mB = _mm256_and_si256(mB, _mm256_cmpgt_epi32(v1B, C1v));
}
accA = _mm256_sub_epi32(accA, mA);
accB = _mm256_sub_epi32(accB, mB);
}
}
// per-query tails: each member finishes [j, QN[.]) on its own
{
int jj = j;
__m256i acc = accA;
for (; jj + 8 <= nA; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nA) {
int n = nA;
__m256i msk = _mm256_cmpgt_epi32(kvA, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0A, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1A, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2A, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[QID[q]] += hsum256(acc);
}
{
int jj = j;
__m256i acc = accB;
for (; jj + 8 <= nB; jj += 8) {
__m256i msk = _mm256_cmpgt_epi32(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
acc = _mm256_sub_epi32(acc, msk);
}
if (jj < nB) {
int n = nB;
__m256i msk = _mm256_cmpgt_epi32(kvB, _mm256_load_si256((const __m256i *)(IK + jj)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0B, _mm256_load_si256((const __m256i *)(IC0 + jj))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1B, _mm256_load_si256((const __m256i *)(IC1 + jj))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2B, _mm256_load_si256((const __m256i *)(IC2 + jj))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - jj]));
acc = _mm256_sub_epi32(acc, msk);
}
OUT[QID[q + 1]] += hsum256(acc);
}
}
#endif
for (; q < nq; q++) {
const int n = (int)QN[q];
__m256i kv = _mm256_set1_epi32((int)QK[q]);
__m256i v0 = _mm256_set1_epi32((int)QC0[q]);
__m256i v1 = _mm256_set1_epi32((int)QC1[q]);
__m256i v2 = _mm256_set1_epi32((int)QC2[q]);
__m256i acc = _mm256_setzero_si256();
int j = 0;
switch (nd) {
case 0:
for (; j + 8 <= n; j += 8)
acc = _mm256_sub_epi32(acc, _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j))));
break;
case 1:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
case 2:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
default:
for (; j + 8 <= n; j += 8) {
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC2 + j))));
acc = _mm256_sub_epi32(acc, msk);
}
break;
}
if (j < n) { // partial block, masked by a table lookup
__m256i msk = _mm256_cmpgt_epi32(kv, _mm256_load_si256((const __m256i *)(IK + j)));
if (nd > 0) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v0, _mm256_load_si256((const __m256i *)(IC0 + j))));
if (nd > 1) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v1, _mm256_load_si256((const __m256i *)(IC1 + j))));
if (nd > 2) msk = _mm256_and_si256(msk, _mm256_cmpgt_epi32(v2, _mm256_load_si256((const __m256i *)(IC2 + j))));
msk = _mm256_and_si256(msk, _mm256_loadu_si256((const __m256i *)VMASK[n - j]));
acc = _mm256_sub_epi32(acc, msk);
}
u32 c = hsum256(acc);
OUT[QID[q]] += c;
}
if (!needOut) return;
#if LEAFNB
{
const int NB = (LFAUTO && m <= LEAFNB) ? (m < 8 ? 8 : m) : LFB;
u32 cnt[LEAFNB + 2];
for (int b = 0; b <= NB; b++) cnt[b] = 0;
for (int t = 0; t < m; t++) cnt[(int)(((u32)(seq[t] >> 32)) >> (LFSH + LFD)) + 1]++;
for (int b = 0; b < NB; b++) cnt[b + 1] += cnt[b];
for (int t = 0; t < m; t++) { u64 v = seq[t]; out[cnt[(int)(((u32)(v >> 32)) >> (LFSH + LFD))]++] = v; }
int lo = 0;
for (int b = 0; b < NB; b++) {
int hi = (int)cnt[b];
if (hi - lo > 1)
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
int q = t - 1;
while (q >= lo && v < out[q]) { out[q + 1] = out[q]; q--; }
out[q + 1] = v;
}
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
}
/* ================== ARM BC1: bucketed 16-bit-lane base case for nd == 1 ==============
See BC1_NOTES.md. Coordinates are in [0,n) <= [0,2^17); write c = h*65536 + lo.
c_i < c_q <=> h_i < h_q OR (h_i == h_q AND lo_i < lo_q).
Bucket the inserts by (h_key, h_c0) -> 4 buckets; for a query in bucket (a,b) the
comparison of each dim is either free (h_i < h_q), a 16-bit lane compare (equal), or
an exclusion (h_i > h_q). lo is stored XOR 0x8000 so the SIGNED vpcmpgtw implements
the unsigned compare. Exact: no approximations anywhere. */
static const u16 BC1_TAIL[17][16] = {
{0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0u},
{0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu,0xFFFFu}
};
#define BC1_JUNK (KG_BASE + 8)
static u16 BC1_IK[4][BC1_JUNK + 1];
static u16 BC1_IC[4][BC1_JUNK + 1];
static u32 BC1_cn[4];
static u32 BC1_QK[BC1_JUNK + 1];
static u32 BC1_QC[BC1_JUNK + 1];
static u32 BC1_QID[BC1_JUNK + 1];
static u32 BC1_QAB[BC1_JUNK + 1];
static u16 BC1_N0[BC1_JUNK + 1], BC1_N1[BC1_JUNK + 1], BC1_N2[BC1_JUNK + 1], BC1_N3[BC1_JUNK + 1];
static inline u32 bc1_hsum(__m256i acc) {
__m256i s = _mm256_madd_epi16(acc, _mm256_set1_epi16(1));
__m128i t = _mm_add_epi32(_mm256_castsi256_si128(s), _mm256_extracti128_si256(s, 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);
}
/* one dim, 16 lanes */
static inline u32 bc1_one(const u16 *a, int len, __m256i vq) {
__m256i acc = _mm256_setzero_si256();
int j = 0;
for (; j + 16 <= len; j += 16) {
__m256i v = _mm256_loadu_si256((const __m256i *)(a + j));
acc = _mm256_sub_epi16(acc, _mm256_cmpgt_epi16(vq, v));
}
if (j < len) {
__m256i v = _mm256_loadu_si256((const __m256i *)(a + j));
__m256i m = _mm256_and_si256(_mm256_cmpgt_epi16(vq, v),
_mm256_loadu_si256((const __m256i *)BC1_TAIL[len - j]));
acc = _mm256_sub_epi16(acc, m);
}
return bc1_hsum(acc);
}
/* two dims, 16 lanes */
static inline u32 bc1_both(const u16 *ak, const u16 *ac, int len, __m256i vk, __m256i vc) {
__m256i acc = _mm256_setzero_si256();
int j = 0;
for (; j + 16 <= len; j += 16) {
__m256i xk = _mm256_loadu_si256((const __m256i *)(ak + j));
__m256i xc = _mm256_loadu_si256((const __m256i *)(ac + j));
__m256i m = _mm256_and_si256(_mm256_cmpgt_epi16(vk, xk), _mm256_cmpgt_epi16(vc, xc));
acc = _mm256_sub_epi16(acc, m);
}
if (j < len) {
__m256i xk = _mm256_loadu_si256((const __m256i *)(ak + j));
__m256i xc = _mm256_loadu_si256((const __m256i *)(ac + j));
__m256i m = _mm256_and_si256(_mm256_cmpgt_epi16(vk, xk), _mm256_cmpgt_epi16(vc, xc));
m = _mm256_and_si256(m, _mm256_loadu_si256((const __m256i *)BC1_TAIL[len - j]));
acc = _mm256_sub_epi16(acc, m);
}
return bc1_hsum(acc);
}
static void base_case_nd1(u64 *seq, int m, int d, u64 *out, int needOut) {
const u32 *X0 = XD[d + 2];
BC1_cn[0] = BC1_cn[1] = BC1_cn[2] = BC1_cn[3] = 0;
int q = 0;
for (int t = 0; t < m; t++) {
u64 v = seq[t];
int id = eid(v);
u32 k = (u32)(v >> 32);
u32 c0 = X0[id];
int isins = (int)((v >> 31) & 1u);
int B = (int)(((k >> 16) & 1u) * 2u + ((c0 >> 16) & 1u));
int ii = isins ? (int)BC1_cn[B] : BC1_JUNK;
BC1_IK[B][ii] = (u16)(k ^ 0x8000u);
BC1_IC[B][ii] = (u16)(c0 ^ 0x8000u);
BC1_cn[B] += (u32)isins;
int qi = isins ? BC1_JUNK : q;
BC1_QK[qi] = k; BC1_QC[qi] = c0; BC1_QID[qi] = (u32)id; BC1_QAB[qi] = (u32)B;
BC1_N0[qi] = (u16)BC1_cn[0]; BC1_N1[qi] = (u16)BC1_cn[1];
BC1_N2[qi] = (u16)BC1_cn[2]; BC1_N3[qi] = (u16)BC1_cn[3];
q += 1 - isins;
}
for (int i = 0; i < q; i++) {
u32 kq = BC1_QK[i], cq = BC1_QC[i];
__m256i vk = _mm256_set1_epi16((short)(u16)(kq ^ 0x8000u));
__m256i vc = _mm256_set1_epi16((short)(u16)(cq ^ 0x8000u));
u32 ab = BC1_QAB[i], cnt;
int n0 = BC1_N0[i], n1 = BC1_N1[i], n2 = BC1_N2[i], n3 = BC1_N3[i];
if (ab == 0) {
cnt = bc1_both(BC1_IK[0], BC1_IC[0], n0, vk, vc);
} else if (ab == 1) { /* a=0,b=1 */
cnt = bc1_one(BC1_IK[0], n0, vk) + bc1_both(BC1_IK[1], BC1_IC[1], n1, vk, vc);
} else if (ab == 2) { /* a=1,b=0 */
cnt = bc1_one(BC1_IC[0], n0, vc) + bc1_both(BC1_IK[2], BC1_IC[2], n2, vk, vc);
} else { /* a=1,b=1 */
cnt = (u32)n0 + bc1_one(BC1_IC[1], n1, vc) + bc1_one(BC1_IK[2], n2, vk)
+ bc1_both(BC1_IK[3], BC1_IC[3], n3, vk, vc);
}
OUT[BC1_QID[i]] += cnt;
}
if (!needOut) return;
#if LEAFNB
{
const int NB = (LFAUTO && m <= LEAFNB) ? (m < 8 ? 8 : m) : LFB;
u32 cnt[LEAFNB + 2];
for (int b = 0; b <= NB; b++) cnt[b] = 0;
for (int t = 0; t < m; t++) cnt[(int)(((u32)(seq[t] >> 32)) >> (LFSH + LFD)) + 1]++;
for (int b = 0; b < NB; b++) cnt[b + 1] += cnt[b];
for (int t = 0; t < m; t++) { u64 v = seq[t]; out[cnt[(int)(((u32)(v >> 32)) >> (LFSH + LFD))]++] = v; }
int lo = 0;
for (int b = 0; b < NB; b++) {
int hi = (int)cnt[b];
if (hi - lo > 1)
for (int t = lo + 1; t < hi; t++) {
u64 v = out[t];
int q2 = t - 1;
while (q2 >= lo && v < out[q2]) { out[q2 + 1] = out[q2]; q2--; }
out[q2 + 1] = v;
}
lo = hi;
}
}
#else
for (int t = 1; t < m; t++) {
u64 v = seq[t];
int p = t - 1;
while (p >= 0 && v < seq[p]) { seq[p + 1] = seq[p]; p--; }
seq[p + 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_nd1(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;
}
}
static void solve(u64 *seq, int m, int d, u64 *out, int needOut) {
if (m <= 1) return;
if (d >= K - 1) {
int cnt = 0;
for (int t = 0; t < m; t++) {
u64 v = seq[t];
if ((v >> 31) & 1) cnt++; // insert (bit31 = 1)
else OUT[eid(v)] += (u32)cnt; // query
}
return;
}
if (m <= KG_BASE) { base_case(seq, m, d, out, needOut); return; }
(void)0;
int h = m >> 1, r = m - h;
u64 *save = ABtop;
u64 *retA = ABtop; ABtop += h;
u64 *retB = ABtop; ABtop += r;
solve(seq, h, d, retA, 1);
solve(seq + h, r, d, retB, 1);
u64 *cross = ABtop; ABtop += m + 1;
const int rekey = (d + 2 <= K - 1);
const u32 *XN = XD[rekey ? d + 2 : 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[eid(retA[ia + KG_PF])], 0, 1);
if (!(ib & (KG_PF - 1)) && ib + KG_PF < r) __builtin_prefetch(&XN[eid(retB[ib + KG_PF])], 0, 1);
#endif
u64 ea = retA[ia], eb = retB[ib];
int takeA = (ea < eb);
ia += takeA; ib += 1 - takeA;
u64 v = takeA ? ea : eb;
if (needOut) out[t] = v;
t++;
int role = (int)((v >> 31) & 1);
#if KG_COND
if (keep_test(takeA, role)) {
cross[cx++] = rekey ? (((u64)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 < 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 < 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 == K - 3) {
u32 *dirt = (u32 *)ABtop; ABtop += (cx + 3) / 2 + 1;
ctr_solve(cross, cx, dirt);
} else {
u64 *cres = ABtop; ABtop += cx;
solve(cross, cx, d + 1, cres, 0);
}
ABtop = save;
}
static u32 *g_ord, *g_x;
static void sort0(int n) {
const u32 *X0 = XD[0];
u32 *cnt = g_x;
memset(cnt, 0, sizeof(u32) * (size_t)(n + 1));
for (int i = 0; i < n; i++) cnt[X0[i] + 1]++;
for (int i = 0; i < n; i++) cnt[i + 1] += cnt[i];
for (int i = 0; i < n; i++) g_ord[cnt[X0[i]]++] = (u32)i;
}
void run(int n, const u32 **x, int k, u32 *out) {
LFSH = 0; LFD = 0;
K = k;
RSZN = n;
for (int d = 0; d < k; d++) XD[d] = x[d];
OUT = out;
ABtop = AB;
if (n <= 1) return;
FUSE_OK = 0;
if (K >= 4 && n > 2 && (u32)n < CT_MAXV) {
const u32 *XL = XD[K - 1];
memset(CT0, 0, (size_t)n);
int ok = 1;
for (int i = 0; i < n; i++) if (++CT0[XL[i]] > 127) { ok = 0; break; }
memset(CT0, 0, (size_t)n);
FUSE_OK = ok;
}
{
int bit = 0;
while ((n >> bit) >= LEAFNB && bit < 30) bit++;
LFSH = bit; LFD = LF_D;
LFB = (int)((u32)n >> (LFSH + LFD)) + 1;
if (LFB > LEAFNB + 1) LFB = LEAFNB + 1;
}
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++;
for (int t = g; t < g1; t++) ev[m++] = mk(X1[g_ord[t]], 1, (int)g_ord[t]);
for (int t = g; t < g1; t++) ev[m++] = mk(X1[g_ord[t]], 0, (int)g_ord[t]);
g = g1;
}
ABtop = AB + m;
u64 *res = ABtop; ABtop += m;
solve(ev, m, 0, res, 0);
}
} // namespace kg
static u64 g_arena[KG_ARENA];
static u32 g_ord_st[300005], g_x_st[300005];
static int g_ready = 0;
static void setup() {
if (g_ready) return;
g_ready = 1;
kg::AB = g_arena;
kg::g_ord = g_ord_st;
kg::g_x = g_x_st;
}
void count_4d(int n, const unsigned *x[4], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run(n, (const unsigned **)x, 4, out);
}
void count_5d(int n, const unsigned *x[5], unsigned *out) {
setup();
memset(out, 0, sizeof(unsigned) * (size_t)n);
kg::run(n, (const unsigned **)x, 5, out);
}
#pragma GCC pop_options
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 86.772 ms | 9 MB + 468 KB | Accepted | Score: 100 | 显示更多 |