代码

// ===== [ml1y6_] mmml1k **同码带尺（对照锚）**（单变量：折叠拼写，机器码逐字节不变）=====
//
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1 提交 #113923 <https://duck.ac/submission/113923>
//     （110.774978 ms = 本题现役最好件 `mine`）—— 本文件正文 = 该提交的正文（经 `work/ml1y5_base.cpp`
//     承载，`.text` md5 `fd2e18f552d1a81848763a031f5dead3`、74,799 B，与该提交逐字节相同）。
// [2] duck.ac 用户 saffah_codex_6s_agg2 提交 #110273 <https://duck.ac/submission/110273>
//     （111.019302 ms = 本题实时 T；严格支线阈 `0.99*T+1µs` = 109.910109 ms）—— 只用于标定缺口。
// [3] 本仓 `problems/mmml1k/work/ml1y5_pricing.md` §一 —— 在档 4 发同码锚（#116437/#116447/#116457）的
//     批内带 = 9.176 µs；本席（ml1y6_）据此把本批窗口钉在 `mine` 上（(1108) 配对设计）。
// 合规：duck.ac 提交正文按站点规则公开可见；本文件按 RULES.md 第 3 节署名作者账号与原地址。
// ======================
// ===== 思路 =====
// （ml1y6_ 席）**用途 = 对照锚**：本批（ml1y6_）只发了一条真臂（K1，见 `work/ml1y6_k1.cpp`），
// 它是**同码带尺**：唯一改动 = 一个宏的折叠拼写（KC 之外的第二个：`NC`），常量折叠等价 ⇒ 机器码逐字节不变
// （本地闸门已核 `.text` md5 = `fd2e18f552d1a81848763a031f5dead3`，四档 n=128/256/512/1024 逐位同）。
// ⇒ 它与 `mine` 的差 = **当批窗口 + 带幅**，用来判 K1 的读数（111.030071）里有多少是窗口上浮。
// ★ 本席在档交班（ml1y5_）："带锚 = 臂的同码孪生 ≥3；对照锚 = 基座的同码孪生 ≥1" ⇒ 本批补 ≥1 对照锚 ✓
// ============================================================================================

/* 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 0xC
#define NC (1 << 7)
#define UU 1
#define CUTL (2 * 64)

static TE Apanel[(size_t)MC * KC * 1 + 64] __attribute__((aligned(4096)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(4096)));
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) {
  register __m256i R0 __asm__("ymm1");
  register __m256i R1 __asm__("ymm3");
  register __m256i R2 __asm__("ymm5");
  register __m256i R3 __asm__("ymm7");
  register __m256i R4 __asm__("ymm9");
  register __m256i R5 __asm__("ymm11");
  __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"
    "prefetcht0 0x40(%[ap],%[idx],1)\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"
    "vpsrlq $32, %%ymm3, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
    "vpsllq $32, %%ymm3, %%ymm3\n\t"
    "vpaddq %%ymm2, %%ymm3, %%ymm3\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"
    "vpsrlq $32, %%ymm7, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
    "vpsllq $32, %%ymm7, %%ymm7\n\t"
    "vpaddq %%ymm6, %%ymm7, %%ymm7\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"
    "vpsrlq $32, %%ymm11, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm11, %%ymm11\n\t"
    "vpsllq $32, %%ymm11, %%ymm11\n\t"
    "vpaddq %%ymm10, %%ymm11, %%ymm11\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx),
      [R0] "=x"(R0), [R1] "=x"(R1), [R2] "=x"(R2), [R3] "=x"(R3), [R4] "=x"(R4), [R5] "=x"(R5)
    :
    : "ymm0",  "ymm2",  "ymm4",  "ymm6", 
      "ymm8",  "ymm10",  "ymm12", "ymm13", "ymm14", "ymm15",
      "cc", "memory");

  __m256i D0,D1,D2,D3,D4,D5;
  if(c2 && op2 != 0){
    D0=LDQ((const __m256i*)(c2+(size_t)0*ldc2));
    D1=LDQ((const __m256i*)(c2+(size_t)1*ldc2));
    D2=LDQ((const __m256i*)(c2+(size_t)2*ldc2));
    D3=LDQ((const __m256i*)(c2+(size_t)3*ldc2));
    D4=LDQ((const __m256i*)(c2+(size_t)4*ldc2));
    D5=LDQ((const __m256i*)(c2+(size_t)5*ldc2));
  }
    if (op1 == 0) {
      STQ((__m256i*)(c1 + (size_t)0*ldc1), R0);
      STQ((__m256i*)(c1 + (size_t)1*ldc1), R1);
      STQ((__m256i*)(c1 + (size_t)2*ldc1), R2);
      STQ((__m256i*)(c1 + (size_t)3*ldc1), R3);
      STQ((__m256i*)(c1 + (size_t)4*ldc1), R4);
      STQ((__m256i*)(c1 + (size_t)5*ldc1), R5);
    } else {
      __m256i E1_0=LDQ((const __m256i*)(c1+(size_t)0*ldc1));
      __m256i E1_1=LDQ((const __m256i*)(c1+(size_t)1*ldc1));
      __m256i E1_2=LDQ((const __m256i*)(c1+(size_t)2*ldc1));
      __m256i E1_3=LDQ((const __m256i*)(c1+(size_t)3*ldc1));
      __m256i E1_4=LDQ((const __m256i*)(c1+(size_t)4*ldc1));
      __m256i E1_5=LDQ((const __m256i*)(c1+(size_t)5*ldc1));
      if (op1 == 1) {
        STQ((__m256i*)(c1+(size_t)0*ldc1), _mm256_add_epi64(E1_0,R0));
        STQ((__m256i*)(c1+(size_t)1*ldc1), _mm256_add_epi64(E1_1,R1));
        STQ((__m256i*)(c1+(size_t)2*ldc1), _mm256_add_epi64(E1_2,R2));
        STQ((__m256i*)(c1+(size_t)3*ldc1), _mm256_add_epi64(E1_3,R3));
        STQ((__m256i*)(c1+(size_t)4*ldc1), _mm256_add_epi64(E1_4,R4));
        STQ((__m256i*)(c1+(size_t)5*ldc1), _mm256_add_epi64(E1_5,R5));
      } else {
        STQ((__m256i*)(c1+(size_t)0*ldc1), _mm256_sub_epi64(E1_0,R0));
        STQ((__m256i*)(c1+(size_t)1*ldc1), _mm256_sub_epi64(E1_1,R1));
        STQ((__m256i*)(c1+(size_t)2*ldc1), _mm256_sub_epi64(E1_2,R2));
        STQ((__m256i*)(c1+(size_t)3*ldc1), _mm256_sub_epi64(E1_3,R3));
        STQ((__m256i*)(c1+(size_t)4*ldc1), _mm256_sub_epi64(E1_4,R4));
        STQ((__m256i*)(c1+(size_t)5*ldc1), _mm256_sub_epi64(E1_5,R5));
      }
    }
  if(c2){
    if(op2==0){
      STQ((__m256i*)(c2+(size_t)0*ldc2),R0);
      STQ((__m256i*)(c2+(size_t)1*ldc2),R1);
      STQ((__m256i*)(c2+(size_t)2*ldc2),R2);
      STQ((__m256i*)(c2+(size_t)3*ldc2),R3);
      STQ((__m256i*)(c2+(size_t)4*ldc2),R4);
      STQ((__m256i*)(c2+(size_t)5*ldc2),R5);
    }else if(op2==1){
      STQ((__m256i*)(c2+(size_t)0*ldc2),_mm256_add_epi64(D0,R0));
      STQ((__m256i*)(c2+(size_t)1*ldc2),_mm256_add_epi64(D1,R1));
      STQ((__m256i*)(c2+(size_t)2*ldc2),_mm256_add_epi64(D2,R2));
      STQ((__m256i*)(c2+(size_t)3*ldc2),_mm256_add_epi64(D3,R3));
      STQ((__m256i*)(c2+(size_t)4*ldc2),_mm256_add_epi64(D4,R4));
      STQ((__m256i*)(c2+(size_t)5*ldc2),_mm256_add_epi64(D5,R5));
    }else{
      STQ((__m256i*)(c2+(size_t)0*ldc2),_mm256_sub_epi64(D0,R0));
      STQ((__m256i*)(c2+(size_t)1*ldc2),_mm256_sub_epi64(D1,R1));
      STQ((__m256i*)(c2+(size_t)2*ldc2),_mm256_sub_epi64(D2,R2));
      STQ((__m256i*)(c2+(size_t)3*ldc2),_mm256_sub_epi64(D3,R3));
      STQ((__m256i*)(c2+(size_t)4*ldc2),_mm256_sub_epi64(D4,R4));
      STQ((__m256i*)(c2+(size_t)5*ldc2),_mm256_sub_epi64(D5,R5));
    }
  }
  return;
  } 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"
    "prefetcht0 64(%[ap],%[idx],1)\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"
    "prefetcht0 64(%[ap],%[idx],1)\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"
    "prefetcht0 64(%[ap],%[idx],1)\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;
  /* [l1y_] nt 特化：把 nt 判决提到 r 环之外，热循环体内 0 条 nt 测试。
     nt>=3 的通用路在下面**原样保留**（冷分支），以满足一般 n 的正确性。 */
  if (nt == 1) {
    const long long *s0b = P + (size_t)o->c[0] + pc;
    const long r0 = o->r[0];
    for (int r = 0; r < MC; r++) {
      int k = 0;
      if (r < mc) {
        const long long *s0 = s0b + (size_t)(r0 + ic + r) * ld;
        for (; k + 4 <= kc; k += 4) { STU(p, LDU(s0 + k)); p += 4; }
        for (; k < kc; k++) *p++ = s0[k];
      }
      for (; k < KC; k++) *p++ = 0;
    }
    return;
  }
  if (nt == 2) {
    const long long *s0b = P + (size_t)o->c[0] + pc;
    const long long *s1b = P + (size_t)o->c[1] + pc;
    const long r0 = o->r[0], r1 = o->r[1];
    if (o->s[1] > 0) {
      for (int r = 0; r < MC; r++) {
        int k = 0;
        if (r < mc) {
          const long long *s0 = s0b + (size_t)(r0 + ic + r) * ld;
          const long long *s1 = s1b + (size_t)(r1 + ic + r) * ld;
          for (; k + 4 <= kc; k += 4)
            { STU(p, _mm256_add_epi64(LDU(s0 + k), LDU(s1 + k))); p += 4; }
          for (; k < kc; k++) *p++ = s0[k] + s1[k];
        }
        for (; k < KC; k++) *p++ = 0;
      }
    } else {
      for (int r = 0; r < MC; r++) {
        int k = 0;
        if (r < mc) {
          const long long *s0 = s0b + (size_t)(r0 + ic + r) * ld;
          const long long *s1 = s1b + (size_t)(r1 + ic + r) * ld;
          for (; k + 4 <= kc; k += 4)
            { STU(p, _mm256_sub_epi64(LDU(s0 + k), LDU(s1 + k))); p += 4; }
          for (; k < kc; k++) *p++ = s0[k] - s1[k];
        }
        for (; k < KC; k++) *p++ = 0;
      }
    }
    return;
  }
  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;
  if (nc == NC && kc == KC && nt <= 2) {
    const long long *s0 = P + (size_t)o->r[0] * ld + o->c[0] + jc;
    if (nt == 1) {
      for (int k = 0; k < KC; k += 2) {
        const long long *r0 = s0 + (size_t)(pc + k) * ld;
        const long long *r0b = r0 + ld;
        long long *qb = Bpanel + (size_t)k * 4;
        for (int g = 0; g < NC / 8; ++g) {
          long long *q = qb + (size_t)g * 2 * KC * 4;
          STU(q, LDU(r0 + g * 8));
          STU(q + 4, LDU(r0b + g * 8));
          STU(q + (size_t)KC * 4, LDU(r0 + g * 8 + 4));
          STU(q + (size_t)KC * 4 + 4, LDU(r0b + g * 8 + 4));
        }
      }
    } else {
      const long long *s1 = P + (size_t)o->r[1] * ld + o->c[1] + jc;
      if (o->s[1] > 0) {
        for (int k = 0; k < KC; k += 2) {
          const long long *r0 = s0 + (size_t)(pc + k) * ld;
          const long long *r1 = s1 + (size_t)(pc + k) * ld;
          long long *qb = Bpanel + (size_t)k * 4;
          for (int g = 0; g < NC / 8; ++g) {
            long long *q = qb + (size_t)g * 2 * KC * 4;
            STU(q, _mm256_add_epi64(LDU(r0 + g*8), LDU(r1 + g*8)));
            STU(q + 4, _mm256_add_epi64(LDU(r0 + ld + g*8), LDU(r1 + ld + g*8)));
            STU(q + (size_t)KC*4, _mm256_add_epi64(LDU(r0 + g*8+4), LDU(r1 + g*8+4)));
            STU(q + (size_t)KC*4 + 4, _mm256_add_epi64(LDU(r0 + ld + g*8+4), LDU(r1 + ld + g*8+4)));
          }
        }
      } else {
        for (int k = 0; k < KC; k += 2) {
          const long long *r0 = s0 + (size_t)(pc + k) * ld;
          const long long *r1 = s1 + (size_t)(pc + k) * ld;
          long long *qb = Bpanel + (size_t)k * 4;
          for (int g = 0; g < NC / 8; ++g) {
            long long *q = qb + (size_t)g * 2 * KC * 4;
            STU(q, _mm256_sub_epi64(LDU(r0 + g*8), LDU(r1 + g*8)));
            STU(q + 4, _mm256_sub_epi64(LDU(r0 + ld + g*8), LDU(r1 + ld + g*8)));
            STU(q + (size_t)KC*4, _mm256_sub_epi64(LDU(r0 + g*8+4), LDU(r1 + g*8+4)));
            STU(q + (size_t)KC*4 + 4, _mm256_sub_epi64(LDU(r0 + ld + g*8+4), LDU(r1 + ld + g*8+4)));
          }
        }
      }
    }
    return;
  }
  /* 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 * 32)
static long long wg_buf_store[9 * (512 * 512 + 3072 + 48)] __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(1 << 21)));
__attribute__((constructor)) static void sm_thp_mid_init(void) {
  uintptr_t a = (uintptr_t)wg_mid, e = a + (uintptr_t)sizeof(wg_mid);
  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 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+32,*M3=M2+sz+3072+32,*M4=M3+sz+3072+32;
  long long *M5=M4+sz+3072+32,*M6=M5+sz+3072+32,*M7=M6+sz+3072+32,*P=M7+sz+3072+32,*Q=P+sz+3072+32;
  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

评测结果