代码

// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104259 <https://duck.ac/submission/104259>
//     （119.726 ms）—— 本文件的引擎与顶层（mikro / packA / packB / leaf / mmxR / Opd /
//     `winograd_top` / `wg_add` / `wg_mutate` / `wg_combine`）来自它，逐字节相同。
//     它自列的数学来源：Brice Boyer et al.
//     <https://hal.science/hal-00163141v3/file/winoschedule.pdf>
//     （只用了七乘积方程与输出重组的数学思想；该文无许可证声明，只参考思想、未复制代码）。
// [2] 本账号（saffah_cc_v41_agg1）提交 #104587 <https://duck.ac/submission/104587>（119.691 ms）
//     —— 引擎基座。
// [3] 本账号（saffah_cc_v41_agg1）提交 #104960 <https://duck.ac/submission/104960>（119.371 ms）
//     —— `wg_combine` 的 4 条 C 写改 `vmovntdq` + `sfence`。
// [4] 本账号（saffah_cc_v41_agg1）提交 #104968 <https://duck.ac/submission/104968>（119.162 ms）
//     —— `wg_buf` 2 MB 对齐 + `madvise(MADV_HUGEPAGE)`。**本文件的直接基座**。
// [5] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104974 <https://duck.ac/submission/104974>
//     —— **本发所要回答的那条读数**：给 `wg_combine` 的 7 条 M 读流加"每行一次"的
//     kickstart 预取 ⇒ −0.19%。该发在**未带 NT/THP 的旧基座**上测得 119.459242（对
//     `ml7_skew` 的 119.691479）。它证明的机理是"**该趟的流跟踪跟不上**"。
// ======================
// ===== 思路 =====
// 现状：mine = 119.162383（#104968），T = 119.726395，严支线 0.99*T+1 = 118.530131
//       ⇒ 还需 **0.632252 ms（0.53%）**。
//
// ★ 唯一改动：**把 7 个 M 缓冲改成"按行走"的交错布局，使 `wg_combine` 的读从 7 条并行流
//   变成 1 条顺序流**（算术一字未动）。
//
// 起因（判据链，不是猜）：
//  1. `#104974` 给这 7 条 M 读流各加一条"每行一次"的 kickstart 预取 ⇒ **−0.19%**。
//     软件预取只在**硬件跟踪跟不上**时才兑现（§2.18.831；而本题所有流的步长都是 2 的幂，
//     按"步长非 2 的幂"的判据本来**不该**兑现）⇒ 唯一自洽的解释是**流数**：
//     `wg_combine` 同时有 **7 条 M 读流 + 4 条 C 写流 = 11 条流**。
//  2. 本队自己的相位账把这一趟标定为 22 MB @ ~8.7–10 GB/s（**低效**）。
//  3. ⇒ **减少流数**才是结构解，而不是给每条流补预取（补一条预取只拉 1 条 64 B 行，
//     而行长 4 KB）。
//
// 布局（只动 M 区，P/Q 原位不动）：
//   改前：`M_k[i][j] = wg_buf + k*(h*h+8) + i*h + j`  （7 个独立的 h×h 矩阵）
//   改后：`M_k[i][j] = wg_buf + i*(7*h) + k*h + j`    （**7 个 M 的第 i 行连成一段**）
//   ⇒ 第 i 行 = 连续 7*h 个元素（28 KB），内层 j 循环看到的是**一条一直向前走的流**。
//
// **叶侧一行代码都不用改**：`mmxR(512)` 的四个象限偏移是 `d, d+h2, d+h2*ld, d+h2*ld+h2`
// （h2 = n/2 = 256），所以只要把 `PROD` 传的 `ld` 从 `h` 改成 `7*h`，象限偏移就正好落在交错行上
// （256*3584 = 917504 = 256 行 × 28672 B = M_k 的第 256 行）。总占用不变（7*h*h 个元素，
// 仍在 `wg_buf` 内）。叶的写地址仍全部 32 B 对齐（4096 与 28672 都是 64 的倍数）。
// 未动：mikro 的汇编体、packA/packB、`wg_add`/`wg_mutate`、P/Q、Strassen/Winograd 分解合并式、
//       MR/KC/MC/NC/CUTL、[3] 的 NT 写、[4] 的 THP。
// 正确性：与 [4] 在同一进程内逐元素对拍 n=128/256/512/1024 int64，**bad=0**（work/sm_ilv_vfy.cpp）。
// ================
/* mmml: C = A*B mod 2^64, long long, FUNCTION-style C++ linkage.
   5-ALU-uop mikro (vpmuludq + vpaddq + vpmulld + vpaddd per k x row x 4 cols, no p5 uop)
   plus Strassen with fused operand combination and dual-destination stores
   (no extra memory).  MR=6 KC=256 MC=24 NC=256 CUTL=256 U=1 BSW=0 EXP=0 */
#include <immintrin.h>
#include <stddef.h>
#include <stdint.h>
#include <sys/mman.h>
#ifndef MADV_HUGEPAGE
#define MADV_HUGEPAGE 14
#endif

/* force single-instruction 32-byte loads/stores (gcc 9.3 -O2 splits some
   _mm256_loadu_si256 into vmovdqu xmm + vinserti128 = 2 extra uops) */
__attribute__((target("avx2"))) static inline __m256i LDU(const void *p) {
  __m256i v; __asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i *)p)); return v;
}
__attribute__((target("avx2"))) static inline __m256i LDA(const void *p) {
  __m256i v; __asm__("vmovdqa %1, %0" : "=x"(v) : "m"(*(const __m256i *)p)); return v;
}
__attribute__((target("avx2"))) static inline void STU(void *p, __m256i v) {
  __asm__("vmovdqu %0, %1" ::"x"(v), "m"(*(__m256i *)p) : "memory");
}

/* sm_nt: 32-byte non-temporal store (write-allocate bypass) */
static inline void STNT(void *p, __m256i v) {
  __asm__ volatile("vmovntdq %0, %1" ::"x"(v), "m"(*(__m256i *)p) : "memory");
}

/* [1] #101958 的 mikro 函数体里用它自己的 32 字节存取宏，逐字节带回这里 */
#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")
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC push_options
#pragma GCC target("avx2")

typedef long long TE;

#define MR 6
#define KC 256
#define MC 24
#define NC 256
#define UU 1
#define CUTL 256

static TE Apanel[(size_t)MC * KC * 1 + 64] __attribute__((aligned(64)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(64)));
static __m256i AC[12] __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) {
  long idx = -(cnt << 3);
  long long *apE = (long long *)ap + cnt;
  long long *bpE = (long long *)bp + (cnt << 2);
  if (rows == 6) {
  __asm__ volatile(
    "vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
    "vmovdqa %%ymm0, %%ymm1\n\t"
    "vmovdqa %%ymm0, %%ymm2\n\t"
    "vmovdqa %%ymm0, %%ymm3\n\t"
    "vmovdqa %%ymm0, %%ymm4\n\t"
    "vmovdqa %%ymm0, %%ymm5\n\t"
    "vmovdqa %%ymm0, %%ymm6\n\t"
    "vmovdqa %%ymm0, %%ymm7\n\t"
    "vmovdqa %%ymm0, %%ymm8\n\t"
    "vmovdqa %%ymm0, %%ymm9\n\t"
    "vmovdqa %%ymm0, %%ymm10\n\t"
    "vmovdqa %%ymm0, %%ymm11\n\t"
    "test %[idx], %[idx]\n\t"
    "jz 2f\n\t"
    "vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
    "1:\n\t"
    "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
    "vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
    "vbroadcastsd 4096(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
    "vbroadcastsd 6144(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
    "vbroadcastsd 8192(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm9, %%ymm9\n\t"
    "vbroadcastsd 10240(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
    "vbroadcastsd 4(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "add $8, %[idx]\n\t"
    "jnz 1b\n\t"
    "2:\n\t"
    "vpsrlq $32, %%ymm1, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
    "vpsllq $32, %%ymm1, %%ymm1\n\t"
    "vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
    "vmovdqa %%ymm1, 0(%[ac])\n\t"
    "vpsrlq $32, %%ymm3, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
    "vpsllq $32, %%ymm3, %%ymm3\n\t"
    "vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
    "vmovdqa %%ymm3, 32(%[ac])\n\t"
    "vpsrlq $32, %%ymm5, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
    "vpsllq $32, %%ymm5, %%ymm5\n\t"
    "vpaddq %%ymm4, %%ymm5, %%ymm5\n\t"
    "vmovdqa %%ymm5, 64(%[ac])\n\t"
    "vpsrlq $32, %%ymm7, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
    "vpsllq $32, %%ymm7, %%ymm7\n\t"
    "vpaddq %%ymm6, %%ymm7, %%ymm7\n\t"
    "vmovdqa %%ymm7, 96(%[ac])\n\t"
    "vpsrlq $32, %%ymm9, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
    "vpsllq $32, %%ymm9, %%ymm9\n\t"
    "vpaddq %%ymm8, %%ymm9, %%ymm9\n\t"
    "vmovdqa %%ymm9, 128(%[ac])\n\t"
    "vpsrlq $32, %%ymm11, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm11, %%ymm11\n\t"
    "vpsllq $32, %%ymm11, %%ymm11\n\t"
    "vpaddq %%ymm10, %%ymm11, %%ymm11\n\t"
    "vmovdqa %%ymm11, 160(%[ac])\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
    : [ac] "r"(&AC[0])
    : "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7",
      "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15",
      "cc", "memory");
  } else
  if (rows == 4) {
  __asm__ volatile(
    "vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
    "vmovdqa %%ymm0, %%ymm1\n\t"
    "vmovdqa %%ymm0, %%ymm2\n\t"
    "vmovdqa %%ymm0, %%ymm3\n\t"
    "vmovdqa %%ymm0, %%ymm4\n\t"
    "vmovdqa %%ymm0, %%ymm5\n\t"
    "vmovdqa %%ymm0, %%ymm6\n\t"
    "vmovdqa %%ymm0, %%ymm7\n\t"
    "test %[idx], %[idx]\n\t"
    "jz 2f\n\t"
    "vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
    "1:\n\t"
    "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
    "vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
    "vbroadcastsd 4096(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
    "vbroadcastsd 6144(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
    "vbroadcastsd 4(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "add $8, %[idx]\n\t"
    "jnz 1b\n\t"
    "2:\n\t"
    "vpsrlq $32, %%ymm1, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
    "vpsllq $32, %%ymm1, %%ymm1\n\t"
    "vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
    "vmovdqa %%ymm1, 0(%[ac])\n\t"
    "vpsrlq $32, %%ymm3, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
    "vpsllq $32, %%ymm3, %%ymm3\n\t"
    "vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
    "vmovdqa %%ymm3, 32(%[ac])\n\t"
    "vpsrlq $32, %%ymm5, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
    "vpsllq $32, %%ymm5, %%ymm5\n\t"
    "vpaddq %%ymm4, %%ymm5, %%ymm5\n\t"
    "vmovdqa %%ymm5, 64(%[ac])\n\t"
    "vpsrlq $32, %%ymm7, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
    "vpsllq $32, %%ymm7, %%ymm7\n\t"
    "vpaddq %%ymm6, %%ymm7, %%ymm7\n\t"
    "vmovdqa %%ymm7, 96(%[ac])\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
    : [ac] "r"(&AC[0])
    : "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7",
      "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15",
      "cc", "memory");
  } else
  if (rows == 2) {
  __asm__ volatile(
    "vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
    "vmovdqa %%ymm0, %%ymm1\n\t"
    "vmovdqa %%ymm0, %%ymm2\n\t"
    "vmovdqa %%ymm0, %%ymm3\n\t"
    "test %[idx], %[idx]\n\t"
    "jz 2f\n\t"
    "vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
    "1:\n\t"
    "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
    "vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
    "vbroadcastsd 4(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "add $8, %[idx]\n\t"
    "jnz 1b\n\t"
    "2:\n\t"
    "vpsrlq $32, %%ymm1, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
    "vpsllq $32, %%ymm1, %%ymm1\n\t"
    "vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
    "vmovdqa %%ymm1, 0(%[ac])\n\t"
    "vpsrlq $32, %%ymm3, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
    "vpsllq $32, %%ymm3, %%ymm3\n\t"
    "vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
    "vmovdqa %%ymm3, 32(%[ac])\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
    : [ac] "r"(&AC[0])
    : "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7",
      "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15",
      "cc", "memory");
  } else
  {
  __asm__ volatile(
    "vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
    "vmovdqa %%ymm0, %%ymm1\n\t"
    "vmovdqa %%ymm0, %%ymm2\n\t"
    "vmovdqa %%ymm0, %%ymm3\n\t"
    "vmovdqa %%ymm0, %%ymm4\n\t"
    "vmovdqa %%ymm0, %%ymm5\n\t"
    "vmovdqa %%ymm0, %%ymm6\n\t"
    "vmovdqa %%ymm0, %%ymm7\n\t"
    "vmovdqa %%ymm0, %%ymm8\n\t"
    "vmovdqa %%ymm0, %%ymm9\n\t"
    "vmovdqa %%ymm0, %%ymm10\n\t"
    "vmovdqa %%ymm0, %%ymm11\n\t"
    "test %[idx], %[idx]\n\t"
    "jz 2f\n\t"
    "vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
    "1:\n\t"
    "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
    "vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
    "vbroadcastsd 4096(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
    "vbroadcastsd 6144(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
    "vbroadcastsd 8192(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm9, %%ymm9\n\t"
    "vbroadcastsd 10240(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
    "vbroadcastsd 4(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "add $8, %[idx]\n\t"
    "jnz 1b\n\t"
    "2:\n\t"
    "vpsrlq $32, %%ymm1, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
    "vpsllq $32, %%ymm1, %%ymm1\n\t"
    "vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
    "vmovdqa %%ymm1, 0(%[ac])\n\t"
    "vpsrlq $32, %%ymm3, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
    "vpsllq $32, %%ymm3, %%ymm3\n\t"
    "vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
    "vmovdqa %%ymm3, 32(%[ac])\n\t"
    "vpsrlq $32, %%ymm5, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
    "vpsllq $32, %%ymm5, %%ymm5\n\t"
    "vpaddq %%ymm4, %%ymm5, %%ymm5\n\t"
    "vmovdqa %%ymm5, 64(%[ac])\n\t"
    "vpsrlq $32, %%ymm7, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
    "vpsllq $32, %%ymm7, %%ymm7\n\t"
    "vpaddq %%ymm6, %%ymm7, %%ymm7\n\t"
    "vmovdqa %%ymm7, 96(%[ac])\n\t"
    "vpsrlq $32, %%ymm9, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
    "vpsllq $32, %%ymm9, %%ymm9\n\t"
    "vpaddq %%ymm8, %%ymm9, %%ymm9\n\t"
    "vmovdqa %%ymm9, 128(%[ac])\n\t"
    "vpsrlq $32, %%ymm11, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm11, %%ymm11\n\t"
    "vpsllq $32, %%ymm11, %%ymm11\n\t"
    "vpaddq %%ymm10, %%ymm11, %%ymm11\n\t"
    "vmovdqa %%ymm11, 160(%[ac])\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
    : [ac] "r"(&AC[0])
    : "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7",
      "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15",
      "cc", "memory");
  }
  if (rows > 0) {
    /* batch: all C reads first, then all stores.  c1 and c2 are 4 MB apart for
       the dual-destination M products, i.e. congruent mod 4096, so an
       interleaved load/store pairs them into the 4K-aliasing stall. */
    __m256i E1[6], E2[6];
    if (op1) for (int r = 0; r < rows; r++) E1[r] = LDQ((const __m256i *)(c1 + (size_t)r * ldc1));
    if (c2 && op2) for (int r = 0; r < rows; r++) E2[r] = LDQ((const __m256i *)(c2 + (size_t)r * ldc2));
    for (int r = 0; r < rows; r++) {
      __m256i res = AC[r];
      long long *cp = c1 + (size_t)r * ldc1;
      __m256i v = (op1 == 0) ? res : (op1 == 1 ? _mm256_add_epi64(E1[r], res) : _mm256_sub_epi64(E1[r], res));
      STQ((__m256i *)cp, v);
      if (c2) { cp = c2 + (size_t)r * ldc2;
        v = (op2 == 0) ? res : (op2 == 1 ? _mm256_add_epi64(E2[r], res) : _mm256_sub_epi64(E2[r], res));
        STQ((__m256i *)cp, v); }
    }
  }
}

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) {
  long long *p = Apanel;
  const int nt = o->nt;
  const long long *P = o->p; const int ld = o->ld;
  for (int r = 0; r < MC; r++) {
    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 = LDU(s0 + k);
          if (nt == 2) {
            __m256i w = LDU(s1 + k);
            v = (sgn1 > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
          }
        STU(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 {
        /* nt >= 3 (mmxR 2-layer composition): same 4-wide vector path as nt == 2,
           just more terms.  Signs are per-term constants of the Opd. */
        const long long *sa = P + (size_t)(o->r[1] + ic + r) * ld + o->c[1] + pc;
        const long long *sb = P + (size_t)(o->r[2] + ic + r) * ld + o->c[2] + pc;
        const long long *sc = (nt >= 4) ? P + (size_t)(o->r[3] + ic + r) * ld + o->c[3] + pc : 0;
        const int ga = o->s[1] > 0, gb = o->s[2] > 0, gc = (nt >= 4) ? (o->s[3] > 0) : 1;
        for (; k + 4 <= kc; k += 4) {
          __m256i v = LDU(s0 + k);
          __m256i wa = LDU(sa + k);
          v = ga ? _mm256_add_epi64(v, wa) : _mm256_sub_epi64(v, wa);
          __m256i wb = LDU(sb + k);
          v = gb ? _mm256_add_epi64(v, wb) : _mm256_sub_epi64(v, wb);
          if (nt >= 4) {
            __m256i wc = LDU(sc + k);
            v = gc ? _mm256_add_epi64(v, wc) : _mm256_sub_epi64(v, wc);
          }
          STU(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);
        }
      }
    }
    for (; k < KC; k++) { *p++ = 0; }
  }
}

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;
  /* 8-column fast path: read a full 64-byte cache line of B per row per pass.
     The 4-column path used only 32 B of each 64 B line while the row stride is
     ld*8 bytes, so each line was fetched twice.  Output layout is unchanged. */
  int jb0 = 0;
  if (nt >= 3) {
    /* mmxR 2-layer composition: 3-4 term sources, same layout as the nt == 2 path */
    const long long *s0 = P + (size_t)o->r[0] * ld + o->c[0] + jc;
    const long long *s1 = P + (size_t)o->r[1] * ld + o->c[1] + jc;
    const long long *s2 = P + (size_t)o->r[2] * ld + o->c[2] + jc;
    const long long *s3 = (nt >= 4) ? P + (size_t)o->r[3] * ld + o->c[3] + jc : 0;
    const int g1 = o->s[1] > 0, g2 = o->s[2] > 0, g3 = (nt >= 4) ? (o->s[3] > 0) : 1;
    int ngrp = 0;
    for (int jb = 0; jb + 8 <= nc; jb += 8) ngrp++;
    int k = 0;
    for (; k + 2 <= kc; k += 2) {
      const long long *r0 = s0 + (size_t)(pc + k) * ld, *r0b = s0 + (size_t)(pc + k + 1) * ld;
      const long long *r1 = s1 + (size_t)(pc + k) * ld, *r1b = s1 + (size_t)(pc + k + 1) * ld;
      const long long *r2 = s2 + (size_t)(pc + k) * ld, *r2b = s2 + (size_t)(pc + k + 1) * ld;
      const long long *r3 = (nt >= 4) ? s3 + (size_t)(pc + k) * ld : 0;
      const long long *r3b = (nt >= 4) ? s3 + (size_t)(pc + k + 1) * ld : 0;
      long long *qb = Bpanel + (size_t)k * 4;
      for (int g = 0; g < ngrp; g++) {
        long long *q0 = qb + (size_t)g * 2 * KC * 4;
        __m256i a0 = LDU(r0 + g * 8),  a1 = LDU(r0 + g * 8 + 4);
        __m256i c0 = LDU(r0b + g * 8), c1 = LDU(r0b + g * 8 + 4);
        { __m256i b0 = LDU(r1 + g * 8),  b1 = LDU(r1 + g * 8 + 4);
          __m256i d0 = LDU(r1b + g * 8), d1 = LDU(r1b + g * 8 + 4);
          if (g1) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
                    c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
          else    { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
                    c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); } }
        { __m256i b0 = LDU(r2 + g * 8),  b1 = LDU(r2 + g * 8 + 4);
          __m256i d0 = LDU(r2b + g * 8), d1 = LDU(r2b + g * 8 + 4);
          if (g2) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
                    c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
          else    { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
                    c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); } }
        if (nt >= 4) {
          __m256i b0 = LDU(r3 + g * 8),  b1 = LDU(r3 + g * 8 + 4);
          __m256i d0 = LDU(r3b + g * 8), d1 = LDU(r3b + g * 8 + 4);
          if (g3) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
                    c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
          else    { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
                    c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); }
        }
        STU(q0, a0);
        STU(q0 + 4, c0);
        STU(q0 + (size_t)KC * 4, a1);
        STU(q0 + (size_t)KC * 4 + 4, c1);
      }
    }
    for (; k < kc; k++) {
      const long long *r0 = s0 + (size_t)(pc + k) * ld;
      const long long *r1 = s1 + (size_t)(pc + k) * ld;
      const long long *r2 = s2 + (size_t)(pc + k) * ld;
      const long long *r3 = (nt >= 4) ? s3 + (size_t)(pc + k) * ld : 0;
      long long *qb = Bpanel + (size_t)k * 4;
      for (int g = 0; g < ngrp; g++) {
        const int jb = g * 8;
        long long *q0 = qb + (size_t)(jb / 4) * KC * 4;
        __m256i a0 = LDU(r0 + jb), a1 = LDU(r0 + jb + 4);
        { __m256i b0 = LDU(r1 + jb), b1 = LDU(r1 + jb + 4);
          if (g1) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1); }
          else    { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1); } }
        { __m256i b0 = LDU(r2 + jb), b1 = LDU(r2 + jb + 4);
          if (g2) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1); }
          else    { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1); } }
        if (nt >= 4) { __m256i b0 = LDU(r3 + jb), b1 = LDU(r3 + jb + 4);
          if (g3) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1); }
          else    { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1); } }
        STU(q0, a0);
        STU(q0 + (size_t)KC * 4, a1);
      }
    }
    jb0 = ngrp * 8;
    for (int k = kc; k < KC; k++) {
      long long *qb = Bpanel + (size_t)k * 4;
      for (int g = 0; g < ngrp; g++)
        { long long *q0 = qb + (size_t)g * 8 * KC; STU(q0, _mm256_setzero_si256()); STU(q0 + (size_t)KC * 4, _mm256_setzero_si256()); }
    }
  }
  if (nt <= 2) {
    const long long *s0 = P + (size_t)o->r[0] * ld + o->c[0] + jc;
    const long long *s1 = (nt == 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc : 0;
    int ngrp = 0;
    for (int jb = 0; jb + 8 <= nc; jb += 8) ngrp++;
    int k = 0;
    /* 两行一组：同一块里 k 与 k+1 的记录相邻 32 B = 同一条 64 B line，
       两个 32 B 写背靠背 ⇒ 整行写。布局/输出逐字节不变。 */
    for (; k + 2 <= kc; k += 2) {
      const long long *r0  = s0 + (size_t)(pc + k) * ld;
      const long long *r0b = s0 + (size_t)(pc + k + 1) * ld;
      const long long *r1  = (nt == 2) ? s1 + (size_t)(pc + k) * ld : 0;
      const long long *r1b = (nt == 2) ? s1 + (size_t)(pc + k + 1) * ld : 0;
      long long *qb = Bpanel + (size_t)k * 4;
      for (int g = 0; g < ngrp; g++) {
        long long *q0 = qb + (size_t)g * 2 * KC * 4;
        __m256i a0 = LDU(r0 + g * 8),  a1 = LDU(r0 + g * 8 + 4);
        __m256i c0 = LDU(r0b + g * 8), c1 = LDU(r0b + g * 8 + 4);
        if (nt == 2) {
          __m256i b0 = LDU(r1 + g * 8),  b1 = LDU(r1 + g * 8 + 4);
          __m256i d0 = LDU(r1b + g * 8), d1 = LDU(r1b + g * 8 + 4);
          if (o->s[1] > 0) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
                             c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
          else             { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
                             c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); }
        }
        STU(q0, a0);
        STU(q0 + 4, c0);
        STU(q0 + (size_t)KC * 4, a1);
        STU(q0 + (size_t)KC * 4 + 4, c1);
      }
    }
    for (; k < kc; k++) {
      const long long *r0 = s0 + (size_t)(pc + k) * ld;
      const long long *r1 = (nt == 2) ? s1 + (size_t)(pc + k) * ld : 0;
      long long *qb = Bpanel + (size_t)k * 4;
      for (int g = 0; g < ngrp; g++) {
        const int jb = g * 8;
        long long *q0 = qb + (size_t)(jb / 4) * KC * 4;
        __m256i a0 = LDU(r0 + jb), a1 = LDU(r0 + jb + 4);
        if (nt == 2) {
          __m256i b0 = LDU(r1 + jb), b1 = LDU(r1 + jb + 4);
          if (o->s[1] > 0) { a0 = _mm256_add_epi64(a0, b0); a1 = _mm256_add_epi64(a1, b1); }
          else { a0 = _mm256_sub_epi64(a0, b0); a1 = _mm256_sub_epi64(a1, b1); }
        }
        STU(q0, a0);
        STU(q0 + (size_t)KC * 4, a1);
      }
    }

    jb0 = ngrp * 8;
    for (int k = kc; k < KC; k++) {
      long long *qb = Bpanel + (size_t)k * 4;
      for (int g = 0; g < ngrp; g++)
        { long long *q0 = qb + (size_t)g * 8 * KC; STU(q0, _mm256_setzero_si256()); STU(q0 + (size_t)KC * 4, _mm256_setzero_si256()); }
    }
  }
  if (nt >= 3) jb0 = NC;
  for (int jb = jb0; jb < NC; 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 = LDU(sq + (size_t)(pc + k) * ld);
        if (nt == 2) {
          __m256i w = LDU(sq1 + (size_t)(pc + k) * ld);
          v = (o->s[1] > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
        }
        STU(q + (size_t)k * 4, v);
      }
      for (int k = kc; k < KC; k++) STU(q + (size_t)k * 4, _mm256_setzero_si256());
    } 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 = LDU(t);
        STU(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 * KC * 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;
      }
    }
}

/* ------------------------------------------------------------------ *
 * 本发：单目的地递归 Strassen（移植自同族 mmml4k 的 `mmxR`）。
 * 旧 `mmx` 的守卫里有 `d1`（第二目的地）⇒ **M1..M5 一进来就落 leaf**，
 * 于是"递归"与"双目的地"绑死：本题 n=1024 实际只有 1 层 Strassen。
 * `mmxR` 只有一个目的地，第二目的地用 `cmb`（先算进 scratch 再加减）实现
 * ⇒ 每一层都能继续递归。`Opd` 的 nt <= 4 正好够：顶层 nt=1 → 512 层 nt=2 →
 * 256 层 nt=4 ⇒ 在 256 层落 leaf（守卫 nt > 2）。
 * ------------------------------------------------------------------ */

/* scratch：n=1024 / CUTL=256 时，顶层 mmxR(1024) 需要 T 与 T2 = T + 512*512
   ⇒ 2 * 512^2 个 long。n<=1024 时递归深度最多用到这一对。 */
static long long g_scratch[2 * (512 * 512) + 4096] __attribute__((aligned(64)));

/* dst op= sign * T , (h x ldt) -> (h x ldd) , sign = +-1 , op in {0,1,2} */
static void cmb(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op) {
  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]; }
    }
  }
}

/* single-destination recursive Strassen: d op= Ao*Bo , T = scratch block (>= 2*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 <= CUTL || (n & 1) || Ao->nt > 2 || Bo->nt > 2) {
    leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); 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 (h <= CUTL) { \
      /* sub-call IS a leaf: hand it BOTH destination quadrants directly \
         (one mikro pass, no scratch, no cmb) */ \
      leaf(h, &(Ma), &(Mb), Q[qi], ld, oa, Q[qj], ld, ob); \
    } else 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
}


/* ------------------------------------------------------------------ *
 * 顶层改成 **物化版 Winograd 7 乘积**：
 *   每个乘积的操作数是**单项**（nt=1）的子矩阵/物化和式 ⇒ 512 级递归里的 256 叶
 *   看到的 nt 从 <=4 降到 <=2（pack 的源读项数减半），而 7 个和式（P/Q）在 4 MB 缓冲里
 *   **增量改写**、每个只算一次；最后 7 个乘积在**一趟向量**里合并进 C（不再有 cmb 的 RMW 链）。
 * ------------------------------------------------------------------ */
/* sm_thp: 2 MB-aligned backing store + MADV_HUGEPAGE, so the 18.9 MB of Winograd
   buffers costs ~9 faults instead of ~4700. */
#define SM_WG_ELEMS (9 * 512 * 512 + 9 * 16)
static long long wg_buf_store[SM_WG_ELEMS] __attribute__((aligned(1 << 21)));
static long long *wg_buf = wg_buf_store;
__attribute__((constructor)) static void sm_thp_init(void) {
  uintptr_t a = (uintptr_t)wg_buf_store;
  uintptr_t e = a + (uintptr_t)sizeof(wg_buf_store);
  uintptr_t ha = (a + ((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 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=LDU(x+j), vy=LDU(y+j);
      STU(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=LDU(d+j),vy=LDU(y+j);
      STU(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*(7*h),  *b=M2+(size_t)i*(7*h),  *c=M3+(size_t)i*(7*h);
    const long long *d=M4+(size_t)i*(7*h),  *e=M5+(size_t)i*(7*h),  *f=M6+(size_t)i*(7*h);
    const long long *g=M7+(size_t)i*(7*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=LDU(a+j),p2=LDU(b+j),p3=LDU(c+j),p4=LDU(d+j);
      __m256i p5=LDU(e+j),p6=LDU(f+j),p7=LDU(g+j);
      __m256i u2=_mm256_add_epi64(p1,p6),u3=_mm256_add_epi64(u2,p7);
      __m256i u4=_mm256_add_epi64(u2,p5);
      STNT(c11+j,_mm256_add_epi64(p1,p2));
      STNT(c12+j,_mm256_add_epi64(u4,p3));
      STNT(c21+j,_mm256_sub_epi64(u3,p4));
      STNT(c22+j,_mm256_add_epi64(u3,p5));
    }
  }
}
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;
  /* sm_ilv: M_k[i][j] = wg_buf + i*(7*h) + k*h + j  -- row i of all 7 M's contiguous */
  long long *M1=wg_buf,*M2=M1+h,*M3=M2+h,*M4=M3+h;
  long long *M5=M4+h,*M6=M5+h,*M7=M6+h;              /* bases; wg_combine strides by 7*h */
  long long *P=wg_buf+(size_t)7*sz+8,*Q=P+sz+8;
  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) do { Opd ao=o1((X),(LX)),bo=o1((Y),(LY)); mmxR(h,&ao,&bo,(Z),7*h,0,g_scratch); } while(0)
  PROD(A11,B11,M1,n,n);
  PROD(A12,B21,M2,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);
  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);
  wg_mutate(h,P,A12,n,1);                 /* S4 = A12 - S2 */
  PROD(P,B22,M3,h,n);
  wg_mutate(h,Q,B21,n,0);                 /* T4 = T2 - B21 */
  PROD(A22,Q,M4,n,h);
  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);
#undef PROD
  wg_combine(h,C,n,M1,M2,M3,M4,M5,M6,M7);
  _mm_sfence();
}

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);
  else if (n == 1024) 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; }

评测结果