代码

// ===== REFERENCES =====
// [0] duck.ac 用户 saffah_codex_6s_agg2，提交 #104643
//     <https://duck.ac/submission/104643>（925.515 ms，mmml2k 实时 rank 1）：**本文件正文
//     逐字节取自该提交**（除本文件头部这两段注释与下面【思路】里写明的那一处常量改写）。
//     他这一发 = 本账号 #103453 的正文 + `CUTL 512→256` + `Opd r/c/s[4]→[8]` +
//     `materialize512`（把 nt>2 的虚拟 512 操作数物化成稠密矩阵，让 7 个孙代 packer 只读一个稠密源）。
//     该件自带的引用块一并保留。
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #103453
//     <https://duck.ac/submission/103453>（938.594 ms）：本账号在 mmml2k 的上一发最好件；
//     [0] 的正文包含它的全部内容（k-外层 packB / t-外层 packA / mikro 剥皮 / prefetcht0 3072）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104836
//     <https://duck.ac/submission/104836>（6142.958 ms，mmml4k）：**本条刀法的出处与先例** ——
//     在 mmml4k 的同族引擎上，把 B 面板的**组距**从 `KC*4` 改成 `kc*4` 一发改到 **−7.27%**。
//     本发把同一条刀**移植**到本题（mmml2k）的 [0] 正文上；未取用其任何代码，只移植方法。
// [3] 本账号在 mmml4k 的其它取证件 #104518 / #104634 <https://duck.ac/submission/104634>：
//     同族"物化"结构，用于说明 [0] 的 `materialize512` 与本账号那两发是同一族手法（旁证，未取码）。
// ======================
// ===== 思路 =====
// 【本发唯一改动】B 面板的**组距** `KC*4` → `kc*4`，共 **9 个读写点**（代码段内逐点改写、
//   锚点用精确计数 + assertion，绝不静默失效）：
//     `packB` 主 16 列快路径 q0/q1、两处 `jb2` 配对路径的 q0/q1、尾循环的 q/q1、
//     以及 `leaf` 里传给 `mikro` 的 `bp`。写入与读出用的是**同一个 kc** ⇒ 自洽。
//   ★ `kc == KC` 的调用**逐字节不变**，所以这是零风险的布局等价改写。
//
// 机理（**预取跨度账**，不是省字节）：`mikro` 每次调用从 B 面板**连续**读 `kc*32 B`
//   （k 循环步进 32 B），而组的**地址跨度**由 `KC*4` 元素决定 = `KC*32` 字节 = **16 KB**。
//   [0] 把叶从 512 降到了 **256**（`CUTL 512→256` + `materialize512`）⇒ 每次调用的
//   `kc = min(KC=512, 256) = 256` ⇒ **每组只用到前 8 KB，后 8 KB 是地址空洞**；
//   而汇编里的 `prefetcht0 3072(%[bp])`（mp 每 k 前进 32 B ⇒ 前瞻 96 次迭代）在 k ≥ 160 之后
//   就把目标指进那个**从不被读的空洞**；同时每次调用的前 96 次迭代又只剩硬件预取器 cover。
//   改成 `kc*4` 后组跨度 = 它被读的长度 ⇒ 前后两段都在有用区间（后半段预取到的是**下一组**，
//   正是下一次 `jg` 调用要读的数据）。
// 先例幅度（同族、同一机关、判题机实测）：**本账号 #104836 在 mmml4k 上 = −7.27%**
//   （那边是 2401 个 256 叶；本题是 343 个 256 叶，且 mikro 占比 90% ⇒ 作用面同量级）。
//   ⚠ 本发**不引用**那个幅度作为预期，本发的幅度由本次正式提交定价。
//
// 为什么本题之前的其它搬运都失败、而本条不同：本题 mikro 的**指令侧已到 p01 端口墙**
//   （本 notes 的端口账：11.88 cyc/k 步 vs 理论 12.0），所以任何 mikro 微优化无利可图；
//   本刀**不动一条指令**，只动访存跨度 ⇒ 打的是剩下的停等余额。
//
// 正确性：`work/mml2_chk2.py 2048`（跨版本全矩阵 FNV）—— 基座 `#103453`（`mdl_m2_pf3072.cpp`）、
//   [0] 原文、本发**三者 FNV 逐位相同 `4180d8ca30183001`** ⇒ 输出逐字节不变。
// 本次为**单变量移植实验**：给"B 面板组距/预取跨度"这条轴在 mmml2k 的 [0] 新结构上定价。
// ================
#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
#define KB_LK kc
#pragma GCC target("avx2")

typedef long long TE;

#define MR 6
#define KC 512
#define MC 48
#define NC 512
#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)));
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) {
  if (rows == 6) {
  __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"
    "add $8, %[ap]\n\t"
    "add $32, %[bp]\n\t"
    "dec %[cnt]\n\t"
    "jle 2f\n\t"
    "vmovdqu 0(%[bp]), %%ymm12\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 3072(%[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"
    "vmovdqa %%ymm0, 0(%[ac])\n\t"
    "vmovdqa %%ymm1, 32(%[ac])\n\t"
    "vmovdqa %%ymm2, 64(%[ac])\n\t"
    "vmovdqa %%ymm3, 96(%[ac])\n\t"
    "vmovdqa %%ymm4, 128(%[ac])\n\t"
    "vmovdqa %%ymm5, 160(%[ac])\n\t"
    "vmovdqa %%ymm6, 192(%[ac])\n\t"
    "vmovdqa %%ymm7, 224(%[ac])\n\t"
    "vmovdqa %%ymm8, 256(%[ac])\n\t"
    "vmovdqa %%ymm9, 288(%[ac])\n\t"
    "vmovdqa %%ymm10, 320(%[ac])\n\t"
    "vmovdqa %%ymm11, 352(%[ac])\n\t"
    ""
    : [ap] "+r"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt)
    : [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(
    "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"
    "vmovdqu 0(%[bp]), %%ymm12\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"
    "prefetcht0 3072(%[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"
    "vmovdqa %%ymm0, 0(%[ac])\n\t"
    "vmovdqa %%ymm1, 32(%[ac])\n\t"
    "vmovdqa %%ymm2, 64(%[ac])\n\t"
    "vmovdqa %%ymm3, 96(%[ac])\n\t"
    "vmovdqa %%ymm4, 128(%[ac])\n\t"
    "vmovdqa %%ymm5, 160(%[ac])\n\t"
    "vmovdqa %%ymm6, 192(%[ac])\n\t"
    "vmovdqa %%ymm7, 224(%[ac])\n\t"
    ""
    : [ap] "+r"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt)
    : [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(
    "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"
    "vmovdqu 0(%[bp]), %%ymm12\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"
    "prefetcht0 3072(%[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"
    "vmovdqa %%ymm0, 0(%[ac])\n\t"
    "vmovdqa %%ymm1, 32(%[ac])\n\t"
    "vmovdqa %%ymm2, 64(%[ac])\n\t"
    "vmovdqa %%ymm3, 96(%[ac])\n\t"
    ""
    : [ap] "+r"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt)
    : [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 %[cnt], %[cnt]\n\t"
    "jle 2f\n\t"
    "vmovdqu 0(%[bp]), %%ymm12\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 3072(%[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"
    "vmovdqa %%ymm0, 0(%[ac])\n\t"
    "vmovdqa %%ymm1, 32(%[ac])\n\t"
    "vmovdqa %%ymm2, 64(%[ac])\n\t"
    "vmovdqa %%ymm3, 96(%[ac])\n\t"
    "vmovdqa %%ymm4, 128(%[ac])\n\t"
    "vmovdqa %%ymm5, 160(%[ac])\n\t"
    "vmovdqa %%ymm6, 192(%[ac])\n\t"
    "vmovdqa %%ymm7, 224(%[ac])\n\t"
    "vmovdqa %%ymm8, 256(%[ac])\n\t"
    "vmovdqa %%ymm9, 288(%[ac])\n\t"
    "vmovdqa %%ymm10, 320(%[ac])\n\t"
    "vmovdqa %%ymm11, 352(%[ac])\n\t"
    ""
    : [ap] "+r"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt)
    : [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) {
    __m256i lo = AC[0], hi = AC[1];
    __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 = AC[2], hi = AC[3];
    __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 = AC[4], hi = AC[5];
    __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 = AC[6], hi = AC[7];
    __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 = AC[8], hi = AC[9];
    __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 = AC[10], hi = AC[11];
    __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[8], c[8], s[8];
};
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 {
        /* ===== mdl_: nt >= 3 的 packA 改 t 外层（源侧连续读）=====
           原序 `for k { for t }` 每步只在**每个源行**上搬 32 B，而 nt 条流以 32 B 步长交错
           —— 与 packB 修复前的形态完全相同（那次改 k 外层换来 −2.61%）。
           现在对固定 (r, t) 把 `pc..pc+kc` 整段**连续读完**（kc=512 ⇒ 4 KB 连续）。
           累加器直接用 A 面板本身（它在 L1/L2 常驻，MC*KCP*8 = 197 KB）⇒ 不需要临时缓冲、
           也不多一趟搬运：t=0 覆盖写，t>0 面板自身做 RMW。内核读 A 面板的方式完全不变。 */
        long long *pp = p;
        for (int t = 0; t < nt; t++) {
          const long long *so = P + (size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc;
          int sgn = (o->s[t] < 0);
          int kk = 0;
          if (t == 0) {
            for (; kk + 4 <= kc; kk += 4) {
              __m256i w = LDQ((const __m256i *)(so + kk));
              if (sgn) w = _mm256_sub_epi64(_mm256_setzero_si256(), w);
              STQ((__m256i *)(pp + kk), w);
            }
            for (; kk < kc; kk++) pp[kk] = sgn ? -so[kk] : so[kk];
          } else {
            for (; kk + 4 <= kc; kk += 4) {
              __m256i w = LDQ((const __m256i *)(so + kk));
              __m256i v = LDQ((const __m256i *)(pp + kk));
              STQ((__m256i *)(pp + kk), sgn ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w));
            }
            for (; kk < kc; kk++) pp[kk] = sgn ? pp[kk] - so[kk] : pp[kk] + so[kk];
          }
        }
        k = kc;
      }
    }
    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;
  /* PACK8: 8 columns per k.  The 4-column form loads 32 B at a stride of
     ld*8 (4 KB or more), so every load is a distinct, half-used line and the
     hardware prefetcher cannot see the stream.  Reading jb and jb+4 together
     makes each k one FULL 64 B line, halving the line count and the DRAM
     traffic.  All 8 columns come from one line because jc+jb is a multiple of
     8 (jb steps by 8) and P is 4096-byte aligned.  Byte-identical: the two
     stores write exactly the values the 4-column loop wrote to the same
     addresses. */
  if (nt <= 2) {
    for (; jb + 8 <= nc; jb += 8) {
      long long *q0 = Bpanel + (size_t)(jb / 4) * kc * 4;
      long long *q1 = q0 + (size_t)kc * 4;
      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;
      int k = 0;
      for (; k + 2 <= kc; k += 2) {
        const long long *a0 = sq + (size_t)(pc + k) * ld;
        const long long *a1 = a0 + ld;
        __m256i v0 = LDQ((const __m256i *)a0), u0 = LDQ((const __m256i *)(a0 + 4));
        __m256i v1 = LDQ((const __m256i *)a1), u1 = LDQ((const __m256i *)(a1 + 4));
        if (nt == 2) {
          const long long *b0 = sq1 + (size_t)(pc + k) * ld;
          const long long *b1 = b0 + ld;
          __m256i w0 = LDQ((const __m256i *)b0), x0 = LDQ((const __m256i *)(b0 + 4));
          __m256i w1 = LDQ((const __m256i *)b1), x1 = LDQ((const __m256i *)(b1 + 4));
          if (o->s[1] > 0) { v0 = _mm256_add_epi64(v0, w0); u0 = _mm256_add_epi64(u0, x0);
                             v1 = _mm256_add_epi64(v1, w1); u1 = _mm256_add_epi64(u1, x1); }
          else             { v0 = _mm256_sub_epi64(v0, w0); u0 = _mm256_sub_epi64(u0, x0);
                             v1 = _mm256_sub_epi64(v1, w1); u1 = _mm256_sub_epi64(u1, x1); }
        }
        _mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), v0);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), u0);
        _mm256_store_si256((__m256i *)(q0 + (size_t)(k + 1) * 4), v1);
        _mm256_store_si256((__m256i *)(q1 + (size_t)(k + 1) * 4), u1);
      }
      for (; k < kc; k++) {
        const long long *a = sq + (size_t)(pc + k) * ld;
        __m256i v = LDQ((const __m256i *)a);
        __m256i u = LDQ((const __m256i *)(a + 4));
        if (nt == 2) {
          const long long *b = sq1 + (size_t)(pc + k) * ld;
          __m256i w = LDQ((const __m256i *)b);
          __m256i x = LDQ((const __m256i *)(b + 4));
          if (o->s[1] > 0) { v = _mm256_add_epi64(v, w); u = _mm256_add_epi64(u, x); }
          else             { v = _mm256_sub_epi64(v, w); u = _mm256_sub_epi64(u, x); }
        }
        _mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), v);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), u);
      }
    }
  }
  /* ===== mdl_: nt >= 3 的 packB 改成 k 外层（源侧连续读）=====
     原序 `for kb { for k { for tt } }` 对固定 (tt, k) 只读源矩阵的 64 B（8 列），
     下一个 k 就跳到 ld*8 = 16 KB 之外 ⇒ 每条流的步长等于一个 4 KB 页的 4 倍。
     改成 k 外层后，对固定 (tt, k) 会连续读完 jc..jc+nc-1 的**整行**（nc=512 ⇒ 4 KB 连续），
     跨 k 才跳 ld*8 ⇒ 步数不变但每步搬运量 ×64。
     ★ 写入侧补偿：把相邻两个 k 配对，使面板存储仍是整条 64 B line
     （q0 收 k 与 k+1 的 v0/u0，q1 收 v1/u1）——否则 32 B 的散写会退回 RFO 半行写。
     目标地址与 stock 逐字节一致（列组 g 在 Bpanel + g*KC*4），输出因此不变。 */
  if (nt >= 3 && (nc & 7) == 0) {
    int k = 0;
    for (; k + 2 <= kc; k += 2) {
      for (int jb2 = 0; jb2 < nc; jb2 += 8) {
        long long *q0 = Bpanel + (size_t)(jb2 / 4) * kc * 4;
        long long *q1 = q0 + (size_t)kc * 4;
        __m256i v0 = _mm256_setzero_si256(), v1 = _mm256_setzero_si256();
        __m256i u0 = _mm256_setzero_si256(), u1 = _mm256_setzero_si256();
        for (int tt = 0; tt < nt; tt++) {
          const long long *sp = P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb2;
          __m256i w0 = LDQ((const __m256i *)sp), w1 = LDQ((const __m256i *)(sp + 4));
          __m256i x0 = LDQ((const __m256i *)(sp + ld)), x1 = LDQ((const __m256i *)(sp + ld + 4));
          if (o->s[tt] < 0) { v0=_mm256_sub_epi64(v0,w0); v1=_mm256_sub_epi64(v1,w1);
                              u0=_mm256_sub_epi64(u0,x0); u1=_mm256_sub_epi64(u1,x1); }
          else              { v0=_mm256_add_epi64(v0,w0); v1=_mm256_add_epi64(v1,w1);
                              u0=_mm256_add_epi64(u0,x0); u1=_mm256_add_epi64(u1,x1); }
        }
        _mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), v0);
        _mm256_store_si256((__m256i *)(q0 + (size_t)(k + 1) * 4), u0);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), v1);
        _mm256_store_si256((__m256i *)(q1 + (size_t)(k + 1) * 4), u1);
      }
    }
    for (; k < kc; k++) {
      for (int jb2 = 0; jb2 < nc; jb2 += 8) {
        long long *q0 = Bpanel + (size_t)(jb2 / 4) * kc * 4;
        long long *q1 = q0 + (size_t)kc * 4;
        __m256i v0 = _mm256_setzero_si256(), v1 = _mm256_setzero_si256();
        for (int tt = 0; tt < nt; tt++) {
          const long long *sp = P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb2;
          __m256i w0 = LDQ((const __m256i *)sp), w1 = LDQ((const __m256i *)(sp + 4));
          if (o->s[tt] < 0) { v0=_mm256_sub_epi64(v0,w0); v1=_mm256_sub_epi64(v1,w1); }
          else              { v0=_mm256_add_epi64(v0,w0); v1=_mm256_add_epi64(v1,w1); }
        }
        _mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), v0);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), v1);
      }
    }
    jb = jbn;
  }
  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) {
      int wide = (jb + 8 <= nc);
      if (wide) {
        long long *q1 = q + (size_t)kc * 4;
        for (int kb = 0; kb < kc; kb += (KB_LK)) {
          int ke = kb + (KB_LK); if (ke > kc) ke = kc;
          for (int k = kb; k < ke; k++) {
            __m256i v0 = _mm256_setzero_si256(), v1 = _mm256_setzero_si256();
            for (int tt = 0; tt < nt; tt++) {
              const long long *sp = P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb;
              __m256i w0 = LDQ((const __m256i *)sp), w1 = LDQ((const __m256i *)(sp + 4));
              if (o->s[tt] < 0) { v0 = _mm256_sub_epi64(v0, w0); v1 = _mm256_sub_epi64(v1, w1); }
              else              { v0 = _mm256_add_epi64(v0, w0); v1 = _mm256_add_epi64(v1, w1); }
            }
            _mm256_store_si256((__m256i *)(q + (size_t)k * 4), v0);
            _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), v1);
          }
        }
        jb += 4;
      } else {
      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]; }
    }
  }
}

/* A virtual 512 operand with three or four terms is reused by seven
   grandchildren. Form it once so their packers read one dense source. */
static long long matA512[512u*512u] __attribute__((aligned(64)));
static long long matB512[512u*512u] __attribute__((aligned(64)));
static void materialize512(const Opd *o, long long *dst) {
  for (int i=0; i<512; ++i) for (int j=0; j<512; j+=4) {
    __m256i v = _mm256_setzero_si256();
    for (int t=0; t<o->nt; ++t) {
      const long long *p = o->p + (size_t)(o->r[t]+i)*o->ld + o->c[t]+j;
      __m256i w = LDQ((const __m256i*)p);
      v = o->s[t] < 0 ? _mm256_sub_epi64(v,w) : _mm256_add_epi64(v,w);
    }
    STQ((__m256i*)(dst+(size_t)i*512+j),v);
  }
}

/* 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 == 512 && (Ao->nt > 2 || Bo->nt > 2)) {
    Opd da=*Ao, db=*Bo;
    if (Ao->nt > 2) { materialize512(Ao,matA512); da=o1(matA512,512); }
    if (Bo->nt > 2) { materialize512(Bo,matB512); db=o1(matB512,512); }
    mmxR(n,&da,&db,d,ld,op,T);
    return;
  }
  if (n <= CUTL || (n & 1) || Ao->nt > 4 || Bo->nt > 4) {
    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 (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
}

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; }
  mmxR(n, &Ao, &Bo, C, n, 0, g_scratch);
}
#pragma GCC pop_options

评测结果