/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/126078: retained the scalar ninth-threshold queue representation.
* saffah_codex_6a_agg3, https://duck.ac/submission/125974: adapted six-dimensional17-bit threshold compression into low16-bit lanes plus high bits in the query handle.
* saffah_cc_v41_agg1, https://duck.ac/submission/122417: retained original engine and all inherited references.
* Idea: The eight remaining thresholds are17-bit values for this100000-point
* task. Pack their low16bits into16bytes and their high bits into unused bits
* of the existing query index, once per query. Replay reconstructs the exact
* full-width vector with one variable shift and OR. This halves the combined
* threshold storage; other dimensional instantiations retain the old path.
* Purpose: Official experiment reducing ambiguity-queue footprint and traffic.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/122417: copied signature/min-fold engine and its deferred ambiguity queue; retained original references.
* Idea: A nine-dimensional ambiguous candidate only uses lane zero of its
* second threshold vector. Save that lane in a four-byte queue and compare it
* directly on replay instead of storing/loading 32 bytes of mostly neutral
* padding. Other dimensions keep the full vector path.
* Purpose: Official experiment removing dead high-vector queue traffic.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #120285
// <https://duck.ac/submission/120285>（判题机 Time = 100.554973 ms，本题当前榜首）
// —— 本文件**以该提交的正文逐字复刻为基座**（含其文件头自带的英文引用/思路区，原样保留）。
// 所用内容：其"签名面建趟（8 点分块字节转置 + 按 d 编译期特化的 switch 分派）"整块，
// 以及它从本账号 #111378 / #114343 继承的全部引擎（min-fold 池、AVM 记录合并、桶缓存等）。
// 原提交正文无独立许可证声明；duck.ac 提交正文按站点规则公开可见（题面"你提交的代码将会被公开"）。
// [2] duck.ac 本账号（saffah_cc_v41_agg1），提交 #111378
// <https://duck.ac/submission/111378>（101.785093 ms）= 本文件改动前的我方最好件，
// 被 [1] 逐字节继承（#120246）+ 一处真改动（#120285）。本文件是它的下游。
// [3] duck.ac 本账号（saffah_cc_v41_agg1），提交 #114343
// <https://duck.ac/submission/114343>（1014b）—— [1] 的建趟刀本体出处（8 点字节转置）。
// [4] 其余上游（#109755 saffah_codex_6s_agg2 / #109677 / #107496 / #107259 与
// Granlund–Montgomery PLDI 1994）沿 [1] 文件头内的英文引用区保留，未作改动。
// 合规：按 /home/yjp/duck.ac/RULES.md 第 3 节署名各作者账号、原提交地址与所用内容。
// ======================
// ===== 思路 =====
// 【席位 b15y_ · 第一发 [B15Y-COPY]（已发 #122326 = 100.551281 ms Accepted）】
// §2.18.349 第一步：把对手最新件 #120285 逐字复刻为基座（正文与 #120285 逐字节相同）。
// 实测 100.551281 ms（对手 T = 100.554973），即"纯抄 = 平手"；距严支线 99.550423 仍缺 1.0009 ms
// ⇒ 按 §2.18.349 继续"抄 + 叠刀"。
// 【本发（第二发）· [B15Y-CPL] 编译期平面索引 —— 本次唯一自变量】
// · 观测：签名面 `s_sig` 为每个扫描维 d 都存了 k 份平面副本，而 d 那一趟**只用** e != d 的
// k-1 份；旧读取循环因此写成 `for e in 0..k-1 { if (e == d) continue; ... }`，而 d 是运行期
// 变量 ⇒ 每 32 候选的循环体里留下 **8 个 cmp+jcc**，8 条平面 load/sub/or 被切成带分支的
// 锯齿链（判题机 gcc9 汇编逐条核过：.L3879 族每条 vmovdqu 后跟 `cmpq $N,%r15; je .Lxxxx`）。
// · 改法：把 dim d 的签名块改成**紧凑存放** —— 只存 e != d 的 k-1 份，顺序 t = e-(e>d)；
// 读取侧 `for t in 0..k-2` 直接 t*MF_SPL，**无分支、无额外 load**；qsb[] 同序重排（每
// (query,dim) 的 setup 里做一次 e = t + (t >= d) 的映射，不进热循环）。因 e == d 那一份在
// 旧代码里被跳过、qsb[d] 恒为 0xFF（subs_epu8(x,0xFF) ≡ 0 恒成立），该改写**逐位等价**。
// · 顺带把 stride 由 k*k*MF_SPL 收到 k*(k-1)*MF_SPL（签名面 8.1 MB → 7.2 MB）；单独量该收缩
// 为 0，故本发全部收益来自索引改写（见下）。
// · 闸门：判题机同款编译器 ref/gcc9/g9.sh + 判题机路径驱动 work/b15x_tmp/rc_1015b_jmain.cpp，
// n=1e5 mod=0 校验和 **19836411**（历史生产点）✓；判题机 probe 侧 FNV
// **h=9900767564265552992** 与基座逐位相同 ✓。
// · 判题机定价（tools/probe.py，同窗口）：基座 102.636842/102.639822 ⇒ 本件
// 102.407230/102.408947/102.483354 ⇒ **−0.2296 ms（−0.224%）**（探针 σ≈0.012% ⇒ 18σ）。
// · 同窗口一并判负（未采纳）：MF_PFB3 16 ⇒ 102.990（+0.34%）· 24 ⇒ 103.211（+0.56%）·
// Q=510/1020 ⇒ 104.92/104.90（+2.2%）· 足迹收缩并入版 ⇒ +0.009（噪声内）。
// ================
// 【席位 b15y_ · [B15Y-COPY]】§2.18.349 第一步：把对手最新件逐字复刻为基座。
// · 目的（试验性）：定位我方与榜首的差距，作为后续叠刀的起点与对照臂。
// · 本文件除本引用/思路区外，正文与 #120285 **逐字节相同**（md5 见 notes）。
// · 闸门：判题机同款编译器 + 判题机路径驱动，n=1e5 mod=0 校验和 19836411（历史生产点）。
// ================
// 【席位 b15b_ · [B15B-TAILB] 精修热环尾部块改走"块环自己的基址" —— 本次唯一自变量】
// · 观测（本席把该相位按 rdtscp 拆成 setup / 块环 / 尾部三段，n=1e5 判题机同形，795k 次扫描）：
// 块环 4,505,013 次 @ ~14.9 拍/块 · 尾部 **770,416 次 @ ~24.9 拍** · 扫描 setup ~26 拍
// （sumR = 156,340,565 候选；尾部只处理 12.1M 候选 = 7.7% 的字节却吃掉 ~18% 的相位）。
// · 根因（判题机 gcc9 汇编逐条核过）：尾部为了取"以 hi 结尾的 32 B 窗"写了 `base = hi-32`
// + 左移掩码，与块环的 `base = lo+32*q` 是**两个不同的基址** ⇒ gcc9 在尾部重新生成 8 条
// `lea <disp>(%rax),%rXX` 并把它们**溢出到栈**再逐条重载（`mov 0x2f0(%rsp),%rbx` 族），
// 即每次尾部多付 ~30 条 uop + 一次冷代码路径。
// · 改法：尾部**继续块环自己的 aj**（`base = aj, lk = 0`，掩码改为低 r 位）——窗口块数
// 完全不变（仍是 ceil(R/32)），但 8 条平面地址保持块环的 `base+disp+aj` 形式，gcc9 不再
// 物化/溢出指针（新汇编中 `lea 0x186e0(%rax)` 族已消失）。
// · 越读安全：MF_SPL = N + 64（= 100064 > 100000），尾部最多多读 31 B，落在平面内 64 B 的
// slack 里；读到的那部分由掩码清零，**逐位等价**。
// · 闸门：判题机同款编译器 ref/gcc9/g9.sh + 判题机路径驱动 rc_1015b_jmain.cpp：
// n=1e5/mod=0 cs **19836411**、n=1e5/mod=1000 cs 18684868、n=65536 cs 8309240、
// n=150000 cs 44789182 —— 四组全与在册件 #122360 逐位相同 ✓；
// 判题机 probe 侧 FNV **h=9900767564265552992** 与在册件逐位相同 ✓。
// · 判题机定价（tools/probe.py --reps 3，同窗口）：在册件 #122360 形
// 102.469805/102.440494/102.433609 ⇒ 本件 101.590396/101.591466/101.537996
// ⇒ **−0.8957 ms（−0.874%）**（同件散布 0.05%，两窗不相交 ⇒ 约 30σ）。
// · 同批判负（未采纳）：ASCR and_popcount 主循环 8 词展开（逐位等价，但本机 3/4 组更慢，
// 判题机前景不明，弃）；vpermd+vpshufb 阈值广播改 vpbroadcastb 栈暂存（逐位等价，
// 本机 setup −4% 但块环 +9%，净负，弃）；MF_PFB3 = 0（本机 min-of-6 190.6 → 248.1 ms，
// 大幅转负 ⇒ 该族预取是净收益，**同时否证"本环与访存无关"**）。
// ================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/111378:
* copied the accepted signature/min-fold engine and preserved its original credits.
* saffah_cc_v41_agg1, https://duck.ac/submission/114343:
* adapted the eight-point blocked byte-transpose idea for signature-plane output.
* Both public sources displayed no separate software license notices.
* Additional idea: Dispatch the omitted dimension once outside the blocked loop,
* specializing each plane emitter so its per-block dimension tests disappear.
* Idea: Replace scalar signature-plane scattering with an eight-record SSE byte
* transpose and 64-bit stores, retaining the existing compact cached signatures
* and excluding the ordering dimension exactly as in the original emitter.
* Purpose: Official experiment reducing plane-build store instructions and indexing.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #109755
// <https://duck.ac/submission/109755>（判题机 Time = 103.797941 ms）
// —— 本文件**以该提交的正文为基座**（其自带的整条引用/思路注释区在下方原样保留）。
// [2] duck.ac 本账号（saffah_cc_v41_agg1）提交 #109782
// <https://duck.ac/submission/109782>（103.770878 ms）= #109755 的逐字节副本 + 引用头，
// 即本文件在改动之前的形态，也是本次判题机读数的对照臂。
// [3] duck.ac 题目页 <https://duck.ac/problem/1015b>（公开约束），仅用于读取公开约束，未复制代码。
// 合规：duck.ac 提交正文按站点规则公开可见（题面"你提交的代码将会被公开"），可直接取用；
// 原提交正文无独立许可证声明。本文件按 /home/yjp/duck.ac/RULES.md 第 3 节署名作者账号与原提交地址。
// ======================
// ===== 思路 =====
// 【席位 a15y_（本席新工作）· [A15Y-AVM]】消掉 pass B 每 (查询,维) 的**第三条随机流**
// （与 1015a #111338 同一把刀，那边实测 −2.16%）。三条流 = s_pm 记录（6.4 MB 随机 64 B）、
// s_gt[d]（400 KB）、以及 per-(ds,d) 粗向量数组 avpool[]（3.2 MB 随机 32 B）。
// · 关键观测：`Avm = Av | notexlo`（exmask = {e < d} \ {ds}）⇒ 粗向量 Av[e] **只在 lane
// e < d ∧ e != ds 上可观测**；记录 64 B 的 **lane 9..15 从不被读取**（lane 8 = coord[8] 由
// qrh 读；9..15 被 khi8hv 置 0x7FFFFFFF，k=9 时不可达）。
// · 改法：把 m[e] = min(coord[e], Av[e]) 直接写进记录的空闲 lane 9+e（e = 0..6，**同一 cache
// line、零额外足迹**），m[7] 写进 s_m7[]（仅在 d == 8 那一趟需要，故按 d==8 条件取用）。
// 读侧用一条 vpermd 把 lane 9..15 对齐到 lane 0..6。因 min(x, min(x,a)) == min(x,a)，
// 该改写**逐位等价**（不是近似）。三条随机流 → 两条，avpool 的建趟也随之消失。
// · 闸门：判题机同款编译器 ref/gcc9/g9.sh + 判题机路径驱动 work/rc_1015b_jmain.cpp；
// 逐位对拍 n = 100000 / 20000 / 2000 × mode 0 / 3 = **6/6 SAME**
// （n=1e5 mode0 checksum 19836411，与历史生产点一致）。
// · 正控（验证闸门读数真变）：变异体 pc1（去掉 lane 0..6 的粗夹）cs 759915 → 760125；
// pc2（去掉 lane 7 的粗夹）759915 → 759920 ⇒ 合并进去的 lane 与 s_m7 **确被读取** ✓
// ================
// References:
// [1] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/109677 : copied the accepted 1015b min-fold engine, including ping-pong pass-B order buffers, build and output prefetches, compact signature rows, bucket cache, and current parameters. The public submission gives attribution to its ancestors; it displayed no separate software license notice.
// [2] Our account saffah_codex_6s_agg2, https://duck.ac/submission/107496 : compact signature-row preprocessing retained within [1].
// [3] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/107259 : prior min-fold, byte signature, and deferred fallback engine inherited through [1].
// [4] Torbjoern Granlund and Peter L. Montgomery, "Division by Invariant Integers using Multiplication," PLDI 1994, https://dl.acm.org/doi/10.1145/178243.178249 : arithmetic principle inherited through [3]; no paper text or software is copied.
// The cited duck.ac submissions displayed no separate software license notice. Each reused source or idea, author, and original link is listed above.
// Approach: For the unfiltered production shape, compute rank-boundary and coarse bucket values while scattering the coordinate order. This removes a later full coordinate pass; the filtered shape retains its original build loop.
// Purpose: Experimentally price a fused coordinate-order and bucket-construction pass on duck.ac.
#define MF_Q 255
#define MF_PFB 20
#define MF_PFB1 64
#define MF_PFB2 64
#define MF_PFB3 8
#define MF_PFB_OUT 4
#define MF_PFA 2
#define MAXQ 256
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing,align-functions=64")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#define MF_SPL ((size_t)100064)
#include <immintrin.h>
static inline __m256i LDU(const void *p) {
__m256i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline __m256i LDA(const void *p) { // p 必须 32 字节对齐
__m256i v;
__asm__("vmovdqa %1, %0" : "=x"(v) : "m"(*(const __m256i *)p));
return v;
}
static inline void STU(void *p, __m256i v) {
__asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline void STA(void *p, __m256i v) { // p 必须 32 字节对齐
__asm__ volatile("vmovdqa %1, %0" : "+m"(*(__m256i *)p) : "x"(v));
}
static inline __m128i LDQ(const void *p) {
__m128i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m128i *)p));
return v;
}
#define LDU_M(p) ({ __m256i _v; __asm__("vmovdqu %1, %0" : "=x"(_v) : "m"(*(const __m256i *)(p))); _v; })
#define LDA_M(p) ({ __m256i _v; __asm__("vmovdqa %1, %0" : "=x"(_v) : "m"(*(const __m256i *)(p))); _v; })
#define STU_M(p, v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i *)(p)) : "x"(_vv)); } while (0)
#define STA_M(p, v) do { __m256i _vv = (__m256i)(v); __asm__ volatile("vmovdqa %1, %0" : "+m"(*(__m256i *)(p)) : "x"(_vv)); } while (0)
extern "C" void *malloc(unsigned long);
typedef unsigned u32;
typedef unsigned char u8;
typedef unsigned short u16;
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
#define MAXK 9
#define STRIDE 16
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 95
#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 1
#endif
// [a15y AVM] 1 = merge the per-(ds,d) coarse vector Av[] into the s_pm record's free
// high lanes 9..15 (m[0..6]) + s_m7[] (m[7]): removes one of the three random
// per-(query,dim) streams (s_gt / s_pm / avpool). k == 9 only.
#define MF_AVM 1
#ifndef MF_NOSDMP
#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 u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u32 s_m7[MAXN]; // [a15y AVM] m[7] = min(coord[7], BVAL(7, bucket[7]))
static u8 s_bd[(size_t)MAXN * MAXK];
static u32 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u16 s_outlo[MAXN];
static inline void add_out16(u32 i, u32 add) { s_outlo[i] = (u16)((u32)s_outlo[i] + add); }
static u8 s_bkt[MAXN];
static u8 s_bks[MAXN]; // [K3a] sort keys scattered into post-sort order
static u8 s_bka[MAXN]; // n17r: permuted d1-bucket cache (post-sort order)
#define AMBD_CAP 4096
static u32 s_abi[AMBD_CAP];
static const u32 *s_abp[AMBD_CAP];
static __m256i s_abqT[AMBD_CAP], s_abqTh[AMBD_CAP];
static u32 s_abq9[AMBD_CAP];
alignas(64) static __m128i s_abqLow[AMBD_CAP];
static u32 s_nid[MAXN];
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN]; // s_posp[p*MAXK+d] = rank of p in dim d
static u64 s_facc[(size_t)(MAXK - 1) * (((MAXN) + 63) / 64 + 16)];
static u32 s_bnde[(size_t)MAXK * (MAXQ + 2)]; // per-dim coarse-cell boundary positions
static 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)));
static u32 s_maxr[(size_t)MAXK * MAXK * MAXQ];
static u8 *s_sig = 0;
static u8 a15b_sigrow[(size_t)MF_SPL * STRIDE] __attribute__((aligned(64)));
static int g_sig_on = 0;
static int g_sig_merged = 0;
static u32 g_sigsh = 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; // row stride of the transposed coordinate rows (64 B for k>8)
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
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, acc = z;
u32 w = 0, cnt = 0, total = 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);
acc = _mm256_add_epi8(acc, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
if (++cnt == 31) { tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z)); acc = z; cnt = 0; }
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc, z));
total = hsum256(tot);
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 1; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v);
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 1; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & ((1ull << rem) - 1));
}
return total;
}
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 int k = KK;
const u32 ASTRL = (KK <= 9) ? 8u : 16u;
#if MF_AVM
const bool MERGE = (KK == 9); // record carries Av in lanes 9..15 (m[0..6]) + s_m7
#else
const bool MERGE = false;
#endif
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_PFB1 > 0
if (qi + MF_PFB1 < q1) _mm_prefetch((const char *)(s_bd + (size_t)s_qord[qi + MF_PFB1] * MAXK + d1), _MM_HINT_T0);
#endif
{ u32 n17r_bk = s_bd[(size_t)s_qord[qi] * MAXK + d1];
s_bkt[qi] = n17r_bk; cntv[n17r_bk + 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]; u32 n17r_p = cntv[s_bkt[qi] + 1]++;
s_tmp[n17r_p] = id; s_bka[n17r_p] = s_bkt[qi]; }
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] * MAXK + dd] + 1]++;
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) { u32 id = s_qord[qi]; s_tmp[cntv[s_bd[(size_t)id * MAXK + dd] + 1]++] = id; }
for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = s_tmp[qi];
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = s_qord[qi + PF_DIST];
const u8 *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
#if MF_ASCR
{
u32 qi = q0;
while (qi < q1) {
u32 b1 = s_bka[qi];
u32 qe = qi + 1;
while (qe < q1 && s_bka[qe] == b1) qe++;
const u64 *bp1 = g_plen_on ? (s_bs + s_pofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
for (u32 qq = qi; qq < qe; qq++) {
u32 i = s_qord[qq];
const u8 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
if (e == (int)e1) bp[e] = bp1;
else { int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W); }
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
s_outlo[i] = (u16)cnt;
}
qi = qe;
}
}
#else
for (u32 qi = q0; qi < q1; qi++) {
u32 i = s_qord[qi];
const u8 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = s_qord[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = s_qord[qi + PF_DIST];
const u8 *bdn = s_bd + (size_t)j * MAXK;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const char *pp = (const char *)(s_bs + ((size_t)e * (Q + 1) + bdn[d]) * W);
_mm_prefetch(pp, _MM_HINT_T0);
_mm_prefetch(pp + 64, _MM_HINT_T0);
_mm_prefetch(pp + 128, _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
if (R && !g_skip_and) {
const u64 *bp[MAXK];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
bp[e] = g_plen_on ? (s_bs + s_pofs[(size_t)e * MAXQ + bd[d]])
: (s_bs + ((size_t)e * (Q + 1) + bd[d]) * W);
#ifdef DBG_PLEN
if (g_plen_on) {
u32 need = (R >> 6) + ((R & 63) ? 1 : 0);
if (need > s_plen[(size_t)e * MAXQ + bd[d]])
printf("DBG ds=%d e=%d d=%d b=%u R=%u need=%u len=%u\n", ds, e, d, bd[d], R, need, s_plen[(size_t)e * MAXQ + bd[d]]);
}
#endif
}
cnt = and_popcount_range<KK - 1>(bp, R >> 6, R & 63u);
}
s_outlo[i] = (u16)cnt;
}
#endif
if (g_skip_fringe) return;
static u32 *avpool = 0; static size_t avcap = 0;
{
size_t needav = (size_t)N * STRIDE * 4 + 256;
if (needav > avcap) { avpool = (u32 *)pool_alloc(needav); avcap = avpool ? needav : 0; }
if (avpool && !MERGE) for (u32 qi = q0; qi < q1; qi++) {
u32 ii = s_qord[qi];
const u8 *bdi = s_bd + (size_t)ii * MAXK;
u32 *av = avpool + (size_t)ii * ASTRL;
for (int e = 0; e < k && e < (int)ASTRL; e++) av[e] = (e == (int)ds) ? 0x7FFFFFFFu : BVAL(e, bdi[e]);
for (int e = k; e < (int)ASTRL; e++) av[e] = 0x7FFFFFFFu;
}
}
u32 *qcur = s_qord, *qdst = s_tmp;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
#if !MF_FR_NOSORT
{
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) {
#if MF_PFB2 > 0
if (qi + MF_PFB2 < q1) _mm_prefetch((const char *)(s_bd + (size_t)qcur[qi + MF_PFB2] * MAXK + d), _MM_HINT_T0);
#endif
u32 bb = s_bd[(size_t)qcur[qi] * MAXK + d];
s_bkt[qi] = bb;
cntv[bb + 1]++;
}
{ u32 ac = q0; for (u32 v = 0; v <= Q; v++) { u32 c = cntv[v]; cntv[v] = ac; ac += c; } }
for (u32 qi = q0; qi < q1; qi++) {
u32 id = qcur[qi];
u32 p = cntv[s_bkt[qi] + 1]++;
qdst[p] = id;
s_bks[p] = s_bkt[qi]; // [K3a] carry the key through the permutation
}
{ u32 *_t = qcur; qcur = qdst; qdst = _t; }
}
#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 __m256i idxm = _mm256_setr_epi32(1, 2, 3, 4, 5, 6, 7, 0); // [a15y AVM]
int amb_n = 0;
const u32 *const odd = s_ord[d];
const u8 *const splane = s_sig + (size_t)d * k * MF_SPL;
for (u32 qi = q0; qi < q1; qi++) {
u32 i = qcur[qi];
#if MF_PFB3 > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = qcur[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);
if (!MERGE) _mm_prefetch((const char *)(avpool + (size_t)j * ASTRL), _MM_HINT_T0);
else if (d == 8) _mm_prefetch((const char *)(s_m7 + j), _MM_HINT_T0);
{
u32 qo = qi + MF_PFB_OUT;
if (qo < q1) {
u32 jo = qcur[qo];
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(s_outlo + jo)));
} else {
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(s_outlo + j)));
}
}
}
}
#endif
const u32 lo = FPOS(d, s_bks[qi]); // [K3a] sequential key, not a cold s_bd line
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(LDA((const __m256i *)rec), khi8v);
__m256i Avm;
if (MERGE) {
// [a15y AVM] h = [coord8, m0 .. m6]; permute -> [m0..m6, *]. Lane 7 is only
// ever unmasked on the d == 8 sweep, so the blend below is exact.
__m256i h = LDU256((rec + 8));
__m256i msh = _mm256_permutevar8x32_epi32(h, idxm);
if (d == 8) msh = _mm256_blend_epi32(msh, _mm256_set1_epi32((int)s_m7[i]), 0x80);
Avm = _mm256_or_si256(msh, notexlo);
} else {
Avm = _mm256_or_si256(LDU256((avpool + (size_t)i * ASTRL)), notexlo);
}
__m256i qT = _mm256_min_epi32(qr, Avm);
__m256i qTh;
if (k > 9) {
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_or_si256(LDU256((avpool + (size_t)i * ASTRL + 8)), notexloh);
qTh = _mm256_min_epi32(qrh, Avmh);
} else if (k > 8) {
qTh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
} else {
qTh = all8;
}
u32 add = 0;
const u32 *c0 = s_sdmp[d] ? (s_sdmp[d] + (size_t)lo * g_KS) : 0;
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const 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) * g_KS) : (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;
}
}
add_out16(i, add);
continue;
}
if (g_sig_on && (hi - lo) >= 2) {
__m256i qsb[16];
{
const __m256i zero8 = _mm256_setzero_si256();
const __m256i shv = _mm256_set1_epi32((int)g_sigsh);
const __m256i one8 = _mm256_set1_epi8(1);
__m256i sw = _mm256_subs_epu8(_mm256_srlv_epi32(qT, shv), one8);
__m256i swh = _mm256_subs_epu8(_mm256_srlv_epi32(qTh, shv), one8);
// [B15Y-CPL] compact plane index: dim d's signature block stores the
// k-1 planes for dims != d back to back (t = e - (e > d)), so the scan
// below needs no per-plane `e == d` test and no branch at all.
for (int t = 0; t < k - 1; t++) {
int e = (t >= (int)d) ? t + 1 : t;
__m256i src = (e < 8) ? sw : swh;
__m256i bc = _mm256_permutevar8x32_epi32(src, _mm256_set1_epi32(e & 7));
qsb[t] = _mm256_shuffle_epi8(bc, zero8);
}
}
__m128i queue_low;u32 queue_i=i;
if constexpr(KK==9){
__m256i lows=_mm256_and_si256(qT,_mm256_set1_epi32(65535));
queue_low=_mm_packus_epi32(_mm256_castsi256_si128(lows),_mm256_extracti128_si256(lows,1));
queue_i|=(u32)_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_slli_epi32(qT,15)))<<20;
}
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
__m256i ac = _mm256_setzero_si256();
{
for (int t = 0; t < k - 1; t++) {
__m256i cv = LDU256((splane + (size_t)t * MF_SPL + aj));
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[t]));
}
}
u32 decm = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_setzero_si256()));
u32 amb = (u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)));
aacc += (u32)__builtin_popcount(decm);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = s_pm + (size_t)odd[aj + (u32)t] * STRIDE;
if (amb_n < AMBD_CAP) {
s_abi[amb_n] = queue_i;
s_abp[amb_n] = cp;
if constexpr(KK==9) _mm_store_si128(&s_abqLow[amb_n],queue_low);
else _mm256_store_si256(&s_abqT[amb_n],qT);
if constexpr(KK==9) s_abq9[amb_n]=(u32)_mm_cvtsi128_si32(_mm256_castsi256_si128(qTh));
else if(k>8) _mm256_store_si256(&s_abqTh[amb_n],qTh);
amb_n++;
} else {
__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) {
// [B15B-TAILB] continue the block loop's own base instead of
// re-deriving hi-32: same number of 32-byte windows, but the
// eight plane addresses keep the loop's base+disp+aj form, so
// gcc9 no longer materialises and spills eight fresh pointers
// here. MF_SPL = N+64 supplies the slack this may read past hi.
u32 base = aj, lk = 0u;
u32 kmask = (1u << r) - 1u;
__m256i ac = _mm256_setzero_si256();
{
for (int t = 0; t < k - 1; t++) {
__m256i cv = LDU256((splane + (size_t)t * MF_SPL + base));
ac = _mm256_or_si256(ac, _mm256_subs_epu8(cv, qsb[t]));
}
}
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 = s_pm + (size_t)odd[base + (u32)t] * STRIDE;
if (amb_n < AMBD_CAP) {
s_abi[amb_n] = queue_i;
s_abp[amb_n] = cp;
if constexpr(KK==9) _mm_store_si128(&s_abqLow[amb_n],queue_low);
else _mm256_store_si256(&s_abqT[amb_n],qT);
if constexpr(KK==9) s_abq9[amb_n]=(u32)_mm_cvtsi128_si32(_mm256_castsi256_si128(qTh));
else if(k>8) _mm256_store_si256(&s_abqTh[amb_n],qTh);
amb_n++;
} else {
__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 + g_KS));
__m256i p2 = LDU256((c + 2 * g_KS));
__m256i p3 = LDU256((c + 3 * g_KS));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8);
c += 4 * g_KS;
}
for (; j < hi; j++, c += g_KS) {
__m256i pv = LDU256(c);
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
add = a0 + a1 + a2 + a3;
} else {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + g_KS));
__m256i p2 = LDU256((c + 2 * g_KS));
__m256i p3 = LDU256((c + 3 * g_KS));
__m256i h0 = LDU256((c + 8));
__m256i h1 = LDU256((c + g_KS + 8));
__m256i h2 = LDU256((c + 2 * g_KS + 8));
__m256i h3 = LDU256((c + 3 * g_KS + 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 * g_KS;
}
for (; j < hi; j++, c += g_KS) {
__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 = LDA((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDA((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
add_out16(i, add);
}
for (int z = 0; z < amb_n; z++) {
if (z + 12 < amb_n) __builtin_prefetch(s_abp[z + 12], 0, 0);
const u32 *cp = s_abp[z];
__m256i pv = LDU256(cp);
__m256i qv;
if constexpr(KK==9){
__m256i low=_mm256_cvtepu16_epi32(_mm_load_si128(&s_abqLow[z]));
__m256i high=_mm256_set1_epi32(s_abi[z]>>20);
high=_mm256_and_si256(_mm256_sllv_epi32(high,_mm256_setr_epi32(16,15,14,13,12,11,10,9)),_mm256_set1_epi32(65536));
qv=_mm256_or_si256(low,high);
}else qv=_mm256_load_si256(&s_abqT[z]);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qv, pv), all8);
if constexpr(KK==9) { ok &= cp[8] < s_abq9[z]; }
else if (k > 8) {
__m256i pvh = LDU256(cp + 8);
__m256i qvh = _mm256_load_si256(&s_abqTh[z]);
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qvh, pvh), all8);
}
add_out16(s_abi[z] & (KK==9 ? 0xfffffu:0xffffffffu),ok);
}
}
if (qcur != s_qord) { for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = qcur[qi]; }
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int k) {
__builtin_memset(s_outlo, 0, (size_t)N * sizeof(u16));
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++; // ceil(sqrt(n)) without libm
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}; // per-array fill, k=4..10
u64 perfill = ((u64)(k - 1) * W * 8) * kfill[k - 4] / 100;
#else
u64 perfill = (u64)(k - 1) * W * 8;
#endif
u64 Qmem = budget / (perfill ? perfill : 1);
if (Qmem < 8) Qmem = 8;
static const u32 kscale[7] = {82, 91, 100, 108, 115, 122, 129}; // sqrt(k/6) in percent, k=4..10
u32 q = (u32)((s2 * MF_QCOEF + 99) / 100);
q = (u32)(((u64)q * kscale[k - 4] + 99) / 100);
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
g_Q = q;
g_B = (N + q - 1) / q;
}
const u32 Q = g_Q, B = g_B;
const u64 bmag = ((1ull << 40) + B - 1) / B;
#if MF_PLEN
g_plen_on = 1; // independent of MF_WTRIM: lengths only, stride stays W
#endif
#if MF_WTRIM
u32 Wpool = W; // MF_WTRIM is KNOWN-BROKEN, see the warning above
#else
const u32 Wpool = W;
#endif
#if MF_FILTAUTO
if (!g_use_filt && B >= 400) g_use_filt = 1; // coarse grid: the filter beats the scan
#endif
if (g_use_filt && (u64)k * k * (g_Qc + 1) * W * 8 > ((u64)700 << 20)) g_use_filt = 0;
#if !MF_FILTAUTO
g_use_filt = 0;
#endif
{
int lo8[8], hi8[8];
for (int e = 0; e < 8; e++) lo8[e] = (e < k) ? 0 : 0x7FFFFFFF;
for (int e = 0; e < 8; e++) hi8[e] = (8 + e < k) ? 0 : 0x7FFFFFFF;
khi8v = LDU256(lo8);
khi8hv = LDU256(hi8);
}
const int K1 = k - 1;
for (int d = 0; d < k; d++) {
const unsigned *xd = x[d];
u32 *od = s_ord[d];
for (u32 i = 0; i <= N; i++) s_cntv[i] = 0;
{ const u32 B15Z_PD = 24;
for (u32 i = 0; i < N; i++) {
if (i + B15Z_PD < N) { u32 vf = xd[i + B15Z_PD]; if (vf >= N) vf = 0;
__asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(s_cntv + vf))); }
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_tmp[i] = s; s += c; } }
#if MF_FB2
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i;
s_posp[(size_t)i * MAXK + d] = q;
}
#else
if (g_use_filt) {
u32 *pp = s_pos[d];
for (u32 i = 0; i < N; i++) { u32 v = xd[i]; if (v >= N) v = 0; u32 q = s_tmp[v]++; od[q] = i; pp[i] = q; }
} else {
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
od[s_tmp[v]++] = i;
s_gt[d][i] = s_cntv[v];
u32 num = s_cntv[v + 1] - 1u; if (num > N) num = N;
u32 bb = (u32)(((u64)num * bmag) >> 40);
s_bd[(size_t)i * MAXK + d] = bb > Q - 1 ? Q - 1 : bb;
}
}
#endif
if (g_use_filt || MF_FB2) { u32 *gt = s_gt[d];
for (u32 i = 0; i < N; i++) {
u32 v = xd[i]; if (v >= N) v = 0;
gt[i] = s_cntv[v];
u32 num = s_cntv[v + 1] - 1u; if (num > N) num = N;
u32 bb = (u32)(((u64)num * bmag) >> 40);
s_bd[(size_t)i * MAXK + d] = bb > Q - 1 ? Q - 1 : bb;
} }
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++) {
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;
#if MF_AVM
// [a15y AVM] Av[e] = BVAL(e, bucket_e(i)) is only ever observed on lanes e < d
// (exmask), so for k == 9 the lanes e = 0..6 live in the record's never-read high
// lanes 9..15 (lane 8 = coord[8] is read by qr_h; 9..15 are killed by khi8hv) and
// e = 7 goes to s_m7[] (needed only on the d == 8 sweep). min(x, min(x,a)) == min(x,a).
if (k == 9) {
for (int e = 0; e < 7; e++) {
u32 av = s_bval[(size_t)e * MAXQ + s_bd[(size_t)i * MAXK + e]], q = p[e];
p[9 + e] = (av < q) ? av : q;
}
{ u32 av = s_bval[(size_t)7 * MAXQ + s_bd[(size_t)i * MAXK + 7]], q = p[7];
s_m7[i] = (av < q) ? av : q; }
}
#endif
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; } // smallest sh with (N-1)>>sh <= 255
g_sigsh = sh;
size_t need = (size_t)k * k * MF_SPL;
static u8 *sigpool = 0; static size_t sigcap = 0;
if (need > ((size_t)220 << 20) || (size_t)N + 32 > MF_SPL) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
g_sig_on = 1;
g_sig_merged = 0;
}
}
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), LDA((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_storeu_si256((__m256i *)(dst + w), LDU256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), LDA((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
{
size_t 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;
if (g_plen_on) {
const u8 *bdv = s_bd + (size_t)i * MAXK;
u32 *mr = s_maxr + ((size_t)ds * MAXK) * MAXQ;
u32 bb = best;
for (int d = 0; d < k; d++) {
if (d == ds) continue;
u32 b = bdv[d];
if (bb > mr[(size_t)d * MAXQ + b]) mr[(size_t)d * MAXQ + b] = bb;
}
}
}
for (u32 i = 0; i < N; i++) s_qord[i] = i;
{
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 (s_sig && k == 9 && (size_t)N <= MF_SPL) {
for (u32 i = 0; i < N; i++) {
u8 *dst = a15b_sigrow + (size_t)i * STRIDE;
for (int e = 0; e < 9; e++) {
u32 z_ = x[e][i] >> g_sigsh;
dst[e] = (u8)(z_ ? z_ : 1u);
}
}
}
#if MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 8; // 32-byte aligned candidate rows for k<=8
#endif
size_t per = (size_t)N * KS * 4;
size_t lim = (size_t)150 << 20;
if (k >= 9) lim = (size_t)470 << 20; // k=9/10 need the rows: the random s_pm fallback is fatal
if (per * k <= lim) {
static u32 *pool = 0; static size_t cap = 0;
if (per * k > cap) { pool = (u32 *)pool_alloc(per * k + 256); cap = per * k; }
if (pool) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb = s_sig + (size_t)d * k * MF_SPL;
for (int e = 0; e < k; e++) if (e != (int)d) { u32 z_ = row[e] >> g_sigsh; sb[(size_t)(e - (e > (int)d ? 1 : 0)) * MF_SPL + j] = (u8)(z_ ? z_ : 1u); }
}
}
if (s_sig) g_sig_merged = 1;
}
}
}
}
if (s_sig && !g_sig_merged && k==9) {
for(int d=0;d<k;d++) {
const u32 *od=s_ord[d];
u8 *base=s_sig+(size_t)d*k*MF_SPL;
u32 j=0;
switch(d) {
case 0: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 1: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 2: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 3: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 4: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 5: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 6: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),_mm_srli_si128(v3,8));
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 7: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),v4);
}
} break;
case 8: {
for(;j+8<=N;j+=8) {
if(j+64<N)_mm_prefetch((const char *)(a15b_sigrow+(size_t)od[j+64]*16),_MM_HINT_T0);
__m128i r0=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+0]*16));
__m128i r1=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+1]*16));
__m128i r2=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+2]*16));
__m128i r3=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+3]*16));
__m128i r4=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+4]*16));
__m128i r5=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+5]*16));
__m128i r6=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+6]*16));
__m128i r7=_mm_loadu_si128((const __m128i *)(a15b_sigrow+(size_t)od[j+7]*16));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
__m128i h0=_mm_unpackhi_epi8(r0,r1);
__m128i h1=_mm_unpackhi_epi8(r2,r3);
__m128i h2=_mm_unpackhi_epi8(r4,r5);
__m128i h3=_mm_unpackhi_epi8(r6,r7);
__m128i hu0=_mm_unpacklo_epi16(h0,h1),hu1=_mm_unpacklo_epi16(h2,h3);
__m128i v4=_mm_unpacklo_epi32(hu0,hu1);
_mm_storel_epi64((__m128i *)(base+(size_t)0*MF_SPL+j),v0);
_mm_storel_epi64((__m128i *)(base+(size_t)1*MF_SPL+j),_mm_srli_si128(v0,8));
_mm_storel_epi64((__m128i *)(base+(size_t)2*MF_SPL+j),v1);
_mm_storel_epi64((__m128i *)(base+(size_t)3*MF_SPL+j),_mm_srli_si128(v1,8));
_mm_storel_epi64((__m128i *)(base+(size_t)4*MF_SPL+j),v2);
_mm_storel_epi64((__m128i *)(base+(size_t)5*MF_SPL+j),_mm_srli_si128(v2,8));
_mm_storel_epi64((__m128i *)(base+(size_t)6*MF_SPL+j),v3);
_mm_storel_epi64((__m128i *)(base+(size_t)7*MF_SPL+j),_mm_srli_si128(v3,8));
}
} break;
}
for(;j<N;j++) {
const u8 *row=a15b_sigrow+(size_t)od[j]*16;
for(int e=0;e<k;e++)if(e!=d)base[(size_t)(e-(e>d?1:0))*MF_SPL+j]=row[e];
}
}
g_sig_merged=1;
}
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * k * MF_SPL;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB < N) {
if (k == 9) _mm_prefetch((const char *)(a15b_sigrow + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
else _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB] * STRIDE), _MM_HINT_T0);
}
#endif
if (k == 9) {
const u8 *row = a15b_sigrow + (size_t)od[j] * STRIDE;
for (int e = 0; e < 9; e++) if (e != d) base[(size_t)(e - (e > d ? 1 : 0)) * MF_SPL + j] = row[e];
} else {
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) if (e != d) { u32 z_ = row[e] >> g_sigsh; base[(size_t)(e - (e > d ? 1 : 0)) * MF_SPL + 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 = s_gt[ds][s_qord[qi]]; if (rr > mx) mx = rr; }
u32 wp = ((mx >> 6) + 8u) & ~7u;
Wpool = (wp < W) ? wp : W;
}
#endif
#if MF_PLEN
if (g_plen_on) {
u64 tot = 0;
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const u32 *mr = s_maxr + ((size_t)ds * MAXK + d) * MAXQ;
u32 *po = s_pofs + (size_t)e * MAXQ;
u32 *pl = s_plen + (size_t)e * MAXQ;
for (u32 b = 0; b <= Q; b++) {
u32 RR = mr[b];
u32 len = ((RR >> 6) + 8u) & ~7u;
if (len > Wpool) len = Wpool;
po[b] = (u32)tot; pl[b] = len; tot += len;
}
}
if (tot > (u64)K1 * (Q + 1) * W) { g_plen_on = 0; } // never; safety
}
#endif
{
const u32 *od = s_ord[ds];
{ const u32 B15Z_PD = 24;
for (u32 j = 0; j < N; j++) {
if (j + B15Z_PD < N) __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(s_nid + od[j + B15Z_PD])));
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;
{ const u32 B15Z_PDB = 12;
while (ptr < pend) {
if (ptr + B15Z_PDB < pend) __asm__ __volatile__("prefetcht0 %0" :: "m"(*(const char *)(s_nid + od[ptr + B15Z_PDB])));
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), LDA((const __m256i *)(cur + w)));
for (; w < wlen; w++) dst[w] = cur[w];
}
#else
{ u32 w = 0;
for (; w + 4 <= wlen; w += 4) _mm256_storeu_si256((__m256i *)(dst + w), LDU256((const __m256i *)(cur + w)));
for (; 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;
}
}
for (u32 i = 0; i < N; i++) out[i] = (u32)s_outlo[i];
}
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 4); }
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 5); }
void count_6d(int n, const unsigned *x[6], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 6); }
void count_7d(int n, const unsigned *x[7], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 7); }
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
void count_9d(int n, const unsigned *x[9], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 9); }
void count_10d(int n, const unsigned *x[10], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 10); }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 99.195 ms | 33 MB + 760 KB | Accepted | Score: 100 | 显示更多 |