代码

/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126423: retained the accepted once-per-panel two-source dispatch and intrinsic B-panel construction, with all inherited notices below.
Idea: Also resolve one-term operands before the B-packing loops. Their intrinsic builder reads each quadrant once, applies the optional whole-operand negation, and forms the same seven packed products. The successful two-term paths and complete generic fallback remain intact.
Purpose: Expose the pure-copy operand shape to the same scheduling freedom that improved two-term packing, without changing panel geometry or any assembly.
*/
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126417: retained the accepted signed-pair B packing and all inherited source notices.
Idea: Extend the once-per-panel dispatch to two positive operand terms. Their known two-source sum now uses the same ordinary intrinsic construction as the successful signed-pair path, keeping term count and sign checks entirely outside the packing loops. Negative and one-term paths retain their established behavior.
Purpose: Remove the remaining positive two-term packing's generic control and load/store scheduling barriers, with identical seven-panel layout and modulo-64-bit values.
*/
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126042: retained the accepted long-integer engine, static column drivers, panel layouts and all inherited notices below.
Idea: For a two-term B operand containing a negative term, choose the source order and add/subtract operation once before panel construction. Form each quadrant directly from its two sources, then emit the same seven shuffled B panels with ordinary AVX2 intrinsics. The double-negative case keeps an explicit modulo-64-bit negation; all other operand shapes retain their existing path.
Purpose: Remove repeated generic term-loop, sign tests and zero accumulation in signed two-term B packing, preserving every matrix product and existing assembly region.
*/
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/123578: copied accepted indexed long-integer microkernel, packing/Strassen engine and the documented fixed-row/fixed-operation specialization idea. All inherited notices retained.
No separate software license notice is shown by these public sources.
Idea: Implement dedicated single-destination 64-product column drivers: select rows4/6 and output store/add/sub once outside all16 column microtiles. Their assembly has exactly the required output registers, no zero-count or runtime row dispatch, and no executed padding before the loop. The original first-product initialization, indexed multiply schedule, redundant lookahead broadcast, B-load pipeline and full modulo-64-bit arithmetic are preserved. Other call shapes retain the generic kernel.
Purpose: Test whether removing repeated dispatch and epilogue branches yields a real official improvement.
*/
// ===== [ml1g_] mmml1k 席 ===== **刀 K-A：mikro 热环的「入口对齐 + 不执行填充」**
//
// ===== REFERENCES（依 RULES.md §3；引用他人代码列原作者账号 + 原提交 URL + 所用内容）=====
// [1] duck.ac 用户 **saffah_cc_v41_260924**，提交 **#122376**
//     <https://duck.ac/submission/122376>（Accepted，106.718334 ms = 本题榜面 T）——
//     本件正文的**最底层基座**（逐字节复制；Strassen 递归、mikro 汇编环、面板布局、
//     buildA_blk/buildB/cmb2/cmb3/cmb*_blk 全部保留）。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#123056**（106.205045 ms =
//     我方现役 mine）<https://duck.ac/submission/123056> —— 本件正文的**直接基座**（逐字节
//     复制）。相对 [1] 它只加了本账号 ml1d_ 席的「冷写 RFO」刀（`matop` n>=512 与 `cmb3`
//     dst>=256^2 用 `vmovntdq`，见 ml1d_ 席 notes 节）。**本件不动那一刀的任何一个字节。**
// [3] 仪器（本账号自有，非引用）：`tools/exact.py`（判题机 ns 级真值/T 与严支）、
//     `tools/probe.py`（立即体验通道 = 判题机 in-process 整进程时间）、`ref/gcc9/g9.sh`
//     （判题同款 gcc 9.3，`-O2 -std=c++17 -static -U_FORTIFY_SOURCE`）。
//
// ===== 新思路（ml1g_ 席）=====
// 【真靶】开棒现读 `tools/exact.py mmml1k`：`mine #123056 = 106.205045` · `T #122376 =
//   106.718334` · 严支 **0.99·T+1µs = 105.652151** ⇒ 缺 **0.552894 ms（0.521%）· RED**。
//
// 【本刀的唯一改动】`mikro` 的汇编里，在「已剥出的第一次迭代」与热环标号 `1:` 之间原本
//   有 **56 字节 `nop` 填充**（作者在原始引擎里用来把热环对齐）。它在**落穿路径上**：
//   每次调用（n=1024 路径上 422 576 次）都要发射这 7 条 nop。
//   本刀把它换成 **`jmp 1f` + `.p2align 6`** 两个指令级的改动：
//     ① `jmp 1f`  → 落穿不再经过填充（7 条 nop uop → 1 条 jmp uop）；
//     ② `.p2align 6` → 由汇编器把热环起点对齐到 **64 字节边界**，不再依赖「函数起始
//        地址 mod 64」这个**会被任何 .text 变更平移的量**（本席实测：现役件里 mikro 起始
//        mod64 = 0、热环起点 mod64 = **3**；本刀后热环起点 mod64 = **0**，且与函数落点无关）。
//   **语义零变化**：`jmp 1f`/`.p2align` 不触碰任何寄存器、内存或标志位；热环体、A/B 面板
//   读图案、所有存指令一字未动。**逐位等价闸门**：判题同款 gcc 9.3，n = 3,4,5,7,8,9,11,
//   15,16,17,23,31,32,33,47,63,64,65,96,127,128,129,192,255,256,257,384,512,700,1024
//   共 30 档，FNV 与 SUM **全档逐位同基座** ✓（n=1024 的 FNV = 10686382741757082234 与
//   在册 ck 一致）。
//
// 【判据 / 定价】判题机 in-process 三臂篮（`work/ml1g__tmp/ml1g_arms3s.cpp`，驱动钉在引擎
//   之前；同二进制内三臂逐位同 ck）读到本刀族 −0.41 ~ −0.54 ms；但本席另用**对齐受控**的
//   三臂篮（`ml1g_arms6s.cpp`，所有臂的 mikro 都 `aligned(64)` + 环前 `.p2align 6`，实测
//   四份 mikro 的热环**全部落 mod64 = 0**、环长 201 B）复测，**同一批内 A0(现役对齐+填充)
//   = 106.887 ms** ⇒ 本刀（=A0 的代码 + jmp + .p2align 6）的板面价**必须用板面判**（探针
//   装置对「改 .text 体积」的候选有 ~0.9% 的布局敏感性，本席两篮读数符号相反 ⇒ 不可用）。
//   ⇒ **本件是「板面定价件」**：配合同批锚 `ml1g_anchor.cpp`（= #123056 逐字节 + `MR 6`
//   → `(6)` 的词法 no-op，机码与 #123056 逐字节同）分离刀价与窗口。
//
// 【实验性提交声明】本件为**实验性提交**：目的是在**判题机板面**（ns 级、已证窗口 ±0.06 ms）
//   上给「热环入口对齐 + 去掉执行填充」这一个单点定号。若板面 ≤ 现 mine 则本件即新 mine；
//   若 > 现 mine 则本件不进榜面最快件，`mine` 不受影响。
// ===== [ml1d_] mmml1k ===== **刀：冷写的 RFO —— 按作用域律只铺 n≥512 级（非临时存 `vmovntdq`）** =====
//
// ===== REFERENCES（依 RULES.md §3；引用他人代码列原作者账号 + 原提交 URL + 所用内容）=====
// [1] duck.ac 用户 **saffah_cc_v41_260924**，提交 **#122376**
//     <https://duck.ac/submission/122376>（Accepted，106.718334 ms = 本题榜面 T）——
//     **本件正文的基座**（逐字节复制；其正文里的一切注释、宏、面板布局、Strassen 递归、
//     mikro 汇编环、buildA_blk/buildB/cmb2/cmb3/cmb*_blk 全部保留）。其自述刀：
//     `l1k6_rv swp2`（B 面板存预交换形态、环内两条相隔 4 B 的载入替掉 `vpshufd`）、
//     `l1k15_rv fus3`（叶相位三目标融合）、`l1k7_rv rba`（A 面板按 MR 行分块重建）。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#122780**（106.625010 ms = 我方
//     现役 mine；机器码与 [1] 逐字节同）—— 本件的**对照件 / 同批锚**（见 §闸门）。
// [3] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），**mmml2k** 提交 **#123023**
//     <https://duck.ac/submission/123023>（749.966163 ms，−4.94 ms）—— **刀法来源**（本账号
//     自有的另一题，非逐字复制）：该件给出「**冷写的 RFO 是本族唯一被板面兑现的写侧杠杆**」
//     与**作用域判据**（该级平面在消费者读到它之前是否已被逐出；实测 8 MB/2 MB 级付钱、
//     512 KB 级 `vmovntdq` 反而 +6.5 ms）。**本件把该判据按 1k 的尺寸逐级重算后只铺在付钱的那一级。**
// [4] 仪器（本账号自有，非引用）：`tools/exact.py`（判题机 ns 级真值）、`tools/probe.py`
//     （立即体验通道 = 判题机 in-process 整进程时间）、本席生成器
//     `work/ml1d__tmp/{ml1d_gen.py,ml1d_gen2.py,ml1d_mksub.py}`；闸门 = 判题同款 gcc 9.3
//     `-O2 -std=c++17`。**本件的每一支臂的实测读数出自本席本轮的探针（见 §定价）。**
// ===== 新思路 =====（ml1d_ 席的工作）
// 【真靶】开棒现读 `tools/exact.py mmml1k` + `tools/subhist.py mmml1k` + 榜面逐行复核：
//   `mine #122780 = 106.625010` · `T #122376 = 106.718334` · 严支 **0.99·T+1µs = 105.652151**
//   ⇒ 缺 **0.972859 ms（0.912%）**。（派单里的 108.884787/107.483462/1.401325 是 **ml1y7_ 时期
//   的两代前口径**，与判题机现读不符；以现读为准。）
// 【判据：谁的钱】判题机 = i3-8100，L1d 32 KB / L2 256 KB / **L3 6 MB**。NT 只对
//   「**写后不再原样读回、且写出前已被逐出**」的块付钱（本族在册判词：`#104960` NT 写 C = −0.27% ✓、
//   `ntadd` NT 写 M 目的地 = +0.53% ✗，差别就在消费者还读不读这份数据）。
//   本席用计数器把本题**全部写侧账量清**（`work/ml1d__tmp/ml1d_cnt4.cpp`，n=1024 实跑）：
//     `matop`   80 呼 / **55.000 MB 纯存**（n=512 十呼 = 20 MB；n=256 七十呼 = 35 MB）
//     `cmb3`  2736 呼 / 16.531 MB 纯存，其中 dst ≥ 256×256 只有 **8.875 MB**（其余 2695 呼是叶级热瓦片）
//     `cmb2`   162 呼 / 63.000 MB、`cmb2_blk` 10780 呼 / 61.250 MB —— **全是 RMW（先读后写，
//              行已独占、本来就无 RFO）⇒ NT 在结构上不可用，不在本刀域内**
// 【作用域律的 1k 版（逐级数出来的）】
//   * n=**512** 的 `matop` dst = **2 MB**：写完还要先写另一侧 2 MB 才轮到消费者，之后 7 个 256-子节点
//     各自再物化 1 MB + 各跑 7 个叶（buildB 每叶 229 KB 写）⇒ 写—读间隔 ≈ 6~7 MB > L3 6 MB
//     ⇒ **逐出成立 = 冷** ⇒ 每写白付一次 RFO（20 MB）。
//   * n=**256** 的 `matop` dst = **512 KB**：其 7 个 128-孙节点的 L3 足迹 ≈ 512 KB + 512 KB
//     + 230 KB（B 面板回收复用）≈ **1.3 MB ≪ 6 MB** ⇒ **驻留 = 热**（消费者 122 MB 读全在 L3）
//     ⇒ **NT 会把那 122 MB 打到 DRAM**。**这是本刀唯一的域边界，本席用两支探针臂把它量成了数。**
// 【本件改动】（只改**存指令形态**，值/次序/字节数/面板布局一字未动）
//   ① `matop2`：`n >= 512` 时用 `vmovntdq`（消 20 MB 写-分配 RFO），`n == 256` 保持 `vmovdqu`；
//   ② `cmb3`：dst 行数 ≥ 256 且 dst 32 B 对齐时用 `vmovntdq`（再消 8.875 MB）。
//   `vmovntdq` 需 32 B 对齐：`matop` 的 dst = `matbuf[0/1]`（64 B 对齐，行距 n*8 = 2048/4096 B）✓；
//   `cmb3` 在**入口**做一次性对齐 + `ldd` 检查，不满足的调用原样走可缓存存 ✓。单线程同核
//   「同地址 load-after-store」由 x86 内存模型保证（2k 的 #123023 已在 n=2048 全元素验证过同一条）。
// 【定价（判题机 in-process，本席本轮两支探针，A0=基座同二进制同机码）】
//   探针#1：A0 107.302 · A1(`matop` n≥512 NT) **106.826 → −0.476** · A2(再加 n=256 NT) 108.925
//           → **+1.623** · A3(A1+cmb3≥256 NT) **106.629 → −0.673**
//   探针#2：A0 106.985 · A1 **106.482 → −0.503** · A2(+cmb3≥256) 106.486 → −0.499 ·
//           A3(A2+cmb3≥128+目的地 `prefetchw`) 106.713 → −0.272
//   ⇒ ★ **两发同号：`matop` n≥512 的 NT = −0.48/−0.50 ms（本刀的全部确定性收益）**；
//     `cmb3` ≥256 = −0.197/+0.004（噪声内，取了低风险的一半）；**n=256 级 NT = +2.10 ms（判负，
//     与 mmml2k 的 512 KB 级 +6.5 ms 同形）**；目的地 `prefetchw` = +0.23（判负，未采纳）。
//   四臂 FNV 校验和逐位相同（`10686382741757082234`）✓。
// 【诚实交底】本刀只买到 ≈ **−0.50 ms**，而缺口是 0.972859 ms ⇒ **本件发射后预期仍 RED（约缺 0.47 ms）**，
//   发射的目的是**把已兑现的收益落袋 + 给下一手留下「哪一级付钱」的判题机实测边界**，
//   不是为了冲线。**未采纳项（n=256 NT / prefetchw）已判负，别再重走。**
// 【闸门】判题同款 gcc 9.3、`-O2 -std=c++17`：n = 128/256/512/1024 四档 FNV 与 [1]/[2] 逐位相同 ✓
// ===== [ml1b_] mmml1k ===== **同批基座锚（诊断件，非候选）** =====
//
// ===== REFERENCES（依 RULES.md §3；引用他人代码列原作者账号 + 原提交 URL + 所用内容）=====
// [1] duck.ac 用户 **saffah_cc_v41_260924**，提交 **#122376**
//     <https://duck.ac/submission/122376>（Accepted，106.718334 ms = 本题榜面 T）——
//     **本件正文的基座**（逐字节复制 + 本头部与 §思路 所述改动；其正文里的一切注释、宏、
//     面板布局、Strassen 递归、mikro 汇编环、buildA_blk/buildB/cmb* 全部保留）。
//     其自述刀：`l1k6_rv swp2`（B 面板存预交换形态、环内两条相隔 4 B 的载入替掉 `vpshufd`）、
//     `l1k15_rv fus3`（叶相位三目标融合）、`l1k7_rv rba`（A 面板按 MR 行分块重建）。
// [2] duck.ac 本账号 **saffah_cc_v41_agg1**，提交 **#122723**（106.759485 ms = 我方现役 mine，
//     机器码与 [1] 逐字节同，`.text` md5 `a420b4f863263b3d`）—— 本件的**对照件**（闸门 + 探针标定）。
// [3] 仪器（本账号自有，非引用）：`tools/exact.py`（判题机 ns 级真值）、`tools/probe.py`
//     （立即体验通道 = 判题机 in-process 整进程 Time (ms)）、本席生成器
//     `work/ml1b__tmp/gen1.py`（从 [1] 的正文机器化生成本件）、闸门 = 判题同款 gcc 9.3 `-O2`。
// ===== 思路 =====（ml1b_ 席的新工作）
// 【真靶】T = #122376 = 106.718334 ⇒ **真靶 = 0.99·T+1µs = 105.652151**；mine = 106.759485
//   ⇒ 缺 **1.107334 ms（1.037%）**。派单方向 =「64 级叶路径专用化 / 去通用代码」（mmml4k 的
//   对手链实测真价：定长打包器 + store-only 专用驱动器买 −29.9/−23.4 ms）。
// 【判题机相位表（ml1a_ 装置，基座口径）】引擎 106.612 = mikro 81.37（76.3%）· matop 8.41 ·
//   buildB 5.09 · cmb2 4.85 · cmb3 1.46 · cmb2_blk 1.05 · cmb3_blk 0.30 · 其余 4.1。
// 【本刀的单点性】mikro 是 `static inline`，其**汇编环逐字节未动**（34 指令 / 39 uop / 24 MAC，
//   环体字节数与指令次序完全一致），改的只有**每次调用的调用方/收尾侧**：
//   ① 形参 10 个 → 5 个：原形参表 `(cnt,ap,bp,c1,ldc1,op1,c2,ldc2,op2,rows)` 里 **后 4 个走栈**
//      （第 7~10 参）⇒ 每次调用都有 4 处栈写 + 栈读（`movl 16(%rbp),%r10d` 等）；
//   ② 尾声的**逐行守卫链**：`if (rows > 0) … if (rows > 5)` = 6 组 `cmpl/je`；
//   ③ 尾声的**逐行 c2 判空**：`testq %r9,%r9 / je` × 6；
//   ④ 序言的 `cmp $4, %[rows] / je` 形状分派。
//   现改为 `template<int ROWS,int OP> mikro64s(cnt,ap,bp,c1,ldc1)`：**ROWS∈{6,4}、OP∈{0,1,2}
//   全部编译期**，尾声是 6（或 4）段直落代码（无分支、无守卫、无栈参）。
//   **调用方 = leaflevel3 → leafP_blk，其 7 个调用点全部满足 `C2==0 ∧ n==KC==64 ∧ rows∈{6,4}
//   ∧ op1∈{0,1,2}`**（L=0,1,2 存储；L=5,6 加；L=3,4 先落 T 再 cmb2_blk，也是 C2=0）；不满足该
//   条件的调用（其它 n、双目标）**原样走旧通用 mikro**，行为逐位不变。
// 【为什么这不是「环内减 uop」】环体一字未动 ⇒ 不触碰 ml4j_/ml4k_ 判死的「环体字节数不增」判据，
//   也不动端口配比；改掉的是**每 64 k 步一次**的固定开销（422 576 次/全程）。
// 【定价（判题机 in-process，同批；probe 双胞胎臂噪声 0.0086%）】基座 **112.475784 / 112.522063**
//   ⇒ 本件 **111.699004 / 111.805461** = **−0.7768 / −0.7166 ms**（两发同号，均值 ≈ **−0.75 ms**）；
//   两件 stdout 校验和逐位同（`ck=10686382741757082234`）✓。按在册「减指令族 1:1 传导」
//   （a14ah_：探针 −1.63% ↔ 板面 −1.54%）预期板面同量级。
// 【闸门】判题同款 gcc 9.3、`-O2 -std=c++17`：**30 档** n = 3,5,6,7,9,10,11,13,15,17,20,24,27,
//   33,39,40,64,100,128,129,130,200,256,260,400,511,512,700,1023,1024 的 FNV/SUM 与 [1]/[2]
//   **逐位全同** ✓（含全部 n≥2·CUTL 的叶路径档）。`.text`：86 246 B → 87 974 B（+1.7 KB）。
// 【未动】mikro 汇编环、面板布局与尺寸、`matop`/`buildB`/`buildA_blk`/`cmb2`/`cmb3`/`cmb*_blk`、
//   全部宏常量（MR/KC/MC/NC/UU/CUTL/KCP）、Strassen 递归结构、C 输出寻址 ✓
// ============================================================================================
#pragma GCC optimize("O3","unroll-loops")
/* mmml: C = A*B mod 2^64, long long, FUNCTION-style C++ linkage.
5-ALU-uop mikro (vpmuludq + vpaddq + vpmulld + vpaddd per k x row x 4 cols, no p5 uop)
l1k6_rv `swp2`: the B panel is STORED PRE-SWAPPED and the mikro reads BOTH operands out
of it with two 32-byte loads 4 bytes apart -- the p5 `vpshufd` is deleted, one L1 load
replaces it, no second panel and no extra memory.
plus Strassen with fused operand combination and dual-destination stores
(no extra memory).  MR=6 KC=512 MC=48 NC=512 CUTL=512 U=1 BSW=0 EXP=0 */
#include <immintrin.h>

/* Single-instruction 256-bit unaligned access, NO alignment precondition.
g++-9 -O2 with no -march expands the unaligned 256-bit intrinsics into 2-3 uops, one
of them a p5-only vinsertf128/vextractf128; these emit the single instruction the
intrinsic is named for.  AT&T order is `vmovdqu src,dst`: the STORE is %0,%1 (ymm
first, memory second).  A REVERSED store operand silently emits a LOAD instead --
that is how this was first written and it produced wrong output with no diagnostic.
volatile + "memory" is REQUIRED: a non-volatile asm with only an "m" input is pure
and gets eliminated. */
#define LDQ(p) ({ __m256i _ldq_v; __asm__("vmovdqu %1,%0" : "=x"(_ldq_v) : "m"(*(const __m256i *)(p)) : "memory"); _ldq_v; })
#define STQ(p, v) __asm__ volatile("vmovdqu %0,%1" :: "x"(v), "m"(*(__m256i *)(p)) : "memory")
/* SWQ == `vpshufd $0xB1`: exchange the two 32-bit halves inside every 64-bit lane.
   The B panel is STORED pre-swapped so the mikro can get both of the operands it needs
   out of the SAME panel with two loads 4 bytes apart (see the mikro). */
#define SWQ(v) _mm256_shuffle_epi32((v), 0xB1)
#include <stddef.h>
#pragma GCC push_options
#pragma GCC target("avx2")

#define STQNT(p, v) __asm__ volatile("vmovntdq %0,%1" :: "x"(v), "m"(*(__m256i *)(p)) : "memory")
static void cmb3(long long *d, int ldd, const long long *a, const long long *b, const long long *c, int ld, int rows, int ncol);
typedef long long TE;

#define MR (6)
#define KC 64
#define MC ( 12 )
#define NC 512
#define UU 1
#define CUTL 64

/* A-panel row stride PADDED by 8 uint64 (64 B).  With the natural stride KC*8 = 4096 B
every one of the 6 rows of a micro-tile maps to the SAME L1 set (4096/64 = 64 sets
apart, 64 sets in the cache) -- a textbook conflict.  +8 uint64 makes the stride
4160 B = 65 lines, so consecutive rows land in consecutive sets. */
#define KCP (KC + 12)
static TE Apanel[(size_t)MC * KCP + 64] __attribute__((aligned(64)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(64)));

/* op: 0 = store, 1 = add, 2 = sub;  c2 may be 0 (single destination) */
static inline void mikro(long cnt, const long long *ap, const long long *bp,
long long *c1, int ldc1, int op1,
long long *c2, int ldc2, int op2, int rows) {
__m256i AC[12];
long idx = -(cnt << 3);
long long *apE = (long long *)ap + cnt;
long long *bpE = (long long *)bp + (cnt << 2);
    __asm__ volatile(
    "test %[idx], %[idx]\n\t"
    "jz 8f\n\t"
    "vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu 4(%[bp],%[idx],4), %%ymm13\n\t"
    "cmp $4, %[rows]\n\t"
    "je 6f\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c0]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c1]\n\t"
    "vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c2]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c3]\n\t"
    "vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c4]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c5]\n\t"
    "vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c6]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c7]\n\t"
    "vbroadcastsd 2432(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c8]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c9]\n\t"
    "vbroadcastsd 3040(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c10]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c11]\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
    "add $8, %[idx]\n\t"
    "jz 2f\n\t"
    ".byte 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00\n\t"
    "1:\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c1], %[c1]\n\t"
    "vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c0], %[c0]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c3], %[c3]\n\t"
    "vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c2], %[c2]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c5], %[c5]\n\t"
    "vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c4], %[c4]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c7], %[c7]\n\t"
    "vbroadcastsd 2432(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c6], %[c6]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c9], %[c9]\n\t"
    "vbroadcastsd 3040(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c8], %[c8]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c11], %[c11]\n\t"
    "vbroadcastsd 8(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c10], %[c10]\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
    "add $8, %[idx]\n\t"
    "jnz 1b\n\t"
    "jmp 2f\n\t"
    "6:\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c0]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c1]\n\t"
    "vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c2]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c3]\n\t"
    "vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c4]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c5]\n\t"
    "vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %[c6]\n\t"
    "vpmulld %%ymm12, %%ymm14, %[c7]\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
    "add $8, %[idx]\n\t"
    "jz 2f\n\t"
    "5:\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c1], %[c1]\n\t"
    "vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c0], %[c0]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c3], %[c3]\n\t"
    "vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c2], %[c2]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c5], %[c5]\n\t"
    "vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %[c4], %[c4]\n\t"
    "vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
    "vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %[c7], %[c7]\n\t"
    "vpaddq %%ymm15, %[c6], %[c6]\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
    "add $8, %[idx]\n\t"
    "jnz 5b\n\t"
    "jmp 2f\n\t"
    "8:\n\t"
    "vpxor %[c0], %[c0], %[c0]\n\t"
    "vmovdqa %[c0], %[c1]\n\t"
    "vmovdqa %[c0], %[c2]\n\t"
    "vmovdqa %[c0], %[c3]\n\t"
    "vmovdqa %[c0], %[c4]\n\t"
    "vmovdqa %[c0], %[c5]\n\t"
    "vmovdqa %[c0], %[c6]\n\t"
    "vmovdqa %[c0], %[c7]\n\t"
    "vmovdqa %[c0], %[c8]\n\t"
    "vmovdqa %[c0], %[c9]\n\t"
    "vmovdqa %[c0], %[c10]\n\t"
    "vmovdqa %[c0], %[c11]\n\t"
    "2:\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx), [c0] "=x"(AC[0]), [c1] "=x"(AC[1]), [c2] "=x"(AC[2]), [c3] "=x"(AC[3]), [c4] "=x"(AC[4]), [c5] "=x"(AC[5]), [c6] "=x"(AC[6]), [c7] "=x"(AC[7]), [c8] "=x"(AC[8]), [c9] "=x"(AC[9]), [c10] "=x"(AC[10]), [c11] "=x"(AC[11])
    : [rows] "r"((long)rows)
    : "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory");
if (rows > 0) {
__m256i lo = AC[0], hi = AC[1];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)0 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)0 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 1) {
__m256i lo = AC[2], hi = AC[3];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)1 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)1 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 2) {
__m256i lo = AC[4], hi = AC[5];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)2 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)2 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 3) {
__m256i lo = AC[6], hi = AC[7];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)3 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)3 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 4) {
__m256i lo = AC[8], hi = AC[9];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)4 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)4 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 5) {
__m256i lo = AC[10], hi = AC[11];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)5 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)5 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
}

template<int R,int OP> static inline __attribute__((always_inline)) void mikro64_static(const long long* ap,const long long* bp,long long*c,int ldc){
__m256i AC[2*R]; long idx=-512; const long long*apE=ap+64,*bpE=bp+256;
if constexpr(R==6){
 __asm__ volatile(
  "vmovdqu (%[bp],%[idx],4), %%ymm12\n"
  "\tvmovdqu 4(%[bp],%[idx],4), %%ymm13\n"
  "\tvbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c0]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c1]\n"
  "\tvbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c2]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c3]\n"
  "\tvbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c4]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c5]\n"
  "\tvbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c6]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c7]\n"
  "\tvbroadcastsd 2432(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c8]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c9]\n"
  "\tvbroadcastsd 3040(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c10]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c11]\n"
  "\tvmovdqu 32(%[bp],%[idx],4), %%ymm12\n"
  "\tvmovdqu 36(%[bp],%[idx],4), %%ymm13\n"
  "\tadd $8, %[idx]\n"
  "\t1:\n"
  "\tvbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c1], %[c1]\n"
  "\tvbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c0], %[c0]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c3], %[c3]\n"
  "\tvbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c2], %[c2]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c5], %[c5]\n"
  "\tvbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c4], %[c4]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c7], %[c7]\n"
  "\tvbroadcastsd 2432(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c6], %[c6]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c9], %[c9]\n"
  "\tvbroadcastsd 3040(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c8], %[c8]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c11], %[c11]\n"
  "\tvbroadcastsd 8(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c10], %[c10]\n"
  "\tvmovdqu 32(%[bp],%[idx],4), %%ymm12\n"
  "\tvmovdqu 36(%[bp],%[idx],4), %%ymm13\n"
  "\tadd $8, %[idx]\n"
  "\tjnz 1b\n"
  "\t"
 : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx), [c0] "=x"(AC[0]), [c1] "=x"(AC[1]), [c2] "=x"(AC[2]), [c3] "=x"(AC[3]), [c4] "=x"(AC[4]), [c5] "=x"(AC[5]), [c6] "=x"(AC[6]), [c7] "=x"(AC[7]), [c8] "=x"(AC[8]), [c9] "=x"(AC[9]), [c10] "=x"(AC[10]), [c11] "=x"(AC[11])
 : : "ymm12","ymm13","ymm14","ymm15","cc","memory");
}
else{
 __asm__ volatile(
  "vmovdqu (%[bp],%[idx],4), %%ymm12\n"
  "\tvmovdqu 4(%[bp],%[idx],4), %%ymm13\n"
  "\t\tvbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c0]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c1]\n"
  "\tvbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c2]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c3]\n"
  "\tvbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c4]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c5]\n"
  "\tvbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %[c6]\n"
  "\tvpmulld %%ymm12, %%ymm14, %[c7]\n"
  "\tvmovdqu 32(%[bp],%[idx],4), %%ymm12\n"
  "\tvmovdqu 36(%[bp],%[idx],4), %%ymm13\n"
  "\tadd $8, %[idx]\n"
  "\t5:\n"
  "\tvbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c1], %[c1]\n"
  "\tvbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c0], %[c0]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c3], %[c3]\n"
  "\tvbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c2], %[c2]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c5], %[c5]\n"
  "\tvbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n"
  "\tvpaddq %%ymm15, %[c4], %[c4]\n"
  "\tvpmuludq %%ymm14, %%ymm13, %%ymm15\n"
  "\tvpmulld %%ymm12, %%ymm14, %%ymm14\n"
  "\tvpaddd %%ymm14, %[c7], %[c7]\n"
  "\tvpaddq %%ymm15, %[c6], %[c6]\n"
  "\tvmovdqu 32(%[bp],%[idx],4), %%ymm12\n"
  "\tvmovdqu 36(%[bp],%[idx],4), %%ymm13\n"
  "\tadd $8, %[idx]\n"
  "\tjnz 5b\n"
  "\t"
 : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx), [c0] "=x"(AC[0]), [c1] "=x"(AC[1]), [c2] "=x"(AC[2]), [c3] "=x"(AC[3]), [c4] "=x"(AC[4]), [c5] "=x"(AC[5]), [c6] "=x"(AC[6]), [c7] "=x"(AC[7])
 : : "ymm12","ymm13","ymm14","ymm15","cc","memory");
}
#pragma GCC unroll 6
for(int r=0;r<R;r++){
__m256i hi=AC[2*r+1],lo=AC[2*r];
__m256i y=_mm256_add_epi32(hi,_mm256_srli_epi64(hi,32));
__m256i z=_mm256_add_epi64(lo,_mm256_slli_epi64(y,32));
long long*d=c+(size_t)r*ldc;
if constexpr(OP==1) z=_mm256_add_epi64(LDQ((const __m256i*)d),z);
if constexpr(OP==2) z=_mm256_sub_epi64(LDQ((const __m256i*)d),z);
STQ((__m256i*)d,z);
}
}
template<int R,int OP> static inline void columns64_static(const long long*ap,const long long*pb,long long*c,int ld){
for(int j=0;j<64;j+=4) mikro64_static<R,OP>(ap,pb+(size_t)(j>>2)*64*4,c+j,ld);
}

struct Opd {                 /* value(i,j) = sum_t s[t]*p[(r[t]+i)*ld + c[t]+j]  (nt <= 4) */
const long long *p; int ld; int nt; int r[8], c[8], s[8];
};
static inline Opd oq(const Opd *o, int qi, int qj, int h) {
Opd d = *o;
for (int t = 0; t < o->nt; t++) { d.r[t] += qi * h; d.c[t] += qj * h; }
return d;
}
static inline Opd oadd(const Opd *a, const Opd *b, int sgn) {
Opd d = *a;
for (int t = 0; t < b->nt; t++) { d.r[d.nt] = b->r[t]; d.c[d.nt] = b->c[t]; d.s[d.nt] = b->s[t] < 0 ? -sgn : sgn; d.nt++; }
return d;
}
static inline Opd o1(const long long *p, int ld) {
Opd d; d.p = p; d.ld = ld; d.nt = 1; d.r[0] = 0; d.c[0] = 0; d.s[0] = 1; return d;
}

/* A panel: rows [ic, ic+mc), k in [pc, pc+kc) of the operand */
static void packA(const Opd *o, int ic, int mc, int pc, int kc) {
const int nt = o->nt;
const long long *P = o->p; const int ld = o->ld;
for (int r = 0; r < MC; r++) {
long long *p = Apanel + (size_t)r * KCP;
int k = 0;
if (r < mc) {
const long long *s0 = P + (size_t)(o->r[0] + ic + r) * ld + o->c[0] + pc;
if (nt == 1 || nt == 2) {
const long long *s1 = (nt == 2) ? P + (size_t)(o->r[1] + ic + r) * ld + o->c[1] + pc : 0;
int sgn1 = (nt == 2) ? o->s[1] : 1;
for (; k + 4 <= kc; k += 4) {
__m256i v = LDQ((const __m256i *)(s0 + k));
if (nt == 2) {
__m256i w = LDQ((const __m256i *)(s1 + k));
v = (sgn1 > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
}
STQ((__m256i *)p, v);
p += 4;
}
for (; k < kc; k++) { long long v = s0[k];
if (nt == 2) v = (sgn1 > 0) ? v + s1[k] : v - s1[k];
*p++ = (v); }
} else {
for (; k + 16 <= kc; k += 16) {
__m256i v0 = _mm256_setzero_si256(), v1 = v0, v2 = v0, v3 = v0;
for (int t = 0; t < nt; t++) {
const long long *sp = P + (size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k;
__m256i w0 = LDQ((const __m256i *)sp);
__m256i w1 = LDQ((const __m256i *)(sp + 4));
__m256i w2 = LDQ((const __m256i *)(sp + 8));
__m256i w3 = LDQ((const __m256i *)(sp + 12));
if (o->s[t] < 0) {
v0 = _mm256_sub_epi64(v0, w0); v1 = _mm256_sub_epi64(v1, w1);
v2 = _mm256_sub_epi64(v2, w2); v3 = _mm256_sub_epi64(v3, w3); }
else {
v0 = _mm256_add_epi64(v0, w0); v1 = _mm256_add_epi64(v1, w1);
v2 = _mm256_add_epi64(v2, w2); v3 = _mm256_add_epi64(v3, w3); }
}
STQ((__m256i *)p, v0);
STQ((__m256i *)(p + 4), v1);
STQ((__m256i *)(p + 8), v2);
STQ((__m256i *)(p + 12), v3);
p += 16;
}
for (; k + 4 <= kc; k += 4) {
__m256i v = _mm256_setzero_si256();
for (int t = 0; t < nt; t++) {
__m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k));
v = (o->s[t] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
}
STQ((__m256i *)p, v);
p += 4;
}
for (; k < kc; k++) {
long long v = 0;
for (int t = 0; t < nt; t++)
v += o->s[t] < 0 ? -P[(size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k]
:  P[(size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k];
*p++ = (v);
}
}
}
p += KC - k;   /* k >= kc is never read: mikro runs kc steps */
}
}

static void packB(const Opd *o, int jc, int nc, int pc, int kc) {
const int nt = o->nt;
const long long *P = o->p; const int ld = o->ld;
const int jbn = (nc + 3) & ~3;   /* mikro reads only jg < nc, step 4 */
int jb = 0;
/* PACK8: 8 columns per k.  The 4-column form loads 32 B at a stride of
ld*8 (4 KB or more), so every load is a distinct, half-used line and the
hardware prefetcher cannot see the stream.  Reading jb and jb+4 together
makes each k one FULL 64 B line, halving the line count and the DRAM
traffic.  All 8 columns come from one line because jc+jb is a multiple of
8 (jb steps by 8) and P is 4096-byte aligned.  Byte-identical: the two
stores write exactly the values the 4-column loop wrote to the same
addresses. */
if (nt <= 2) {
for (; jb + 8 <= nc; jb += 8) {
long long *q0 = Bpanel + (size_t)(jb / 4) * KC * 4;
long long *q1 = q0 + (size_t)KC * 4;
const long long *sq  = P + (size_t)o->r[0] * ld + o->c[0] + jc + jb;
const long long *sq1 = (nt >= 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc + jb : 0;
int k = 0;
for (; k + 2 <= kc; k += 2) {
const long long *a0 = sq + (size_t)(pc + k) * ld;
const long long *a1 = a0 + ld;
if (k + 3 < kc) { _mm_prefetch((const char *)(a0 + (size_t)3 * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(a1 + (size_t)3 * ld), _MM_HINT_T0); }
__m256i v0 = LDQ((const __m256i *)a0), u0 = LDQ((const __m256i *)(a0 + 4));
__m256i v1 = LDQ((const __m256i *)a1), u1 = LDQ((const __m256i *)(a1 + 4));
if (nt == 2) {
const long long *b0 = sq1 + (size_t)(pc + k) * ld;
const long long *b1 = b0 + ld;
__m256i w0 = LDQ((const __m256i *)b0), x0 = LDQ((const __m256i *)(b0 + 4));
__m256i w1 = LDQ((const __m256i *)b1), x1 = LDQ((const __m256i *)(b1 + 4));
if (o->s[1] > 0) { v0 = _mm256_add_epi64(v0, w0); u0 = _mm256_add_epi64(u0, x0);
v1 = _mm256_add_epi64(v1, w1); u1 = _mm256_add_epi64(u1, x1); }
else             { v0 = _mm256_sub_epi64(v0, w0); u0 = _mm256_sub_epi64(u0, x0);
v1 = _mm256_sub_epi64(v1, w1); u1 = _mm256_sub_epi64(u1, x1); }
}
_mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), SWQ(v0));
_mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), SWQ(u0));
_mm256_store_si256((__m256i *)(q0 + (size_t)(k + 1) * 4), SWQ(v1));
_mm256_store_si256((__m256i *)(q1 + (size_t)(k + 1) * 4), SWQ(u1));
}
for (; k < kc; k++) {
const long long *a = sq + (size_t)(pc + k) * ld;
__m256i v = LDQ((const __m256i *)a);
__m256i u = LDQ((const __m256i *)(a + 4));
if (nt == 2) {
const long long *b = sq1 + (size_t)(pc + k) * ld;
__m256i w = LDQ((const __m256i *)b);
__m256i x = LDQ((const __m256i *)(b + 4));
if (o->s[1] > 0) { v = _mm256_add_epi64(v, w); u = _mm256_add_epi64(u, x); }
else             { v = _mm256_sub_epi64(v, w); u = _mm256_sub_epi64(u, x); }
}
_mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), SWQ(v));
_mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), SWQ(u));
}
}
}
for (; jb < jbn; jb += 4) {
long long *q = Bpanel + (size_t)(jb / 4) * KC * 4;
int full = (jb + 4 <= nc);
int lim = nc - jb; if (lim > 4) lim = 4; if (lim < 0) lim = 0;
const long long *sq = P + (size_t)o->r[0] * ld + o->c[0] + jc + jb;
const long long *sq1 = (nt >= 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc + jb : 0;
if (full && nt <= 2) {
for (int k = 0; k < kc; k++) {
__m256i v = LDQ((const __m256i *)(sq + (size_t)(pc + k) * ld));
if (nt == 2) {
__m256i w = LDQ((const __m256i *)(sq1 + (size_t)(pc + k) * ld));
v = (o->s[1] > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
}
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(v));
}
} else if (full) {
int k = 0;
for (; k + 8 <= kc; k += 8) {
__m256i v0 = _mm256_setzero_si256(), v1 = v0, v2 = v0, v3 = v0;
__m256i v4 = v0, v5 = v0, v6 = v0, v7 = v0;
for (int tt = 0; tt < nt; tt++) {
const long long *sp = P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb;
__m256i w0 = LDQ((const __m256i *)sp);
__m256i w1 = LDQ((const __m256i *)(sp + ld));
__m256i w2 = LDQ((const __m256i *)(sp + 2 * ld));
__m256i w3 = LDQ((const __m256i *)(sp + 3 * ld));
__m256i w4 = LDQ((const __m256i *)(sp + 4 * ld));
__m256i w5 = LDQ((const __m256i *)(sp + 5 * ld));
__m256i w6 = LDQ((const __m256i *)(sp + 6 * ld));
__m256i w7 = LDQ((const __m256i *)(sp + 7 * ld));
if (o->s[tt] < 0) {
v0 = _mm256_sub_epi64(v0, w0); v1 = _mm256_sub_epi64(v1, w1);
v2 = _mm256_sub_epi64(v2, w2); v3 = _mm256_sub_epi64(v3, w3);
v4 = _mm256_sub_epi64(v4, w4); v5 = _mm256_sub_epi64(v5, w5);
v6 = _mm256_sub_epi64(v6, w6); v7 = _mm256_sub_epi64(v7, w7);
} else {
v0 = _mm256_add_epi64(v0, w0); v1 = _mm256_add_epi64(v1, w1);
v2 = _mm256_add_epi64(v2, w2); v3 = _mm256_add_epi64(v3, w3);
v4 = _mm256_add_epi64(v4, w4); v5 = _mm256_add_epi64(v5, w5);
v6 = _mm256_add_epi64(v6, w6); v7 = _mm256_add_epi64(v7, w7);
}
}
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(v0));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 1) * 4), SWQ(v1));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 2) * 4), SWQ(v2));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 3) * 4), SWQ(v3));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 4) * 4), SWQ(v4));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 5) * 4), SWQ(v5));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 6) * 4), SWQ(v6));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 7) * 4), SWQ(v7));
}
for (; k < kc; k++) {
__m256i v = _mm256_setzero_si256();
for (int tt = 0; tt < nt; tt++) {
__m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb));
v = (o->s[tt] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
}
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(v));
}
} else {
for (int k = 0; k < KC; k++) {
long long t[4] = {0, 0, 0, 0};
if (k < kc) {
int ncol = full ? 4 : lim;
for (int c = 0; c < ncol; c++) {
long long v = 0;
for (int tt = 0; tt < nt; tt++)
v += o->s[tt] < 0 ? -P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c]
:  P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c];
t[c] = v;
}
}
__m256i v = LDQ((const __m256i *)t);
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(v));
}
}
}
}


static inline int effop(int base, int term) {   /* term: 1 = +M, 2 = -M ; result: 0 store 1 add 2 sub */
if (base == 0) return term == 1 ? 0 : 2;
if (base == 1) return term == 2 ? 2 : 1;
return term == 2 ? 1 : 2;
}

/* C1 op1 (+- second destination C2 op2, either may be absent) = Ao * Bo, block size n */
static void leaf(int n, const Opd *Ao, const Opd *Bo,
long long *C1, int ldc1, int op1, long long *C2, int ldc2, int op2) {
int nb4 = n - (n & 3);
for (int jc = 0; jc < nb4; jc += NC) {
int nc = (nb4 - jc < NC) ? (nb4 - jc) : NC;
for (int pc = 0; pc < n; pc += KC) {
int kc = (n - pc < KC) ? (n - pc) : KC;
packB(Bo, jc, nc, pc, kc);
int e1 = (pc == 0) ? op1 : (op1 == 2 ? 2 : 1);
int e2 = (pc == 0) ? op2 : (op2 == 2 ? 2 : 1);
for (int ic = 0; ic < n; ic += MC) {
int mc = (n - ic < MC) ? (n - ic) : MC;
packA(Ao, ic, mc, pc, kc);
for (int ir = 0; ir < mc; ir += MR) {
int rows = (mc - ir < MR) ? (mc - ir) : MR;
const long long *ap = Apanel + (size_t)ir * KCP * 1;
for (int jg = 0; jg < nc; jg += 4) {
const long long *bp = Bpanel + (size_t)(jg / 4) * KC * 4;
mikro(kc / UU, ap, bp, C1 + (size_t)(ic + ir) * ldc1 + jc + jg, ldc1, e1,
C2 ? C2 + (size_t)(ic + ir) * ldc2 + jc + jg : 0, ldc2, e2, rows);
}
}
}
}
}
for (int j = nb4; j < n; j++)
for (int i = 0; i < n; i++) {
unsigned long long acc = 0;
for (int k = 0; k < n; k++) {
long long av = 0, bv = 0;
for (int t = 0; t < Ao->nt; t++)
av += Ao->s[t] < 0 ? -Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k]
:  Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k];
for (int t = 0; t < Bo->nt; t++)
bv += Bo->s[t] < 0 ? -Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j]
:  Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j];
acc += (unsigned long long)av * (unsigned long long)bv;
}
long long *cp = C1 + (size_t)i * ldc1 + j;
if (op1 == 0) *cp = (long long)acc;
else if (op1 == 1) *cp += (long long)acc;
else *cp -= (long long)acc;
if (C2) {
long long *cq = C2 + (size_t)i * ldc2 + j;
if (op2 == 0) *cq = (long long)acc;
else if (op2 == 1) *cq += (long long)acc;
else *cq -= (long long)acc;
}
}
}

/* recursive Strassen: C1 op1 op= Ao*Bo  and  C2 op2 op= Ao*Bo */
/* Scratch: one h*h product block per recursion level, carved from the tail of the
buffer.  Only the levels that actually recurse consume space, so a walk of the
recursion computes exactly the bytes matrix_multiply() hands us below. */

/* ===================== FUSED LEAF-OPERAND PANEL BUILDER =====================
Replaces, at a node of size 2*CUTL (whose 7 children are leaves), the pair
[ matop(node operand) ] + [ packA x7 / packB x7 inside the leaves ]
with ONE pass per side that reads the node operand's four quadrants once and
writes all 7 leaf-operand panels in the mikro's own layout. */
#define QN 64                    /* quadrant size == CUTL */
#define PN 7
typedef char pnl_cutl_is_64[(CUTL == 64) ? 1 : -1];
typedef char pnl_mr_is_6[(MR == 6) ? 1 : -1];
typedef char pnl_kc_is_kcp_off[(KCP == KC + 12) ? 1 : -1];

static long long PAN_A_RAW[PN][(size_t)MR*KCP + 32 + 8*1] __attribute__((aligned(64)));
static long long PAN_B_RAW[PN][(size_t)16*KC*4 + 32 + 8*0] __attribute__((aligned(64)));
static long long *const PAN_A[PN] = {
PAN_A_RAW[0]+8*1,
PAN_A_RAW[1]+8*1,
PAN_A_RAW[2]+8*1,
PAN_A_RAW[3]+8*1,
PAN_A_RAW[4]+8*1,
PAN_A_RAW[5]+8*1,
PAN_A_RAW[6]+8*1 };
static long long *const PAN_B[PN] = {
PAN_B_RAW[0]+8*0,
PAN_B_RAW[1]+8*0,
PAN_B_RAW[2]+8*0,
PAN_B_RAW[3]+8*0,
PAN_B_RAW[4]+8*0,
PAN_B_RAW[5]+8*0,
PAN_B_RAW[6]+8*0 };


/* A side: rows i in [0,64), k in [0,64).  Each row is 64 elements = 8 whole
64B lines, so the read pattern is the same as matop2's plane read. */
/* ---- HOISTED BUILDERS ------------------------------------------------------------------
`LDQ`/`STQ` carry a "memory" clobber, so gcc may NOT keep any memory-derived value
in a register across them.  The shipped builders therefore re-load `o->r[t]`,
`o->c[t]`, `o->s[t]` and re-multiply `o->r[t]*ld` INSIDE the k-loop (verified in the
g++-9 .s: `movslq 20(%rax),%rsi / imulq %r14,%rsi / movl 84(%rax),%edi / js` appear
in the innermost body).  Hoisting them into locals is the same loads, the same
stores and the same addresses -- only loop-invariant arithmetic moves out.
nt==2 (this row's dominant arity) is specialised; a hoisted generic path covers
every other arity and sign pattern.  Measured on the judge: buildB -6.09 %,
buildA null; both byte-exact at all four sign combinations. */
static void buildA(const Opd *o) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
if (nt == 2) {
const size_t qld = (size_t)QN * ld;
const long long *B0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *B1 = P + (size_t)o->r[1] * ld + o->c[1];
const int n0 = (o->s[0] < 0), n1 = (o->s[1] < 0);
const long long *pA = n0 ? B1 : B0, *pB = n0 ? B0 : B1;
const int sub = (n0 != n1), negall = (n0 && n1);
for (int i = 0; i < QN; i++) {
const size_t pi = (size_t)(i % MR);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
const long long *f0 = pA + (size_t)i * ld, *f1 = pB + (size_t)i * ld;
__m256i a[2][4], w, z = _mm256_setzero_si256();
for (int k = 0; k < QN; k += 8) {
/* BOTH 4-blocks' loads issue before either block's stores: the shipped loop
cannot do this (the "memory" clobber forbids gcc the reorder). */
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
a[b][0] = LDQ((const __m256i *)(f0 + kk));
w = LDQ((const __m256i *)(f1 + kk));
a[b][0] = sub ? _mm256_sub_epi64(a[b][0], w) : _mm256_add_epi64(a[b][0], w);
a[b][1] = LDQ((const __m256i *)(f0 + QN + kk));
w = LDQ((const __m256i *)(f1 + QN + kk));
a[b][1] = sub ? _mm256_sub_epi64(a[b][1], w) : _mm256_add_epi64(a[b][1], w);
a[b][2] = LDQ((const __m256i *)(f0 + qld + kk));
w = LDQ((const __m256i *)(f1 + qld + kk));
a[b][2] = sub ? _mm256_sub_epi64(a[b][2], w) : _mm256_add_epi64(a[b][2], w);
a[b][3] = LDQ((const __m256i *)(f0 + qld + QN + kk));
w = LDQ((const __m256i *)(f1 + qld + QN + kk));
a[b][3] = sub ? _mm256_sub_epi64(a[b][3], w) : _mm256_add_epi64(a[b][3], w);
if (negall) { for (int q = 0; q < 4; q++) a[b][q] = _mm256_sub_epi64(z, a[b][q]); }
}
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
STQ((__m256i *)(p0 + kk), _mm256_add_epi64(a[b][0], a[b][3]));  /* M1 A11+A22 */
STQ((__m256i *)(p1 + kk), _mm256_add_epi64(a[b][2], a[b][3]));  /* M2 A21+A22 */
STQ((__m256i *)(p2 + kk), a[b][0]);                             /* M3 A11     */
STQ((__m256i *)(p3 + kk), a[b][3]);                             /* M4 A22     */
STQ((__m256i *)(p4 + kk), _mm256_add_epi64(a[b][0], a[b][1]));  /* M5 A11+A12 */
STQ((__m256i *)(p5 + kk), _mm256_sub_epi64(a[b][2], a[b][0]));  /* M6 A21-A11 */
STQ((__m256i *)(p6 + kk), _mm256_sub_epi64(a[b][1], a[b][3]));  /* M7 A12-A22 */
}
}
}
return;
}
{ /* generic hoisted path: identical arithmetic, pointers and signs hoisted */
const long long *Bp[8]; int sg[8];
for (int t = 0; t < nt; t++) { Bp[t] = P + (size_t)o->r[t] * ld + o->c[t]; sg[t] = (o->s[t] < 0); }
const size_t qld = (size_t)QN * ld;
for (int i = 0; i < QN; i++) {
const size_t pi = (size_t)(i % MR);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
for (int k = 0; k < QN; k += 4) {
__m256i a11, a12, a21, a22;
if (nt == 1 && !sg[0]) {
a11 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + k));
a12 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + QN + k));
a21 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + k));
a22 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + QN + k));
} else {
a11 = _mm256_setzero_si256(); a12 = a11; a21 = a11; a22 = a11;
for (int t = 0; t < nt; t++) {
const long long *b = Bp[t] + (size_t)i * ld;
const int neg = sg[t]; __m256i w;
w = LDQ((const __m256i *)(b + k));
a11 = neg ? _mm256_sub_epi64(a11, w) : _mm256_add_epi64(a11, w);
w = LDQ((const __m256i *)(b + QN + k));
a12 = neg ? _mm256_sub_epi64(a12, w) : _mm256_add_epi64(a12, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + k));
a21 = neg ? _mm256_sub_epi64(a21, w) : _mm256_add_epi64(a21, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + QN + k));
a22 = neg ? _mm256_sub_epi64(a22, w) : _mm256_add_epi64(a22, w);
}
}
STQ((__m256i *)(p0 + k), _mm256_add_epi64(a11, a22));   /* M1  A11+A22 */
STQ((__m256i *)(p1 + k), _mm256_add_epi64(a21, a22));   /* M2  A21+A22 */
STQ((__m256i *)(p2 + k), a11);                          /* M3  A11     */
STQ((__m256i *)(p3 + k), a22);                          /* M4  A22     */
STQ((__m256i *)(p4 + k), _mm256_add_epi64(a11, a12));   /* M5  A11+A12 */
STQ((__m256i *)(p5 + k), _mm256_sub_epi64(a21, a11));   /* M6  A21-A11 */
STQ((__m256i *)(p6 + k), _mm256_sub_epi64(a12, a22));   /* M7  A12-A22 */
}
}
}

}


// One sign dispatch per B panel; every inner result remains modulo 2^64.
template<bool Subtract, bool Negate, bool Single=false>
static void buildB_signed_pair(const long long *s0,const long long *s1,int ld) {
const size_t qld=(size_t)QN*ld;
for(int jb=0;jb<QN;jb+=8) {
const size_t g=(size_t)(jb>>2)*KC*4;
long long *p0=PAN_B[0]+g,*p1=PAN_B[1]+g,*p2=PAN_B[2]+g;
long long *p3=PAN_B[3]+g,*p4=PAN_B[4]+g,*p5=PAN_B[5]+g,*p6=PAN_B[6]+g;
for(int k=0;k<QN;k++) {
const long long *x=s0+(size_t)k*ld+jb,*y=s1+(size_t)k*ld+jb;
auto pair=[&](size_t off) {
__m256i a=_mm256_loadu_si256((const __m256i*)(x+off));
__m256i v;
if constexpr(Single)v=a;
else {
 __m256i b=_mm256_loadu_si256((const __m256i*)(y+off));
 if constexpr(Subtract)v=_mm256_sub_epi64(a,b);
 else v=_mm256_add_epi64(a,b);
}
if constexpr(Negate)v=_mm256_sub_epi64(_mm256_setzero_si256(),v);
return v;
};
__m256i b11a=pair(0),b11b=pair(4),b12a=pair(QN),b12b=pair(QN+4);
__m256i b21a=pair(qld),b21b=pair(qld+4),b22a=pair(qld+QN),b22b=pair(qld+QN+4);
size_t ka=(size_t)k*4,kb=ka+(size_t)KC*4;
auto put=[](long long *p,__m256i v){_mm256_storeu_si256((__m256i*)p,_mm256_shuffle_epi32(v,0xB1));};
put(p0+ka,_mm256_add_epi64(b11a,b22a));put(p0+kb,_mm256_add_epi64(b11b,b22b));
put(p1+ka,b11a);put(p1+kb,b11b);
put(p2+ka,_mm256_sub_epi64(b12a,b22a));put(p2+kb,_mm256_sub_epi64(b12b,b22b));
put(p3+ka,_mm256_sub_epi64(b21a,b11a));put(p3+kb,_mm256_sub_epi64(b21b,b11b));
put(p4+ka,b22a);put(p4+kb,b22b);
put(p5+ka,_mm256_add_epi64(b11a,b12a));put(p5+kb,_mm256_add_epi64(b11b,b12b));
put(p6+ka,_mm256_add_epi64(b21a,b22a));put(p6+kb,_mm256_add_epi64(b21b,b22b));
}
}
}

static void buildB(const Opd *o) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
if(nt==1){
 const long long *src=P+(size_t)o->r[0]*ld+o->c[0];
 if(o->s[0]<0)buildB_signed_pair<false,true,true>(src,src,ld);
 else buildB_signed_pair<false,false,true>(src,src,ld);
 return;
}
if(nt==2) {
const long long *s0=P+(size_t)o->r[0]*ld+o->c[0];
const long long *s1=P+(size_t)o->r[1]*ld+o->c[1];
if(o->s[0]<0 && o->s[1]<0)buildB_signed_pair<false,true>(s0,s1,ld);
else if(o->s[0]<0)buildB_signed_pair<true,false>(s1,s0,ld);
else if(o->s[1]<0)buildB_signed_pair<true,false>(s0,s1,ld);
else buildB_signed_pair<false,false>(s0,s1,ld);
return;
}
const size_t qld = (size_t)QN * ld;
const long long *Bp[8]; int sg[8]; int allpos = 1;
for (int t = 0; t < nt; t++) {
Bp[t] = P + (size_t)o->r[t] * ld + o->c[t]; sg[t] = (o->s[t] < 0);
if (sg[t]) allpos = 0;
}
for (int jb = 0; jb < QN; jb += 8) {
const size_t g = (size_t)(jb >> 2) * KC * 4;
long long *p0 = PAN_B[0] + g, *p1 = PAN_B[1] + g, *p2 = PAN_B[2] + g,
*p3 = PAN_B[3] + g, *p4 = PAN_B[4] + g, *p5 = PAN_B[5] + g,
*p6 = PAN_B[6] + g;
for (int k = 0; k < QN; k++) {
__m256i b11a, b11b, b12a, b12b, b21a, b21b, b22a, b22b;
if (nt == 1 && allpos) {
const long long *r0p = Bp[0] + (size_t)k * ld;
const long long *r1p = r0p + qld;
b11a = LDQ((const __m256i *)(r0p + jb));        b11b = LDQ((const __m256i *)(r0p + jb + 4));
b12a = LDQ((const __m256i *)(r0p + QN + jb));   b12b = LDQ((const __m256i *)(r0p + QN + jb + 4));
b21a = LDQ((const __m256i *)(r1p + jb));        b21b = LDQ((const __m256i *)(r1p + jb + 4));
b22a = LDQ((const __m256i *)(r1p + QN + jb));   b22b = LDQ((const __m256i *)(r1p + QN + jb + 4));
} else if (nt == 2 && allpos) {
/* pairwise combine straight off the two base pointers: the shipped generic
loop re-derives Bp[0]/Bp[1] and re-tests the signs every (k,t). */
const long long *r0p = Bp[0] + (size_t)k * ld, *q0p = Bp[1] + (size_t)k * ld;
const long long *r1p = r0p + qld,            *q1p = q0p + qld;
__m256i w, x;
b11a = LDQ((const __m256i *)(r0p + jb));      x = LDQ((const __m256i *)(q0p + jb));
b11b = LDQ((const __m256i *)(r0p + jb + 4));  w = LDQ((const __m256i *)(q0p + jb + 4));
b11a = _mm256_add_epi64(b11a, x);             b11b = _mm256_add_epi64(b11b, w);
b12a = LDQ((const __m256i *)(r0p + QN + jb)); x = LDQ((const __m256i *)(q0p + QN + jb));
b12b = LDQ((const __m256i *)(r0p + QN + jb + 4)); w = LDQ((const __m256i *)(q0p + QN + jb + 4));
b12a = _mm256_add_epi64(b12a, x);             b12b = _mm256_add_epi64(b12b, w);
b21a = LDQ((const __m256i *)(r1p + jb));      x = LDQ((const __m256i *)(q1p + jb));
b21b = LDQ((const __m256i *)(r1p + jb + 4));  w = LDQ((const __m256i *)(q1p + jb + 4));
b21a = _mm256_add_epi64(b21a, x);             b21b = _mm256_add_epi64(b21b, w);
b22a = LDQ((const __m256i *)(r1p + QN + jb)); x = LDQ((const __m256i *)(q1p + QN + jb));
b22b = LDQ((const __m256i *)(r1p + QN + jb + 4)); w = LDQ((const __m256i *)(q1p + QN + jb + 4));
b22a = _mm256_add_epi64(b22a, x);             b22b = _mm256_add_epi64(b22b, w);
} else {
b11a = b11b = b12a = b12b = b21a = b21b = b22a = b22b = _mm256_setzero_si256();
for (int t = 0; t < nt; t++) {
const long long *r0p = Bp[t] + (size_t)k * ld;
const long long *r1p = r0p + qld;
const int neg = sg[t];
__m256i w;
w = LDQ((const __m256i *)(r0p + jb));         b11a = neg ? _mm256_sub_epi64(b11a, w) : _mm256_add_epi64(b11a, w);
w = LDQ((const __m256i *)(r0p + jb + 4));     b11b = neg ? _mm256_sub_epi64(b11b, w) : _mm256_add_epi64(b11b, w);
w = LDQ((const __m256i *)(r0p + QN + jb));    b12a = neg ? _mm256_sub_epi64(b12a, w) : _mm256_add_epi64(b12a, w);
w = LDQ((const __m256i *)(r0p + QN + jb + 4));b12b = neg ? _mm256_sub_epi64(b12b, w) : _mm256_add_epi64(b12b, w);
w = LDQ((const __m256i *)(r1p + jb));         b21a = neg ? _mm256_sub_epi64(b21a, w) : _mm256_add_epi64(b21a, w);
w = LDQ((const __m256i *)(r1p + jb + 4));     b21b = neg ? _mm256_sub_epi64(b21b, w) : _mm256_add_epi64(b21b, w);
w = LDQ((const __m256i *)(r1p + QN + jb));    b22a = neg ? _mm256_sub_epi64(b22a, w) : _mm256_add_epi64(b22a, w);
w = LDQ((const __m256i *)(r1p + QN + jb + 4));b22b = neg ? _mm256_sub_epi64(b22b, w) : _mm256_add_epi64(b22b, w);
}
}
const size_t ka = (size_t)k * 4, kb = ka + (size_t)KC * 4;
STQ((__m256i *)(p0 + ka), SWQ(_mm256_add_epi64(b11a, b22a)));
STQ((__m256i *)(p0 + kb), SWQ(_mm256_add_epi64(b11b, b22b)));
STQ((__m256i *)(p1 + ka), SWQ(b11a));
STQ((__m256i *)(p1 + kb), SWQ(b11b));
STQ((__m256i *)(p2 + ka), SWQ(_mm256_sub_epi64(b12a, b22a)));
STQ((__m256i *)(p2 + kb), SWQ(_mm256_sub_epi64(b12b, b22b)));
STQ((__m256i *)(p3 + ka), SWQ(_mm256_sub_epi64(b21a, b11a)));
STQ((__m256i *)(p3 + kb), SWQ(_mm256_sub_epi64(b21b, b11b)));
STQ((__m256i *)(p4 + ka), SWQ(b22a));
STQ((__m256i *)(p4 + kb), SWQ(b22b));
STQ((__m256i *)(p5 + ka), SWQ(_mm256_add_epi64(b11a, b12a)));
STQ((__m256i *)(p5 + kb), SWQ(_mm256_add_epi64(b11b, b12b)));
STQ((__m256i *)(p6 + ka), SWQ(_mm256_add_epi64(b21a, b22a)));
STQ((__m256i *)(p6 + kb), SWQ(_mm256_add_epi64(b21b, b22b)));
}
}
}

/* the leaf, driven from prebuilt panels.  Loop nest, order, mikro call
sequence and destination addressing are the shipped leaf's; packA / packB
are replaced by the panel pointers. */
static void leafP(int n, int L, long long *C1, int ldc1, int op1,
long long *C2, int ldc2, int op2) {
const long long *pa = PAN_A[L];
const long long *pb = PAN_B[L];
for (int ic = 0; ic < n; ic += MC) {
int mc = (n - ic < MC) ? (n - ic) : MC;
for (int ir = 0; ir < mc; ir += MR) {
int rows = (mc - ir < MR) ? (mc - ir) : MR;
const long long *ap = pa + (size_t)(ic + ir) * KCP;
long long *c1p = C1 + (size_t)(ic + ir) * ldc1;
long long *c2p = C2 ? C2 + (size_t)(ic + ir) * ldc2 : 0;
for (int jg = 0; jg < n; jg += 4) {
mikro(n, ap, pb + (size_t)(jg >> 2) * KC * 4, c1p + jg, ldc1, op1,
c2p ? c2p + jg : 0, ldc2, op2, rows);
}
}
}
}

static void leafP_blk(int n, int L, int ic0, int ic1, long long *C1, int ldc1, int op1,
long long *C2, int ldc2, int op2) {
const long long *pa = PAN_A[L];
const long long *pb = PAN_B[L];
for (int ic = ic0; ic < ic1; ic += MC) {
int mc = (n - ic < MC) ? (n - ic) : MC;
if (mc > ic1 - ic) mc = ic1 - ic;
for (int ir = 0; ir < mc; ir += MR) {
int rows = (mc - ir < MR) ? (mc - ir) : MR;
const long long *ap = pa + (size_t)(ic + ir - ic0) * KCP;
long long *c1p = C1 + (size_t)(ic + ir) * ldc1;
long long *c2p = C2 ? C2 + (size_t)(ic + ir) * ldc2 : 0;
if(n==64 && !C2 && (rows==4||rows==6)){
if(rows==6){
if(op1==0) columns64_static<6,0>(ap,pb,c1p,ldc1);
else if(op1==1) columns64_static<6,1>(ap,pb,c1p,ldc1);
else columns64_static<6,2>(ap,pb,c1p,ldc1);
}else{
if(op1==0) columns64_static<4,0>(ap,pb,c1p,ldc1);
else if(op1==1) columns64_static<4,1>(ap,pb,c1p,ldc1);
else columns64_static<4,2>(ap,pb,c1p,ldc1);
}
}else {
for (int jg = 0; jg < n; jg += 4) {
mikro(n, ap, pb + (size_t)(jg >> 2) * KC * 4, c1p + jg, ldc1, op1,
c2p ? c2p + jg : 0, ldc2, op2, rows);
}
}
}
}
}


static long long g_scratch[6u * 1024u * 1024u] __attribute__((aligned(64)));

/* dst op= sign * T   (op: 0 store, 1 add, 2 sub) */
static void cmb(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op) {
__m256i sgn = _mm256_set1_epi64x(sign < 0 ? -1 : 0);
for (int i = 0; i < h; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c = dst + (size_t)i * ldd;
int j = 0;
/* ARM cmbu16: the j-loop's inner step raised 4 -> 16 elements so FOUR loads
are issued before the first store.  `LDQ` carries a "memory" clobber, so
the shipped 4-element step cannot be software-pipelined by gcc; here the
pipelining is explicit.  Same loads, same stores, same addresses. */
if (op == 0) {
for (; j + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)),     v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
if (sign < 0) { __m256i z = _mm256_setzero_si256();
v0 = _mm256_sub_epi64(z, v0); v1 = _mm256_sub_epi64(z, v1);
v2 = _mm256_sub_epi64(z, v2); v3 = _mm256_sub_epi64(z, v3); }
STQ((__m256i *)(c + j), v0);      STQ((__m256i *)(c + j + 4), v1);
STQ((__m256i *)(c + j + 8), v2);  STQ((__m256i *)(c + j + 12), v3);
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
if (sign < 0) v = _mm256_sub_epi64(_mm256_setzero_si256(), v);
STQ((__m256i *)(c + j), v);
}
for (; j < h; j++) c[j] = sign < 0 ? -t[j] : t[j];
} else {
for (; j + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)),     v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
__m256i w0 = LDQ((const __m256i *)(c + j)),     w1 = LDQ((const __m256i *)(c + j + 4));
__m256i w2 = LDQ((const __m256i *)(c + j + 8)), w3 = LDQ((const __m256i *)(c + j + 12));
if ((op == 1) != (sign < 0)) {
STQ((__m256i *)(c + j),      _mm256_add_epi64(w0, v0));
STQ((__m256i *)(c + j + 4),  _mm256_add_epi64(w1, v1));
STQ((__m256i *)(c + j + 8),  _mm256_add_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_add_epi64(w3, v3));
} else {
STQ((__m256i *)(c + j),      _mm256_sub_epi64(w0, v0));
STQ((__m256i *)(c + j + 4),  _mm256_sub_epi64(w1, v1));
STQ((__m256i *)(c + j + 8),  _mm256_sub_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_sub_epi64(w3, v3));
}
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
__m256i w = LDQ((const __m256i *)(c + j));
if ((op == 1) != (sign < 0)) STQ((__m256i *)(c + j), _mm256_add_epi64(w, v));
else                         STQ((__m256i *)(c + j), _mm256_sub_epi64(w, v));
}
for (; j < h; j++) { if ((op == 1) != (sign < 0)) c[j] += t[j]; else c[j] -= t[j]; }
}
}
}

/* fused two-destination combine: reads T ONCE for both destinations.
Semantics identical to cmb(d1,s1,o1) then cmb(d2,s2,o2) on the same T. */
static void cmb2(long long *d1,int ld1,int op1,int s1,
long long *d2,int ld2,int op2,int s2,
const long long *T,int ldt,int h) {
for (int i = 0; i < h; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c1 = d1 + (size_t)i * ld1, *c2 = d2 + (size_t)i * ld2;
int j = 0;
for (; j + 16 <= h; j += 16) {
__m256i v0=LDQ((const __m256i*)(t+j)),   v1=LDQ((const __m256i*)(t+j+4));
__m256i v2=LDQ((const __m256i*)(t+j+8)), v3=LDQ((const __m256i*)(t+j+12));
__m256i a0=LDQ((const __m256i*)(c1+j)),   a1=LDQ((const __m256i*)(c1+j+4));
__m256i a2=LDQ((const __m256i*)(c1+j+8)), a3=LDQ((const __m256i*)(c1+j+12));
__m256i b0=LDQ((const __m256i*)(c2+j)),   b1=LDQ((const __m256i*)(c2+j+4));
__m256i b2=LDQ((const __m256i*)(c2+j+8)), b3=LDQ((const __m256i*)(c2+j+12));
if ((op1 == 1) != (s1 < 0)) { STQ((__m256i*)(c1+j),_mm256_add_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_add_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_add_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_add_epi64(a3,v3)); }
else { STQ((__m256i*)(c1+j),_mm256_sub_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_sub_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_sub_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_sub_epi64(a3,v3)); }
if ((op2 == 1) != (s2 < 0)) { STQ((__m256i*)(c2+j),_mm256_add_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_add_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_add_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_add_epi64(b3,v3)); }
else { STQ((__m256i*)(c2+j),_mm256_sub_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_sub_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_sub_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_sub_epi64(b3,v3)); }
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
__m256i a = LDQ((const __m256i *)(c1 + j));
STQ((__m256i *)(c1 + j), ((op1 == 1) != (s1 < 0)) ? _mm256_add_epi64(a,v)
: _mm256_sub_epi64(a,v));
__m256i b = LDQ((const __m256i *)(c2 + j));
STQ((__m256i *)(c2 + j), ((op2 == 1) != (s2 < 0)) ? _mm256_add_epi64(b,v)
: _mm256_sub_epi64(b,v));
}
for (; j < h; j++) {
long long v = t[j];
if ((op1 == 1) != (s1 < 0)) c1[j] += v; else c1[j] -= v;
if ((op2 == 1) != (s2 < 0)) c2[j] += v; else c2[j] -= v;
}
}
}


/* l1k15_rv `fus3`: the three products whose SECOND destination is C22 (M1, M2, M3)
land in three quadrants that are ALREADY LIVE (C11, C21, C12) -- so the three
sequential combine passes over C22 collapse into ONE pass.  Source bytes are
identical; the destination is written ONCE instead of three times.
Semantics == cmb(C22,C11,+1); cmb(C22,C21,-1); cmb(C22,C12,+1)   [op==0] */
static void cmb3_p(long long *d, int ldd,
const long long *a, const long long *b, const long long *c, int ld,
int rows, int ncol) {
for (int i = 0; i < rows; i++) {
const long long *pa = a + (size_t)i * ld, *pb = b + (size_t)i * ld, *pc = c + (size_t)i * ld;
long long *pd = d + (size_t)i * ldd;
int j = 0;
for (; j + 16 <= ncol; j += 16) {
__m256i A0=LDQ((const __m256i*)(pa+j)),   A1=LDQ((const __m256i*)(pa+j+4));
__m256i A2=LDQ((const __m256i*)(pa+j+8)), A3=LDQ((const __m256i*)(pa+j+12));
__m256i B0=LDQ((const __m256i*)(pb+j)),   B1=LDQ((const __m256i*)(pb+j+4));
__m256i B2=LDQ((const __m256i*)(pb+j+8)), B3=LDQ((const __m256i*)(pb+j+12));
__m256i C0=LDQ((const __m256i*)(pc+j)),   C1=LDQ((const __m256i*)(pc+j+4));
__m256i C2=LDQ((const __m256i*)(pc+j+8)), C3=LDQ((const __m256i*)(pc+j+12));
STQ((__m256i*)(pd+j),    _mm256_sub_epi64(_mm256_add_epi64(A0,C0),B0));
STQ((__m256i*)(pd+j+4),  _mm256_sub_epi64(_mm256_add_epi64(A1,C1),B1));
STQ((__m256i*)(pd+j+8),  _mm256_sub_epi64(_mm256_add_epi64(A2,C2),B2));
STQ((__m256i*)(pd+j+12), _mm256_sub_epi64(_mm256_add_epi64(A3,C3),B3));
}
for (; j + 4 <= ncol; j += 4)
STQ((__m256i*)(pd+j), _mm256_sub_epi64(
_mm256_add_epi64(LDQ((const __m256i*)(pa+j)), LDQ((const __m256i*)(pc+j))),
LDQ((const __m256i*)(pb+j))));
for (; j < ncol; j++) pd[j] = pa[j] - pb[j] + pc[j];
}
}

static void cmb3_ntx(long long *d, int ldd,
const long long *a, const long long *b, const long long *c, int ld,
int rows, int ncol) {
for (int i = 0; i < rows; i++) {
const long long *pa = a + (size_t)i * ld, *pb = b + (size_t)i * ld, *pc = c + (size_t)i * ld;
long long *pd = d + (size_t)i * ldd;
int j = 0;
for (; j + 16 <= ncol; j += 16) {
__m256i A0=LDQ((const __m256i*)(pa+j)),   A1=LDQ((const __m256i*)(pa+j+4));
__m256i A2=LDQ((const __m256i*)(pa+j+8)), A3=LDQ((const __m256i*)(pa+j+12));
__m256i B0=LDQ((const __m256i*)(pb+j)),   B1=LDQ((const __m256i*)(pb+j+4));
__m256i B2=LDQ((const __m256i*)(pb+j+8)), B3=LDQ((const __m256i*)(pb+j+12));
__m256i C0=LDQ((const __m256i*)(pc+j)),   C1=LDQ((const __m256i*)(pc+j+4));
__m256i C2=LDQ((const __m256i*)(pc+j+8)), C3=LDQ((const __m256i*)(pc+j+12));
STQNT((__m256i*)(pd+j),    _mm256_sub_epi64(_mm256_add_epi64(A0,C0),B0));
STQNT((__m256i*)(pd+j+4),  _mm256_sub_epi64(_mm256_add_epi64(A1,C1),B1));
STQNT((__m256i*)(pd+j+8),  _mm256_sub_epi64(_mm256_add_epi64(A2,C2),B2));
STQNT((__m256i*)(pd+j+12), _mm256_sub_epi64(_mm256_add_epi64(A3,C3),B3));
}
for (; j + 4 <= ncol; j += 4)
STQNT((__m256i*)(pd+j), _mm256_sub_epi64(
_mm256_add_epi64(LDQ((const __m256i*)(pa+j)), LDQ((const __m256i*)(pc+j))),
LDQ((const __m256i*)(pb+j))));
for (; j < ncol; j++) pd[j] = pa[j] - pb[j] + pc[j];
}
}

static void cmb3_blk(long long *d, int ldd, const long long *a, const long long *b,
const long long *c, int ld, int h, int i0, int i1) {
if (i1 <= i0) return;
cmb3(d + (size_t)i0 * ldd, ldd, a + (size_t)i0 * ld, b + (size_t)i0 * ld,
c + (size_t)i0 * ld, ld, i1 - i0, h);
}

/* single-destination recursive Strassen: d op= Ao*Bo, with T a scratch block of h*h */
/* materialise one level's operands ONCE per parent call, so the LEAF's packing reads 2
contiguous sources per element instead of 4 strided ones.  Only applied at the level whose
children are leaves (n == 2*CUTL): the cost is 2*n^2 element-writes per parent, and it is
amortised over that parent's 7 leaf calls, which read 4x more elements than it writes. */
static long long matbufA[3][262144 + 64], matbufB[3][262144 + 64];
/* lane m142: THE FUSION'S MAXIMUM FORM, built as a measurement.
`matop`'s t-loop over nt is a RUNTIME loop with a sign branch per term; the dominant
arity on this row is nt==2 (m101's firing table: "by arity: all nt=2").  A FULLY
UNROLLED nt==2 path removes exactly the work a "fuse the materialisation into the
combine" construction could remove -- the per-element CONTROL -- while leaving every
byte of traffic identical.  If THIS does not pay, no pass-fusion can: the fusion
changes no traffic either (2 dense reads + 1 dense write, either way).  Correct by
construction == the generic path's arithmetic. */
static void matop2_p(const Opd *o, int n, long long *dst) {
const long long *P = o->p; const int ld = o->ld;
const long long *p0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *p1 = P + (size_t)o->r[1] * ld + o->c[1];
const int s1neg = (o->s[1] < 0);
int i = 0;
/* TEST ARM: the i-loop unrolled by 2.  Two destination rows are in flight,
so the row-jump latency that serialises the shipped loop is exposed once
per TWO rows instead of once per row.  Every load, store and address is
the same as the shipped function's; only the ORDER changes.
The unrolled path is taken only when n % 16 == 0, because the shipped
tail loop writes up to 3 elements PAST the row end and rows are written
in order there; the real call sizes are 128/256/512 and are all % 16 == 0. */
if ((n & 15) == 0) for (; i + 2 <= n; i += 2) {
long long *dpa = dst + (size_t)i * n;
long long *dpb = dpa + n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
const long long *c = a + ld, *e = b + ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 4) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 4) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)),   x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)),   y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)),   z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQ((__m256i *)(dpa+j),    _mm256_add_epi64(w0,y0));
STQ((__m256i *)(dpa+j+4),  _mm256_add_epi64(w1,y1));
STQ((__m256i *)(dpa+j+8),  _mm256_add_epi64(w2,y2));
STQ((__m256i *)(dpa+j+12), _mm256_add_epi64(w3,y3));
STQ((__m256i *)(dpb+j),    _mm256_add_epi64(x0,z0));
STQ((__m256i *)(dpb+j+4),  _mm256_add_epi64(x1,z1));
STQ((__m256i *)(dpb+j+8),  _mm256_add_epi64(x2,z2));
STQ((__m256i *)(dpb+j+12), _mm256_add_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQ((__m256i *)(dpa+j), _mm256_add_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQ((__m256i *)(dpb+j), _mm256_add_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)),   x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)),   y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)),   z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQ((__m256i *)(dpa+j),    _mm256_sub_epi64(w0,y0));
STQ((__m256i *)(dpa+j+4),  _mm256_sub_epi64(w1,y1));
STQ((__m256i *)(dpa+j+8),  _mm256_sub_epi64(w2,y2));
STQ((__m256i *)(dpa+j+12), _mm256_sub_epi64(w3,y3));
STQ((__m256i *)(dpb+j),    _mm256_sub_epi64(x0,z0));
STQ((__m256i *)(dpb+j+4),  _mm256_sub_epi64(x1,z1));
STQ((__m256i *)(dpb+j+8),  _mm256_sub_epi64(x2,z2));
STQ((__m256i *)(dpb+j+12), _mm256_sub_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQ((__m256i *)(dpa+j), _mm256_sub_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQ((__m256i *)(dpb+j), _mm256_sub_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
}
}
for (; i < n; i++) {
long long *dp = dst + (size_t)i * n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 2) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 2) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_add_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_add_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_add_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_add_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQ((__m256i *)(dp+j), w0); STQ((__m256i *)(dp+j+4), w1);
STQ((__m256i *)(dp+j+8), w2); STQ((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQ((__m256i *)(dp+j), _mm256_add_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_sub_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_sub_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_sub_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_sub_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQ((__m256i *)(dp+j), w0); STQ((__m256i *)(dp+j+4), w1);
STQ((__m256i *)(dp+j+8), w2); STQ((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQ((__m256i *)(dp+j), _mm256_sub_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
}
}
}

static void matop2_ntx(const Opd *o, int n, long long *dst) {
const long long *P = o->p; const int ld = o->ld;
const long long *p0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *p1 = P + (size_t)o->r[1] * ld + o->c[1];
const int s1neg = (o->s[1] < 0);
int i = 0;
/* TEST ARM: the i-loop unrolled by 2.  Two destination rows are in flight,
so the row-jump latency that serialises the shipped loop is exposed once
per TWO rows instead of once per row.  Every load, store and address is
the same as the shipped function's; only the ORDER changes.
The unrolled path is taken only when n % 16 == 0, because the shipped
tail loop writes up to 3 elements PAST the row end and rows are written
in order there; the real call sizes are 128/256/512 and are all % 16 == 0. */
if ((n & 15) == 0) for (; i + 2 <= n; i += 2) {
long long *dpa = dst + (size_t)i * n;
long long *dpb = dpa + n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
const long long *c = a + ld, *e = b + ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 4) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 4) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)),   x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)),   y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)),   z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQNT((__m256i *)(dpa+j),    _mm256_add_epi64(w0,y0));
STQNT((__m256i *)(dpa+j+4),  _mm256_add_epi64(w1,y1));
STQNT((__m256i *)(dpa+j+8),  _mm256_add_epi64(w2,y2));
STQNT((__m256i *)(dpa+j+12), _mm256_add_epi64(w3,y3));
STQNT((__m256i *)(dpb+j),    _mm256_add_epi64(x0,z0));
STQNT((__m256i *)(dpb+j+4),  _mm256_add_epi64(x1,z1));
STQNT((__m256i *)(dpb+j+8),  _mm256_add_epi64(x2,z2));
STQNT((__m256i *)(dpb+j+12), _mm256_add_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQNT((__m256i *)(dpa+j), _mm256_add_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQNT((__m256i *)(dpb+j), _mm256_add_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)),   x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)),   y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)),   z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQNT((__m256i *)(dpa+j),    _mm256_sub_epi64(w0,y0));
STQNT((__m256i *)(dpa+j+4),  _mm256_sub_epi64(w1,y1));
STQNT((__m256i *)(dpa+j+8),  _mm256_sub_epi64(w2,y2));
STQNT((__m256i *)(dpa+j+12), _mm256_sub_epi64(w3,y3));
STQNT((__m256i *)(dpb+j),    _mm256_sub_epi64(x0,z0));
STQNT((__m256i *)(dpb+j+4),  _mm256_sub_epi64(x1,z1));
STQNT((__m256i *)(dpb+j+8),  _mm256_sub_epi64(x2,z2));
STQNT((__m256i *)(dpb+j+12), _mm256_sub_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQNT((__m256i *)(dpa+j), _mm256_sub_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQNT((__m256i *)(dpb+j), _mm256_sub_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
}
}
for (; i < n; i++) {
long long *dp = dst + (size_t)i * n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 2) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 2) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_add_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_add_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_add_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_add_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQNT((__m256i *)(dp+j), w0); STQNT((__m256i *)(dp+j+4), w1);
STQNT((__m256i *)(dp+j+8), w2); STQNT((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQNT((__m256i *)(dp+j), _mm256_add_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)),   w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_sub_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_sub_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_sub_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_sub_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQNT((__m256i *)(dp+j), w0); STQNT((__m256i *)(dp+j+4), w1);
STQNT((__m256i *)(dp+j+8), w2); STQNT((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQNT((__m256i *)(dp+j), _mm256_sub_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
}
}
}


/* [ml1d_] NT cold-write knife.  matop's destination is a *fresh* n*n plane written
wholesale and consumed only after >L3 of intervening traffic (see head).  For n>=512
that plane is 2 MB and provably evicted before its consumer -> the write-allocate
(RFO) read is pure waste -> vmovntdq.  For n==256 the plane is 512 KB and stays
L3-resident -> NT would only destroy the consumer's hits (MEASURED +2.10 ms) -> keep
the cacheable store.  cmb3's destination is never read back at all -> NT for the
large planes; the 2695 sub-64 leaf tiles stay cacheable. */
static void matop2(const Opd *o, int n, long long *dst) {
if (n >= 512) { matop2_ntx(o, n, dst); return; }
matop2_p(o, n, dst);
}
static void cmb3(long long *d, int ldd, const long long *a, const long long *b,
const long long *c, int ld, int rows, int ncol) {
if (rows >= 256 && ((unsigned long)d & 31u) == 0 && (ldd & 3) == 0) {
cmb3_ntx(d, ldd, a, b, c, ld, rows, ncol); return;
}
cmb3_p(d, ldd, a, b, c, ld, rows, ncol);
}

static void matop(const Opd *o, int n, long long *dst) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
if (nt == 2) { matop2(o, n, dst); return; }
for (int i = 0; i < n; i++) {
long long *dp = dst + (size_t)i * n;
int j = 0;
for (; j + 16 <= n; j += 16) {
__m256i v0=_mm256_setzero_si256(),v1=v0,v2=v0,v3=v0;
for (int t = 0; t < nt; t++) {
const long long *sp = P + (size_t)(o->r[t] + i) * ld + o->c[t] + j;
__m256i w0=LDQ((const __m256i *)sp), w1=LDQ((const __m256i *)(sp+4));
__m256i w2=LDQ((const __m256i *)(sp+8)), w3=LDQ((const __m256i *)(sp+12));
if (o->s[t] < 0) { v0=_mm256_sub_epi64(v0,w0); v1=_mm256_sub_epi64(v1,w1);
v2=_mm256_sub_epi64(v2,w2); v3=_mm256_sub_epi64(v3,w3); }
else             { v0=_mm256_add_epi64(v0,w0); v1=_mm256_add_epi64(v1,w1);
v2=_mm256_add_epi64(v2,w2); v3=_mm256_add_epi64(v3,w3); }
}
STQ((__m256i *)(dp+j), v0); STQ((__m256i *)(dp+j+4), v1);
STQ((__m256i *)(dp+j+8), v2); STQ((__m256i *)(dp+j+12), v3);
}
for (; j < n; j += 4) {
__m256i v = _mm256_setzero_si256();
for (int t = 0; t < nt; t++) {
__m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[t] + i) * ld + o->c[t] + j));
v = (o->s[t] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
}
STQ((__m256i *)(dp + j), v);
}
}
}
/* the strassen node whose 7 children are leaves, driven from one fused pass */
static void leaflevel(int n, const Opd *Ao, const Opd *Bo, long long *d, int ld, int op,
long long *T) {
const int h = n >> 1;
long long *Q[4];
Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
int wr[4] = {0, 0, 0, 0};
buildA(Ao); buildB(Bo);
#define MKL(qi, si, qj, sj, L) do { \
int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
int ob = wr[qj] ? effop(1, (sj)) : effop(op, (sj)); \
wr[qi] = 1; wr[qj] = 1; \
if (oa == 0) { \
leafP(h, (L), Q[qi], ld, 0, 0, 0, 0); \
cmb(Q[qj], ld, Q[qi], ld, h, (sj), ob); \
} else { \
leafP(h, (L), T, h, 0, 0, 0, 0); \
cmb2(Q[qi], ld, oa, (si), Q[qj], ld, ob, (sj), T, h, h); \
} \
} while (0)
#define MKL1(qi, si, L) do { \
int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
wr[qi] = 1; \
leafP(h, (L), Q[qi], ld, oa, 0, 0, 0); \
} while (0)
{ MKL(0, 1, 3, 1, 0); }   /* M1 = (A11+A22)(B11+B22) -> +C11, +C22 */
{ MKL(2, 1, 3, 2, 1); }   /* M2 = (A21+A22)B11       -> +C21, -C22 */
{ MKL(1, 1, 3, 1, 2); }   /* M3 = A11(B12-B22)       -> +C12, +C22 */
{ MKL(0, 1, 2, 1, 3); }   /* M4 = A22(B21-B11)       -> +C11, +C21 */
{ MKL(0, 2, 1, 1, 4); }   /* M5 = (A11+A12)B22       -> -C11, +C12 */
{ MKL1(3, 1, 5); }        /* M6 = (A21-A11)(B11+B12) -> +C22      */
{ MKL1(0, 1, 6); }        /* M7 = (A12-A22)(B21+B22) -> +C11      */
#undef MKL
#undef MKL1
}

static void buildA_blk(const Opd *o, int i0, int i1) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
if (nt == 2) {
const size_t qld = (size_t)QN * ld;
const long long *B0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *B1 = P + (size_t)o->r[1] * ld + o->c[1];
const int n0 = (o->s[0] < 0), n1 = (o->s[1] < 0);
const long long *pA = n0 ? B1 : B0, *pB = n0 ? B0 : B1;
const int sub = (n0 != n1), negall = (n0 && n1);
for (int i = i0; i < i1; i++) {
const size_t pi = (size_t)(i - i0);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
const long long *f0 = pA + (size_t)i * ld, *f1 = pB + (size_t)i * ld;
__m256i a[2][4], w, z = _mm256_setzero_si256();
for (int k = 0; k < QN; k += 8) {
/* BOTH 4-blocks' loads issue before either block's stores: the shipped loop
cannot do this (the "memory" clobber forbids gcc the reorder). */
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
a[b][0] = LDQ((const __m256i *)(f0 + kk));
w = LDQ((const __m256i *)(f1 + kk));
a[b][0] = sub ? _mm256_sub_epi64(a[b][0], w) : _mm256_add_epi64(a[b][0], w);
a[b][1] = LDQ((const __m256i *)(f0 + QN + kk));
w = LDQ((const __m256i *)(f1 + QN + kk));
a[b][1] = sub ? _mm256_sub_epi64(a[b][1], w) : _mm256_add_epi64(a[b][1], w);
a[b][2] = LDQ((const __m256i *)(f0 + qld + kk));
w = LDQ((const __m256i *)(f1 + qld + kk));
a[b][2] = sub ? _mm256_sub_epi64(a[b][2], w) : _mm256_add_epi64(a[b][2], w);
a[b][3] = LDQ((const __m256i *)(f0 + qld + QN + kk));
w = LDQ((const __m256i *)(f1 + qld + QN + kk));
a[b][3] = sub ? _mm256_sub_epi64(a[b][3], w) : _mm256_add_epi64(a[b][3], w);
if (negall) { for (int q = 0; q < 4; q++) a[b][q] = _mm256_sub_epi64(z, a[b][q]); }
}
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
STQ((__m256i *)(p0 + kk), _mm256_add_epi64(a[b][0], a[b][3]));  /* M1 A11+A22 */
STQ((__m256i *)(p1 + kk), _mm256_add_epi64(a[b][2], a[b][3]));  /* M2 A21+A22 */
STQ((__m256i *)(p2 + kk), a[b][0]);                             /* M3 A11     */
STQ((__m256i *)(p3 + kk), a[b][3]);                             /* M4 A22     */
STQ((__m256i *)(p4 + kk), _mm256_add_epi64(a[b][0], a[b][1]));  /* M5 A11+A12 */
STQ((__m256i *)(p5 + kk), _mm256_sub_epi64(a[b][2], a[b][0]));  /* M6 A21-A11 */
STQ((__m256i *)(p6 + kk), _mm256_sub_epi64(a[b][1], a[b][3]));  /* M7 A12-A22 */
}
}
}
return;
}
{ /* generic hoisted path: identical arithmetic, pointers and signs hoisted */
const long long *Bp[8]; int sg[8];
for (int t = 0; t < nt; t++) { Bp[t] = P + (size_t)o->r[t] * ld + o->c[t]; sg[t] = (o->s[t] < 0); }
const size_t qld = (size_t)QN * ld;
for (int i = i0; i < i1; i++) {
const size_t pi = (size_t)(i - i0);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
for (int k = 0; k < QN; k += 4) {
__m256i a11, a12, a21, a22;
if (nt == 1 && !sg[0]) {
a11 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + k));
a12 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + QN + k));
a21 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + k));
a22 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + QN + k));
} else {
a11 = _mm256_setzero_si256(); a12 = a11; a21 = a11; a22 = a11;
for (int t = 0; t < nt; t++) {
const long long *b = Bp[t] + (size_t)i * ld;
const int neg = sg[t]; __m256i w;
w = LDQ((const __m256i *)(b + k));
a11 = neg ? _mm256_sub_epi64(a11, w) : _mm256_add_epi64(a11, w);
w = LDQ((const __m256i *)(b + QN + k));
a12 = neg ? _mm256_sub_epi64(a12, w) : _mm256_add_epi64(a12, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + k));
a21 = neg ? _mm256_sub_epi64(a21, w) : _mm256_add_epi64(a21, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + QN + k));
a22 = neg ? _mm256_sub_epi64(a22, w) : _mm256_add_epi64(a22, w);
}
}
STQ((__m256i *)(p0 + k), _mm256_add_epi64(a11, a22));   /* M1  A11+A22 */
STQ((__m256i *)(p1 + k), _mm256_add_epi64(a21, a22));   /* M2  A21+A22 */
STQ((__m256i *)(p2 + k), a11);                          /* M3  A11     */
STQ((__m256i *)(p3 + k), a22);                          /* M4  A22     */
STQ((__m256i *)(p4 + k), _mm256_add_epi64(a11, a12));   /* M5  A11+A12 */
STQ((__m256i *)(p5 + k), _mm256_sub_epi64(a21, a11));   /* M6  A21-A11 */
STQ((__m256i *)(p6 + k), _mm256_sub_epi64(a12, a22));   /* M7  A12-A22 */
}
}
}

}


static void cmb_blk(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op, int i0, int i1) {
__m256i sgn = _mm256_set1_epi64x(sign < 0 ? -1 : 0);
for (int i = i0; i < i1; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c = dst + (size_t)i * ldd;
int j = 0;
/* ARM cmbu16: the j-loop's inner step raised 4 -> 16 elements so FOUR loads
are issued before the first store.  `LDQ` carries a "memory" clobber, so
the shipped 4-element step cannot be software-pipelined by gcc; here the
pipelining is explicit.  Same loads, same stores, same addresses. */
if (op == 0) {
for (; j + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)),     v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
if (sign < 0) { __m256i z = _mm256_setzero_si256();
v0 = _mm256_sub_epi64(z, v0); v1 = _mm256_sub_epi64(z, v1);
v2 = _mm256_sub_epi64(z, v2); v3 = _mm256_sub_epi64(z, v3); }
STQ((__m256i *)(c + j), v0);      STQ((__m256i *)(c + j + 4), v1);
STQ((__m256i *)(c + j + 8), v2);  STQ((__m256i *)(c + j + 12), v3);
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
if (sign < 0) v = _mm256_sub_epi64(_mm256_setzero_si256(), v);
STQ((__m256i *)(c + j), v);
}
for (; j < h; j++) c[j] = sign < 0 ? -t[j] : t[j];
} else {
for (; j + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)),     v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
__m256i w0 = LDQ((const __m256i *)(c + j)),     w1 = LDQ((const __m256i *)(c + j + 4));
__m256i w2 = LDQ((const __m256i *)(c + j + 8)), w3 = LDQ((const __m256i *)(c + j + 12));
if ((op == 1) != (sign < 0)) {
STQ((__m256i *)(c + j),      _mm256_add_epi64(w0, v0));
STQ((__m256i *)(c + j + 4),  _mm256_add_epi64(w1, v1));
STQ((__m256i *)(c + j + 8),  _mm256_add_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_add_epi64(w3, v3));
} else {
STQ((__m256i *)(c + j),      _mm256_sub_epi64(w0, v0));
STQ((__m256i *)(c + j + 4),  _mm256_sub_epi64(w1, v1));
STQ((__m256i *)(c + j + 8),  _mm256_sub_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_sub_epi64(w3, v3));
}
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
__m256i w = LDQ((const __m256i *)(c + j));
if ((op == 1) != (sign < 0)) STQ((__m256i *)(c + j), _mm256_add_epi64(w, v));
else                         STQ((__m256i *)(c + j), _mm256_sub_epi64(w, v));
}
for (; j < h; j++) { if ((op == 1) != (sign < 0)) c[j] += t[j]; else c[j] -= t[j]; }
}
}
}


static void cmb2_blk(long long *d1,int ld1,int op1,int s1,
long long *d2,int ld2,int op2,int s2,
const long long *T,int ldt,int h, int i0, int i1) {
for (int i = i0; i < i1; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c1 = d1 + (size_t)i * ld1, *c2 = d2 + (size_t)i * ld2;
int j = 0;
for (; j + 16 <= h; j += 16) {
__m256i v0=LDQ((const __m256i*)(t+j)),   v1=LDQ((const __m256i*)(t+j+4));
__m256i v2=LDQ((const __m256i*)(t+j+8)), v3=LDQ((const __m256i*)(t+j+12));
__m256i a0=LDQ((const __m256i*)(c1+j)),   a1=LDQ((const __m256i*)(c1+j+4));
__m256i a2=LDQ((const __m256i*)(c1+j+8)), a3=LDQ((const __m256i*)(c1+j+12));
__m256i b0=LDQ((const __m256i*)(c2+j)),   b1=LDQ((const __m256i*)(c2+j+4));
__m256i b2=LDQ((const __m256i*)(c2+j+8)), b3=LDQ((const __m256i*)(c2+j+12));
if ((op1 == 1) != (s1 < 0)) { STQ((__m256i*)(c1+j),_mm256_add_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_add_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_add_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_add_epi64(a3,v3)); }
else { STQ((__m256i*)(c1+j),_mm256_sub_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_sub_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_sub_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_sub_epi64(a3,v3)); }
if ((op2 == 1) != (s2 < 0)) { STQ((__m256i*)(c2+j),_mm256_add_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_add_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_add_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_add_epi64(b3,v3)); }
else { STQ((__m256i*)(c2+j),_mm256_sub_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_sub_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_sub_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_sub_epi64(b3,v3)); }
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
__m256i a = LDQ((const __m256i *)(c1 + j));
STQ((__m256i *)(c1 + j), ((op1 == 1) != (s1 < 0)) ? _mm256_add_epi64(a,v)
: _mm256_sub_epi64(a,v));
__m256i b = LDQ((const __m256i *)(c2 + j));
STQ((__m256i *)(c2 + j), ((op2 == 1) != (s2 < 0)) ? _mm256_add_epi64(b,v)
: _mm256_sub_epi64(b,v));
}
for (; j < h; j++) {
long long v = t[j];
if ((op1 == 1) != (s1 < 0)) c1[j] += v; else c1[j] -= v;
if ((op2 == 1) != (s2 < 0)) c2[j] += v; else c2[j] -= v;
}
}
}



/* row-blocked leaflevel: the A-side build for rows [ic0,ic1) of all 7 panels is
done immediately before those rows are consumed.  Every operand row is still
read exactly once, the instruction stream of each phase is unchanged.
The per-product destination operators (oa/ob) are functions of the ORIGINAL
sequential wr[] evolution, so they are computed ONCE here, before the loop. */
static void leaflevel3(int n, const Opd *Ao, const Opd *Bo, long long *d, int ld, int op,
long long *T) {
const int h = n >> 1;
long long *Q[4];
Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
int OA[7], OB[7];
{
int wr[4] = {0, 0, 0, 0};
const int qi[7] = {0, 2, 1, 0, 0, 3, 0};
const int si[7] = {1, 1, 1, 1, 2, 1, 1};
const int qj[7] = {3, 3, 3, 2, 1, -1, -1};
const int sj[7] = {1, 2, 1, 1, 1, 0, 0};
for (int L = 0; L < 7; L++) {
OA[L] = wr[qi[L]] ? effop(1, si[L]) : effop(op, si[L]);
wr[qi[L]] = 1;
if (qj[L] >= 0) { OB[L] = wr[qj[L]] ? effop(1, sj[L]) : effop(op, sj[L]); wr[qj[L]] = 1; }
}
}
buildB(Bo);
/* ROW-BLOCKED A BUILD (l1k7_rv `rba`): the A panel depends on the ROW index and
NOT on the column group, and buildA_blk reads every operand row exactly ONCE
across the blocks, so blocking the build by MR rows costs no extra read -- while
the live A window drops from 7 x 64 x 608 B = 268 KB to 7 x MR x 512 B = 21.5 KB.
B stays a single up-front pass (its k is the contraction index and every mikro
call reads all KC of it).  Only the ORDER of the mikro calls changes. */
{
static const int DQ1[7] = {0, 2, 1, 0, 0, 3, 0};
static const int DQ2[7] = {3, 3, 3, 2, 1, -1, -1};
static const int SIx[7] = {1, 1, 1, 1, 2, 1, 1};
static const int SJx[7] = {1, 2, 1, 1, 1, 0, 0};
for (int ic0 = 0; ic0 < h; ic0 += MR) {
int ic1 = ic0 + MR; if (ic1 > h) ic1 = h;
buildA_blk(Ao, ic0, ic1);
if (op == 0) {
/* FUS3 at the leaf: M1->Q[0], M2->Q[2], M3->Q[1], then ONE C22 pass. */
leafP_blk(h, 0, ic0, ic1, Q[0], ld, 0, 0, 0, 0);
leafP_blk(h, 1, ic0, ic1, Q[2], ld, 0, 0, 0, 0);
leafP_blk(h, 2, ic0, ic1, Q[1], ld, 0, 0, 0, 0);
cmb3_blk(Q[3], ld, Q[0], Q[2], Q[1], ld, h, ic0, ic1);
for (int L = 3; L < 7; L++) {
long long *D1 = Q[DQ1[L]];
long long *D2 = (DQ2[L] >= 0) ? Q[DQ2[L]] : 0;
const int oa = OA[L], ob = OB[L];
if (D2 == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, oa, 0, 0, 0);
} else if (oa == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, 0, 0, 0, 0);
cmb_blk(D2, ld, D1, ld, h, SJx[L], ob, ic0, ic1);
} else {
leafP_blk(h, L, ic0, ic1, T, h, 0, 0, 0, 0);
cmb2_blk(D1, ld, oa, SIx[L], D2, ld, ob, SJx[L], T, h, h, ic0, ic1);
}
}
} else {
for (int L = 0; L < 7; L++) {
long long *D1 = Q[DQ1[L]];
long long *D2 = (DQ2[L] >= 0) ? Q[DQ2[L]] : 0;
const int oa = OA[L], ob = OB[L];
if (D2 == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, oa, 0, 0, 0);
} else if (oa == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, 0, 0, 0, 0);
cmb_blk(D2, ld, D1, ld, h, SJx[L], ob, ic0, ic1);
} else {
leafP_blk(h, L, ic0, ic1, T, h, 0, 0, 0, 0);
cmb2_blk(D1, ld, oa, SIx[L], D2, ld, ob, SJx[L], T, h, h, ic0, ic1);
}
}
}
}
}
}

static void mmxR(int n, const Opd *Ao0, const Opd *Bo0, long long *d, int ld, int op,
long long *T) {
const Opd *Ao = Ao0, *Bo = Bo0; Opd A1, B1;
if (n <= CUTL || (n & 1) || Ao->nt > 4 || Bo->nt > 4) {
leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); return; }
if (n == 2 * CUTL) {
{ leaflevel3(n, Ao, Bo, d, ld, op, T); return; }
leaflevel(n, Ao, Bo, d, ld, op, T); return; }
if (n <= 8 * CUTL && (Ao->nt > 1 || Bo->nt > 1)) {
int Lm = (n > 4 * CUTL) ? 0 : (n > 2 * CUTL) ? 1 : 2;
if (Ao->nt > 1) { matop(Ao, n, matbufA[Lm]); A1 = o1(matbufA[Lm], n); Ao = &A1; }
if (Bo->nt > 1) { matop(Bo, n, matbufB[Lm]); B1 = o1(matbufB[Lm], n); Bo = &B1; }
}
const int h = n >> 1;
Opd A11 = oq(Ao, 0, 0, h), A12 = oq(Ao, 0, 1, h), A21 = oq(Ao, 1, 0, h), A22 = oq(Ao, 1, 1, h);
Opd B11 = oq(Bo, 0, 0, h), B12 = oq(Bo, 0, 1, h), B21 = oq(Bo, 1, 0, h), B22 = oq(Bo, 1, 1, h);
long long *Q[4];
Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
long long *T2 = T + (size_t)h * h;
int wr[4] = {0, 0, 0, 0};
#define MK(qi, si, qj, sj, Ma, Mb) do { \
int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
int ob = wr[qj] ? effop(1, (sj)) : effop(op, (sj)); \
wr[qi] = 1; wr[qj] = 1; \
if (oa == 0) { \
mmxR(h, &(Ma), &(Mb), Q[qi], ld, 0, T); \
cmb(Q[qj], ld, Q[qi], ld, h, (sj), ob); \
} else { \
mmxR(h, &(Ma), &(Mb), T, h, 0, T2); \
cmb2(Q[qi], ld, oa, (si), Q[qj], ld, ob, (sj), T, h, h); \
} \
} while (0) \

#define MK1(qi, si, Ma, Mb) do { \
int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
wr[qi] = 1; \
mmxR(h, &(Ma), &(Mb), Q[qi], ld, oa, T); \
} while (0)
if (op == 0) {
/* FUS3: M1, M2, M3 all write C22.  Their three sources (C11, C21, C12) are
live simultaneously, so one pass replaces three.  Same mmxR call order. */
wr[0] = 1; wr[3] = 1;
{ Opd Ma = oadd(&A11, &A22, 1), Mb = oadd(&B11, &B22, 1);
mmxR(h, &Ma, &Mb, Q[0], ld, 0, T); }
wr[2] = 1;
{ Opd Ma = oadd(&A21, &A22, 1), Mb = B11;
mmxR(h, &Ma, &Mb, Q[2], ld, 0, T); }
wr[1] = 1;
{ Opd Ma = A11, Mb = oadd(&B12, &B22, -1);
mmxR(h, &Ma, &Mb, Q[1], ld, 0, T); }
cmb3(Q[3], ld, Q[0], Q[2], Q[1], ld, h, h);
} else {
{ /* M1 = (A11+A22)(B11+B22) -> +C11, +C22 */
Opd Ma = oadd(&A11, &A22, 1), Mb = oadd(&B11, &B22, 1);
MK(0, 1, 3, 1, Ma, Mb);
}
{ /* M2 = (A21+A22)B11 -> +C21, -C22 */
Opd Ma = oadd(&A21, &A22, 1), Mb = B11;
MK(2, 1, 3, 2, Ma, Mb);
}
{ /* M3 = A11(B12-B22) -> +C12, +C22 */
Opd Ma = A11, Mb = oadd(&B12, &B22, -1);
MK(1, 1, 3, 1, Ma, Mb);
}
}
{ /* M4 = A22(B21-B11) -> +C11, +C21 */
Opd Ma = A22, Mb = oadd(&B21, &B11, -1);
MK(0, 1, 2, 1, Ma, Mb);
}
{ /* M5 = (A11+A12)B22 -> -C11, +C12 */
Opd Ma = oadd(&A11, &A12, 1), Mb = B22;
MK(0, 2, 1, 1, Ma, Mb);
}
{ /* M6 = (A21-A11)(B11+B12) -> +C22 */
Opd Ma = oadd(&A21, &A11, -1), Mb = oadd(&B11, &B12, 1);
MK1(3, 1, Ma, Mb);
}
{ /* M7 = (A12-A22)(B21+B22) -> +C11 */
Opd Ma = oadd(&A12, &A22, -1), Mb = oadd(&B21, &B22, 1);
MK1(0, 1, Ma, Mb);
}
#undef MK
#undef MK1
}

void matrix_multiply(int n, const long long *A, const long long *B, long long *C) {
Opd Ao = o1(A, n), Bo = o1(B, n);
if (n <= CUTL || (n & 1)) { leaf(n, &Ao, &Bo, C, n, 0, 0, 0, 0); return; }
mmxR(n, &Ao, &Bo, C, n, 0, g_scratch);
}
#pragma GCC pop_options

评测结果