代码

/*
References:
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/104753: copied our accepted mmml4k top-level Winograd implementation, retaining all inherited credits and license notes below.
- Brice Boyer et al., "Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm", https://hal.science/hal-00163141v3/file/winoschedule.pdf: the same seven-product mathematical equations, independently applied to a second level; no code or text copied.
Approach:
- Apply the existing output-backed Winograd schedule again inside each 2048 product. A separate 40MB child workspace preserves the parent operands and products while the 1024 matrix calls use the original recursive Strassen engine. Explicit row strides support parent output quadrants and materialized inputs.
Purpose:
- Experimental judged speed and correctness measurement of two Winograd levels on mmml4k.
*/
/*
References:
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/104714: copied our latest accepted mmml4k engine, including original source credits and license notes below.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/104259: adapted the top-level Winograd seven-product schedule and modular AVX2 operand preparation from our accepted mmml1k solution.
- Brice Boyer et al., "Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm", https://hal.science/hal-00163141v3/file/winoschedule.pdf: mathematical equations only; no paper code or text copied.
Approach:
- At n=4096, precompute top-level Winograd sums and run seven 2048 products using the inherited recursive engine. Four products write directly to output quadrants; three use temporary planes. Recombine by loading all seven products before each four-quadrant store, retaining modular 64-bit arithmetic while limiting extra memory.
Purpose:
- Experimental judged measurement of a memory-conscious top-level Winograd schedule, following full-size output equivalence checks.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/104652: copied its fastest accepted mmml4k source, including the 512-block materialization and all inherited credits below. Its public source has no separate license notice.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/104263: our earlier 1024-block selective materialization, inherited through #104652.
Approach:
- Also materialize a 1024-block operand when it contains two signed source terms. A single dense operand lets more 512 children retain the extra Strassen recursion to 256, reducing multiplication count. Keep the newer 512-block materialization and all modular AVX2 arithmetic unchanged.
Purpose:
- Experimental judged measurement of the expanded 1024 operand materialization threshold on mmml4k.
*/
// ===== REFERENCES =====
// [0] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104518
//     <https://duck.ac/submission/104518>（6757.731 ms，本题实时 rank 1）：本件的**直接基座**，
//     正文逐字节取自它。
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104263
//     <https://duck.ac/submission/104263>（6765.498 ms）：#104518 的来源。该提交首创
//     "在 n=1024 处把 4 项惰性操作数物化成连续矩阵、从而让递归再下推两层" 这一结构，
//     `materialize1024` 函数体**直接复制**自它（本发把它泛化成 `materializeN(o,dst,n)`
//     并在 n=512 复用）。其自带的引擎血统引用（#103986 / #103919）一并继承。
// [2] 本账号提交 #104263 的血统说明见 [1]；本发**没有**引用任何其它他人代码。
// ======================
// ===== 思路 =====
// 【本发唯一改动】把 [1] 在 **n=1024** 处的那把"物化"刀**推广到 n=512**：
//   `mmxR` 的守卫（`n<=CUTL || n&1 || nt>2`）**之后**再加一个钩子：
//     若 `n == 512 && (Ao->nt > 1 || Bo->nt > 1)` ⇒ 先把惰性操作数物化成 512² 连续矩阵
//     （`materializeN`，目的地是**新加的 2 MB 静态缓冲** `material_C/D`），描述符退化成 `o1()`，
//     再递归。守卫已过 ⇒ 本次**本来就会递归** ⇒ **叶结构（叶数/叶尺寸/叶的位置）一字不变**。
//
// 为什么这一步应该赚（这是"pack 流量"账，不是 MAC 账 —— MAC 一字未变）：
//   * `leaf(n)` 里 `packA`/`packB` 每个元素的**源读项数 = 操作数的 nt**。
//     n=512 这一层的操作数 nt ∈ {1,2}，递归到 h=256 时子问题描述符被 `oadd` 加成 nt ∈ {2,4}
//     ⇒ **256 叶的 pack 源读是 2~4 项**。
//   * 物化后 nt 掉到 1 ⇒ 256 叶的描述符 nt **只有 2**（`oadd` 两个单项）⇒ **pack 源读项数减半**。
//     按 `work/l5_tree.py` 的树叶枚举（41 个 1024 级钩子结点 → 287 个 n=512 结点 → 每个 7 个 256 叶），
//     每个 256 叶的 pack 流量 5.2 MB → 3.1 MB（含写），而物化本身只花 6 MB/结点
//     ⇒ 每个 n=512 结点净省 ≈ 8.7 MB ⇒ 全题 ≈ −2.5 GB 的串行搬运。
//   * 依据：本发进场前已用判题机把每档叶定价（`work/l5_probe_gen.py`：`leaf(1024)=533.96M`、
//     `leaf(512)=71.12M`、`leaf(256)=11.02M` ticks），并用 `41×leaf(1024)+56×leaf(512)` 对
//     #103166 的 7252.095 ms 标定出 **F = 3.568e9 ticks/s**；据此反解出 #104518 的树叶
//     = `176×leaf(512) + 1169×leaf(256)`，其中 **256 叶实测 9.92 M ticks vs mikro 端口下界 8.39 M**
//     ⇒ 多出来的 1.53 M ≈ 该叶 pack 的 5.2 MB / (4 B/cyc) ⇒ **pack 在本题是串行的**，正是本发打的地方。
//   * 端口账（说明为什么只能从 pack 下手）：mikro 每迭代 24 条 p0/p1 类向量 uop ⇒ **12 周期硬下界**，
//     实测 leaf(1024) 11.94（100%）、leaf(512) 12.72（94%）、leaf(256) 14.19（85%）⇒ 指令侧无余量。
//
// 恒等性与安全性：
//   * 语义恒等：mod 2^64 下加法可结合 ⇒ 先求和再乘 ≡ 先乘再各自求和，**输出逐位不变**。
//   * ★ **缓冲不能复用**：1024 级钩子递归期间仍在读 `material_A/B`，若 512 级写到同一块会**就地覆盖**
//     ⇒ 本发新开 `material_C/D`（各 512²×8 = 2 MB）。（n=512 结点内不会再触发更深的物化 ⇒ 一对足够。）
//   * 正确性：`work/ml4_vfy.cpp`（mod 2^64 逐元素全矩阵对拍）**n=512 / 1024 / 2048 全部 bad=0**；
//     本机 n=1024 全矩阵耗时 304.4 ms → 266.7 ms（**−12.4%**，本机数字仅作方向参考，
//     本机对结构类改动会反号，判题机值由本发定价）。
// 本次为**单变量结构提交**：给"把物化下推一级以砍 pack 源读"这条轴定价。
// ================
#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>
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC push_options
#pragma GCC target("avx2")

typedef long long TE;

#define MR 6
#define KC 512
#define MC 12
#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 + 8)
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 4160(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm2\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd 8320(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm4\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd 12480(%[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 4160(%[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 8320(%[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 12480(%[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 4160(%[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 4160(%[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 4160(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm2\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm3\n\t"
        "vbroadcastsd 8320(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm4\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm5\n\t"
        "vbroadcastsd 12480(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm6\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm7\n\t"
        "vbroadcastsd 16640(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm8\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm9\n\t"
        "vbroadcastsd 20800(%[ap]), %%ymm14\n\t"
        "vpmuludq %%ymm14, %%ymm12, %%ymm10\n\t"
        "vpmulld %%ymm13, %%ymm14, %%ymm11\n\t"
        "prefetcht0 4096(%[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 4160(%[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 8320(%[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 12480(%[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 16640(%[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 20800(%[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 4096(%[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(64)));
static long long material_B[1024u * 1024u] __attribute__((aligned(64)));
/* 512 级物化的目的地：必须与 material_A/B 分开 —— 1024 级钩子的递归仍在读 material_A/B。 */
static long long material_C[512u * 512u] __attribute__((aligned(64)));
static long long material_D[512u * 512u] __attribute__((aligned(64)));
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;
    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 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,*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 PROD2(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); } 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(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,*M7=M6+sz,*P=M7+sz,*Q=P+sz;
  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(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 mmxR(n, &Ao, &Bo, C, n, 0, g_scratch);
}
#pragma GCC pop_options

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

评测结果