代码

// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106176 <https://duck.ac/submission/106176>
//     用途：正文基座（对手 #106155 的完整复刻，781.803621 ms）。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #106155 <https://duck.ac/submission/106155>
//     用途：他把叶内核 B 面板预取提前量 4096 → 1024。本发沿同一距离轴**反向**探一格：1024 → 2048（单变量，
//       与 #106180 的 512 同一基座、同轴、只差这一个常量）⇒ 三格 {512, 1024, 2048} 构成一次探峰。
//       本族纪律（协调者 §2.18.990）：距离类符号只信判题机、单边可信（只信负号），探针读 0 ≠ 没收益。
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106162 <https://duck.ac/submission/106162>
//     用途：**本文件的正文基座**（`work/m1c_sp256.cpp` = #106080 + 对手 #106134 的三处组间距；782.297607 ms）。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #106155 <https://duck.ac/submission/106155>
//     用途：逐字复刻他相对 [1] 的**全部**增量（逐行 diff `ref2/rival_106134.cpp → ref2/rival_106155.cpp`）：
//       `#define KCP (KC + 32)` → **`(KC + 16)`**，并把 mikro 汇编里**全部 18 处** A 行偏移按行号同步减
//       `128 B/行`（多重集精确对应：1280×6→1152、2560×4→2304、3840×4→3456、5120×2→4608、6400×2→5760）
//       ⇒ 即"A 面板行距 pad 与其从属的全部 A 行取数偏移必须同步"（与 `mmml1k` 上 `KC` 减半必须同步减半
//       A 行偏移、以及对手 #105565 的手法完全同型）。
// ======================
// ===== 思路 =====
// **单机制（A 面板行距 pad −16 元素，及其从属的全部汇编行偏移）**：`KCP = KC + 16`、18 处 `vbroadcastsd`
// 偏移按行号 −128 B/行。正文其余一字未动。
// 预期：他这一发在**同版引擎**上拿到 780.126072（−0.28% vs [1]）⇒ 本发复现同量级；
//   且 `mine` 每压低 1 ms，窗口门槛 `mine/1.005` 同步降低 ⇒ **双向有效**（既逼近严支，又加宽窗口）。
// 正确性：纯地址改动 ⇒ 双引擎应逐位相同（`l1n_vfy2.py` n=512/1024/2048）。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #106080 <https://duck.ac/submission/106080>
//     用途：**本文件的正文基座**（`work/m2x_v1_mini.cpp`，六级物化 Winograd，786.401305 ms）。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #106134 <https://duck.ac/submission/106134>
//     用途：逐字复刻他相对 [1] 的**全部**增量。他的抬头原文："retained the accepted six-level
//       Winograd engine"（即 [1]）+ "Space the three upper Winograd scratch-plane groups by 256
//       long long elements and shift the nested leaf workspace past the expanded middle group."
//       代码差异恰好 4 行（3 处组间距 + 1 处 wleaf 偏移）。
// ======================
// ===== 思路 =====
// **单机制（组间距 + 与之耦合的工作区偏移）**：三级上层的 `M5..Q` 组间距 `+32 → +256` 长整型元素
// （`wg_buf1`/`wg_buf2`/`wg_buf` 各一行），并把叶子工作区 `wleaf` 由 `+128` 移到 `+1024`（越过被撑大的那组；
// 这两步**必须耦合**，否则两区重叠）。叶子层自己的 `M5=wleaf,*M6=M5+sz,...`（无 pad）**未动** ✓。
// ★ 本发是**红态弱支配**下的"武装窗口"件：现 `mine=786.401305 / T=782.3078`，宽支门槛 `1.005*T=786.219317`
//   ⇒ **只差 0.181988 ms（0.023%）**。只要把它压过去，窗口即武装；此后对手任何新最好件（他们现就在 782.31）
//   都会当场把我们转绿。对手正是用这同一组改动在**我们这版引擎**上拿到 782.3078（−0.52%）⇒ 最强的先验。
// 正确性：纯地址改动（算术与全部存数表达式不变）⇒ 双引擎应逐位相同；容量已核（wleaf 移位后与 M5..Q 组不重叠）。
// ================
// ===== 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 _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")
#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"
        "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
        "vbroadcastsd 0(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm0\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm1\n\t"
        "vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm2\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm4\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm6\n\t"
        "vpmulld %%ymm13, %%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"
        "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
        "vbroadcastsd 0(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%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"
        "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
        "vbroadcastsd 0(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm0\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm1\n\t"
        "vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm2\n\t"
        "vpmulld %%ymm13, %%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"
        "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
        "vbroadcastsd 0(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%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"
        "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
        "vbroadcastsd 0(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm0\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm1\n\t"
        "vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm2\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm4\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm6\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm7\n\t"
        "vbroadcastsd 4608(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm8\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm9\n\t"
        "vbroadcastsd 5760(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm10\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm11\n\t"
        "prefetcht0 2048(%[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"
        "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
        "vbroadcastsd 0(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vbroadcastsd 1152(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vbroadcastsd 2304(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vbroadcastsd 3456(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vbroadcastsd 4608(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
        "vpaddd %%ymm14, %%ymm9, %%ymm9\n\t"
        "vbroadcastsd 5760(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "prefetcht0 2048(%[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;
  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 + 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;
  /* ---- 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;
      for (int k = 0; k < kc; k++) {
        const long long *r0 = s0 + (size_t)(pc + k) * ld;
        __m256i v0 = LDQ((const __m256i *)(r0));
        __m256i v1 = LDQ((const __m256i *)(r0 + 4));
        __m256i v2 = LDQ((const __m256i *)(r0 + 8));
        __m256i v3 = LDQ((const __m256i *)(r0 + 12));
        if (nt == 2) {
          const long long *r1 = s1 + (size_t)(pc + k) * ld;
          __m256i w0 = LDQ((const __m256i *)(r1));
          __m256i w1 = LDQ((const __m256i *)(r1 + 4));
          __m256i w2 = LDQ((const __m256i *)(r1 + 8));
          __m256i w3 = LDQ((const __m256i *)(r1 + 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 *)(q0 + (size_t)k * 4), v0);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), v1);
        _mm256_store_si256((__m256i *)(q2 + (size_t)k * 4), v2);
        _mm256_store_si256((__m256i *)(q3 + (size_t)k * 4), v3);
      }
    } 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), v0);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), v1);
        _mm256_store_si256((__m256i *)(q2 + (size_t)k * 4), v2);
        _mm256_store_si256((__m256i *)(q3 + (size_t)k * 4), v3);
      }
    }
  }
  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), v);
      }
    } 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), 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), 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. */
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)));

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) + 1024;   /* wg_buf1 slack: grandchild uses 5*256^2 */
static long long wg_buf_mini[5*64*64 + 4*256] __attribute__((aligned(2097152)));
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,*P=M7+sz+256,*Q=P+sz+256;
  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 PRM(X,Y,Z,LX,LY,LZ) do { Opd ao=o1((X),(LX)),bo=o1((Y),(LY)); mmxR(h,&ao,&bo,(Z),(LZ),0,g_scratch+((size_t)1<<20)); } while(0)
  PRM(A11,B11,M1,lda,ldb,ldc);
  PRM(A12,B21,M2,lda,ldb,ldc);
  wg_add(h,P,A21,lda,A22,lda,+1);
  wg_add(h,Q,B12,ldb,B11,ldb,-1);
  PRM(P,Q,M5,h,h,h);
  wg_mutate(h,P,A11,lda,0);
  wg_mutate(h,Q,B22,ldb,1);
  PRM(P,Q,M6,h,h,h);
  wg_mutate(h,P,A12,lda,1);
  PRM(P,B22,M3,h,ldb,ldc);
  wg_mutate(h,Q,B21,ldb,0);
  PRM(A22,Q,M4,lda,h,ldc);
  wg_add(h,P,A11,lda,A21,lda,-1);
  wg_add(h,Q,B22,ldb,B12,ldb,-1);
  PRM(P,Q,M7,h,h,h);
#undef PRM
  wg_combine_cache(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
}

/* 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,*M7=M6+sz,*P=M7+sz,*Q=P+sz;
  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+256,*M7=M6+sz+256,*P=M7+sz+256,*Q=P+sz+256;
  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;
  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_add(h,P,A21,lda,A22,lda,+1);
  wg_add(h,Q,B12,ldb,B11,ldb,-1);
  PROD2(P,Q,M5,h,h,h);
  wg_mutate(h,P,A11,lda,0);
  wg_mutate(h,Q,B22,ldb,1);
  PROD2(P,Q,M6,h,h,h);
  wg_mutate(h,P,A12,lda,1);
  PROD2(P,B22,M3,h,ldb,ldc);
  wg_mutate(h,Q,B21,ldb,0);
  PROD2(A22,Q,M4,lda,h,ldc);
  wg_add(h,P,A11,lda,A21,lda,-1);
  wg_add(h,Q,B22,ldb,B12,ldb,-1);
  PROD2(P,Q,M7,h,h,h);
#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+256,*M7=M6+sz+256,*P=M7+sz+256,*Q=P+sz+256;
  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_add(h,P,A21,n,A22,n,+1);             // S1
  wg_add(h,Q,B12,n,B11,n,-1);             // T1
  PROD(P,Q,M5,h,h,h);
  wg_mutate(h,P,A11,n,0);                 // S2 = S1 - A11
  wg_mutate(h,Q,B22,n,1);                 // T2 = B22 - T1
  PROD(P,Q,M6,h,h,h);
  wg_mutate(h,P,A12,n,1);                 // S4 = A12 - S2
  PROD(P,B22,M3,h,n,n);
  wg_mutate(h,Q,B21,n,0);                 // T4 = T2 - B21
  PROD(A22,Q,M4,n,h,n);
  wg_add(h,P,A11,n,A21,n,-1);             // S3
  wg_add(h,Q,B22,n,B12,n,-1);             // T3
  PROD(P,Q,M7,h,h,h);
#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

__attribute__((weak)) int main() { return 0; }

评测结果