代码

/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140864 supplies this engine; https://duck.ac/submission/140845 introduced the edge gate and https://duck.ac/submission/140865 supplies streamed byte-prefix snapshots. Inherited credits/licenses below retained; no separate license.
Idea/Purpose: copy64-byte-aligned prefix rows with two256-bit non-temporal stores per cacheline, then fence once before consumers. Exact bytes and geometry stay unchanged; test reduced write allocation. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140857 supplies this nine-dimensional final-edge gate and inherited notices below. Our https://duck.ac/submission/140856 supplies the positively measured independent inverse-edge halves, adapted here without its other ten-dimensional behavior. No separate license declared.
Idea/Purpose: after the combined empty-block rejection, reconstruct only each nonempty64-bit half of the final inverse edge. Each omitted half contributes exactly zero; avoid its64-byte row read and two movemasks. Preserve all ordinary bitmap loads, prefix masks and main scans. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140707 supplies this nine-dimensional engine and inherited source/license notices below. Our https://duck.ac/submission/140845 supplies the successfully measured ten-dimensional final-edge survivor gate, adapted here to this engine's seven-plane intersection. Public submissions declared no separate license.
Idea: intersect the final partial128-point block's ordinary planes with its prefix masks before reading its inverse-edge permutation. If both survivor words are zero, the exact final contribution is zero and no edge loads/comparisons/popcounts are needed. Otherwise retain the original edge arithmetic. Main scans, exact fringe checks, query order and mapping remain unchanged.
Purpose: test whether avoided random edge access benefits the still-unmet nine-dimensional family as observed in ten dimensions. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140641 supplies the current nine-dimensional engine and all inherited source/license notices below. Our https://duck.ac/submission/140413 establishes that the ninth threshold is the original coordinate because no later fringe can consume its exclusion bound. No separate license was declared.
Idea: in nine-dimensional signature ambiguity checks, only lane0 of the high vector is live. Compare that exact query threshold with candidate coordinate8 as unsigned scalars, replacing a32-byte candidate load, vector comparison, bad-lane mask and vector test. Keep every fringe, all ambiguity bits, low-eight checks and strict ninth-coordinate condition; ten-dimensional and nonsignature fallback paths remain unchanged.
Purpose: reduce exact ambiguity work without omitting a coordinate comparison or changing signature filtering. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140596 supplies this accepted nine-dimensional engine and all inherited notices; https://duck.ac/submission/140638 supplies our scalar final-edge-mask combination from the ten-dimensional sibling. No separate license was declared.
 * Idea/Purpose: retain scalar inverse-edge and prefix masks for direct combination with the two extracted64-bit bitmap intersections at the final popcounts. This avoids scalar-to-XMM mask construction and two vector ANDs while preserving every selected bit and all traversal geometry. Test the same mechanism at nine dimensions.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140516 supplies the accepted nine-coordinate inverse-edge engine and all inherited notices. saffah_v41f1017_agg4 https://duck.ac/submission/137126 hoists fringe query-shuffle indices once per dimension in its ten-coordinate engine; adapt that explicit invariant placement here. Its public source and this base declare no separate license; inherited citations are retained.
 * Idea: the eight byte-shuffle selection indices depend only on the current fringe dimension d, not on the query. Construct their vectors once before the query loop and reuse those same masks for every query signature broadcast. Threshold values, plane order, query records and all comparisons remain identical.
 * Purpose: official experiment removing repeated index/broadcast construction from the nine-coordinate query loop without changing data geometry or count semantics.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140413 supplies the accepted nine-coordinate compact-bound/online-anchor engine and full inherited notices. Our https://duck.ac/submission/140499 introduced inverse edge permutations and https://duck.ac/submission/140510 adapted them to128-point blocks; this copies that128-bit helper and construction while preserving this problem's online-anchor schedule. No separate licenses were declared.
 * Idea: store the exact original anchor offset at each permuted bitmap position, using128 bytes per block instead of129 pairs of u64 masks. Four byte comparisons and movemasks reconstruct the requested partial-block mask. Original offsets and cut are0..127, absent partial-block positions are127, and all full-word counting, bounds, query order and prefix updates stay unchanged.
 * Purpose: test whether the smaller edge construction and resident footprint improve this distinct nine-coordinate engine, with full output and boundary validation before online measurement.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/139851 supplies this accepted online-anchor nine-dimensional engine and its inherited source/license notices below. saffah_codex_6a_agg3 https://duck.ac/submission/140133 establishes the unused last-coordinate bound in the nine-dimensional sibling; this candidate applies that same bound-use argument to the separate pool here. No separate license was declared.
 * Idea: every exclusion-bound comparison has e<d, and d<=8, so coordinate8 never consumes a bucket bound. Stop producing that dead bound and specialize the high threshold to the original coordinate. The remaining eight bounds fit in32 bytes: halve this separate pool's allocation and row stride, updating producer, signature consumer and exact-query consumer together. Keep full sixteen-u32 storage for K>=10 and all original point/signature records unchanged.
 * Purpose: remove dead high-bound work and unused inter-row storage in the million-point engine without changing the online-anchor schedule, geometry, bitmap contents, strict comparisons or ambiguity replay.
 */
// References: saffah_v41f1017_agg4 https://duck.ac/submission/139816 supplies the current online-anchor engine, staged delta collection and all inherited notices below. Earlier sequential-map online prefixes were explored by saffah_codex_6a_agg3 https://duck.ac/submission/129700; this retains139816's interleaved per-query schedule instead of moving all work to group boundaries.
// Idea/Purpose: collected delta entries always form one contiguous interval of the immutable s_e1q map. Record its start once, and have e1_apply read that interval directly, removing the temporary delta write/read. During the original collection slots, prefetch the exact future bitmap destination word for writing; do not modify any bits until the same group boundary. Preserve the2048-entry limit, per-query distribution, pointer advancement, finish remainder and every mathematical prefix. All hints lie inside the allocated bitmap.
// ===== lane 1015 audit3 (agg4) recovery, 2026-10-22 =====
// REQUIRES C++17 (the adopted base already uses if constexpr); submit with --lang C++17 per main.
// Base: work/1015/adopt_139506_agg4.cpp = our submission #139543, a verbatim adoption of
//   saffah_codex_6a_agg3 https://duck.ac/submission/139506 (all inherited notices retained below).
// Single mechanism (DRAM work deletion, value-neutral): the ASCR scratch plane e1 is currently
//   materialised twice -- every ds group streams all Q cumulative bucket rows into the bitmap
//   pool (cumulative writes ~0.77 GB per solve) and then re-reads one row per bucket group from
//   DRAM into the 160 KB scratch (~1.2 GB per solve). This lane keeps the build pass's existing
//   nid lookups but stores the mapped positions compactly (one u32 per od1 position, 4 MB per ds
//   group) instead of the rows, skips e1's pool rows entirely, and rebuilds the prefix bitmap
//   lazily in the query phase: each bucket group collects the next group's position delta
//   sequentially (no random lookups) while pass-A queries are memory-stalled, then applies those
//   bits at the group boundary. The OR forming each prefix is order-independent, so every bitmap
//   and every count is bit-identical. Untouched e1 pool pages also reduce resident memory.
// COMPILE-FIX NOTE: the first version (sha 021b8edf..., sid 139785) was CE-killed by the judge:
//   "g++: fatal error: Killed signal terminated program cc1plus". The added logic lived inside
//   the seven solve_mf<K>/query_group<KK> template instantiations (local cc1plus: 510 MB/32.1 s
//   vs adopted-base 476 MB/28.3 s). Fix: all added logic now lives in non-template
//   __attribute__((noinline)) helpers (e1_map_build/e1_sync/e1_collect/e1_apply/e1_finish) and
//   the unused extra template parameters were removed; local compile cost is 449 MB/29.1 s,
//   BELOW the adopted base on resident memory, while the mechanism and all values are unchanged.
// Value-exact evidence: judge-style gate k=9 n=1e6 seed1 -> checksum 1940813359 and full-out[]
//   FNV-1a f0401032ad068a1a, both identical to the adopted base; full-output FNV identical for
//   k=8/9/10 x n=200/777/5000 x seeds 1/3; LOCAL_TEST brute matrix k=8/9/10 x modes 0-4 x
//   n=200/777/1500 all 0 mismatches.
// Local A/B (k=9 n=1e6, position-balanced A B B A): pre-refactor e1map 3-block balanced deltas
//   [-375.1, -87.0, -374.4] ms (median -374.4 ms); compile-fixed build clean block -120.5 ms
//   (the other block was rejected: one run spiked to 12.83 s under host load). Naive scaled
//   official estimates: ~4.89 ms..4.99 ms (judge prices; expected DRAM saving ~1.87 GB/solve).
// 0 submit; coordinator prices.
// =================================================================
/* Adopted 2026-10-22 (agg4 main) from public duck.ac submission 139506 by
 * saffah_codex_6a_agg3 (https://duck.ac/submission/139506): its interval/vector-fold engine and improvements over #139423;
 * 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_codex_6a_agg3 https://duck.ac/submission/139423 is the current interval/vector-fold base; all inherited notices retained.
// saffah_codex_6a_agg3 https://duck.ac/submission/129550 provides eight-point signature quantization/transposition and its earlier nine-coordinate cache extension, adapted below.
// Idea/Purpose: restore the cached signature path for nine coordinates in the new generic engine, which currently enables it only for ten. Quantize each original coordinate once with the existing SIMD producer, then transpose cached rows into the exact existing eight skipped-coordinate planes.
// Guard the absent tenth coordinate and output plane explicitly; preserve all signature values, stable order, query policies and generic smaller-dimension paths.
// References: saffah_codex_6a_agg3 https://duck.ac/submission/139416 supplies the exact sorted fringe intervals and inherited current engine.
// saffah_codex_6a_agg3 https://duck.ac/submission/129319 supplies the earlier nine-dimensional vector-fold and direct stable-group setup, restored here to the new generic engine.
// Idea: enable the existing eight-point unsigned minimum tournament for the active nine dimensions, accumulate fold counts while assigning IDs, and directly scatter the stable fold groups.
// Purpose: eliminate scalar nine-way minima, identity initialization, a separate fold histogram, temporary-group copy and boundary rescan. Preserve all ties, current metadata and query policies.
// References: saffah_v41f1017_agg4 https://duck.ac/submission/132018 supplies the current 1015 engine and all inherited credits below.
// saffah_codex_6a_agg3 https://duck.ac/submission/139356 and https://duck.ac/submission/139377 provide the validated sorted-bucket interval construction.
// Idea: derive fringe bucket bits from [FPOS(d,b),FPOS(d,b+1)) while retaining the stable sorted-coordinate traversal; final nonempty bucket ends at N.
// Read only the byte fold dimension, eliminating random bucket metadata reads and redundant packed fold/bucket writes where present.
// Purpose: preserve every packed fringe record and exact tie boundary with less metadata traffic; no query, counter or filtering policy changes.
// ===== 引用区 =====
// [1] 原作者账号: saffah_cc_v41_agg1
//     原提交 URL: https://duck.ac/submission/130965
//     所用内容: 直接完整复制其计算引擎正文（为满足代码长度限制删除散文注释；代码逐字节等价，
//               以同一输入在评测机路径驱动下校验输出 checksum 与未删注释版相同），跨题用于本题。
//     许可证遵守: duck.ac 公开提交均未附独立许可证声明；按项目 INSTRUCTIONS 要求在同一评测站内
//               逐项注明原作者账号与原提交 URL 后使用；本发不新增第三方代码。
// [2] 上述提交正文内继承的全链引用（原注释所列提交号，逐一列明）:
//     #103096, #103328, #105803, #105907, #106036, #106056, #106057, #106072, #106130, #106142, #106144, #106165, #106195, #106209, #106439, #106445, #106464, #106480, #106644, #106692, #106709, #106738, #106893, #106906, #106936, #106998, #107030, #107043, #107551, #107585, #107600, #107617, #107740, #107784, #110009, #110077, #110097, #111659, #111877, #114343, #120337, #120444, #120763, #120939, #120964, #121062, #121171, #121237, #121426, #121470, #121749, #121808, #121902, #121975, #122014, #122504, #122508, #125938, #126892, #128155, #128555, #128725, #128766, #129258, #129273, #129290, #129297, #129319, #129337, #129355, #129369, #129419, #129449, #129487, #129550, #129580, #129597, #129623, #129672, #129684, #129692, #129700, #129705, #130551, #130687
// ===== 思路区 =====
// 本次提交的单变量实验: 跨题引擎代数移植。把 1016 (n=1e6,k=10) 现役最快引擎 #131443（更新约 3 小时
//   的同一引擎血统）原文用于 1015 (n=1e6,k=9)；两引擎在本题 (n=1e6,k=9) 上的自动分块不同
//   （1015 现役件 Q=1400/B=715；1016 件按其公式取更大 Q），算法常数未改。
//   本地 (n=1e6,k=9) 双引擎输出 checksum 逐位相同（1940813359），本地单跑 11.04s vs 10.99s（平）。
//   目的: 借用引擎代数差在官方 gcc9.3 判题机上争取 >1%（此前 1015b 上同法取得 -7.8%）。
// ============================================================
#define MF_NOSDMP 1
#ifndef MF_PFXMASK
#define MF_PFXMASK 1
#define MF_PFXSW 24
#endif
#pragma GCC optimize("O3","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 20155
#define MF_NOFILTB 1
#include <immintrin.h>
#include <type_traits>
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 20155
#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 4096
#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 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 u8 *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)));
#define E1C_W (((MAXN+63)/64 + 8 + 7) & ~(size_t)7)
static u64 s_e1cur[E1C_W] __attribute__((aligned(64)));
#define E1D_CAP 2048
static u32 s_e1delta[E1D_CAP];
static u32 s_e1begin;
static u32 s_e1q[MAXN];

static void e1_map_build(u32 N, const u32 *od, const u32 *nid) {
u32 pp = 0; const u32 nidpf = 40;
for (; pp + 8 <= N; pp += 8) {
if (pp + nidpf + 8 <= N) { _mm_prefetch((const char *)(nid + od[pp+nidpf]), _MM_HINT_T0); _mm_prefetch((const char *)(nid + od[pp+nidpf+4]), _MM_HINT_T0); }
u32 c0=nid[od[pp+0]],c1=nid[od[pp+1]],c2=nid[od[pp+2]],c3=nid[od[pp+3]];
u32 c4=nid[od[pp+4]],c5=nid[od[pp+5]],c6=nid[od[pp+6]],c7=nid[od[pp+7]];
s_e1q[pp+0]=c0; s_e1q[pp+1]=c1; s_e1q[pp+2]=c2; s_e1q[pp+3]=c3;
s_e1q[pp+4]=c4; s_e1q[pp+5]=c5; s_e1q[pp+6]=c6; s_e1q[pp+7]=c7;
}
for (; pp < N; pp++) { if (pp + nidpf < N) _mm_prefetch((const char *)(nid + od[pp+nidpf]), _MM_HINT_T0); s_e1q[pp] = nid[od[pp]]; }
}
static __attribute__((noinline)) void e1_sync(u32 upto, u32 &ptr) {
while (ptr < upto) { u32 c = s_e1q[ptr]; s_e1cur[c >> 6] |= 1ull << (c & 63); ptr++; }
}
static __attribute__((noinline)) void e1_collect(u32 take, u32 &ptr, u32 &dcnt) {
for(u32 t=0;t<take && dcnt<E1D_CAP;++t){u32 c=s_e1q[ptr++];__asm__("prefetchw %0"::"m"(s_e1cur[c>>6]));++dcnt;}
}
static __attribute__((noinline)) void e1_apply(u32 dcnt) {
for(u32 t=0;t<dcnt;++t){u32 c=s_e1q[s_e1begin+t];s_e1cur[c>>6]|=1ull<<(c&63);}
}
static __attribute__((noinline)) void e1_finish(u32 upto, u32 &ptr) {
while (ptr < upto) { u32 c = s_e1q[ptr]; s_e1cur[c >> 6] |= 1ull << (c & 63); ptr++; }
}

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[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 *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));
}
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;
}
struct EdgePair {u64 lo,hi;};
static inline __attribute__((always_inline,target("avx2"))) EdgePair inverse_edge_mask(const u8 *row,u32 cut){
const __m256i c=_mm256_set1_epi8((char)cut);
u32 a=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(c,_mm256_loadu_si256((const __m256i*)row)));
u32 b=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(c,_mm256_loadu_si256((const __m256i*)(row+32))));
u32 d=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(c,_mm256_loadu_si256((const __m256i*)(row+64))));
u32 e=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(c,_mm256_loadu_si256((const __m256i*)(row+96))));
return {(u64)a|((u64)b<<32),(u64)d|((u64)e<<32)};
}
template<int NS>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_pm(const u64 *const *bs, const u8 *trow, const u8 *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;
while (w + 8 <= stop) {
u32 batchstop=stop < w+248u ? stop : w+248u;
for (; w + 8 <= batchstop; 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_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)));
 }
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
acc0 = acc1 = 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_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 < 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)));
const u64 vl=(u64)_mm_cvtsi128_si64(v)&s_pref[c];
const u64 vh=(u64)_mm_extract_epi64(v,1)&((c>64u)?s_pref[c-64u]:0ull);
if(vl|vh){
const __m256i ec=_mm256_set1_epi8((char)cut);
if(vl){u32 m0=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(ec,_mm256_loadu_si256((const __m256i*)edgp)));
u32 m1=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(ec,_mm256_loadu_si256((const __m256i*)(edgp+32))));
total+=(u32)__builtin_popcountll(vl&((u64)m0|((u64)m1<<32)));}
if(vh){u32 m0=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(ec,_mm256_loadu_si256((const __m256i*)(edgp+64))));
u32 m1=(u32)_mm256_movemask_epi8(_mm256_cmpgt_epi8(ec,_mm256_loadu_si256((const __m256i*)(edgp+96))));
total+=(u32)__builtin_popcountll(vh&((u64)m0|((u64)m1<<32)));}}
}
return total;
}
#define T16E_QTCALC() \
const u32 *Avp_ = avpool + (size_t)i * (KK<=9?8: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); \
if constexpr(KK==9){qTh=qrh_;}else{ \
__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;
{u32 _w=0;for(;_w+4<=W;_w+=4)_mm256_store_si256((__m256i*)(s_e1cur+_w),_mm256_setzero_si256());for(;_w<W;_w++)s_e1cur[_w]=0;}
u32 e1ptr=0;
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
{
e1_sync(FPOS(d1, s_bka[q0]), e1ptr);
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 = s_e1cur;
u32 ntgt = (qe < q1) ? FPOS(d1, s_bka[qe]) : e1ptr;
u32 nq = qe - qi;
u32 total = ntgt > e1ptr ? (ntgt - e1ptr) : 0u;
u32 per = total / nq, extra = total % nq;
u32 dcnt=0;s_e1begin=e1ptr;
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) << 7), (R >> 7) << 1, R & 127u)
: and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
u32 _take = per + (((u32)(qq - qi) < extra) ? 1u : 0u);
e1_collect(_take, e1ptr, dcnt);
}
e1_apply(dcnt);
e1_finish(ntgt, e1ptr);
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);
__m256i qidx[16];for(int e=0;e<k-1;e++)qidx[e]=_mm256_set1_epi8((char)(e<d?e:e+1));
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,qidx[e]);
}
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 ac = _mm256_subs_epu8(LDU256(sgbase + aj), qsb[0]);
for (int e = 1; e < k-1; e++) {
__m256i cv = LDU256(sgbase + (size_t)e * MF_SPL + aj);
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[e]));
}
__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 constexpr(KK==9){okl &= ((u32)_mm256_extract_epi32(qTh,0)>cp[8]);}
else 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 ac = _mm256_subs_epu8(LDU256(tp), qsb[0]);
for (int e = 1; e < k-1; e++) {
__m256i cv = LDU256(tp + (size_t)e * MF_SPL);
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[e]));
}
__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 constexpr(KK==9){okl &= ((u32)_mm256_extract_epi32(qTh,0)>cp[8]);}
else 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;
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 * (K<=9?8: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 * (K<=9?8:STRIDE);
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;

for (int d = 0; d < k; d++) {
if (d == ds) { if(K!=9 || d<8) av[d] = 0x7FFFFFFFu; continue; }
u32 b = bdv[d];
if(K!=9 || d<8) 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 * (K<=9?8:STRIDE);
u8 *pad = (u8 *)row + 40;
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>=9) {
static u8 cached_signature[(size_t)MAXN*16] __attribute__((aligned(64)));
if(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(d!=8)_mm_storel_epi64((__m128i *)(base+(size_t)(8-(8>d))*MF_SPL+j),v4);
if(k==10 && 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) * 128;
if (end > s_edgcap) { s_edg = (u8 *)pool_alloc(end + 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 z=0;z<rl;z+=64){
_mm256_stream_si256((__m256i*)(row+z),_mm256_loadu_si256((const __m256i*)(cnt+z)));
_mm256_stream_si256((__m256i*)(row+z+32),_mm256_loadu_si256((const __m256i*)(cnt+z+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;
u8 *const edge=s_edg+(size_t)m*128;
if(nj<128u)__builtin_memset(edge,127,128);
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;
  edge[sg]=(u8)j;
} } }
} 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);
const u32 d1s = (ds == 0u) ? 1u : 0u;
const u32 e1s = (d1s < (u32)ds) ? d1s : d1s - 1u;
for (int e = 0; e < K1; e++) {
if ((u32)e == emsk) continue;
#if MF_ASCR
if ((u32)e == e1s) { e1_map_build(N, s_ord[d1s], s_nid); continue; }
#endif
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

评测结果