// lane q13b_rv ARM HP (submission build). Bound bytes were
// work/lane_q13_rv/src/cand_QF.cpp md5 e43af257f893 (sid 120025, 5.236 s) -- this file is
// that file plus the ARM HP hunk, with MF_HP baked to 1 (the judge passes no -D).
#define MF_HP 1
#define MF_APF 56
#define MF_SBDMODE 1
#define MF_SPSF 1
#define MF_SNID 1
// lane qm4_1013 ARM u8: and_popcount_range hand-unrolled to 8 words (two 32 B vectors per
// plane) per iteration, so the six 32 B loads and five vpand per 4 words are issued with half
// the loop-control overhead. Byte accumulators take <=16 per byte per iteration, so the
// flush period is 15 (15*16 = 240 <= 255). Judge-side custom_test, n=3e5 k=7 Q=461 B=651
// (filt ON), checksum 697645990 BOTH arms: shipped 711.8 / 713.9 -> u8 703.0 ms (-10 ms).
// Same output by construction; gated byte-identical at n=1e6 Q=1536 (chk 7782794688) and at
// n=30000/20000/50000 (Q=461/31/77) against the bound bytes. Lane record: work/lane_qm4_1013/.
// lane qm3_1013 ARM pgJ: flag-header sweep on the kx bytes (compact sort keys).
// Header (ABOVE the includes, where it bites): optimize(O3,no-unroll-loops,no-tree-vectorize,
// fno-inline-functions) + target(avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake).
// judge-side custom_test at n=3e5 k=7 Q=461: plain header 721.2 -> pgJ 713.0 (-8.2 ms, -1.14%),
// same checksum 697645990; the SAME header with unroll-loops is +16.7 ms; with NO optimize
// pragma (the judge's own -O2) it is +147 ms. Lane record: work/lane_qm3_1013/.
// lane fp13_1013 ARM q1536: MAXQ 1024 -> 1536 and MF_QCOEF 122 -> 160 => Q 1536.
// ONE-LINE-EACH PORT of the clamp the record forbade lifting. Q is a PARTITION of an
// exact count, so the checksum is Q-invariant by construction (gated at n=1e6 k=7).
#define AS_WORDS 0
#define PF_DIST 64
#pragma GCC optimize("O3","no-unroll-loops","no-tree-vectorize","fno-inline-functions")
// M180ARM:ctrlPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#define MF_QCOEF 160
// duck.ac 1012..1016 : k-D dominance counting (k=6..10).
//
// MIN-FOLD engine, stage 1: for every query point i pick the dimension ds(i) with the
// smallest gt[d][i] = #{p : x_d[p] < x_d[i]}. Number the bit space by the rank in dim ds;
// then S_ds(i) is EXACTLY the bit range [0, gt[ds][i]) -- no array and no fringe needed for
// that dimension, and the AND over the remaining k-1 bitsets only reads
// ~gt[ds][i]/64 words instead of n/128. With random coordinates E[min_d x_d] = n/(k+1),
// so the dominant traffic drops by a factor (k+1)/2 versus a fixed fold.
//
// Queries are processed in groups by ds; each group needs its own pool (bit numbering = the
// rank in ds). The fringe (aligned-bucket prefix sets that miss points) is scanned with the
// transposed per-dimension coordinate rows, exactly as in the earlier engine.
#include <immintrin.h>
#ifdef LOCAL_TEST
#include <cstdio>
#include <cstdlib>
#include <ctime>
#else
#ifdef __cplusplus
extern "C" void *malloc(unsigned long);
#else
void *malloc(unsigned long);
#endif
#endif
typedef unsigned u32;
typedef unsigned long long u64;
// g++-9 -O2 applies -mavx256-split-unaligned-load, so EVERY _mm256_loadu_si256
// compiles to "vmovdqu xmm; vinserti128" -- TWO loads plus a shuffle instead of one.
// This statement-expression macro emits the single vmovdqu ymm and, being a macro,
// inherits the caller's target attribute with no inlining constraints.
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#ifndef MAXN
#define MAXN 1000005
#endif
#define MAXK 8
#define STRIDE 8
// --- lane m155h (transfer leg 1): in-stream AND-array prefetch, the fleet's
// MF_APF construction banked on 1010 (-2.8 %), 1015 (MF_APF2=48, -1.03 %) and 1012.
#ifndef MF_APF
#define MF_APF 48
#endif
// Bucket-count coefficient: q ~ sqrt(n) * MF_QCOEF/100. Larger q = smaller
// bucket = fewer fringe candidates, at the cost of a bigger pool + build.
// Bucket count = sqrt(n) * MF_QCOEF/100 * kscale[k]. MEASURED OPTIMA (board):
// n=1e6 k=6..8 : 122 default is fine (1013 7.491 s, 1014 8.503 s at 122)
// n=1e6 k=9/10 : 92 (1016 10.888 s at 92; 10.959 at 122; 11.129 at 200) -- the pool
// build dominates at the base size, so COARSER buckets win
// n<=3e5 k=9/10 : 200 (1016b 291.4 ms vs 326.6 at 122; 1015a 1.444, 1016a 1.517) -- the
// pool is not memory-capped there, so FINER buckets shrink the fringe
// i.e. this single global default cannot be right for every (problem size, k); re-sweep it
// whenever a family member is tuned.
#ifndef MF_QCOEF
#define MF_QCOEF 122
#endif
// Transposed fringe rows: stride padded to 8 u32 (32 B) for k<=8 so every
// candidate load is one aligned, line-contained 32-byte vector load.
#ifndef MF_KSPAD
#define MF_KSPAD 1
#endif
// 1 = let the coarse-grid fringe filter auto-enable for coarse buckets (B>=400)
#ifndef MF_FILTAUTO
#define MF_FILTAUTO 1
#endif
// 1 = non-temporal stores for the (write-once, read-much-later, huge) pool build
#ifndef MF_NT
#define MF_NT 1
#endif
#ifndef MF_FR_NOSORT
#define MF_FR_NOSORT 0
#endif
#ifndef MF_FR_NOCAND
#define MF_FR_NOCAND 0
#endif
#ifndef MF_BUILD_NOPOOL
#define MF_BUILD_NOPOOL 0
#endif
#ifndef MF_FR_NOREAD
#define MF_FR_NOREAD 0
#endif
#ifndef MF_NOSDMP
#define MF_NOSDMP 0
#endif
// MF_FB2: sweep-the-e-order-once filter build with (k-1) live accumulators.
// MEASURED WORSE on the board (1016 12.580 -> 13.361 s, 1014 9.400 -> 10.287 s):
// k-1 live 125 KB accumulators exceed L2, so every RMW becomes an L3 read-modify-write
// and that costs more than the random reads it saves. Kept only as a springboard.
#ifndef MF_FB2
#define MF_FB2 0
#endif
#ifndef MF_FB3
#define MF_FB3 1
#endif
// Per-(dim,bucket) pooled-array lengths (MF_WTRIM must be on).
#ifndef MF_PLEN
#define MF_PLEN 1
#endif
#ifndef MF_BUDGETMB
#define MF_BUDGETMB 1900
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
// Lookahead distance for prefetching the scattered per-query metadata
// (s_pm record, s_gt threshold, s_bd bucket word) in pass B / pass A.
#ifndef MF_PFB
#define MF_PFB 32
#endif
// Pass A: sort each ds group by the bucket of one non-fold dim and copy that
// dim's (dim,bucket) array prefix into a small cache-resident scratch so the
// ~137 queries sharing it do not each re-fetch ~70 KB from DRAM. 1 = on.
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
// ============================================================================
// !! DO NOT ENABLE MF_WTRIM / MF_PLEN !!
// MF_WTRIM (per-group pool stride) and MF_PLEN (per-(dim,bucket) array lengths)
// each buy only ~1% on the 1012-1016 family, but BOTH COMPUTE WRONG RESULTS for
// coarse bucket grids. Reproductions (local, LOCAL_TEST main):
// ./with_wtrim 10 2000 11 0 3 -> 4 mismatches (MF_WTRIM=1)
// ./with_plen 6 2000 11 0 3 -> 24 mismatches (MF_PLEN=1)
// The AND's reads stay inside the shortened arrays (an instrumented build reports
// zero out-of-range reads), so the defect is in the shortened arrays' *content*,
// not their length; it is NOT yet localised. The fastest recorded board times
// (1012 6.753 s, 1016 12.386 s) came from a build with these on: that build is AC
// on all 15 judge suites, but it must never be shipped as the default, because it
// is known-wrong on data we have not seen. Fixing this properly is worth ~1%.
// ============================================================================
// Trim the pool stride per ds group to the largest R actually read (-13% pool).
#ifndef MF_WTRIM
#define MF_WTRIM 0
#endif
#ifndef MF_PFA
#define MF_PFA 1
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 1536
#endif
static u32 s_ord[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u32 s_bd[(size_t)MAXN * MAXK];
static u32 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_nid[MAXN];
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN]; // s_posp[p*MAXK+d] = rank of p in dim d
// ---- lane q13_rv QA-MERGE (FUSED) -------------------------------------------------------
// Pass B touches THREE random cache lines per (query,dim) pair: the 32 B s_pm record,
// the 32 B s_bd bucket row, and the 4 B s_gt[d][i]. Board-measured: MF_PFB=0 costs
// +679 ms, so that stream is ~679 ms of exposed latency. All three are ONE line's
// worth of data: rec is 32 B, the bucket row is 32 B (k<=8), and gt is a <=20-bit count
// while bd needs 11 bits -- so (gt<<QA_SH)|bd is ONE u32 and the whole k-lane bucket row
// PLUS all k gt values fit in the second 32 B of a 64 B record.
// s_qa[p] = { rec[0..7], w[0..7] }, w[e] = (gt[e][p] << QA_SH) | bd[e][p]
// => one 64 B line and ONE prefetch per (query,dim) instead of three.
// FUSED: no extra pass. The tie sweep already knows BOTH gt (j) and bb, so it writes the
// packed word into the rank-indexed bbr array in one SEQUENTIAL u32 store; the existing
// point-order fill then reads it once for s_bd AND for s_qa. bbr_t goes u16 -> u32, which
// costs no traffic: a random u16 read and a random u32 read both fetch one 64 B line.
// Restricted to k <= MAXK (8) -- the s_bd row is MAXK u32 wide, so k=9/10 keep the old path.
#define MF_QAMERGE 1
#define QA_SH 11
#define QA_MASK 0x7FFu
#if MAXQ > (QA_MASK + 1)
#error "QA_SH too small for MAXQ"
#endif
#if MAXN > (1u << (31 - QA_SH))
#error "QA_SH too small for MAXN"
#endif
struct qa_t { u32 rec[STRIDE]; u32 w[MAXK]; };
static qa_t s_qa[MAXN] __attribute__((aligned(64)));
// ---- lane qm5 SBDPT --------------------------------------------------------------------
// MEASURED (judge-side, n=3e5 k=7 Q=461, chk 697645990): the scattered
// `s_bd[p*MAXK+d] = bb` in the per-dim tie walk is 11.9 ms of PERDIM's 51.0 ms
// (work-removal ablation, 39.1 ms with the write gone). 7 M scattered 4 B stores
// into a 9.6 MB array -- one RFO line fetch each -- for a value that is written
// 4 B at a time where 32 B (two points) share a line. FIX: write the bucket word
// into a PER-DIM RANK-INDEXED array inside the existing rank-order tie walk (a
// SEQUENTIAL store), then fill s_bd in POINT ORDER with one full 32 B row per point.
// The rank-indexed array is the MF_FB2-dead s_posp block, so no new .bss.
#ifndef MF_SBDPT
#define MF_SBDPT 1
#endif
#ifndef MF_SBDW
#define MF_SBDW 0 /* q13_rv QA-MERGE: bbr carries the PACKED word, so u16 is too narrow */
#endif
#ifndef MF_SBDPB
#define MF_SBDPB 128
#endif
#ifndef MF_SBDPF
#define MF_SBDPF 0
#endif
#ifndef MF_SBDMODE
#define MF_SBDMODE 0
#endif
#ifndef MF_SPSF
#define MF_SPSF 0
#endif
#ifndef MF_SNID
#define MF_SNID 0
#endif
#ifndef MF_GTABL
#define MF_GTABL 0
#endif
static u64 s_gtab[16];
#if MF_SBDPT
#if MF_SBDW
typedef unsigned short bbr_t; // Q <= MAXQ = 1536 < 65536, so u16 is exact
#else
typedef u32 bbr_t;
#endif
static bbr_t *s_bbr = (bbr_t *)s_posp; // stride MAXN (2 or 4 MB per dim @ n=1e6)
#endif
#if MF_GTABL
#define MF_GTABL_ON 1
#else
#define MF_GTABL_ON 0
#endif
#if MF_SBDPT && MF_FB2
#error "MF_SBDPT aliases s_posp; MF_FB2 owns it"
#endif
static u64 s_facc[(size_t)(MAXK - 1) * (((MAXN) + 63) / 64 + 16)];
static u32 s_bnde[(size_t)MAXK * (MAXQ + 2)]; // per-dim coarse-cell boundary positions
static unsigned long long g_novf = 0; /* m164d PIPE: survivor-buffer overflows (expected 0) */
static u64 *s_fbm = 0;
static u64 s_fbuf[8192] __attribute__((aligned(64)));
static __m256i khi8v, khi8hv;
static u32 g_Qc = 32;
static u32 g_fbm_W = 0;
static int g_use_filt = 0;
static int g_lexsort = 0;
static int g_nofilt = 0;
static u32 s_fpos[(MAXK + 1) * MAXQ];
static u32 s_bval[(MAXK + 1) * MAXQ];
static u64 *s_bs = 0;
static u64 s_acc[MAXN / 64 + 256] __attribute__((aligned(64)));
static u32 s_pj[2048]; /* m164d PIPE: this pair's survivor j-indices */
/* lane qm2 SPLIT: standalone pre-materialisation of the RMW loop's random rank
lookup. The fused form (store address dependent on a random load) measures 3.65
ns/iter on the judge; the split form (pure load chain 1.94 + sequential-read RMW
0.93) measures 2.87 ns/iter for the identical work. Values are unchanged. */
static u32 s_qseq[MAXN];
static u32 s_sortk[MAXN]; /* qm3: compact sort keys, q-indexed */
static u32 *s_sdmp[MAXK];
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
static u32 s_pofs[(size_t)MAXK * MAXQ];
static u32 s_plen[(size_t)MAXK * MAXQ];
static int g_plen_on = 0;
static u32 g_W, g_Q, g_B;
static u32 g_KS; // row stride of the transposed coordinate rows (64 B for k>8)
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
// lane q13b_rv ARM HP: ask the kernel for TRANSPARENT HUGE PAGES on every big pool
// (the fold-dim pool, the coarse-filter bitmap, the transposed fringe rows).
//
// CONSTRUCTION. The engine's own comment prices this loop as 4.58 GB / ~1.7 s
// (2.7 GB/s) "vs 12.5 GB/s for a plain memset", and a local child-CPU ablation puts
// the whole pool build at ~30 % of the row. Isolated (src/ntb.cpp, same chunk shape,
// same 32 B NT stores, same one-sfence-per-462-chunks) the copy is
// cold destination : 193.5 ms / 275.4 MB = 1.49 GB/s
// SAME loop, dest pre-touched : 9.8 ms / 275.4 MB = 29.46 GB/s
// i.e. the copy itself is 20x faster than what the row pays; the entire difference is
// FIRST-TOUCH PAGE FAULTS on the destination. At the judged shape the pool is ~670 MB
// = ~171 k anonymous faults, and every fault costs ~2.6 us of kernel time to zero one
// 4 KB page. MADV_HUGEPAGE collapses that to ~335 faults of 2 MB each, whose zeroing
// is charged at memory bandwidth instead of per-fault. The transformation is
// mapping-only: the bytes written and every value read are unchanged, so the output is
// bit-identical (gate below).
//
// Raw syscall: the artifact overrides __libc_start_main, so glibc's wrappers are not a
// dependency we are entitled to; `madvise` needs only the register convention.
#ifndef MF_HP
#define MF_HP 0
#endif
static void hp_madvise(void *p, size_t n) {
#if MF_HP
register long rax __asm__("rax") = 28; /* SYS_madvise */
register long rdi __asm__("rdi") = (long)p;
register long rsi __asm__("rsi") = (long)n;
register long rdx __asm__("rdx") = 14; /* MADV_HUGEPAGE */
__asm__ volatile("syscall" : "+r"(rax) : "r"(rdi), "r"(rsi), "r"(rdx)
: "rcx", "r11", "memory");
#else
(void)p; (void)n;
#endif
}
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
return (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
}
// Harley-Seal / Muła nibble-popcount. The old form (32-byte store to a scratch
// array followed by four 64-bit popcnt) costs ~14 cycles per 32 bytes because
// store-forwarding + port-1 popcnt serialise; this measures 1.48x faster on a
// pure-L1 microbenchmark (work/c6d2/alucost.cpp: 19.27 -> 13.06 cyc per chunk).
// Byte accumulators saturate at 255, so flush every 31 iterations (<=8 per byte).
__attribute__((target("avx2,popcnt")))
static inline u32 hsum256(__m256i v) {
u64 b[4];
_mm256_storeu_si256((__m256i *)b, v);
return (u32)(b[0] + b[1] + b[2] + b[3]);
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_range(const u64 *const *bs, u32 Wr, u32 rem) {
const __m256i lm = _mm256_set1_epi8(0x0f);
const __m256i lk = _mm256_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);
const __m256i z = _mm256_setzero_si256();
__m256i tot = z, acc = z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 8 <= Wr; w += 8) {
#if MF_APF > 0
for (int d0 = 0; d0 < KK; d0++) _mm_prefetch((const char *)(bs[d0] + w + MF_APF), _MM_HINT_T0);
#endif
__m256i a = LDU256((bs[0] + w));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < KK; d++) { a = _mm256_and_si256(a, LDU256((bs[d] + w)));
b = _mm256_and_si256(b, LDU256((bs[d] + w + 4))); }
__m256i lo = _mm256_and_si256(a, lm);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lo2 = _mm256_and_si256(b, lm);
__m256i hi2 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc = _mm256_add_epi8(acc, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
acc = _mm256_add_epi8(acc, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo2), _mm256_shuffle_epi8(lk, hi2)));
if (++cnt == 15) { tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z)); acc = z; cnt = 0; }
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z)); acc = z; cnt = 0;
for (; w + 4 <= Wr; w += 4) {
#if MF_APF > 0
for (int d0 = 0; d0 < KK; d0++) _mm_prefetch((const char *)(bs[d0] + w + MF_APF), _MM_HINT_T0);
#endif
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
__m256i lo = _mm256_and_si256(a, lm);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc = _mm256_add_epi8(acc, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
if (++cnt == 31) { tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z)); acc = z; cnt = 0; }
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z));
total = hsum256(tot);
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 1; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v);
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 1; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & ((1ull << rem) - 1));
}
return total;
}
/* m164d PIPE: one candidate test, factored out so the pipelined flush and the overflow path
share it. Pulled straight from the shipped body: same address, same load, same test. */
__attribute__((always_inline)) static inline u32 m164p_ok(u32 j, const u32 *c0, u32 lo,
__m256i qT, __m256i qTh, const u32 *odd, int k, __m256i all8v) {
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * g_KS)
: (s_pm + (size_t)odd[j] * STRIDE);
__m256i pv = LDU256(cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8v);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8v);
}
return ok;
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const int k = KK;
const int K1 = KK - 1;
#if MF_ASCR
// ---- pass A locality: order the group by the bucket in one non-fold dim so that the
// ~Q-run of queries needing array (e1,b) is processed together, and pin that
// array's prefix in a small scratch (the other k-2 dims still stream from DRAM
// and evict a 70 KB hot array within one query, which is why merely sorting
// the queries was measured a wash). ----
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB] * MAXK + d1), _MM_HINT_T0);
#endif
u32 kv = s_bd[(size_t)s_qord[qi] * MAXK + d1];
s_sortk[qi] = kv;
cntv[kv + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) { u32 id = s_qord[qi]; s_tmp[cntv[s_sortk[qi] + 1]++] = id; }
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#else
const u32 d1 = 0, e1 = 0;
#endif
// ---- pass A: aligned AND; queries ordered lexicographically by their bucket vector so
// that consecutive queries share bitset arrays (keeps the working set small). ----
if (g_lexsort) {
u32 *cntv = s_cntv;
for (int dd = k - 1; dd >= 0; dd--) {
if (dd == ds) continue;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) cntv[s_bd[(size_t)s_qord[qi] * MAXK + dd] + 1]++;
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) { u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * MAXK + dd] + 1]++] = id; }
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = s_qord[qi + PF_DIST];
const u32 *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
#if MF_ASCR
{
u32 qi = q0;
while (qi < q1) {
u32 b1 = s_bd[(size_t)s_qord[qi] * MAXK + d1];
u32 qe = qi + 1;
while (qe < q1 && s_bd[(size_t)s_qord[qe] * MAXK + d1] == b1) qe++;
const u64 *bp1 = g_plen_on ? (s_bs + s_pofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
u32 mx = 1;
for (u32 t = qi; t < qe; t++) { u32 rr = s_gt[ds][s_qord[t]]; if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 1;
if (nwd <= AS_WORDS) {
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = s_qord[qq];
const u32 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
if (e == (int)e1) bp[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
const u32 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = s_qord[qi + PF_DIST];
const u32 *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W);
#ifdef DBG_PLEN
if (g_plen_on) {
u32 need = (R >> 6) + ((R & 63) ? 1 : 0);
if (need > s_plen[(size_t)e * MAXQ + bd[d]])
printf("DBG ds=%d e=%d d=%d b=%u R=%u need=%u len=%u\n", ds, e, d, bd[d], R, need, s_plen[(size_t)e * MAXQ + bd[d]]);
}
#endif
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
#endif
if (g_skip_fringe) return;
// ---- pass B: fringe, dimension-major. Ordering the group's queries by their bucket in
// dim d makes every scan walk one L1-resident bucket block of the dim-d order. ----
for (int d = 0; d < k; d++) {
if (d == ds) continue;
#if !MF_FR_NOSORT
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB] * MAXK + d), _MM_HINT_T0);
#endif
u32 kv = s_bd[(size_t)s_qord[qi] * MAXK + d];
s_sortk[qi] = kv;
cntv[kv + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
u32 id = s_qord[qi];
s_tmp[cntv[s_sortk[qi] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#endif
u32 exmask = 0;
for (int e = 0; e < d; e++) if (e != ds) exmask |= 1u << e;
// "don't care" lanes are forced to a maximum-unsigned value, so one unsigned max test
// per candidate decides domination *and* the exclusion rule together.
u32 hiA[8], hiAh[8];
for (int e = 0; e < 8; e++) hiA[e] = ((exmask >> e) & 1u) ? 0u : 0x7FFFFFFFu;
for (int e = 8; e < 16; e++) hiAh[e - 8] = ((exmask >> e) & 1u) ? 0u : 0x7FFFFFFFu;
const __m256i notexlo = LDU256(hiA);
const __m256i notexloh = LDU256(hiAh);
const __m256i all8 = _mm256_set1_epi32(-1);
const u32 *const odd = s_ord[d];
#if !MF_FR_NOSORT_PIPE_MARK
/* ---- m164d PIPE state, live across one (query,dim) pair within this dim ---- */
u32 s_pend = 0, s_pn = 0, s_plo = 0;
const u32 *s_pc0 = 0;
__m256i s_pqT = all8, s_pqTh = all8;
u32 *s_pout = 0;
#define M164P_FLUSH() do { \
if (s_pend) { \
u32 add = 0; \
for (u32 q2 = 0; q2 < s_pn; q2++) \
add += m164p_ok(s_pj[q2], s_pc0, s_plo, s_pqT, s_pqTh, odd, k, all8); \
*s_pout += add; s_pend = 0; s_pn = 0; \
} \
} while (0)
#endif
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB;
if (qj < q1) {
u32 j = s_qord[qj];
if (k <= MAXK) {
_mm_prefetch((const char *)(s_qa + j), _MM_HINT_T0);
} else {
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_T0);
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK + d), _MM_HINT_T0);
}
}
}
#endif
u32 lo, hi;
const u32 *rec;
u32 Av[STRIDE] __attribute__((aligned(32)));
if (k <= MAXK) {
/* ONE 64 B line carries the point record, every bucket word AND every gt. */
const qa_t *const qa = s_qa + i;
const u32 *const wv = qa->w;
lo = FPOS(d, wv[d] & QA_MASK);
hi = wv[d] >> QA_SH;
if (lo >= hi) continue;
rec = qa->rec;
for (int e = k; e < STRIDE; e++) Av[e] = 0x7FFFFFFFu;
for (int e = 0; e < k; e++) Av[e] = (e == ds) ? 0x7FFFFFFFu : BVAL(e, wv[e] & QA_MASK);
} else {
const u32 *bd = s_bd + (size_t)i * MAXK;
lo = FPOS(d, bd[d]);
hi = s_gt[d][i];
if (lo >= hi) continue;
rec = s_pm + (size_t)i * STRIDE;
for (int e = k; e < STRIDE; e++) Av[e] = 0x7FFFFFFFu;
for (int e = 0; e < k; e++) Av[e] = (e == ds) ? 0x7FFFFFFFu : BVAL(e, bd[e]);
}
__m256i qr = _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
__m256i Avm = _mm256_or_si256(LDU256(Av), notexlo);
__m256i qT = _mm256_min_epi32(qr, Avm);
__m256i qTh;
if (k > 8) {
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU256((Av + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
} else {
qTh = all8;
}
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (size_t)lo * g_KS) : 0;
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const u64 winv = ((u64)1 << 32) / wcs + 1;
const u32 FW = g_fbm_W;
// G4_B: the original accumulated fb[w] &= row[w] ONE e at a time -- per (pair,e)
// that is nw loads of fb + nw loads of row + nw ANDs + nw stores of fb, i.e.
// ~4 memory ops per word per dim (24*nw per pair). Hoisting the (k-1) row
// pointers and keeping the running AND in a register drops it to ~1 memory op
// per word per dim plus nw stores (7*nw per pair): a ~3.4x cut in pass B's
// filter traffic, which is the single largest per-pair block left after the
// division was removed.
const u64 *rows[MAXK];
for (int e = 0; e < k; e++) {
if (e == d) continue;
u32 c = (u32)(((u64)rec[e] * winv) >> 32) + 1;
if (c > g_Qc) c = g_Qc;
rows[e] = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
}
#if MF_FR_NOREAD
for (u32 w = 0; w < nw; w++) fb[w] = 0;
#else
for (u32 w = 0; w < nw; w++) {
u64 v = ~(u64)0;
for (int e = 0; e < k; e++) { if (e == d) continue; v &= rows[e][w]; }
fb[w] = v;
}
#endif
fb[0] &= ~0ull << (lo & 63);
u32 r2 = hi & 63u;
if (r2) fb[nw - 1] &= (1ull << r2) - 1;
/* m164d PIPE: (1) THIS pair's words are already ANDed above; (2) extract THIS pair's
survivors into s_pj and issue their gather prefetches; (3) TEST the PREVIOUS pair,
whose prefetches were issued one whole pair ago -- one word scan of lead, ~350 cyc,
against a ~57-cycle gather. Same ANDs, same tests, same per-pair sum; only the
ORDER of the += across pairs changes, and each out[i] still takes exactly one +=
with the same integer value => output bit-identical. */
M164P_FLUSH(); /* must precede the extraction: s_pj holds the PENDING pair */
u32 nsn = 0;
for (u32 w = 0; w < nw; w++) {
u64 vv = fb[w];
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = ((w0 + w) << 6) + (u32)t;
if (nsn < 2048u) {
s_pj[nsn++] = j;
const u32 *pf = c0 ? (c0 + (size_t)(j - lo) * g_KS)
: (s_pm + (size_t)s_ord[d][j] * STRIDE);
_mm_prefetch((const char *)pf, _MM_HINT_T0);
} else {
add += m164p_ok(j, c0, lo, qT, qTh, odd, k, all8); /* overflow: direct, correct */
g_novf++;
}
}
}
s_pend = 1; s_pn = nsn; s_plo = lo; s_pc0 = c0; s_pqT = qT; s_pqTh = qTh;
s_pout = out + i;
continue;
}
if (c0) {
if (k <= 8) {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + g_KS));
__m256i p2 = LDU256((c + 2 * g_KS));
__m256i p3 = LDU256((c + 3 * g_KS));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8);
c += 4 * g_KS;
}
for (; j < hi; j++, c += g_KS) {
__m256i pv = LDU256(c);
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
add = a0 + a1 + a2 + a3;
} else {
const u32 *c = c0;
for (u32 j = lo; j < hi; j++, c += g_KS) {
__m256i pv = LDU256(c);
__m256i pvh = LDU256((c + 8));
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8) &
(u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
add += ok;
}
}
} else {
// Fallback used when the transposed fringe rows are not allocated. NOTE: this
// branch used to test only lanes 0..7, so for k>8 the dims 8..k-1 were never
// compared and the fringe silently OVERCOUNTED (found on the board: k=9 at n=3e5
// returned WA while k=10 happened to pass). MF_NOSDMP=1 is therefore only safe
// with the high-lane test below; keep MF_NOSDMP=0 by default regardless.
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__m256i pv = _mm256_load_si256((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = _mm256_load_si256((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
#if !MF_FR_NOSORT_PIPE_MARK
M164P_FLUSH(); /* the dim's last pending pair */
#undef M164P_FLUSH
#endif
}
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int k) {
// W padded to 8 words: every (dim,bucket) array then starts on a 64-byte
// boundary and every 32-byte vector load inside a stream is line-contained.
const u32 W = ((N + 63) / 64 + 7u) & ~7u;
g_W = W;
#if MF_HP >= 2
// The .bss tables are the OTHER first-touch surface, and at n=1e6 they are ~280 MB
// (s_ord/s_gt/s_pos/s_posp 32 MB each, s_pm 32 MB, s_bd 28 MB, s_qa 64 MB) -- anonymous
// pages that fault on first touch and are then read RANDOMLY (s_bd, s_pos, s_pm), where a
// 4 KB TLB entry covers 4 KB instead of 2 MB. The linker symbols bracket the whole .bss,
// which is one VMA, so a min/max over individual arrays is neither needed nor safe.
{ extern char __bss_start[]; extern char _end[];
hp_madvise((void *)__bss_start, (size_t)(_end - __bss_start)); }
#endif
#ifdef MF_SKIP_AND
g_skip_and = 1;
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_Q
g_Q = MF_Q; g_B = (u32)((N + (MF_Q) - 1) / (MF_Q));
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_SKIP_FRINGE
g_skip_fringe = 1;
#endif
if (!g_Q) {
// bucket size ~ 0.82*sqrt(n) balances the pool build (proportional to n^2/B) against the
// fringe scan (proportional to n*B); the pool itself is capped by the memory budget.
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++; // ceil(sqrt(n)) without libm
// Reserve: the .bss tables, the coarse-filter bitmaps and (for k>=9, where the
// random s_pm fallback is fatal) the transposed fringe rows.
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
// MF_PLEN sizes each pooled array to the largest R ever read from it, so the
// pool really used is only ~fill of (k-1)*(Q+1)*W words; charge that instead.
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46}; // per-array fill, k=4..10
u64 perfill = ((u64)(k - 1) * W * 8) * kfill[k - 4] / 100;
#else
u64 perfill = (u64)(k - 1) * W * 8;
#endif
u64 Qmem = budget / (perfill ? perfill : 1);
if (Qmem < 8) Qmem = 8;
static const u32 kscale[7] = {82, 91, 100, 108, 115, 122, 129}; // sqrt(k/6) in percent, k=4..10
u32 q = (u32)((s2 * MF_QCOEF + 99) / 100);
q = (u32)(((u64)q * kscale[k - 4] + 99) / 100);
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
g_Q = q;
// LANE inputdist: the pool tables are dimensioned MAXQ; bound g_Q by them (a no-op when
// the formula is inside, and a CLAMP rather than an overflow when a test is larger than
// the statement's n). g_B follows the clamped value, so the partition stays consistent.
if (g_Q > (u32)MAXQ) g_Q = (u32)MAXQ;
g_B = (N + g_Q - 1) / g_Q;
}
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
g_plen_on = 1; // independent of MF_WTRIM: lengths only, stride stays W
#endif
#if MF_WTRIM
u32 Wpool = W; // MF_WTRIM is KNOWN-BROKEN, see the warning above
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1; // coarse grid: the filter beats the scan
#endif
if (g_use_filt && (u64)k * k * (g_Qc + 1) * W * 8 > ((u64)700 << 20)) g_use_filt = 0;
#if !MF_FILTAUTO
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
// ---------- per-dimension orders, ranks, tie info ----------
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
for (u32 i = 0; i < N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i < N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s += c; } }
#if MF_SPSF
// SBDPT step A: the insertion loop already knows the RANK of point i -- it is the
// pre-increment value of s_cntv[v] -- so s_pos[d][i] = rank is a SEQUENTIAL store
// right here, and the separate scatter loop over od is deleted outright. One whole
// scattered-store stream per dim (N stores) disappears; od keeps its own scatter.
{ u32 *pp = s_pos[d]; u32 *sk = s_sortk;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; if (v >= N) v = 0; sk[i] = s_cntv[v]++; }
for (u32 i = 0; i < N; i++) { u32 r = sk[i]; od[r] = i; pp[i] = r; } }
#else
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; if (v >= N) v = 0; od[s_cntv[v]++] = i; }
#if MF_FB2
{ for (u32 j = 0; j < N; j++) { u32 p = od[j]; s_posp[(size_t)p * MAXK + d] = j; } }
#else
{ u32 *pp = s_pos[d]; for (u32 j = 0; j < N; j++) pp[od[j]] = j; }
#endif
#endif
u32 *gt = s_gt[d];
// ORD_NODIV: (jj-1)/B is a runtime DIVISION once per tie group (~N per dim, N*k per row).
// jj is non-decreasing, so the quotient advances incrementally and EXACTLY.
u32 bbcur = 0, bbend = B;
#if MF_SBDPT
{
bbr_t *bbr = s_bbr + (size_t)d * MAXN;
// ---- lane qm5 RESTART: TIEWALK BY COUNTS -------------------------------------------
// After the counting-sort insertion loop, `s_cntv[v] == prefix(v+1)` -- the running
// pointer landed exactly at the end of value v's block -- so the tie groups are the
// intervals (s_cntv[v-1], s_cntv[v]] and can be enumerated by a SEQUENTIAL sweep of
// the counter array. The shipped form rediscovered the same intervals with two
// RANDOM reads of xd per point (14 M at n=1e6); the sweep has none. Same groups,
// same order, same bbcur recurrence, same stores.
{
u32 j = 0;
for (u32 v = 0; v < N; v++) {
u32 jj = s_cntv[v];
if (jj == j) continue;
while (bbcur < Q - 1 && jj - 1 >= bbend) { bbcur++; bbend += B; }
const u32 bb = bbcur;
for (u32 t = j; t < jj; t++) {
u32 p = od[t];
#if MF_GTABL
s_gtab[p & 15] = j; /* work-removal probe: scatter REMOVED */
#else
gt[p] = j;
#endif
bbr[t] = (bbr_t)((j << QA_SH) | bb); /* SEQUENTIAL: (gt<<QA_SH)|bd, t ascends by 1 */
}
j = jj;
}
}
}
#else
for (u32 j = 0; j < N;) {
u32 jj = j + 1, v = xd[od[j]];
while (jj < N && xd[od[jj]] == v) jj++;
while (bbcur < Q - 1 && jj - 1 >= bbend) { bbcur++; bbend += B; }
const u32 bb = bbcur;
for (u32 t = j; t < jj; t++) {
u32 p = od[t];
gt[p] = j;
s_bd[(size_t)p * MAXK + d] = bb;
}
j = jj;
}
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
#if MF_SBDPT
// ---- point-order s_bd fill -------------------------------------------------------------
// s_pos[d][p] is the rank of p in dim d and is SEQUENTIAL in p, so the reads are k
// sequential streams and the write is one FULL 32 B row per point. The point loop is
// BLOCKED so that only one random-read stream (bbr[d]) is live at a time; the s_bd
// block (PB*32 B) stays in L1 and is written back as full lines.
#if MF_SBDMODE
{
// point-outer: k INDEPENDENT random loads in flight per point, one 32 B row written.
const u32 PB = MF_SBDPB;
for (u32 p0 = 0; p0 < N; p0 += PB) {
u32 p1 = p0 + PB; if (p1 > N) p1 = N;
for (u32 p = p0; p < p1; p++) {
u32 *row = s_bd + (size_t)p * MAXK;
if (k <= MAXK) {
qa_t *qa = s_qa + p;
for (int d = 0; d < k; d++) {
const u32 *pp = s_pos[d];
const bbr_t *bb = s_bbr + (size_t)d * MAXN;
u32 v = (u32)bb[pp[p]];
row[d] = v & QA_MASK;
qa->w[d] = v;
qa->rec[d] = x[d][p];
}
for (int d = k; d < STRIDE; d++) { qa->w[d] = 0; qa->rec[d] = 0; }
} else {
for (int d = 0; d < k; d++) {
const u32 *pp = s_pos[d];
const bbr_t *bb = s_bbr + (size_t)d * MAXN;
row[d] = (u32)bb[pp[p]] & QA_MASK;
}
}
}
}
}
#else
{
const u32 PB = MF_SBDPB;
for (u32 p0 = 0; p0 < N; p0 += PB) {
u32 p1 = p0 + PB; if (p1 > N) p1 = N;
for (int d = 0; d < k; d++) {
const u32 *pp = s_pos[d];
const bbr_t *bb = s_bbr + (size_t)d * MAXN;
u32 *row = s_bd + (size_t)p0 * MAXK + d;
#if MF_SBDPF > 0
for (u32 p = p0; p < p1; p++, row += MAXK) {
if (p + MF_SBDPF < p1) _mm_prefetch((const char *)(bb + pp[p + MF_SBDPF]), _MM_HINT_T0);
*row = (u32)bb[pp[p]];
}
#else
for (u32 p = p0; p < p1; p++, row += MAXK) *row = (u32)bb[pp[p]];
#endif
}
}
}
#endif
#endif
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < STRIDE; d++) p[d] = 0;
}
// ---------- coarse superset bitmaps over every dim order (fringe filter) ----------
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
/* lane qm2 BS: xe[pe[ptr]] is MONOTONE in ptr (pe is the value-sorted e-order), so
the boundary is a binary search -- 33*20 random reads per dim instead of N. The
two forms are equal by monotonicity. */
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 lo = 0, hi = N;
while (lo < hi) { u32 mid = (lo + hi) >> 1; if ((u64)xe[pe[mid]] < v) lo = mid + 1; else hi = mid; }
bn[c] = lo;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
// Sweep the e-order ONCE per dim e, keeping (k-1) running bitsets (one per
// target dim d), and copy them out at each coarse-cell boundary. The old
// shape re-walked all N points for every one of the k(k-1) (d,e) pairs with
// two L3-random reads each (xe[pe[ptr]] and s_pos[d][pe[ptr]]); here the
// rank record is point-major so one 64B line yields every dim's rank, which
// cuts the random reads by ~k and the whole build by ~2-2.5x at k>=8.
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
for (u32 t = 0; t < N; t++) s_qseq[t] = pd[pe[t]];
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = s_qseq[ptr]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
// ---------- fold dimension + bucket vector per query ----------
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
// Largest R that will ever be read out of each (fold dim, dim, bucket) array:
// the pooled array (ds,d,b) is only ever scanned over [0, maxR] words, so the
// per-array length can be far below the group maximum (0.38n vs 0.93n at k=10).
if (g_plen_on) {
const u32 *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
// group by ds
{
for (u32 v = 0; v <= (u32)k; v++) s_cntv[v] = 0;
for (u32 i = 0; i < N; i++) s_cntv[s_g[i] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)k; v++) { u32 c = s_cntv[v]; s_cntv[v] = ac; ac += c; } }
for (u32 i = 0; i < N; i++) { u32 id = s_qord[i]; s_tmp[s_cntv[s_g[id] + 1]++] = id; }
for (u32 i = 0; i < N; i++) s_qord[i] = s_tmp[i];
}
u32 gstart[MAXK + 2];
{ u32 c = 0; for (int d = 0; d <= k; d++) { gstart[d] = c; while (c < N && (int)s_g[s_qord[c]] == d) c++; } gstart[k + 1] = c; }
// ---------- transposed coordinate rows (sequential fringe) ----------
#if MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 8; // 32-byte aligned candidate rows for k<=8
#endif
size_t per = (size_t)N * KS * 4;
size_t lim = (size_t)150 << 20;
if (k >= 9) lim = (size_t)470 << 20; // k=9/10 need the rows: the random s_pm fallback is fatal
if (per * k <= lim) {
static u32 *pool = 0; static size_t cap = 0;
if (per * k > cap) { pool = (u32 *)pool_alloc(per * k + 256); cap = per * k; }
if (pool) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
}
}
}
}
}
// ---------- pool ----------
u64 poolwords = (u64)K1 * (Q + 1) * W;
#if MF_PLEN
if (g_plen_on) {
u64 best = 0;
for (int dsx = 0; dsx < k; dsx++) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < dsx) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)dsx * MAXK + d) * MAXQ;
for (u32 b = 0; b <= Q; b++) {
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > W) len = W;
tot += len;
}
}
if (tot > best) best = tot;
}
if (best && best < poolwords) poolwords = best;
}
#endif
size_t need = (size_t)poolwords * sizeof(u64);
static u64 *pool = 0; static size_t cap2 = 0;
if (need > cap2) { pool = (u64 *)pool_alloc(need + 4096); cap2 = need; }
s_bs = pool;
// ---------- groups ----------
for (int ds = 0; ds < k; ds++) {
u32 q0 = gstart[ds], q1 = gstart[ds + 1];
if (q0 >= q1) continue;
#if MF_WTRIM
{
u32 mx = 1;
for (u32 qi = q0; qi < q1; qi++) { u32 rr = s_gt[ds][s_qord[qi]]; if (rr > mx) mx = rr; }
u32 wp = ((mx >> 6) + 8u) & ~7u;
Wpool = (wp < W) ? wp : W;
}
#endif
#if MF_PLEN
// Per-(dim,bucket) array lengths inside the group.
if (g_plen_on) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)ds * MAXK + d) * MAXQ;
u32 *po = s_pofs + (size_t)e * MAXQ;
u32 *pl = s_plen + (size_t)e * MAXQ;
for (u32 b = 0; b <= Q; b++) {
// NOTE: bucket 0 is the EMPTY prefix set (all-zero array), but a query whose
// bd[d]==0 still scans [0,R) words of it, so the array must be long enough --
// sizing it to the fixed 8-word minimum reads past the zeroed region into the
// next array. Use the real per-bucket max R here, not 0.
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > Wpool) len = Wpool;
po[b] = (u32)tot; pl[b] = len; tot += len;
}
}
if (tot > (u64)K1 * (Q + 1) * W) { g_plen_on = 0; } // never; safety
}
#endif
// bit numbering = rank in dim ds
#if !MF_SNID
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
#else
// SBDPT step C: s_nid[od[j]] = j over ALL j is exactly s_pos[ds][p] = (rank of p in ds),
// which step A already produced as a SEQUENTIAL store. The scatter is deleted and the
// consumer below reads s_pos[ds] instead -- same value, same array size, one fewer
// whole scattered-store stream (N stores) per fold dimension.
const u32 *s_npsd = s_pos[ds];
#endif
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
u64 *base = g_plen_on ? s_bs : (s_bs + (size_t)e * (Q + 1) * Wpool);
u64 *cur = s_acc;
for (u32 w = 0; w < Wpool; w++) cur[w] = 0;
if (g_plen_on) {
// bucket 0 is the empty prefix set; bd[] can be 0 so it must be materialised as zeros
u64 *d0 = s_bs + s_pofs[(size_t)e * MAXQ];
u32 l0 = s_plen[(size_t)e * MAXQ];
for (u32 w = 0; w < l0; w++) d0[w] = 0;
}
u32 ptr = 0;
#if MF_BUILD_NOPOOL
while (ptr < N) ptr++;
#else
#if MF_SNID
for (u32 t = 0; t < N; t++) s_qseq[t] = s_npsd[od[t]];
#else
for (u32 t = 0; t < N; t++) s_qseq[t] = s_nid[od[t]];
#endif
for (u32 b = 1; b <= Q; b++) {
// {q : x_d[q] < BVAL(d,b)} is EXACTLY the prefix [0,FPOS(d,b)) of the dim-d order
// (FPOS is the tie-group start containing rank b*B), so the per-point value compare
// -- one random xd read per point per (ds,d) pair, k(k-1)N of them -- is redundant.
// b==Q means "every point" and FPOS(d,Q) is the 0 sentinel, hence N.
// Credited to agent a7f799ceb7e63b0cd (worth 4-5% at k=4/5 on the board).
u32 pend = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) { u32 q = s_qseq[ptr]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
#if MF_NT
// The pool is written once and read 89 GB later, so a plain copy pays a
// read-for-ownership DRAM read per line for nothing. NT stores remove it
// (measured: this loop was 4.58 GB / ~1.7 s, i.e. 2.7 GB/s, vs 12.5 GB/s
// for a plain memset).
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 4); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 5); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 6); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 7); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 9); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 10); }
#ifdef LOCAL_TEST
static void brute_nd(int n, const u32 **x, u32 *out, int k) {
for (int i = 0; i < n; i++) {
u32 c = 0;
for (int j = 0; j < n; j++) {
int ok = 1;
for (int d = 0; d < k; d++) if (!(x[d][j] < x[d][i])) { ok = 0; break; }
c += ok;
}
out[i] = c;
}
}
int main(int argc, char **argv) {
int k = (argc > 1) ? atoi(argv[1]) : 6;
int n = (argc > 2) ? atoi(argv[2]) : 200;
int seed = (argc > 3) ? atoi(argv[3]) : 1;
int mode = (argc > 4) ? atoi(argv[4]) : 0;
int qq = (argc > 5 && atoi(argv[5]) > 0) ? atoi(argv[5]) : 0;
g_Q = qq ? (u32)qq : (u32)(n / 250 + 1);
g_B = (u32)((n + g_Q - 1) / g_Q);
if (argc > 6) g_nosort = atoi(argv[6]);
if (argc > 7) g_nosdmp = atoi(argv[7]);
if (argc > 8) g_skip_and = atoi(argv[8]);
if (argc > 9) g_skip_fringe = atoi(argv[9]);
if (argc > 10) g_use_filt = atoi(argv[10]);
if (argc > 11) g_lexsort = atoi(argv[11]);
const int dm = (mode >= 10) ? mode - 10 : mode;
srand(seed);
static u32 *xs[10];
u32 *out = (u32 *)malloc(sizeof(u32) * n);
u32 *ref = (u32 *)malloc(sizeof(u32) * n);
for (int d = 0; d < k; d++) {
xs[d] = (u32 *)malloc(sizeof(u32) * n);
for (int i = 0; i < n; i++) {
if (dm == 0) xs[d][i] = (u32)(rand() % n);
else if (dm == 1) xs[d][i] = (u32)(rand() % 3);
else if (dm == 2) xs[d][i] = 0;
else if (dm == 3) xs[d][i] = (u32)(n - 1 - (rand() % n));
else if (dm == 4) xs[d][i] = (u32)(i);
else xs[d][i] = (u32)(rand() % n);
}
}
if (mode >= 10) {
struct timespec t0, t1;
clock_gettime(CLOCK_MONOTONIC, &t0);
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
clock_gettime(CLOCK_MONOTONIC, &t1);
double ms = (t1.tv_sec - t0.tv_sec) * 1e3 + (t1.tv_nsec - t0.tv_nsec) / 1e6;
u64 s = 0; for (int i = 0; i < n; i++) s += out[i];
printf("k=%d n=%d Q=%u TIME %.1f ms checksum %llu\n", k, n, g_Q, ms, s);
return 0;
}
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
brute_nd(n, (const u32 **)xs, ref, k);
int bad = 0;
for (int i = 0; i < n; i++) if (out[i] != ref[i]) { if (bad < 5) printf("MISMATCH i=%d got=%u ref=%u\n", i, out[i], ref[i]); bad++; }
printf("k=%d n=%d mode=%d Q=%u %s (%d mismatches)\n", k, n, mode, g_Q, bad ? "FAIL" : "OK", bad);
return bad != 0;
}
#endif