代码

/* Adopted 2026-11-04 by account saffah_v41f1017_agg4 (campaign seat agg4) from public
 * duck.ac submission 141698 (https://duck.ac/submission/141698, account
 * saffah_codex_6a_agg3, Accepted 336.892992 ms, C++17, problem 1014a). The complete
 * source and every inherited citation/license notice below are retained verbatim.
 * That public duck.ac submission declares no separate license; it is credited by
 * account name and submission URL. Used: the whole submission as the base.
 * Changed, two value-neutral edits in the BANDFUSE far loop:
 *  (1) Next-query mask-row prefetch: for each far record the seven 128-byte band
 *      table rows of the next record are prefetched to L1 at distance one query
 *      (both 64-byte lines each, fourteen T0 prefetches). The inherited code kept
 *      this prefetch only under !MF_BANDFUSE, where it targeted the deleted snapshot
 *      pool; the band-local table is L2-resident and the far loop's dependent
 *      rec-load -> mask-load -> serial-AND chain still exposes that latency.
 *  (2) Removed the far loop's zero-intersection early exit (testz + continue after
 *      the first four dims); the remaining three masks are now always ANDed. This is
 *      output-identical: if the four-dim partial AND is zero, ANDing further masks
 *      keeps it zero, so its popcount contribution is zero either way. The branch
 *      fired in only 0.8% of lane-steps and cost more issue slots than it saved.
 * No other change; the body carries our own https://duck.ac/submission/141115
 * band-local engine plus the rival's carried-prefix construction and Q 228->254.
 * Purpose/approach: hide L2 mask-row latency and delete a near-dead branch in the
 * dominant far phase; rebuilt with the judge compiler and verified value-exact
 * before official timing (348 shapes byte-identical, registered anchor n=300000
 * seed=1 mode=0 checksum 351867771). */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/141687 supplies the locally built1024-band engine with carried prefix vectors; saffah_v41f1017_agg4 https://duck.ac/submission/141115 supplies its band-local algorithm. Our older https://duck.ac/submission/140321 used254 buckets with a full snapshot pool; that obsolete pool is not used here. All inherited notices remain; public sources declare no separate license.
Idea: after local-band construction removed the full snapshot pool, increase the geometry coefficient36 to40, yielding254 buckets at300000 points. The seven-plane live table remains under256KiB (228480 bytes), while finer buckets shorten fringe scans. Retain carried prefixes,1024-rank bands, exact ties and every counting operation.
Purpose: measure the changed construction/fringe crossover in the new local-table architecture; full-output and boundary validation precedes official timing. */
/* References: saffah_v41f1017_agg4 https://duck.ac/submission/141115 supplies the complete1024-rank band-local eight-dimensional bitmap engine, used in full below; its lineage includes saffah_codex_6a_agg3 https://duck.ac/submission/140330. All inherited source/license notices remain. These public Duck sources declare no separate license.
Idea: hold the four cumulative256-bit vectors in registers across each bucket's prefix construction. The zeroed row0 receives no inserted points because insertion indexes bucket+1, so the initial carries are zero; each later row is the exact OR of its delta and the retained predecessor. Eliminate four predecessor-row reloads per bucket, retaining every written slice and all query/fringe arithmetic.
Purpose: evaluate genuine load deletion on the newly adopted local-band producer; validate complete outputs, ties, guards, hashes and independent samples before official timing. */
/* Adopted 2026-10-31 by campaign seat agg4 from public duck.ac submission 140330
 * (https://duck.ac/submission/140330, saffah_codex_6a_agg3, Accepted 404.616890 ms,
 * C++17, problem 1014a); the complete source and every inherited citation/license
 * notice below are retained verbatim. The cited public submission is a duck.ac site
 * submission and declares no separate license; it is credited by account and
 * submission URL. Used: the whole submission as the base of this file; the only
 * change added here is the BANDFUSE cut described below.
 *
 * BANDFUSE cut (single mechanism, value-exact): delete the band-snapshot pool
 * producer/consumer round trip of the far phase. The pool s_bs (per (dim,bucket)
 * prefix bitmap rows, rebuilt for each of the 8 anchor groups) is no longer
 * materialised with a scattered RMW pass plus a streamed NT copy and re-read as
 * band slices. Instead, for each 1024-rank band and each anchor group, the
 * 7x(Q+1) band-slice table (7*(228+1)*16 words = 200.4 KB, <=256 KB) is rebuilt
 * in s_btab directly from the band's 1024 points: each point's dim bucket scatters
 * one bit into row bucket+1, then one prefix-OR pass per dim makes row b the exact
 * {bucket < b} slice that the old pool row slice held. The far and partial-band
 * consumers read the same bit-identical slices from the L2-resident table; the
 * fused path also drops the 14 per-iteration next-query pool prefetches (their
 * target pool no longer exists; rows are L2-resident), following the proven
 * BANDFUSE1024 artifact of work/1015b (sid 141104). All answers, strict tie
 * boundaries, signature screens, fold selection, record order, fringe traversal
 * and output formulas are unchanged. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140321 supplies bounded-byte survivor counting and every inherited notice below. saffah_v41f1017_agg4 https://duck.ac/submission/139842 supplies coefficient36 for the dimension-major eight-dimensional engine; its original source and attribution are already inherited below. No separate license was declared; notices are retained.
 * Idea: combine the newly cheaper signature scan with the cited228-bucket geometry rather than the254-bucket variant in our prior submission. Reducing prefix snapshot construction trades for more signature work, whose per-block SAD was just eliminated. Change only the coefficient40 to36; retain the original correct tie boundaries, live prefix construction, byte-overflow chunking and exact query predicates.
 * Purpose: evaluate the geometry balance after a real signature-cost change, combining the external best geometry with our new reductions instead of repeating either historical implementation.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140318 supplies the current static eight-dimensional fringe engine, including its original per-vector survivor SAD and all inherited source/license notices below. No separate license was declared; attribution remains intact.
 * Idea: count each definite-survivor byte as one by subtracting its all-ones equality mask into a byte accumulator. Process at most255 vectors per chunk, so every lane is exactly0..255 with no wrap, then widen that chunk once using SAD. Arbitrarily long ranges use further chunks. This removes SAD and 64-bit addition from every signature iteration, plus the final multiply that previously divided255-valued counts. Ambiguity classification/replay, tail masks and geometric boundaries stay identical.
 * Purpose: experimentally replace repeated horizontal reductions with bounded vertical accumulation; no timing claims rely on local execution.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/139895 supplies this eight-dimensional global-fringe engine and all inherited notices below. Our https://duck.ac/submission/140277 supplies the statically inactive first-dimension bound elimination in the seven-dimensional sibling. No separate license was declared; original notices remain intact.
 * Idea: for compile-time fringe dimension zero, exmask is zero and the bound vector is ORed with INT_MAX in every lane. Coordinates are nonnegative ranks below300000, so min(qr,INT_MAX) equals qr. Omit this dimension's unused bound vector read and minimum. Keep its queries, bitmaps, signatures, exact comparisons and every other dimension unchanged.
 * Purpose: experimentally reduce dead query-setup work in this distinct eight-dimensional bitmap-band engine.
 */
// References: saffah_v41f1017_agg4 https://duck.ac/submission/139842 supplies its36-coefficient geometry and lookahead on saffah_codex_6a_agg3 https://duck.ac/submission/139740's dimension-major fringe engine; inherited credits remain below. Earlier40-coefficient geometry on the different per-anchor fringe engine was explored by saffah_codex_6a_agg3 in the chain of https://duck.ac/submission/121022; that old traversal is not reintroduced.
// Idea/Purpose: test coefficient40 specifically after replacing per-anchor fringe sorting with a single monotone dimension traversal and clustered band records. Relative to139842's36, this increases actual live buckets from228 to254 at300000 points, shortening all same-bucket signature scans while increasing prefix snapshots. Keep the inherited255 cap, every strict tie boundary, all cached signatures, exact ambiguity checks, hints and output formulas unchanged. This is an official geometry tradeoff for the new traversal, not an identical old-grid resubmission.
/* Further optimized 2026-10-22 from the adopted public submission below.
 * Named citation: saffah_codex_6a_agg3, https://duck.ac/submission/139740
 * (Accepted 419.514288 ms, C++17, 1014a); its complete source and every
 * inherited reference/license notice below are retained verbatim. The cited
 * public submission states no separate license. This file is not that
 * submission: it makes three value-preserving
 * tuning changes, all checked through the site's C++17 judge-machine
 * experience channel and by byte-exact output comparison:
 *  1. MF_QCOEF 32 -> 36 (automatic Q 203 -> 228 at official n=300000).
 *     The global dimension-major fringe supplied by the cited lineage scans
 *     about 4*N*N/Q same-bucket signature slots, while pool construction grows
 *     with Q; after that traversal replaced per-dimension fringe sorting, the
 *     inherited coefficient is no longer at the judge-machine balance point.
 *  2. MF_PFB3 4 -> 8: that same traversal walks ranks monotonically inside each
 *     tie-aware bucket, so its candidate lookahead is retuned.
 *  3. Delete the inherited candidate prefetch of s_bd[j] in
 *     dimension_major_fringe. The new traversal derives each bucket's lower
 *     bound from FPOS(d,bucket) and never reads s_bd for a candidate j, so this
 *     random 300000-row prefetch was dead work.
 * Answers, tie handling, ordering, counters, thresholds, all other prefetches,
 * all references and all inherited notices are unchanged. */
/* Adopted 2026-10-22 from public duck.ac submission 139740 by saffah_codex_6a_agg3
 * (https://duck.ac/submission/139740, Accepted 419.514288 ms, C++17, 1014a); complete
 * source and every inherited reference/license notice below are retained verbatim.
 * This adoption changes no code; public submission, no separate license stated. */
// References: saffah_codex_6a_agg3 https://duck.ac/submission/139703 supplies first-plane-clustered sparse complete bands and all inherited notices below. saffah_v41f1017_agg4 https://duck.ac/submission/139705 supplies global coordinate-order fringe traversal on1016a.
// Idea/Purpose: eliminate this8D engine's repeated per-anchor/per-dimension histogram, scatter and copy-back of fringe query IDs. After full-band answers finish, traverse each existing coordinate order by bucket, skip queries anchored in that dimension, and reuse the same compile-time dimension signature kernels. Set each point's previously ignored anchor Av lane toINT_MAX once so dimension-only masks preserve every threshold. Other Av lanes, all prefix records, exact replay and output values stay unchanged.
// This changes traversal and removes intermediary sorting work; it is separate from the prior sparse-band gate. Stable coordinate ranks delimit the same tie-aware buckets, with the last nonempty bucket ending at N.
// References: saffah_codex_6a_agg3 https://duck.ac/submission/139655 supplies first-plane-clustered rank bands with all inherited credits; https://duck.ac/submission/139356 and its original cited1016 sources provide first-four bitmap early rejection.
// Idea/Purpose: test the same exact zero-intersection gate in this different seven-stream1024-bit band engine. Compute the original first four streams for each512-bit slice; when their entire intersection is zero, later ANDs and popcount are provably zero, so skip those demand loads and arithmetic. No stream order or prefetch address changes.
// Preserve partial own-band queries and all record order, counters, full-mask semantics and output. This is a sparse-data branch/load tradeoff measured only by the official submission, not an assumed speedup.
// References: saffah_codex_6a_agg3 https://duck.ac/submission/139584 supplies the current eight-coordinate live-prefix/fold-maximum engine, with all inherited notices below; https://duck.ac/submission/134015 supplies first-plane bucket clustering within rank bands on1016a.
// Idea/Purpose: adapt that source-reuse ordering to this1024-bit band engine. Count-sort by (R/1024,first non-fold bucket) so adjacent queries reuse the same first bitmap slice; retain exact R and point ID in every record.
// Every query consumes complete earlier bands and its exact partial own band, including zero partials at1024 boundaries. Reuse the dead s_bandacc array for the smaller histogram; ceil(MAXN/1024)*(Q+1)<=ceil(1000005/1024)*256<MAXN, and reset accumulators after scattering.
// This also replaces the full N-key initial rank histogram with the joint band/bucket key domain. No bitmap contents, fringe traversal, pool layout or numerical formulas change.
// References: saffah_v41f1017_agg4 https://duck.ac/submission/134673 supplies the current 1014a engine; all inherited attribution/license notices remain below.
// saffah_codex_6a_agg3 https://duck.ac/submission/139559 provides retained stable bitmap positions; https://duck.ac/submission/139573 and https://duck.ac/submission/139566 provide live-prefix-only construction and fold maximum retention.
// Idea/Purpose: actual buckets always lie in [0,Q), so stop constructing the unconsumed terminal snapshot Q. Preserve allocation offsets, all live snapshots, arithmetic, ordering and fallback behavior.
// Count each fold maximum while assigning points, replacing only the subsequent Wpool rank rescan. Store the dimension-only s_g in bytes (values0..9) to reduce group-building traffic; no packed rank representation is changed.
// ===== 引用区 =====
// [1] 本账号 saffah_v41f1017_agg4，提交 #132635（472.482 ms）—— 本件正文的底盘：#131038 引擎
//     + 旗标消融组合；其自带的完整嵌套引用链（#131038 及下列提交号）原样保留。
// [2] 对手账号 saffah_codex_6a_agg3，公开提交 #133921 <https://duck.ac/submission/133921>
//     （469.014 ms）—— 所用内容：setup 相融合（每 8 点 AVX2 min_d(rank<<4|d) 折叠选择、
//     组内计数、一次前缀、cursor 原地稳定散布）。本账号已于 #133949 移植并官方实测 468.986 ms。
// [3] 对手账号 saffah_codex_6a_agg3，公开提交 #133933 <https://duck.ac/submission/133933>
//     （467.724 ms，2026-10-18 03:10:02）—— **本发唯一新增改动来源**：签名缓存的 SIMD 量化
//     与字节转置（每 8 点、8 维：mullo+srl23+max1+pshufb+unpack 转置、4×16B 存，取代
//     64 次标量乘/存；尾块标量兜底）。
// [4] 继承链（#131038 正文原列）: #104600, #105554, #108160, #108803, #108836, #109857,
//     #111877, #112420, #114041, #114076, #114343, #120185, #120231, #120336, #120475,
//     #120507, #120554, #120788, #120919, #120922, #120949, #120953, #121022, #121062,
//     #121092, #121378, #121423, #121504, #121549, #121604, #121720, #122354, #122473,
//     #126043, #126050, #126060, #126875, #127013, #127053, #128242, #128317, #128389,
//     #128441, #128447, #128501, #128525, #128537, #128686
// [5] 本账号 saffah_v41f1017_agg4，提交 #133969 <https://duck.ac/submission/133969>
//     （467.718 ms，2026-10-18 03:18:17）—— 本件底盘与本发唯一变量载体（far 相预取补丁
//     直接落在其 far 向量循环内）；其自带引用链 [1]–[4] 原样保留，无新增外部来源。
// LICENSE：均为 duck.ac 站内公开提交，无第三方许可证声明；按「账号/提交地址/所用内容」具名引用。
// ===== 思路区 =====
// 【本发性质】移植件（在 #133949 之上叠加单变量，输出逐位不变）。
//   签名缓存构建原为每点每维一次标量 64 位乘 + 移位 + 字节存（N*8 = 2.4e6 次）。
//   改为每 8 点一组：8 条 _mm256_mullo_epi32（g_sigmul 乘、src>>23、max 1）经 pshufb 取字节，
//   再用 unpacklo/hi 链做 8×8 字节转置，4 条 16B 对齐存落 cached_signature；
//   尾块（<8 点）保留标量 sig_nz。算术逐位等价（乘积 < 256*2^23 < 2^32 ⇒ mullo 等于 64 位低半）。
// 【值闸门】判题同款 gcc 9.3.0：本地 n=200 × seed{1,2} × mode{0,1,3,4} 对暴力 oracle 全 OK；
//   n=300000 mode=10 checksum=351867771 与底盘 #132635 逐位相同；同机 A/B min 923→922 ms。
// 【判靶】严支线 464.325 ms（0.99T+1µs，T=467.724）；本发预计 ~467.7 ms（取回对手 −0.27%）。
// 【本发性质】自研优化（在 #133969 之上叠加单变量：far 相下一查询掩码切片预取；输出逐位不变）。
//   far 相账（现役源 in-engine）：9,619,861 query-band 迭代 × 896B（7 维 × 2 lane × 64B），
//   实测 57 拍/迭代 = 4.1× 载入端口地板；15.2B/拍 只合 54.6GB/s 而孤核 DRAM 峰 ~6B/拍
//   ⇒ 逐查询随机桶偏移下的载入延迟/MLP 界（上一席已证非 DRAM、非端口界）。
//   改动：far 向量循环内用 rec 下一行（rec + (t+1)*(KK+1)）的 7 个维桶偏移 rq[1+e]，
//   对 s_bs + (off8<<3) + rq[1+e] 的 ±两条 64B 行各发 _mm_prefetch(T0)（每迭代 14 条），
//   即把下一查询的 14 行请求提前一个迭代（~57 拍）放进在飞窗口；plen 布局下桶偏移均为
//   64B 对齐（len 8 字对齐），128B 切片恰为两行，无覆盖缺口。语义零改动（只读 rec 下一行）。
// 【值闸门】判题同款 gcc 9.3.0 本地：n=300000 全变体 checksum 355887373 与底盘逐位同；
//   判题机 experience checksum 356972483 与底盘 ph300 逐位同（判题 libc 的 rand 序列与本地
//   厂商链不同，故判题侧 checksum 数值与本地不同，两侧各自恒等）。
// 【判题机 A/B】（experience，同机同源，两轮交错）base 473.81/473.80/473.81 ms
//   vs 本发 d1 462.54/462.56 ms = −11.26 ms（−2.38%）；引擎自计 query 相
//   1,399.5M→1,358.8M ticks（−40.7M）；距离 2 = −10.7 ms、距离 4 = −5.3 ms、
//   单行预取 = +1.3 ms（负）⇒ 「双行 + 距离 1」唯一正点，剂量曲线单调。
// 【判靶】严支线 463.048 ms（0.99T+1µs，T=467.724 #133933）；本发预计 ~456.4 ms（−11.3 ms 移植量）。
// ============================================================
#define MF_NOSDMP 1
#define MF_QCOEF 40
#define MF_BANDFUSE 1
#pragma GCC optimize("O3,no-strict-aliasing","live-range-shrinkage","no-tree-forwprop","no-tree-ch","unroll-all-loops")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,arch=skylake")
#undef __k8
#undef __k8__
#undef __code_model_small__
#undef __MMX_WITH_SSE__
#define __skylake 1
#define __skylake__ 1
#define __tune_skylake__ 1
#define __code_model_small__ 1
#define __SSE3__ 1
#define __SSSE3__ 1
#define __SSE4_1__ 1
#define __SSE4_2__ 1
#define __AES__ 1
#define __PCLMUL__ 1
#define __AVX__ 1
#define __AVX2__ 1
#define __SGX__ 1
#define __FMA__ 1
#define __BMI__ 1
#define __BMI2__ 1
#define __LZCNT__ 1
#define __CRC32__ 1
#define __POPCNT__ 1
#define __FSGSBASE__ 1
#define __RDRND__ 1
#define __F16C__ 1
#define __RDSEED__ 1
#define __PRFCHW__ 1
#define __ADX__ 1
#define __XSAVE__ 1
#define __XSAVEOPT__ 1
#define __CLFLUSHOPT__ 1
#define __XSAVEC__ 1
#define __XSAVES__ 1
#define __MMX_WITH_SSE__ 1
#define __LAHF_SAHF__ 1
#define __MOVBE__ 1
#define MF_NOFILTB 1
#include <immintrin.h>
#include <type_traits>
#include <type_traits>
#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 char u8;
typedef unsigned long long u64;
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#define AND256S(a,pz,off,disp) ({ __m256i _v; __asm__("vpand " #disp "(%2,%3), %1, %0" : "=x"(_v) : "x"(a), "r"(pz), "r"(off)); _v; })
#define AND256M(a,p) ({ __m256i _v; __asm__("vpand %2, %1, %0" : "=x"(_v) : "x"(a), "m"(*(const __m256i *)(p))); _v; })
#define LDU8U16(p) ({ __m256i _v; \
__asm__("vpmovzxbd %1, %0" : "=x"(_v) : "m"(*(const u64 *)(p))); _v; })
#ifndef MAXN
#define MAXN 1000005
#endif
#define MAXK 8
#define STRIDE 16
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 42
#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 1900
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
#define MF_PFB 32
#ifndef MF_PFB1
#ifndef MF_PFOR
#define MF_PFOR 192
#endif
#define MF_PFB1 128
#endif
#ifndef MF_PFB2
#define MF_PFB2 128
#endif
#ifndef MF_PFB3
#define MF_PFB3 8
#endif
#ifndef MF_PFB4
#define MF_PFB4 32
#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_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 256
#endif
#define MF_SIGPL ((((size_t)MAXN + 63u) & ~(size_t)63u) + 0x840u)
static u32 s_ord[MAXK][MAXN];
static u8 s_bd[(size_t)MAXN * MAXK];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u32 s_gt[MAXK][MAXN];
static u8 s_bdv[(size_t)MAXN * MAXK];
static u8 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_bandrec[(size_t)MAXN*(MAXK+1)],s_bandacc[MAXN];
static __m128i s_bandvec[MAXN];
static u32 s_nid[MAXN];
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;
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)));
#if MF_BANDFUSE
static u64 s_btab[(size_t)(MAXK - 1) * (MAXQ + 1) * 16 + 64] __attribute__((aligned(64)));
#endif
static u32 *s_sdmp[MAXK];
static u32 *inc_anchor=0;static size_t inc_cap=0;static bool inc_on=false;
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
#define MAXRS 256u
static u32 s_maxr[(size_t)MAXK * MAXK * MAXRS];
static u8 *s_sig = 0;
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 0;
alignas(32) static const unsigned char q8_bcast_mask[4][32] = {
{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4},
{8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8},
{12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12},
};
static u32 g_sigmul = 1;
static inline u8 sig_nz(u32 v) {u32 z=(u32)(((u64)v*g_sigmul)>>23u);return (u8)(z?z: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)]
extern "C" int madvise(void *, unsigned long, int);
#define MF_MADV_HUGEPAGE 14
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
void *p = (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
madvise(p, bytes, MF_MADV_HUGEPAGE);
return p;
}
__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 + 16 <= (Wr); w += 16) {
{ __asm__ __volatile__("nop;nop;nop;nop;nop;nop" :: "m"(*(const char(*)[32])bs[0]) : "ymm15"); }
_mm_prefetch((const char *)(bs[1] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[1] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 40 + 8), _MM_HINT_T0);
__m256i a, b, ta, tb;
a = _mm256_and_si256(LDU256((bs[0] + w)), LDU256((bs[1] + w)));
b = _mm256_and_si256(LDU256((bs[0] + w + 4)), LDU256((bs[1] + w + 4)));
ta = _mm256_and_si256(LDU256((bs[2] + w)), LDU256((bs[3] + w)));
tb = _mm256_and_si256(LDU256((bs[2] + w + 4)), LDU256((bs[3] + w + 4)));
a = _mm256_and_si256(a, ta);
b = _mm256_and_si256(b, tb);
ta = _mm256_and_si256(LDU256((bs[4] + w)), LDU256((bs[5] + w)));
tb = _mm256_and_si256(LDU256((bs[4] + w + 4)), LDU256((bs[5] + w + 4)));
a = _mm256_and_si256(a, ta);
b = _mm256_and_si256(b, tb);
a = _mm256_and_si256(a, LDU256((bs[6] + w)));
b = _mm256_and_si256(b, LDU256((bs[6] + w + 4)));
__m256i c, d4, tc, td;
c = _mm256_and_si256(LDU256((bs[0] + w + 8)), LDU256((bs[1] + w + 8)));
d4 = _mm256_and_si256(LDU256((bs[0] + w + 12)), LDU256((bs[1] + w + 12)));
tc = _mm256_and_si256(LDU256((bs[2] + w + 8)), LDU256((bs[3] + w + 8)));
td = _mm256_and_si256(LDU256((bs[2] + w + 12)), LDU256((bs[3] + w + 12)));
c = _mm256_and_si256(c, tc);
d4 = _mm256_and_si256(d4, td);
tc = _mm256_and_si256(LDU256((bs[4] + w + 8)), LDU256((bs[5] + w + 8)));
td = _mm256_and_si256(LDU256((bs[4] + w + 12)), LDU256((bs[5] + w + 12)));
c = _mm256_and_si256(c, tc);
d4 = _mm256_and_si256(d4, td);
c = _mm256_and_si256(c, LDU256((bs[6] + w + 8)));
d4 = _mm256_and_si256(d4, LDU256((bs[6] + w + 12)));
if (KK - 1 > 7) {
for (int d = 8; d < KK; d++) {
a = _mm256_and_si256(a, LDU256((bs[d] + w)));
b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
c = _mm256_and_si256(c, LDU256((bs[d] + w + 8)));
d4 = _mm256_and_si256(d4, LDU256((bs[d] + w + 12)));
}
}
{
__m256i lc = _mm256_and_si256(c, lm), hc = _mm256_and_si256(_mm256_srli_epi16(c, 4), lm);
__m256i ld = _mm256_and_si256(d4, lm), hd = _mm256_and_si256(_mm256_srli_epi16(d4, 4), lm);
__m256i la = _mm256_and_si256(a, lm), ha = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lb = _mm256_and_si256(b, lm), hb = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
__m256i x0 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, la), _mm256_shuffle_epi8(lk, ha));
__m256i x1 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lb), _mm256_shuffle_epi8(lk, hb));
__m256i x2 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lc), _mm256_shuffle_epi8(lk, hc));
__m256i x3 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, ld), _mm256_shuffle_epi8(lk, hd));
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(x0, x1));
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(x2, x3));
if (++cnt == 15) {
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 KK>
__attribute__((target("avx2,popcnt")))
__attribute__((flatten))
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, *tmp = s_tmp;
if (!g_skip_and) {
const u32 nq = q1 - q0;
u32 *rec = s_bandrec;
u32 max_query_rank=0;
{
 const u32 first_dim=ds?0u:1u,stride=Q+1u,keys=(((N-1u)>>10)+1u)*stride;
 u32 *cv=s_bandacc;
 for(u32 key=0;key<=keys;++key)cv[key]=0;
 for(u32 t=0;t<nq;++t){u32 i=s_qord[q0+t],R=s_gt[ds][i];
  ++cv[(R>>10)*stride+s_bd[(size_t)i*MAXK+first_dim]+1];
  if(R>max_query_rank)max_query_rank=R;
 }
 {u32 sum=0;for(u32 key=0;key<=keys;++key){u32 c=cv[key];cv[key]=sum;sum+=c;}}
 for(u32 t=0;t<nq;++t){u32 i=s_qord[q0+t],R=s_gt[ds][i];const u8 *bd=s_bd+(size_t)i*MAXK;
  u32 key=(R>>10)*stride+bd[first_dim],slot=cv[key+1]++;u32 *rp=rec+(size_t)slot*(KK+1);
  rp[0]=R;rp[KK]=i;
  for(int e=0;e<K1;++e){int d=e<ds?e:e+1;
#if MF_BANDFUSE
   rp[1+e]=(u32)(((size_t)e*(Q+1)+bd[d])*16u)<<3;
#else
   rp[1+e]=(g_plen_on?s_pofs[(size_t)e*MAXQ+bd[d]]:(u32)(((size_t)e*(Q+1)+bd[d])*W))<<3;
#endif
  }
 }
}
{
u32 *acc = s_bandacc;
for (u32 t = 0; t < nq; t++){ acc[t] = 0; s_bandvec[t]=_mm_setzero_si128(); }
if (nq) {
const u32 maxR = max_query_rank;
const u32 maxband = (maxR + 1023u) >> 10;
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 zz = _mm256_setzero_si256();
u32 ps = 0;
for (u32 bk = 0; bk <= maxband; bk++) {
const u32 lim = bk << 10;
while (ps < nq && (rec[(size_t)ps*(KK+1)]>>10)<bk) ps++;
if (ps >= nq) break;
u32 psp_ = ps;
while (psp_ < nq && (rec[(size_t)psp_*(KK+1)]>>10)==bk) psp_++;
const size_t off8 = (size_t)bk << 4;
#if MF_BANDFUSE
{
const size_t TN = (size_t)(KK - 1) * (Q + 1) * 16u;
for (size_t w = 0; w < TN; w += 4) _mm256_store_si256((__m256i *)(s_btab + w), zz);
u32 jend = lim + 1024u; if (jend > max_query_rank) jend = max_query_rank;
const u32 *const odds = s_ord[ds];
for (u32 j = lim; j < jend; j++) {
const u32 id = odds[j];
const u8 *const bdp = s_bd + (size_t)id * MAXK;
const u32 r = j - lim, wi = r >> 6;
const u64 one = (u64)1 << (r & 63);
for (int e = 0; e < KK - 1; e++) { const int d = (e < (int)ds) ? e : e + 1;
s_btab[((size_t)e * (Q + 1) + bdp[d] + 1) * 16 + wi] |= one; }
}
for (int e = 0; e < KK - 1; e++) {
u64 *const eB = s_btab + (size_t)e * (Q + 1) * 16;
__m256i carry0=zz,carry1=zz,carry2=zz,carry3=zz;
for (u32 b = 1; b <= Q; b++) {
u64 *cp = eB + (size_t)b * 16;
carry0=_mm256_or_si256(carry0,LDU256(cp));carry1=_mm256_or_si256(carry1,LDU256(cp+4));
carry2=_mm256_or_si256(carry2,LDU256(cp+8));carry3=_mm256_or_si256(carry3,LDU256(cp+12));
_mm256_store_si256((__m256i*)cp,carry0);_mm256_store_si256((__m256i*)(cp+4),carry1);
_mm256_store_si256((__m256i*)(cp+8),carry2);_mm256_store_si256((__m256i*)(cp+12),carry3);
}
}
}
#endif
for (u32 t = psp_; t < nq; t++) {
const u32 *rp = rec + (size_t)t * (KK+1);
{ const u32 tq = t + 1;
  if (tq < nq) {
    const u32 *rq = rec + (size_t)tq * (KK+1);
#if MF_BANDFUSE
    const char *pb = (const char *)s_btab;
#else
    const char *pb = (const char *)s_bs + (size_t)(off8 << 3);
#endif
    for (int e = 0; e < K1; e++) {
      const char *pn = pb + rq[1 + e];
      _mm_prefetch(pn, _MM_HINT_T0);
      _mm_prefetch(pn + 64, _MM_HINT_T0);
    }
  }
}
__m256i bytes=zz;
for (u32 lane=0;lane<16;lane+=8) {
#if MF_BANDFUSE
const char *pz = (const char *)(s_btab + lane);
#else
const char *pz = (const char *)(s_bs + off8 + lane);
#endif
__m256i a0 = LDU256(pz + rp[1]), a1 = LDU256(pz + rp[1] + 32);
for(int e=1;e<K1 && e<4;++e){const size_t oe=(size_t)rp[1+e];a0=AND256S(a0,pz,oe,0);a1=AND256S(a1,pz,oe,32);}
if constexpr(K1>4){
 for(int e=4;e<K1;++e){const size_t oe=(size_t)rp[1+e];a0=AND256S(a0,pz,oe,0);a1=AND256S(a1,pz,oe,32);}
}
__m256i lo0 = _mm256_and_si256(a0, lm), hi0 = _mm256_and_si256(_mm256_srli_epi16(a0, 4), lm);
__m256i lo1 = _mm256_and_si256(a1, lm), hi1 = _mm256_and_si256(_mm256_srli_epi16(a1, 4), lm);
__m256i s0 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo0), _mm256_shuffle_epi8(lk, hi0));
__m256i s1 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo1), _mm256_shuffle_epi8(lk, hi1));
bytes = _mm256_add_epi8(bytes,_mm256_add_epi8(s0,s1));
}
__m256i sum64=_mm256_sad_epu8(bytes,zz);
__m128i lanes=_mm_add_epi32(_mm256_castsi256_si128(sum64),_mm256_extracti128_si256(sum64,1));
s_bandvec[t]=_mm_add_epi32(s_bandvec[t],lanes);
}
for (u32 t = ps; t < psp_; t++) {
const u32 *rp = rec + (size_t)t * (KK+1);
const u32 R = rp[0];
const u64 *bp[MAXK];
#if MF_BANDFUSE
for (int e = 0; e < K1; e++) bp[e] = (const u64 *)((const char *)s_btab + (size_t)rp[1 + e]);
#else
for (int e = 0; e < K1; e++) bp[e] = (const u64 *)((const char *)s_bs + (off8 << 3) + (size_t)rp[1 + e]);
#endif
acc[t] += and_popcount_range<KK - 1>(bp, (R >> 6) - (u32)off8, R & 63u);
}
}
for (u32 t = 0; t < nq; t++){__m128i lanes=s_bandvec[t];out[rec[(size_t)t * (KK+1) + KK]] = acc[t]+(u32)_mm_cvtsi128_si32(lanes)+(u32)_mm_extract_epi32(lanes,2);}
}
}
}
}

template<int KK>
__attribute__((target("avx2,popcnt")))
static void dimension_major_fringe(u32 N,u32*out,u32 B){
if(g_skip_fringe)return;
const int k=KK;const u32 KS_C=KK<=8?8u:(u32)KK;
auto fringe_dim=[&](auto dim){ constexpr int d=decltype(dim)::value;
u32 exmask = 0;
for(int e=0;e<d;++e)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];
const u8 *const pbase_sig = s_sig + (size_t)d * (k - 1) * MF_SIGPL;
const u32 *const sdmpd = s_sdmp[d];
const u32 last_bucket=(N-1u)/B;
for(u32 bucket=0;bucket<=last_bucket;++bucket){
const u32 lo=FPOS(d,bucket);
const u32 end=bucket==last_bucket?N:FPOS(d,bucket+1);
for(u32 rank=lo;rank<end;++rank){
const u32 i=odd[rank];
if(s_g[i]==d)continue;
#if MF_PFB > 0
{
u32 qj=rank+MF_PFB3;
if(qj<end){
u32 j=odd[qj];
__asm__ __volatile__("prefetchw %0" ::"m"(*(const char *)(&out[j])));
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_T0);
}
}
#endif
const u32 hi = s_gt[d][i];
if (lo >= hi) continue;
const u32 *rec = s_pm + (size_t)i * STRIDE;
__m256i qr = (KK == 8) ? _mm256_load_si256((const __m256i *)rec)
: _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
__m256i Avraw = _mm256_load_si256((const __m256i *)(rec + 8));
__m256i Avm = _mm256_or_si256(Avraw, notexlo);
__m256i qT;
if constexpr (d == 0) qT = qr;
else qT = _mm256_min_epi32(qr, Avm);
__m256i qTh;
if (k > 8) {
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_set1_epi32(0x7FFFFFFF);
qTh = _mm256_min_epi32(qrh, Avmh);
} else {
qTh = all8;
}
u32 add = 0;
const u32 *c0 = sdmpd ? (sdmpd + (size_t)lo * KS_C) : 0;
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const 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_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (g_sig_on && (hi - lo) >= 2) {
__m256i qsv[16];
{
const __m256i mulv = _mm256_set1_epi32((int)g_sigmul);
const __m256i maxv = _mm256_set1_epi32(0x7FFFFFFF);
__m256i sv = _mm256_srli_epi32(_mm256_mullo_epi32(qT, mulv), 23);
if (k > 8) sv = _mm256_or_si256(sv, _mm256_cmpeq_epi32(qT, maxv));
if (k == 8) sv = _mm256_subs_epu8(sv, _mm256_set1_epi8(1));
const __m256i lohalf = _mm256_broadcastsi128_si256(_mm256_castsi256_si128(sv));
const __m256i hihalf = _mm256_broadcastsi128_si256(_mm256_extracti128_si256(sv, 1));
if (k == 8) {
switch (d) {
case 0:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 1:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 2:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 3:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 4:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 5:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 6:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 7:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
break;
}
} else {
int t = 0;
for (int e = 0; e < k; e++) {
if (e == (int)d) continue;
const __m256i src = (e < 4) ? lohalf : hihalf;
qsv[t++] = _mm256_subs_epu8(
_mm256_shuffle_epi8(src, _mm256_load_si256((const __m256i *)q8_bcast_mask[e & 3])),
_mm256_set1_epi8(1));
}
}
}
u32 aj = lo, aacc = 0, aexact = 0;
__m256i dec64 = _mm256_setzero_si256();
if (hi - lo >= 32u) { const u32 ajlim=hi-32u;
while (aj <= ajlim) {
const u32 chunklim = (ajlim < aj + 254u * 32u ? ajlim : aj + 254u * 32u);
__m256i dec8 = _mm256_setzero_si256();
do {
__m256i ac0 = _mm256_setzero_si256(), ac1 = ac0;
for (int t = 0; t < KK - 1; t++) {
__m256i cv = LDU256((const u8 *)(pbase_sig + (size_t)t * MF_SIGPL + aj));
__m256i df = _mm256_subs_epu8(cv, qsv[t]);
if (t & 1) ac1 = _mm256_or_si256(ac1, df);
else ac0 = _mm256_or_si256(ac0, df);
}
__m256i ac = _mm256_or_si256(ac0, ac1);
dec8 = _mm256_sub_epi8(dec8, _mm256_cmpeq_epi8(ac, _mm256_setzero_si256()));
u32 amb=(u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)));
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
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_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
aj += 32u;
} while (aj <= chunklim);
dec64 = _mm256_add_epi64(dec64, _mm256_sad_epu8(dec8, _mm256_setzero_si256()));
}
}
{
__m128i s64 = _mm_add_epi64(_mm256_castsi256_si128(dec64), _mm256_extracti128_si256(dec64, 1));
static_assert((u64)MAXN*255u < (u64(1)<<32),"SAD total fits32 bits");
s64 = _mm_add_epi32(s64, _mm_srli_si128(s64, 8));
const u32 tot = (u32)_mm_cvtsi128_si32(s64);
aacc += tot;
}
{
u32 r = hi - aj;
if (r) {
const u32 base = aj, lk = 0u;
u32 kmask = (1u << r) - 1u;
__m256i ac0 = _mm256_setzero_si256(), ac1 = ac0;
for (int t = 0; t < KK - 1; t++) {
__m256i cv = LDU256((const u8 *)(pbase_sig + (size_t)t * MF_SIGPL + base));
__m256i df = _mm256_subs_epu8(cv, qsv[t]);
if (t & 1) ac1 = _mm256_or_si256(ac1, df);
else ac0 = _mm256_or_si256(ac0, df);
}
__m256i ac = _mm256_or_si256(ac0, ac1);
u32 decm=((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_setzero_si256()))) & kmask;
u32 amb=((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(decm);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
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_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
}
add = aacc;
} else if (c0) {
if (k <= 8) {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + 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_cmpgt_epi32(qTh, h0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h3), all8);
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_cmpgt_epi32(qTh, pvh), all8);
}
add = a0 + a1 + a2 + a3;
}
} else {
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__m256i pv = _mm256_load_si256((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = _mm256_load_si256((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
}
};
if constexpr(KK>0) fringe_dim(std::integral_constant<int,0>{});
if constexpr(KK>1) fringe_dim(std::integral_constant<int,1>{});
if constexpr(KK>2) fringe_dim(std::integral_constant<int,2>{});
if constexpr(KK>3) fringe_dim(std::integral_constant<int,3>{});
if constexpr(KK>4) fringe_dim(std::integral_constant<int,4>{});
if constexpr(KK>5) fringe_dim(std::integral_constant<int,5>{});
if constexpr(KK>6) fringe_dim(std::integral_constant<int,6>{});
if constexpr(KK>7) fringe_dim(std::integral_constant<int,7>{});
if constexpr(KK>8) fringe_dim(std::integral_constant<int,8>{});
if constexpr(KK>9) fringe_dim(std::integral_constant<int,9>{});
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int 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 - 1) * (u64)MF_SIGPL;
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;
g_Q = q;
g_B = (N + q - 1) / q;
}
if (g_Q > 255u) { g_Q = 255u; g_B = (N + g_Q - 1) / g_Q; }
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
#if !MF_BANDFUSE
g_plen_on = 1;
#endif
#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
#if MF_NOFILTB
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
static u8 hist8[MAXN+1] __attribute__((aligned(64)));
for(u32 i=0;i<=N;i++)hist8[i]=0;
for(u32 i=0;i<N;i++){u32 v=xd[i];hist8[v<N?v:0]++;}
    u32 sum=0, pi=0;
    const __m256i last16=_mm256_set1_epi16(0x0f0e);
    for(;pi+16<=N+1;pi+=16){
      __m256i raw=_mm256_cvtepu8_epi16(_mm_loadu_si128((const __m128i*)(hist8+pi)));
      __m256i p=_mm256_add_epi16(raw,_mm256_slli_si256(raw,2));
      p=_mm256_add_epi16(p,_mm256_slli_si256(p,4));
      p=_mm256_add_epi16(p,_mm256_slli_si256(p,8));
      p=_mm256_add_epi16(p,_mm256_shuffle_epi8(_mm256_permute2x128_si256(p,p,0x08),last16));
      __m256i ex=_mm256_sub_epi16(p,raw),base=_mm256_set1_epi32(sum);
      _mm256_storeu_si256((__m256i*)(s_tmp+pi),_mm256_add_epi32(base,_mm256_cvtepu16_epi32(_mm256_castsi256_si128(ex))));
      _mm256_storeu_si256((__m256i*)(s_tmp+pi+8),_mm256_add_epi32(base,_mm256_cvtepu16_epi32(_mm256_extracti128_si256(ex,1))));
      sum+=(u32)_mm_extract_epi16(_mm256_extracti128_si256(p,1),7);
    }
    for(;pi<=N;pi++){u32 c=hist8[pi];s_tmp[pi]=sum;sum+=c;}
bool compact=(sum==N);
if(!compact){
  for(u32 i=0;i<=N;i++)s_cntv[i]=0;
  for(u32 i=0;i<N;i++){u32 v=xd[i];s_cntv[v<N?v:0]++;}
  sum=0;for(u32 i=0;i<=N;i++){u32 c=s_cntv[i];s_cntv[i]=sum;s_tmp[i]=sum;sum+=c;}
}
u32 *gt=s_gt[d],*pp=s_pos[d];
const u64 BMAG=(((u64)1<<41)/B)+1u;
auto scatter=[&](auto tag){
  for(u32 i=0;i<N;i++){
    if constexpr(decltype(tag)::value){
      if(i+32<N){u32 v=xd[i+32];if(v>=N)v=0;__builtin_prefetch(s_tmp+v,0,3);}
      if(i+16<N){u32 v=xd[i+16];if(v>=N)v=0;u32 w=s_tmp[v];__builtin_prefetch(od+(w&0xfffffu)+(w>>20),1,3);}
    }
    u32 v=xd[i];if(v>=N)v=0;
    u32 st,q,num;
    if constexpr(decltype(tag)::value){
      u32 w=s_tmp[v];st=w&0xfffffu;q=st+(w>>20);s_tmp[v]=w+(1u<<20);
      num=(s_tmp[v+1]&0xfffffu)-1u;
    }else{st=s_cntv[v];q=s_tmp[v]++;num=s_cntv[v+1]-1u;}
    od[q]=i;pp[i]=q;gt[i]=st;
#if MF_FB2
    s_posp[(size_t)i*MAXK+d]=q;
#endif
    u32 bb=(u32)(((u64)num*BMAG)>>41);
    s_bdv[(size_t)d*N+i]=bb>Q-1?Q-1:bb;
  }
};
if(compact)scatter(std::true_type{});else scatter(std::false_type{});
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++) {
u8 *dst = s_bd + (size_t)i * MAXK;
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) { dst[d] = s_bdv[(size_t)d * N + i]; p[d] = x[d][i]; }
{
__m256i idx = _mm256_add_epi32(_mm256_cvtepu8_epi32(_mm_loadl_epi64((const __m128i *)dst)),
_mm256_setr_epi32(0, MAXQ, 2*MAXQ, 3*MAXQ,
4*MAXQ, 5*MAXQ, 6*MAXQ, 7*MAXQ));
_mm256_store_si256((__m256i *)(p + 8),
_mm256_i32gather_epi32((const int *)s_bval, idx, 4));
if (k < 8) for (int d = k; d < 8; d++) p[8 + d] = 0x7fffffffu;
}
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
{
u32 M = (u32)(((u64)255 << 23) / (u64)(N - 1)) + 1u;
if (M < 1u) M = 1u;
while (M > 1u && ((u32)(((u64)(N - 1)) * M >> 23)) > 255u) M--;
g_sigmul = M;
}
size_t need = (size_t)k * (k - 1) * (size_t)MF_SIGPL;
static u8 *sigpool = 0; static size_t sigcap = 0;
if (need > ((size_t)220 << 20)) { 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), _mm256_load_si256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
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), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
#if !MF_BANDFUSE
{
size_t mrn = (size_t)MAXK * MAXK * MAXRS;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
#endif
u32 gstart[MAXK+2] = {};u32 fold_maximum[MAXK];for(int d=0;d<MAXK;++d)fold_maximum[d]=1;
for(u32 base=0;base<N;base+=8) {
 const u32 end=base+8<N?base+8:N;
 alignas(32) u32 folds[8];
 if(end-base==8) {
  __m256i v=_mm256_slli_epi32(_mm256_loadu_si256((const __m256i*)(s_gt[0]+base)),4);
  for(int d=1;d<k;++d) {
   __m256i q=_mm256_or_si256(_mm256_slli_epi32(_mm256_loadu_si256((const __m256i*)(s_gt[d]+base)),4),_mm256_set1_epi32(d));
   v=_mm256_min_epu32(v,q);
  }
  _mm256_store_si256((__m256i*)folds,v);
 } else {
  for(u32 i=base;i<end;++i) {
   u32 v=s_gt[0][i]<<4;
   for(int d=1;d<k;++d) {u32 q=(s_gt[d][i]<<4)|(u32)d;if(q<v)v=q;}
   folds[i-base]=v;
  }
 }
#pragma GCC unroll 1
 for(u32 i=base;i<end;++i) {
  const u32 best=folds[i-base]>>4;
  const int ds=(int)(folds[i-base]&15u);
  ++gstart[ds+1];if(best>fold_maximum[ds])fold_maximum[ds]=best;
s_g[i] = (u32)ds;
s_pm[(size_t)i*STRIDE+8+ds]=0x7fffffffu;
if (g_plen_on) {
const u8 *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXRS;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXRS + b]) mr[(size_t)d * MAXRS + b] = bb;
}
}
 }
}
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;}
#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;
if (KSx == (u32)k && (u32)k == (u32)STRIDE) {
u8 *sb = s_sig ? (s_sig + (size_t)d * (k - 1) * MF_SIGPL) : 0;
u32 j = 0;
if (sb && k == 8) {
const __m256i sigmul = _mm256_set1_epi32((int)g_sigmul);
const __m256i one = _mm256_set1_epi32(1);
for (; j + 8 <= N; j += 8) {
#if MF_PFB > 0
if (j + 0 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 0 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row0 = s_pm + (size_t)od[j + 0] * STRIDE;
__m256i rv0 = _mm256_loadu_si256((const __m256i *)row0);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 0) * KSx), rv0);
__m256i r0 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv0, sigmul), 23), one);
#if MF_PFB > 0
if (j + 1 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 1 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row1 = s_pm + (size_t)od[j + 1] * STRIDE;
__m256i rv1 = _mm256_loadu_si256((const __m256i *)row1);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 1) * KSx), rv1);
__m256i r1 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv1, sigmul), 23), one);
#if MF_PFB > 0
if (j + 2 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 2 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row2 = s_pm + (size_t)od[j + 2] * STRIDE;
__m256i rv2 = _mm256_loadu_si256((const __m256i *)row2);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 2) * KSx), rv2);
__m256i r2 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv2, sigmul), 23), one);
#if MF_PFB > 0
if (j + 3 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 3 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row3 = s_pm + (size_t)od[j + 3] * STRIDE;
__m256i rv3 = _mm256_loadu_si256((const __m256i *)row3);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 3) * KSx), rv3);
__m256i r3 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv3, sigmul), 23), one);
#if MF_PFB > 0
if (j + 4 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 4 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row4 = s_pm + (size_t)od[j + 4] * STRIDE;
__m256i rv4 = _mm256_loadu_si256((const __m256i *)row4);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 4) * KSx), rv4);
__m256i r4 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv4, sigmul), 23), one);
#if MF_PFB > 0
if (j + 5 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 5 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row5 = s_pm + (size_t)od[j + 5] * STRIDE;
__m256i rv5 = _mm256_loadu_si256((const __m256i *)row5);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 5) * KSx), rv5);
__m256i r5 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv5, sigmul), 23), one);
#if MF_PFB > 0
if (j + 6 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 6 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row6 = s_pm + (size_t)od[j + 6] * STRIDE;
__m256i rv6 = _mm256_loadu_si256((const __m256i *)row6);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 6) * KSx), rv6);
__m256i r6 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv6, sigmul), 23), one);
#if MF_PFB > 0
if (j + 7 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 7 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row7 = s_pm + (size_t)od[j + 7] * STRIDE;
__m256i rv7 = _mm256_loadu_si256((const __m256i *)row7);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 7) * KSx), rv7);
__m256i r7 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv7, sigmul), 23), one);
__m256i t0 = _mm256_unpacklo_epi32(r0, r1);
__m256i t1 = _mm256_unpackhi_epi32(r0, r1);
__m256i t2 = _mm256_unpacklo_epi32(r2, r3);
__m256i t3 = _mm256_unpackhi_epi32(r2, r3);
__m256i t4 = _mm256_unpacklo_epi32(r4, r5);
__m256i t5 = _mm256_unpackhi_epi32(r4, r5);
__m256i t6 = _mm256_unpacklo_epi32(r6, r7);
__m256i t7 = _mm256_unpackhi_epi32(r6, r7);
__m256i u0 = _mm256_unpacklo_epi64(t0, t2);
__m256i u1 = _mm256_unpackhi_epi64(t0, t2);
__m256i u2 = _mm256_unpacklo_epi64(t1, t3);
__m256i u3 = _mm256_unpackhi_epi64(t1, t3);
__m256i u4 = _mm256_unpacklo_epi64(t4, t6);
__m256i u5 = _mm256_unpackhi_epi64(t4, t6);
__m256i u6 = _mm256_unpacklo_epi64(t5, t7);
__m256i u7 = _mm256_unpackhi_epi64(t5, t7);
if ((int)d != 0) {
__m256i v = _mm256_permute2x128_si256(u0, u4, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((0 < (int)d) ? 0 : 0 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 1) {
__m256i v = _mm256_permute2x128_si256(u1, u5, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((1 < (int)d) ? 1 : 1 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 2) {
__m256i v = _mm256_permute2x128_si256(u2, u6, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((2 < (int)d) ? 2 : 2 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 3) {
__m256i v = _mm256_permute2x128_si256(u3, u7, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((3 < (int)d) ? 3 : 3 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 4) {
__m256i v = _mm256_permute2x128_si256(u0, u4, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((4 < (int)d) ? 4 : 4 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 5) {
__m256i v = _mm256_permute2x128_si256(u1, u5, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((5 < (int)d) ? 5 : 5 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 6) {
__m256i v = _mm256_permute2x128_si256(u2, u6, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((6 < (int)d) ? 6 : 6 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 7) {
__m256i v = _mm256_permute2x128_si256(u3, u7, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((7 < (int)d) ? 7 : 7 - 1) * MF_SIGPL + j,
&bytes, 8);
}
}
}
for (; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
__m256i rv = _mm256_loadu_si256((const __m256i *)row);
_mm256_stream_si256((__m256i *)dp, rv);
if (sb) {
__m256i sg = _mm256_srli_epi32(
_mm256_mullo_epi32(rv, _mm256_set1_epi32((int)g_sigmul)), 23);
sg = _mm256_max_epu32(sg, _mm256_set1_epi32(1));
u32 tmpk[MAXK];
_mm256_storeu_si256((__m256i *)tmpk, sg);
for (int e = 0; e < k; e++) { if (e != (int)d) sb[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + j] = (u8)tmpk[e]; }
}
}
} else
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB4] * 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 *sb2 = s_sig + (size_t)d * (k - 1) * MF_SIGPL;
for (int e = 0; e < k; e++) { if (e == (int)d) continue; sb2[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + j] = sig_nz(row[e]); }
}
}
if (s_sig) g_sig_merged = 1;
}
_mm_sfence();
}
}
}
if (s_sig && !g_sig_merged && k == 8) {
static u8 cached_signature[(size_t)MAXN * 8] __attribute__((aligned(64)));
const __m256i mul=_mm256_set1_epi32((int)g_sigmul),one=_mm256_set1_epi32(1);
const __m256i pick=_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);
u32 si=0;
for(;si+8<=N;si+=8) {
 __m128i c[8];
#pragma GCC unroll 8
 for(int d=0;d<8;++d) {
  __m256i q=_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(x[d]+si)),mul);
  q=_mm256_max_epu32(_mm256_srli_epi32(q,23),one);
  q=_mm256_shuffle_epi8(q,pick);
  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);
 _mm_store_si128((__m128i*)(cached_signature+(size_t)si*8),_mm_unpacklo_epi32(u0,u2));
 _mm_store_si128((__m128i*)(cached_signature+(size_t)si*8+16),_mm_unpackhi_epi32(u0,u2));
 _mm_store_si128((__m128i*)(cached_signature+(size_t)si*8+32),_mm_unpacklo_epi32(u1,u3));
 _mm_store_si128((__m128i*)(cached_signature+(size_t)si*8+48),_mm_unpackhi_epi32(u1,u3));
}
for(;si<N;++si)for(int d=0;d<8;++d)cached_signature[(size_t)si*8+d]=sig_nz(x[d][si]);
for (int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*(k-1)*MF_SIGPL;
u32 j=0;
for (;j+8<=N;j+=8) {
__m128i r0=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+0]*8));
__m128i r1=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+1]*8));
__m128i r2=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+2]*8));
__m128i r3=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+3]*8));
__m128i r4=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+4]*8));
__m128i r5=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+5]*8));
__m128i r6=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+6]*8));
__m128i r7=_mm_loadl_epi64((const __m128i *)(cached_signature+(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 < d ? 0 : 0-1)*MF_SIGPL+j),v0);
if(d!=1) _mm_storel_epi64((__m128i *)(base+(size_t)(1 < d ? 1 : 1-1)*MF_SIGPL+j),_mm_srli_si128(v0,8));
if(d!=2) _mm_storel_epi64((__m128i *)(base+(size_t)(2 < d ? 2 : 2-1)*MF_SIGPL+j),v1);
if(d!=3) _mm_storel_epi64((__m128i *)(base+(size_t)(3 < d ? 3 : 3-1)*MF_SIGPL+j),_mm_srli_si128(v1,8));
if(d!=4) _mm_storel_epi64((__m128i *)(base+(size_t)(4 < d ? 4 : 4-1)*MF_SIGPL+j),v2);
if(d!=5) _mm_storel_epi64((__m128i *)(base+(size_t)(5 < d ? 5 : 5-1)*MF_SIGPL+j),_mm_srli_si128(v2,8));
if(d!=6) _mm_storel_epi64((__m128i *)(base+(size_t)(6 < d ? 6 : 6-1)*MF_SIGPL+j),v3);
if(d!=7) _mm_storel_epi64((__m128i *)(base+(size_t)(7 < d ? 7 : 7-1)*MF_SIGPL+j),_mm_srli_si128(v3,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*8;
for(int e=0;e<k;e++)if(e!=d)base[(size_t)(e<d?e:e-1)*MF_SIGPL+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_SIGPL;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) { if (e == (int)d) continue; base[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + j] = sig_nz(row[e]); }
}
}
}
#if !MF_BANDFUSE
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) * MAXRS;
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;
{size_t bytes=(size_t)N*sizeof(u32)+4096;if(bytes>inc_cap){inc_anchor=(u32*)pool_alloc(bytes);inc_cap=inc_anchor?bytes:0;}inc_on=inc_anchor&&(u64)N<=(u64)AS_WORDS*64;}
#endif
for (int ds = 0; ds < k; ds++) {
u32 q0 = gstart[ds], q1 = gstart[ds + 1];
if (q0 >= q1) continue;
#if MF_WTRIM
{
const u32 mx=fold_maximum[ds];
u32 wp = ((mx >> 6) + 8u) & ~7u;
Wpool = (wp < W) ? wp : W;
}
#endif
#if !MF_BANDFUSE
#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) * MAXRS;
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
const u32 *const spp = s_pos[ds];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
u64 *base = g_plen_on ? s_bs : (s_bs + (size_t)e * (Q + 1) * Wpool);
u64 *cur = s_acc;
for (u32 w = 0; w < Wpool; w++) cur[w] = 0;
if (g_plen_on) {
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 = FPOS(d,b);
while (ptr + 12 <= pend) {
if (ptr+32 < pend) _mm_prefetch((const char *)(spp + od[ptr+32]), _MM_HINT_T0);
u32 q = spp[od[ptr+0]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+1]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+2]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+3]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+4]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+5]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+6]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+7]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+8]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+9]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+10]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+11]]; cur[q >> 6] |= 1ull << (q & 63);
ptr += 12; }
while (ptr < pend) { u32 q = spp[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
#if MF_NT
{
u32 w = 0;
for (; w + 4 <= wlen; w += 4)
_mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
for (u32 w = 0; w < wlen; w++) dst[w] = cur[w];
#endif
}
_mm_sfence();
#endif
}
#endif
query_group<8>(N, out, ds, q0, q1, Q, Wpool);
}
dimension_major_fringe<8>(N,out,B);
}
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
#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

#pragma GCC optimize("peel-loops")
#pragma GCC optimize("no-jump-tables")

评测结果