// ===== REFERENCES（a14t_ 席，2026-10-02）=====
// [1] duck.ac 用户 saffah_codex_6a_agg3，提交 **#121092** <https://duck.ac/submission/121092>
// （**556.002216 ms** = 现役 T，2026-10-02 20:55:43）。本文件正文 = 该提交正文**逐字节** + 我方
// #121378 的 2 键 ASCR + 本席唯一改动（见下"思路"）。其实质保留的三处真改动：
// (a) `MF_QCOEF 36 -> 32`（桶数 203）；(b) **增量锚位图**（ASCR 趟内共享平面 (e1,b1) 随 b1
// 递增只 OR 新进点，省掉该维整个池 NT 拷贝）；(c) `fringe_dim` 的 `if constexpr` 维常量特化。
// [2] duck.ac 本账号 saffah_cc_v41_agg1，提交 **#121378** <https://duck.ac/submission/121378>
// （**554.427224 ms** = 本席开工时的 mine）= [1] 的正文 + a14u_ 席的 **2 键 ASCR**
// （d1 桶计数排序之前加一趟 d2 稳定计数排序 ⇒ 桶组内按 (b1,b2) 有序），本件逐字节保留。
// [3] 本账号 saffah_cc_v41_agg1，提交 #120919 / #120554 / #120336 / #114076 / #114041：在册六旗
// （`no-tree-reassoc` / `live-range-shrinkage` / `no-tree-forwprop` / `no-tree-ch` /
// `unroll-all-loops` / `no-sched-critical-path-heuristic`）与 `MAXQ=256` 下界的来源。
// [4] duck.ac 用户 saffah_codex_6s_agg2，提交 #105554 / #104600 / #108803 / #108836：正文中
// ASCR bucket-key cache、scatter load-chase 预取等刀的来源（正文内逐处署名，继承链见正文头部）。
// [5] duck.ac/problem/1014a。以上均为公开提交，正文按站点规则公开可见；各原提交无独立许可证声明，
// 按 /home/yjp/duck.ac/RULES.md 第 3 节署名引用账号与原提交地址。
// ======================
// ===== 思路 =====（a14t_ 席，2026-10-02）
// 【本发唯一变量 = **头部竞技场（`MF_HEADH`）+ 读长分类分趟**】正文 `solve_mf` 里新增一块紧凑
// `s_bsH` 竞技场：每个平面 (e,b) 的前 min(len,544) 个词按**自身长度紧密**拷贝到该区（总 ~6.2 MB，
// 平面间距由池的 len≈1600 词压到 544 词 ⇒ 读地址密度 ×2.9）。随后把每个 ds 趟的查询按**读长**分成
// 两段先后处理（`R < 544*64` 的短读查询在前），`g_ppool/g_ppofs` 令短读段只从竞技场取平面、长读段
// 照旧读原池。答案逐位不变（见闸门）。
// 【为什么 —— 本席判题机 rdtsc 单变量定价（`tools/probe.py`，同指令数、故意算错的消融件）】
// 基座 1991.9 Mcyc（与 a14u_ 的 1992.8/1992.7 复现一致，重复性 0.016%）：
// · **只把 pass A 的读地址指针整体压缩 2×**（写路径不动、读长度不动、指令数不动）= **1830.0 Mcyc（−8.13%）**；
// **压缩 4× = 1609.9 Mcyc（−19.2%）**；**读+写一起压缩 2× = 1834.9 Mcyc（−7.88%）** ⇒ 效应全在**读侧**，
// 与写体积无关。
// · **只按读长分两趟、地址不动 = −0.06%（零效应）** ⇒ 单纯换处理次序无用，必须真把读地址变密。
// ⇒ 结论：pass A 的瓶颈是**池的读地址足迹**（18 MB 池 ≫ 判题机 6 MB L3），不是 L1-miss 数、不是写体积。
// 【本件实测】判题机 **1976.6 Mcyc（−0.77%）**，ck 与基座逐字相同（352704819）。
// 【闸门】与 [2] 正文逐位对拍：**16 组全同**（n=3e5 × mode{10,11,13} × Q{自主,200,228,255} +
// n=2000 × mode{0..4} + n∈{2e4,5e4,1e5}），n≤2e4 组跑暴力 oracle **0 mismatch**；生产 checksum 351867771 不变。
// 【同族判负（本席判题机读数，勿重走）】把整池拆成 [头部 544 词][尾部 len-544 词] 两个紧致竞技场、
// 让长读查询走两段读（省掉竞技场拷贝）= **+2.5%**（分段读本身贵 ~3.3%，已由"分量法"定位：两件的
// 短读段完全相同，差别只在长读段）；池按 `(plane,b)` 词级交错/转置 = 单调放大行数，未试即否。
// 【目的】把 pass A 的读地址足迹压小，作为本席的最优实现交板面。
// ======================
#define MF_NOSDMP 1
#define MF_QCOEF 32
#pragma GCC optimize("O3,no-strict-aliasing","no-tree-reassoc","live-range-shrinkage","no-tree-forwprop","no-tree-ch","unroll-all-loops","no-sched-critical-path-heuristic","no-gcse")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=skylake")
#define MF_NOFILTB 1
#include <immintrin.h>
#include <type_traits>
#ifdef LOCAL_TEST
#include <cstdio>
#include <cstdlib>
#include <ctime>
#else
#ifdef __cplusplus
extern "C" void *malloc(unsigned long);
#else
void *malloc(unsigned long);
#endif
#endif
typedef unsigned u32;
typedef unsigned char u8;
typedef unsigned long long u64;
#define LDU256(p) ({ __m256i _v256; __asm__("vmovdqu %1, %0" : "=x"(_v256) : "m"(*(const __m256i *)(p))); _v256; })
#define LDU8U16(p) ({ __m256i _v; \
__asm__("vpmovzxbd %1, %0" : "=x"(_v) : "m"(*(const u64 *)(p))); _v; })
#ifndef MAXN
#define MAXN 1000005
#endif
#define MAXK 8
#define STRIDE 16
#ifndef MF_APF
#define MF_APF 0
#endif
#ifndef MF_QCOEF
#define MF_QCOEF 42
#endif
#ifndef MF_KSPAD
#define MF_KSPAD 1
#endif
#ifndef MF_FILTAUTO
#define MF_FILTAUTO 1
#endif
#ifndef MF_NT
#define MF_NT 1
#endif
#ifndef MF_FR_NOSORT
#define MF_FR_NOSORT 0
#endif
#ifndef MF_FR_NOCAND
#define MF_FR_NOCAND 0
#endif
#ifndef MF_BUILD_NOPOOL
#define MF_BUILD_NOPOOL 0
#endif
#ifndef MF_FR_NOREAD
#define MF_FR_NOREAD 0
#endif
#ifndef MF_NOSDMP
#define MF_NOSDMP 0
#endif
#ifndef MF_FB2
#define MF_FB2 0
#endif
#ifndef MF_FB3
#define MF_FB3 1
#endif
#ifndef MF_PLEN
#define MF_PLEN 1
#endif
#ifndef MF_BUDGETMB
#define MF_BUDGETMB 1900
#endif
#ifndef MF_NOFILTB
#define MF_NOFILTB 0
#endif
#define MF_PFB 32
#ifndef MF_PFB1
#ifndef MF_PFOR
#define MF_PFOR 192
#endif
#define MF_PFB1 128
#endif
#ifndef MF_PFB2
#define MF_PFB2 128
#endif
#ifndef MF_PFB3
#define MF_PFB3 4
#endif
#ifndef MF_PFB4
#define MF_PFB4 32
#endif
#ifndef MF_HEADH
#define MF_HEADH 544
#endif
#ifndef MF_2KEY
#define MF_2KEY 1
#endif
#ifndef MF_ASCR
#define MF_ASCR 1
#endif
#ifndef MF_WTRIM
#define MF_WTRIM 0
#endif
#define AS_WORDS 20480
#ifndef MF_PFA
#define MF_PFA 1
#endif
#ifndef MF_FR_NOMETA
#define MF_FR_NOMETA 0
#endif
#ifndef MAXQ
#define MAXQ 256
#endif
#define MF_SIGPL ((((size_t)MAXN + 63u) & ~(size_t)63u) + 0x840u)
static u32 s_ord[MAXK][MAXN];
static u32 s_gt[MAXK][MAXN];
static u32 s_pm[(size_t)MAXN * STRIDE] __attribute__((aligned(64)));
static u8 s_bd[(size_t)MAXN * MAXK];
static u8 s_bdv[(size_t)MAXN * MAXK];
static u32 s_g[MAXN];
static u32 s_qord[MAXN], s_tmp[MAXN], s_cntv[MAXN];
static u32 s_nid[MAXN];
static u32 s_pos[MAXK][MAXN];
static u32 s_posp[(size_t)MAXK * MAXN];
static u64 s_facc[(size_t)(MAXK - 1) * (((MAXN) + 63) / 64 + 16)];
static u32 s_bnde[(size_t)MAXK * (MAXQ + 2)];
static u64 *s_fbm = 0;
static u64 s_fbuf[8192] __attribute__((aligned(64)));
static __m256i khi8v, khi8hv;
static u32 g_Qc = 32;
static u32 g_fbm_W = 0;
static int g_use_filt = 0;
static int g_lexsort = 0;
static int g_nofilt = 0;
static u32 s_fpos[(MAXK + 1) * MAXQ];
static u32 s_bval[(MAXK + 1) * MAXQ];
static u64 *s_bs = 0;
static u64 *s_bsH = 0;
static const u64 *g_ppool = 0;
static const u32 *g_ppofs = 0;
static u32 s_pofsH[(size_t)MAXK * MAXQ];
static u32 s_plenH[(size_t)MAXK * MAXQ];
static u64 s_acc[MAXN / 64 + 256] __attribute__((aligned(64)));
static u32 *s_sdmp[MAXK];
static u32 *inc_anchor=0;static size_t inc_cap=0;static bool inc_on=false;
static u64 s_ascr[AS_WORDS] __attribute__((aligned(64)));
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;
alignas(32) static const unsigned char q8_bcast_mask[4][32] = {
{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4},
{8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8},
{12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12},
};
static u32 g_sigmul = 1;
static inline u8 sig_nz(u32 v) {u32 z=(u32)(((u64)v*g_sigmul)>>23u);return (u8)(z?z:1u);}
static u32 s_pofs[(size_t)MAXK * MAXQ];
static u32 s_plen[(size_t)MAXK * MAXQ];
static int g_plen_on = 0;
static u32 g_W, g_Q, g_B;
static u32 g_KS;
static int g_nosort = 0;
static int g_nosdmp = 0;
static int g_skip_and = 0, g_skip_fringe = 0;
#define FPOS(d, b) s_fpos[(size_t)(d) * MAXQ + (b)]
#define BVAL(d, b) s_bval[(size_t)(d) * MAXQ + (b)]
extern "C" int madvise(void *, unsigned long, int);
#define MF_MADV_HUGEPAGE 14
static void *pool_alloc(size_t bytes) {
void *raw = malloc(bytes + (2u << 20));
if (!raw) return 0;
void *p = (void *)(((size_t)raw + (2u << 20) - 1) & ~(size_t)((2u << 20) - 1));
madvise(p, bytes, MF_MADV_HUGEPAGE);
return p;
}
__attribute__((target("avx2,popcnt")))
static inline u32 hsum256(__m256i v) {
u64 b[4];
_mm256_storeu_si256((__m256i *)b, v);
return (u32)(b[0] + b[1] + b[2] + b[3]);
}
template<int KK>
__attribute__((target("avx2,popcnt")))
static u32 and_popcount_range(const u64 *const *bs, u32 Wr, u32 rem) {
const __m256i lm = _mm256_set1_epi8(0x0f);
const __m256i lk = _mm256_setr_epi8(0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);
const __m256i z = _mm256_setzero_si256();
__m256i tot = z, acc0 = z, acc1 = z;
u32 w = 0, cnt = 0, total = 0;
for (; w + 16 <= (Wr); w += 16) {
{ __asm__ __volatile__("nop;nop;nop;nop;nop;nop" :: "m"(*(const char(*)[32])bs[0]) : "ymm15"); }
_mm_prefetch((const char *)(bs[1] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[1] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[2] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[3] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[4] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[5] + w + 40 + 8), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 40), _MM_HINT_T0);
_mm_prefetch((const char *)(bs[6] + w + 40 + 8), _MM_HINT_T0);
__m256i a, b, ta, tb;
a = _mm256_and_si256(LDU256((bs[0] + w)), LDU256((bs[1] + w)));
b = _mm256_and_si256(LDU256((bs[0] + w + 4)), LDU256((bs[1] + w + 4)));
ta = _mm256_and_si256(LDU256((bs[2] + w)), LDU256((bs[3] + w)));
tb = _mm256_and_si256(LDU256((bs[2] + w + 4)), LDU256((bs[3] + w + 4)));
a = _mm256_and_si256(a, ta);
b = _mm256_and_si256(b, tb);
ta = _mm256_and_si256(LDU256((bs[4] + w)), LDU256((bs[5] + w)));
tb = _mm256_and_si256(LDU256((bs[4] + w + 4)), LDU256((bs[5] + w + 4)));
a = _mm256_and_si256(a, ta);
b = _mm256_and_si256(b, tb);
a = _mm256_and_si256(a, LDU256((bs[6] + w)));
b = _mm256_and_si256(b, LDU256((bs[6] + w + 4)));
__m256i c, d4, tc, td;
c = _mm256_and_si256(LDU256((bs[0] + w + 8)), LDU256((bs[1] + w + 8)));
d4 = _mm256_and_si256(LDU256((bs[0] + w + 12)), LDU256((bs[1] + w + 12)));
tc = _mm256_and_si256(LDU256((bs[2] + w + 8)), LDU256((bs[3] + w + 8)));
td = _mm256_and_si256(LDU256((bs[2] + w + 12)), LDU256((bs[3] + w + 12)));
c = _mm256_and_si256(c, tc);
d4 = _mm256_and_si256(d4, td);
tc = _mm256_and_si256(LDU256((bs[4] + w + 8)), LDU256((bs[5] + w + 8)));
td = _mm256_and_si256(LDU256((bs[4] + w + 12)), LDU256((bs[5] + w + 12)));
c = _mm256_and_si256(c, tc);
d4 = _mm256_and_si256(d4, td);
c = _mm256_and_si256(c, LDU256((bs[6] + w + 8)));
d4 = _mm256_and_si256(d4, LDU256((bs[6] + w + 12)));
if (KK - 1 > 7) {
for (int d = 8; d < KK; d++) {
a = _mm256_and_si256(a, LDU256((bs[d] + w)));
b = _mm256_and_si256(b, LDU256((bs[d] + w + 4)));
c = _mm256_and_si256(c, LDU256((bs[d] + w + 8)));
d4 = _mm256_and_si256(d4, LDU256((bs[d] + w + 12)));
}
}
{
__m256i lc = _mm256_and_si256(c, lm), hc = _mm256_and_si256(_mm256_srli_epi16(c, 4), lm);
__m256i ld = _mm256_and_si256(d4, lm), hd = _mm256_and_si256(_mm256_srli_epi16(d4, 4), lm);
__m256i la = _mm256_and_si256(a, lm), ha = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
__m256i lb = _mm256_and_si256(b, lm), hb = _mm256_and_si256(_mm256_srli_epi16(b, 4), lm);
__m256i x0 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, la), _mm256_shuffle_epi8(lk, ha));
__m256i x1 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lb), _mm256_shuffle_epi8(lk, hb));
__m256i x2 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, lc), _mm256_shuffle_epi8(lk, hc));
__m256i x3 = _mm256_add_epi8(_mm256_shuffle_epi8(lk, ld), _mm256_shuffle_epi8(lk, hd));
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(x0, x1));
acc1 = _mm256_add_epi8(acc1, _mm256_add_epi8(x2, x3));
if (++cnt == 15) {
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
acc0 = acc1 = z; cnt = 0;
}
}
}
for (; w + 4 <= Wr; w += 4) {
__m256i a = LDU256((bs[0] + w));
for (int d = 1; d < KK; d++) a = _mm256_and_si256(a, LDU256((bs[d] + w)));
__m256i lo = _mm256_and_si256(a, lm);
__m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), lm);
acc0 = _mm256_add_epi8(acc0, _mm256_add_epi8(_mm256_shuffle_epi8(lk, lo), _mm256_shuffle_epi8(lk, hi)));
}
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc0, z));
tot = _mm256_add_epi64(tot, _mm256_sad_epu8(acc1, z));
total = hsum256(tot);
for (; w < Wr; w++) {
u64 v = bs[0][w];
for (int d = 1; d < KK; d++) v &= bs[d][w];
total += (u32)__builtin_popcountll(v);
}
if (rem) {
u64 v = bs[0][Wr];
for (int d = 1; d < KK; d++) v &= bs[d][Wr];
total += (u32)__builtin_popcountll(v & ((1ull << rem) - 1));
}
return total;
}
template<int KK>
__attribute__((target("avx2,popcnt")))
__attribute__((flatten))
static void query_group(u32 N, u32 *out, u32 ds, u32 q0, u32 q1, u32 Q, u32 W) {
const u32 KS_C = (KK <= 8) ? 8u : (u32)KK;
const int k = KK;
const int K1 = KK - 1;
u32 *qo = s_qord, *tmp = s_tmp;
#if MF_ASCR
const u32 d1 = (ds == 0) ? 1u : 0u;
const u32 e1 = (d1 < ds) ? d1 : d1 - 1;
#if MF_2KEY
{
const u32 d2 = (ds <= 1) ? 2u : 1u;
u32 *cntv = s_cntv;
for (u32 v = 0; v <= Q; v++) cntv[v] = 0;
for (u32 qi = q0; qi < q1; qi++) { cntv[s_bd[(size_t)qo[qi] * MAXK + d2] + 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++) {
const u32 id = qo[qi];
tmp[cntv[s_bd[(size_t)id * MAXK + d2] + 1]++] = id;
}
u32 *sw = qo; qo = tmp; tmp = sw;
}
#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_PFB1 < q1) _mm_prefetch((const char *)(s_bd + (size_t)qo[qi + MF_PFB1] * MAXK + d1), _MM_HINT_T0);
#endif
const u32 id = qo[qi];
const u32 bb = s_bd[(size_t)id * MAXK + d1];
qo[qi] = id | (bb << 20);
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++) {
const u32 packed = qo[qi];
const u32 bb = packed >> 20;
const u32 id = packed & ((1u << 20) - 1u);
tmp[cntv[bb + 1]++] = id;
}
u32 *sw = qo; qo = tmp; tmp = sw;
}
#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)qo[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 = qo[qi]; tmp[cntv[s_bd[(size_t)id * MAXK + dd] + 1]++] = id; }
u32 *sw = qo; qo = tmp; tmp = sw;
}
}
#ifdef PF_DIST
for (u32 qi = q0; qi < q0 + PF_DIST && qi + PF_DIST < q1; qi++) {
u32 j = qo[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 used=0;if(inc_on)for(u32 w=0;w<((N+63u)>>6);w++)s_ascr[w]=0;
u32 qi = q0;
while (qi < q1) {
u32 b1 = s_bd[(size_t)qo[qi] * MAXK + d1];
u32 qe = qi + 1;
while (qe < q1 && s_bd[(size_t)qo[qe] * MAXK + d1] == b1) qe++;
const u64 *bp1 = g_plen_on ? (g_ppool + g_ppofs[(size_t)e1 * MAXQ + b1])
: (s_bs + ((size_t)e1 * (Q + 1) + b1) * W);
if(inc_on){u32 end=b1<Q?FPOS(d1,b1):(u32)N;while(used<end){u32 r=inc_anchor[used++];s_ascr[r>>6]|=1ull<<(r&63);}bp1=s_ascr;}else{
u32 mx = 1;
for (u32 t = qi; t < qe; t++) { u32 rr = s_gt[ds][qo[t]]; if (rr > mx) mx = rr; }
u32 nwd = (mx >> 6) + 1;
if (nwd <= AS_WORDS) {
const u64 *src = bp1;
for (u32 w = 0; w < nwd; w++) s_ascr[w] = src[w];
bp1 = s_ascr;
}
}
for (u32 qq = qi; qq < qe; qq++) {
u32 i = qo[qq];
const u8 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qq + MF_PFA;
if (qj < qe) {
u32 j = qo[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
const u32 R = s_gt[ds][i];
u32 cnt = 0;
u32 zmsk;
{
const __m256i bdv = LDU8U16(bd);
u32 zm = (u32)_mm256_movemask_ps(_mm256_castsi256_ps(
_mm256_cmpeq_epi32(bdv, _mm256_setzero_si256())));
zm &= ~(1u << (u32)ds);
if (k < 8) zm &= (1u << (u32)k) - 1u;
zmsk = zm;
}
if (R && !g_skip_and && !zmsk) {
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 ? (g_ppool + g_ppofs[(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 = qo[qi];
const u8 *bd = s_bd + (size_t)i * MAXK;
#if MF_PFA > 0
{
u32 qj = qi + MF_PFA;
if (qj < q1) {
u32 j = qo[qj];
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[ds] + j), _MM_HINT_T0);
}
}
#endif
#ifdef PF_DIST
if (qi + PF_DIST < q1) {
u32 j = qo[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 ? (g_ppool + g_ppofs[(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) { if (qo != s_qord) for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = qo[qi]; return; }
auto fringe_dim=[&](auto dim){ constexpr int d=decltype(dim)::value;
if (d == ds) return;
#if !MF_FR_NOSORT
#if MF_ASCR
if (d != (int)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)qo[qi + MF_PFB2] * MAXK + d), _MM_HINT_T0);
#endif
const u32 id = qo[qi];
const u32 bb = s_bd[(size_t)id * MAXK + d];
qo[qi] = id | (bb << 20);
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++) {
const u32 packed = qo[qi];
const u32 bb = packed >> 20;
const u32 id = packed & ((1u << 20) - 1u);
tmp[cntv[bb + 1]++] = id;
}
u32 *sw = qo; qo = tmp; tmp = sw;
}
#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];
const u8 *const pbase_sig = s_sig + (size_t)d * (k - 1) * MF_SIGPL;
const u32 *const sdmpd = s_sdmp[d];
for (u32 qi = q0; qi < q1; qi++) {
u32 i = qo[qi];
#if MF_PFB > 0
{
u32 qj = qi + MF_PFB3;
if (qj < q1) {
u32 j = qo[qj];
__asm__ __volatile__("prefetchw %0" ::"m"(*(const char *)(&out[j])));
_mm_prefetch((const char *)(s_pm + (size_t)j * STRIDE), _MM_HINT_T0);
_mm_prefetch((const char *)(s_gt[d] + j), _MM_HINT_T0);
_mm_prefetch((const char *)(s_bd + (size_t)j * MAXK + d), _MM_HINT_T0);
}
}
#endif
const u8 *bd = s_bd + (size_t)i * MAXK;
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 = (KK == 8) ? _mm256_load_si256((const __m256i *)rec)
: _mm256_or_si256(_mm256_load_si256((const __m256i *)rec), khi8v);
__m256i Avraw = _mm256_load_si256((const __m256i *)(rec + 8));
__m256i Avm = _mm256_or_si256(Avraw, notexlo);
__m256i qT = _mm256_min_epi32(qr, Avm);
__m256i qTh;
if (k > 8) {
__m256i qrh = _mm256_or_si256(LDU256((rec + 8)), khi8hv);
__m256i Avmh = _mm256_set1_epi32(0x7FFFFFFF);
qTh = _mm256_min_epi32(qrh, Avmh);
} else {
qTh = all8;
}
u32 add = 0;
const u32 *c0 = sdmpd ? (sdmpd + (size_t)lo * KS_C) : 0;
if (g_use_filt) {
const u32 w0 = lo >> 6;
const u32 nw = ((hi - 1) >> 6) - w0 + 1;
u64 *fb = s_fbuf;
const u64 wcs = (u64)N / g_Qc + 1;
const u32 FW = g_fbm_W;
int first = 1;
#if MF_FR_NOREAD
first = 0;
for (u32 w = 0; w < nw; w++) fb[w] = 0;
#endif
for (int e = 0; e < k; e++) {
if (e == d) continue;
u32 c = (u32)((u64)rec[e] / wcs) + 1;
if (c > g_Qc) c = g_Qc;
const u64 *row = s_fbm + ((size_t)(d * k + e) * (g_Qc + 1) + c) * FW + w0;
if (first) { for (u32 w = 0; w < nw; w++) fb[w] = row[w]; first = 0; }
else { for (u32 w = 0; w < nw; w++) fb[w] &= row[w]; }
}
fb[0] &= ~0ull << (lo & 63);
u32 r2 = hi & 63u;
if (r2) fb[nw - 1] &= (1ull << r2) - 1;
for (u32 w = 0; w < nw; w++) {
u64 vv = fb[w];
#if MF_FR_NOCAND
vv = 0;
#endif
while (vv) {
int t = __builtin_ctzll(vv);
vv &= vv - 1;
u32 j = ((w0 + w) << 6) + (u32)t;
const u32 *cp = c0 ? (c0 + (size_t)(j - lo) * KS_C) : (s_pm + (size_t)s_ord[d][j] * STRIDE);
__m256i pv = LDU256(cp);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
continue;
}
if (g_sig_on && (hi - lo) >= 2) {
__m256i qsv[16];
{
const __m256i mulv = _mm256_set1_epi32((int)g_sigmul);
const __m256i maxv = _mm256_set1_epi32(0x7FFFFFFF);
__m256i sv = _mm256_srli_epi32(_mm256_mullo_epi32(qT, mulv), 23);
if (k > 8) sv = _mm256_or_si256(sv, _mm256_cmpeq_epi32(qT, maxv));
if (k == 8) sv = _mm256_subs_epu8(sv, _mm256_set1_epi8(1));
const __m256i lohalf = _mm256_broadcastsi128_si256(_mm256_castsi256_si128(sv));
const __m256i hihalf = _mm256_broadcastsi128_si256(_mm256_extracti128_si256(sv, 1));
if (k == 8) {
switch (d) {
case 0:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 1:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 2:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 3:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 4:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 5:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 6:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
break;
case 7:
qsv[0] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[1] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[2] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
qsv[3] = _mm256_shuffle_epi8(lohalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[3]));
qsv[4] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[0]));
qsv[5] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[1]));
qsv[6] = _mm256_shuffle_epi8(hihalf, _mm256_load_si256((const __m256i *)q8_bcast_mask[2]));
break;
}
} else {
int t = 0;
for (int e = 0; e < k; e++) {
if (e == (int)d) continue;
const __m256i src = (e < 4) ? lohalf : hihalf;
qsv[t++] = _mm256_subs_epu8(
_mm256_shuffle_epi8(src, _mm256_load_si256((const __m256i *)q8_bcast_mask[e & 3])),
_mm256_set1_epi8(1));
}
}
}
u32 aj = lo, aacc = 0, aexact = 0;
for (; aj + 32 <= hi; aj += 32) {
__m256i ac0 = _mm256_setzero_si256(), ac1 = ac0;
for (int t = 0; t < KK - 1; t++) {
__m256i cv = LDU256((const u8 *)(pbase_sig + (size_t)t * MF_SIGPL + aj));
__m256i df = _mm256_subs_epu8(cv, qsv[t]);
if (t & 1) ac1 = _mm256_or_si256(ac1, df);
else ac0 = _mm256_or_si256(ac0, df);
}
__m256i ac = _mm256_or_si256(ac0, ac1);
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 = 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) {
const u32 base = aj, lk = 0u;
u32 kmask = (1u << r) - 1u;
__m256i ac0 = _mm256_setzero_si256(), ac1 = ac0;
for (int t = 0; t < KK - 1; t++) {
__m256i cv = LDU256((const u8 *)(pbase_sig + (size_t)t * MF_SIGPL + base));
__m256i df = _mm256_subs_epu8(cv, qsv[t]);
if (t & 1) ac1 = _mm256_or_si256(ac1, df);
else ac0 = _mm256_or_si256(ac0, df);
}
__m256i ac = _mm256_or_si256(ac0, ac1);
u32 decm=((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_setzero_si256()))) & kmask;
u32 amb=((u32)_mm256_movemask_epi8(_mm256_cmpeq_epi8(ac, _mm256_set1_epi8(1)))) & kmask;
aacc += (u32)__builtin_popcount(decm);
while (amb) {
int t = __builtin_ctz(amb); amb &= amb - 1;
const u32 *cp = c0 ? (c0 + (size_t)(base + (u32)t - lo) * KS_C) : (s_pm + (size_t)odd[base + (u32)t] * STRIDE);
__m256i pv = LDU256(cp);
u32 okl = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = LDU256((cp + 8));
okl &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
aacc += okl;
}
}
}
add = aacc;
} else if (c0) {
if (k <= 8) {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8);
c += 4 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__m256i pv = LDU256(c);
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
}
add = a0 + a1 + a2 + a3;
} else {
const u32 *c = c0;
u32 a0 = 0, a1 = 0, a2 = 0, a3 = 0;
u32 j = lo;
for (; j + 4 <= hi; j += 4) {
__m256i p0 = LDU256((c));
__m256i p1 = LDU256((c + KS_C));
__m256i p2 = LDU256((c + 2 * KS_C));
__m256i p3 = LDU256((c + 3 * KS_C));
__m256i h0 = LDU256((c + 8));
__m256i h1 = LDU256((c + KS_C + 8));
__m256i h2 = LDU256((c + 2 * KS_C + 8));
__m256i h3 = LDU256((c + 3 * KS_C + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p0), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h0), all8);
a1 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p1), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h1), all8);
a2 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p2), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h2), all8);
a3 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, p3), all8) & (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, h3), all8);
c += 4 * KS_C;
}
for (; j < hi; j++, c += KS_C) {
__m256i pv = LDU256(c);
__m256i pvh = LDU256((c + 8));
a0 += (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8) &
(u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add = a0 + a1 + a2 + a3;
}
} else {
for (u32 j = lo; j < hi; j++) {
const u32 *c = s_pm + (size_t)odd[j] * STRIDE;
__m256i pv = _mm256_load_si256((const __m256i *)c);
u32 ok = (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qT, pv), all8);
if (k > 8) {
__m256i pvh = _mm256_load_si256((const __m256i *)(c + 8));
ok &= (u32)_mm256_testc_si256(_mm256_cmpgt_epi32(qTh, pvh), all8);
}
add += ok;
}
}
out[i] += add;
}
};
if constexpr(KK>0) fringe_dim(std::integral_constant<int,0>{});
if constexpr(KK>1) fringe_dim(std::integral_constant<int,1>{});
if constexpr(KK>2) fringe_dim(std::integral_constant<int,2>{});
if constexpr(KK>3) fringe_dim(std::integral_constant<int,3>{});
if constexpr(KK>4) fringe_dim(std::integral_constant<int,4>{});
if constexpr(KK>5) fringe_dim(std::integral_constant<int,5>{});
if constexpr(KK>6) fringe_dim(std::integral_constant<int,6>{});
if constexpr(KK>7) fringe_dim(std::integral_constant<int,7>{});
if constexpr(KK>8) fringe_dim(std::integral_constant<int,8>{});
if constexpr(KK>9) fringe_dim(std::integral_constant<int,9>{});
if (qo != s_qord) for (u32 qi = q0; qi < q1; qi++) s_qord[qi] = qo[qi];
}
__attribute__((target("avx2")))
static void solve_mf(u32 N, const unsigned **x, u32 *out, int k) {
const u32 W = ((N + 63) / 64 + 7u) & ~7u;
g_W = W;
#ifdef MF_SKIP_AND
g_skip_and = 1;
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_Q
g_Q = MF_Q; g_B = (u32)((N + (MF_Q) - 1) / (MF_Q));
#endif
#ifdef MF_FILT
g_use_filt = 1;
#endif
#ifdef MF_NOFILT
g_use_filt = 0;
#endif
#ifdef MF_SKIP_FRINGE
g_skip_fringe = 1;
#endif
if (!g_Q) {
u64 r = (u64)N, s2 = 0;
while (s2 * s2 < r) s2++;
u64 overhead = (u64)(5 * k + 21) * N * 4;
overhead += (u64)k * k * 4 * N;
overhead += (u64)k * (k - 1) * (u64)MF_SIGPL;
overhead += (u64)N * ((k <= 8) ? 8u : (u32)k) * 4 * k;
u64 budget = ((u64)MF_BUDGETMB << 20);
if (budget > overhead + ((u64)64 << 20)) budget -= overhead; else budget = (u64)64 << 20;
#if MF_PLEN
static const u32 kfill[7] = {72, 67, 62, 58, 54, 50, 46};
u64 perfill = ((u64)(k - 1) * W * 8) * kfill[k - 4] / 100;
#else
u64 perfill = (u64)(k - 1) * W * 8;
#endif
u64 Qmem = budget / (perfill ? perfill : 1);
if (Qmem < 8) Qmem = 8;
static const u32 kscale[7] = {82, 91, 100, 108, 115, 122, 129};
u32 q = (u32)((s2 * MF_QCOEF + 99) / 100);
q = (u32)(((u64)q * kscale[k - 4] + 99) / 100);
if ((u64)q > Qmem) q = (u32)Qmem;
if (q > N) q = N;
g_Q = q;
g_B = (N + q - 1) / q;
}
if (g_Q > 255u) { g_Q = 255u; g_B = (N + g_Q - 1) / g_Q; }
const u32 Q = g_Q, B = g_B;
#if MF_PLEN
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];
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_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
{ u32 *gt = s_gt[d];
const u64 BMAG = (((u64)1 << 41) / B) + 1u;
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;
gt[i] = s_cntv[v];
#endif
u32 bb = (u32)((((u64)(s_cntv[v + 1] - 1u)) * (u64)BMAG) >> 41);
s_bdv[(size_t)d * N + i] = 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++) {
u8 *dst = s_bd + (size_t)i * MAXK;
u32 *p = s_pm + (size_t)i * STRIDE;
for (int d = 0; d < k; d++) { dst[d] = s_bdv[(size_t)d * N + i]; p[d] = x[d][i]; }
{
__m256i idx = _mm256_add_epi32(_mm256_cvtepu8_epi32(_mm_loadl_epi64((const __m128i *)dst)),
_mm256_setr_epi32(0, MAXQ, 2*MAXQ, 3*MAXQ,
4*MAXQ, 5*MAXQ, 6*MAXQ, 7*MAXQ));
_mm256_store_si256((__m256i *)(p + 8),
_mm256_i32gather_epi32((const int *)s_bval, idx, 4));
if (k < 8) for (int d = k; d < 8; d++) p[8 + d] = 0x7fffffffu;
}
}
{
u32 nn = N - 1; u32 sh = 0;
while (nn > 255u) { nn >>= 1; sh++; }
g_sigsh = sh;
{
u32 M = (u32)(((u64)255 << 23) / (u64)(N - 1)) + 1u;
if (M < 1u) M = 1u;
while (M > 1u && ((u32)(((u64)(N - 1)) * M >> 23)) > 255u) M--;
g_sigmul = M;
}
size_t need = (size_t)k * (k - 1) * (size_t)MF_SIGPL;
static u8 *sigpool = 0; static size_t sigcap = 0;
if (need > ((size_t)220 << 20)) { sigpool = 0; sigcap = 0; }
else if (need > sigcap) { sigpool = (u8 *)pool_alloc(need); sigcap = sigpool ? need : 0; }
if (sigpool) {
s_sig = sigpool;
g_sig_on = 1;
g_sig_merged = 0;
}
}
#if MF_NOFILTB
g_use_filt = 0;
#endif
#if MF_FB3
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
u32 ptr = 0;
u32 *bn = s_bnde + (size_t)e * (MAXQ + 2);
for (u32 c = 0; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (v && ptr < N && (u64)xe[pe[ptr]] < v) ptr++;
bn[c] = ptr;
}
}
}
#endif
if (g_use_filt) {
u32 Qc = g_Qc;
u64 wc = (u64)N / Qc + 1;
size_t need2 = (size_t)k * k * (Qc + 1) * W * sizeof(u64);
static u64 *fpool = 0; static size_t fcap = 0;
if (need2 > fcap) { fpool = (u64 *)pool_alloc(need2 + 4096); fcap = need2; }
if (fpool) {
s_fbm = fpool; g_fbm_W = W;
#if MF_FB2
for (int d = 0; d < k; d++)
for (int e = 0; e < k; e++) {
if (e == d) continue;
u64 *b0 = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) b0[w] = 0;
}
for (int e = 0; e < k; e++) {
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
int dmap[MAXK], km1 = 0;
for (int d = 0; d < k; d++) if (d != e) dmap[km1++] = d;
for (int i = 0; i < km1; i++) { u64 *a = s_facc + (size_t)i * W; for (u32 w = 0; w < W; w++) a[w] = 0; }
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
while (ptr < N && (u64)xe[pe[ptr]] < v) {
const u32 *pv = s_posp + (size_t)pe[ptr] * MAXK;
for (int i = 0; i < km1; i++) {
u32 q = pv[dmap[i]];
s_facc[(size_t)i * W + (q >> 6)] |= 1ull << (q & 63);
}
ptr++;
}
for (int i = 0; i < km1; i++) {
int d = dmap[i];
u64 *dst = s_fbm + ((size_t)(d * k + e) * (Qc + 1) + c) * W;
const u64 *src = s_facc + (size_t)i * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(src + w)));
for (; w < W; w++) dst[w] = src[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = src[w];
#endif
}
}
}
#else
for (int d = 0; d < k; d++) {
for (int e = 0; e < k; e++) {
if (e == d) continue;
const unsigned *xe = x[e];
const u32 *pe = s_ord[e];
const u32 *pd = s_pos[d];
u64 *base = s_fbm + (size_t)(d * k + e) * (Qc + 1) * W;
for (u32 w = 0; w < W; w++) base[w] = 0;
u64 *cur = s_acc;
for (u32 w = 0; w < W; w++) cur[w] = 0;
u32 ptr = 0;
for (u32 c = 1; c <= Qc; c++) {
#if MF_FB3
u32 pend = s_bnde[(size_t)e * (MAXQ + 2) + c];
#else
u64 v = (u64)c * wc;
if (v > 0xFFFFFFFFull) v = 0xFFFFFFFFull;
u32 pend = ptr;
{ const unsigned *xe_ = xe; const u32 *pe_ = pe;
while (pend < N && (u64)xe_[pe_[pend]] < v) pend++; }
#endif
while (ptr < pend) { u32 q = pd[pe[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = base + (size_t)c * W;
#if MF_NT
{ u32 w = 0;
for (; w + 4 <= W; w += 4) _mm256_stream_si256((__m256i *)(dst + w), _mm256_load_si256((const __m256i *)(cur + w)));
for (; w < W; w++) dst[w] = cur[w]; }
#else
for (u32 w = 0; w < W; w++) dst[w] = cur[w];
#endif
}
}
}
#endif
} else { g_use_filt = 0; }
_mm_sfence();
}
{
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 MF_NOSDMP
g_nosdmp = 1;
#endif
if (!g_nosdmp) {
u32 KS = k;
#if MF_KSPAD
if (k <= 8) KS = 8;
#endif
size_t per = (size_t)N * KS * 4;
size_t lim = (size_t)150 << 20;
if (k >= 9) lim = (size_t)470 << 20;
if (per * k <= lim) {
static u32 *pool = 0; static size_t cap = 0;
if (per * k > cap) { pool = (u32 *)pool_alloc(per * k + 256); cap = per * k; }
if (pool) {
g_KS = KS;
for (int d = 0; d < k; d++) {
s_sdmp[d] = pool + (size_t)d * N * g_KS;
const u32 *od = s_ord[d];
u32 *dst = s_sdmp[d];
const u32 KSx = g_KS;
if (KSx == (u32)k && (u32)k == (u32)STRIDE) {
u8 *sb = s_sig ? (s_sig + (size_t)d * (k - 1) * MF_SIGPL) : 0;
u32 j = 0;
if (sb && k == 8) {
const __m256i sigmul = _mm256_set1_epi32((int)g_sigmul);
const __m256i one = _mm256_set1_epi32(1);
for (; j + 8 <= N; j += 8) {
#if MF_PFB > 0
if (j + 0 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 0 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row0 = s_pm + (size_t)od[j + 0] * STRIDE;
__m256i rv0 = _mm256_loadu_si256((const __m256i *)row0);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 0) * KSx), rv0);
__m256i r0 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv0, sigmul), 23), one);
#if MF_PFB > 0
if (j + 1 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 1 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row1 = s_pm + (size_t)od[j + 1] * STRIDE;
__m256i rv1 = _mm256_loadu_si256((const __m256i *)row1);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 1) * KSx), rv1);
__m256i r1 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv1, sigmul), 23), one);
#if MF_PFB > 0
if (j + 2 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 2 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row2 = s_pm + (size_t)od[j + 2] * STRIDE;
__m256i rv2 = _mm256_loadu_si256((const __m256i *)row2);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 2) * KSx), rv2);
__m256i r2 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv2, sigmul), 23), one);
#if MF_PFB > 0
if (j + 3 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 3 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row3 = s_pm + (size_t)od[j + 3] * STRIDE;
__m256i rv3 = _mm256_loadu_si256((const __m256i *)row3);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 3) * KSx), rv3);
__m256i r3 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv3, sigmul), 23), one);
#if MF_PFB > 0
if (j + 4 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 4 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row4 = s_pm + (size_t)od[j + 4] * STRIDE;
__m256i rv4 = _mm256_loadu_si256((const __m256i *)row4);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 4) * KSx), rv4);
__m256i r4 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv4, sigmul), 23), one);
#if MF_PFB > 0
if (j + 5 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 5 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row5 = s_pm + (size_t)od[j + 5] * STRIDE;
__m256i rv5 = _mm256_loadu_si256((const __m256i *)row5);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 5) * KSx), rv5);
__m256i r5 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv5, sigmul), 23), one);
#if MF_PFB > 0
if (j + 6 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 6 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row6 = s_pm + (size_t)od[j + 6] * STRIDE;
__m256i rv6 = _mm256_loadu_si256((const __m256i *)row6);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 6) * KSx), rv6);
__m256i r6 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv6, sigmul), 23), one);
#if MF_PFB > 0
if (j + 7 + MF_PFB4 < N)
_mm_prefetch((const char *)(s_pm + (size_t)od[j + 7 + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row7 = s_pm + (size_t)od[j + 7] * STRIDE;
__m256i rv7 = _mm256_loadu_si256((const __m256i *)row7);
_mm256_stream_si256((__m256i *)(dst + (size_t)(j + 7) * KSx), rv7);
__m256i r7 = _mm256_max_epu32(
_mm256_srli_epi32(_mm256_mullo_epi32(rv7, sigmul), 23), one);
__m256i t0 = _mm256_unpacklo_epi32(r0, r1);
__m256i t1 = _mm256_unpackhi_epi32(r0, r1);
__m256i t2 = _mm256_unpacklo_epi32(r2, r3);
__m256i t3 = _mm256_unpackhi_epi32(r2, r3);
__m256i t4 = _mm256_unpacklo_epi32(r4, r5);
__m256i t5 = _mm256_unpackhi_epi32(r4, r5);
__m256i t6 = _mm256_unpacklo_epi32(r6, r7);
__m256i t7 = _mm256_unpackhi_epi32(r6, r7);
__m256i u0 = _mm256_unpacklo_epi64(t0, t2);
__m256i u1 = _mm256_unpackhi_epi64(t0, t2);
__m256i u2 = _mm256_unpacklo_epi64(t1, t3);
__m256i u3 = _mm256_unpackhi_epi64(t1, t3);
__m256i u4 = _mm256_unpacklo_epi64(t4, t6);
__m256i u5 = _mm256_unpackhi_epi64(t4, t6);
__m256i u6 = _mm256_unpacklo_epi64(t5, t7);
__m256i u7 = _mm256_unpackhi_epi64(t5, t7);
if ((int)d != 0) {
__m256i v = _mm256_permute2x128_si256(u0, u4, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((0 < (int)d) ? 0 : 0 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 1) {
__m256i v = _mm256_permute2x128_si256(u1, u5, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((1 < (int)d) ? 1 : 1 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 2) {
__m256i v = _mm256_permute2x128_si256(u2, u6, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((2 < (int)d) ? 2 : 2 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 3) {
__m256i v = _mm256_permute2x128_si256(u3, u7, 0x20);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((3 < (int)d) ? 3 : 3 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 4) {
__m256i v = _mm256_permute2x128_si256(u0, u4, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((4 < (int)d) ? 4 : 4 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 5) {
__m256i v = _mm256_permute2x128_si256(u1, u5, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((5 < (int)d) ? 5 : 5 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 6) {
__m256i v = _mm256_permute2x128_si256(u2, u6, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((6 < (int)d) ? 6 : 6 - 1) * MF_SIGPL + j,
&bytes, 8);
}
if ((int)d != 7) {
__m256i v = _mm256_permute2x128_si256(u3, u7, 0x31);
__m128i lo = _mm256_castsi256_si128(v);
__m128i hi = _mm256_extracti128_si256(v, 1);
__m128i p16 = _mm_packus_epi32(lo, hi);
__m128i p8 = _mm_packus_epi16(p16, p16);
u64 bytes = (u64)_mm_cvtsi128_si64(p8);
__builtin_memcpy(sb + (size_t)((7 < (int)d) ? 7 : 7 - 1) * MF_SIGPL + j,
&bytes, 8);
}
}
}
for (; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
__m256i rv = _mm256_loadu_si256((const __m256i *)row);
_mm256_stream_si256((__m256i *)dp, rv);
if (sb) {
__m256i sg = _mm256_srli_epi32(
_mm256_mullo_epi32(rv, _mm256_set1_epi32((int)g_sigmul)), 23);
sg = _mm256_max_epu32(sg, _mm256_set1_epi32(1));
u32 tmpk[MAXK];
_mm256_storeu_si256((__m256i *)tmpk, sg);
for (int e = 0; e < k; e++) { if (e != (int)d) sb[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + j] = (u8)tmpk[e]; }
}
}
} else
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
u32 *dp = dst + (size_t)j * KSx;
for (u32 e = 0; e < KSx; e++) dp[e] = (e < (u32)k) ? row[e] : 0u;
if (s_sig) {
u8 *sb2 = s_sig + (size_t)d * (k - 1) * MF_SIGPL;
for (int e = 0; e < k; e++) { if (e == (int)d) continue; sb2[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + j] = sig_nz(row[e]); }
}
}
if (s_sig) g_sig_merged = 1;
}
_mm_sfence();
}
}
}
if (s_sig && !g_sig_merged && k == 8) {
static u8 cached_signature[(size_t)MAXN * 8] __attribute__((aligned(64)));
for (u32 i = 0; i < N; i++)
for (int e = 0; e < 8; e++) {
u32 v = (u32)sig_nz(x[e][i]);
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-1)*MF_SIGPL;
u32 j=0;
for (;j+8<=N;j+=8) {
__m128i r0=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+0]*8));
__m128i r1=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+1]*8));
__m128i r2=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+2]*8));
__m128i r3=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+3]*8));
__m128i r4=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+4]*8));
__m128i r5=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+5]*8));
__m128i r6=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+6]*8));
__m128i r7=_mm_loadl_epi64((const __m128i *)(cached_signature+(size_t)od[j+7]*8));
__m128i t0=_mm_unpacklo_epi8(r0,r1);
__m128i t1=_mm_unpacklo_epi8(r2,r3);
__m128i t2=_mm_unpacklo_epi8(r4,r5);
__m128i t3=_mm_unpacklo_epi8(r6,r7);
__m128i u0=_mm_unpacklo_epi16(t0,t1);
__m128i u1=_mm_unpackhi_epi16(t0,t1);
__m128i u2=_mm_unpacklo_epi16(t2,t3);
__m128i u3=_mm_unpackhi_epi16(t2,t3);
__m128i v0=_mm_unpacklo_epi32(u0,u2);
__m128i v1=_mm_unpackhi_epi32(u0,u2);
__m128i v2=_mm_unpacklo_epi32(u1,u3);
__m128i v3=_mm_unpackhi_epi32(u1,u3);
if(d!=0) _mm_storel_epi64((__m128i *)(base+(size_t)(0 < d ? 0 : 0-1)*MF_SIGPL+j),v0);
if(d!=1) _mm_storel_epi64((__m128i *)(base+(size_t)(1 < d ? 1 : 1-1)*MF_SIGPL+j),_mm_srli_si128(v0,8));
if(d!=2) _mm_storel_epi64((__m128i *)(base+(size_t)(2 < d ? 2 : 2-1)*MF_SIGPL+j),v1);
if(d!=3) _mm_storel_epi64((__m128i *)(base+(size_t)(3 < d ? 3 : 3-1)*MF_SIGPL+j),_mm_srli_si128(v1,8));
if(d!=4) _mm_storel_epi64((__m128i *)(base+(size_t)(4 < d ? 4 : 4-1)*MF_SIGPL+j),v2);
if(d!=5) _mm_storel_epi64((__m128i *)(base+(size_t)(5 < d ? 5 : 5-1)*MF_SIGPL+j),_mm_srli_si128(v2,8));
if(d!=6) _mm_storel_epi64((__m128i *)(base+(size_t)(6 < d ? 6 : 6-1)*MF_SIGPL+j),v3);
if(d!=7) _mm_storel_epi64((__m128i *)(base+(size_t)(7 < d ? 7 : 7-1)*MF_SIGPL+j),_mm_srli_si128(v3,8));
}
for(;j<N;j++) {
const u8 *row=cached_signature+(size_t)od[j]*8;
for(int e=0;e<k;e++)if(e!=d)base[(size_t)(e<d?e:e-1)*MF_SIGPL+j]=row[e];
}
}
g_sig_merged=1;
}
if (s_sig && !g_sig_merged) {
for (int d = 0; d < k; d++) {
const u32 *od = s_ord[d];
u8 *base = s_sig + (size_t)d * (k - 1) * MF_SIGPL;
for (u32 j = 0; j < N; j++) {
#if MF_PFB > 0
if (j + MF_PFB4 < N) _mm_prefetch((const char *)(s_pm + (size_t)od[j + MF_PFB4] * STRIDE), _MM_HINT_T0);
#endif
const u32 *row = s_pm + (size_t)od[j] * STRIDE;
for (int e = 0; e < k; e++) { if (e == (int)d) continue; base[(size_t)((e < (int)d) ? e : e - 1) * MF_SIGPL + j] = sig_nz(row[e]); }
}
}
}
u64 poolwords = (u64)K1 * (Q + 1) * W;
#if MF_PLEN
if (g_plen_on) {
u64 best = 0, bestH = 0;
for (int dsx = 0; dsx < k; dsx++) {
u64 tot = 0, toth = 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;
toth += (len > MF_HEADH) ? (u32)MF_HEADH : len;
}
}
if (tot > best) best = tot;
if (toth > bestH) bestH = toth;
}
if (best && best < poolwords) poolwords = best;
{
static u64 *hpool = 0; static size_t hcap = 0;
size_t needH = (size_t)bestH * sizeof(u64) + 8192;
if (needH > hcap) { hpool = (u64 *)pool_alloc(needH); hcap = hpool ? needH : 0; }
s_bsH = hpool;
if (!s_bsH) g_plen_on = 0;
}
}
#endif
size_t need = (size_t)poolwords * sizeof(u64);
static u64 *pool = 0; static size_t cap2 = 0;
if (need > cap2) { pool = (u64 *)pool_alloc(need + 4096); cap2 = need; }
s_bs = pool;
{size_t bytes=(size_t)N*sizeof(u32)+4096;if(bytes>inc_cap){inc_anchor=(u32*)pool_alloc(bytes);inc_cap=inc_anchor?bytes:0;}inc_on=inc_anchor&&(u64)N<=(u64)AS_WORDS*64;}
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, toth = 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;
u32 *poH = s_pofsH + (size_t)e * MAXQ;
u32 *plH = s_plenH + (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;
{ u32 lh = (len > (u32)MF_HEADH) ? (u32)MF_HEADH : len;
poH[b] = (u32)toth; plH[b] = lh; toth += lh; }
}
}
if (tot > (u64)K1 * (Q + 1) * W) { g_plen_on = 0; }
}
#endif
const u32 *const spp = s_pos[ds];
for (int e = 0; e < K1; e++) {
int d = (e < ds) ? e : e + 1;
const unsigned *xd = x[d];
const u32 *od = s_ord[d];
if(inc_on&&d==((ds==0)?1:0)){for(u32 j=0;j<N;j++)inc_anchor[j]=spp[od[j]];continue;}
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 + 12 <= pend) {
u32 q = spp[od[ptr+0]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+1]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+2]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+3]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+4]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+5]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+6]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+7]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+8]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+9]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+10]]; cur[q >> 6] |= 1ull << (q & 63);
q = spp[od[ptr+11]]; cur[q >> 6] |= 1ull << (q & 63);
ptr += 12; }
while (ptr < pend) { u32 q = spp[od[ptr]]; cur[q >> 6] |= 1ull << (q & 63); ptr++; }
u64 *dst = g_plen_on ? (base + s_pofs[(size_t)e * MAXQ + b])
: (base + (size_t)b * Wpool);
u32 wlen = g_plen_on ? s_plen[(size_t)e * MAXQ + b] : Wpool;
if (g_plen_on && s_bsH) {
u32 lh = s_plenH[(size_t)e * MAXQ + b];
u64 *dh = s_bsH + s_pofsH[(size_t)e * MAXQ + b];
for (u32 w = 0; w < lh; w++) dh[w] = cur[w];
}
#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
}
{
const u32 R0 = (u32)MF_HEADH * 64u;
u32 lo = q0, hi = q1;
while (lo < hi) {
if (s_gt[ds][s_qord[lo]] < R0) lo++;
else { hi--; u32 t = s_qord[lo]; s_qord[lo] = s_qord[hi]; s_qord[hi] = t; }
}
if (s_bsH) {
g_ppool = s_bsH; g_ppofs = s_pofsH;
if (lo > q0) query_group<8>(N, out, ds, q0, lo, Q, Wpool);
g_ppool = s_bs; g_ppofs = s_pofs;
if (q1 > lo) query_group<8>(N, out, ds, lo, q1, Q, Wpool);
} else {
g_ppool = s_bs; g_ppofs = s_pofs;
query_group<8>(N, out, ds, q0, q1, Q, Wpool);
}
}
}
}
void count_8d(int n, const unsigned *x[8], unsigned *out) { solve_mf((u32)n, (const unsigned **)x, out, 8); }
#ifdef LOCAL_TEST
static void brute_nd(int n, const u32 **x, u32 *out, int k) {
for (int i = 0; i < n; i++) {
u32 c = 0;
for (int j = 0; j < n; j++) {
int ok = 1;
for (int d = 0; d < k; d++) if (!(x[d][j] < x[d][i])) { ok = 0; break; }
c += ok;
}
out[i] = c;
}
}
int main(int argc, char **argv) {
int k = (argc > 1) ? atoi(argv[1]) : 6;
int n = (argc > 2) ? atoi(argv[2]) : 200;
int seed = (argc > 3) ? atoi(argv[3]) : 1;
int mode = (argc > 4) ? atoi(argv[4]) : 0;
int qq = (argc > 5 && atoi(argv[5]) > 0) ? atoi(argv[5]) : 0;
g_Q = qq ? (u32)qq : (u32)(n / 250 + 1);
g_B = (u32)((n + g_Q - 1) / g_Q);
if (argc > 6) g_nosort = atoi(argv[6]);
if (argc > 7) g_nosdmp = atoi(argv[7]);
if (argc > 8) g_skip_and = atoi(argv[8]);
if (argc > 9) g_skip_fringe = atoi(argv[9]);
if (argc > 10) g_use_filt = atoi(argv[10]);
if (argc > 11) g_lexsort = atoi(argv[11]);
const int dm = (mode >= 10) ? mode - 10 : mode;
srand(seed);
static u32 *xs[10];
u32 *out = (u32 *)malloc(sizeof(u32) * n);
u32 *ref = (u32 *)malloc(sizeof(u32) * n);
for (int d = 0; d < k; d++) {
xs[d] = (u32 *)malloc(sizeof(u32) * n);
for (int i = 0; i < n; i++) {
if (dm == 0) xs[d][i] = (u32)(rand() % n);
else if (dm == 1) xs[d][i] = (u32)(rand() % 3);
else if (dm == 2) xs[d][i] = 0;
else if (dm == 3) xs[d][i] = (u32)(n - 1 - (rand() % n));
else if (dm == 4) xs[d][i] = (u32)(i);
else xs[d][i] = (u32)(rand() % n);
}
}
if (mode >= 10) {
struct timespec t0, t1;
clock_gettime(CLOCK_MONOTONIC, &t0);
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
clock_gettime(CLOCK_MONOTONIC, &t1);
double ms = (t1.tv_sec - t0.tv_sec) * 1e3 + (t1.tv_nsec - t0.tv_nsec) / 1e6;
u64 s = 0; for (int i = 0; i < n; i++) s += out[i];
printf("k=%d n=%d Q=%u TIME %.1f ms checksum %llu\n", k, n, g_Q, ms, s);
return 0;
}
switch (k) {
case 4: count_4d(n, (const unsigned **)xs, out); break;
case 5: count_5d(n, (const unsigned **)xs, out); break;
case 6: count_6d(n, (const unsigned **)xs, out); break;
case 7: count_7d(n, (const unsigned **)xs, out); break;
case 8: count_8d(n, (const unsigned **)xs, out); break;
case 9: count_9d(n, (const unsigned **)xs, out); break;
default: count_10d(n, (const unsigned **)xs, out); break;
}
brute_nd(n, (const u32 **)xs, ref, k);
int bad = 0;
for (int i = 0; i < n; i++) if (out[i] != ref[i]) { if (bad < 5) printf("MISMATCH i=%d got=%u ref=%u\n", i, out[i], ref[i]); bad++; }
printf("k=%d n=%d mode=%d Q=%u %s (%d mismatches)\n", k, n, mode, g_Q, bad ? "FAIL" : "OK", bad);
return bad != 0;
}
#endif
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 550.499 ms | 98 MB + 348 KB | Accepted | Score: 100 | 显示更多 |