代码

// ===== REFERENCES =====
// [1] 本账号（saffah_cc_v41_agg1）提交 #105839 <https://duck.ac/submission/105839>
//     （112.541260 ms，我方最好 = 红态现行最好件）—— 本文件 = 该件正文 + 一处结构改动。
// [2] duck.ac 用户 saffah_codex_6s_agg2 #105791 <https://duck.ac/submission/105791>
//     （112.512637，T）—— 同族对手（NC->nc 一行）。
// [3] 同族 mmml4k 车道（本会话 `dm_`）的"再物化一级"形态与判据（§2.18.930：
//     刀只在"最内层之上还有一个**未分解**的乘积级"时付钱）。
// ======================
// ===== 思路 =====
// 现状：mine = 112.541260（#105839），T = 112.512637，严支线 0.99*T+1us = 111.388511
//       ⇒ 还缺 1.152749 ms（1.024%）。
//
// ★ 唯一结构改动：**把最内层的 128 级乘积也物化** —— 叶 128 -> 7 x 叶 64。
//   现状链条：winograd_top(1024) -> winograd_mid(512) -> mmxR(256) -> 7 x **leaf(128)**
//             （mmxR 的 MK 宏在 h=128 == CUTL 时直接落"双目的地叶"，128 级**裸算、未分解**）
//   改后：                   ... -> mmxR(256) -> 7 x **st128** -> 7 x **leaf(64)**
//   ⇒ 每个 128 级乘积从 128^3 变成 7 x 64^3 ⇒ **mikro MAC x 7/8 = -12.5%**；而 mikro 占
//   判题机 in-situ 相位账的 **87.3%（346.9 M TSC / 397.4 M）**、且这 346.9 M **全部**来自这
//   343 个 128 级叶 ⇒ 理论上界 = **-10.9%**（远超 1.02% 缺口）。
//
// ★ 为什么必须是"物化"形态（不是"降 CUTL"）：降 CUTL 会让 64 级叶拿到 nt=4 的**四重和式**
//   描述符（pack 源读 x2、且走慢路径）；物化形态先把操作数落成单矩阵（nt=1），再由物化操作数
//   造**二重**和式（nt<=2，pack 走向量快路径）。同族 `l1n_` 的读数：递归形态 +3.4% vs 物化形态 -4.10%。
//
// ★ 实现（见正文 sm1_l4 段）：
//   1) `wg_mat` 把 Opd（nt<=4 的和式）显式落成 128x128 连续矩阵 MA/MB；
//   2) `wg_qadd` 由 MA/MB 的 64x64 象限造 10 个二重和式（sa1/a2/a5/a6/a7、sb1/b3/b4/b6/b7）；
//   3) 7 个 leaf(64) 乘积落进 128x128 临时区 P 的四个子象限，写/加/减语义按 Strassen 映射
//      **预先排定**（代数已逐项验证：C11=M1+M4-M5+M7, C12=M3+M5, C21=M2+M4, C22=M1-M2+M3+M6，
//       其中 M1=(a11+a22)(b11+b22), M2=(a21+a22)b11, M3=a11(b12-b22), M4=a22(b21-b11),
//             M5=(a11+a12)b22, M6=(a21-a11)(b11+b12), M7=(a12-a22)(b21+b22)）；
//   4) 两趟 `cmb` 把 P 加到两个父目的地（带 op/sign），**与父级的 op 语义逐位一致**。
//   5) 另：leaf 内 `mikro(KC / UU, ...)` -> `mikro(kc / UU, ...)`。对现有叶 kc == KC == 128
//      **恒等**（零行为变化），是新叶（k=64）的必要条件：否则 mikro 会白跑 128 步里的后 64 步。
//
// ★ 常驻：st128 的 720 KB 临时区**全部从 `g_scratch` 的死尾切**（§证明见正文），净增常驻 = 0。
//
// ★ 本发同时含**第二条**（kc 维死区，与第一条同向、且对现有叶恒等）：
//   packA 的 `for (; k < KC; k++) *p++ = 0;` -> `p += (KC - k);`；packB 的四处
//   `kc..KC` 零填充块整块删除。理由：mikro 现在按 `cnt = kc/UU` 走 k 步 ⇒ **k >= kc 的
//   面板区永不被读**（A 面板块行距 KC*8 不变、B 分组基址不变 ⇒ 布局不动）。
//   对现有叶 kc == KC == 128 ⇒ **逐字节恒等**（零行为变化）；对新叶（kc=64）去掉约
//   205 MB/run 的死 store（这是本地读数从 +22% 收到 +6% 的主因）。
//
// 正确性闸门：与 [1] 同进程逐元素对拍 n=64/128/256/512/1024 **bad=0**
//   （Winograd/Strassen 在 mod 2^64 下精确；写/加/减语义逐项代数验证）。
//
// 本地读数（**仅作方向参考；本机 load ≈ 14.5/16、同侪二进制占核 ⇒ 本族 probe 与提交会反号**）：
//   交错 A/B，各 15 次：base min 371.8 M / median ≈ 380 M；l4b min 395.4 M / median ≈ 414 M
//   ⇒ 本地 **+6~9%（更慢）**。判题机是唯一裁判 ⇒ 本发 = **定价件/排除发**。
// 未动：mikro 汇编体、packA/packB、面板几何、KC/MC/NC/CUTL、Winograd 各层与间距、预取、THP。
// ================
/* 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 128
#define MC 12
#define NC 256
#define UU 1
#define CUTL 128

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 1024(%[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 2048(%[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 3072(%[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 4096(%[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 5120(%[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 1024(%[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 2048(%[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 3072(%[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 1024(%[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 1024(%[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 2048(%[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 3072(%[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 4096(%[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 5120(%[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);
        }
      }
    }
    p += (KC - k);   /* sm1_l4: k>=kc 的面板区 mikro 不读，删死 store（kc==KC 时恒等） */
  }
}

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;
  }
  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;
  }
  if (nt >= 3) jb0 = NC;
  /* tail: only the `nc` columns actually used by this call.  The old bound was
     the BUFFER capacity NC (256) while nc = 128 for a 128-wide leaf => 32 dead
     4-column groups (KC*4*4 B = 131 KB of zero stores) per call x 343 calls
     = ~45 MB/run written into panel blocks the mikro never reads. */
  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);
      }
    } 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]; }
    }
  }
}


/* ================= sm1_l4: 128 级物化（新叶 = 64） =================
   临时区 = 2*128*128（MA/MB）+ 10*64*64（十个和式）+ 128*128（乘积临时区 P）
          = 90112 元素 = 720 KB，**从 g_scratch 的死尾切**。
   死尾证明（n<=1024，即本题与闸门全部 n）：
     g_scratch 仅供 mmxR 的 T/T2 使用；n=1024 走 winograd_top（不用 g_scratch），
     winograd_mid 的 MPROD 只把 g_scratch 当 mmxR(256) 的 T 传入，而 mmxR(256) 的
     h=128 分支（本发起的 st128）**不读 T**（旧代码那里是 leaf(h) 双目的地直落，同样不读 T）；
     n=512 时最深的触碰是 mmxR(512) 的 T（256*256）+T2 = 2*65536 = 131072 元素。
     => 实际最大触碰 index <= 131072 << 528384-90112-64 = 438208。安全。 */
#define SM1_L4_ELEMS (2*128*128 + 10*64*64 + 128*128)
static long long *const wg_st = g_scratch + (2 * (512*512) + 4096 - SM1_L4_ELEMS - 64);

/* dst(H x H, ld = H) = 描述符 A 的显式和（nt <= 4，向量化） */
static void wg_mat(int H, long long *dst, const Opd *A) {
  const int nt = A->nt, ld = A->ld;
  const long long *p = A->p;
  for (int i = 0; i < H; i++) {
    long long *dr = dst + (size_t)i * H;
    for (int j = 0; j < H; j += 4) {
      __m256i v = LDU(p + (size_t)(A->r[0] + i) * ld + A->c[0] + j);
      if (A->s[0] < 0) v = _mm256_sub_epi64(_mm256_setzero_si256(), v);
      for (int t = 1; t < nt; t++) {
        __m256i w = LDU(p + (size_t)(A->r[t] + i) * ld + A->c[t] + j);
        v = (A->s[t] > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
      }
      STU(dr + j, v);
    }
  }
}
/* dst(64x64, ld=64) = a (+-) b，a/b 是 128x128 物化操作数的象限（ld=128） */
static void wg_qadd(long long *dst, const long long *a, const long long *b, int sg) {
  for (int i = 0; i < 64; i++) {
    const long long *x = a + (size_t)i * 128, *y = b + (size_t)i * 128;
    long long *z = dst + (size_t)i * 64;
    for (int j = 0; j < 64; j += 4) {
      __m256i vx = LDU(x + j), vy = LDU(y + j);
      STU(z + j, sg > 0 ? _mm256_add_epi64(vx, vy) : _mm256_sub_epi64(vx, vy));
    }
  }
}
/* 128 级乘积 Ao*Bo -> 加进 (C1,op1) 与 (C2,op2)。C2==0 表示单目的地。 */
static void st128(const Opd *Ao, const Opd *Bo,
                  long long *C1, int ld1, int op1, long long *C2, int ld2, int op2) {
  long long *MA = wg_st, *MB = MA + 128*128;
  long long *S = MB + 128*128;          /* 10 x (64x64) */
  long long *P = S + 10*64*64;          /* 128x128 乘积临时区 */
  wg_mat(128, MA, Ao); wg_mat(128, MB, Bo);
  const long long *a11 = MA, *a12 = MA + 64, *a21 = MA + 64*128, *a22 = MA + 64*128 + 64;
  const long long *b11 = MB, *b12 = MB + 64, *b21 = MB + 64*128, *b22 = MB + 64*128 + 64;
  long long *sa1=S, *sa2=S+4096, *sa5=S+8192, *sa6=S+12288, *sa7=S+16384;
  long long *sb1=S+20480, *sb3=S+24576, *sb4=S+28672, *sb6=S+32768, *sb7=S+36864;
  wg_qadd(sa1, a11, a22, +1); wg_qadd(sa2, a21, a22, +1); wg_qadd(sa5, a11, a12, +1);
  wg_qadd(sa6, a21, a11, -1); wg_qadd(sa7, a12, a22, -1);
  wg_qadd(sb1, b11, b22, +1); wg_qadd(sb3, b12, b22, -1); wg_qadd(sb4, b21, b11, -1);
  wg_qadd(sb6, b11, b12, +1); wg_qadd(sb7, b21, b22, +1);
  Opd q_a11 = o1(a11, 128), q_a22 = o1(a22, 128), q_b11 = o1(b11, 128), q_b22 = o1(b22, 128);
  Opd q_sa1 = o1(sa1, 64), q_sa2 = o1(sa2, 64), q_sa5 = o1(sa5, 64);
  Opd q_sa6 = o1(sa6, 64), q_sa7 = o1(sa7, 64);
  Opd q_sb1 = o1(sb1, 64), q_sb3 = o1(sb3, 64), q_sb4 = o1(sb4, 64);
  Opd q_sb6 = o1(sb6, 64), q_sb7 = o1(sb7, 64);
  long long *p11 = P, *p12 = P + 64, *p21 = P + 64*128, *p22 = P + 64*128 + 64;
#define L4(X_, Y_, D1, O1, D2, O2) \
  leaf(64, &(X_), &(Y_), (long long *)(D1), 128, (O1), (long long *)(D2), 128, (O2))
  L4(q_sa1, q_sb1, p11, 0, p22, 0);   /* M1 -> +C11 +C22 */
  L4(q_sa2, q_b11, p21, 0, p22, 2);   /* M2 -> +C21 -C22 */
  L4(q_a11, q_sb3, p12, 0, p22, 1);   /* M3 -> +C12 +C22 */
  L4(q_a22, q_sb4, p11, 1, p21, 1);   /* M4 -> +C11 +C21 */
  L4(q_sa5, q_b22, p11, 2, p12, 1);   /* M5 -> -C11 +C12 */
  L4(q_sa6, q_sb6, p22, 1, 0, 0);     /* M6 -> +C22 */
  L4(q_sa7, q_sb7, p11, 1, 0, 0);     /* M7 -> +C11 */
#undef L4
  cmb(C1, ld1, P, 128, 128, +1, op1);
  if (C2) cmb(C2, ld2, P, 128, 128, +1, op2);
}
/* ================= /sm1_l4 ================= */

/* 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 (h == CUTL) { \
      /* sm1_l4: 128 级物化 —— 新叶 64（mikro MAC x7/8） */ \
      st128(&(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; \
    if (h == CUTL) { \
      st128(&(Ma), &(Mb), Q[qi], ld, oa, 0, 0, 0); \
    } else 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 * 32)
static long long wg_buf_store[9 * (512 * 512 + 3072 + 16)] __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*h,*b=M2+(size_t)i*h,*c=M3+(size_t)i*h;
    const long long *d=M4+(size_t)i*h,*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;
    /* 7 条 M 读流各"每行踢一脚"：行内 64 B 步长由硬件跟踪，行间 4 KB 的跳越超出
       L2 流预取器的能力 ⇒ 每行补一条非临时预取（M 只被读一次 ⇒ nta 不污染 cache）。
       代价 7 uop / 行（行体 512 元素）≈ 0.014 uop/元素。 */
    if (i + 4 < h) {
      long long *pa=(long long*)a+(size_t)4*h, *pb=(long long*)b+(size_t)4*h;
      long long *pc=(long long*)c+(size_t)4*h, *pd=(long long*)d+(size_t)4*h;
      long long *pe=(long long*)e+(size_t)4*h, *pf=(long long*)f+(size_t)4*h;
      long long *pg=(long long*)g+(size_t)4*h;
      __builtin_prefetch(pa,0,0); __builtin_prefetch(pb,0,0); __builtin_prefetch(pc,0,0);
      __builtin_prefetch(pd,0,0); __builtin_prefetch(pe,0,0); __builtin_prefetch(pf,0,0);
      __builtin_prefetch(pg,0,0);
    }
    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));
    }
  }
}

/* ---- sm_mid: materialised Winograd at the 512 level (a further layer) ---- */
static long long wg_mid[9 * (256 * 256 + 3072 + 16)] __attribute__((aligned(64)));
static void winograd_mid(int n, const long long *X, int ldx, const long long *Y, int ldy,
                         long long *Z, int ldz) {
  const int g = n >> 1; const size_t z2 = (size_t)g * g;
  long long *m1=wg_mid,*m2=m1+z2+3072,*m3=m2+z2+3072,*m4=m3+z2+3072;
  long long *m5=m4+z2+3072,*m6=m5+z2+3072,*m7=m6+z2+3072,*pp=m7+z2+3072,*qq=pp+z2+3072;
  const long long *X11=X,*X12=X+g,*X21=X+(size_t)g*ldx,*X22=X21+g;
  const long long *Y11=Y,*Y12=Y+g,*Y21=Y+(size_t)g*ldy,*Y22=Y21+g;
#define MPROD(A_,LNA_,B_,LNB_,Z_) do { Opd ao=o1((A_),(LNA_)),bo=o1((B_),(LNB_)); mmxR(g,&ao,&bo,(Z_),g,0,g_scratch); } while(0)
  MPROD(X11,ldx,Y11,ldy,m1);
  MPROD(X12,ldx,Y21,ldy,m2);
  wg_add(g,pp,X21,ldx,X22,ldx,+1);   wg_add(g,qq,Y12,ldy,Y11,ldy,-1);   MPROD(pp,g,qq,g,m5);
  wg_mutate(g,pp,X11,ldx,0);         wg_mutate(g,qq,Y22,ldy,1);         MPROD(pp,g,qq,g,m6);
  wg_mutate(g,pp,X12,ldx,1);         MPROD(pp,g,Y22,ldy,m3);
  wg_mutate(g,qq,Y21,ldy,0);         MPROD(X22,ldx,qq,g,m4);
  wg_add(g,pp,X11,ldx,X21,ldx,-1);   wg_add(g,qq,Y22,ldy,Y12,ldy,-1);   MPROD(pp,g,qq,g,m7);
#undef MPROD
  wg_combine(g,Z,ldz,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=wg_buf,*M2=M1+sz+3072,*M3=M2+sz+3072,*M4=M3+sz+3072;
  long long *M5=M4+sz+3072,*M6=M5+sz+3072,*M7=M6+sz+3072,*P=M7+sz+3072,*Q=P+sz+3072;
  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 { winograd_mid(h,(X),(LX),(Y),(LY),(Z),h); } 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; }

评测结果