/* ===== [f2xc_ 席] mmmf2k：**screen 相位的两把纯 ALU 刀**（无访存宽度/布局改动，逐位等价）
* 刀① `bmoments` 内环：3 级链式前缀和 (c0->c1->c2) 拆成**两条互不依赖的半流**（同一条
* 循环里 512 列上半 + 512 列下半各一组累加器），末尾用闭式合并恢复整行的 (c0,c1,c2)：
* c1 = c1a+c1b+M*c0a，c2 = c2a+c2b+M*c1a+M(M+1)/2*c0a（M=64 组）。链长减半、uop 不变、
* 访存形态一字未动。判题机 tryrun −165,004 拍（−0.0682%）。
* 刀② screen 的 `expected0/1` 互积环：`vpmulld`（8 个 int32 积/2 uop）换成 **`vpmaddwd`**
* （16 个 int16 积/1 uop）。做法：把 b 行矩心按 16 位拆成 v = 65536*h + l（l=(int16)v，
* h 与 l 都恰好落在 int16 内），a 侧用 **原始无符号字节**（+80 偏置折成常量 −80*sum_k r*[k]），
* 于是 Σ(a−80)·r = Σ a·l + 65536*Σ a·h − 80*Σ r。每个部分和都 < 2^31（逐项验证），
* 末尾在 int64 里合成再截断到 32 位 ⇒ 与原来的 int32 累加**逐位相同**。p01 压力减半，
* 实测本机该相位 −25%，判题机 tryrun −522,762 拍（−0.2160%）。
* 两刀合体（本发）判题机 tryrun **−660,976 拍 = −0.2731%**（同批两轮同号，chk 全等）。
* ===== REFERENCES =====
* [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #121727 <https://duck.ac/submission/121727>
* （判题机 Time = 71.356208 ms）= 正文出处；其后继承的全部注释/引用区原样保留在下方。
* [2] duck.ac 本账号 saffah_cc_v41_agg1：**#123937（70.800284 ms）= 本件基座**（前席 f2xb_：
* pack 环 bsum 的 p5 去化 —— 本席的两刀都**不碰**该处）· #123875（前席 f2zf_：mikro A 面 T1
* 预取）· #123395（前席 f2ze_）· #121816 · #123869。
* 合规：duck.ac 提交正文按站点公开可见、可直接取用；上述各件均无独立许可证声明。
* 按 RULES 第 3 节署名作者账号与原提交地址；本注释块之后保留 [1]/[2] 继承的全部注释区。
* ======================
* ======================
* ===== 思路 =====
* 【单变量组（本发 = 两把同族 ALU 刀，均不改访存宽度/存类型/布局）】
* ① `bmoments`（screen 的 b 行矩心）：`for(j=n-8;j>=0;j-=8){x=load;c0+=x;c1+=c0;c2+=c1;}` 的
* 3 级链每 8 列深 3 拍 ⇒ 整行 128 组 = 384 拍下限。改为 `h=((n>>1)&~7)`，同一条循环里
* 上半 (j=n-8..h) 与下半 (j=h-8..0) 各持一组 (c0,c1,c2)（6 个累加器 + 2 个载入 = 8 个
* ymm，寄存器够），互不依赖 ⇒ 每 8 列链深仍是 3，但两条链并行。
* 末尾**闭式合并**（纯整数结合律，无任何近似）：设上半有 M 组、下半有 Q 组，则
* 整行的 c0 = c0a+c0b；c1 = c1a+c1b+M*c0a；c2 = c2a+c2b+M*c1a+M(M+1)/2*c0a。
* （推导：下半的剩余步数恰为 Q-t = 它自己的局部倒计数；上半每组多计入 M。）
* 合并用 2 条 `vpmulld` 广播常量 + 6 条 `vpaddd`，每行 1 次，可忽略。
* **上界**：n=1024 时 |c2| ≤ 127*127*129/2*… 实测量级 ~2.9e7 << 2^31 ⇒ 无回绕；
* 即便回绕，int32 加法可结合可交换 ⇒ 与原名逐位相同。
* ② screen 的 `expected0/1`：原式为每行 1024 次的 int32 互积 `Σ_k (a[i][k]-80)*r0[k]`
* （`vpmulld` 在 SKX 上是 2 uop/p01、10 拍延迟）。改为 **int16 互积**：
* - 预拆 b 行矩心（每叶 1 次、1024 次迭代）：`int v=r0[k]; int lo=(short)(v&0xffff);
* r0l[k]=(int16)lo; r0h[k]=(int16)((v-lo)>>16);` ⇒ v = 65536*h + l，**h、l 都在
* int16 内**（|r0| ≤ 1024*80 = 81920 ⇒ h ∈ [-2,1]；|r1| ≤ 41.9e6 ⇒ h ∈ [-640,640]）。
* - 主环：`_mm256_cvtepu8_epi16` 一次吃 16 个 a 字节（**用原始无符号字节，不减去 80**），
* 用 Σ(a−80)·r = Σ a·l + 65536·Σ a·h − 80·(Σ_k r) 把偏置折成一个**常量**；
* 每个累加器每 16 列只吃 1 条 `vpmaddwd`（16 个积）+ 1 条 `vpaddd`。
* - **溢出闸门（逐项）**：vpmaddwd 的两两积和 ≤ 2*255*32767 = 1.67e7 < 2^31 ✓；
* 累加器 64 次 ⇒ ≤ 1.07e9 < 2^31 ✓（h 支路更小）；**全链路无回绕**。末尾
* `(uint32)((int64)(uint32)hsum32(x0) + 65536*(int64)(uint32)hsum32(x1) − 80*(int64)RS0)`
* —— 32 位模运算下与旧式 int32 累加完全同余 ⇒ 逐位相同。
* - 寄存器：**一行一趟、4 个累加器**（先试过"两行一趟 8 累加器"⇒ GCC 把 4 个累加器
* spill 到栈上（`vpaddd 160(%rsp)`），本机只值 −14%；改成一行一趟后 −25%）。
* 【判题机判据（`tools/tryrun.py`，同批 4 臂轮转 x2 轮、warm rdtsc 取 min、chk 全等 854faeacca6e9c9c）】
* 基座(#123937 正文) 242,009,446 / 242,027,790 → ① 刀 241,844,442（−0.0682%）→
* ② 刀 241,486,684（−0.2160%）→ **①+② 241,348,470 / 241,392,128（−660,976 拍 = −0.2731%）**。
* 【本机逐位闸门（4 种子 x 2 分布，含 ±1 端点与 repair 计数）】
* 854faeacca6e9c9c/175 · 37a50e8457f6257a/15360 · efa28c3c62f2153b/100 · 90f92986f1dd0380/15360 ·
* fc5cc61502503ce6/205 · f697ce823c9fd021/15360 · 28803d4dbba50635/145 · b38c6d610629e798/15360
* —— 与基座 #123937 逐项全等（本席另核：`bmoments` 的 r0/r1/r2 三数组、panel、bsum 逐字节不变）。
* 【与已判死族的关系】不改 mikro/combine/pack 一行；不改任何访存宽度或存类型（本席只动 ALU 与
* 两处栈上数组的**布局无关**拆分）；不新增循环外访存；不引入 AVX-512；提交件不含任何 main 函数（本席已核：提交件只有库函数与 `matrix_multiply`，无驱动程序）。
* ================
*/
/* ===== [f2xb_ 席] mmmf2k：**pack 环 bsum 的 p5 去化** —— 列和改在 int16 域累加
* 32 个 k-step 后才加宽（`vpmovsxwd` 2 条/k-step → 0.125 条/k-step），pack 环 p5 40 → 32。
* **bp 逐字节不变**（转置环一行未动）；**mikro 环（含 #123875 的 A 面 T1 预取）一行未动**。
* ===== REFERENCES =====
* [1] duck.ac 用户 saffah_codex_6a_agg3，提交 #121727 <https://duck.ac/submission/121727>
* （判题机 Time = 71.356208 ms）= 正文出处；其后继承的全部注释/引用区原样保留在下方。
* [2] duck.ac 本账号 saffah_cc_v41_agg1：**#123875（70.835434 ms）= 本件基座**（前席 f2zf_：
* mikro A 面「下一个 ic 块」T1 预取）· #123395（前席 f2ze_：pack 4 条 T0 预取）· #121816
* · #123869 · 以及本席前两发 #123935 / #123936（见下「两发定位」）。
* 合规：duck.ac 提交正文按站点公开可见、可直接取用；上述各件均无独立许可证声明。
* 按 RULES 第 3 节署名作者账号与原提交地址；本注释块之后保留 [1]/[2] 继承的全部注释区。
* ======================
* ======================
* ===== 思路 ===== （本席第三发 = 前两发的去污染重做；唯一变量 = bsum 的 p5 去化）
* 【★★★ 前两发为什么判负：基座文件取错了（全队可引的过程教训）】
* 本席先前拿 `work/f2zf__tmp/f2zf_a0.cpp` 当「#123875 正文」用 —— 实测它其实是 **#123395 的正文**
* （**缺 f2zf_ 那条 A 面预取**，逐行 diff 只有 3 处差别：缺 1 条 `_MM_HINT_T1` 预取、多 2 处
* `f2ze_rep` 计数插桩）。⇒ #123935（71.190757）/ #123936（71.168929）两发的板面读数**都落在
* #123395 的 71.190138 附近**（±21 µs），而不是落在基座 #123875 的 70.835434 附近 ——
* **两发并没有「反号」，它们只是少了那条 −0.355 ms 的预取**。本席已用 `tools/tryrun.py` 复核：
* 取真正的 #123875 正文（`f2zf_sub3.cpp` 去头）作锚，读数 **242,212,668 拍**，与 f2zf_ 记录在案的
* #123875 读数 242,220,974 一致 ✓（而错基座在判题机上读 243,561,082 = #123395 的档位）。
* **教训：臂文件要与「当前最好件的发放件」逐行 diff 校验，不能信族内先辈留下的 *_a0.cpp。**
* 【本发改动（唯一）】`leaf()` 的 pack 转置环里，每个 16 列组的列和 `wa`/`wb` 不再每 k-step 就
* `vpmovsxwd` 加宽（2 条 p5 + 2 条 vpaddd），而是用 `vpaddw`（p015）攒在 int16 域，每 32 个
* k-step 才 flush 一次。**上界 32 x 4 行 x 128 = 16384 < 32767 ⇒ 不可能溢出**；整数加法可结合
* ⇒ 结果逐位相同；bp 与 demás 面板逐字节不变 ⇒ mikro/combine/repair 全未受影响。
* 【判题机判据（`tools/tryrun.py`，同批两臂轮转 x3 轮、warm rdtsc 取 min、chk 全等）】
* 真基座(#123875 正文) 242,212,668 / 242,274,478 / 242,300,666 → **本件 241,967,820 / 242,002,104 /
* 242,040,830 = −244,848 拍 = −0.1011%**（3 轮同号）。**这才是只含 bsum 那一个变量的读数。**
* 【传导率（本席三发合起来才看得出的数）】同一族里只有 **ALU/uop** 的子族有部分传导：
* bsum 刀在板面（去污染前的那一发）实测 −21.2 µs ↔ tryrun −87 µs ⇒ **传导 ≈ 34%**；
* 而「访存宽度/存取型」子族（16B 存→32B 存、WB→NT）两发都是 **0 传导**（NT 臂 tryrun +0.435%、
* 板面同号；ymm/32B 存那一半在板面把 bsum 的收益吃掉）。本发只保留 **纯 ALU** 的那一半。
* 【闸门（本机）】4 种子 x 2 分布（均匀[-1,1] / ±1 端点）C 的 FNV 与 repair 计数逐项全等：
* 854faeacca6e9c9c/175 · 37a50e8457f6257a/15360 · efa28c3c62f2153b/100 · 90f92986f1dd0380/15360 ·
* fc5cc61502503ce6/205 · f697ce823c9fd021/15360 · f05aaf0406da4b2d/115 · a10f7dae53d489cb/15360。
* 本件不改任何量化/偏置 ⇒ 误差与基座恒等（无需 f2ze_err.cpp）。
* ================
*/
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing,web,ira-loop-pressure")
#pragma GCC target("avx2,fma")
#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>
#include <cstring>
static constexpr int MR=4, MC=208, NR=16, QA=40, QB=40;
static int8_t *aq,*bq;
alignas(4096) static unsigned char workspace[28*1024*1024];
alignas(64) static uint8_t bp[1024*1024]; /* m2k_: the A panel `ap` is gone */
alignas(64) static int bsum[1024];
struct Mat { const int8_t *x; int ld; const int8_t *y=nullptr; int ly=0; int sign=0; };
// 取 16 个 int8；x±y 都在 int8 域完成（各源 |v|<=46 => 和 |v|<=92 <=127，不溢出）
static inline __m128i ld(Mat a,int i,int j) {
__m128i x=_mm_loadu_si128((const __m128i*)(a.x+i*a.ld+j));
if(a.sign) {__m128i y=_mm_loadu_si128((const __m128i*)(a.y+i*a.ly+j));
x=a.sign>0?_mm_add_epi8(x,y):_mm_sub_epi8(x,y);}
return x;
}
static inline __m128i load4(Mat a,int i,int k) {
__m128i x=_mm_loadl_epi64((const __m128i*)(a.x+i*a.ld+k));
if(a.sign) {
__m128i y=_mm_loadl_epi64((const __m128i*)(a.y+i*a.ly+k));
x=a.sign>0?_mm_add_epi16(x,y):_mm_sub_epi16(x,y);
}
return x;
}
static __attribute__((noinline)) void kernel(const uint8_t *a,const uint8_t *b,uint32_t *out,const int *sums) {
ptrdiff_t stride=28672;
int k=256;
asm volatile (
"vpxor %%ymm14, %%ymm14, %%ymm14\n\t"
"vmovdqa (%[sums]), %%ymm11\n\t"
"vmovdqa 32(%[sums]), %%ymm12\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm0\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm1\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm2\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm3\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm4\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm5\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm6\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm7\n\t"
"vpcmpeqw %%ymm10, %%ymm10, %%ymm10\n\t"
"vpsrlw $15, %%ymm10, %%ymm10\n\t"
".p2align 5\n\t"
"1:\n\t"
"vmovdqa 0(%[b]), %%ymm11\n\t"
"vmovdqa 32(%[b]), %%ymm12\n\t"
"vpbroadcastd 0(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 1024(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 2048(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 3072(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 64(%[b]), %%ymm11\n\t"
"vmovdqa 96(%[b]), %%ymm12\n\t"
"vpbroadcastd 4(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 1028(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 2052(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 3076(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 128(%[b]), %%ymm11\n\t"
"vmovdqa 160(%[b]), %%ymm12\n\t"
"vpbroadcastd 8(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 1032(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 2056(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 3080(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 192(%[b]), %%ymm11\n\t"
"vmovdqa 224(%[b]), %%ymm12\n\t"
"vpbroadcastd 12(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 1036(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 2060(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm8\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm9\n\t"
"vpbroadcastd 3084(%[a]), %%ymm13\n\t"
"vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
"vpmaddubsw %%ymm12, %%ymm13, %%ymm15\n\t"
"vpmaddwd %%ymm10, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"add $16, %[a]\n\t"
"add $256, %[b]\n\t"
"sub $4, %[k]\n\t"
"jnz 1b\n\t"
"vmovntdq %%ymm0, (%[out])\n\t"
"vmovntdq %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntdq %%ymm2, (%[out])\n\t"
"vmovntdq %%ymm3, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntdq %%ymm4, (%[out])\n\t"
"vmovntdq %%ymm5, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntdq %%ymm6, (%[out])\n\t"
"vmovntdq %%ymm7, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [sums] "r"(sums)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
alignas(64) static uint32_t observed_moments[1024][16];
alignas(64) static int moment_weights[1024];
static __attribute__((noinline)) void kernel_pair2(const uint8_t *a,const uint8_t *b,uint32_t *out,const int *sums,uint32_t*mom,const int*weight) {
ptrdiff_t stride=28672;
int k=256;
asm volatile (
"vpxor %%ymm14, %%ymm14, %%ymm14\n\t"
"vmovdqa (%[sums]), %%ymm11\n\t"
"vmovdqa 32(%[sums]), %%ymm12\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm0\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm1\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm2\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm3\n\t"
"vpcmpeqw %%ymm15, %%ymm15, %%ymm15\n\t"
"vpsrlw $15, %%ymm15, %%ymm15\n\t"
".p2align 5\n\t"
"1:\n\t"
"vpbroadcastd 0(%[a]), %%ymm10\n\t"
"vpmaddubsw 0(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 32(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 4(%[a]), %%ymm10\n\t"
"vpmaddubsw 64(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1024(%[a]), %%ymm10\n\t"
"vpmaddubsw 0(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 32(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1028(%[a]), %%ymm10\n\t"
"vpmaddubsw 64(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 8(%[a]), %%ymm10\n\t"
"vpmaddubsw 128(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 160(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 12(%[a]), %%ymm10\n\t"
"vpmaddubsw 192(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1032(%[a]), %%ymm10\n\t"
"vpmaddubsw 128(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 160(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1036(%[a]), %%ymm10\n\t"
"vpmaddubsw 192(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 16(%[a]), %%ymm10\n\t"
"vpmaddubsw 256(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 288(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 20(%[a]), %%ymm10\n\t"
"vpmaddubsw 320(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1040(%[a]), %%ymm10\n\t"
"vpmaddubsw 256(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 288(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1044(%[a]), %%ymm10\n\t"
"vpmaddubsw 320(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 24(%[a]), %%ymm10\n\t"
"vpmaddubsw 384(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 416(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 28(%[a]), %%ymm10\n\t"
"vpmaddubsw 448(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1048(%[a]), %%ymm10\n\t"
"vpmaddubsw 384(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 416(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1052(%[a]), %%ymm10\n\t"
"vpmaddubsw 448(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 32(%[a]), %%ymm10\n\t"
"vpmaddubsw 512(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 544(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 36(%[a]), %%ymm10\n\t"
"vpmaddubsw 576(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1056(%[a]), %%ymm10\n\t"
"vpmaddubsw 512(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 544(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1060(%[a]), %%ymm10\n\t"
"vpmaddubsw 576(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 40(%[a]), %%ymm10\n\t"
"vpmaddubsw 640(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 672(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 44(%[a]), %%ymm10\n\t"
"vpmaddubsw 704(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1064(%[a]), %%ymm10\n\t"
"vpmaddubsw 640(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 672(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1068(%[a]), %%ymm10\n\t"
"vpmaddubsw 704(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 48(%[a]), %%ymm10\n\t"
"vpmaddubsw 768(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 800(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 52(%[a]), %%ymm10\n\t"
"vpmaddubsw 832(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1072(%[a]), %%ymm10\n\t"
"vpmaddubsw 768(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 800(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1076(%[a]), %%ymm10\n\t"
"vpmaddubsw 832(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 56(%[a]), %%ymm10\n\t"
"vpmaddubsw 896(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 928(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 60(%[a]), %%ymm10\n\t"
"vpmaddubsw 960(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1080(%[a]), %%ymm10\n\t"
"vpmaddubsw 896(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 928(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1084(%[a]), %%ymm10\n\t"
"vpmaddubsw 960(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm5\n\t"
"vpaddw %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"add $64, %[a]\n\t" "add $1024, %[b]\n\t" "sub $16, %[k]\n\t" "jnz 1b\n\t"
"vmovntdq %%ymm0, (%[out])\n\t"
"vmovntdq %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntdq %%ymm2, (%[out])\n\t"
"vmovntdq %%ymm3, 32(%[out])\n\t"
"test %[mom], %[mom]\n\t" "jz 3f\n\t"
"vmovdqa (%[weight]), %%ymm4\n\t"
"vmovdqa 32(%[weight]), %%ymm5\n\t"
"vpaddd %%ymm0, %%ymm1, %%ymm6\n\t"
"vpaddd 0(%[mom]), %%ymm6, %%ymm6\n\t"
"vmovdqa %%ymm6, 0(%[mom])\n\t"
"vpmulld %%ymm4, %%ymm0, %%ymm0\n\t"
"vpmulld %%ymm5, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm0, %%ymm1, %%ymm6\n\t"
"vpaddd 32(%[mom]), %%ymm6, %%ymm6\n\t"
"vmovdqa %%ymm6, 32(%[mom])\n\t"
"vpaddd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpaddd 64(%[mom]), %%ymm6, %%ymm6\n\t"
"vmovdqa %%ymm6, 64(%[mom])\n\t"
"vpmulld %%ymm4, %%ymm2, %%ymm2\n\t"
"vpmulld %%ymm5, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpaddd 96(%[mom]), %%ymm6, %%ymm6\n\t"
"vmovdqa %%ymm6, 96(%[mom])\n\t"
"3:\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [sums] "r"(sums), [mom] "r"(mom), [weight] "r"(weight)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
// Checksums of AB with weights 1,j,j^2; these are exact integer identities.
/* [f2mx] Exact row moments of a b panel: r0[k]=sum_j b[k][j], r1[k]=sum_j j*b[k][j],
r2[k]=sum_j j*j*b[k][j]. These are the closed forms of the scalar triple recurrence
the previous version ran (x=sum b, y=sum (j+1)b, z=sum (j+1)(j+2)/2 b => r0=x,
r1=y-x, r2=2z-3y+x), so every value -- and therefore every repair decision -- is
bit-identical; only the instruction count changes (that loop was 68% of the repair's
instructions: 14.16M of 20.85M, measured with callgrind on the copied file).
32 columns per iteration: the local index t (0..31) is handled with FIXED int16
weights via vpmaddwd, and the block number q (j = 32q + t) enters polynomially
(j*j = 1024q*q + 64q*t + t*t), so no per-block weight vector is built. The two
16-column halves of the block reuse one another: for t = 16+u,
madd(hi,w1) = sum u*b, and madd(hi,w1) + 16*madd(hi,one) = sum t*b. */
static __attribute__((noinline)) void kernel_pair2s(const uint8_t *a,const uint8_t *b,uint32_t *out,const int *sums) {
ptrdiff_t stride=28672;
int k=256;
asm volatile (
"vpxor %%ymm14, %%ymm14, %%ymm14\n\t"
"vmovdqa (%[sums]), %%ymm11\n\t"
"vmovdqa 32(%[sums]), %%ymm12\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm0\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm1\n\t"
"vpsubd %%ymm11, %%ymm14, %%ymm2\n\t"
"vpsubd %%ymm12, %%ymm14, %%ymm3\n\t"
"vpcmpeqw %%ymm15, %%ymm15, %%ymm15\n\t"
"vpsrlw $15, %%ymm15, %%ymm15\n\t"
".p2align 5\n\t"
"1:\n\t"
"vpbroadcastd 0(%[a]), %%ymm10\n\t"
"vpmaddubsw 0(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 32(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 4(%[a]), %%ymm10\n\t"
"vpmaddubsw 64(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1024(%[a]), %%ymm10\n\t"
"vpmaddubsw 0(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 32(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1028(%[a]), %%ymm10\n\t"
"vpmaddubsw 64(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 8(%[a]), %%ymm10\n\t"
"vpmaddubsw 128(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 160(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 12(%[a]), %%ymm10\n\t"
"vpmaddubsw 192(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1032(%[a]), %%ymm10\n\t"
"vpmaddubsw 128(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 160(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1036(%[a]), %%ymm10\n\t"
"vpmaddubsw 192(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 16(%[a]), %%ymm10\n\t"
"vpmaddubsw 256(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 288(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 20(%[a]), %%ymm10\n\t"
"vpmaddubsw 320(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1040(%[a]), %%ymm10\n\t"
"vpmaddubsw 256(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 288(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1044(%[a]), %%ymm10\n\t"
"vpmaddubsw 320(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 24(%[a]), %%ymm10\n\t"
"vpmaddubsw 384(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 416(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 28(%[a]), %%ymm10\n\t"
"vpmaddubsw 448(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1048(%[a]), %%ymm10\n\t"
"vpmaddubsw 384(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 416(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1052(%[a]), %%ymm10\n\t"
"vpmaddubsw 448(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 32(%[a]), %%ymm10\n\t"
"vpmaddubsw 512(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 544(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 36(%[a]), %%ymm10\n\t"
"vpmaddubsw 576(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1056(%[a]), %%ymm10\n\t"
"vpmaddubsw 512(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 544(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1060(%[a]), %%ymm10\n\t"
"vpmaddubsw 576(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 40(%[a]), %%ymm10\n\t"
"vpmaddubsw 640(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 672(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 44(%[a]), %%ymm10\n\t"
"vpmaddubsw 704(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1064(%[a]), %%ymm10\n\t"
"vpmaddubsw 640(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 672(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1068(%[a]), %%ymm10\n\t"
"vpmaddubsw 704(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 48(%[a]), %%ymm10\n\t"
"vpmaddubsw 768(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 800(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 52(%[a]), %%ymm10\n\t"
"vpmaddubsw 832(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1072(%[a]), %%ymm10\n\t"
"vpmaddubsw 768(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 800(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1076(%[a]), %%ymm10\n\t"
"vpmaddubsw 832(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"vpbroadcastd 56(%[a]), %%ymm10\n\t"
"vpmaddubsw 896(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 928(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 60(%[a]), %%ymm10\n\t"
"vpmaddubsw 960(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm4, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm5, %%ymm1, %%ymm1\n\t"
"vpbroadcastd 1080(%[a]), %%ymm10\n\t"
"vpmaddubsw 896(%[b]), %%ymm10, %%ymm8\n\t"
"vpmaddubsw 928(%[b]), %%ymm10, %%ymm9\n\t"
"vpbroadcastd 1084(%[a]), %%ymm10\n\t"
"vpmaddubsw 960(%[b]), %%ymm10, %%ymm4\n\t"
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm5\n\t"
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpmaddwd %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmaddwd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm8, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm4, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm5, %%ymm3, %%ymm3\n\t"
"add $64, %[a]\n\t" "add $1024, %[b]\n\t" "sub $16, %[k]\n\t" "jnz 1b\n\t"
"vmovntdq %%ymm0, (%[out])\n\t"
"vmovntdq %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntdq %%ymm2, (%[out])\n\t"
"vmovntdq %%ymm3, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [sums] "r"(sums)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
// Checksums of AB with weights 1,j,j^2; these are exact integer identities.
/* [f2mx] Exact row moments of a b panel: r0[k]=sum_j b[k][j], r1[k]=sum_j j*b[k][j],
r2[k]=sum_j j*j*b[k][j]. These are the closed forms of the scalar triple recurrence
the previous version ran (x=sum b, y=sum (j+1)b, z=sum (j+1)(j+2)/2 b => r0=x,
r1=y-x, r2=2z-3y+x), so every value -- and therefore every repair decision -- is
bit-identical; only the instruction count changes (that loop was 68% of the repair's
instructions: 14.16M of 20.85M, measured with callgrind on the copied file).
32 columns per iteration: the local index t (0..31) is handled with FIXED int16
weights via vpmaddwd, and the block number q (j = 32q + t) enters polynomially
(j*j = 1024q*q + 64q*t + t*t), so no per-block weight vector is built. The two
16-column halves of the block reuse one another: for t = 16+u,
madd(hi,w1) = sum u*b, and madd(hi,w1) + 16*madd(hi,one) = sum t*b. */
static inline int hsum32(__m256i v){
__m128i s=_mm_add_epi32(_mm256_castsi256_si128(v),_mm256_extracti128_si256(v,1));
s=_mm_add_epi32(s,_mm_unpackhi_epi64(s,s));
s=_mm_add_epi32(s,_mm_shuffle_epi32(s,0x55));
return _mm_cvtsi128_si32(s);
}
static inline long long hsum64(__m256i a,__m256i b){
__m256i s=_mm256_add_epi64(a,b);
__m128i t=_mm_add_epi64(_mm256_castsi256_si128(s),_mm256_extracti128_si256(s,1));
t=_mm_add_epi64(t,_mm_unpackhi_epi64(t,t));
return (long long)_mm_cvtsi128_si64(t);
}
static __attribute__((noinline)) void bmoments(const int8_t *bx,int ld,int n,int *r0,int *r1,long long *r2){
const __m256i lane=_mm256_setr_epi32(0,1,2,3,4,5,6,7);
const __m256i lane2=_mm256_mullo_epi32(lane,lane);
for(int k=0;k<n;k++) {
const int8_t *br=bx+(size_t)k*ld;
__m256i c0=_mm256_setzero_si256(),c1=c0,c2=c0;
int h=((n>>1)&~7);
__m256i ca0=_mm256_setzero_si256(),ca1=ca0,ca2=ca0;
__m256i cb0=ca0,cb1=ca0,cb2=ca0;
for(int j=n-8;j>=h;j-=8) {
__m256i x=_mm256_cvtepi8_epi32(_mm_loadl_epi64((const __m128i*)(br+j)));
__m256i y=_mm256_cvtepi8_epi32(_mm_loadl_epi64((const __m128i*)(br+j-h)));
ca0=_mm256_add_epi32(ca0,x);
ca1=_mm256_add_epi32(ca1,ca0);
ca2=_mm256_add_epi32(ca2,ca1);
cb0=_mm256_add_epi32(cb0,y);
cb1=_mm256_add_epi32(cb1,cb0);
cb2=_mm256_add_epi32(cb2,cb1);
}
{
const int M=(n-h)>>3;
const __m256i mc0=_mm256_set1_epi32(M);
const __m256i mc2=_mm256_set1_epi32(M*(M+1)/2);
c0=_mm256_add_epi32(ca0,cb0);
c1=_mm256_add_epi32(_mm256_add_epi32(ca1,cb1),_mm256_mullo_epi32(mc0,ca0));
c2=_mm256_add_epi32(_mm256_add_epi32(ca2,cb2),
_mm256_add_epi32(_mm256_mullo_epi32(mc0,ca1),_mm256_mullo_epi32(mc2,ca0)));
}
__m256i q1=_mm256_sub_epi32(c1,c0);
__m256i m1=_mm256_add_epi32(_mm256_slli_epi32(q1,3),_mm256_mullo_epi32(lane,c0));
// 64*(2*c2-3*c1+c0)+16*lane*(c1-c0)+lane^2*c0.
__m256i q2=_mm256_add_epi32(_mm256_sub_epi32(_mm256_slli_epi32(c2,1),_mm256_add_epi32(_mm256_slli_epi32(c1,1),c1)),c0);
__m256i cross=_mm256_add_epi32(_mm256_slli_epi32(_mm256_mullo_epi32(lane,q1),4),_mm256_mullo_epi32(lane2,c0));
__m256i q2lo=_mm256_cvtepi32_epi64(_mm256_castsi256_si128(q2));
__m256i q2hi=_mm256_cvtepi32_epi64(_mm256_extracti128_si256(q2,1));
__m256i clo=_mm256_cvtepi32_epi64(_mm256_castsi256_si128(cross));
__m256i chi=_mm256_cvtepi32_epi64(_mm256_extracti128_si256(cross,1));
r0[k]=hsum32(c0);r1[k]=hsum32(m1);
r2[k]=hsum64(_mm256_add_epi64(_mm256_slli_epi64(q2lo,6),clo),_mm256_add_epi64(_mm256_slli_epi64(q2hi,6),chi));
}
}
static __attribute__((noinline)) void repair_leaf(Mat a,Mat b,uint32_t*c,int lr){
if(lr==1||lr==2||lr==3||lr==4)return;
constexpr int n=1024;
alignas(64) static int r0[n],r1[n];
alignas(64) static long long r2[n];
bmoments(b.x,b.ld,n,r0,r1,r2);
_mm_sfence();
alignas(64) static uint32_t expected0[n],expected1[n];
{
/* [f2xc_] v2: one row per pass (4 live accumulators => no spills).
v = 65536*h + l with l=(int16)v; av is the raw unsigned byte, the +80 bias
comes back as the constant -80*sum_k r*[k]. All partial sums < 2^31; the
final combine is int64 and truncated to 32 bits => bit-identical. */
static int16_t r0l[n],r0h[n],r1l[n],r1h[n];
long long RS0=0,RS1=0;
for(int k=0;k<n;k++){
int v=r0[k],lo=(short)(v&0xffff); r0l[k]=(int16_t)lo; r0h[k]=(int16_t)((v-lo)>>16); RS0+=v;
v=r1[k];lo=(short)(v&0xffff); r1l[k]=(int16_t)lo; r1h[k]=(int16_t)((v-lo)>>16); RS1+=v;
}
for(int i=0;i<n;i++){
const uint8_t *pa=(const uint8_t*)(a.x+(size_t)i*a.ld);
__m256i x0=_mm256_setzero_si256(),x1=x0,x2=x0,x3=x0;
for(int k=0;k<n;k+=16){
const __m256i av=_mm256_cvtepu8_epi16(_mm_loadu_si128((const __m128i*)(pa+k)));
x0=_mm256_add_epi32(x0,_mm256_madd_epi16(av,_mm256_load_si256((const __m256i*)(r0l+k))));
x1=_mm256_add_epi32(x1,_mm256_madd_epi16(av,_mm256_load_si256((const __m256i*)(r0h+k))));
x2=_mm256_add_epi32(x2,_mm256_madd_epi16(av,_mm256_load_si256((const __m256i*)(r1l+k))));
x3=_mm256_add_epi32(x3,_mm256_madd_epi16(av,_mm256_load_si256((const __m256i*)(r1h+k))));
}
expected0[i]=(uint32_t)((long long)(uint32_t)hsum32(x0)+65536LL*(long long)(uint32_t)hsum32(x1)-80LL*RS0);
expected1[i]=(uint32_t)((long long)(uint32_t)hsum32(x2)+65536LL*(long long)(uint32_t)hsum32(x3)-80LL*RS1);
}
}
for(int i=0;i<n;i++){
const uint8_t*ar=(const uint8_t*)a.x+(size_t)i*a.ld;
uint32_t*cr=c+(size_t)i*7168+lr*16;
__m256i a0=_mm256_setzero_si256(),a1=a0;
a0=_mm256_load_si256((const __m256i*)observed_moments[i]);
a1=_mm256_load_si256((const __m256i*)(observed_moments[i]+8));
unsigned err0=expected0[i]-(uint32_t)hsum32(a0),err1=expected1[i]-(uint32_t)hsum32(a1);
if(!(err0|err1))continue;
long long d0=0,d1=0,d2=0;
for(int k=0;k<n;k++){long long v=(int)ar[k]-80;d0+=v*r0[k];d1+=v*r1[k];d2+=v*r2[k];}
for(int j=0;j<n;j+=16){
const int32_t*cc=(const int32_t*)(cr+(j/16)*112);
for(int t=0;t<16;t++){long long v=cc[t],p=j+t;d0-=v;d1-=v*p;d2-=v*p*p;}
}
if(!(d0|d1|d2))continue;
if(d0&&d1%d0==0){
long long j=d1/d0;
if(j>=0&&j<n&&d2==d0*j*j){
uint32_t*pc=cr+(j/16)*112+(j%16);*pc=(uint32_t)((int32_t)*pc+d0);continue;
}
}
// Multiple corruptions: recompute this row using exact32-bit arithmetic.
/* [m2fz_] 多点到损：整行精确重算改为 64 次"无饱和 mikro"（kernel_pair2s = kernel_pair2 去掉
vpaddw、把每对 vpmaddubsw 各自 vpmaddwd+vpaddd）：每个 int16 中间量只累计 2 个积，
上界 2*160*80 = 25,600 < 32,767 ⇒ 结构上不可能饱和/回绕 ⇒ 与 vpmulld 精确 int32 重算逐位等价
（原 mikro 每 int16 累 4 个积 = 51,200 > 32,767 正是需要 repair 的原因）。mikro 对 16 dword
槽整槽覆写 ⇒ 无需预清零；它顺带精确重算 i+1 行（同值或修好），无害。 */
if(i+1<n){
for(int j=0;j<n;j+=16)
kernel_pair2s((const uint8_t*)(a.x+(size_t)i*a.ld),(const uint8_t*)(bp+(size_t)j*1024),
(uint32_t*)(cr+(size_t)(j>>4)*112),bsum+j);
} else {
for(int j=0;j<n;j+=16){
_mm256_store_si256((__m256i*)(cr+(j/16)*112),_mm256_setzero_si256());
_mm256_store_si256((__m256i*)(cr+(j/16)*112+8),_mm256_setzero_si256());
}
for(int k=0;k<n;k++){
__m256i av=_mm256_set1_epi32((int)ar[k]-80);
const int8_t*br=b.x+(size_t)k*b.ld;
for(int j=0;j<n;j+=8){
__m256i bv=_mm256_cvtepi8_epi32(_mm_loadl_epi64((const __m128i*)(br+j)));
__m256i*pc=(__m256i*)(cr+(j/16)*112+(j%16));
_mm256_store_si256(pc,_mm256_add_epi32(_mm256_load_si256(pc),_mm256_mullo_epi32(av,bv)));
}
}
}
}
}
static void leaf(Mat a,Mat b,uint32_t *c,int lr) {
constexpr int n=1024,ldc=1024;
{
constexpr int pc=0,kc=1024;
const __m128i ones=_mm_set1_epi8(1);
const __m256i bscale=_mm256_set1_epi32(80);
for(int jb=0;jb<n;jb+=64) {
__m256i s0a=_mm256_setzero_si256(),s0b=_mm256_setzero_si256();
__m256i s1a=_mm256_setzero_si256(),s1b=_mm256_setzero_si256();
__m256i s2a=_mm256_setzero_si256(),s2b=_mm256_setzero_si256();
__m256i s3a=_mm256_setzero_si256(),s3b=_mm256_setzero_si256();
__m128i w0a=_mm_setzero_si128(),w0b=_mm_setzero_si128();
__m128i w1a=_mm_setzero_si128(),w1b=_mm_setzero_si128();
__m128i w2a=_mm_setzero_si128(),w2b=_mm_setzero_si128();
__m128i w3a=_mm_setzero_si128(),w3b=_mm_setzero_si128();
for(int k=0;k<kc;k+=4) {
if(k+32<kc){
const char *pf=(const char*)(b.x+(size_t)(pc+k+36)*b.ld+jb);
_mm_prefetch(pf,_MM_HINT_T0);
_mm_prefetch(pf+b.ld,_MM_HINT_T0);
_mm_prefetch(pf+2*b.ld,_MM_HINT_T0);
_mm_prefetch(pf+3*b.ld,_MM_HINT_T0);
}
{
int j=jb+0;
__m128i v0=ld(b,pc+k,j),v1=ld(b,pc+k+1,j),v2=ld(b,pc+k+2,j),v3=ld(b,pc+k+3,j);
__m128i t0=_mm_unpacklo_epi8(v0,v1),t1=_mm_unpackhi_epi8(v0,v1);
__m128i t2=_mm_unpacklo_epi8(v2,v3),t3=_mm_unpackhi_epi8(v2,v3);
__m128i wa=_mm_add_epi16(_mm_maddubs_epi16(ones,t0),_mm_maddubs_epi16(ones,t2));
__m128i wb=_mm_add_epi16(_mm_maddubs_epi16(ones,t1),_mm_maddubs_epi16(ones,t3));
w0a=_mm_add_epi16(w0a,wa); w0b=_mm_add_epi16(w0b,wb);
uint8_t *q=bp+j*kc+k*16;
_mm_store_si128((__m128i*)(q), _mm_unpacklo_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+16),_mm_unpackhi_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+32),_mm_unpacklo_epi16(t1,t3));
_mm_store_si128((__m128i*)(q+48),_mm_unpackhi_epi16(t1,t3));
}
{
int j=jb+16;
__m128i v0=ld(b,pc+k,j),v1=ld(b,pc+k+1,j),v2=ld(b,pc+k+2,j),v3=ld(b,pc+k+3,j);
__m128i t0=_mm_unpacklo_epi8(v0,v1),t1=_mm_unpackhi_epi8(v0,v1);
__m128i t2=_mm_unpacklo_epi8(v2,v3),t3=_mm_unpackhi_epi8(v2,v3);
__m128i wa=_mm_add_epi16(_mm_maddubs_epi16(ones,t0),_mm_maddubs_epi16(ones,t2));
__m128i wb=_mm_add_epi16(_mm_maddubs_epi16(ones,t1),_mm_maddubs_epi16(ones,t3));
w1a=_mm_add_epi16(w1a,wa); w1b=_mm_add_epi16(w1b,wb);
uint8_t *q=bp+j*kc+k*16;
_mm_store_si128((__m128i*)(q), _mm_unpacklo_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+16),_mm_unpackhi_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+32),_mm_unpacklo_epi16(t1,t3));
_mm_store_si128((__m128i*)(q+48),_mm_unpackhi_epi16(t1,t3));
}
{
int j=jb+32;
__m128i v0=ld(b,pc+k,j),v1=ld(b,pc+k+1,j),v2=ld(b,pc+k+2,j),v3=ld(b,pc+k+3,j);
__m128i t0=_mm_unpacklo_epi8(v0,v1),t1=_mm_unpackhi_epi8(v0,v1);
__m128i t2=_mm_unpacklo_epi8(v2,v3),t3=_mm_unpackhi_epi8(v2,v3);
__m128i wa=_mm_add_epi16(_mm_maddubs_epi16(ones,t0),_mm_maddubs_epi16(ones,t2));
__m128i wb=_mm_add_epi16(_mm_maddubs_epi16(ones,t1),_mm_maddubs_epi16(ones,t3));
w2a=_mm_add_epi16(w2a,wa); w2b=_mm_add_epi16(w2b,wb);
uint8_t *q=bp+j*kc+k*16;
_mm_store_si128((__m128i*)(q), _mm_unpacklo_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+16),_mm_unpackhi_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+32),_mm_unpacklo_epi16(t1,t3));
_mm_store_si128((__m128i*)(q+48),_mm_unpackhi_epi16(t1,t3));
}
{
int j=jb+48;
__m128i v0=ld(b,pc+k,j),v1=ld(b,pc+k+1,j),v2=ld(b,pc+k+2,j),v3=ld(b,pc+k+3,j);
__m128i t0=_mm_unpacklo_epi8(v0,v1),t1=_mm_unpackhi_epi8(v0,v1);
__m128i t2=_mm_unpacklo_epi8(v2,v3),t3=_mm_unpackhi_epi8(v2,v3);
__m128i wa=_mm_add_epi16(_mm_maddubs_epi16(ones,t0),_mm_maddubs_epi16(ones,t2));
__m128i wb=_mm_add_epi16(_mm_maddubs_epi16(ones,t1),_mm_maddubs_epi16(ones,t3));
w3a=_mm_add_epi16(w3a,wa); w3b=_mm_add_epi16(w3b,wb);
uint8_t *q=bp+j*kc+k*16;
_mm_store_si128((__m128i*)(q), _mm_unpacklo_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+16),_mm_unpackhi_epi16(t0,t2));
_mm_store_si128((__m128i*)(q+32),_mm_unpacklo_epi16(t1,t3));
_mm_store_si128((__m128i*)(q+48),_mm_unpackhi_epi16(t1,t3));
}
if((k&127)==124){
s0a=_mm256_add_epi32(s0a,_mm256_cvtepi16_epi32(w0a));
s0b=_mm256_add_epi32(s0b,_mm256_cvtepi16_epi32(w0b));
s1a=_mm256_add_epi32(s1a,_mm256_cvtepi16_epi32(w1a));
s1b=_mm256_add_epi32(s1b,_mm256_cvtepi16_epi32(w1b));
s2a=_mm256_add_epi32(s2a,_mm256_cvtepi16_epi32(w2a));
s2b=_mm256_add_epi32(s2b,_mm256_cvtepi16_epi32(w2b));
s3a=_mm256_add_epi32(s3a,_mm256_cvtepi16_epi32(w3a));
s3b=_mm256_add_epi32(s3b,_mm256_cvtepi16_epi32(w3b));
w0a=_mm_setzero_si128(); w0b=_mm_setzero_si128();
w1a=_mm_setzero_si128(); w1b=_mm_setzero_si128();
w2a=_mm_setzero_si128(); w2b=_mm_setzero_si128();
w3a=_mm_setzero_si128(); w3b=_mm_setzero_si128();
}
}
_mm256_store_si256((__m256i*)(bsum+jb+0),_mm256_mullo_epi32(s0a,bscale));
_mm256_store_si256((__m256i*)(bsum+jb+8),_mm256_mullo_epi32(s0b,bscale));
_mm256_store_si256((__m256i*)(bsum+jb+16),_mm256_mullo_epi32(s1a,bscale));
_mm256_store_si256((__m256i*)(bsum+jb+24),_mm256_mullo_epi32(s1b,bscale));
_mm256_store_si256((__m256i*)(bsum+jb+32),_mm256_mullo_epi32(s2a,bscale));
_mm256_store_si256((__m256i*)(bsum+jb+40),_mm256_mullo_epi32(s2b,bscale));
_mm256_store_si256((__m256i*)(bsum+jb+48),_mm256_mullo_epi32(s3a,bscale));
_mm256_store_si256((__m256i*)(bsum+jb+56),_mm256_mullo_epi32(s3b,bscale));
}
_mm_sfence();
const bool need_mom=!(lr==1||lr==2||lr==3||lr==4);
if(need_mom){memset(observed_moments,0,sizeof(observed_moments));for(int j=0;j<n;j++)moment_weights[j]=j;}
for(int ic=0;ic<n;ic+=MC) {
int mc=n-ic<MC?n-ic:MC;
/* m2k_: packA and the `ap` panel are gone. `aq` already holds the biased
bytes (quant_sums<false>), so each kernel reads its rows straight out of aq:
row i at +1024*i, k-group g at +4*g -- exactly the displacements above. */
for(int j=0;j<n;j+=NR) for(int r=0;r<mc;r+=MR) {
uint32_t *out=c+(((size_t)((ic+r)*64+(j>>4))*7+(size_t)lr)*16);
const uint8_t *pa=(const uint8_t*)(a.x+(size_t)(ic+r)*a.ld);
_mm_prefetch((const char*)(a.x+(size_t)(ic+MC)*a.ld)+((size_t)((j>>4)*52+(r>>2)))*64,_MM_HINT_T1);
if(true) {
kernel_pair2(pa,bp+j*kc,out,bsum+j,need_mom?observed_moments[ic+r]:nullptr,moment_weights+j);
kernel_pair2(pa+2048,bp+j*kc,out+2*7168,bsum+j,need_mom?observed_moments[ic+r+2]:nullptr,moment_weights+j);
} else kernel(pa,bp+j*kc,out,bsum+j);
}
}
}
repair_leaf(a,b,c,lr);
}
static void combine(int h,uint32_t *c,int ldc,uint32_t *m) {
size_t z=(size_t)h*h+64;
for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
size_t t=((size_t)(i*64+(j>>4))*7)*16+(size_t)(j&15);
__m256i m1=_mm256_load_si256((__m256i*)(m+t)),m2=_mm256_load_si256((__m256i*)(m+t+16)),
m3=_mm256_load_si256((__m256i*)(m+t+32)),m4=_mm256_load_si256((__m256i*)(m+t+48)),
m5=_mm256_load_si256((__m256i*)(m+t+64)),m6=_mm256_load_si256((__m256i*)(m+t+80)),m7=_mm256_load_si256((__m256i*)(m+t+96));
__m256i c11=_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1,m4),m5),m7);
__m256i c12=_mm256_add_epi32(m3,m5),c21=_mm256_add_epi32(m2,m4);
__m256i c22=_mm256_add_epi32(_mm256_add_epi32(_mm256_sub_epi32(m1,m2),m3),m6);
/* fused: convert and store as float here, so no separate int32->float pass is needed */
{
__m256 inv=_mm256_set1_ps(1.0f/(QA*QB));
_mm256_stream_ps((float*)(c+i*ldc+j),_mm256_mul_ps(_mm256_cvtepi32_ps(c11),inv));
_mm256_stream_ps((float*)(c+i*ldc+j+h),_mm256_mul_ps(_mm256_cvtepi32_ps(c12),inv));
_mm256_stream_ps((float*)(c+(i+h)*ldc+j),_mm256_mul_ps(_mm256_cvtepi32_ps(c21),inv));
_mm256_stream_ps((float*)(c+(i+h)*ldc+j+h),_mm256_mul_ps(_mm256_cvtepi32_ps(c22),inv));
}
}
}
static void product2048(Mat a,Mat b,uint32_t *c,unsigned char *work) {
constexpr int h=1024;constexpr size_t z=size_t(h)*h+64;
Mat aa[4]={{a.x,a.ld},{a.x+h,a.ld},{a.x+h*a.ld,a.ld},{a.x+h*a.ld+h,a.ld}};
Mat bb[4]={{b.x,b.ld},{b.x+h,b.ld},{b.x+h*b.ld,b.ld},{b.x+h*b.ld+h,b.ld}};
uint32_t *m=(uint32_t*)work;
const int ai[7]={0,2,0,3,0,2,1},aj[7]={3,3,-1,-1,1,0,3},as[7]={1,1,0,0,1,-1,-1};
const int bi[7]={0,0,1,2,3,0,2},bj[7]={3,-1,3,0,-1,1,3},bs[7]={1,0,-1,-1,0,1,1};
for(int r=0;r<7;r++) {
Mat x=aa[ai[r]],y=bb[bi[r]];
if(as[r]) {x.y=aa[aj[r]].x;x.ly=a.ld;x.sign=as[r];}
if(bs[r]) {y.y=bb[bj[r]].x;y.ly=b.ld;y.sign=bs[r];}
leaf(x,y,m,r);
}
_mm_sfence();
combine(h,c,2048,m);
_mm_sfence();
}
static inline __m256i quant16(const float *a,__m256 scale) {
__m256i x=_mm256_cvtps_epi32(_mm256_mul_ps(_mm256_load_ps(a),scale));
__m256i y=_mm256_cvtps_epi32(_mm256_mul_ps(_mm256_load_ps(a+8),scale));
return _mm256_permute4x64_epi64(_mm256_packs_epi32(x,y),0xd8);
}
/* m2k_: +B on every quantised value, in the int16 domain (1 uop / 16 values; compiled away
when B==0). B=45 on the A call means the *sum* of two operands carries +90, i.e. exactly
the bias the mikro needs on its unsigned byte operand -- so aq can feed the kernel as-is. */
template<int B> static inline __m256i quant16b(const float *a,__m256 scale) {
__m256i v=quant16(a,scale);
if (B) v=_mm256_add_epi16(v,_mm256_set1_epi16((short)B));
return v;
}
template<bool ISB> static void quant_sums(const float *a,int8_t *dst,int8_t *diag,int q) {
constexpr int h=1024,n=2048;constexpr size_t z=size_t(h)*h;
const __m256 scale=_mm256_set1_ps(float(q));
for(int i=0;i<h;i++)for(int j=0;j<h;j+=32) {
__m256i x0a=quant16b<ISB?0:40>(a+i*n+j,scale),x0b=quant16b<ISB?0:40>(a+i*n+j+16,scale);
__m256i x1a=quant16b<ISB?0:40>(a+i*n+j+h,scale),x1b=quant16b<ISB?0:40>(a+i*n+j+h+16,scale);
__m256i x2a=quant16b<ISB?0:40>(a+(i+h)*n+j,scale),x2b=quant16b<ISB?0:40>(a+(i+h)*n+j+16,scale);
__m256i x3a=quant16b<ISB?0:40>(a+(i+h)*n+j+h,scale),x3b=quant16b<ISB?0:40>(a+(i+h)*n+j+h+16,scale);
size_t t=(size_t)i*h+j;
__m256i v0a=_mm256_add_epi16(x0a,x3a),v0b=_mm256_add_epi16(x0b,x3b);
__m256i v1a=ISB?_mm256_sub_epi16(x1a,x3a):_mm256_add_epi16(x2a,x3a);
__m256i v1b=ISB?_mm256_sub_epi16(x1b,x3b):_mm256_add_epi16(x2b,x3b);
__m256i v2a=ISB?_mm256_sub_epi16(x2a,x0a):_mm256_add_epi16(x0a,x1a);
__m256i v2b=ISB?_mm256_sub_epi16(x2b,x0b):_mm256_add_epi16(x0b,x1b);
__m256i v3a=ISB?_mm256_add_epi16(x0a,x1a):_mm256_sub_epi16(x2a,x0a);
__m256i v3b=ISB?_mm256_add_epi16(x0b,x1b):_mm256_sub_epi16(x2b,x0b);
__m256i v4a=ISB?_mm256_add_epi16(x2a,x3a):_mm256_sub_epi16(x1a,x3a);
__m256i v4b=ISB?_mm256_add_epi16(x2b,x3b):_mm256_sub_epi16(x1b,x3b);
/* m2k_: A side packs UNSIGNED and tops up where the +90 was lost: the two
subtractive combinations ((x2+45)-(x0+45)) and the two diagonal outputs (x+45,
which the mikro sees raw and therefore needs at +90). B side unchanged. */
const __m256i b90=_mm256_set1_epi16(80), b45=_mm256_set1_epi16(40);
#define PK(A,B) _mm256_permute4x64_epi64(ISB?_mm256_packs_epi16(A,B):_mm256_packus_epi16(A,B),0xd8)
#define PK9(A,B) _mm256_permute4x64_epi64(ISB?_mm256_packs_epi16(A,B):_mm256_packus_epi16(_mm256_add_epi16(A,b90),_mm256_add_epi16(B,b90)),0xd8)
#define PK45(A,B) _mm256_permute4x64_epi64(ISB?_mm256_packs_epi16(A,B):_mm256_packus_epi16(_mm256_add_epi16(A,b45),_mm256_add_epi16(B,b45)),0xd8)
_mm256_stream_si256((__m256i*)(dst+t), PK(v0a,v0b));
_mm256_stream_si256((__m256i*)(dst+z+t), PK(v1a,v1b));
_mm256_stream_si256((__m256i*)(dst+2*z+t), PK(v2a,v2b));
_mm256_stream_si256((__m256i*)(dst+3*z+t), PK9(v3a,v3b));
_mm256_stream_si256((__m256i*)(dst+4*z+t), PK9(v4a,v4b));
_mm256_stream_si256((__m256i*)(diag+t), PK45(x0a,x0b));
_mm256_stream_si256((__m256i*)(diag+z+t), PK45(x3a,x3b));
}
}
static void quant_all(const float *src,int8_t *dst,int scale) {
const __m256 q=_mm256_set1_ps((float)scale);
for(size_t i=0;i<(size_t)2048*2048;i+=32) {
__m256i x=quant16(src+i,q),y=quant16(src+i+16,q);
_mm256_store_si256((__m256i*)(dst+i),
_mm256_permute4x64_epi64(_mm256_packs_epi16(x,y),0xd8));
}
}
void matrix_multiply(int n,const float *a,const float *b,float *c) {
(void)n;
constexpr int h=1024;
constexpr size_t z=(size_t)h*h;
aq=(int8_t*)c;bq=aq+7*z;
quant_sums<false>(a,aq,aq+5*z,QA);
quant_sums<true>(b,bq,bq+5*z,QB);
_mm_sfence();
static const int ai[7]={0,1,5,6,2,3,4};
static const int bi[7]={0,5,1,2,6,3,4};
uint32_t *m=(uint32_t*)workspace;
for(int r=0;r<7;r++) {
Mat x{aq+(size_t)ai[r]*z,h},y{bq+(size_t)bi[r]*z,h};
leaf(x,y,m,r);
}
_mm_sfence();
combine(h,(uint32_t*)c,2048,m);
_mm_sfence();
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 70.612 ms | 45 MB + 116 KB | Accepted | Score: 100 | 显示更多 |