代码

/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/126504: final-four-quadrant fusion with old-product read-ahead and safe four-row tail-A retirement; inherited numerical and source credits below.
 * Idea: M5 is never materialized by this fused algorithm, so remove its empty slot from every16-column product record. Store M1/M2/M3/M4/M6/M7 in six contiguous cache lines, reducing the arena28MiB to24MiB. Remap all main, fallback repair, two-row repair, read-ahead and final reconstruction addresses consistently; logical product IDs and repair decisions remain unchanged.
 * Purpose: Remove4MiB of unused physical page footprint and improve dense six-product locality without changing any arithmetic or final output schedule.
 */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/126501 : the proven maximal
 * tail-A retirement schedule and four-quadrant final product engine.
 * All earlier author, source and license notices remain below.
 * Idea: Before each final4-row/16-column dot-product kernel, prefetch the six
 * prior product cache lines for each of its four rows. The multiplication
 * supplies lead time before the epilogue first reads these otherwise cold
 * strided product rows. Exactly24 valid64-byte lines are hinted, omitting M5;
 * no extra data is read semantically and no arithmetic or assembly changes.
 * Experiment purpose: restore timely old-product availability after fusion
 * changed reconstruction from a linear pass to a column-panel traversal.
 */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/126496 : direct tail A,
 * four-quadrant final product, and all inherited credits retained below.
 * Idea: Choose each final direct-input panel by the exact overwrite bound,
 * not a guessed tile size. For a panel beginning at a and ending at b,
 * it is safe without copying when 8*b-7161<a. After four208-row panels,
 * rows832..999 (168 rows) and1000..1019 (20 rows) meet this bound. Copy only
 * the last four rows1020..1023 (4KiB) into the retired moment buffer. This
 * preserves more large-panel work while minimizing saved input traffic;
 * every kernel and all integer/float arithmetic are unchanged.
 * Experiment purpose: test a mathematically derived monotone retirement
 * schedule using the same direct-tail production and existing scratch.
 */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/126496 : direct tail A
 * production and retired moment staging, preserving all prior credits below.
 * Idea: Keep the first four 208-row final-product panels on their original
 * contiguous A source, then process rows832..959 directly as one128-row panel.
 * Their tail output overwrites at most A row511, already retired. Only rows
 * 960..1023 need a64KiB save to observed_moments before the last panel. This
 * removes15 of16 panel copies and retains the original large-panel schedule
 * for832 rows. All final arithmetic and assembly remain unchanged; bounded
 * read-ahead points only into the remaining original tail A plane.
 * Experiment purpose: reduce saved-input traffic and preserve cache reuse
 * while maintaining the proven monotone output/input overlap invariant.
 */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/126488 : four-quadrant final
 * kernel and inherited engine; full earlier attribution is retained below.
 * Idea: The A plane numbered 2 is consumed only by final product M5. Produce
 * it directly at C+15 MiB during A quantization instead of first at C+2 MiB,
 * avoiding its subsequent 1 MiB relocation. B's plane2 and all other planes
 * retain their original locations. The existing retired-moment panel copies
 * and monotone overwrite proof still protect the final contiguous A input.
 * Experiment purpose: remove a full temporary-plane copy without changing
 * any quantized value, arithmetic instruction or final kernel.
 */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/126488 : retain final M5,
 * the four-quadrant epilogue, all quantization and checksum/repair arithmetic.
 * All inherited references, author notices and licenses remain below.
 * Idea: Store final A in the last 1 MiB of C, beyond the 14 MiB quantized
 * inputs, and copy each final 64-row panel to the retired 64 KiB observed
 * moment buffer before computing it. No new array or memory pages are needed;
 * the existing kernel retains contiguous 1024-byte A rows. Outputs touch this
 * A region only for logical rows i>=896 and overwrite A rows
 * [8*(i-896),8*(i-896)+7], whose upper bound is at most i for i<=1023.
 * Thus only completed or currently saved A panels are overwritten, never a
 * future panel. Final B is completely packed before any such output.
 * Experiment purpose: avoid the private A page footprint without sparse A
 * addressing; only the final product uses the 64-row staging loop.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126480: final-M5 scheduling and direct top-half output, retaining all engine credits below.
Idea: Extend the final epilogue to all four quadrants. Load each of M1/M2/M3/M4/M6/M7 once and reuse it for both dependent results, removing the bottom-half pass that duplicated three complete product reads. Integer sum order, signed conversion and floating scale match the original combine. All C-backed inputs are dead after final B packing and the private A copy.
Purpose: Eliminate the M5 round trip without repeated product reads. An identical assembly-core string macro and removal of uncalled helpers keep the source within its size limit; they do not alter those existing instructions. */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/126462 : retain the validated
 * lane-local quantization, seeded moments, seven products and repair engine.
 * Earlier authors, sources and license notices are preserved below.
 * Idea: Evaluate product M5 last. Save its 1 MiB A plane outside C, then pack
 * its complete B plane before writing any final result. After B is packed, load all six prior products once and consume the final
 * M5 registers to produce all four output quadrants in one epilogue.
 * This avoids materializing and reading M5's 4 MiB product. The other six
 * products and every applicable checksum/repair remain unchanged. Experiment
 * purpose: replace one product round trip with a single input-plane copy.
 */
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126457 : the accepted lane-local intermediate quantization idea, combined here with the accepted observed-seeded residuals.
- saffah_codex_6a_agg3, https://duck.ac/submission/126453 : copied the current quad-byte row moments and exact quantized engine, retaining all inherited citations below.
Idea:
Keep quantized 16-bit values in the natural lane-local order produced by
packs_epi32 throughout the elementwise Strassen sums/differences. Correct the
order only once after final byte packing, using the fixed dword permutation
[0,4,1,5,2,6,3,7]. This removes eight intermediate cross-lane shuffles per
32-element group. Rounding, saturation, biases, stores, and final byte order
are exactly the same; all existing assembly is untouched.
Experimental purpose: reduce quantization shuffle-port work and dependencies.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126448 : current accepted residual-reduction and exact-repair engine; copied with all inherited references below.
Idea:
Initialize the two low-factor dot accumulators with the negated observed
moments. Distributivity modulo 2^32 makes their final reductions identical
to subtracting the observed vectors after the dot loops. This moves the two
observed loads/subtractions away from the reduction critical path while
using the same live accumulators; all repair decisions remain identical.
Experimental purpose: reduce late dependency latency in the ordinary checksum
screen. Existing matrix kernels and assembly remain untouched.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126443 : copied the current residual-reduction engine and complete checksum repair, retaining all inherited citations below.
Idea:
Compute B row moments using 32-byte blocks split into eight four-byte groups.
Unsigned local weights {1}, {0,1,2,3}, {0,1,4,9} times signed input bytes
produce exact local moments with nonsaturating maddubs/maddwd. Six prefix sums
then recover global zeroth/first/second moments by polynomial identities;
second-moment high terms are widened before multiplication by 1024.
Experimental purpose: halve the serial prefix iteration count without changing
matrix arithmetic, precision, checks or repair; existing assembly is untouched.
*/
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126412: retained the accepted quantized engine, split16 moments, observed-moment recording and complete repair; all inherited credits follow.
Idea: Linearity modulo 2^32 lets the expected split low/high dot accumulators and the observed lane sums combine before horizontal reduction. Form low+(high<<16)-observed as eight vector lanes, reduce once for each of the two residual moments, and then subtract the unchanged scalar input bias. Their exact zero tests and all full-width repair remain unchanged.
Purpose: Replace repeated horizontal sums of individual expected/observed terms with two final residual reductions. Matrix kernels and every existing assembly region are unchanged.
*/
/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/126408: retained the accepted inlined paired tiles, invariant recording dispatch, negative bias and complete quantized/repair engine; inherited credits remain below.
 * Idea: Load the already-negative B column bias directly into the first row accumulators, then copy it to the remaining rows. Eliminate the two temporary bias-to-first-row moves from each hot tile and the exact repair kernel. Every multiplication starts from the identical values and all loop arithmetic is unchanged.
 * Purpose: Reduce microtile initialization instructions without altering precision or repair behavior.
 */
/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/126398: retained leaf-level moment dispatch, negative-bias paired kernels, full repair and all inherited credits.
 * Idea: Inline the two hot four-row kernels directly into their leaf loops, while explicitly keeping one inner-row loop body. This removes the repeated call boundary, argument setup and intermediate vzeroupper, without multiplying the large assembly body by compiler row-loop unrolling. The rare exact-repair kernel remains outlined; all assembly arithmetic and output predicates are unchanged.
 * Purpose: Reduce hot microtile call overhead with controlled instruction footprint and bit-identical output.
 */
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/126387: retained negative-bias four-row kernels, moments, exact repair and all inherited references.
* Idea: Select whether the current leaf needs moment recording once outside both microtile loops. Specialized four-row kernels then either always record or contain no moment epilogue, eliminating the per-call null branch and unused moment/weight arguments. Remove the already-unreachable legacy kernel without changing any active arithmetic or repair.
* Purpose: Reduce invariant dispatch and call setup across all microtiles while preserving every output bit.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/126377: retained four-row paired kernel, complete quantized algorithm, checksum repair and inherited citations.
* saffah_codex_6a_agg3 https://duck.ac/submission/126365: reused the idea of storing the required signed initial correction directly instead of repeating subtraction at each microkernel entry.
* Idea: Materialize negative80-times B column sums once during packing. Every matrix and repair kernel then copies this already-correct initial accumulator bias instead of subtracting positive sums from zero for each output row. All integer values entering the multiplication loops are unchanged.
* Purpose: Remove repeated vector bias subtraction and zero setup from thousands of kernel calls while preserving exact output bits.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/125990: retained the accepted quantized Strassen engine, narrow moment checks, exact repair, and inherited citations.
* Idea: Process the four adjacent output rows in one paired-byte microkernel with eight independent accumulators. Share bias initialization, weight loads and the optional moment branch across all four rows; keep the identical paired signed-word wrap and per-row reduction order. Reuse the existing complete exact-repair path without changes.
* Purpose: Remove duplicate two-row call setup and epilogue work while preserving every output bit.
*/
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125917: exact repair skipping and the accepted quantized matrix engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/123940: inherited split16-bit moment dot products and independent B moment setup; complete source notices retained below.
No separate software license notice was present in those public submissions.
Idea: Check once while preparing B moments whether every zeroth row moment fits signed16bits. In that case its split high word is identically zero, so dispatch to a three-dot-product checksum loop that omits the fourth zero product. The general four-product loop remains for any larger row sum. First moments, full-width repair and all quantized arithmetic remain exact.
Purpose: Remove a provably zero dot-product stream on bounded row sums with a complete exact fallback.
*/
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/123940: current accepted quantized matrix engine and split-16-bit checksum dot products, copied and adapted.
- saffah_codex_6a_agg3, https://duck.ac/submission/121727: original engine ancestry, explicitly cited by the reference implementation.
The public submissions declare no separate software license. All inherited comments, including their full references and credit notices, follow this reference block before the new idea and code.
*/
/* ===== [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）。
* ================
*/
/* m2k_: the A panel `ap` is gone */
// 取 16 个 int8；x±y 都在 int8 域完成（各源 |v|<=46 => 和 |v|<=92 <=127，不溢出）
// 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. */
// 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. */
// 64*(2*c2-3*c1+c0)+16*lane*(c1-c0)+lane^2*c0.
// 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 行（同值或修好），无害。 */
// both consecutive rows are now exact
/* 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. */
/* fused: convert and store as float here, so no separate int32->float pass is needed */
/* 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. */
/* 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. */
/* Idea: Compute each expected checksum immediately before testing that row, eliminating the expected0/1 arrays and permitting work to be skipped. The exact pair repair already recomputes both i and i+1, so advance past i+1 after that repair; its old observed checksum is stale and cannot require another correction. Every unskipped row still uses the complete original checksum/repair logic. Experiment purpose: avoid duplicate verification of rows already exactly repaired, and retain the row operand in cache for exceptional repair. */
#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[24*1024*1024];
alignas(64) static uint8_t bp[1024*1024];   
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; };
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;
}
alignas(64) static uint32_t observed_moments[1024][16];
alignas(64) static int moment_weights[1024];
#define F2_NOMOM_CORE \
"vmovdqa (%[sums]), %%ymm0\n\t" \
"vmovdqa 32(%[sums]), %%ymm1\n\t" \
"vmovdqa %%ymm0, %%ymm2\n\t" \
"vmovdqa %%ymm1, %%ymm3\n\t" \
"vmovdqa %%ymm0, %%ymm4\n\t" \
"vmovdqa %%ymm1, %%ymm5\n\t" \
"vmovdqa %%ymm0, %%ymm6\n\t" \
"vmovdqa %%ymm1, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2048(%[a]), %%ymm10\n\t" \
"vpmaddubsw 0(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 32(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2052(%[a]), %%ymm10\n\t" \
"vpmaddubsw 64(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3072(%[a]), %%ymm10\n\t" \
"vpmaddubsw 0(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 32(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3076(%[a]), %%ymm10\n\t" \
"vpmaddubsw 64(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 96(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2056(%[a]), %%ymm10\n\t" \
"vpmaddubsw 128(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 160(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2060(%[a]), %%ymm10\n\t" \
"vpmaddubsw 192(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3080(%[a]), %%ymm10\n\t" \
"vpmaddubsw 128(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 160(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3084(%[a]), %%ymm10\n\t" \
"vpmaddubsw 192(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 224(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2064(%[a]), %%ymm10\n\t" \
"vpmaddubsw 256(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 288(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2068(%[a]), %%ymm10\n\t" \
"vpmaddubsw 320(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3088(%[a]), %%ymm10\n\t" \
"vpmaddubsw 256(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 288(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3092(%[a]), %%ymm10\n\t" \
"vpmaddubsw 320(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 352(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2072(%[a]), %%ymm10\n\t" \
"vpmaddubsw 384(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 416(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2076(%[a]), %%ymm10\n\t" \
"vpmaddubsw 448(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3096(%[a]), %%ymm10\n\t" \
"vpmaddubsw 384(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 416(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3100(%[a]), %%ymm10\n\t" \
"vpmaddubsw 448(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 480(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2080(%[a]), %%ymm10\n\t" \
"vpmaddubsw 512(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 544(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2084(%[a]), %%ymm10\n\t" \
"vpmaddubsw 576(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3104(%[a]), %%ymm10\n\t" \
"vpmaddubsw 512(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 544(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3108(%[a]), %%ymm10\n\t" \
"vpmaddubsw 576(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 608(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2088(%[a]), %%ymm10\n\t" \
"vpmaddubsw 640(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 672(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2092(%[a]), %%ymm10\n\t" \
"vpmaddubsw 704(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3112(%[a]), %%ymm10\n\t" \
"vpmaddubsw 640(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 672(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3116(%[a]), %%ymm10\n\t" \
"vpmaddubsw 704(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 736(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2096(%[a]), %%ymm10\n\t" \
"vpmaddubsw 768(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 800(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2100(%[a]), %%ymm10\n\t" \
"vpmaddubsw 832(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3120(%[a]), %%ymm10\n\t" \
"vpmaddubsw 768(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 800(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3124(%[a]), %%ymm10\n\t" \
"vpmaddubsw 832(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 864(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\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, %%ymm11\n\t" \
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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, %%ymm11\n\t" \
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%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 2104(%[a]), %%ymm10\n\t" \
"vpmaddubsw 896(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 928(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 2108(%[a]), %%ymm10\n\t" \
"vpmaddubsw 960(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t" \
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t" \
"vpbroadcastd 3128(%[a]), %%ymm10\n\t" \
"vpmaddubsw 896(%[b]), %%ymm10, %%ymm8\n\t" \
"vpmaddubsw 928(%[b]), %%ymm10, %%ymm9\n\t" \
"vpbroadcastd 3132(%[a]), %%ymm10\n\t" \
"vpmaddubsw 960(%[b]), %%ymm10, %%ymm11\n\t" \
"vpmaddubsw 992(%[b]), %%ymm10, %%ymm12\n\t" \
"vpaddw %%ymm11, %%ymm8, %%ymm8\n\t" \
"vpaddw %%ymm12, %%ymm9, %%ymm9\n\t" \
"vpmaddwd %%ymm15, %%ymm8, %%ymm8\n\t" \
"vpmaddwd %%ymm15, %%ymm9, %%ymm9\n\t" \
"vpaddd %%ymm8, %%ymm6, %%ymm6\n\t" \
"vpaddd %%ymm9, %%ymm7, %%ymm7\n\t" \
"add $64, %[a]\n\t" \
"add $1024, %[b]\n\t" \
"sub $16, %[k]\n\t" \
"jnz 1b\n\t"
static __attribute__((always_inline)) inline void kernel_pair4(const uint8_t*a,const uint8_t*b,uint32_t*out,const int*sums,uint32_t*mom,const int*weight){
int k=256;
asm volatile(
F2_NOMOM_CORE
"vmovntdq %%ymm0, 0(%[out])\n\t"
"vmovntdq %%ymm1, 32(%[out])\n\t"
"vmovntdq %%ymm2, 24576(%[out])\n\t"
"vmovntdq %%ymm3, 24608(%[out])\n\t"
"vmovntdq %%ymm4, 49152(%[out])\n\t"
"vmovntdq %%ymm5, 49184(%[out])\n\t"
"vmovntdq %%ymm6, 73728(%[out])\n\t"
"vmovntdq %%ymm7, 73760(%[out])\n\t"
"vmovdqa (%[weight]), %%ymm8\n\t"
"vmovdqa 32(%[weight]), %%ymm9\n\t"
"vpaddd %%ymm0, %%ymm1, %%ymm10\n\t"
"vpaddd 0(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 0(%[mom])\n\t"
"vpmulld %%ymm8, %%ymm0, %%ymm0\n\t"
"vpmulld %%ymm9, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm0, %%ymm1, %%ymm10\n\t"
"vpaddd 32(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 32(%[mom])\n\t"
"vpaddd %%ymm2, %%ymm3, %%ymm10\n\t"
"vpaddd 64(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 64(%[mom])\n\t"
"vpmulld %%ymm8, %%ymm2, %%ymm2\n\t"
"vpmulld %%ymm9, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm2, %%ymm3, %%ymm10\n\t"
"vpaddd 96(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 96(%[mom])\n\t"
"vpaddd %%ymm4, %%ymm5, %%ymm10\n\t"
"vpaddd 128(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 128(%[mom])\n\t"
"vpmulld %%ymm8, %%ymm4, %%ymm4\n\t"
"vpmulld %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm4, %%ymm5, %%ymm10\n\t"
"vpaddd 160(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 160(%[mom])\n\t"
"vpaddd %%ymm6, %%ymm7, %%ymm10\n\t"
"vpaddd 192(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 192(%[mom])\n\t"
"vpmulld %%ymm8, %%ymm6, %%ymm6\n\t"
"vpmulld %%ymm9, %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm6, %%ymm7, %%ymm10\n\t"
"vpaddd 224(%[mom]), %%ymm10, %%ymm10\n\t"
"vmovdqa %%ymm10, 224(%[mom])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [k] "+&r"(k)
: [out] "r"(out), [sums] "r"(sums), [mom] "r"(mom), [weight] "r"(weight)
: "memory","cc","ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
static __attribute__((always_inline)) inline void kernel_pair4_nomom(const uint8_t*a,const uint8_t*b,uint32_t*out,const int*sums){
int k=256;
asm volatile(
F2_NOMOM_CORE
"vmovntdq %%ymm0, 0(%[out])\n\t"
"vmovntdq %%ymm1, 32(%[out])\n\t"
"vmovntdq %%ymm2, 24576(%[out])\n\t"
"vmovntdq %%ymm3, 24608(%[out])\n\t"
"vmovntdq %%ymm4, 49152(%[out])\n\t"
"vmovntdq %%ymm5, 49184(%[out])\n\t"
"vmovntdq %%ymm6, 73728(%[out])\n\t"
"vmovntdq %%ymm7, 73760(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [k] "+&r"(k)
: [out] "r"(out), [sums] "r"(sums)
: "memory","cc","ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
static __attribute__((noinline)) void kernel_pair2s(const uint8_t *a,const uint8_t *b,uint32_t *out,const int *sums) {
ptrdiff_t stride=24576;
int k=256;
asm volatile (
"vmovdqa (%[sums]), %%ymm0\n\t"
"vmovdqa 32(%[sums]), %%ymm1\n\t"
"vmovdqa %%ymm0, %%ymm2\n\t"
"vmovdqa %%ymm1, %%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");
}
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 ones8=_mm256_set1_epi8(1),ones16=_mm256_set1_epi16(1);
const __m256i w1=_mm256_set1_epi32(0x03020100),w2=_mm256_set1_epi32(0x09040100);
const __m256i lane=_mm256_setr_epi32(0,4,8,12,16,20,24,28),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,d0=c0,d1=c0,e0=c0;
 for(int j=n-32;j>=0;j-=32){
  const __m256i b=_mm256_loadu_si256((const __m256i*)(br+j));
  const __m256i x0=_mm256_madd_epi16(_mm256_maddubs_epi16(ones8,b),ones16);
  const __m256i x1=_mm256_madd_epi16(_mm256_maddubs_epi16(w1,b),ones16);
  const __m256i x2=_mm256_madd_epi16(_mm256_maddubs_epi16(w2,b),ones16);
  c0=_mm256_add_epi32(c0,x0);c1=_mm256_add_epi32(c1,c0);c2=_mm256_add_epi32(c2,c1);
  d0=_mm256_add_epi32(d0,x1);d1=_mm256_add_epi32(d1,d0);e0=_mm256_add_epi32(e0,x2);
 }
 const __m256i q1=_mm256_sub_epi32(c1,c0);
 const __m256i q2=_mm256_add_epi32(_mm256_sub_epi32(_mm256_slli_epi32(c2,1),_mm256_add_epi32(_mm256_slli_epi32(c1,1),c1)),c0);
 const __m256i m1=_mm256_add_epi32(_mm256_slli_epi32(q1,5),_mm256_add_epi32(_mm256_mullo_epi32(lane,c0),d0));
 const __m256i cross=_mm256_add_epi32(_mm256_slli_epi32(_mm256_add_epi32(_mm256_mullo_epi32(lane,q1),_mm256_sub_epi32(d1,d0)),6),_mm256_add_epi32(_mm256_mullo_epi32(lane2,c0),_mm256_add_epi32(_mm256_slli_epi32(_mm256_mullo_epi32(lane,d0),1),e0)));
 __m256i qlo=_mm256_cvtepi32_epi64(_mm256_castsi256_si128(q2)),qhi=_mm256_cvtepi32_epi64(_mm256_extracti128_si256(q2,1));
 __m256i clo=_mm256_cvtepi32_epi64(_mm256_castsi256_si128(cross)),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(qlo,10),clo),_mm256_add_epi64(_mm256_slli_epi64(qhi,10),chi));
}
}

template<bool NarrowR0>
static __attribute__((always_inline)) inline void moment_errors(const uint8_t *ar,const int16_t*r0l,const int16_t*r0h,const int16_t*r1l,const int16_t*r1h,long long RS0,long long RS1,const uint32_t *observed,uint32_t&error0,uint32_t&error1){
constexpr int n=1024;
const __m256i zero=_mm256_setzero_si256();
__m256i x0=_mm256_sub_epi32(zero,_mm256_load_si256((const __m256i*)observed)),x1=zero;
__m256i x2=_mm256_sub_epi32(zero,_mm256_load_si256((const __m256i*)(observed+8))),x3=zero;
for(int k=0;k<n;k+=16){
const __m256i av=_mm256_cvtepu8_epi16(_mm_loadu_si128((const __m128i*)(ar+k)));
x0=_mm256_add_epi32(x0,_mm256_madd_epi16(av,_mm256_load_si256((const __m256i*)(r0l+k))));
if constexpr(!NarrowR0) 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))));
}
__m256i sum0=x0;
if constexpr(!NarrowR0)sum0=_mm256_add_epi32(sum0,_mm256_slli_epi32(x1,16));
const __m256i sum1=_mm256_add_epi32(x2,_mm256_slli_epi32(x3,16));
error0=(uint32_t)hsum32(sum0)-(uint32_t)(80LL*RS0);
error1=(uint32_t)hsum32(sum1)-(uint32_t)(80LL*RS1);
}
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();
static int16_t r0l[n],r0h[n],r1l[n],r1h[n];
long long RS0=0,RS1=0; bool narrow_r0=true;
for(int k=0;k<n;k++){
int v=r0[k],lo=(short)(v&0xffff); narrow_r0 &= (v==lo); 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*ar=(const uint8_t*)a.x+(size_t)i*a.ld;
uint32_t err0,err1;
if(narrow_r0) moment_errors<true>(ar,r0l,r0h,r1l,r1h,RS0,RS1,observed_moments[i],err0,err1);
else moment_errors<false>(ar,r0l,r0h,r1l,r1h,RS0,RS1,observed_moments[i],err0,err1);
uint32_t*cr=c+(size_t)i*6144+(lr-(lr>4))*16;
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)*96);
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)*96+(j%16);*pc=(uint32_t)((int32_t)*pc+d0);continue;
}
}
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)*96),bsum+j);
++i; 
} else {
for(int j=0;j<n;j+=16){
_mm256_store_si256((__m256i*)(cr+(j/16)*96),_mm256_setzero_si256());
_mm256_store_si256((__m256i*)(cr+(j/16)*96+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)*96+(j%16));
_mm256_store_si256(pc,_mm256_add_epi32(_mm256_load_si256(pc),_mm256_mullo_epi32(av,bv)));
}
}
}
}
}
static __attribute__((always_inline)) inline void kernel_pair4_final(const uint8_t*a,const uint8_t*b,uint32_t*out,const int*sums,float*dst){
const float final_scale=1.0f/(QA*QB);
int k=256;
asm volatile(
F2_NOMOM_CORE
"vbroadcastss %[scale], %%ymm15\n\t"
"vmovdqa -256(%[out]), %%ymm8\n\t"
"vmovdqa -192(%[out]), %%ymm9\n\t"
"vmovdqa -128(%[out]), %%ymm10\n\t"
"vmovdqa -64(%[out]), %%ymm11\n\t"
"vmovdqa 0(%[out]), %%ymm12\n\t"
"vmovdqa 64(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm0, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm0, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 0(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 4096(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8388608(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8392704(%[dst])\n\t"
"vmovdqa -224(%[out]), %%ymm8\n\t"
"vmovdqa -160(%[out]), %%ymm9\n\t"
"vmovdqa -96(%[out]), %%ymm10\n\t"
"vmovdqa -32(%[out]), %%ymm11\n\t"
"vmovdqa 32(%[out]), %%ymm12\n\t"
"vmovdqa 96(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm1, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm1, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 32(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 4128(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8388640(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8392736(%[dst])\n\t"
"vmovdqa 24320(%[out]), %%ymm8\n\t"
"vmovdqa 24384(%[out]), %%ymm9\n\t"
"vmovdqa 24448(%[out]), %%ymm10\n\t"
"vmovdqa 24512(%[out]), %%ymm11\n\t"
"vmovdqa 24576(%[out]), %%ymm12\n\t"
"vmovdqa 24640(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm2, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm2, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 8192(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 12288(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8396800(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8400896(%[dst])\n\t"
"vmovdqa 24352(%[out]), %%ymm8\n\t"
"vmovdqa 24416(%[out]), %%ymm9\n\t"
"vmovdqa 24480(%[out]), %%ymm10\n\t"
"vmovdqa 24544(%[out]), %%ymm11\n\t"
"vmovdqa 24608(%[out]), %%ymm12\n\t"
"vmovdqa 24672(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm3, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm3, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 8224(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 12320(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8396832(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8400928(%[dst])\n\t"
"vmovdqa 48896(%[out]), %%ymm8\n\t"
"vmovdqa 48960(%[out]), %%ymm9\n\t"
"vmovdqa 49024(%[out]), %%ymm10\n\t"
"vmovdqa 49088(%[out]), %%ymm11\n\t"
"vmovdqa 49152(%[out]), %%ymm12\n\t"
"vmovdqa 49216(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm4, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 16384(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 20480(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8404992(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8409088(%[dst])\n\t"
"vmovdqa 48928(%[out]), %%ymm8\n\t"
"vmovdqa 48992(%[out]), %%ymm9\n\t"
"vmovdqa 49056(%[out]), %%ymm10\n\t"
"vmovdqa 49120(%[out]), %%ymm11\n\t"
"vmovdqa 49184(%[out]), %%ymm12\n\t"
"vmovdqa 49248(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm5, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 16416(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 20512(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8405024(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8409120(%[dst])\n\t"
"vmovdqa 73472(%[out]), %%ymm8\n\t"
"vmovdqa 73536(%[out]), %%ymm9\n\t"
"vmovdqa 73600(%[out]), %%ymm10\n\t"
"vmovdqa 73664(%[out]), %%ymm11\n\t"
"vmovdqa 73728(%[out]), %%ymm12\n\t"
"vmovdqa 73792(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm6, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 24576(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 28672(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8413184(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8417280(%[dst])\n\t"
"vmovdqa 73504(%[out]), %%ymm8\n\t"
"vmovdqa 73568(%[out]), %%ymm9\n\t"
"vmovdqa 73632(%[out]), %%ymm10\n\t"
"vmovdqa 73696(%[out]), %%ymm11\n\t"
"vmovdqa 73760(%[out]), %%ymm12\n\t"
"vmovdqa 73824(%[out]), %%ymm13\n\t"
"vpaddd %%ymm11, %%ymm9, %%ymm14\n\t"
"vpsubd %%ymm9, %%ymm8, %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm12, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm11, %%ymm8, %%ymm8\n\t"
"vpsubd %%ymm7, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm13, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm7, %%ymm10, %%ymm10\n\t"
"vcvtdq2ps %%ymm8, %%ymm8\n\t"
"vmulps %%ymm15, %%ymm8, %%ymm8\n\t"
"vmovntps %%ymm8, 24608(%[dst])\n\t"
"vcvtdq2ps %%ymm10, %%ymm10\n\t"
"vmulps %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovntps %%ymm10, 28704(%[dst])\n\t"
"vcvtdq2ps %%ymm14, %%ymm14\n\t"
"vmulps %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntps %%ymm14, 8413216(%[dst])\n\t"
"vcvtdq2ps %%ymm9, %%ymm9\n\t"
"vmulps %%ymm15, %%ymm9, %%ymm9\n\t"
"vmovntps %%ymm9, 8417312(%[dst])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [k] "+&r"(k)
: [out] "r"(out), [sums] "r"(sums), [dst] "r"(dst), [scale] "m"(final_scale)
: "memory","cc","ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}

template<bool FINAL=false> static void leaf(Mat a,Mat b,uint32_t *c,int lr,float *dst=nullptr) {
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();
if constexpr(FINAL) {
// The final B is already packed. C-tail A remains live only until its panel.
static_assert(sizeof(observed_moments)==64*1024,"retired final panel capacity");
for(int ic=0;ic<n;) {
const int mc=ic<832?MC:(ic==832?168:(ic==1000?20:4));
const uint8_t *panel=(const uint8_t*)(a.x+(size_t)ic*a.ld);
if(ic==1020) {
memcpy(observed_moments,panel,4*1024);
panel=(const uint8_t*)observed_moments;
}
const int pf_columns=ic<624?1024:(ic==624?944:(ic==832?144:(ic==1000?192:0)));
for(int j=0;j<n;j+=NR) {
#pragma GCC unroll 1
for(int r=0;r<mc;r+=MR) {
uint32_t *out=c+(((size_t)((ic+r)*64+(j>>4))*6+4)*16);
const uint8_t *pa=panel+(size_t)r*1024;
if(j<pf_columns) _mm_prefetch((const char*)(a.x+(size_t)(ic+mc)*a.ld)+((size_t)((j>>4)*(mc/MR)+(r>>2)))*64,_MM_HINT_T1);
#pragma GCC unroll 4
for(int oldrow=0;oldrow<4;oldrow++) {
const char *old=(const char*)out+(size_t)oldrow*6144*4;
_mm_prefetch(old-256,_MM_HINT_T0);
_mm_prefetch(old-192,_MM_HINT_T0);
_mm_prefetch(old-128,_MM_HINT_T0);
_mm_prefetch(old-64,_MM_HINT_T0);
_mm_prefetch(old,_MM_HINT_T0);
_mm_prefetch(old+64,_MM_HINT_T0);
}
kernel_pair4_final(pa,bp+j*kc,out,bsum+j,dst+(size_t)(ic+r)*2048+j);
}
}
ic+=mc;
}
}else{
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;
if(need_mom){
for(int j=0;j<n;j+=NR) {
#pragma GCC unroll 1
for(int r=0;r<mc;r+=MR) {
uint32_t *out=c+(((size_t)((ic+r)*64+(j>>4))*6+(size_t)(lr-(lr>4)))*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);
{
kernel_pair4(pa,bp+j*kc,out,bsum+j,observed_moments[ic+r],moment_weights+j);
}
}
}
}else{
for(int j=0;j<n;j+=NR) {
#pragma GCC unroll 1
for(int r=0;r<mc;r+=MR) {
uint32_t *out=c+(((size_t)((ic+r)*64+(j>>4))*6+(size_t)(lr-(lr>4)))*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);
{
kernel_pair4_nomom(pa,bp+j*kc,out,bsum+j);
}
}
}
}
}
}
}
if constexpr(!FINAL) repair_leaf(a,b,c,lr);
}
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);
}
template<int B> static inline __m256i quant16b(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));
__m256i v=_mm256_packs_epi32(x,y);
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));
const __m256i byte_order=_mm256_setr_epi32(0,4,1,5,2,6,3,7);
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);
const __m256i b90=_mm256_set1_epi16(80), b45=_mm256_set1_epi16(40);
#define PK(A,B)   _mm256_permutevar8x32_epi32(ISB?_mm256_packs_epi16(A,B):_mm256_packus_epi16(A,B),byte_order)
#define PK9(A,B)  _mm256_permutevar8x32_epi32(ISB?_mm256_packs_epi16(A,B):_mm256_packus_epi16(_mm256_add_epi16(A,b90),_mm256_add_epi16(B,b90)),byte_order)
#define PK45(A,B) _mm256_permutevar8x32_epi32(ISB?_mm256_packs_epi16(A,B):_mm256_packus_epi16(_mm256_add_epi16(A,b45),_mm256_add_epi16(B,b45)),byte_order)
_mm256_stream_si256((__m256i*)(dst+t),     PK(v0a,v0b));
_mm256_stream_si256((__m256i*)(dst+z+t),   PK(v1a,v1b));
_mm256_stream_si256((__m256i*)(dst+(ISB?2:15)*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));
}
}

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;
static const int order[6]={0,1,2,3,5,6};
for(int k=0;k<6;k++) {
int r=order[k];
Mat x{aq+(size_t)ai[r]*z,h},y{bq+(size_t)bi[r]*z,h};
leaf(x,y,m,r);
}
// Final A was produced beyond the other quantized inputs. Its tail location
// is retired monotonically only after each current panel has been copied.
int8_t *tail_a=(int8_t*)c+15*z;
leaf<true>(Mat{tail_a,h},Mat{bq+6*z,h},m,4,c);
_mm_sfence();
}

评测结果