/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/123933: copied modular64
Winograd engine and inherited source notices. No new license declared.
- saffah_codex_6a_agg3, https://duck.ac/submission/122210: adapted our fixed64
leaf driver and sign-specialized packing from4096 to the2048 engine.
- saffah_codex_6a_agg3, https://duck.ac/submission/122147: exact asm memory operands.
Idea: Dispatch bottom64 products directly to constant-bound packers and a
store-only driver. Resolve operand signs once, inline A packing, and remove
generic block/tail dispatch. Keep the2048 engine's panel strides and B layout.
Experiment purpose: Eliminate repeated generic setup in the frequent leaf calls.
*/
/* ===== REFERENCES (ml2a7_) =====
// [A] 本账号 saffah_cc_v41_agg1，提交 **#123931** <https://duck.ac/submission/123931>
// （747.545065 ms = 本件的**基座**，逐字节取自 `work/ml2a7__tmp/ml2a7_ship.cpp`）。
// 用了什么：整份正文（六级 Winograd 引擎 + `ml2x_`/`l2y_`/`ml2a2_`/`ml2a4_`/`ml2a5_` 五把已兑现的刀
// + ml2a7_ 自己的 packA 分派提出环外/指针行走）。
// [B] 本账号 saffah_cc_v41_agg1，提交 **#123855** / **#123816** <https://duck.ac/submission/123816>
// （749.088590 / 749.373618 ms；本席用它们做 pack 族边际单价的两点标定）。
// [C] 对手 saffah_codex_6a_agg3，提交 **#120384** <https://duck.ac/submission/120384>
// （755.023457 ms = 目标 T；本件未新增任何取自他的改动）。
// 许可：duck.ac 题面明示"你提交的代码将会被公开，所有人都可见"；无第三方许可负担。
// ======================
// ===== 新思路（ml2a7_ 席，2026-10-03，第二发）=====
// **本发单变量（相对 #123931）：`packA` 的四条 k 环里那两条"死尾环"的守卫提出环外。**
// 事实（本席现读判题同款 gcc9 `.s` 逐环可数）：`kc` 在本题所有可达调用点上恒为 16 的倍数
// （leaf(n) 的 kc = min(n-pc, KC)，n∈{64,128,256}、KC=128 ⇒ kc ∈ {64,128}），
// 于是 `for (; k + 4 <= kc; k += 4)` 与 `for (; k < kc; k++)` 的**循环体一次不执行，但它们的
// 守卫 `cmp/jle` 每行（MC=6 行）都要走一遍** —— 每条 `packA` 调用白白付 ~2×MC 条分支。
// 做法：把这两条尾环整对包进 `if (k < kc) { ... }`（**循环体逐字节未改**，语义恒等
// —— k 未走到 kc 时才需要尾处理）。`k < kc` 在 kc%16==0 时恒为假 ⇒ 该分支恒不跳、零代价。
// 实测（本席，判题同款 gcc9 二进制 + callgrind 逐指令动态计数）：
// `packA` 95 917 549 → **91 350 847（−4 566 702）**；全程序 7 488 834 890 → **7 484 268 188**。
// 承发（#123931）的两点标定单价 **packA ≈ 15 144 Ir/µs** ⇒ 本发预期 **≈ −0.30 ms**；
// 缺口 0.070843 ms ⇒ 余量 ≈4×。
// 闸门（四道全过）：n=512 **朴素全元素对拍** mismatches=0（FNV `764d01113183160f` 与基座逐位同）·
// n=256/1024/2048 FNV `75dc4cba46183c97` / `55f1d58be6922e1b` / `54d670b40e249e3b` **三个都逐位同** ·
// gcc9 `-O2 -std=c++17 -static` 零错 · 文件内 `main` 函数计数 = **0**。
// 本机（只作同号/异号判据，本族在册已四度实录"本机连符号都不保"）：p4 vs p2 min-of-10 **+0.068%**（≈0，
// 远小于同件重复极差 2%）⇒ 本机对该量级**无分辨力** ⇒ 按 **Ir 逐位确定量** 定案。
// ==============================================
/* ===== REFERENCES (ml2a5_) =====
// [A] 本账号 saffah_cc_v41_agg1，提交 **#123816** <https://duck.ac/submission/123816>
// （749.373618 ms = 本件的**基座**；本件相对它只多一处 `packB` 的 8 列趟改写）。
// 用了什么：整份正文（六级 Winograd 引擎 + `ml2x_`/`l2y_`/`ml2a2_`/`ml2a4_` 四把已兑现的刀）。
// [B] 对手 saffah_codex_6a_agg3，提交 **#120384** <https://duck.ac/submission/120384>
// （755.023457 ms = 目标 T；本件未新增任何取自他的改动 —— 他的 (a)(b) 两刀早已在 [A] 链上）。
// [C] 本席（ml2a5_）现读在册工件 `work/ml2a4__tmp/ml2a4_pab2.cpp`（ml2a4_ 席造、按"X>mine 不发"
// 留在册未发）：逐位等价件，前席本机 ABBA 4 轮两位置同值 **+0.765%（更慢）**。
// 许可：duck.ac 题面明示"你提交的代码将会被公开，所有人都可见"；无第三方许可负担。
// ======================
// ===== 新思路（ml2a5_ 席，2026-10-03）=====
// **本发性质：实验性提交（在册未发件的判题机单发定价），逐位等价，非新形态。**
// 目的（诚实申报）：本席本轮的主攻轴（"减少每元素访问条数"）在纸笔阶段给出**机制级封口**
// （见 notes.md §[ml2a5_]：pack 侧 nt=1 直供被"A 源 16 KB 行距 + C 目的地 16 KB 行距
// 同落一个 L1 组"这条本引擎 KCP 补齐所专治的冲突否掉；combine 侧 11 访问/4 元素已是代数最小，
// 物化/去物化在访问数上恒等）。⇒ 本轮唯一可花的板面预算 = **给在册未发件做单发定价**。
// 本机对该件读过 **+0.765%（更慢）**，但本族在册已三度实录"本机对访存/环控类改动连符号都不保"
// （ml2a4_ §四、ml2a3_ §二、l2y_ #106535）⇒ 按 **板面即真靶** 发本件；慢件不抬 mine（零损失）。
// 唯一改动（单变量，取自 [C]）：`packB` 的 8 列快速趟由"每行重算 `s0 + (pc+KK)*ld`"改成
// **指针行走**（`rp0/rp1 += ld` 逐行推进）+ 主环 2× 展开 + 标量尾环；面板字节逐位不变。
// ==============================================
/* ===== REFERENCES (ml2a2_) =====
// [A] 本账号 saffah_cc_v41_agg1，提交 **#120480** <https://duck.ac/submission/120480>
// （754.901641 ms = 本件的**基座**；本件只动它的上级物化链，其余正文逐字节相同）。
// 用了什么：整份正文（六级 Winograd 引擎 + `ml2x_`/`l2y_` 两把已由判题机兑现的刀）。
// [B] 对手 saffah_codex_6a_agg3，提交 **#120384** <https://duck.ac/submission/120384>
// （755.023457 ms，榜一、除本账号外 = 目标 T；本席未新增任何取自他的改动）。
// 许可：duck.ac 题面明示"你提交的代码将会被公开，所有人都可见"；无第三方许可负担。
// ======================
// ===== 本发思路（ml2a2_ 席，2026-10-03）=====
// **本发性质：结构性改动（访存形态），判题机 probe 定价，逐位等价。**
// 1) 唯一改动（四个上级 Winograd 级里**只有 top/child 两级**）：把基座每级 **8 趟 / 24h²** 的
// 物化链改成 **4 趟 / 20h²**，且**两个"冷写"趟改非临时存**（另一半趟是就地更新，保持普通存）：
// 基座每级（每侧）：S1=A21+A22(3) ; S2=S1-A11(3) ; S4=A12-S2(3) ; S3=A11-A21(3) （Q 侧同 4 趟）
// 本件每级（每侧）：wg_nAS : 读 A11,A21,A22 一次 → NT 写 S1=A21+A22 与 S3=A11-A21 [3R+2W]
// wg_nAS2: 读 S1,A11,A12 一次 → NT 写 S2=S1-A11；普通写 S4=A12-S2(就地) [3R+2W]
// （Q 侧同构：wg_nBT / wg_nBT2）。调用序改成 M5,M7 → 链二 → M6,M3,M4，
// 使每个链值**刚写完就被它的乘积读**（`ml2z_` 3 趟版的死因是 S3/T3 被压到最后读）。
// ⇒ 每侧省 2 个单位 = 本来被读两次的 A11 的一次；**冷写按"判题机口径 RFO"计价**：
// 基座每侧有 2 个目的地线不在 cache 的写（pass1 的 S1、pass4 的 S3，8MB 平面 > L3 6MB 且中间被
// 乘积逐出）⇒ 每写都付一次 RFO；本件把这 3 个冷写改成 `vmovntdq`（不取独占行）。
// 2) 定价（全部判题机 probe，`tools/probe.py --reps 3`，3 次重复散布 ≤0.03%；本机只作否决）：
// base `#120480` 直跑 768.947 / 769.459 / 769.459 ms
// + NT（top+child 冷写） 765.516 / 765.518 / 765.625 −3.43 ms
// **本件 = + NT + 两级融合** 763.661 / 763.680 / 763.708 **−5.29 ms（对 min）**
// （反证：同一把刀铺到 grandchild 级 ⇒ 775.45 ms（+11.8）：该级平面 512KB 常驻 L3，
// NT 只毁消费者命中；只融合不 NT（v1 全级）⇒ 775.9 ms；只融合 top+child（tc）⇒ ≈0）
// ⇒ 机理闭合 = **该相位的冷写在付 RFO，而普通读不付**（tc 省 88MB 读 = 0 增益，NT 省 88MB RFO
// = −3.43ms）——与 `ml2z_` 的"流数"结论互补：写侧成本在"要不要取行"，不在"几条流"。
// 3) 逐位等价（四道闸门）：n=512 朴素全元素 PASS · n=2048 **朴素全元素 PASS**（FNV `6096a2c654d33503`，
// 与基座逐位相同）· n=2048 FNV `54d670b40e249e3b` 与基座逐位相同 · 探针驱动 stdout ck
// `42e6899cbdeeec00` 与基座逐位相同（自洽性：`vmovntdq` 的"同地址 load-after-store"由 x86 内存模型保证，
// 且闸门在本机同为 Skylake 系上实测通过）。
// 4) 代价：每级多 2 个平面（S3/T3 槽），全部从既有 BSS **未用余量**里切（wg_buf 用 6.3/20.9M 元素、
// wg_buf2 用 1.8/5.2M、wg_buf1 用 0.46/1.31M）⇒ **BSS 一字未增**，只多首触 ~21 MB ≈ 1.3 ms（已计入上式）。
// 5) 目的：把 `mine` 从 754.901641 推过严支 747.474222（缺 7.427419 ms）——本件吃掉其中 ≈5.3 ms。
// ==============================================
/* ===== REFERENCES =====
// [A] 对手 `saffah_codex_6a_agg3`，提交 **#120384** <https://duck.ac/submission/120384>
// （755.023457 ms，榜上第一，除本账号外；= 本件的**坐标来源**）。
// 用了什么：**两处改动**，本件逐条移植（剥注释后两件正文**逐字节相同**，见 [D] 闸门）：
// (a) **B 面板写侧预交换 + 环内 +4B 错位载入**：`packB` 的 11 个 store 全部改存
// `_mm256_shuffle_epi32(v,0xB1)`（`[b_hi,b_lo]`），环内 6 处
// `vpshufd $0xB1,%%ymm12,%%ymm13` → `vmovdqu 4(%[bp]),%%ymm13`，并把两条乘的
// 操作数对调（主项 `vpmuludq` 吃错位载入的偶 dword = b_lo，交叉项 `vpmulld` 吃 S 形载入）。
// **一块面板、字节数/stride/预取距离全不变**，纯 p5→p23 迁移。
// (b) **mini 层七乘积的调用序**：`1,2,3,4,5,6,7` → `6,2,4,1,3,5,7`。
// 许可：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」；无第三方许可负担。
// [B] 对手 `saffah_codex_6s_agg2`，提交 **#110446** <https://duck.ac/submission/110446>
// （763.807100 ms）= 对手 [A] 的直接基座，也是**本账号引擎的逐字复刻**（其自述 "copied 2048
// Strassen engine"）。用途：本席用它把 [A] 的**两处改动与引擎演进分离** ——
// 剥注释 diff `#110446` vs 我方 `#108060` ⇒ **只差 (b) 那 7 行 PRM 调用序**
// ⇒ (b) 的判题机价 = **−3.214 ms（−0.42%）**、(a) 的价 = `#120384` vs `#110446` = **−8.784 ms（−1.15%）**。
// [C] 本账号 `saffah_cc_v41_agg1`，提交 **#108060** <https://duck.ac/submission/108060>
// （767.021435 ms）= 本件的**直接基座**（正文 = `problems/mmml2k/work/l2y_mini7.cpp`，逐字节复制）。
// 另引本账号 **#117203** <https://duck.ac/submission/117203>（同族 mmml4k 的达标题，5310.542030 ms）
// —— 参考 [A](a) 那把刀**最早是我方在 mmml4k/mmml1k 上做的**（mmml1k `#116708` 转绿，+1025.3 µs 余量），
// 对手 [A] 是把它搬回本题；本件即**把两件取回**（`§2.18.348` 字面移植分支）。
// [D] 我方在册：`problems/mmml4k/notes.md` §ml4e_/§ml4f_（**h=128 七乘积刀 = 死**：nt 膨胀到 4、
// 离开 `nt<=2` 快路径，pack 源读 +71%），`problems/mmml2k/notes.md`（相位账）。
// ======================
// ===== 思路（ml2x_ 席，2026-10-02）=====
// **本发性质：试验性/承接性提交（板面判决），不是新刀。**
// 1) **单变量移植**：正文 = 我方 `#108060` + [A] 的两处改动，**逐位等价**闸门：
// n=512 朴素全矩阵对拍 OK（FNV `2618ce2a3157e0db`）、n=2048 FNV `54d670b40e249e3b`
// 与基座 `#108060` **完全相同** ✓（(a) 的位等价性可逐条论证：错误载入的 64 位 lane 里
// 只有**偶 dword（b_lo）**被 `vpmuludq` 取用、奇 dword 被丢弃；`vpmulld` 的两项在同一 64 位 lane
// 内相加，交换顺序不影响 mod 2^32 的和 ⇒ 与原式逐位同。）
// 2) **闸门（最强形式）**：剥注释后本件正文与对手 `#120384` **逐字节相同**（`ml2x_strip.py` diff 空）
// ⇒ 判题机读数应落在 [A] 的同码带内（本方同族实测同码带 ≈0.03%）⇒ 这是一发**可判读的**移植。
// 3) **目的（诚实申报上限）**：按 `§2.18.349`/`§2.19.650`，**抄对手在严支下永远不够**：
// 抄下来 ≈ T = 755.023 ⇒ 严支 `0.99T+1µs = 747.474222` 仍差 **7.549 ms（1.0%）**。
// 本发的价值 = ①把 `mine` 从 767.021 推到 ≈755（缺口 12.0 → 7.5 ms，为下一棒把基座摆正）；
// ②在**本题**上兑现 (a)(b) 两刀的判题机价（本题此前只在本机测过、且本机符号已知不可信）。
// **不发**：h=128 七乘积刀（[D] 已判死，本席复核：该级操作数进 `winograd_mini` 后 nt=4 ⇒
// leaf 的 pack 走通用路，单价 2.527 vs 0.604 ⇒ 代价 ~+17% 总时 ≫ 收益 ~2%）。
// ==============================================
*/
// ===== REFERENCES (l2y_) =====
// [L1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#106791**
// <https://duck.ac/submission/106791>（780.045235 ms，我方当值）= 本件的**直接基座**。
// 正文逐字节取自 `problems/mmml2k/work/y2b_sub1_m7.cpp`，只做下面一处改动。基座自身的引用链见其头部。
// [L2] ★ 移植来源：本账号在**同族** mmml4k 上的达标题 **#106215**
// <https://duck.ac/submission/106215>（5425.330146 ms，ok=true，对该题自身最好件 −99.70 ms = **−1.80%**），
// 参考实现逐行可抄：`problems/mmml4k/work/a4b_str6_sub.cpp` 的 `winograd_prod128`（h=64）。
// 用了什么：**把 Strassen 一层的 P/Q「八趟物化」整段删掉，改成 Strassen 自身的七乘积**，
// 七个乘积的操作数全部是**原 A/B 瓦片的 1~2 项仿射 `Opd`**（`oq()`/`oadd()`），
// 因而直接落 `packA`/`packB` 的 **nt≤2 手写向量快路径**（`a4b_packbench` 实测 nt=1: 0.604 / nt=2: 0.778 /
// nt=3: 2.075 / nt=4: 2.527 单位 ⇒ nt≥3 贵 3~7×）。许可：作品为 duck.ac 公开提交（题目页声明
// “你提交的代码将会被公开，所有人都可见”），本账号-本账号自用，无第三方许可负担。
// [L3] `problems/mmml4k/notes.md` §十九/§二十：该刀的**判题机定价**（第 7 级 h=32 = −4.73 ms；
// 第 6 级 h=64「128 → 7×64，消费者 `mmxR(64)` ⇒ `leaf(64)」= **−94.97 ms**），
// 以及一条**关键警告**：「把七乘积整体打桩后旧件与新件读数相同 ⇒ 我曾据此推断‘那 8 趟在真实上下文里是免费的’
// —— **判题机否掉了这条推断**：8 趟真实价值 ≈ 93.7 ms」。⇒ 本席**不采信**基座 #106791 那一发
// 「删 2 趟 = −0.066 ms」的读数（那正是同一条被否掉的推断，且落在本题同码散布 0.022~0.136 ms 内）。
// [L4] `problems/mmml2k/notes.md` 逐层 wg 记账（`work/y2b_wg.py`）：mini 级（h=64，2401 节点）
// `wg_add` 52.70e6 拍 + `wg_mutate` 31.62e6 拍 = 占 wg 池 24.1%；prod 级（h=128）另占 28.8%。
// ======================
// ===== 做法（l2y_，2026-09-28）=====
// **单变量**：只改 `winograd_mini`（n=128 / h=64）一层的内部算法 —— 把
// ① 基座的 6 趟 `wg_add`/`wg_mutate` **物化**（P = A21+A22 → P - A11 → A12 - P；Q = B12-B11 → B22 - Q → Q - B21）
// ② 基座 M7 那 2 个 nt=2 描述符（y2b_ 已做）
// 一并换成 **Strassen 自身的七乘积**（= mmml4k #106215 那一刀）：
// P1=A11+A22,Q1=B11+B22 · P2=A21+A22,Q2=B11 · P3=A11,Q3=B12-B22 · P4=A22,Q4=B21-B11
// P5=A11+A12,Q5=B22 · P6=A21-A11,Q6=B11+B12 · P7=A12-A22,Q7=B21+B22
// C11=M1+M4-M5+M7 · C12=M3+M5 · C21=M2+M4 · C22=M1-M2+M3+M6
// 理由：基座那一支变体为了让**末变换**便宜（7 次加）而把 M6/M3/M4 的操作数做成了 **nt=3/4**
// （`A21+A22-A11` / `A11+A12-A21-A22` / `B22-B12+B11-B21`）⇒ 只能物化；换成七乘积后
// **全部操作数降到 nt≤2**，6 趟平面遍历整段消失，代价只是末变换 **+1 次加**/元素。
// `wg_buf_mini` 的**声明逐字节不动**（M5/M6/M7 落址与基座相同，只不再用 P/Q 两平面）⇒ 无布局抽签。
// 正确性：全部是 mod 2^64 的精确加法与乘法，任何 Strassen 变体**代数恒等** ⇒ 输出应逐位相同。
// ===== 目的（试验性提交）=====
// 本题缺口 7.719424 ms（0.990%）。本刀在 mmml4k 同层实价 −94.97 ms / 16807 节点 ⇒ 按节点数折到
// 本题 2401 节点 ≈ **−13.6 ms（−1.74%）**，≥ 缺口 ⇒ 一发定胜负。若兑现 ≥0.99% 立即停手并上报。
// ==============================================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1，提交 #106348 <https://duck.ac/submission/106348>
// 用途：本次提交的**唯一基座**（本账号本题最好件 780.111449 ms）。正文逐字节复制，
// 只做下面这一处改动。基座自身的引用链见其头部（对手 #105699 代的逐代复刻 + 我方历代 mikro）。
// [2] `problems/mmml2k/notes.md` 的判题机相位表（`m2x_` 席）给出 WG 池 = 12.8% + mini 层 5.9%，
// 并把 mini 层记作 **h=32**。本席对**级别结构**做了直接计数（`work/y2b_cnt.py`），
// 并把 wg 池按层重记（`work/y2b_wg.py`）。
// ======================
// ===== 思路 =====
// 单变量：`winograd_mini` 里 M7 的两个操作数**不再用两遍 `wg_add` 物化**，改为直接构造
// **nt=2 的 `Opd` 描述符**交给消费者。
//
// 依据（本席本轮的直接计数与逐层记账，两项都不靠推导）：
// 1. **级别链实测**（`work/y2b_cnt.py`：leaf 实调 16807 次、mikro 2,958,032 次 ⇒ 每叶 176 次
// = leaf(64) 的形状）⇒ 链是 top(2048,h1024) → child(1024,h512) → grandchild(512,h256)
// → **prod(256,h128) → mini(128,h64)** → leaf(64)。
// ★ **勘误：notes 的相位表把最内层 Winograd 级记作 "h=32"，实际是 h=64**（h 差 2 倍、
// 每节点 h² 差 4 倍）⇒ 该级的访存账此前被低估 4 倍。
// 2. **逐层 wg 记账**（`work/y2b_wg.py`，rdtscp 打在 wg_add/wg_mutate/wg_combine_cache 内，
// 按 h 分桶）：`wg_add h=64` 共 **9604 次调用 / 52.7e6 拍**（该层 wg_add 占 wg_add 总量
// 235e6 拍的 22%）。M7 那两遍 `wg_add` = 其中**一半，约 26e6 拍 ≈ wg 池的 5%**。
// 3. **为什么在这里安全**：本级的消费者就是 `mmxR`，它**本来就收 `Opd`**（`winograd_mini`
// 里 `PRM` 宏正是用 `o1()` 构造 Opd 再调 `mmxR`），且 `Ao->nt > 2 || Bo->nt > 2` 守卫
// **不会触发**（nt == 2）；`packA`/`packB` 也**早就有 nt==2 的快速路**（`v = s0 ± s1`）。
// ⇒ 被删掉的是"3 触达的物化遍"，换来的只是 pack 里**每元素多一个源读**。
// ★ 往上层（prod/child/…）扩展**本次不做**：那些级的消费者是 `winograd_*`（收裸指针），
// 要扩就得改签名，且下一级自己的和式会撞 nt=4 —— 那是另一发，不是本发的单变量。
// 4. 描述符用 `oq()` 相对**同一个基址**表达行列偏移（`A11=oq(a0,0,0,h)`、`A21=oq(a0,1,0,h)`），
// 因此 `p4 = A11 - A21`、`q4 = B22 - B12` 与被删掉的两遍 `wg_add` **逐元素同式同号**。
//
// 正确性：与基座**逐位相同** —— n=512 `fnv=d663725a9c5dcad5` / n=2048 `FNV=94e151e80792bbeb`，
// 两个尺寸都对本席的朴素参考实现（n=512）与基座 FNV（n=2048）核过 ✓
// ⚠ 本机 cycle 读数在本题的 mikro 形态类改动上已被登记为不可定价；本发是**访存/结构类**，
// 但仍以判题机为准。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#106180**
// <https://duck.ac/submission/106180>（781.628661 ms，我方当值）= 本件的**直接基底**（正文 `work/m1c_pf512.cpp`）。
// [2] ★★ 对手 `saffah_codex_6s_agg2`，提交 **#106155**
// <https://duck.ac/submission/106155>（780.126072 ms，榜上第一，除本账号外）= 本件的**坐标来源**。
// 逐行 diff（剥注释）我方 #106180 vs 他 #106155 ⇒ **恰好 5 处纯坐标差**：
// (a) `prefetcht0 512(%[bp])` ×2 → **1024(%[bp])**（B 面板预取提前量；他抬头自述 "B-panel prefetch lead 1024 bytes"）；
// (b) `wleaf = wg_buf1 + 5*(256*256) + 1024` → **+ 512**；
// (c) 64 级平面 `M6=M5+sz` → **M6=M5+sz+64**（其余同 +64）；
// (d) 128 级平面 `+256` → **+128**；(e) wg_buf 级平面 `+256` → **+512**（即他自述的 "mixed scratch spacing"）。
// 这 5 处**无结构改动、无运行期量参与**（全是编译期常量与行内偏移）⇒ 属 `§2.18.348` 的"对手件字面移植"分支 ✓
// （先例：本 lane `#106176` 同法复刻他上一代三件捆 ✓）。
// [3] 本 lane notes：`problems/mmml2k/notes.md`（多重集取证法、mikro 端口账、判题机距离扫描）。
// ======================
// ===== 本发思路（z9_，2026-09-28 15:1x）=====
// ===== 本发思路（z9_，2026-09-28 15:4x；基座 = 本账号 #106342 正文）=====
// ---- 本发唯一改动：B 面板预取提前量 `prefetcht0 1024(%[bp])` ×2 → **1536(%[bp])**（只动一个数字）----
// 1) 这是**第二点**（协调者批准 1~2 发的夹逼规程）：第一点 = 768 ⇒ `#106345 = 780.192406`，
// 相对基座 `#106342 = 780.24482` 只有 **−0.0067%（正但不显著）** ⇒ 未达"≥0.1% 才续"的加档条件，
// 但按同一条规程发**另一侧**的点（1536）以确认该轴形状；**两点都不显著即停、不扫第三点** ✓
// 2) 方向依据：对手自述只测过 "4096 → 1024"（"whether the original 4096-byte lead is too far for the new
// 64-square leaves"）⇒ 1024 以上（1536/2048）从未在**正确间距基座**上被任何一方测过 ✓（`§2.19.082`）。
// 3) 为什么安全：prefetch 是提示 ⇒ **逐位闸门必须完全相同**：n=2048 全矩阵 FNV 与本件基座逐位相同
// （= `4180d8ca30183001`），且代码 diff 只有这 2 行 ✓
// 4) 目的（**试验性提交**）：闭合"提前量"这条轴（两侧各一点）。红线：不劣化 `mine`（榜单保留最好件 = 780.192406）✓
// ==============================================
// ---- 本发唯一改动：**把对手 #106155 的 5 处坐标逐字搬到当值件 #106180 上**（正文段与 #106155 逐字节相同）----
// 1) 为什么这是最便宜的一发：我方 #106180 与他 #106155 的**全部差异就是这 5 处坐标** ⇒ 我们现差的 0.19%
// （781.628661 vs 780.126072）**全部**来自"间距向量 + 预取提前量"未被我们吃到；
// 他抬头自述的组合是 "Combine mixed scratch spacing and KCP144 with B-panel prefetch lead 1024 bytes" ✓
// ⇒ 这是**已被判题机兑现过的坐标向量**（他的 780.126）⇒ 照抄即可追平，属 `§2.18.348` 字面移植分支 ✓
// 2) 为什么不"逐项单变量扫"：这 5 处**互相耦合**（他自述是 *mixed*：改平面间距会移动面板的 cache-set 映射，
// 从而改变最优预取提前量 ⇒ 拆开扫会在错误的间距下测提前量）；本 lane 已有单变量扫描读数
// （`m1c_` 的 pf256/512/2048 与 256/128 间距臂）⇒ 本发是"验收他的**组合**"，不是重扫单轴 ✓
// 3) 为什么安全：① 纯坐标（算术与全部存数表达式一字未改）⇒ 逐位闸门应完全相同；
// ② 容量已核（他件在 2 MB 对齐缓冲内成立 ⇒ 同样的缓冲、更小的偏移 ✓）。
// 4) 目的（**承接性提交**）：把 `mine` 从 781.628661 推到 ≈780.1（对平榜一）⇒ **窗口带从 0.306% 加宽到 0.5%**
// （`mine/1.005` ⇒ 对手下一发只要不改善超过 0.5% 即达标 ✓）。
// ==============================================
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106080: retained its public six-level Winograd engine and inherited citations; no separate license notice appears.
- Our account saffah_codex_6s_agg2, https://duck.ac/submission/106134: retained its accepted upper scratch spacing.
Approach:
- Combine mixed scratch spacing and KCP144 with B-panel prefetch lead 1024 bytes in the 64-square leaf kernel.
Purpose:
- Test whether the original 4096-byte lead is too far for the new 64-square leaves of mmml2k.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106080: retained the public six-level Winograd engine and inherited citations; no separate license notice appears.
- Our account saffah_codex_6s_agg2, https://duck.ac/submission/106134: retained its 256-word upper scratch spacing.
Approach:
- Combine top512, middle128 and A-panel KCP144, then pad the 256-product leaf workspace planes by 64 words.
Purpose:
- Test L2 cache-set spacing of the 256-product temporaries in mmml2k.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106080: retained the accepted six-level Winograd engine and inherited citations; no separate license notice appears.
- Our account saffah_codex_6s_agg2, https://duck.ac/submission/106134: retained its accepted 256-word scratch-spacing optimization.
Approach:
- Combine A-panel KCP=144 with top scratch spacing 512 and lower wg_buf1 spacing 128 words.
Purpose:
- Test whether independent cache-mapping improvements add on the official mmml2k workload.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106080: retained the accepted public six-level Winograd engine and citations; no separate license notice appears.
- Our account saffah_codex_6s_agg2, https://duck.ac/submission/106134: retained 256-word upper scratch spacing.
Approach:
- Use A-panel row stride KC+16, updating all handwritten microkernel broadcast offsets.
Purpose:
- Test L1 cache-set placement of the 64-square leaf operand panel in mmml2k.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/106080: retained the accepted six-level Winograd engine and all inherited citations; no separate license notice appears.
- Our account saffah_codex_6s_agg2, https://duck.ac/submission/106055: inherited sixth-level Winograd idea is credited in parent source.
Approach:
- Space the three upper Winograd scratch-plane groups by 256 long long elements and shift the nested leaf workspace past the expanded middle group.
Purpose:
- Test cache-set placement of temporary products in the official mmml2k environment.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #105976 <https://duck.ac/submission/105976>
// （判题机 820.180 ms）—— 本文件正文段（自 `#include` 起，本发新增的 `winograd_mini` 及其
// 缓冲声明除外）与该提交**逐字节相同**。它件头自述其来源：本账号提交 #105789（见 [2]）。
// 它自列的数学来源：Brice Boyer et al.,
// "Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm",
// <https://hal.science/hal-00163141v3/file/winoschedule.pdf>
// —— 只用七乘积方程与输出重组的数学思想，未复制文本或代码；该文无许可证声明。
// [2] 本账号（saffah_cc_v41_agg1）提交 #105789 <https://duck.ac/submission/105789>（831.898902 ms）
// —— 本文件的引擎与顶层、以及四把刀（最内层物化 Winograd 乘积 256 -> 7x128、`wg_combine` 的
// kickstart 预取、`material_*` 的 2 MB 对齐 + THP、Winograd 平面间隔 32 字）全部来自它。
// 本账号 #106003 <https://duck.ac/submission/106003> 采纳了 [1] 的整份；本账号 #106015
// <https://duck.ac/submission/106015>（819.987544 ms）= 本发的**基座**（单变量 `KC 512 -> 128`）。
// [3] ★★ 本发新增的唯一改动（第五级物化）来自：duck.ac 用户 saffah_codex_6s_agg2 在题 `mmml4k`
// 上的提交 #106055 <https://duck.ac/submission/106055>（5526.277 ms，自其 #106009 的 5718.922
// = **−3.37%**）。该件自述 "Materialize a sixth Winograd level on each 128-square product,
// with 64-square leaf multiplications and an independent scratch buffer"。本文件只取用其
// **思想与结构形态**：(a) 在最内层物化级之下再加一级 Winograd（128 乘积 -> 7 个 64 乘积，新叶 64）；
// (b) 新级的 5 个 scratch 平面放在**独立的 2 MB 对齐数组**（`wg_buf_mini[5*64*64 + 4*256]`）里，
// 不与既有平面缓冲共享 4 KB 相位；(c) 该级平面相互间隔 256 个 long long（2048 B）。
// `winograd_mini` 的代码是我们按本引擎（= 他件的祖先）的对应位置自行重写的：函数体与
// `winograd_prod` 同构（h 减半、平面换缓冲），7 个乘积用 `mmxR(h=64, ..., g_scratch+(1<<20))`。
// [4] 本文件未取用其它他人代码。
// ======================
// ===== 思路 =====
// 本发为**单变量**：在最内层物化级（`winograd_prod`：256 乘积 -> 7x128 叶）之下**再加一级物化
// Winograd**（128 乘积 -> 7x64 叶，**新叶 = 64**）。这是协调者指派的"物化更多级 / 新叶尺寸"轴在
// 本题上的**逐引擎定价**。
//
// 为什么值得量（当前 mine = 819.987544（#106015），T = 820.180234，严支 thr = 811.979432 ⇒ 缺口
// −8.008 ms = −0.98%）：
// * 本引擎此刻**全部剩余乘积工作都在最内层**：2048^3 x (7/8)^4 = 58.6% 的 MAC 里，128 乘积
// （2401 个 = 7^4）就是整棵树的叶子。再物化一级 ⇒ 最内层 MAC x (7/8) ⇒ **全题 MAC −12.5%**
// （按上一级实测的 mikro 单价折算约 −85 ms 量级），代价是新增 2401 x 8 x 64^2 = 78.6 M 元素次
// 变换（新级平面 = 64^2 x 8 = **32 KB/个**，可整块驻留 L2 ⇒ 单位成本应低于上一级的 128^2 平面）。
// * ★★ 家族律"新叶 = 64 ⇒ 判负"来自**别的引擎**（`sm_` 的四点归纳），已被对手在 `mmml4k` 上的
// **判题机实测 −3.37%**（[3]）推翻 ⇒ 口径为**逐引擎定价、禁止外推**。
// ⇒ 闸门①：本题这份引擎 `MC = 6`、叶级账单（MR=6 / MC=6 / KC=128 面板 + 一行 mikro 块）与 [3]
// 所在引擎同形，且**变换完全相同**（都是"叶 128 -> 7x64"）⇒ 可量。
// * 不引入新的"大"常驻缓冲：新级 5 个平面共 5x64^2 + 4x256 = 21504 字 = **168 KB**（远小于上级的
// 656 KB），且声明在所有既有数组**之后** ⇒ **不移动任何既有数组的地址**（不引入布局抽签）。
//
// 正确性（本地闸门，基座 = 我方 #106015）：`mml2_chk2.py` n = 2048 / 1024 / 512 与基座**逐位相同**
// （n=2048 FNV=4180d8ca30183001）。除新增 `wg_buf_mini`/`winograd_mini` 及把 `winograd_prod` 的
// `PR` 宏改为调用 `winograd_mini` 之外，正文一字未改。
// ================
#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 v; __asm__("vmovdqu %1,%0" : "=x"(v) : "m"(*(const __m256i *)(p))); v; })
#define STQ(p,v) __asm__ volatile("vmovdqu %1,%0" : "=m"(*(__m256i*)(p)) : "x"(v))
#include <stddef.h>
#include <stdint.h>
#include <sys/mman.h>
#ifndef MADV_HUGEPAGE
#define MADV_HUGEPAGE 14
#endif
static inline void wg_store_nt(long long *p, __m256i v) {
__asm__ volatile("vmovntdq %0, %1" ::"x"(v), "m"(*(__m256i *)p) : "memory");
}
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC push_options
#pragma GCC target("avx2")
typedef long long TE;
#define MR 6
#define KC 128
#define MC 6
#define NC 1024
#define UU 1
#define CUTL 256
/* 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 + 16)
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) {
register __m256i v0 asm("ymm0");
register __m256i v1 asm("ymm1");
register __m256i v2 asm("ymm2");
register __m256i v3 asm("ymm3");
register __m256i v4 asm("ymm4");
register __m256i v5 asm("ymm5");
register __m256i v6 asm("ymm6");
register __m256i v7 asm("ymm7");
register __m256i v8 asm("ymm8");
register __m256i v9 asm("ymm9");
register __m256i v10 asm("ymm10");
register __m256i v11 asm("ymm11");
if (rows == 4) {
__asm__ volatile(
"test %[cnt], %[cnt]\n\t"
"jle 2f\n\t"
"vmovdqu 0(%[bp]), %%ymm12\n\t"
"vmovdqu 4(%[bp]), %%ymm13\n\t"
"vbroadcastsd 0(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm0\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm1\n\t"
"vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm2\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm3\n\t"
"vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm4\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm5\n\t"
"vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm6\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm7\n\t"
"add $8, %[ap]\n\t"
"add $32, %[bp]\n\t"
"dec %[cnt]\n\t"
"jle 2f\n\t"
"1:\n\t"
"vmovdqu 0(%[bp]), %%ymm12\n\t"
"vmovdqu 4(%[bp]), %%ymm13\n\t"
"vbroadcastsd 0(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
"vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
"vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
"vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
"add $8, %[ap]\n\t"
"add $32, %[bp]\n\t"
"dec %[cnt]\n\t"
"jnz 1b\n\t"
"2:\n\t"
""
: [ap] "+r"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt), [v0] "=x"(v0), [v1] "=x"(v1), [v2] "=x"(v2), [v3] "=x"(v3), [v4] "=x"(v4), [v5] "=x"(v5), [v6] "=x"(v6), [v7] "=x"(v7)
:
: "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory");
} else
if (rows == 2) {
__asm__ volatile(
"test %[cnt], %[cnt]\n\t"
"jle 2f\n\t"
"vmovdqu 0(%[bp]), %%ymm12\n\t"
"vmovdqu 4(%[bp]), %%ymm13\n\t"
"vbroadcastsd 0(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm0\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm1\n\t"
"vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm2\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm3\n\t"
"add $8, %[ap]\n\t"
"add $32, %[bp]\n\t"
"dec %[cnt]\n\t"
"jle 2f\n\t"
"1:\n\t"
"vmovdqu 0(%[bp]), %%ymm12\n\t"
"vmovdqu 4(%[bp]), %%ymm13\n\t"
"vbroadcastsd 0(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
"vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
"add $8, %[ap]\n\t"
"add $32, %[bp]\n\t"
"dec %[cnt]\n\t"
"jnz 1b\n\t"
"2:\n\t"
""
: [ap] "+r"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt), [v0] "=x"(v0), [v1] "=x"(v1), [v2] "=x"(v2), [v3] "=x"(v3)
:
: "ymm4", "ymm5", "ymm6", "ymm7",
"ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory");
} else
__asm__ volatile(
"test %[cnt], %[cnt]\n\t"
"jle 2f\n\t"
"vmovdqu 0(%[bp]), %%ymm12\n\t"
"vmovdqu 4(%[bp]), %%ymm13\n\t"
"vbroadcastsd 0(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm0\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm1\n\t"
"vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm2\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm3\n\t"
"vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm4\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm5\n\t"
"vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm6\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm7\n\t"
"vbroadcastsd 4608(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm8\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm9\n\t"
"vbroadcastsd 5760(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm10\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm11\n\t"
"prefetcht0 1536(%[bp])\n\t"
"vmovdqu 32(%[bp]), %%ymm12\n\t"
"add $8, %[ap]\n\t"
"add $32, %[bp]\n\t"
"dec %[cnt]\n\t"
"jle 2f\n\t"
"1:\n\t"
"vmovdqu 4(%[bp]), %%ymm13\n\t"
"vbroadcastsd 0(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
"vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
"vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
"vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
"vbroadcastsd 4608(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm9, %%ymm9\n\t"
"vbroadcastsd 5760(%[ap]), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
"prefetcht0 1536(%[bp])\n\t"
"vmovdqu 32(%[bp]), %%ymm12\n\t"
"add $8, %[ap]\n\t"
"add $32, %[bp]\n\t"
"dec %[cnt]\n\t"
"jnz 1b\n\t"
"2:\n\t"
""
: [ap] "+r"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt), [v0] "=x"(v0), [v1] "=x"(v1), [v2] "=x"(v2), [v3] "=x"(v3), [v4] "=x"(v4), [v5] "=x"(v5), [v6] "=x"(v6), [v7] "=x"(v7), [v8] "=x"(v8), [v9] "=x"(v9), [v10] "=x"(v10), [v11] "=x"(v11)
:
: "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory");
if (rows > 0) {
__m256i lo = v0, hi = v1;
__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 = v2, hi = v3;
__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 = v4, hi = v5;
__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 = v6, hi = v7;
__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 = v8, hi = v9;
__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 = v10, hi = v11;
__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)); } }
}
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[4], c[4], s[4];
};
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;
/* [ml2a7_] The (nt, s[1]-sign) dispatch and the per-row source base
`P + (r[t]+ic+r)*ld + c[t]+pc` used to be re-evaluated inside the r loop, i.e. MC
times per call; the judge's gcc9 spends ~20 instructions per row on that prologue.
Hoisting it gives one dispatch per call and makes the row addressing a pointer walk.
Every k-loop body below is byte-identical to the original. */
const long long *b0 = P + (size_t)(o->r[0] + ic) * ld + o->c[0] + pc;
const long long *b1 = (nt == 2) ? P + (size_t)(o->r[1] + ic) * ld + o->c[1] + pc : 0;
const int sg1 = (nt == 2) ? o->s[1] : 1;
long long *ap = Apanel;
if (nt == 2 && sg1 > 0) {
for (const long long *s0 = b0, *s1 = b1; ap < Apanel + (size_t)(mc < MC ? mc : MC) * KCP; s0 += ld, s1 += ld, ap += KCP) {
long long *p = ap;
int k = 0;
for (; k + 16 <= kc; k += 16) {
#define MLA4(J) do { \
__m256i a_ = LDQ((const __m256i *)(s0 + k + (J))); \
__m256i b_ = LDQ((const __m256i *)(s1 + k + (J))); \
STQ((__m256i *)(p + (J)), _mm256_add_epi64(a_, b_)); } while (0)
MLA4(0); MLA4(4); MLA4(8); MLA4(12);
#undef MLA4
p += 16;
}
if (k < kc) {
for (; k + 4 <= kc; k += 4) {
__m256i v = LDQ((const __m256i *)(s0 + k));
__m256i w = LDQ((const __m256i *)(s1 + k));
STQ((__m256i *)p, _mm256_add_epi64(v, w));
p += 4;
}
for (; k < kc; k++) *p++ = s0[k] + s1[k];
}
}
} else if (nt == 2) {
for (const long long *s0 = b0, *s1 = b1; ap < Apanel + (size_t)(mc < MC ? mc : MC) * KCP; s0 += ld, s1 += ld, ap += KCP) {
long long *p = ap;
int k = 0;
for (; k + 16 <= kc; k += 16) {
#define MLS4(J) do { \
__m256i a_ = LDQ((const __m256i *)(s0 + k + (J))); \
__m256i b_ = LDQ((const __m256i *)(s1 + k + (J))); \
STQ((__m256i *)(p + (J)), _mm256_sub_epi64(a_, b_)); } while (0)
MLS4(0); MLS4(4); MLS4(8); MLS4(12);
#undef MLS4
p += 16;
}
if (k < kc) {
for (; k + 4 <= kc; k += 4) {
__m256i v = LDQ((const __m256i *)(s0 + k));
__m256i w = LDQ((const __m256i *)(s1 + k));
STQ((__m256i *)p, _mm256_sub_epi64(v, w));
p += 4;
}
for (; k < kc; k++) *p++ = s0[k] - s1[k];
}
}
} else if (nt == 1) {
for (const long long *s0 = b0; ap < Apanel + (size_t)(mc < MC ? mc : MC) * KCP; s0 += ld, ap += KCP) {
long long *p = ap;
int k = 0;
for (; k + 16 <= kc; k += 16) {
#define MLC4(J) STQ((__m256i *)(p + (J)), LDQ((const __m256i *)(s0 + k + (J))))
MLC4(0); MLC4(4); MLC4(8); MLC4(12);
#undef MLC4
p += 16;
}
if (k < kc) {
for (; k + 4 <= kc; k += 4) { STQ((__m256i *)p, LDQ((const __m256i *)(s0 + k))); p += 4; }
for (; k < kc; k++) *p++ = s0[k];
}
}
} else {
for (int r = 0; r < MC && r < mc; r++, ap += KCP) {
long long *p = ap;
int k = 0;
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);
}
}
}
}
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;
/* ---- 8-column pass: the two 4-column groups jb and jb+4 share ONE 64-byte
source line (the row stride is ld*8 bytes, so a 4-column load uses exactly
half a line). The 4-column loop walks all kc rows before returning to the
neighbouring group, i.e. after kc strided accesses, so that line is long
evicted and gets fetched a second time. Taking both groups in one pass
reads a full line per row per visit: source line fetches are halved.
Panel layout is byte-identical (group jb sits at (jb/4)*KC). ---- */
for (; jb + 16 <= nc; jb += 16) {
long long *q0 = Bpanel + (size_t)(jb / 4) * kc * 4;
long long *q1 = q0 + (size_t)kc * 4;
long long *q2 = q1 + (size_t)kc * 4;
long long *q3 = q2 + (size_t)kc * 4;
if (nt <= 2) {
const long long *s0 = P + (size_t)o->r[0] * ld + o->c[0] + jc + jb;
const long long *s1 = (nt == 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc + jb : 0;
const int sg1 = (nt == 2) ? o->s[1] : 1;
{
const long long *rp0 = s0 + (size_t)pc * ld;
const long long *rp1 = (nt == 2) ? s1 + (size_t)pc * ld : 0;
long long *p0 = q0, *p1 = q1, *p2 = q2, *p3 = q3;
int k = 0;
#define PBR(KK, RO, RO1) do { \
__m256i v0 = LDQ((const __m256i *)(RO)); \
__m256i v1 = LDQ((const __m256i *)(RO + 4)); \
__m256i v2 = LDQ((const __m256i *)(RO + 8)); \
__m256i v3 = LDQ((const __m256i *)(RO + 12)); \
if (nt == 2) { \
__m256i w0 = LDQ((const __m256i *)(RO1)); \
__m256i w1 = LDQ((const __m256i *)(RO1 + 4)); \
__m256i w2 = LDQ((const __m256i *)(RO1 + 8)); \
__m256i w3 = LDQ((const __m256i *)(RO1 + 12)); \
if (sg1 > 0) { v0 = _mm256_add_epi64(v0, w0); v1 = _mm256_add_epi64(v1, w1); \
v2 = _mm256_add_epi64(v2, w2); v3 = _mm256_add_epi64(v3, w3); } \
else { v0 = _mm256_sub_epi64(v0, w0); v1 = _mm256_sub_epi64(v1, w1); \
v2 = _mm256_sub_epi64(v2, w2); v3 = _mm256_sub_epi64(v3, w3); } \
} \
_mm256_store_si256((__m256i *)(p0), _mm256_shuffle_epi32(v0, 0xB1)); \
_mm256_store_si256((__m256i *)(p1), _mm256_shuffle_epi32(v1, 0xB1)); \
_mm256_store_si256((__m256i *)(p2), _mm256_shuffle_epi32(v2, 0xB1)); \
_mm256_store_si256((__m256i *)(p3), _mm256_shuffle_epi32(v3, 0xB1)); \
p0 += 4; p1 += 4; p2 += 4; p3 += 4; (void)(KK); } while (0)
for (; k + 1 < kc; k += 2) {
PBR(k, rp0, rp1);
rp0 += ld; rp1 += ld;
PBR(k + 1, rp0, rp1);
rp0 += ld; rp1 += ld;
}
for (; k < kc; k++) {
PBR(k, rp0, rp1);
rp0 += ld; rp1 += ld;
}
#undef PBR
}
} else {
for (int k = 0; k < kc; k++) {
__m256i v0 = _mm256_setzero_si256(), v1 = _mm256_setzero_si256();
__m256i v2 = _mm256_setzero_si256(), v3 = _mm256_setzero_si256();
for (int tt = 0; tt < nt; tt++) {
const long long *r0 = P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb;
__m256i w0 = LDQ((const __m256i *)r0);
__m256i w1 = LDQ((const __m256i *)(r0 + 4));
__m256i w2 = LDQ((const __m256i *)(r0 + 8));
__m256i w3 = LDQ((const __m256i *)(r0 + 12));
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); }
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); }
}
_mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), _mm256_shuffle_epi32(v0, 0xB1));
_mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), _mm256_shuffle_epi32(v1, 0xB1));
_mm256_store_si256((__m256i *)(q2 + (size_t)k * 4), _mm256_shuffle_epi32(v2, 0xB1));
_mm256_store_si256((__m256i *)(q3 + (size_t)k * 4), _mm256_shuffle_epi32(v3, 0xB1));
}
}
}
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), _mm256_shuffle_epi32(v, 0xB1));
}
} else if (full) {
for (int k = 0; 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), _mm256_shuffle_epi32(v, 0xB1));
}
} 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), _mm256_shuffle_epi32(v, 0xB1));
}
}
}
}
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;
}
template<int Sign>
static __attribute__((always_inline)) inline void packA_fixed64(const long long*s0,const long long*s1,int ld,int mc){
long long*d=Apanel;
for(int r=0;r<mc;r++){
#pragma GCC unroll 4
for(int k=0;k<64;k+=4){
__m256i v=LDQ((const __m256i*)(s0+k));
if constexpr(Sign) { __m256i w=LDQ((const __m256i*)(s1+k));v=Sign>0?_mm256_add_epi64(v,w):_mm256_sub_epi64(v,w); }
STQ((__m256i*)(d+k),v);
}
s0+=ld;if constexpr(Sign)s1+=ld;d+=KCP;
}
}
template<int Sign>
static __attribute__((noinline)) void packB_fixed64(const long long*s0,const long long*s1,int ld){
for(int jb=0;jb<64;jb+=16){
long long*q0=Bpanel+(jb/4)*64*4,*q1=q0+256,*q2=q1+256,*q3=q2+256;
const long long*r0=s0+jb,*r1=s1;
if constexpr(Sign)r1+=jb;
#pragma GCC unroll 1
for(int k=0;k<64;k++){
__m256i v0=LDQ((const __m256i*)r0),v1=LDQ((const __m256i*)(r0+4)),v2=LDQ((const __m256i*)(r0+8)),v3=LDQ((const __m256i*)(r0+12));
if constexpr(Sign){
__m256i w0=LDQ((const __m256i*)r1),w1=LDQ((const __m256i*)(r1+4)),w2=LDQ((const __m256i*)(r1+8)),w3=LDQ((const __m256i*)(r1+12));
if constexpr(Sign>0){v0=_mm256_add_epi64(v0,w0);v1=_mm256_add_epi64(v1,w1);v2=_mm256_add_epi64(v2,w2);v3=_mm256_add_epi64(v3,w3);}
else {v0=_mm256_sub_epi64(v0,w0);v1=_mm256_sub_epi64(v1,w1);v2=_mm256_sub_epi64(v2,w2);v3=_mm256_sub_epi64(v3,w3);}
}
_mm256_store_si256((__m256i*)q0,_mm256_shuffle_epi32(v0,0xB1));
_mm256_store_si256((__m256i*)q1,_mm256_shuffle_epi32(v1,0xB1));
_mm256_store_si256((__m256i*)q2,_mm256_shuffle_epi32(v2,0xB1));
_mm256_store_si256((__m256i*)q3,_mm256_shuffle_epi32(v3,0xB1));
q0+=4;q1+=4;q2+=4;q3+=4;r0+=ld;if constexpr(Sign)r1+=ld;
}
}
}
template<int ASign>
static __attribute__((noinline)) void leaf64_by_sign(const Opd*ao,const Opd*bo,long long*c,int ldc){
const long long *a0=ao->p+(size_t)ao->r[0]*ao->ld+ao->c[0];
const long long *a1=ao->nt==2?ao->p+(size_t)ao->r[1]*ao->ld+ao->c[1]:a0;
const long long *b0=bo->p+(size_t)bo->r[0]*bo->ld+bo->c[0];
const long long *b1=bo->nt==2?bo->p+(size_t)bo->r[1]*bo->ld+bo->c[1]:b0;
if(bo->nt==1)packB_fixed64<0>(b0,b1,bo->ld);
else if(bo->s[1]>0)packB_fixed64<1>(b0,b1,bo->ld);else packB_fixed64<-1>(b0,b1,bo->ld);
for(int ic=0;ic<60;ic+=6){
packA_fixed64<ASign>(a0+(size_t)ic*ao->ld,a1+(size_t)ic*ao->ld,ao->ld,6);
for(int j=0;j<64;j+=4) mikro(64,Apanel,Bpanel+(size_t)j*64,c+(size_t)ic*ldc+j,ldc,0,nullptr,0,0,6);
}
packA_fixed64<ASign>(a0+(size_t)60*ao->ld,a1+(size_t)60*ao->ld,ao->ld,4);
for(int j=0;j<64;j+=4) mikro(64,Apanel,Bpanel+(size_t)j*64,c+(size_t)60*ldc+j,ldc,0,nullptr,0,0,4);
}
static void fixed_leaf64_product(const Opd*ao,const Opd*bo,long long*c,int ldc){
if(ao->nt==1)leaf64_by_sign<0>(ao,bo,c,ldc);
else if(ao->s[1]>0)leaf64_by_sign<1>(ao,bo,c,ldc);
else leaf64_by_sign<-1>(ao,bo,c,ldc);
}
/* 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. */
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;
if (op == 0) {
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 + 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]; }
}
}
}
/* Reuse one pair of temporary 1024-square operands across depth-first leaves. */
static long long material_A[1024u * 1024u] __attribute__((aligned(1 << 21)));
__attribute__((constructor)) static void l1n_huge_material_A(void) {
long long *a = material_A; uintptr_t p = (uintptr_t)a;
uintptr_t e = p + sizeof(material_A);
uintptr_t ha = (p + ((uintptr_t)1 << 21) - 1) & ~(((uintptr_t)1 << 21) - 1);
uintptr_t he = e & ~(((uintptr_t)1 << 21) - 1);
if (he > ha) madvise((void *)ha, (size_t)(he - ha), MADV_HUGEPAGE);
}
static long long material_B[1024u * 1024u] __attribute__((aligned(1 << 21)));
__attribute__((constructor)) static void l1n_huge_material_B(void) {
long long *a = material_B; uintptr_t p = (uintptr_t)a;
uintptr_t e = p + sizeof(material_B);
uintptr_t ha = (p + ((uintptr_t)1 << 21) - 1) & ~(((uintptr_t)1 << 21) - 1);
uintptr_t he = e & ~(((uintptr_t)1 << 21) - 1);
if (he > ha) madvise((void *)ha, (size_t)(he - ha), MADV_HUGEPAGE);
}
/* 512 级物化的目的地：必须与 material_A/B 分开 —— 1024 级钩子的递归仍在读 material_A/B。 */
static long long material_C[512u * 512u] __attribute__((aligned(1 << 21)));
__attribute__((constructor)) static void l1n_huge_material_C(void) {
long long *a = material_C; uintptr_t p = (uintptr_t)a;
uintptr_t e = p + sizeof(material_C);
uintptr_t ha = (p + ((uintptr_t)1 << 21) - 1) & ~(((uintptr_t)1 << 21) - 1);
uintptr_t he = e & ~(((uintptr_t)1 << 21) - 1);
if (he > ha) madvise((void *)ha, (size_t)(he - ha), MADV_HUGEPAGE);
}
static long long material_D[512u * 512u] __attribute__((aligned(1 << 21)));
__attribute__((constructor)) static void l1n_huge_material_D(void) {
long long *a = material_D; uintptr_t p = (uintptr_t)a;
uintptr_t e = p + sizeof(material_D);
uintptr_t ha = (p + ((uintptr_t)1 << 21) - 1) & ~(((uintptr_t)1 << 21) - 1);
uintptr_t he = e & ~(((uintptr_t)1 << 21) - 1);
if (he > ha) madvise((void *)ha, (size_t)(he - ha), MADV_HUGEPAGE);
}
static void materializeN(const Opd *o, long long *dst, const int n) {
for (int i = 0; i < n; i++) {
long long *out = dst + (size_t)i * n;
int j = 0;
for (; j + 4 <= n; j += 4) {
__m256i value = _mm256_setzero_si256();
for (int t = 0; t < o->nt; t++) {
const long long *src = o->p + (size_t)(o->r[t] + i) * o->ld + o->c[t] + j;
__m256i part = LDQ((const __m256i *)src);
value = o->s[t] < 0 ? _mm256_sub_epi64(value, part) : _mm256_add_epi64(value, part);
}
STQ((__m256i *)(out + j), value);
}
for (; j < n; j++) {
unsigned long long value = 0;
for (int t = 0; t < o->nt; t++) {
unsigned long long part = (unsigned long long)o->p[(size_t)(o->r[t] + i) * o->ld + o->c[t] + j];
value = o->s[t] < 0 ? value - part : value + part;
}
out[j] = (long long)value;
}
}
}
/* single-destination recursive Strassen: d op= Ao*Bo, with T a scratch block of h*h */
static void mmxR(int n, const Opd *Ao, const Opd *Bo, long long *d, int ld, int op,
long long *T) {
if (n == 1024 && (Ao->nt > 1 || Bo->nt > 1)) {
Opd am, bm; const Opd *a = Ao, *b = Bo;
if (Ao->nt > 1) { materializeN(Ao, material_A, 1024); am = o1(material_A, n); a = &am; }
if (Bo->nt > 1) { materializeN(Bo, material_B, 1024); bm = o1(material_B, n); b = &bm; }
mmxR(n, a, b, d, ld, op, T);
return;
}
if (n <= CUTL || (n & 1) || Ao->nt > 2 || Bo->nt > 2) {
leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); return; }
/* n == 512 且本次**本来就会递归**（守卫已过 ⇒ nt <= 2、n > CUTL）时，先把惰性操作数
物化成 512^2 连续矩阵再递归。语义恒等（mod 2^64 下加法可结合），但 h=256 那一层的描述符
nt 从 {2,4} 掉到 {2} ⇒ **每个 256 叶的 pack 源读项数减半**（pack 流量 ≈ −2 MB/叶）。 */
if (n == 512 && (Ao->nt > 1 || Bo->nt > 1)) {
Opd am, bm; const Opd *a = Ao, *b = Bo;
if (Ao->nt > 1) { materializeN(Ao, material_C, 512); am = o1(material_C, 512); a = &am; }
if (Bo->nt > 1) { materializeN(Bo, material_D, 512); bm = o1(material_D, 512); b = &bm; }
mmxR(n, a, b, d, ld, op, T);
return;
}
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); \
cmb(Q[qi], ld, T, h, h, (si), oa); \
cmb(Q[qj], ld, T, h, h, (sj), ob); \
} \
} 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)
{ /* 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
}
static inline __m256i wg_load(const long long *p) { return LDQ((const __m256i *)p); }
static inline void wg_store(long long *p, __m256i v) { STQ((__m256i *)p, v); }
static long long wg_buf[5 * 2048 * 2048] __attribute__((aligned(64)));
/* ---------------- [ml2a2_/v5] fused two-pass chain, cold writes non-temporal -- */
static void wg_nAS(int h, long long *S1, long long *S3, const long long *A11,
const long long *A21, const long long *A22, int lda) {
for (int i = 0; i < h; i++) {
const long long *c = A11 + (size_t)i*lda, *a = A21 + (size_t)i*lda, *b = A22 + (size_t)i*lda;
long long *s1 = S1 + (size_t)i*h, *s3 = S3 + (size_t)i*h;
for (int j = 0; j < h; j += 4) {
__m256i va = wg_load(a+j), vb = wg_load(b+j), vc = wg_load(c+j);
wg_store_nt(s1+j, _mm256_add_epi64(va, vb)); /* S1 = A21+A22 */
wg_store_nt(s3+j, _mm256_sub_epi64(vc, va)); /* S3 = A11-A21 */
}
}
}
static void wg_nBT(int h, long long *T1, long long *T3, const long long *B11,
const long long *B12, const long long *B22, int ldb) {
for (int i = 0; i < h; i++) {
const long long *c = B11 + (size_t)i*ldb, *a = B12 + (size_t)i*ldb, *b = B22 + (size_t)i*ldb;
long long *t1 = T1 + (size_t)i*h, *t3 = T3 + (size_t)i*h;
for (int j = 0; j < h; j += 4) {
__m256i va = wg_load(a+j), vb = wg_load(b+j), vc = wg_load(c+j);
wg_store_nt(t1+j, _mm256_sub_epi64(va, vc)); /* T1 = B12-B11 */
wg_store_nt(t3+j, _mm256_sub_epi64(vb, va)); /* T3 = B22-B12 */
}
}
}
static void wg_nAS2(int h, long long *S2, long long *S4, const long long *S1,
const long long *A11, const long long *A12, int lda) {
for (int i = 0; i < h; i++) {
const long long *p = S1 + (size_t)i*h, *c = A11 + (size_t)i*lda, *d = A12 + (size_t)i*lda;
long long *s2 = S2 + (size_t)i*h, *s4 = S4 + (size_t)i*h;
for (int j = 0; j < h; j += 4) {
__m256i v = _mm256_sub_epi64(wg_load(p+j), wg_load(c+j)); /* S2 = S1-A11 */
wg_store_nt(s2+j, v);
wg_store(s4+j, _mm256_sub_epi64(wg_load(d+j), v)); /* S4 = A12-S2 (in place) */
}
}
}
static void wg_nBT2(int h, long long *T2, long long *T4, const long long *T1,
const long long *B22, const long long *B21, int ldb) {
for (int i = 0; i < h; i++) {
const long long *p = T1 + (size_t)i*h, *c = B22 + (size_t)i*ldb, *d = B21 + (size_t)i*ldb;
long long *t2 = T2 + (size_t)i*h, *t4 = T4 + (size_t)i*h;
for (int j = 0; j < h; j += 4) {
__m256i v = _mm256_sub_epi64(wg_load(c+j), wg_load(p+j)); /* T2 = B22-T1 */
wg_store_nt(t2+j, v);
wg_store(t4+j, _mm256_sub_epi64(v, wg_load(d+j))); /* T4 = T2-B21 (in place) */
}
}
}
static void wg_add(int h, long long *dst, const long long *a, int lda,
const long long *b, int ldb, int sg) {
for (int i = 0; i < h; i++) {
long long *d = dst + (size_t)i*h;
const long long *x = a + (size_t)i*lda, *y = b + (size_t)i*ldb;
for (int j = 0; j < h; j += 4) {
__m256i vx=wg_load(x+j), vy=wg_load(y+j);
wg_store(d+j, sg > 0 ? _mm256_add_epi64(vx,vy) : _mm256_sub_epi64(vx,vy));
}
}
}
static void wg_mutate(int h, long long *dst, const long long *b, int ldb, int reverse) {
for (int i = 0; i < h; i++) {
long long *d = dst + (size_t)i*h; const long long *y=b+(size_t)i*ldb;
for (int j=0; j<h; j+=4) {
__m256i vd=wg_load(d+j),vy=wg_load(y+j);
wg_store(d+j, reverse ? _mm256_sub_epi64(vy,vd) : _mm256_sub_epi64(vd,vy));
}
}
}
static void wg_combine(int h, long long *C, int ldc, const long long *M1,
const long long *M2, const long long *M3, const long long *M4,
const long long *M5, const long long *M6, const long long *M7) {
for (int i=0; i<h; i++) {
const long long *a=M1+(size_t)i*ldc,*b=M2+(size_t)i*ldc,*c=M3+(size_t)i*ldc;
const long long *d=M4+(size_t)i*ldc,*e=M5+(size_t)i*h,*f=M6+(size_t)i*h,*g=M7+(size_t)i*h;
long long *c11=C+(size_t)i*ldc,*c12=c11+h,*c21=c11+(size_t)h*ldc,*c22=c21+h;
if (i + 4 < h) { /* 每行一脚：行间跳越超出 L2 流预取器能力（M 只读一次 => nta 不污染） */
__builtin_prefetch((const void *)(a+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(b+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(c+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(d+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(e+(size_t)4*h),0,0);
__builtin_prefetch((const void *)(f+(size_t)4*h),0,0);
__builtin_prefetch((const void *)(g+(size_t)4*h),0,0);
}
for (int j=0; j<h; j+=4) {
__m256i p1=wg_load(a+j),p2=wg_load(b+j),p3=wg_load(c+j),p4=wg_load(d+j);
__m256i p5=wg_load(e+j),p6=wg_load(f+j),p7=wg_load(g+j);
__m256i u2=_mm256_add_epi64(p1,p6),u3=_mm256_add_epi64(u2,p7);
__m256i u4=_mm256_add_epi64(u2,p5);
wg_store_nt(c11+j,_mm256_add_epi64(p1,p2));
wg_store_nt(c12+j,_mm256_add_epi64(u4,p3));
wg_store_nt(c21+j,_mm256_sub_epi64(u3,p4));
wg_store_nt(c22+j,_mm256_add_epi64(u3,p5));
}
}
_mm_sfence();
}
static void wg_combine_cache(int h, long long *C, int ldc, const long long *M1,
const long long *M2, const long long *M3, const long long *M4,
const long long *M5, const long long *M6, const long long *M7) {
for (int i=0; i<h; i++) {
const long long *a=M1+(size_t)i*ldc,*b=M2+(size_t)i*ldc,*c=M3+(size_t)i*ldc;
const long long *d=M4+(size_t)i*ldc,*e=M5+(size_t)i*h,*f=M6+(size_t)i*h,*g=M7+(size_t)i*h;
long long *c11=C+(size_t)i*ldc,*c12=c11+h,*c21=c11+(size_t)h*ldc,*c22=c21+h;
if (i + 4 < h) { /* 每行一脚：行间跳越超出 L2 流预取器能力（M 只读一次 => nta 不污染） */
__builtin_prefetch((const void *)(a+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(b+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(c+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(d+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(e+(size_t)4*h),0,0);
__builtin_prefetch((const void *)(f+(size_t)4*h),0,0);
__builtin_prefetch((const void *)(g+(size_t)4*h),0,0);
}
for (int j=0; j<h; j+=4) {
__m256i p1=wg_load(a+j),p2=wg_load(b+j),p3=wg_load(c+j),p4=wg_load(d+j);
__m256i p5=wg_load(e+j),p6=wg_load(f+j),p7=wg_load(g+j);
__m256i u2=_mm256_add_epi64(p1,p6),u3=_mm256_add_epi64(u2,p7);
__m256i u4=_mm256_add_epi64(u2,p5);
wg_store(c11+j,_mm256_add_epi64(p1,p2));
wg_store(c12+j,_mm256_add_epi64(u4,p3));
wg_store(c21+j,_mm256_sub_epi64(u3,p4));
wg_store(c22+j,_mm256_add_epi64(u3,p5));
}
}
}
static long long wg_buf2[5 * 1024 * 1024] __attribute__((aligned(64)));
static long long wg_buf1[5 * 512 * 512] __attribute__((aligned(64)));
static long long *wleaf = wg_buf1 + 5*(256*256) + 512; /* wg_buf1 slack: grandchild uses 5*256^2 */
static long long wg_buf_mini[5*64*64 + 4*256] __attribute__((aligned(2097152)));
__attribute__((constructor)) static void ml2a2v5_chk(void) {
if (wg_buf1 + 5*512*512 < wleaf + 7*((size_t)(128*128) + 64)) __builtin_trap();
}
static long long *const m5_base = wg_buf_mini;
/* 128 乘积 -> 7 个 64 乘积（新叶 64）。形态与 winograd_prod 完全同构，只是 h 减半、
平面换成 wg_buf_mini（独立 2 MB 对齐，与其它平面缓冲不共享 4 KB 相位）。 */
static void winograd_mini(int n, const long long *A, int lda, const long long *B, int ldb,
long long *C, int ldc) {
const int h=n/2; const size_t sz=(size_t)h*h;
long long *M1=C,*M2=C+h,*M3=C+(size_t)h*ldc,*M4=M3+h;
long long *M5=m5_base,*M6=M5+sz+256,*M7=M6+sz+256;
/* [l2y_] Strassen 自身七乘积：全部操作数是原瓦片的 1~2 项仿射 Opd ⇒
`packA`/`packB` 的 nt≤2 快路径；基座的 P/Q 六趟物化（含 M6/M3/M4 的 nt=3/4 操作数）整段消失。
M5/M6/M7 落址与基座逐字节相同（仍不用 P/Q 两平面，`wg_buf_mini` 声明未动）。 */
Opd A0=o1(A,lda), B0=o1(B,ldb);
Opd A11=oq(&A0,0,0,h), A12=oq(&A0,0,1,h), A21=oq(&A0,1,0,h), A22=oq(&A0,1,1,h);
Opd B11=oq(&B0,0,0,h), B12=oq(&B0,0,1,h), B21=oq(&B0,1,0,h), B22=oq(&B0,1,1,h);
Opd P1=oadd(&A11,&A22,1), Q1=oadd(&B11,&B22,1); /* M1=(A11+A22)(B11+B22) */
Opd P2=oadd(&A21,&A22,1), Q2=B11; /* M2=(A21+A22)B11 */
Opd P3=A11, Q3=oadd(&B12,&B22,-1); /* M3=A11(B12-B22) */
Opd P4=A22, Q4=oadd(&B21,&B11,-1); /* M4=A22(B21-B11) */
Opd P5=oadd(&A11,&A12,1), Q5=B22; /* M5=(A11+A12)B22 */
Opd P6=oadd(&A21,&A11,-1),Q6=oadd(&B11,&B12,1); /* M6=(A21-A11)(B11+B12) */
Opd P7=oadd(&A12,&A22,-1),Q7=oadd(&B21,&B22,1); /* M7=(A12-A22)(B21+B22) */
#define PRM(AO,BO,Z,LZ) do { Opd a_=(AO),b_=(BO); \
fixed_leaf64_product(&a_,&b_,(Z),(LZ)); } while(0)
PRM(P6,Q6,M6,h);
PRM(P2,Q2,M2,ldc);
PRM(P4,Q4,M4,ldc);
PRM(P1,Q1,M1,ldc);
PRM(P3,Q3,M3,ldc);
PRM(P5,Q5,M5,h);
PRM(P7,Q7,M7,h);
#undef PRM
for (int i=0;i<h;i++) {
const long long *m1=M1+(size_t)i*ldc,*m2=M2+(size_t)i*ldc;
const long long *m3=M3+(size_t)i*ldc,*m4=M4+(size_t)i*ldc;
const long long *m5=M5+(size_t)i*h,*m6=M6+(size_t)i*h,*m7=M7+(size_t)i*h;
long long *c11=C+(size_t)i*ldc,*c12=c11+h,*c21=c11+(size_t)h*ldc,*c22=c21+h;
if (i + 4 < h) { /* 每行一脚：行间跳越超出 L2 流预取器能力（M 只读一次） */
__builtin_prefetch((const void *)(m1+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(m2+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(m3+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(m4+(size_t)4*ldc),0,0);
__builtin_prefetch((const void *)(m5+(size_t)4*h),0,0);
__builtin_prefetch((const void *)(m6+(size_t)4*h),0,0);
__builtin_prefetch((const void *)(m7+(size_t)4*h),0,0);
}
for (int j=0;j<h;j+=4) {
__m256i p1=wg_load(m1+j),p2=wg_load(m2+j),p3=wg_load(m3+j),p4=wg_load(m4+j);
__m256i p5=wg_load(m5+j),p6=wg_load(m6+j),p7=wg_load(m7+j);
wg_store(c11+j,_mm256_add_epi64(_mm256_sub_epi64(_mm256_add_epi64(p1,p4),p5),p7));
wg_store(c12+j,_mm256_add_epi64(p3,p5));
wg_store(c21+j,_mm256_add_epi64(p2,p4));
wg_store(c22+j,_mm256_add_epi64(_mm256_add_epi64(_mm256_sub_epi64(p1,p2),p3),p6));
}
}
}
/* One more Winograd level at the innermost 256-product: the three tops are already
Winograd, but the products they call are plain `mmxR` leaves. Turning each 256-product
into 7 products of 128 with MATERIALIZED operand sums makes the products' operands nt=1
(pack reads unchanged, MAC -12.5% at that level). Judge-priced precedent: the rival's
#105699 group A on mmml4k, reproduced by our teammate at -1.24%. */
static void winograd_prod(int n, const long long *A, int lda, const long long *B, int ldb,
long long *C, int ldc) {
const int h=n/2; const size_t sz=(size_t)h*h;
long long *M1=C,*M2=C+h,*M3=C+(size_t)h*ldc,*M4=M3+h;
long long *M5=wleaf,*M6=M5+sz+64,*M7=M6+sz+64,*P=M7+sz+64,*Q=P+sz+64;
const long long *A11=A,*A12=A+h,*A21=A+(size_t)h*lda,*A22=A21+h;
const long long *B11=B,*B12=B+h,*B21=B+(size_t)h*ldb,*B22=B21+h;
#define PR(X,Y,Z,LX,LY,LZ) do { winograd_mini(h,(X),(LX),(Y),(LY),(Z),(LZ)); } while(0)
PR(A11,B11,M1,lda,ldb,ldc);
PR(A12,B21,M2,lda,ldb,ldc);
wg_add(h,P,A21,lda,A22,lda,+1);
wg_add(h,Q,B12,ldb,B11,ldb,-1);
PR(P,Q,M5,h,h,h);
wg_mutate(h,P,A11,lda,0);
wg_mutate(h,Q,B22,ldb,1);
PR(P,Q,M6,h,h,h);
wg_mutate(h,P,A12,lda,1);
PR(P,B22,M3,h,ldb,ldc);
wg_mutate(h,Q,B21,ldb,0);
PR(A22,Q,M4,lda,h,ldc);
wg_add(h,P,A11,lda,A21,lda,-1);
wg_add(h,Q,B22,ldb,B12,ldb,-1);
PR(P,Q,M7,h,h,h);
#undef PR
wg_combine_cache(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
}
static void winograd_grandchild(int n, const long long *A, int lda, const long long *B, int ldb,
long long *C, int ldc) {
const int h=n/2; const size_t sz=(size_t)h*h;
long long *M1=C,*M2=C+h,*M3=C+(size_t)h*ldc,*M4=M3+h;
long long *M5=wg_buf1,*M6=M5+sz+128,*M7=M6+sz+128,*P=M7+sz+128,*Q=P+sz+128;
const long long *A11=A,*A12=A+h,*A21=A+(size_t)h*lda,*A22=A21+h;
const long long *B11=B,*B12=B+h,*B21=B+(size_t)h*ldb,*B22=B21+h;
#define PROD3(X,Y,Z,LX,LY,LZ) do { winograd_prod(h,(X),(LX),(Y),(LY),(Z),(LZ)); } while(0)
PROD3(A11,B11,M1,lda,ldb,ldc);
PROD3(A12,B21,M2,lda,ldb,ldc);
wg_add(h,P,A21,lda,A22,lda,+1);
wg_add(h,Q,B12,ldb,B11,ldb,-1);
PROD3(P,Q,M5,h,h,h);
wg_mutate(h,P,A11,lda,0);
wg_mutate(h,Q,B22,ldb,1);
PROD3(P,Q,M6,h,h,h);
wg_mutate(h,P,A12,lda,1);
PROD3(P,B22,M3,h,ldb,ldc);
wg_mutate(h,Q,B21,ldb,0);
PROD3(A22,Q,M4,lda,h,ldc);
wg_add(h,P,A11,lda,A21,lda,-1);
wg_add(h,Q,B22,ldb,B12,ldb,-1);
PROD3(P,Q,M7,h,h,h);
#undef PROD3
wg_combine_cache(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
}
static void winograd_child(int n, const long long *A, int lda, const long long *B, int ldb,
long long *C, int ldc) {
const int h=n/2; const size_t sz=(size_t)h*h;
long long *M1=C,*M2=C+h,*M3=C+(size_t)h*ldc,*M4=M3+h;
long long *M5=wg_buf2,*M6=M5+sz+256,*M7=M6+sz+256,*P=M7+sz+256,*Q=P+sz+256;
long long *S3=Q+sz+256,*T3=S3+sz+256; /* [ml2a2_/v5] two extra chain slots */
const long long *A11=A,*A12=A+h,*A21=A+(size_t)h*lda,*A22=A21+h;
const long long *B11=B,*B12=B+h,*B21=B+(size_t)h*ldb,*B22=B21+h;
#define PROD2(X,Y,Z,LX,LY,LZ) do { winograd_grandchild(h,(X),(LX),(Y),(LY),(Z),(LZ)); } while(0)
PROD2(A11,B11,M1,lda,ldb,ldc);
PROD2(A12,B21,M2,lda,ldb,ldc);
wg_nAS(h,P,S3,A11,A21,A22,lda); /* P=S1=A21+A22 ; S3=A11-A21 (NT) */
wg_nBT(h,Q,T3,B11,B12,B22,ldb); /* Q=T1=B12-B11 ; T3=B22-B12 (NT) */
PROD2(P,Q,M5,h,h,h);
PROD2(S3,T3,M7,h,h,h);
wg_nAS2(h,S3,P,P,A11,A12,lda); /* S3=S2=S1-A11(NT) ; P=S4=A12-S2 */
wg_nBT2(h,T3,Q,Q,B22,B21,ldb); /* T3=T2=B22-T1(NT) ; Q=T4=T2-B21 */
PROD2(S3,T3,M6,h,h,h);
PROD2(P,B22,M3,h,ldb,ldc);
PROD2(A22,Q,M4,lda,h,ldc);
#undef PROD2
wg_combine_cache(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
}
static void winograd_top(int n, const long long *A, const long long *B, long long *C) {
const int h=n/2; const size_t sz=(size_t)h*h;
long long *M1=C,*M2=C+h,*M3=C+(size_t)h*n,*M4=M3+h;
long long *M5=wg_buf,*M6=M5+sz+512,*M7=M6+sz+512,*P=M7+sz+512,*Q=P+sz+512;
long long *S3=Q+sz+512,*T3=S3+sz+512; /* [ml2a2_/v5] two extra chain slots */
const long long *A11=A,*A12=A+h,*A21=A+(size_t)h*n,*A22=A21+h;
const long long *B11=B,*B12=B+h,*B21=B+(size_t)h*n,*B22=B21+h;
#define PROD(X,Y,Z,LX,LY,LZ) do { winograd_child(h,(X),(LX),(Y),(LY),(Z),(LZ)); } while(0)
PROD(A11,B11,M1,n,n,n);
PROD(A12,B21,M2,n,n,n);
wg_nAS(h,P,S3,A11,A21,A22,n); /* P=S1=A21+A22 ; S3=A11-A21 (NT) */
wg_nBT(h,Q,T3,B11,B12,B22,n); /* Q=T1=B12-B11 ; T3=B22-B12 (NT) */
PROD(P,Q,M5,h,h,h);
PROD(S3,T3,M7,h,h,h);
wg_nAS2(h,S3,P,P,A11,A12,n); /* S3=S2=S1-A11(NT) ; P=S4=A12-S2 */
wg_nBT2(h,T3,Q,Q,B22,B21,n); /* T3=T2=B22-T1(NT) ; Q=T4=T2-B21 */
PROD(S3,T3,M6,h,h,h);
PROD(P,B22,M3,h,n,n);
PROD(A22,Q,M4,n,h,n);
#undef PROD
wg_combine_cache(h,C,n,M1,M2,M3,M4,M5,M6,M7);
}
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; }
if (n == 4096) winograd_top(n, A, B, C);
else if (n == 2048) winograd_top(n, A, B, C);
else mmxR(n, &Ao, &Bo, C, n, 0, g_scratch);
}
#pragma GCC pop_options
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 743.398 ms | 105 MB + 300 KB | Accepted | Score: 100 | 显示更多 |