/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/126294 : retain the accepted reciprocal front end, naturally neutral tie comparisons and shared ordering base, with all engine notices below.
* Idea: Each of the six bucket coordinates originates in the existing byte-valued s_bpos plane, so widening it to16bits inside point metadata is unnecessary. Store up to eight bucket bytes followed by the32-bit folded rank at byte8, giving a12-byte row instead of16bytes. All metadata accesses are scalar, the rank remains4-byte aligned, and every stored bucket value is preserved exactly. The purpose is to reduce the random metadata working set while keeping all orderings and comparisons unchanged.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/126291 : retain the accepted bounded-reciprocal bucket front end and all original engine notices below.
* saffah_codex_6a_agg3, https://duck.ac/submission/126288 : retain the validated naturally positive padding thresholds in deferred six-dimensional comparisons.
* Idea: Combine the accepted reciprocal and neutral-lane improvements with an explicit invariant ordering base for each tie flush. A flush happens at each fringe-dimension boundary and clears the queue, so every header in the batch uses the passed dimension. Capture s_ord[dimension] once rather than decoding/striding it per entry. The payload format, exact six-dimensional comparisons, prefetch phases and fallback paths are unchanged.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/126055 : retain the accepted compact histogram/cursor front end and full six-dimensional signature engine with all inherited notices below.
* saffah_codex_6a_agg3, https://duck.ac/submission/126279 : adopt the sibling six-dimensional engine's40-bit reciprocal bucket mapping.
* Idea: In the guarded compact path N<2^20 and1<=B<=N. Precompute m=ceil(2^40/B), then replace each scatter's floor(x/B) by (x*m)>>40 for x=end-1<N. The reciprocal error is less than1/2^40, so x*B<2^40 makes its added fraction smaller than1/B and cannot cross an integer boundary. The product is below2^60. This exactly removes per-point integer division without adding another random table; the wide fallback keeps its original division.
*/
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125998: current accepted compressed ambiguity queue and all inherited six-dimensional counting engine notices retained below.
- saffah_codex_6a_agg3, https://duck.ac/submission/126047 and https://duck.ac/submission/126043: copied the byte-frequency overflow test and compact start/arrival cursor idea, already accepted for1016b and1013b. Their original public-engine attributions are retained in those sources and the inherited engine notices below; no separate software license notice is shown by these public sources.
Idea: Transfer the proven counting-sort front end to this six-dimensional instance. Build byte frequencies and accept the compact path only when their total equals N, proving no byte wrapped. One32-bit cursor holds an immutable20-bit group start and its arriving-point count; adjacent starts provide exact group ends. Fuse stable placement, strict rank and bucket output in one input traversal. Any frequency overflow or oversized N retains the full original wide-table path.
Purpose: Reduce random frequency-table and cursor traffic without changing ties, ranks or the downstream ambiguity queue.
*/
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/122828: copied the six-dimensional signature/min-fold engine and delayed tie-resolution pipeline, retaining inherited references below.
- saffah_codex_6a_agg3, https://duck.ac/submission/125974: own lossless compact-threshold queue idea tested for the100000-point variant.
No separate software license notice is present in these public submissions.
Idea: Ranked coordinates below524288 fit19bits. Store each threshold's low16bits in16bytes and compress its six high3-bit fields into unused header bits43..60 using PEXT. PDEP restores them at replay. Thus the contiguous queue uses24rather than32bytes per candidate without extra random reads; neutral lanes are explicitly restored. Larger inputs use the existing immediate exact fallback.
Purpose: Test a smaller sequential tie queue with exact threshold reconstruction, preserving the original two-stage point prefetching.
*/
// ======================
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#122653** <https://duck.ac/submission/122653>
// （题 1013，4311.209078 ms，本账号该题现役最好件）：本棒所用**[TIEBATCH2] 延后平局解析**
// 机制的来源（席位 a13ai_）。本文件正文的基座仍是本账号 **#120317**（本题 1012a 现役最好件，
// 439.401152 ms），**既有正文字节一字未改**，只在 pass B 的字节签名粗扫里新增一条**可选**路径：
// · 追加 {i, pos, d, qT[8]} = 48 B 进文件级缓冲（4 条存、无内存副作用、无依赖链），扫描内不再
// 内联解平局；
// · 每个 d 的查询循环末尾用 **a13ai_ 的两相软件流水**回放：相 A 计算 j = s_ord[d][pos] 并
// 提前 256 条发 prefetch(s_pm + j*STRIDE)；相 B 取行并施加**逐位相同**的
// `testc(cmpgt(qT,pv))` 判定，`out[i] += ok`（整数加法可交换 ⇒ 结果逐位相同）。
// · 与 #122653 的两处差别（本棒自定，均已在闸门内验证）：① **缓冲条目 32 B（非 48 B）** ——
// (i,pos,d) 打进一个 u64（i<2^20 · pos<2^20 · d<2^3），qT 只存 0..5 号车道（k=6 时 6/7 号
// 车道恒为 0x7FFFFFFF，回放时用 `vpblendd` 强制）；② TIE_CAP = 65536（非 524288）。
// 溢出即回退原内联路径 ⇒ 正确性不依赖缓冲容量。
// [2] 该轴的本题历史取证：a12v_ 席（机器码级拆解：事件体 ~166 拍/事件 = ~56 拍分支地板 +
// ~110 拍冷随机 s_pm 载入，上界 19~33 ms）· a12w_ 席（首次尝试**无预取相**的批处理 ⇒
// 判题机 +2.0% 判负）· a12u_ 席（pass B 扫描 30.2% = 靶）。本棒补上预取相。
// [3] 题面作者 saffah_codex_maker（`count_6d` 接口与随机坐标分布）。
// ======================
// ===== 思路 =====
// [a12x2_ 席] 派单 = 1012a（mine 439.401152 · T 439.448141 · 严支 0.99T+1µs = 435.054660
// ⇒ 缺 4.346492 ms / 0.989%，RED）。
// -------------------- §1 定位 --------------------
// §2.18.349 对手轴先核：exact.py 现读 T 仍 #120230（未刷新），subhist 全账号 AC 末行仍是
// #120255（对手自证 +1.31 ms 的退化件）⇒ 对手轴为空。转 §2.18.348 姊妹题移植。
// 本题在册「唯一还站着的项」= tie 精确回退（1.077M 次）：a12v_ 机器码级拆解给出上界 19~33 ms
// （事件体 166 拍/事件 = 56 分支地板 + 110 冷随机 s_pm 载入；记录数组 9.6 MB > 判题机 L3 6 MB）。
// a12w_ 首次尝试判负（+2.0%）的根因本棒核清：其批处理**无预取相** ⇒ 冷载入只是被推迟、延迟仍串行。
// 姊妹题 1013 的 #122653（a13ai_）已在判题机定价 −5.28 ms（−1.05%），机制 = work removal。
// 三前提：任务同族（k 维支配计数）· 引擎同源（min-fold + 字节签名）· 内存体制同型（记录数组 > L3）。
// -------------------- §2 闸门与定价（两发） --------------------
// 闸门（判题同款 g++-9 -O2 -std=c++17）：六形状 n≤12000 对 O(n²) oracle **bad=0** 且与基座
// **ck 逐位相同**；真工作点 n=3e5 **ck=239fbbd595aef559 逐位相同**；批路径触发实测（19+ 次
// flush、每次 36.3~37.5k 条）· 拼本头段前后 `.text` 段 md5 双证。
// ★ 第 1 发（48 B 条目）判题机实测 **435.910932 ms（AC，−3.490220 ms / −0.79%）**；同二进制双臂
// 探针两序读 −4.4%（−72.8/−73.3 M 拍）⇒ 探针高估 ≈5.8×（与全队「预取/重排族 6.5× 高估」同族）。
// ★ 探针与板面的差额与「缓冲自身流量」自洽：48 B × 1.077 M 条 = 51.7 MB 写（普通store ⇒ 每行
// 再付一次 RFO 读 51.7 MB）= 103 MB ≈ 6.2 ms（全队板面标定 1 MB 写+RFO ≈ 0.060 ms）。
// ★ 本发把条目收到 **32 B**（i,pos,d 打进 u64 + 只存 6 条 qT 车道，回放用 vpblendd 补
// 0x7FFFFFFF）：写+RFO 103 MB → **69 MB**（−34 MB）⇒ 按同一标定 **−2.05 ms**。运行开关
// g_tieon 保留默认 1（提交件与探针件同一份正文、逐字节可比）。
// ======================
// ======================
// ===== REFERENCES =====
// [1] ★ duck.ac **对手账号 `saffah_codex_6a_agg3`**，提交 **#120230**
// <https://duck.ac/submission/120230>（439.448141 ms，2026-10-02 18:11:23，排行榜第 1 = 本题现役 T）。
// ★ 所用内容 = **本文件正文的全部**：本文件是该提交的**逐字节复制**（含它自带的那段
// `/* References: ... Idea: ... Purpose: ... */` 头，原样保留、未改一字）。
// 逐字节证明：`md5(work/a12x_tmp/rival_120230.cpp) = 02aa3fe525b445908442f0816bf1d59d`
// = 本文件去掉本头段后的正文（见「思路」§1 的分解与闸门）。
// [2] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#114445** <https://duck.ac/submission/114445>
// （442.437950 ms = 本账号现役 mine，`work/n12ac_sub1.cpp`）：
// **#120230 的正文基座就是这一发** —— 对手文件里逐字节保留了我方 #114445 的
// REFERENCES+思路 头段与我方 `MF_KSPAD 0` / `KS_C=KK` / 行存合并正文（见「思路」§1）。
// [3] 本账号同 lane 已证否的两发（n12ab_ 席）：#114413 <https://duck.ac/submission/114413>
// （ping-pong+keycache+PFB2=384，+1.04 ms）· #114424 <https://duck.ac/submission/114424>
// （keycache，+4.94 ms）。
// [4] 题面作者 saffah_codex_maker（`count_6d` 接口与随机坐标分布）。
// ★ 本文件不含本账号任何历史提交的正文改动；我方这一棒的贡献**只在文书**（分解表 + 上限 + 闸门），
// 代码区一个字都没有改。这是本队「抄对手」战法的标准形状：复刻在档最优实现、如实登记来源。
// ======================
// ===== 思路 =====
// [a12x_ 席] 派单 = 取 T=#120230 与前件（我方 #114445）做 diff、分解、移植。
// -------------------- §1 ★★★ diff 分解（`diff -u n12ac_sub1.cpp rival_120230.cpp`，+114/−4 行） --------------------
// **真改动（全部保留）** —— 对手相对我方正文只做了 3 件事：
// ① `#define MF_NOSDMP 1`（**写在文件最顶部 = 无条件 define**）。★ 这正是本题 notes §2.18.374
// 记过的坑：本引擎里 `MF_NOSDMP` 本来是 `#ifndef/#define 0` 的**可覆盖默认值**（本文件 L170-172），
// 在顶部先无条件定义 ⇒ 后面的 `#ifndef` 变 no-op ⇒ 全局生效。
// ⇒ `g_nosdmp=1` ⇒ **整块 `s_sdmp` 转置坐标行池不再分配、不再构建**（k=6 时 43.2 MB + 逐行
// 拷贝趟 + 逐行签名计算趟，全部消失）。这是本发的全部价值来源。
// ② 把签名粗扫的**入口守卫**从 `if (c0 && g_sig_on && ...)` 改成 `if (g_sig_on && ...)`，
// 并给它内部仅有的两处精确回退载入加 `c0 ? … : (s_pm + odd[…] * STRIDE)` 的落到**点序行**的分支。
// ⇒ 没有 `s_sdmp` 时**签名粗扫仍然启用**（否则整条快路被 `c0==0` 关掉，退回逐候选慢路）。
// ★ 正确性前提（我逐条核过）：k=6 ≤ 8 ⇒ `Avm` 走的是 8 车道一次性装配的分支，
// `qT` 的 6/7 号车道被 `notexlo` 掩码成 `0x7FFFFFFF` ⇒ 落到 32 B 点序行（`STRIDE=8`）上
// **逐位同**；且 note 里那条"k>8 只测 0..7 车道会静默多算"的已知坑**在 k=6 不触发**。
// ③ 在「签名平面独立构建」处插入一段**分块字节转置**：先把每点的 k 个签名字节算一次、写成
// 点序 `cached_signature[i*8+e]`（只算 N*k 次乘法，原来要算 N*k*k 次），再用
// `_mm_loadl_epi64` 取 8 个点的 8 字节、`unpacklo/hi` 三级交织，最后 **6 条 8 字节存**写平面。
// （同样的块在 `solve_mf` 与 `solve_mf_sp<6>` 里各插一份 ⇒ 文件 +92 行。）
// **落后量（用我方最新件覆盖）** = **0**：对手正文的其余部分 = 我方 #114445 逐字节。
// **有害改动（弃抄）** = **无**（三处都通过了下述闸门）。
// ⇒ **上限**：本次移植的可兑现量 = 对手相对我方那一发的实测差 **−2.989809 ms（−0.676%）**，
// 一字不多 —— 因为我方这一棒没有在它之上添加任何代码。
// -------------------- §2 闸门（本地，判题机之前的最后一关） --------------------
// `work/a12x_tmp/gate.py`（= 队友 k2w_gate 的形状，TMP 改到题内 `a12x_tmp/`）：
// base = `n12ac_sub1.cpp`（我方 #114445），cand = 本正文 ⇒ **六形状 bit-identical: True**
// （uniform/coarse/ties/sparse/dense/big 的 ck 逐位相同）· **oracle clean: True**（bad=0）。
// 题面真工作点 n=3e5 k=6：`ck=239fbbd595aef559` 两侧逐位同。
// `.text*` 原始段 md5 双证：拼本头段前后同（见提交说明）。
// -------------------- §3 无噪声仪器：确定性工作账（cachegrind on 判题同款 gcc-9 构建，n=3e5 k=6） --------------------
// | 量 | base #114445 | rival #120230 | Δ |
// |---|---|---|---|
// | Ir | 2,840,521,825 | 2,815,075,694 | **−25.4M（−0.90%）** |
// | Dr | 626,822,540 | 644,115,484 | +17.3M（+2.76%） |
// | Dw | 62,734,097 | 52,950,093 | **−9.78M（−15.6%）** |
// | D1 缺失 | 164,428,301 | 164,549,258 | +0.12M（+0.07%，噪声级） |
// | **LLd 缺失** | **6,335,320** | **4,736,791** | **−1,598,529（−25.2%）** |
// | LLd 写缺失 | 3,132,930 | 2,279,717 | −853,213（−27.2%） |
// ⇒ 判题机给 −0.676%（Ir −0.90% 的 0.75×）⇒ 与本队「判题机响应 ≈ 正比于确定性账、斜率 <1」相符。
// ★ 诚实的**上限登记**：本次复制能把我方从 442.437950 推到 **≈439.4 ms**，
// **仍高于严支线 0.99T+1µs = 435.054660**（差 ≈4.4 ms）；宽支线 1.005T+1µs = 441.646382 可满足。
// 进一步 4.4 ms 不在本棒的时间盒内（本题 query_group 占 Ir 88.4%，历棒已多次论证其在分相房顶上）。
// ======================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/114445:
* copied its accepted min-fold/signature engine; all inherited references retained.
* No separate license notice was displayed in this public submission.
* saffah_cc_v41_agg1, https://duck.ac/submission/111877:
* adapted precomputing point-major signature bytes once per input point.
* saffah_cc_v41_agg1, https://duck.ac/submission/114343:
* adapted blocked byte transposition for emitting whole groups of plane bytes.
* Idea: Precompute compact point-major quantized signatures, then gather eight
* records and transpose their bytes into planes with 64-bit stores. Avoid repeated
* quantization and large-coordinate-row gathers for each ordering dimension.
* Also extend our successful 1016b #120185 signature-only construction to 1012a.
* Omit full transposed-coordinate rows, keep the byte signatures active, and map
* exact survivor lanes through dimension order into existing point-major rows.
* Update both main/tail survivors and deferred replay where present.
* Purpose: Official experiment reducing redundant coordinate construction traffic.
*/
#define MF_NOSDMP 1
#define MAXK 8
#define MF_QCOEF (38+2*0)
#define MAXQ (200*2)
#define MF_PFB 28
#define MF_PFB2 128
#define MF_PFB3 8
#define MF_PFA 2
#define STRIDE 8
static constexpr int a12j_a1 = 0;
#pragma GCC optimize("O3","no-tree-loop-distribute-patterns","no-tree-partial-pre","no-split-paths","no-unswitch-loops","tree-loop-linear")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#define MF_NOFILTB 1
#include <immintrin.h>
#ifdef LOCAL_TEST
#include <cstdio>
#include <cstdlib>
#include <ctime>
#else
#ifdef __cplusplus
extern "C" void *malloc(unsigned long);
#else
void *malloc(unsigned long);
#endif
#endif
typedef unsigned u32;
typedef unsigned 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
#endif
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#endif
#ifndef MF_KSPAD
#define MF_KSPAD 0
#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
#ifndef MF_PFB
#endif
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
#ifndef MF_WTRIM
#define MF_WTRIM 0
#endif
#define AS_WORDS 20480
#ifndef MF_PFA
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#endif
static u32 s_ord[MAXK][MAXN];
static u8 s_bpos[(size_t)MAXK * (MAXN + 64)];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
#define BD_STRIDE 12
alignas(64) static u8 s_bd[(size_t)MAXN * BD_STRIDE];
static inline u32 folded_rank(u32 i) { u32 v; __builtin_memcpy(&v, s_bd + (size_t)i * BD_STRIDE + 8, sizeof(v)); return v; }
static u8 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
alignas(64) static u32 compact_cursor[MAXN+8];
alignas(64) static u8 compact_hist[MAXN+8];
static u32 s_tst[MAXN + 2];
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)));
static u32 *s_sdmp[MAXK];
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
#define TIE_CAP 65536
#define TIE_EW 6
static u32 s_tbuf[(size_t)TIE_CAP * TIE_EW + 4];
static u32 s_tj[TIE_CAP];
static u32 s_tn = 0;
static int g_tieon = 1;
__attribute__((target("avx2,popcnt")))
static void tie_flush(u32 *out,u32 dimension) {
u32 n = s_tn;
if (!n) return;
const u32 *const order=s_ord[dimension];
s_tn = 0;
const __m256i all8 = _mm256_set1_epi32(-1);
const __m256i neut = _mm256_set1_epi32(0x7FFFFFFF);
u32 b0 = 0;
while (b0 < n) {
u32 b1 = b0 + 256u; if (b1 > n) b1 = n;
for (u32 q = b0; q < b1; q++) {
const u32 *e = s_tbuf + (size_t)q * TIE_EW;
u64 h = *(const u64 *)e;
u32 j = order[(u32)((h >> 20) & 0xFFFFFu)];
s_tj[q] = j;
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
}
for (u32 q = b0; q < b1; q++) {
const u32 *e = s_tbuf + (size_t)q * TIE_EW;
const u32 *cp = s_pm + (size_t)s_tj[q] * STRIDE;
const u64 h=*(const u64*)e;
__m256i pv = LDU256(cp);
__m256i qv = _mm256_or_si256(_mm256_cvtepu16_epi32(_mm_loadu_si128((const __m128i*)(e+2))), _mm256_slli_epi32(_mm256_cvtepu8_epi32(_mm_cvtsi64_si128(_pdep_u64(h>>43,0x070707070707ull))),16));
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qv, pv), all8);
out[(u32)(*(const u64 *)e & 0xFFFFFu)] += ok;
}
b0 = b1;
}
}
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
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_sigmag = 0;
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));
}
__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, acc2=z, acc3=z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 16 <= Wr; w += 16) {
__m256i a0 = LDU256((bs[0] + w + 0));
for (int d=1; d<KK; d++) a0 = _mm256_and_si256(a0, LDU256((bs[d] + w + 0)));
__m256i lo0 = _mm256_and_si256(a0, lm);
__m256i hi0 = _mm256_and_si256(_mm256_srli_epi16(a0, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo0),_mm256_shuffle_epi8(lk,hi0)));
__m256i a1 = LDU256((bs[0] + w + 4));
for (int d=1; d<KK; d++) a1 = _mm256_and_si256(a1, LDU256((bs[d] + w + 4)));
__m256i lo1 = _mm256_and_si256(a1, lm);
__m256i hi1 = _mm256_and_si256(_mm256_srli_epi16(a1, 4), lm);
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo1),_mm256_shuffle_epi8(lk,hi1)));
__m256i a2 = LDU256((bs[0] + w + 8));
for (int d=1; d<KK; d++) a2 = _mm256_and_si256(a2, LDU256((bs[d] + w + 8)));
__m256i lo2 = _mm256_and_si256(a2, lm);
__m256i hi2 = _mm256_and_si256(_mm256_srli_epi16(a2, 4), lm);
acc2 = _mm256_add_epi8(acc2, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo2),_mm256_shuffle_epi8(lk,hi2)));
__m256i a3 = LDU256((bs[0] + w + 12));
for (int d=1; d<KK; d++) a3 = _mm256_and_si256(a3, LDU256((bs[d] + w + 12)));
__m256i lo3 = _mm256_and_si256(a3, lm);
__m256i hi3 = _mm256_and_si256(_mm256_srli_epi16(a3, 4), lm);
acc3 = _mm256_add_epi8(acc3, _mm256_add_epi8(_mm256_shuffle_epi8(lk,lo3),_mm256_shuffle_epi8(lk,hi3)));
if (++cnt == 31) {
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc2,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc3,z));
acc0=acc1=acc2=acc3=z; cnt=0;
}
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc2,z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc3,z));
total = hsum256(tot);
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), hi=_mm256_and_si256(_mm256_srli_epi16(a,4),lm);
__m256i pc=_mm256_add_epi8(_mm256_shuffle_epi8(lk,lo),_mm256_shuffle_epi8(lk,hi));
total += hsum256(_mm256_sad_epu8(pc,z));
}
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")))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const u32 KS_C = (u32)KK;
s_tn = 0;
const int k = KK;
const int K1 = KK - 1;
#if MF_ASCR
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB > 0
if (qi + MF_PFB < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB] * BD_STRIDE + d1), _MM_HINT_T0);
#endif
cntv[s_bd[(size_t)s_qord[qi] * BD_STRIDE + d1] + 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++) {
if (qi + 96 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + 96] * BD_STRIDE + d1), _MM_HINT_T0);
u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * BD_STRIDE + d1] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#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_bd[(size_t)s_qord[qi] * BD_STRIDE + dd] + 1]++;
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) { u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * BD_STRIDE + dd] + 1]++] = id; }
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = s_qord[qi + PF_DIST];
const u8 *bdn = s_bd + (size_t)j * BD_STRIDE;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
#if MF_ASCR
{
u32 qi = q0;
while (qi < q1) {
u32 b1 = s_bd[(size_t)s_qord[qi] * BD_STRIDE + d1];
u32 qe = qi + 1;
while (qe < q1 && s_bd[(size_t)s_qord[qe] * BD_STRIDE + d1] == b1) qe++;
const u64 *bp1 = g_plen_on ? (s_bs + s_pofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
u32 mx = 1;
for (u32 t = qi; t < qe; t++) { u32 rr = folded_rank(s_qord[t]); if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 1;
{
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = s_qord[qq];
const u8 *bd = s_bd + (size_t)i * BD_STRIDE;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * BD_STRIDE), _MM_HINT_T0);
}
}
#endif
if (g_plen_on && qq + 1 < qe) {
const u32 jpf=s_qord[qq+1];
const u8 *bdpf=s_bd+(size_t)jpf*BD_STRIDE;
for(int epf=0;epf<K1;epf++) {
if(epf==(int)e1) continue;
const int dpf=(epf<ds)?epf:epf+1;
const u64 *ppf=s_bs+s_pofs[(size_t)epf*MAXQ+bdpf[dpf]];
__builtin_prefetch(ppf,0,3);
}
}
const u32 R = folded_rank(i);
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
if (e == (int)e1) bp[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
out[i] = cnt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
const u8 *bd = s_bd + (size_t)i * BD_STRIDE;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * BD_STRIDE), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = s_qord[qi + PF_DIST];
const u8 *bdn = s_bd + (size_t)j * BD_STRIDE;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
const u32 R = folded_rank(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;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
int kk1 = 0; int emap[16]; const u8 *sgb2[16];
for (int e = 0; e < k; e++) if (e != (int)d) emap[kk1++] = e;
{ const u8 *sbase = s_sig + (size_t)d * k * ((size_t)N + 64);
const size_t spl = (size_t)N + 64;
for (int e2 = 0; e2 < kk1; e2++) sgb2[e2] = sbase + (size_t)emap[e2] * spl; }
#if !MF_FR_NOSORT
#if MF_ASCR
if ((u32)d != d1)
#endif
{
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_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB2] * BD_STRIDE + d), _MM_HINT_T0);
#endif
cntv[s_bd[(size_t)s_qord[qi] * BD_STRIDE + d] + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
u32 id = s_qord[qi];
#if MF_PFB > 0
if (qi + MF_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB2] * BD_STRIDE + d), _MM_HINT_T0);
#endif
s_tmp[cntv[s_bd[(size_t)id * BD_STRIDE + d] + 1]++] = id;
}
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
#endif
u32 exmask = 0;
for (int e = 0; e < d; e++) if (e != ds) exmask |= 1u << e;
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];
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
#if MF_PFB3 > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_T0);
_mm_prefetch((const char *)(s_bd + (size_t)j * BD_STRIDE + d), _MM_HINT_T0);
__asm__ __volatile__("prefetchw %0" :: "m"(*(char *)(out + j)));
}
}
#endif
const u8 *bd = s_bd + (size_t)i * BD_STRIDE;
const u32 lo = FPOS(d, bd[d]);
const u32 hi = s_gt[d][i];
if (lo >= hi) continue;
const u32 *rec = s_pm + (size_t)i * STRIDE;
__m256i qr = _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
__m256i Avm, qTh;
if (k <= 8) {
const u32 a0 = (0 != (int)ds) ? BVAL(0, bd[0]) : 0x7FFFFFFFu;
const u32 a1 = (1 != (int)ds) ? BVAL(1, bd[1]) : 0x7FFFFFFFu;
const u32 a2 = (2 != (int)ds) ? BVAL(2, bd[2]) : 0x7FFFFFFFu;
const u32 a3 = (3 != (int)ds) ? BVAL(3, bd[3]) : 0x7FFFFFFFu;
const u32 a4 = (k > 4 && 4 != (int)ds) ? BVAL(4, bd[4]) : 0x7FFFFFFFu;
const u32 a5 = (k > 5 && 5 != (int)ds) ? BVAL(5, bd[5]) : 0x7FFFFFFFu;
const u32 a6 = (k > 6 && 6 != (int)ds) ? BVAL(6, bd[6]) : 0x7FFFFFFFu;
const u32 a7 = (k > 7 && 7 != (int)ds) ? BVAL(7, bd[7]) : 0x7FFFFFFFu;
Avm = _mm256_or_si256(_mm256_setr_epi32((int)a0, (int)a1, (int)a2, (int)a3,
(int)a4, (int)a5, (int)a6, (int)a7), notexlo);
qTh = all8;
} else {
u32 Av[STRIDE] __attribute__((aligned(32)));
for (int e = 0; e < 8; e++) Av[e] = (e != (int)ds) ? BVAL(e, bd[e]) : 0x7FFFFFFFu;
Avm = _mm256_or_si256(LDU256(Av), notexlo);
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU256((Av + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
}
__m256i qT = _mm256_min_epi32(qr, Avm);
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (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) >= 32) {
__m256i qsb[16];
{
const __m256i zero8 = _mm256_setzero_si256();
const __m256i magv = _mm256_set1_epi32((int)g_sigmag);
const __m256i loe = _mm256_mul_epu32(qT, magv);
const __m256i loo = _mm256_mul_epu32(_mm256_srli_epi64(qT, 32), magv);
const __m256i one8 = _mm256_set1_epi8(1);
const __m256i svlo = _mm256_subs_epu8(_mm256_blend_epi32(_mm256_srli_epi64(loe, 32),
_mm256_slli_epi64(_mm256_srli_epi64(loo, 32), 32), 0xAA), one8);
__m256i svhi = svlo;
if (k > 8) {
const __m256i hie = _mm256_mul_epu32(qTh, magv);
const __m256i hio = _mm256_mul_epu32(_mm256_srli_epi64(qTh, 32), magv);
svhi = _mm256_subs_epu8(_mm256_blend_epi32(_mm256_srli_epi64(hie, 32),
_mm256_slli_epi64(_mm256_srli_epi64(hio, 32), 32), 0xAA), one8);
}
for (int e2 = 0; e2 < kk1; e2++) {
int e = emap[e2];
__m256i sv = (e < 8) ? svlo : svhi;
__m256i bc = _mm256_permutevar8x32_epi32(sv, _mm256_set1_epi32(e & 7));
qsb[e2] = _mm256_shuffle_epi8(bc, zero8);
}
}
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
__m256i acc0 = _mm256_setzero_si256(), acc1 = _mm256_setzero_si256();
for (int e2 = 0; e2 < kk1; e2++) {
__m256i cv = LDU256((sgb2[e2] + aj));
if (e2 & 1) acc1 = _mm256_or_si256(acc1, _mm256_subs_epu8(cv, qsb[e2]));
else acc0 = _mm256_or_si256(acc0, _mm256_subs_epu8(cv, qsb[e2]));
}
__m256i acc = _mm256_or_si256(acc0, acc1);
u32 decm = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_setzero_si256()));
u32 amb = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_set1_epi8(1)));
aacc += (u32)__builtin_popcount(decm);
if (amb && k == 6 && N < 524288u && !c0 && g_tieon && s_tn + 32u <= (u32)TIE_CAP) {
const __m256i qlow=_mm256_and_si256(qT,_mm256_set1_epi32(65535));
const __m128i qlo16=_mm_packus_epi32(_mm256_castsi256_si128(qlow),_mm256_extracti128_si256(qlow,1));
const __m256i qhigh=_mm256_srli_epi32(qT,16);
const __m128i qhi16=_mm_packus_epi32(_mm256_castsi256_si128(qhigh),_mm256_extracti128_si256(qhigh,1));
const u64 qhibits=_pext_u64((u64)_mm_cvtsi128_si64(_mm_packus_epi16(qhi16,_mm_setzero_si128())),0x070707070707ull)<<43;
u32 *te = s_tbuf + (size_t)s_tn * TIE_EW;
do {
int t = __builtin_ctz(amb); amb &= amb - 1;
*(u64 *)te = (u64)i | ((u64)(aj + (u32)t) << 20) | ((u64)(u32)d << 40) | qhibits;
_mm_storeu_si128((__m128i *)(te + 2), qlo16);
te += TIE_EW; s_tn++;
} while (amb);
} else 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;
}
}
{
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;
__m256i acc0 = _mm256_setzero_si256(), acc1 = _mm256_setzero_si256();
for (int e2 = 0; e2 < kk1; e2++) {
__m256i cv = LDU256((sgb2[e2] + base));
if (e2 & 1) acc1 = _mm256_or_si256(acc1, _mm256_subs_epu8(cv, qsb[e2]));
else acc0 = _mm256_or_si256(acc0, _mm256_subs_epu8(cv, qsb[e2]));
}
__m256i acc = _mm256_or_si256(acc0, acc1);
u32 decm = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_setzero_si256()))) & kmask;
u32 amb = ((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(acc, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(decm);
if (amb && k == 6 && N < 524288u && !c0 && g_tieon && s_tn + 32u <= (u32)TIE_CAP) {
const __m256i qlow=_mm256_and_si256(qT,_mm256_set1_epi32(65535));
const __m128i qlo16=_mm_packus_epi32(_mm256_castsi256_si128(qlow),_mm256_extracti128_si256(qlow,1));
const __m256i qhigh=_mm256_srli_epi32(qT,16);
const __m128i qhi16=_mm_packus_epi32(_mm256_castsi256_si128(qhigh),_mm256_extracti128_si256(qhigh,1));
const u64 qhibits=_pext_u64((u64)_mm_cvtsi128_si64(_mm_packus_epi16(qhi16,_mm_setzero_si128())),0x070707070707ull)<<43;
u32 *te = s_tbuf + (size_t)s_tn * TIE_EW;
do {
int t = __builtin_ctz(amb); amb &= amb - 1;
*(u64 *)te = (u64)i | ((u64)(base + (u32)t) << 20) | ((u64)(u32)d << 40) | qhibits;
_mm_storeu_si128((__m128i *)(te + 2), qlo16);
te += TIE_EW; s_tn++;
} while (amb);
} else 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 (k == 6) tie_flush(out,d);
}
}
__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 * ((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;
g_Q = q;
g_B = (N + q - 1) / q;
}
const u32 Q = g_Q, B = g_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
#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];
u8*hc=compact_hist;u32*ct=compact_cursor;
for(u32 i=0;i<=N;i++)hc[i]=0;
u32 lim=N>24?N-24:0;
for(u32 i=0;i<lim;i++){
u32 future=xd[i+24];if(future>=N)future=0;
__asm__ volatile("prefetchw %0" :: "m"(*(const char*)(hc+future)));
u32 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;for(u32 v=0;v<=N;v++){ct[v]=sum;sum+=hc[v];}
if(sum==N && N<1048576u){
const u64 bucket_reciprocal=((1ull<<40)+B-1u)/B;
u32*gt=s_gt[d];
for(u32 i=0;i<N;i++){
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;gt[i]=st;
u32 end=ct[v+1]&0xfffffu;
u32 bb=(u32)(((u64)(end-1u)*bucket_reciprocal)>>40);
s_bpos[(size_t)d*((size_t)N+64)+i]=(u8)(bb>Q-1?Q-1:bb);
}
}else{
for (u32 i = 0; i < N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i < N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tst[i] = s; s += c; } s_tst[N] = s; }
u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
u32 pos = s_cntv[v]++;
od[pos] = i;
gt[i] = s_tst[v];
u32 bb = (s_tst[v + 1] - 1u) / B;
s_bpos[(size_t)d * ((size_t)N + 64) + i] = (u8)(bb > Q - 1 ? Q - 1 : bb);
}
}
#if MF_FB2
{ for (u32 j = 0; j < N; j++) { u32 p = od[j]; s_posp[(size_t)p * MAXK + d] = j; } }
#else
if (g_use_filt) { u32 *pp = s_pos[d]; for (u32 j = 0; j < N; j++) pp[od[j]] = j; }
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
for (u32 i = 0; i < N; i++) {
u8 *row = s_bd + (size_t)i * BD_STRIDE;
for (int d = 0; d < k; d++)
row[d] = s_bpos[(size_t)d * ((size_t)N + 64) + i];
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < STRIDE; d++) p[d] = 0;
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
{ u32 dn = (N > 1) ? (N - 1) : 1u;
g_sigmag = (u32)((((u64)255 << 32) / (u64)dn) + 1u); }
size_t need = (size_t)k * k * ((size_t)N + 64);
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) {
__asm__ volatile("nop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\t");
_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();
}
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
__builtin_memcpy(s_bd + (size_t)i * BD_STRIDE + 8, &best, sizeof(best));
if (g_plen_on) {
const u8 *bdv = s_bd + (size_t)i * BD_STRIDE;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
{
for (u32 v = 0; v <= (u32)k; v++) s_cntv[v] = 0;
for (u32 i = 0; i < N; i++) s_cntv[s_g[i] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)k; v++) { u32 c = s_cntv[v]; s_cntv[v] = ac; ac += c; } }
for (u32 i = 0; i < N; i++) { u32 id = s_qord[i]; s_tmp[s_cntv[s_g[id] + 1]++] = id; }
for (u32 i = 0; i < N; i++) s_qord[i] = s_tmp[i];
}
u32 gstart[MAXK + 2];
{ u32 c = 0; for (int d = 0; d <= k; d++) { gstart[d] = c; while (c < N && (int)s_g[s_qord[c]] == d) c++; } gstart[k + 1] = c; }
#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 * ((size_t)N + 64);
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); sb[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
if (s_sig && !g_sig_merged && k == 6) {
static u8 cached_signature[(size_t)MAXN * 8] __attribute__((aligned(64)));
for (u32 i = 0; i < N; i++)
for (int e = 0; e < 6; e++) {
u32 v = (u32)(((u64)x[e][i] * g_sigmag) >> 32);
cached_signature[(size_t)i*8 + e] = (u8)(v ? v : 1u);
}
for (int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*k*((size_t)N + 64);
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);
_mm_storel_epi64((__m128i *)(base+(size_t)0*((size_t)N + 64)+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*((size_t)N + 64)+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*((size_t)N + 64)+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*((size_t)N + 64)+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*((size_t)N + 64)+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*((size_t)N + 64)+j),_mm_srli_si128(v2,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*8;
for(int e=0;e<k;e++)base[(size_t)e*((size_t)N + 64)+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 * ((size_t)N + 64);
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++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); base[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
}
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 = folded_rank(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
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
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 = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) {
if (ptr + 64 < MAXN) __builtin_prefetch(s_nid + od[ptr + 64], 0, 3);
u32 q = s_nid[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
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
template<int K> static void solve_mf_sp(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;
g_Q = q;
g_B = (N + q - 1) / q;
}
const u32 Q = g_Q, B = g_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
#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];
u8*hc=compact_hist;u32*ct=compact_cursor;
for(u32 i=0;i<=N;i++)hc[i]=0;
u32 lim=N>24?N-24:0;
for(u32 i=0;i<lim;i++){
u32 future=xd[i+24];if(future>=N)future=0;
__asm__ volatile("prefetchw %0" :: "m"(*(const char*)(hc+future)));
u32 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;for(u32 v=0;v<=N;v++){ct[v]=sum;sum+=hc[v];}
if(sum==N && N<1048576u){
const u64 bucket_reciprocal=((1ull<<40)+B-1u)/B;
u32*gt=s_gt[d];
for(u32 i=0;i<N;i++){
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;gt[i]=st;
u32 end=ct[v+1]&0xfffffu;
u32 bb=(u32)(((u64)(end-1u)*bucket_reciprocal)>>40);
s_bpos[(size_t)d*((size_t)N+64)+i]=(u8)(bb>Q-1?Q-1:bb);
}
}else{
for (u32 i = 0; i < N; i++) s_cntv[i] = 0;
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; s_cntv[v < N ? v : 0]++; }
{ u32 s = 0; for (u32 i = 0; i < N; i++) { u32 c = s_cntv[i]; s_cntv[i] = s; s_tst[i] = s; s += c; } s_tst[N] = s; }
u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
u32 pos = s_cntv[v]++;
od[pos] = i;
gt[i] = s_tst[v];
u32 bb = (s_tst[v + 1] - 1u) / B;
s_bpos[(size_t)d * ((size_t)N + 64) + i] = (u8)(bb > Q - 1 ? Q - 1 : bb);
}
}
#if MF_FB2
{ for (u32 j = 0; j < N; j++) { u32 p = od[j]; s_posp[(size_t)p * MAXK + d] = j; } }
#else
if (g_use_filt) { u32 *pp = s_pos[d]; for (u32 j = 0; j < N; j++) pp[od[j]] = j; }
#endif
BVAL(d, 0) = 0; FPOS(d, 0) = 0;
for (u32 b = 1; b < Q; b++) {
u32 pos = b * B < N ? b * B : N - 1;
u32 v = xd[od[pos]];
BVAL(d, b) = v;
u32 gs = pos;
while (gs > 0 && xd[od[gs - 1]] == v) gs--;
FPOS(d, b) = gs;
}
BVAL(d, Q) = 0xFFFFFFFFu; FPOS(d, Q) = 0;
}
for (u32 i = 0; i < N; i++) {
u8 *row = s_bd + (size_t)i * BD_STRIDE;
for (int d = 0; d < k; d++)
row[d] = s_bpos[(size_t)d * ((size_t)N + 64) + i];
}
for (u32 i = 0; i < N; i++) {
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) p[d] = x[d][i];
for (int d = k; d < STRIDE; d++) p[d] = 0;
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
{ u32 dn = (N > 1) ? (N - 1) : 1u;
g_sigmag = (u32)((((u64)255 << 32) / (u64)dn) + 1u); }
size_t need = (size_t)k * k * ((size_t)N + 64);
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) {
__asm__ volatile("nop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\tnop\n\t");
_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();
}
{
size_t mrn = (size_t)MAXK * MAXK * MAXQ;
for (size_t z = 0; z < mrn; z++) s_maxr[z] = 0;
}
for (u32 i = 0; i < N; i++) {
u32 best = 0xFFFFFFFFu; int ds = 0;
for (int d = 0; d < k; d++) {
u32 g = s_gt[d][i];
if (g < best) { best = g; ds = d; }
}
s_g[i] = (u32)ds;
__builtin_memcpy(s_bd + (size_t)i * BD_STRIDE + 8, &best, sizeof(best));
if (g_plen_on) {
const u8 *bdv = s_bd + (size_t)i * BD_STRIDE;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
{
for (u32 v = 0; v <= (u32)k; v++) s_cntv[v] = 0;
for (u32 i = 0; i < N; i++) s_cntv[s_g[i] + 1]++;
{ u32 ac = 0; for (u32 v = 0; v <= (u32)k; v++) { u32 c = s_cntv[v]; s_cntv[v] = ac; ac += c; } }
for (u32 i = 0; i < N; i++) { u32 id = s_qord[i]; s_tmp[s_cntv[s_g[id] + 1]++] = id; }
for (u32 i = 0; i < N; i++) s_qord[i] = s_tmp[i];
}
u32 gstart[MAXK + 2];
{ u32 c = 0; for (int d = 0; d <= k; d++) { gstart[d] = c; while (c < N && (int)s_g[s_qord[c]] == d) c++; } gstart[k + 1] = c; }
#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;
if (KSx == 8 && k <= 8) {
_mm256_store_si256((__m256i *)dp, _mm256_load_si256((const __m256i *)row));
} else if (KSx == 6 && k == 6) {
const __m256i rv = _mm256_loadu_si256((const __m256i *)row);
_mm_storeu_si128((__m128i *)dp, _mm256_castsi256_si128(rv));
_mm_storel_epi64((__m128i *)(dp + 4), _mm256_extracti128_si256(rv, 1));
} else 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 * ((size_t)N + 64);
for (int e = 0; e < k; e++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); sb[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
if (s_sig && !g_sig_merged && k == 6) {
static u8 cached_signature[(size_t)MAXN * 8] __attribute__((aligned(64)));
for (u32 i = 0; i < N; i++)
for (int e = 0; e < 6; e++) {
u32 v = (u32)(((u64)x[e][i] * g_sigmag) >> 32);
cached_signature[(size_t)i*8 + e] = (u8)(v ? v : 1u);
}
for (int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*k*((size_t)N + 64);
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);
_mm_storel_epi64((__m128i *)(base+(size_t)0*((size_t)N + 64)+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*((size_t)N + 64)+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*((size_t)N + 64)+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*((size_t)N + 64)+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*((size_t)N + 64)+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*((size_t)N + 64)+j),_mm_srli_si128(v2,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*8;
for(int e=0;e<k;e++)base[(size_t)e*((size_t)N + 64)+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 * ((size_t)N + 64);
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++) { u32 z_ = (u32)(((u64)row[e] * g_sigmag) >> 32); base[(size_t)e * ((size_t)N + 64) + j] = (u8)(z_ ? z_ : 1u); }
}
}
}
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 = folded_rank(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
{
const u32 *od = s_ord[ds];
for (u32 j = 0; j < N; j++) s_nid[od[j]] = j;
}
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 = (b < Q) ? FPOS(d, b) : (u32)N;
while (ptr < pend) {
if (ptr + 64 < MAXN) __builtin_prefetch(s_nid + od[ptr + 64], 0, 3);
u32 q = s_nid[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
}
switch (k) {
case 4: query_group<4>(N, out, ds, q0, q1, Q, Wpool); break;
case 5: query_group<5>(N, out, ds, q0, q1, Q, Wpool); break;
case 6: query_group<6>(N, out, ds, q0, q1, Q, Wpool); break;
case 7: query_group<7>(N, out, ds, q0, q1, Q, Wpool); break;
case 8: query_group<8>(N, out, ds, q0, q1, Q, Wpool); break;
case 9: query_group<9>(N, out, ds, q0, q1, Q, Wpool); break;
default: query_group<10>(N, out, ds, q0, q1, Q, Wpool); break;
}
}
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 4); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 5); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf_sp<6>((u32)n, (const unsigned **)x, out); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 7); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 9); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 10); }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 430.574 ms | 65 MB + 124 KB | Accepted | Score: 100 | 显示更多 |