/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140756 supplies the current eight-dimensional engine and its inherited source/license notices below. Our seven-dimensional https://duck.ac/submission/140865 supplies the positively measured non-temporal prefix-snapshot copy mechanism. No separate license was declared.
Idea: copy each byte-count snapshot with two256-bit streaming stores per64-byte destination cacheline. This engine's allocator aligns every snapshot pool to2MiB and every row offset/length is a multiple of64; use unaligned source loads so no mutable-counter layout change is needed. Fence once after all snapshots and before any consumer. Keep128-point blocks, all count bytes, query arithmetic, edge masks and geometry unchanged.
Purpose: reduce write allocation and displacement of live prefix counters during construction, independently of query-loop gating. Validate complete outputs, repeated-rank cases, brute-force samples and guards before official timing. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140704 supplies the current eight-dimensional engine, including all inherited source/license notices below. Our seven-dimensional https://duck.ac/submission/126214 supplies immediate SAD reduction into independent64-bit accumulators. No separate license was declared.
Idea: apply immediate widening only to this engine's selective prefix kernel: every nonempty256-bit intersection's nibble counts are immediately SAD-reduced and accumulated in64-bit lanes. Empty512-bit first-four intersections still skip all later loads and count work. Remove the31-iteration batch boundaries, accumulator clearing and reductions executed even for empty batches. Preserve both prefetch horizons, clamped tail, stream order and edge semantics.
Purpose: test whether simpler loop control and reductions proportional to surviving blocks improve the eight-dimensional early-rejection kernel. All counts remain exact without byte overflow. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140531 supplies the current engine and all inherited source/license notices below. No separate license was declared.
Idea: test the previously unmeasured nearby coefficient118, reducing the bitmap snapshot count from coefficient122 while retaining exact fringe correction. Existing physical metadata bounds remain sufficient because the new logical geometry is smaller. All rank, tie, edge, query and fallback semantics remain exact.
Purpose: locate the current construction-versus-fringe tradeoff on the coarser side after the retained byte/prefix and signature improvements; this concrete geometry has no prior result in this engine. Retain only if the official runtime improves. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140219 supplies the accepted eight-coordinate million-point engine and inherited references. Our https://duck.ac/submission/140142 used independent even/odd saturated-difference OR chains in the ten-coordinate300k scanner; adapt that exact reduction to this engine's seven live planes and cached query bytes. Public sources declared no separate license.
* Idea: accumulate the same seven unsigned saturated byte differences in two independent even/odd OR chains, then combine once before zero/one classification. Longest reduction depth decreases while every loaded coordinate, signature, ambiguity mask, scalar tail and exact replay remains identical. Apply to both complete and partial32-point blocks; no geometry, prefix or edge changes are included.
* Purpose: shorten the fringe Boolean dependency in this distinct million-point implementation after validating complete output equivalence.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/139881 supplies the current coalesced-query eight-dimensional engine and all inherited notices below. Our https://duck.ac/submission/140142 supplies exact capacity-bounded cache-line metadata rows in the ten-dimensional300k sibling. No separate license was declared; inherited credits remain intact.
* Idea: for every1<=N<=1000000, ceil(sqrt(N))<=1000 and the unchanged coefficient122/scale115 imply Q<=1403 before the memory cap. Every metadata index is in0..Q inclusive, so round1404 slots to1408 for64-byte-aligned rows. Replace the unused4096-slot physical stride only; runtime Q/B, sorting, bitmaps, signatures and all query geometry are unchanged.
* Purpose: remove unreachable row capacity from metadata initialization and working sets, preserving the exact algorithm and testing the larger-scale locality effect separately from bitmap instruction scheduling.
*/
// References: saffah_v41f1017_agg4 https://duck.ac/submission/139828 supplies the current coalesced8-byte cached-signature row and scratch backing, with its inherited chain including saffah_codex_6a_agg3 https://duck.ac/submission/139516. All attribution and notices remain below.
// Idea/Purpose: in the existing per-point record pass, quantize eight coordinates and eight bound values with two AVX2 multiply/shift groups. Produce both the dense signature row and the query's coordinate-minus-one bytes from the same vector, replacing15 scalar multiplies and individual byte stores. Select low bytes with PSHUFB, preserving unsigned32-bit multiplication wrap and byte truncation even for theINT_MAX anchor sentinel. Store exactly the original8+7 query bytes, preserving all neighbouring padding; fallback dimensions keep their original scalar path. This is different from the removed separate8-point signature pass: both headers are generated where their source row is already live.
// ===== lane 1014 huge-page/TLB layer (agg4), 2026-10-22 =====
// Base: work/1014/agg4_audit3_sigfold_139516.cpp = our submission #139727 (audit3 sigfold
// candidate, sha256 04a96fcf0a435179246e0a24f15f9088dff5e4bf2e01a65db7b9b3d43606c11c); its
// full inherited citation chain is retained verbatim below this block.
// Single-variable change (memory backing only; no arithmetic and no data-layout change):
// (1) pool_alloc now backs every large scratch pool with a 2MB-aligned anonymous mmap and hints
// MADV_HUGEPAGE on it (malloc fallback if mmap fails);
// (2) solve_mf hints MADV_HUGEPAGE on the large static working arrays before first use.
// Reused idea/mechanism citation: the 2MB-aligned MADV_HUGEPAGE scratch-backing mechanism with
// runtime fallback was introduced in this campaign by agg4 main's #139768 candidate
// (work/1005e7/CAND_HPTHP_agg4.cpp, https://duck.ac/submission/139768), itself an adoption of
// pdoom's public submission https://duck.ac/submission/139676 with a modified allocation layer.
// Public duck.ac submissions carry no separate license statement; attribution recorded here.
// Only the mechanism idea was reused; the code below was written for this engine's own arrays.
// Purpose: 1014's hot passes are random-access over ~0.8 GB of tables; 4KB-page dTLB coverage is
// the binding constraint on the local host. Backing those tables with 2MB pages cuts TLB misses.
// Value-exact evidence: judge-style gate k=8 n=1e6 seed=1 -> checksum 3889404381, full-output
// FNV-1a eff6d5432dab9c80 (both identical to #139727's gate values); brute force k=8/9/10 x
// n=200/777/1500 all 0 mismatches; in-process paired A/B diff=0.
// Timing evidence (local, isolated fresh-process alternating runs, n=1e6 seed=1, CPU time):
// base 9625.8/9674.3/11011.6 ms vs variant 8412.0/8486.8/8433.8 ms -> min-vs-min -12.61%.
// Judge caveat: a prior probe recorded MADV_HUGEPAGE returning EPERM in the judge sandbox
// (work/1014/agg4_audit3_EVIDENCE_139516_sigfold.json, negative_findings[0]); if that still
// holds, the hints are no-ops and the mmap path is performance-neutral (glibc malloc already
// uses mmap for these >128KB pools). Peak RSS +19 MB (831 -> 850 MB) at n=1e6.
// 0 submit; coordinator prices.
// ============================================================
// ===== lane 1014 audit3 (agg4) recovery seat, 2026-10-22 =====
// Base: work/1014/adopt_139516_agg4.cpp = our submission #139572, an unmodified adoption of
// saffah_codex_6a_agg3 https://duck.ac/submission/139516 (all inherited notices retained below).
// Single-variable change (work deletion, no functional edit):
// Delete the separate prepare_cached_signatures<8> pass. The 8-bit signatures it materialised
// for point i, dim e are exactly max(sigfull(x[e][i]),1) = pad[e]+1, where pad[e] is already
// computed for the s_pm record in the existing build loop. cand3 therefore packs those bytes
// into one u64 store per point during that loop (dense s_rsig8, 8 B/point) and the s_sig
// transpose for the k=8 path reads the dense 8-byte rows; k=9/10 keep the original prepare
// path byte-for-byte. The s_sig values themselves are unchanged.
// Purpose: delete one O(N) preparation pass (direct in-binary cost 5.23 ms @ n=300k,
// 16.55 ms @ n=1e6 on the audit host) and halve the s_sig random-read footprint (16 MB -> 8 MB).
// Value-exact evidence: judge-style gate k=8 n=1e6 seed=1 -> checksum 3889404381 and full-out[]
// FNV-1a eff6d5432dab9c80, both identical to base (also seed 2/3 @ n=300k and seed 2 @ n=1e6);
// LOCAL_TEST brute-force matrix k=8/9/10 x modes 0-4 x n=200/777/1500 all 0 mismatches.
// 0 submit; coordinator prices. Expected official effect: ~8-18 ms of the 45.71 ms strict gap
// (prepare pass ~7.8 ms official-equivalent; s_sig footprint gain is judge-L3 dependent), i.e.
// estimated official time ~4,513,000-4,523,500 us (4.513-4.523 s), NOT a strict met candidate on its own.
// =================================================================
/* Adopted 2026-10-22 (agg4 main) from public duck.ac submission 139516 by
* saffah_codex_6a_agg3 (https://duck.ac/submission/139516): its interval/vector-fold engine and improvements over #139506/#132595;
* all inherited notices below are retained. License: public submission, no
* separate license stated. Approach/purpose: price the strongest available
* engine as our baseline before further optimization. No functional edits. */
// References: saffah_v41f1017_agg4 https://duck.ac/submission/132595 supplies the current eight-dimensional engine, with inherited notices retained below.
// saffah_codex_6a_agg3 https://duck.ac/submission/139506 and its cited139423/139416/129550 provide sorted bucket intervals, vector-fold stable grouping and cached SIMD signature construction.
// Idea/Purpose: apply these missing preparation mechanisms to the active eight-coordinate path. Replace packed random bucket metadata with monotone intervals; compute unsigned minima in vectors, count and scatter stable groups directly, and quantize each coordinate once before the existing ordered transpose.
// Guard absent coordinates8/9 and their output planes; keep all query/filter algorithms, ties and fallback dimensions intact. This removes redundant setup traffic and scalar quantization without altering numerical representation.
// ===== 引用区 =====
// [1] 原作者账号: saffah_cc_v41_agg1
// 原提交 URL: https://duck.ac/submission/130654
// 所用内容: 直接完整复制其计算引擎正文（为满足代码长度限制删除散文注释；代码逐字节等价，
// 以同一输入在评测机路径驱动下校验输出 checksum 与未删注释版相同），跨题用于本题。
// 许可证遵守: duck.ac 公开提交均未附独立许可证声明；按项目 INSTRUCTIONS 要求在同一评测站内
// 逐项注明原作者账号与原提交 URL 后使用；本发不新增第三方代码。
// [2] 上述提交正文内继承的全链引用（原注释所列提交号，逐一列明）:
// #103361, #103367, #106719, #106747, #106762, #108184, #108329, #108359, #108399, #108401, #108502, #109438, #111338, #121252, #122788, #122903, #123047, #123142, #126143, #126941, #127524, #128321, #128473, #129085, #129575, #130301, #130561, #130562
// ===== 思路区 =====
// 单变量实验（1014）: 1016 引擎 #131443 自动路径 QCOEF=122，优化级 O3→O2（字符串首项），其余逐字节不动。
// 目的: 补齐 0.40%。本地 checksum 同 3889404381。
// ============================================================
#define MF_NOSDMP 1
#ifndef MF_PFXMASK
#define MF_PFXMASK 1
#define MF_PFXSW 24
#endif
#pragma GCC optimize("O2","unroll-loops","rename-registers","no-strict-aliasing","align-functions=32","align-loops=256,no-gcse,live-range-shrinkage,no-if-conversion,no-tree-coalesce-vars,no-tree-partial-pre,no-ivopts")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#undef __code_model_small__
#undef __MMX_WITH_SSE__
#define __tune_skylake__ 1
#define __code_model_small__ 1
#define __SSE3__ 1
#define __SSSE3__ 1
#define __SSE4_1__ 1
#define __SSE4_2__ 1
#define __AVX__ 1
#define __AVX2__ 1
#define __BMI__ 1
#define __BMI2__ 1
#define __LZCNT__ 1
#define __CRC32__ 1
#define __POPCNT__ 1
#define __XSAVE__ 1
#define __MMX_WITH_SSE__ 1
#ifndef MF_FRDPF
#define MF_FRDPF 32
#endif
#ifndef MF_CNTPF
#define MF_CNTPF 64
#endif
#ifndef MF_GQPF
#define MF_GQPF 32
#endif
#define MF_QCOEF 118
#define MF_NOFILTB 1
#include <immintrin.h>
#include <type_traits>
#include <sys/mman.h>
#include <stdint.h>
static inline __m256i LDU(const void *p) {
__m256i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline __m256i LDA(const void *p) {
__m256i v;
__asm__("vmovdqa %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline void STU(void *p, __m256i v) {
__asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline void STA(void *p, __m256i v) {
__asm__ volatile("vmovdqa %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline __m128i LDQ(const void *p) {
__m128i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m128i *)p));
return v;
}
#define LDU_M(p) ({ __m256i _v; __asm__("vmovdqu %1, %0" : "=x"(_v) : "m"(*(const __m256i *)(p))); _v; })
#define LDA_M(p) ({ __m256i _v; __asm__("vmovdqa %1, %0" : "=x"(_v) : "m"(*(const __m256i *)(p))); _v; })
#define STU_M(p,v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i *)(p)) : "x"(_vv)); } while (0)
#define STA_M(p,v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqa %1, %0" : "+m"(*(__m256i *)(p)) : "x"(_vv)); } while (0)
#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 short u16;
typedef unsigned char u8;
typedef unsigned long long u64;
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#ifndef MAXN
#define MAXN 1000005
#define MF_SPL ((size_t)1000064)
#define MF_PFB1 64
#define MF_PFB2 64
#define MF_PFB3 8
#endif
#define MAXK 10
#define STRIDE 16
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 122
#endif
#ifndef MF_KSPAD
#define MF_KSPAD 1
#endif
#ifndef MF_FILTAUTO
#define MF_FILTAUTO 1
#endif
#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
#ifndef MF_FB2
#define MF_FB2 0
#endif
#ifndef MF_FB3
#define MF_FB3 1
#endif
#ifndef MF_PLEN
#define MF_PLEN 1
#endif
#ifndef MF_BUDGETMB
#define MF_BUDGETMB 1940
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
#ifndef MF_PFB
#define MF_PFB 24
#endif
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
#ifndef MF_WTRIM
#define MF_WTRIM 0
#endif
#define AS_WORDS 20480
#ifndef MF_PFA
#define MF_PFA 1
#endif
#ifndef MF_PFAOUT
#define MF_PFAOUT 1
#endif
#ifndef MF_PFAOUTD
#define MF_PFAOUTD 24
#endif
#ifndef MF_SCAPF
#define MF_SCAPF 16
#endif
#define Z16S_PFW(p) __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(p)))
#ifndef MF_ORDPF
#define MF_ORDPF 24
#endif
#ifndef MF_NIDPF
#define MF_NIDPF (12)
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 1408
#endif
static u32 s_ord[MAXK][MAXN];
static u32 s_frord_[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u16 s_bd[(size_t)MAXN * MAXK];
static u16 s_bdT[(size_t)MAXK * MAXN];
static u8 s_g[MAXN];
static u8 s_rsig8[(size_t)MAXN*8] __attribute__((aligned(64)));
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
alignas(64) static u8 compact_hist[MAXN+8];
alignas(64) static u32 compact_cursor[MAXN+8];
static u16 s_bkt[MAXN];
static u16 s_bka[MAXN];
static u32 s_nid[MAXN];
static u32 s_npos[MAXN] __attribute__((aligned(64)));
static u32 s_blkk[64], s_blkp[64];
static u8 *s_tm = 0; static size_t s_tmcap = 0;
static u64 *s_edg = 0; static size_t s_edgcap = 0;
static u64 s_pref[256] __attribute__((aligned(64)));
static int g_dm = -1;
static u32 g_tstride = 0;
static u32 s_tofs[MAXQ + 1], s_rlen[MAXQ + 1];
static u8 s_cntx[16384];
static int g_mask_on = 0;
static int g_mask_ok = 0;
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN];
static u64 s_facc[(size_t)(MAXK - 1) * (((MAXN) + 63) / 64 + 16)];
static u32 s_bnde[(size_t)MAXK * (MAXQ + 2)];
static u64 *s_fbm = 0;
static u64 s_fbuf[8192] __attribute__((aligned(64)));
static __m256i khi8v, khi8hv, kbadhi, kokhi;
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_sdmp[MAXK];
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
static u32 *s_avpool = 0;
static size_t s_avcap = 0;
static u8 *s_sig = 0;
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
static u32 g_sigmul;
static inline u32 sigfull(u32 x){return (x*g_sigmul)>>23;}
template<int C>
static void prepare_cached_signatures(u8 *dst,const unsigned **x,u32 N){
const __m256i mul=_mm256_set1_epi32(g_sigmul),one=_mm256_set1_epi32(1);
const __m256i gather=_mm256_setr_epi8(0,4,8,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,4,8,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1);
const __m128i mask2=_mm_set_epi64x(0,65535);
u32 i=0;
for(;i+8<=N;i+=8){
__m128i c[10];c[8]=c[9]=_mm_setzero_si128();
#pragma GCC unroll 9
for(int d=0;d<C;++d){
__m256i q=_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x[d]+i)),mul);
q=_mm256_max_epu32(_mm256_srli_epi32(q,23),one);
q=_mm256_shuffle_epi8(q,gather);
c[d]=_mm_unpacklo_epi32(_mm256_castsi256_si128(q),_mm256_extracti128_si256(q,1));
}
__m128i t0=_mm_unpacklo_epi8(c[0],c[1]),t1=_mm_unpacklo_epi8(c[2],c[3]);
__m128i t2=_mm_unpacklo_epi8(c[4],c[5]),t3=_mm_unpacklo_epi8(c[6],c[7]);
__m128i u0=_mm_unpacklo_epi16(t0,t1),u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3),u3=_mm_unpackhi_epi16(t2,t3);
__m128i r[4]={_mm_unpacklo_epi32(u0,u2),_mm_unpackhi_epi32(u0,u2),_mm_unpacklo_epi32(u1,u3),_mm_unpackhi_epi32(u1,u3)};
__m128i tail=_mm_unpacklo_epi8(c[8],c[9]);
#pragma GCC unroll 8
for(int p=0;p<8;++p){
__m128i lo=(p&1)?_mm_srli_si128(r[p/2],8):r[p/2];
__m128i hi=_mm_and_si128(_mm_srli_si128(tail,p*2),mask2);
_mm_storeu_si128((__m128i*)(dst+(size_t)(i+p)*16),_mm_unpacklo_epi64(lo,hi));
}
}
for(;i<N;++i)for(int d=0;d<C;++d){u32 v=sigfull(x[d][i]);dst[(size_t)i*16+d]=(u8)(v?v:1u);}
}
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;
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)]
static void hp_hint(void *p, size_t len) {
if (!p || !len) return;
uintptr_t a = (uintptr_t)p & ~(uintptr_t)4095u;
uintptr_t b = ((uintptr_t)p + len + 4095u) & ~(uintptr_t)4095u;
madvise((void *)a, b - a, MADV_HUGEPAGE);
}
static void *pool_alloc(size_t bytes) {
size_t len = bytes + (4u << 20);
void *raw = mmap(0, len, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
if (raw != MAP_FAILED) {
void *al = (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
madvise(al, bytes + (2u << 20), MADV_HUGEPAGE);
return al;
}
raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
return (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
}
static constexpr u32 DEFER_CAP=1u<<19;
alignas(64) static u32 defer_entries[(size_t)DEFER_CAP*9];
static u32 defer_count;
static void flush_deferred(u32*out,const u32*order){
u32 n=defer_count;defer_count=0;
for(u32 b=0;b<n;b+=128){
u32 end=b+128<n?b+128:n;
for(u32 j=b;j<end;j++){
const u32*e=defer_entries+(size_t)j*9;
u32 id=e[1];
__builtin_prefetch(out+e[0],1,3);
_mm_prefetch((const char*)(s_pm+(size_t)id*STRIDE),_MM_HINT_T0);
}
for(u32 j=b;j<end;){
const u32*e=defer_entries+(size_t)j*9;
const u32 query=e[0];
__m256i lo=_mm256_cvtepu16_epi32(_mm_loadu_si128((const __m128i*)(e+2)));
__m256i hi=_mm256_cvtepu8_epi32(_mm_cvtsi64_si128((long long)_pdep_u64(e[6],0x0f0f0f0f0f0f0f0full)));
__m256i th=_mm256_or_si256(lo,_mm256_slli_epi32(hi,16));
__m128i upper=_mm_loadl_epi64((const __m128i*)(e+7));
u32 count=0;
do{
const u32*p=s_pm+(size_t)defer_entries[(size_t)j*9+1]*STRIDE;
__m256i a=_mm256_cmpgt_epi32(th,_mm256_load_si256((const __m256i*)p));
__m128i h=_mm_cmpgt_epi32(upper,_mm_loadl_epi64((const __m128i*)(p+8)));
u32 mask=(u32)_mm256_movemask_epi8(a);
mask&=(u32)_mm_movemask_epi8(h)|0xffffff00u;
count+=(mask==0xffffffffu);
++j;
}while(j<end&&defer_entries[(size_t)j*9]==query);
out[query]+=count;
}
}
}
__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, acc0 = z, acc1 = z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 8 <= Wr; w += 8) {
if (KK > 1) _mm_prefetch((const char *)(bs[1] + w + 32), _MM_HINT_T0);
if (KK > 2) _mm_prefetch((const char *)(bs[2] + w + 32), _MM_HINT_T0);
if (KK > 3) _mm_prefetch((const char *)(bs[3] + w + 32), _MM_HINT_T0);
if (KK > 4) _mm_prefetch((const char *)(bs[4] + w + 32), _MM_HINT_T0);
if (KK > 5) _mm_prefetch((const char *)(bs[5] + w + 32), _MM_HINT_T0);
if (KK > 6) _mm_prefetch((const char *)(bs[6] + w + 32), _MM_HINT_T0);
if (KK > 7) _mm_prefetch((const char *)(bs[7] + w + 32), _MM_HINT_T0);
if (KK > 8) _mm_prefetch((const char *)(bs[8] + w + 32), _MM_HINT_T0);
__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 lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)));
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)));
if (++cnt == 31) {
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
acc0 = acc1 = z; cnt = 0;
}
}
for (; w + 4 <= Wr; w += 4) {
__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);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, 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;
}
template<int NS>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_pm(const u64 *const *bs, const u8 *trow, const u64 *edgp, u32 Wf, u32 cut) {
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();
const __m256i ones=_mm256_set1_epi64x(-1);
const __m128i sub64v=_mm_set1_epi16(0x4000);
__m256i tot = z, acc0 = z, acc1 = z;
u32 w = 0, total = 0;
auto scan_full=[&](auto clamp_tag,u32 stop){
constexpr bool Clamped=decltype(clamp_tag)::value;
for (; w + 8 <= stop; w += 8) {
constexpr int NF=NS<4?NS:4;
const u32 ld32 = Clamped ? ((w + 32u < Wf) ? 32u : 0u) : 32u;
const u32 ld32b = Clamped ? ((w + 64u < Wf) ? 64u : 0u) : 64u;
for (int d = 0; d < NF; d++) { _mm_prefetch((const char *)(bs[d] + w + ld32), _MM_HINT_T0); _mm_prefetch((const char *)(bs[d] + w + ld32b), _MM_HINT_T0); }
_mm_prefetch((const char *)(trow + (w >> 1) + (Clamped ? ((w + 64u < Wf) ? 32u : 0u) : 32u)), _MM_HINT_T0);
__m256i a = LDU256((bs[0] + w));
__m256i b = LDU256((bs[0] + w + 4));
for (int d = 1; d < NF; d++) {
a = _mm256_and_si256(a, LDU256((bs[d] + w)));
b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
}
u32 tb = *(const u32 *)(trow + (w >> 1));
__m128i tv=_mm_cvtsi32_si128((int)tb);
__m128i tdup=_mm_unpacklo_epi8(tv,tv);
__m128i tsub=_mm_subs_epu8(tdup,sub64v);
a = _mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(tsub)),ones));
b = _mm256_and_si256(b,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(_mm_srli_si128(tsub,4))),ones));
if constexpr(NS>4){
__m256i live=_mm256_or_si256(a,b);
if(_mm256_testz_si256(live,live))continue;
for(int d=4;d<NS;d++){
_mm_prefetch((const char*)(bs[d]+w+ld32),_MM_HINT_T0);
a=_mm256_and_si256(a,LDU256(bs[d]+w));
b=_mm256_and_si256(b,LDU256(bs[d]+w+4));
}
}
__m256i lo0 = _mm256_and_si256(a, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lo1 = _mm256_and_si256(b, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0)), z));
acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1)), z));
}
};
scan_full(std::false_type{}, Wf>64u?Wf-64u:0u);
scan_full(std::true_type{},Wf);
for (; w + 4 <= Wf; w += 4) {
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < NS; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
u32 tb = *(const u16 *)(trow + (w >> 1));
__m128i tv=_mm_cvtsi32_si128((int)tb);
__m128i tdup=_mm_unpacklo_epi8(tv,tv);
__m128i tsub=_mm_subs_epu8(tdup,sub64v);
a = _mm256_and_si256(a,_mm256_xor_si256(_mm256_sllv_epi64(ones,_mm256_cvtepu8_epi64(tsub)),ones));
__m256i lo = _mm256_and_si256(a, lm);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(_mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)), z));
}
total = hsum256(_mm256_add_epi64(acc0,acc1));
for (; w < Wf; w++) {
u64 v = bs[0][w];
for (int d = 1; d < NS; d++) v &= bs[d][w];
u32 c = trow[w >> 1];
u32 k = (w & 1u) ? ((c > 64u) ? c - 64u : 0u) : c;
total += (u32)__builtin_popcountll(v & s_pref[k]);
}
{
u32 c = trow[Wf >> 1];
__m128i v = _mm_loadu_si128((const __m128i*)(bs[0] + Wf));
for(int d=1;d<NS;d++) v=_mm_and_si128(v,_mm_loadu_si128((const __m128i*)(bs[d]+Wf)));
v=_mm_and_si128(v,_mm_loadu_si128((const __m128i*)(edgp+(size_t)cut*2u)));
v=_mm_and_si128(v,_mm_set_epi64x((c>64u)?s_pref[c-64u]:0ull,s_pref[c]));
total+=(u32)__builtin_popcountll((u64)_mm_cvtsi128_si64(v));
total+=(u32)__builtin_popcountll((u64)_mm_extract_epi64(v,1));
}
return total;
}
#define T16E_QTCALC() \
const u32 *Avp_ = avpool + (size_t)i * STRIDE; \
__m256i qr_ = _mm256_or_si256(LDA((const __m256i *)rec), khi8v); \
__m256i qT = _mm256_min_epi32(qr_, _mm256_or_si256(LDU256(Avp_), notexlo)); \
__m256i qTh; \
if (k > 8) { \
__m256i qrh_ = _mm256_or_si256(LDU256((rec + 8)), khi8hv); \
__m256i Avmh_ = _mm256_or_si256(LDU256((Avp_ + 8)), notexloh); \
qTh = _mm256_min_epi32(qrh_, Avmh_); \
} else { qTh = all8; }
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 u32 KS_C = (KK <= 8) ? 8u : (u32)KK;
const int k = KK;
const int K1 = KK - 1;
u32 *qo = s_qord, *qt = s_tmp;
#if MF_ASCR
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
const u32 dmx = (u32)g_dm, emx = (dmx < ds) ? dmx : dmx - 1u;
{const u32 d2=(ds<=1u)?2u:1u;u32 *cv=s_cntv;
for(u32 b=0;b<=Q;b++)cv[b]=0;
for(u32 qi=q0;qi<q1;qi++)cv[s_bdT[(size_t)d2*N+qo[qi]]+1]++;
u32 ac=q0;for(u32 b=0;b<=Q;b++){u32 c=cv[b];cv[b]=ac;ac+=c;}
for(u32 qi=q0;qi<q1;qi++){u32 id=qo[qi];qt[cv[s_bdT[(size_t)d2*N+id]+1]++]=id;}
u32*sw=qo;qo=qt;qt=sw;}
{
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_PFB1 < q1) _mm_prefetch((const char *)(s_bdT + (size_t)d1 * N + qo[qi + MF_PFB1]), _MM_HINT_T0);
#endif
{ u32 bb = s_bdT[(size_t)d1 * N + qo[qi]];
s_bkt[qi] = bb; cntv[bb + 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 = qo[qi]; u32 pos = cntv[s_bkt[qi] + 1]++;
qt[pos] = id; s_bka[pos] = s_bkt[qi]; }
{ u32 *swp = qo; qo = qt; qt = swp; }
}
#else
const u32 d1 = 0, e1 = 0;
#endif
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_bdT[(size_t)dd * N + qo[qi]] + 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 = qo[qi]; qt[cntv[s_bdT[(size_t)dd * N + id] + 1]++] = id; }
{ u32 *swp = qo; qo = qt; qt = swp; }
}
for (u32 qi = q0; qi < q1; qi++) s_bka[qi] = s_bdT[(size_t)d1 * N + qo[qi]];
}
#define PF_DIST 4
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = qo[qi + PF_DIST];
const u16 *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 + (g_plen_on ? s_pofs[(size_t)e * MAXQ + bdn[d]]
: ((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_bka[qi];
u32 qe = qi + 1;
while (qe < q1 && s_bka[qe] == 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][qo[t]]; if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 2;
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 = qo[qq];
const u16 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFAOUT
{ u32 oo = qq + MF_PFAOUTD; if (oo < qe) _mm_prefetch((const char *)(out + qo[oo]), _MM_HINT_T0); }
#endif
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = qo[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];
u16 brank[MAXK];
int nbx = 0;
for (int e = 0; e < K1; e++) {
if (g_mask_ok && e == (int)emx) continue;
int rd=(e<ds)?e:e+1; brank[nbx]=bd[rd];
if (e == (int)e1) bp[nbx++] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[nbx++] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
if(g_mask_ok && nbx==8){
__m128i ranks=_mm_loadu_si128((const __m128i*)brank);
const __m128i lanes=_mm_setr_epi16(0,1,2,3,4,5,6,7);
int selected[4];unsigned used=0;
#pragma GCC unroll 4
for(int z=0;z<4;z++){
__m128i v=_mm_minpos_epu16(ranks);
unsigned j=(unsigned)_mm_extract_epi16(v,1);
selected[z]=j;used|=1u<<j;
ranks=_mm_or_si128(ranks,_mm_cmpeq_epi16(lanes,_mm_set1_epi16(j)));
}
const u64* ordered[8];
for(int z=0;z<4;z++)ordered[z]=bp[selected[z]];
unsigned rest=_pext_u32(0xFAC688u,_pdep_u32((~used)&255u,0x249249u)*7u);
for(int z=0;z<4;z++){ordered[4+z]=bp[rest&7u];rest>>=3;}
for(int z=0;z<8;z++)bp[z]=ordered[z];
}else if(g_mask_ok && nbx>=8){
for(int z=0;z<4;z++){
int best=z;
u16 minrank=brank[z];
for(int j=z+1;j<nbx;j++){
const u16 rank=brank[j];
const bool less=rank<minrank;
best=less?j:best;
minrank=less?rank:minrank;
}
const u64*p=bp[z];bp[z]=bp[best];bp[best]=p;
u16 rr=brank[z];brank[z]=brank[best];brank[best]=rr;
}
}
cnt = g_mask_ok
? and_popcount_pm<KK - 2>(bp, s_tm + s_tofs[bd[dmx]], s_edg + (((size_t)(R >> 7) * 129u) << 1), (R >> 7) << 1, R & 127u)
: 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 = qo[qi];
const u16 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = qo[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);
}
}
#if MF_PFAOUT
{ u32 oo = qi + MF_PFAOUTD;
if (oo < q1) _mm_prefetch((const char *)(out + qo[oo]), _MM_HINT_T0); }
#endif
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = qo[qi + PF_DIST];
const u16 *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 + (g_plen_on ? s_pofs[(size_t)e * MAXQ + bdn[d]]
: ((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;
u32 *const avpool = s_avpool;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
const u32 *fo = s_frord_[d];
u32 exmask = 0;
for (int e = 0; e < d; e++) if (e != ds) exmask |= 1u << e;
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];
alignas(16) u8 qmask_bytes[16],dmask_bytes[16];
for(int e=0;e<16;e++){qmask_bytes[e]=(e<d&&e!=(int)ds)?0:255;dmask_bytes[e]=(e==d)?255:0;}
const __m128i qmask=_mm_load_si128((const __m128i*)qmask_bytes),dmask=_mm_load_si128((const __m128i*)dmask_bytes);
for (u32 qi = q0; qi < q1; qi++) {
u32 i = fo[qi] & 0xFFFFFu;
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = fo[qj] & 0xFFFFFu;
_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 *)(out + j), _MM_HINT_T0);
}
}
#endif
const u32 lo = FPOS(d, fo[qi] >> 20);
const u32 hi = s_gt[d][i];
if (lo >= hi) continue;
const u32 *rec = s_pm + (size_t)i * STRIDE;
const u8 *const qpad = ((const u8 *)rec) + 40;
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (size_t)lo * KS_C) : 0;
if (g_use_filt) { T16E_QTCALC();
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 u32 FW = g_fbm_W;
int first = 1;
#if MF_FR_NOREAD
first = 0;
for (u32 w = 0; w < nw; w++) fb[w] = 0;
#endif
for (int e = 0; e < k; e++) {
if (e == d) continue;
u32 c = (u32)((u64)rec[e] / wcs) + 1;
if (c > g_Qc) c = g_Qc;
const u64 *row = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
if (first) { for (u32 w = 0; w < nw; w++) fb[w] = row[w]; first = 0; }
else { for (u32 w = 0; w < nw; w++) fb[w] &= row[w]; }
}
fb[0] &= ~0ull << (lo & 63);
u32 r2 = hi & 63u;
if (r2) fb[nw - 1] &= (1ull << r2) - 1;
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;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * KS_C) : (s_pm + (size_t)s_ord[d][j] * STRIDE);
__m256i pv = LDU256(cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (g_sig_on && (hi - lo) >= 2) {
__m256i qsb[16];
{
__m128i qr8=_mm_loadu_si128((const __m128i*)qpad);
__m128i qb8=_mm_loadl_epi64((const __m128i*)(qpad+k));
if constexpr(KK>=10)qb8=_mm_insert_epi8(qb8,qpad[k+8],8);
__m128i qs=_mm_or_si128(_mm_min_epu8(qr8,_mm_or_si128(qb8,qmask)),dmask);
__m256i qwide=_mm256_broadcastsi128_si256(qs);
for(int e=0;e<k-1;e++)qsb[e]=_mm256_shuffle_epi8(qwide,_mm256_set1_epi8((char)(e<d?e:e+1)));
}
const u8 *const sgbase = s_sig + (size_t)d * (k-1) * MF_SPL;
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
__m256i ac0=_mm256_setzero_si256(),ac1=ac0;
for(int e=0;e<k-1;e++){
const __m256i df=_mm256_subs_epu8(LDU256(sgbase + (size_t)e * MF_SPL + aj),qsb[e]);
if(e&1)ac1=_mm256_or_si256(ac1,df);else ac0=_mm256_or_si256(ac0,df);
}
__m256i ac=_mm256_or_si256(ac0,ac1);
__m256i lt = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = (u32)_mm256_movemask_epi8(lt);
u32 amb = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)));
aacc += (u32)__builtin_popcount(ltm);
if (amb) { T16E_QTCALC();
__m256i lv=_mm256_and_si256(qT,_mm256_set1_epi32(65535));
__m128i packed_low=_mm_packus_epi32(_mm256_castsi256_si128(lv),_mm256_extracti128_si256(lv,1));
__m256i hv=_mm256_and_si256(_mm256_srli_epi32(qT,16),_mm256_set1_epi32(15));
__m128i hw=_mm_packus_epi32(_mm256_castsi256_si128(hv),_mm256_extracti128_si256(hv,1));
__m128i hb=_mm_packus_epi16(hw,hw);
u32 packed_high=(u32)_pext_u64((u64)_mm_cvtsi128_si64(hb),0x0f0f0f0f0f0f0f0full);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
if(KK==10 && !c0 && defer_count<DEFER_CAP){
u32*e=defer_entries+(size_t)defer_count++*9;
e[0]=i;e[1]=odd[aj+(u32)t];
_mm_storeu_si128((__m128i*)(e+2),packed_low);e[6]=packed_high;
_mm_storel_epi64((__m128i*)(e+7),_mm256_castsi256_si128(qTh));
continue;
}
const u32 *cp = c0 ? (c0 + (size_t)(aj + (u32)t - lo) * KS_C) : (s_pm + (size_t)odd[aj + (u32)t] * STRIDE);
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
aacc += okl;
}
}
}
{
u32 r = hi - aj;
if (r) {
u32 base, lk;
if (hi >= 32u) { base = hi - 32u; lk = 32u - r; }
else { base = 0u; lk = lo; }
u32 kmask = ((r >= 32u) ? 0xFFFFFFFFu : ((1u << r) - 1u)) << lk;
const u8 *const tp = sgbase + base;
__m256i ac0=_mm256_setzero_si256(),ac1=ac0;
for(int e=0;e<k-1;e++){
const __m256i df=_mm256_subs_epu8(LDU256(tp + (size_t)e * MF_SPL),qsb[e]);
if(e&1)ac1=_mm256_or_si256(ac1,df);else ac0=_mm256_or_si256(ac0,df);
}
__m256i ac=_mm256_or_si256(ac0,ac1);
__m256i lt = _mm256_cmpeq_epi8(ac, _mm256_setzero_si256());
u32 ltm = ((u32)_mm256_movemask_epi8(lt)) & kmask;
u32 amb = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(ltm);
if (amb) { T16E_QTCALC();
__m256i lv=_mm256_and_si256(qT,_mm256_set1_epi32(65535));
__m128i packed_low=_mm_packus_epi32(_mm256_castsi256_si128(lv),_mm256_extracti128_si256(lv,1));
__m256i hv=_mm256_and_si256(_mm256_srli_epi32(qT,16),_mm256_set1_epi32(15));
__m128i hw=_mm_packus_epi32(_mm256_castsi256_si128(hv),_mm256_extracti128_si256(hv,1));
__m128i hb=_mm_packus_epi16(hw,hw);
u32 packed_high=(u32)_pext_u64((u64)_mm_cvtsi128_si64(hb),0x0f0f0f0f0f0f0f0full);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
if(KK==10 && !c0 && defer_count<DEFER_CAP){
u32*e=defer_entries+(size_t)defer_count++*9;
e[0]=i;e[1]=odd[base+(u32)t];
_mm_storeu_si128((__m128i*)(e+2),packed_low);e[6]=packed_high;
_mm_storel_epi64((__m128i*)(e+7),_mm256_castsi256_si128(qTh));
continue;
}
const u32 *cp = c0 ? (c0 + (size_t)(base + (u32)t - lo) * KS_C) : (s_pm + (size_t)odd[base + (u32)t] * STRIDE);
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
aacc += okl;
}
}
}
}
add = aacc;
} else if (c0) {
T16E_QTCALC();
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 + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
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 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__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;
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 + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
__m256i h0 = LDU256((c + 8));
__m256i h1 = LDU256((c + KS_C + 8));
__m256i h2 = LDU256((c + 2 * KS_C + 8));
__m256i h3 = LDU256((c + 3 * KS_C + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8) & (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h0)), kokhi);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8) & (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h1)), kokhi);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8) & (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h2)), kokhi);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8) & (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, h3)), kokhi);
c += 4 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__m256i pv = LDU256(c);
__m256i pvh = LDU256((c + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8) &
(u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
add = a0 + a1 + a2 + a3;
}
} else {
T16E_QTCALC();
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__m256i pv = LDA((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDA((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_andnot_si256(kbadhi, _mm256_cmpgt_epi32(qTh, pvh)), kokhi);
}
add += ok;
}
}
out[i] += add;
}
if constexpr(KK==10)flush_deferred(out,odd);
}
}
template<int K>
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out) {
const int k = K;
{ hp_hint(s_ord, sizeof(s_ord)); hp_hint(s_frord_, sizeof(s_frord_)); hp_hint(s_gt, sizeof(s_gt));
hp_hint(s_pm, sizeof(s_pm)); hp_hint(s_bd, sizeof(s_bd)); hp_hint(s_bdT, sizeof(s_bdT));
hp_hint(s_pos, sizeof(s_pos)); hp_hint(s_rsig8, sizeof(s_rsig8));
hp_hint(s_qord, sizeof(s_qord)); hp_hint(s_tmp, sizeof(s_tmp)); hp_hint(s_cntv, sizeof(s_cntv)); }
const u32 W = ((N + 63) / 64 + 7u) & ~7u;
g_W = W;
#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) {
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++;
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
overhead += (u64)k * k * ((u64)N + 64);
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;
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46};
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};
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;
if(q>=MAXQ)q=MAXQ-1;
g_Q = q;
g_B = (N + q - 1) / q;
}
for (u32 t = 0; t < 256u; t++) s_pref[t] = (t >= 64u) ? ~0ull : ((1ull << t) - 1ull);
g_mask_on = MF_PFXMASK;
const u32 Q = g_Q, B = g_B;
const u64 bmag = ((1ull << 40) + B - 1) / B;
#if MF_PLEN
g_plen_on = 1;
#endif
#if MF_WTRIM
u32 Wpool = W;
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1;
#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], bad8[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;
for (int e = 0; e < 8; e++) bad8[e] = (8 + e < k) ? 0 : -1;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
kbadhi = LDU256(bad8);
kokhi = _mm256_xor_si256(kbadhi, _mm256_set1_epi32(-1));
}
const int K1 = k - 1;
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
u8*hc=compact_hist;u32*ct=compact_cursor;
for(u32 i=0;i<=N;i++)hc[i]=0;
u32 lim=N>MF_ORDPF?N-MF_ORDPF:0;
for(u32 i=0;i<lim;i++){
u32 v=xd[i+MF_ORDPF];if(v>=N)v=0;Z16S_PFW(hc+v);
v=xd[i];hc[v<N?v:0]++;
}
for(u32 i=lim;i<N;i++){u32 v=xd[i];hc[v<N?v:0]++;}
u32 sum=0;
u32 hv=0;
for(;hv+8<=N+1;hv+=8){
__m256i c=_mm256_cvtepu8_epi32(_mm_loadl_epi64((const __m128i*)(hc+hv)));
__m256i p=_mm256_add_epi32(c,_mm256_slli_si256(c,4));
p=_mm256_add_epi32(p,_mm256_slli_si256(p,8));
__m128i lo=_mm256_castsi256_si128(p),hi=_mm256_extracti128_si256(p,1);
hi=_mm_add_epi32(hi,_mm_shuffle_epi32(lo,255));
__m256i inc=_mm256_inserti128_si256(_mm256_castsi128_si256(lo),hi,1);
__m256i ex=_mm256_add_epi32(_mm256_sub_epi32(inc,c),_mm256_set1_epi32(sum));
_mm256_storeu_si256((__m256i*)(ct+hv),ex);
sum+=(u32)_mm_extract_epi32(hi,3);
}
for(;hv<=N;hv++){u32 c=hc[hv];ct[hv]=sum;sum+=c;}
if(sum==N && N<1048576u){
u32*gt=s_gt[d];u32*pp=s_pos[d];
for(u32 i=0;i<N;i++){
if(i+2*MF_SCAPF<N){u32 v=xd[i+2*MF_SCAPF];if(v>=N)v=0;_mm_prefetch((const char*)(ct+v),_MM_HINT_T0);}
if(i+MF_SCAPF<N){u32 v=xd[i+MF_SCAPF];if(v>=N)v=0;u32 w=ct[v];Z16S_PFW(od+(w&0xfffffu)+(w>>20));}
u32 v=xd[i];if(v>=N)v=0;
u32 w=ct[v],st=w&0xfffffu,pos=st+(w>>20);ct[v]=w+(1u<<20);
od[pos]=i;pp[i]=pos;gt[i]=st;
#if MF_FB2
s_posp[(size_t)i*MAXK+d]=pos;
#endif
u32 num=(ct[v+1]&0xfffffu)-1u;if(num>N)num=N;
u32 bb=(u32)(((u64)num*bmag)>>40);
s_bdT[(size_t)d*N+i]=(u16)(bb>Q-1?Q-1:bb);
}
}else{
for (u32 i = 0; i <= N; i++) s_cntv[i] = 0;
#if MF_ORDPF > 0
{ u32 i = 0;
for (; i + MF_ORDPF < N; i++) {
u32 vn = xd[i + MF_ORDPF];
Z16S_PFW(s_cntv + (vn < N ? vn : 0));
u32 v = xd[i]; s_cntv[v < N ? v : 0]++;
}
for (; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
}
#else
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
#endif
{ u32 s = 0; for (u32 i = 0; i <= N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tmp[i] = s; s += c; } }
#if MF_FB2
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i;
s_posp[(size_t)i * MAXK + d] = q;
s_pos[d][i] = q;
}
#else
{
u32 *pp = s_pos[d];
for (u32 i = 0; i < N; i++) {
if (i + 2*MF_SCAPF < N) { u32 v2 = xd[i + 2*MF_SCAPF]; if (v2 >= N) v2 = 0;
_mm_prefetch((const char *)(s_tmp + v2), _MM_HINT_T0); }
if (i + MF_SCAPF < N) { u32 vf = xd[i + MF_SCAPF]; if (vf >= N) vf = 0;
Z16S_PFW(&od[s_tmp[vf]]); }
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i; pp[i] = q;
}
}
#endif
{ u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
{if (i + MF_CNTPF < N) { u32 vp = xd[i + MF_CNTPF]; if (vp >= N) vp = 0;
_mm_prefetch((const char *)(s_cntv + vp), _MM_HINT_T0); }}
u32 v = xd[i]; if (v >= N) v = 0;
gt[i] = s_cntv[v];
u32 num = s_cntv[v + 1] - 1u; if (num > N) num = N;
u32 bb = (u32)(((u64)num * bmag) >> 40);
{ u16 bbx = (bb > Q - 1 ? Q - 1 : (u16)bb);
s_bdT[(size_t)d * N + i] = bbx; }
} }
}
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;
}
for (u32 i = 0; i < N; i++) {
u16 *row = s_bd + (size_t)i * MAXK;
for (int d = 0; d < k; d++) row[d] = s_bdT[(size_t)d * N + i];
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
g_sigmul=N>1?(u32)((255ULL<<23)/(N-1)):0;
size_t need = (size_t)k * (k-1) * MF_SPL;
static u8 *sigpool = 0; static size_t sigcap = 0;
if (need > ((size_t)220 << 20) || (size_t)N + 32 > MF_SPL) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
g_sig_on = 1;
g_sig_merged = 0;
}
}
#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 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#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
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), LDA((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;
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 = pd[pe[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), LDA((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();
}
{
size_t needav = (size_t)N * STRIDE * 4 + 256;
if (needav > s_avcap) {
s_avpool = (u32 *)pool_alloc(needav);
s_avcap = s_avpool ? needav : 0;
}
}
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
u32 gstart[MAXK+2]={};
{
u32 i=0;
for(;i+8<=N;i+=8){
__m256i v=_mm256_slli_epi32(_mm256_loadu_si256((const __m256i*)(s_gt[0]+i)),4);
for(unsigned d=1;d<(unsigned)k;++d){
__m256i q=_mm256_or_si256(_mm256_slli_epi32(_mm256_loadu_si256((const __m256i*)(s_gt[d]+i)),4),_mm256_set1_epi32(d));
v=_mm256_min_epu32(v,q);
}
_mm256_storeu_si256((__m256i*)(s_tmp+i),v);
}
for(;i<N;++i){u32 v=s_gt[0][i]<<4;for(unsigned d=1;d<(unsigned)k;++d){u32 q=(s_gt[d][i]<<4)|d;if(q<v)v=q;}s_tmp[i]=v;}
}
for (u32 i = 0; i < N; i++) {
const u32 best=s_tmp[i]>>4;const int ds=(int)(s_tmp[i]&15u);
++gstart[ds+1];
s_g[i] = (u32)ds;
const u16 *bdv = s_bd + (size_t)i * MAXK;
u32 *av = s_avpool + (size_t)i * STRIDE;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
for (int d = 0; d < k; d++) {
if (d == ds) { av[d] = 0x7FFFFFFFu; continue; }
u32 b = bdv[d];
av[d] = BVAL(d, b);
if (g_plen_on && best > mr[(size_t)d * MAXQ + b])
mr[(size_t)d * MAXQ + b] = best;
}
}
if (s_sig) {
const u32 sh = g_sigsh;
for (u32 i = 0; i < N; i++) {
u32 *row = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) row[d] = x[d][i];
const u32 *av = s_avpool + (size_t)i * STRIDE;
u8 *pad = (u8 *)row + 40;
if constexpr(K==8){
const __m256i one=_mm256_set1_epi32(1),mul=_mm256_set1_epi32(g_sigmul);
const __m256i perm=_mm256_setr_epi8(0,4,8,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,4,8,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1);
auto bytes8=[&](__m256i q){q=_mm256_shuffle_epi8(q,perm);return _mm_unpacklo_epi32(_mm256_castsi256_si128(q),_mm256_extracti128_si256(q,1));};
__m256i v=_mm256_max_epu32(one,_mm256_srli_epi32(_mm256_mullo_epi32(_mm256_load_si256((const __m256i*)row),mul),23));
_mm_storel_epi64((__m128i*)(s_rsig8+(size_t)i*8),bytes8(v));
_mm_storel_epi64((__m128i*)pad,bytes8(_mm256_sub_epi32(v,one)));
__m256i a=_mm256_max_epu32(one,_mm256_srli_epi32(_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)av),mul),23));
u64 ap=(u64)_mm_cvtsi128_si64(bytes8(_mm256_sub_epi32(a,one)));__builtin_memcpy(pad+8,&ap,7);
}else{
for(int e=0;e<k;++e){u32 v=sigfull(row[e]);pad[e]=(u8)(v?v-1u:0u);}
for(int e=0;e+1<k;++e){u32 v=sigfull(av[e]);pad[k+e]=(u8)(v?v-1u:0u);}
}
}
} else {
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=1;d<=k+1;++d)gstart[d]+=gstart[d-1];
{u32 cursor[MAXK];for(int d=0;d<k;++d)cursor[d]=gstart[d];
for(u32 i=0;i<N;++i)s_qord[cursor[s_g[i]]++]=i;}
for (int d = 0; d < k; d++) {
u32 pos[MAXK];
for (int ds = 0; ds < k; ds++) pos[ds] = gstart[ds];
const u32 last_bucket=(N-1u)/B;u32 r=0;
for(u32 b=0;b<=last_bucket;++b){
const u32 end=b==last_bucket?N:FPOS(d,b+1);
for(;r<end;++r){
if(r+MF_FRDPF<N){u32 ip=s_ord[d][r+MF_FRDPF];_mm_prefetch((const char*)(s_g+ip),_MM_HINT_T0);}
u32 id=s_ord[d][r];u32 dst=pos[s_g[id]]++;
s_frord_[d][dst]=id|(b<<20);
}
}
}
#if MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 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;
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;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * (k-1) * MF_SPL;
for (int e = 0; e < k; e++) if(e!=d) sb[(size_t)(e-(e>d)) * MF_SPL + j] = (u8)((sigfull(row[e])) | ((sigfull(row[e])) == 0));
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
if (s_sig && !g_sig_merged && k>=8) {
if constexpr (K == 8) {
/* cand2: signatures come from the dense 8-byte rows filled during the s_pm build */
for(int d=0;d<8;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*7*MF_SPL;
u32 j=0;
for(;j+8<=N;j+=8) {
if(j+64<N){const u32*op=od+j+64;
_mm_prefetch((const char *)(s_rsig8+(size_t)op[0]*8),_MM_HINT_T0);
_mm_prefetch((const char *)(s_rsig8+(size_t)op[1]*8),_MM_HINT_T0);
_mm_prefetch((const char *)(s_rsig8+(size_t)op[2]*8),_MM_HINT_T0);
_mm_prefetch((const char *)(s_rsig8+(size_t)op[3]*8),_MM_HINT_T0);
_mm_prefetch((const char *)(s_rsig8+(size_t)op[4]*8),_MM_HINT_T0);
_mm_prefetch((const char *)(s_rsig8+(size_t)op[5]*8),_MM_HINT_T0);
_mm_prefetch((const char *)(s_rsig8+(size_t)op[6]*8),_MM_HINT_T0);
_mm_prefetch((const char *)(s_rsig8+(size_t)op[7]*8),_MM_HINT_T0);}
__m128i r0=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+0]*8));
__m128i r1=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+1]*8));
__m128i r2=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+2]*8));
__m128i r3=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+3]*8));
__m128i r4=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+4]*8));
__m128i r5=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+5]*8));
__m128i r6=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+6]*8));
__m128i r7=_mm_loadl_epi64((const __m128i *)(s_rsig8+(size_t)od[j+7]*8));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
if(d!=0)_mm_storel_epi64((__m128i *)(base+(size_t)(0-(0>d))*MF_SPL+j),v0);
if(d!=1)_mm_storel_epi64((__m128i *)(base+(size_t)(1-(1>d))*MF_SPL+j),_mm_srli_si128(v0,8));
if(d!=2)_mm_storel_epi64((__m128i *)(base+(size_t)(2-(2>d))*MF_SPL+j),v1);
if(d!=3)_mm_storel_epi64((__m128i *)(base+(size_t)(3-(3>d))*MF_SPL+j),_mm_srli_si128(v1,8));
if(d!=4)_mm_storel_epi64((__m128i *)(base+(size_t)(4-(4>d))*MF_SPL+j),v2);
if(d!=5)_mm_storel_epi64((__m128i *)(base+(size_t)(5-(5>d))*MF_SPL+j),_mm_srli_si128(v2,8));
if(d!=6)_mm_storel_epi64((__m128i *)(base+(size_t)(6-(6>d))*MF_SPL+j),v3);
if(d!=7)_mm_storel_epi64((__m128i *)(base+(size_t)(7-(7>d))*MF_SPL+j),_mm_srli_si128(v3,8));
}
for(;j<N;j++) {
const u8 *row=s_rsig8+(size_t)od[j]*8;
for(int e=0;e<8;e++)if(e!=(int)d)base[(size_t)(e-(e>d))*MF_SPL+j]=row[e];
}
}
} else {
static u8 cached_signature[(size_t)MAXN*16] __attribute__((aligned(64)));
if constexpr (K == 9) prepare_cached_signatures<9>(cached_signature,x,N); else prepare_cached_signatures<10>(cached_signature,x,N);
for(int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*(k-1)*MF_SPL;
u32 j=0;
for(;j+8<=N;j+=8) {
if(j+64<N){const u32*op=od+j+64;
_mm_prefetch((const char *)(cached_signature+(size_t)op[0]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[1]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[2]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[3]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[4]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[5]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[6]*16),_MM_HINT_T0);
_mm_prefetch((const char *)(cached_signature+(size_t)op[7]*16),_MM_HINT_T0);}
__m128i r0=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(cached_signature+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
if(d!=0)_mm_storel_epi64((__m128i *)(base+(size_t)(0-(0>d))*MF_SPL+j),v0);
if(d!=1)_mm_storel_epi64((__m128i *)(base+(size_t)(1-(1>d))*MF_SPL+j),_mm_srli_si128(v0,8));
if(d!=2)_mm_storel_epi64((__m128i *)(base+(size_t)(2-(2>d))*MF_SPL+j),v1);
if(d!=3)_mm_storel_epi64((__m128i *)(base+(size_t)(3-(3>d))*MF_SPL+j),_mm_srli_si128(v1,8));
if(d!=4)_mm_storel_epi64((__m128i *)(base+(size_t)(4-(4>d))*MF_SPL+j),v2);
if(d!=5)_mm_storel_epi64((__m128i *)(base+(size_t)(5-(5>d))*MF_SPL+j),_mm_srli_si128(v2,8));
if(d!=6)_mm_storel_epi64((__m128i *)(base+(size_t)(6-(6>d))*MF_SPL+j),v3);
if(d!=7)_mm_storel_epi64((__m128i *)(base+(size_t)(7-(7>d))*MF_SPL+j),_mm_srli_si128(v3,8));
if(k>8 && d!=8)_mm_storel_epi64((__m128i *)(base+(size_t)(8-(8>d))*MF_SPL+j),v4);
if(k>9 && d!=9)_mm_storel_epi64((__m128i *)(base+(size_t)(9-(9>d))*MF_SPL+j),_mm_srli_si128(v4,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*16;
for(int e=0;e<k;e++)if(e!=(int)d)base[(size_t)(e-(e>d))*MF_SPL+j]=row[e];
}
}
}
g_sig_merged=1;
}
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * (k-1) * MF_SPL;
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;
for (int e = 0; e < k; e++) if(e!=d) base[(size_t)(e-(e>d)) * MF_SPL + j] = (u8)((sigfull(row[e])) | ((sigfull(row[e])) == 0));
}
}
}
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;
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
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++) {
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; }
}
#endif
g_mask_ok = 0;
if (g_mask_on) {
{ u32 dd0 = (ds == 0) ? 1u : 0u; g_dm = 0; while ((u32)g_dm == (u32)ds || (u32)g_dm == dd0) g_dm++; }
const u32 dmx = (u32)g_dm;
const u32 *const sdsp = s_pos[ds], *const sdmo = s_pos[dmx], *const sdom = s_ord[dmx], *const sdsr = s_ord[ds];
const u32 *const mrx = s_maxr + ((size_t)ds * MAXK + dmx) * MAXQ;
u32 mblk = 1; { u32 acc = 0;
for (u32 b = 0; b <= Q; b++) {
u32 rr = (mrx[b] >> 7) + 1u; if (rr > Wpool) rr = Wpool; if (rr > mblk) mblk = rr;
u32 rl = (rr + 63u) & ~63u; s_tofs[b] = acc; s_rlen[b] = rl; acc += rl; }
if (mblk > Wpool) mblk = Wpool;
{ size_t tnd = acc;
if (tnd > s_tmcap) { s_tm = (u8 *)pool_alloc(tnd + 4096); s_tmcap = s_tm ? tnd : 0; } } }
{ size_t end = ((size_t)mblk + 1) * 129 * 2;
if (end > s_edgcap) { s_edg = (u64 *)pool_alloc(end * 8 + 4096); s_edgcap = s_edg ? end : 0; } }
if (s_tm && s_edg) {
g_mask_ok = 1; g_tstride = mblk;
{ u8 *const cnt = s_cntx;
for (u32 m = 0; m < mblk; m++) cnt[m] = 0;
u32 j2 = 0;
for (u32 b = 0; b <= Q; b++) {
u32 L = (b == 0u) ? 0u : ((b < Q) ? FPOS(dmx, b) : (u32)N);
while (j2 < L) {
{ u32 pf = j2 + MF_PFXSW; if (pf >= N) pf = j2; u32 pp = sdom[pf];
_mm_prefetch((const char *)(sdsp + pp), _MM_HINT_T0); Z16S_PFW(&s_npos[pp]); }
u32 p = sdom[j2]; j2++;
u32 rk = sdsp[p]; u32 blk = rk >> 7;
if (blk < mblk) { u32 c = cnt[blk]; cnt[blk] = c + 1u; s_npos[p] = (blk << 7) | c; }
else s_npos[p] = rk;
}
{ u8 *row = s_tm + s_tofs[b]; u32 rl = s_rlen[b];
for(u32 off=0;off<rl;off+=64){
_mm256_stream_si256((__m256i*)(row+off),_mm256_loadu_si256((const __m256i*)(cnt+off)));
_mm256_stream_si256((__m256i*)(row+off+32),_mm256_loadu_si256((const __m256i*)(cnt+off+32)));
} }
}
_mm_sfence();
}
{ u32 nb = (N + 127u) >> 7; if (mblk > nb) mblk = nb;
for (u32 m = 0; m < mblk; m++) {
u32 b0 = m << 7, nj = N - b0; if (nj > 128u) nj = 128u;
__m128i edge=_mm_setzero_si128();
_mm_storeu_si128((__m128i*)(s_edg+(size_t)m*258),edge);
const __m128i onebit=_mm_set1_epi64x(1),offset=_mm_set_epi64x(64,0);
for (u32 j = 0; j < nj; j++) {
u32 look=b0+j+32;
if(look<N)_mm_prefetch((const char*)(s_npos+sdsr[look]),_MM_HINT_T0);
u32 sg = s_npos[sdsr[b0 + j]] & 127u;
{ const __m128i shifts=_mm_sub_epi64(_mm_set1_epi64x(sg),offset);
edge=_mm_or_si128(edge,_mm_sllv_epi64(onebit,shifts)); }
_mm_storeu_si128((__m128i*)(s_edg+(size_t)m*258+(j+1u)*2u),edge);
}
for (u32 j = nj + 1u; j <= 128u; j++) { _mm_storeu_si128((__m128i*)(s_edg+(size_t)m*258+j*2u),edge); } } }
} else g_mask_on = 0;
}
const u32 *s_nid = g_mask_ok ? s_npos : s_pos[ds];
const u32 emsk = (g_mask_ok ? ((u32)g_dm < (u32)ds ? (u32)g_dm : (u32)g_dm - 1u) : 0xffffffffu);
for (int e = 0; e < K1; e++) {
if ((u32)e == emsk) continue;
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) {
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
for (u32 b = 1; b <= Q; b++) {
u32 pend = (b < Q) ? FPOS(d, b) : (u32)N;
#if MF_NIDPF > 0
{ u32 pp = ptr;
const u32 nidpf = 40;
for (; pp + 8 <= pend; pp += 8) {
if (pp + nidpf + 8 <= pend) {
_mm_prefetch((const char *)(s_nid + od[pp + nidpf]), _MM_HINT_T0);
_mm_prefetch((const char *)(s_nid + od[pp + nidpf + 4]), _MM_HINT_T0);
}
u32 q0 = s_nid[od[pp + 0]], q1 = s_nid[od[pp + 1]], q2 = s_nid[od[pp + 2]], q3 = s_nid[od[pp + 3]];
u32 q4 = s_nid[od[pp + 4]], q5 = s_nid[od[pp + 5]], q6 = s_nid[od[pp + 6]], q7 = s_nid[od[pp + 7]];
cur[q0 >> 6] |= 1ull << (q0 & 63); cur[q1 >> 6] |= 1ull << (q1 & 63);
cur[q2 >> 6] |= 1ull << (q2 & 63); cur[q3 >> 6] |= 1ull << (q3 & 63);
cur[q4 >> 6] |= 1ull << (q4 & 63); cur[q5 >> 6] |= 1ull << (q5 & 63);
cur[q6 >> 6] |= 1ull << (q6 & 63); cur[q7 >> 6] |= 1ull << (q7 & 63);
}
for (; pp + nidpf < pend; pp++) {
_mm_prefetch((const char *)(s_nid + od[pp + nidpf]), _MM_HINT_T0);
u32 q = s_nid[od[pp]]; cur[q >> 6] |= 1ull << (q & 63);
}
for (; pp < pend; pp++) { u32 q = s_nid[od[pp]]; cur[q >> 6] |= 1ull << (q & 63); }
}
ptr = pend;
#else
while (ptr < pend) { u32 q = s_nid[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
#endif
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
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), LDA((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: defer_count=0;
query_group<10>(N, out, ds, q0, q1, Q, Wpool);
break;
}
}
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf<4>((u32)n, (const unsigned **)x, out); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf<5>((u32)n, (const unsigned **)x, out); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf<6>((u32)n, (const unsigned **)x, out); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf<7>((u32)n, (const unsigned **)x, out); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf<8>((u32)n, (const unsigned **)x, out); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf<9>((u32)n, (const unsigned **)x, out); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf<10>((u32)n, (const unsigned **)x, out); }
#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
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 4.503 s | 772 MB + 400 KB | Accepted | Score: 100 | 显示更多 |