代码

/*
References:
- Duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/104048: copied the public matrix multiplication implementation as the base. Its public source contains no license notice.
- Duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/103841: the inherited MC/NC sizing originated in our own submission.
- Duck.ac user saffah_cc_v41_260924, https://duck.ac/submission/102254: the inherited AVX2 microkernel originated here; its public source contains no license notice.
- Brice Boyer et al., "Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm", https://hal.science/hal-00163141v3/file/winoschedule.pdf : used the seven-product equations and output recombination as mathematical ideas; no text or code copied.
Approach:
- Materialize Winograd operands at the 1024 level, then reuse the inherited Strassen leaf structure for each 512 product. This changes the 256 leaves from up to four source terms to at most two, and merges the seven products into C in one vector pass.
Purpose:
- Experimental scored comparison of operand packing cost versus Winograd materialization and output combination.
*/
// ===== REFERENCES =====
// [1] 本账号（saffah_cc_v41_agg1）提交 #103939 <https://duck.ac/submission/103939>
//     （122.819 ms，mem 11328 KB）—— **本文件的引擎正文逐字节来自 [1]**。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #103841 <https://duck.ac/submission/103841>
//     —— [1] 的面板尺寸（MC 48->24、NC 512->256）来自 [2]（[2] 本身是在本账号 #102589
//     的正文上只改这两个 #define 得到的，见 [1] 的 REFERENCES/思路 取证）。
// [3] duck.ac 用户 saffah_cc_v41_260924，提交 #102254 <https://duck.ac/submission/102254>
//     —— mikro() 汇编体与 LDQ/STQ 宏继承自它（见 [1]）。
// ======================
// ===== 思路 =====
// 现状：mine = 122.818854（#103939），T = 122756.363（[2]），严支线 121.529799 ⇒ 还需 1.289 ms。
//
// ★ 本发**只改一处**：`mmxR` 的最后一级（`h == CUTL`）用**双目的地 leaf** 取代
//   「leaf 写一个象限 + cmb 搬到另一个象限」/「写 scratch + 两次 cmb」。
//
//   原 `MK` 宏里，一个乘积 M(qi,qj) 要落到目的地块 d 的两个象限 Q[qi]、Q[qj]：
//     oa == 0 -> `mmxR(h,...,Q[qi],...,0,T)` 再 `cmb(Q[qj], Q[qi], ...)`
//     else    -> `mmxR(h,...,T,...)` 再两次 cmb
//   ★ 而 `leaf(n, Ao, Bo, C1, ldc1, op1, C2, ldc2, op2)` **本来就支持两个目的地**
//     （1 层引擎的 M1..M5 就是这么用的，[3] 的原 `mmx`）。当 `h <= CUTL` 时子调用
//     已经是 leaf，把两个象限一次性交给它即可：**一次 mikro、零 scratch、零 cmb**。
//   （`h > CUTL` 时不能这么做：那时子调用还要再往下递归，一个乘积的两条贡献 × 两个
//     父目的地 = 4 个目的地，超出 leaf 的 2 个。本题 CUTL=256、顶层 n=1024 ⇒ 只有
//     第二级 h=512 的 7 个 mmxR(512) 调用享受本改动。）
//
// ★ 量化（都在**判题机**上、同一二进制内运行期开关轮转，min 各一轮）：
//   `cmb` 调用数 **62 -> 7**、cmb 搬运量 **78.5 MB -> 26.0 MB**（本地 instrument 版普查）；
//   整机 Δ = **-0.24% / -0.15%**（两轮独立 base/dd 交替，方向一致，2/2 为负）。
//   ★ 这是一条**带数字的负结果**：`skip cmb` 消融值 3.9%，但把 cmb 的活搬进 mikro epilogue
//   只收回 0.2% ⇒ **cmb 的价主要在"目的地的 RMW 本身"，不在 scratch 流量**（§2.18.771 的又一例）。
//
// 未动：mikro 内联汇编体、packA/packB 的全部路径、MC/NC/KC/CUTL、eb 式 epilogue、面板布局。
// 正确性：对拍 #103939 引擎 n=128/256/512/1024 逐元素 int64 **bad=0**（work/ml7_vfy_dd）。
// 本次为正式提交：把上一发的 0.2% 落袋（中间改进，无风险）。
// ================
/* 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>

/* 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");
}

/* [1] #101958 的 mikro 函数体里用它自己的 32 字节存取宏，逐字节带回这里 */
#define LDQ(p) ({ __m256i _ldq_v; __asm__("vmovdqu %1,%0" : "=x"(_ldq_v) : "m"(*(const __m256i *)(p)) : "memory"); _ldq_v; })
#define STQ(p, v) __asm__ volatile("vmovdqu %0,%1" :: "x"(v), "m"(*(__m256i *)(p)) : "memory")
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC push_options
#pragma GCC target("avx2")

typedef long long TE;

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

static TE Apanel[(size_t)MC * KC * 1 + 64] __attribute__((aligned(64)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(64)));
static __m256i AC[12] __attribute__((aligned(64)));

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

struct Opd {                 /* value(i,j) = sum_t s[t]*p[(r[t]+i)*ld + c[t]+j]  (nt <= 4) */
  const long long *p; int ld; int nt; int r[4], c[4], s[4];
};
static inline Opd oq(const Opd *o, int qi, int qj, int h) {
  Opd d = *o;
  for (int t = 0; t < o->nt; t++) { d.r[t] += qi * h; d.c[t] += qj * h; }
  return d;
}
static inline Opd oadd(const Opd *a, const Opd *b, int sgn) {
  Opd d = *a;
  for (int t = 0; t < b->nt; t++) { d.r[d.nt] = b->r[t]; d.c[d.nt] = b->c[t]; d.s[d.nt] = b->s[t] < 0 ? -sgn : sgn; d.nt++; }
  return d;
}
static inline Opd o1(const long long *p, int ld) {
  Opd d; d.p = p; d.ld = ld; d.nt = 1; d.r[0] = 0; d.c[0] = 0; d.s[0] = 1; return d;
}

/* A panel: rows [ic, ic+mc), k in [pc, pc+kc) of the operand */

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

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

    jb0 = ngrp * 8;
    for (int k = kc; k < KC; k++) {
      long long *qb = Bpanel + (size_t)k * 4;
      for (int g = 0; g < ngrp; g++)
        { long long *q0 = qb + (size_t)g * 8 * KC; STU(q0, _mm256_setzero_si256()); STU(q0 + (size_t)KC * 4, _mm256_setzero_si256()); }
    }
  }
  if (nt >= 3) jb0 = NC;
  for (int jb = jb0; jb < NC; jb += 4) {
    long long *q = Bpanel + (size_t)(jb / 4) * KC * 4;
    int full = (jb + 4 <= nc);
    int lim = nc - jb; if (lim > 4) lim = 4; if (lim < 0) lim = 0;
    const long long *sq = P + (size_t)o->r[0] * ld + o->c[0] + jc + jb;
    const long long *sq1 = (nt >= 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc + jb : 0;
    if (full && nt <= 2) {
      for (int k = 0; k < kc; k++) {
        __m256i v = LDU(sq + (size_t)(pc + k) * ld);
        if (nt == 2) {
          __m256i w = LDU(sq1 + (size_t)(pc + k) * ld);
          v = (o->s[1] > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
        }
        STU(q + (size_t)k * 4, v);
      }
      for (int k = kc; k < KC; k++) STU(q + (size_t)k * 4, _mm256_setzero_si256());
    } else {
      for (int k = 0; k < KC; k++) {
        long long t[4] = {0, 0, 0, 0};
        if (k < kc) {
          int ncol = full ? 4 : lim;
          for (int c = 0; c < ncol; c++) {
            long long v = 0;
            for (int tt = 0; tt < nt; tt++)
              v += o->s[tt] < 0 ? -P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c]
                                :  P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c];
            t[c] = v;
          }
        }
        __m256i v = LDU(t);
        STU(q + (size_t)k * 4, v);
      }
    }
  }
}


static inline int effop(int base, int term) {   /* term: 1 = +M, 2 = -M ; result: 0 store 1 add 2 sub */
  if (base == 0) return term == 1 ? 0 : 2;
  if (base == 1) return term == 2 ? 2 : 1;
  return term == 2 ? 1 : 2;
}

/* C1 op1 (+- second destination C2 op2, either may be absent) = Ao * Bo, block size n */
static void leaf(int n, const Opd *Ao, const Opd *Bo,
                 long long *C1, int ldc1, int op1, long long *C2, int ldc2, int op2) {
  int nb4 = n - (n & 3);
  for (int jc = 0; jc < nb4; jc += NC) {
    int nc = (nb4 - jc < NC) ? (nb4 - jc) : NC;
    for (int pc = 0; pc < n; pc += KC) {
      int kc = (n - pc < KC) ? (n - pc) : KC;
      packB(Bo, jc, nc, pc, kc);
      int e1 = (pc == 0) ? op1 : (op1 == 2 ? 2 : 1);
      int e2 = (pc == 0) ? op2 : (op2 == 2 ? 2 : 1);
      for (int ic = 0; ic < n; ic += MC) {
        int mc = (n - ic < MC) ? (n - ic) : MC;
        packA(Ao, ic, mc, pc, kc);
        for (int ir = 0; ir < mc; ir += MR) {
          int rows = (mc - ir < MR) ? (mc - ir) : MR;
          const long long *ap = Apanel + (size_t)ir * KC * 1;
          for (int jg = 0; jg < nc; jg += 4) {
            const long long *bp = Bpanel + (size_t)(jg / 4) * KC * 4;
            mikro(KC / UU, ap, bp, C1 + (size_t)(ic + ir) * ldc1 + jc + jg, ldc1, e1,
                  C2 ? C2 + (size_t)(ic + ir) * ldc2 + jc + jg : 0, ldc2, e2, rows);
          }
        }
      }
    }
  }
  for (int j = nb4; j < n; j++)
    for (int i = 0; i < n; i++) {
      unsigned long long acc = 0;
      for (int k = 0; k < n; k++) {
        long long av = 0, bv = 0;
        for (int t = 0; t < Ao->nt; t++)
          av += Ao->s[t] < 0 ? -Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k]
                             :  Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k];
        for (int t = 0; t < Bo->nt; t++)
          bv += Bo->s[t] < 0 ? -Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j]
                             :  Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j];
        acc += (unsigned long long)av * (unsigned long long)bv;
      }
      long long *cp = C1 + (size_t)i * ldc1 + j;
      if (op1 == 0) *cp = (long long)acc;
      else if (op1 == 1) *cp += (long long)acc;
      else *cp -= (long long)acc;
      if (C2) {
        long long *cq = C2 + (size_t)i * ldc2 + j;
        if (op2 == 0) *cq = (long long)acc;
        else if (op2 == 1) *cq += (long long)acc;
        else *cq -= (long long)acc;
      }
    }
}

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

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

/* dst op= sign * T , (h x ldt) -> (h x ldd) , sign = +-1 , op in {0,1,2} */
static void cmb(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op) {
  for (int i = 0; i < h; i++) {
    const long long *t = T + (size_t)i * ldt;
    long long *c = dst + (size_t)i * ldd;
    int j = 0;
    if (op == 0) {
      for (; j + 4 <= h; j += 4) {
        __m256i v = LDQ((const __m256i *)(t + j));
        if (sign < 0) v = _mm256_sub_epi64(_mm256_setzero_si256(), v);
        STQ((__m256i *)(c + j), v);
      }
      for (; j < h; j++) c[j] = sign < 0 ? -t[j] : t[j];
    } else {
      for (; j + 4 <= h; j += 4) {
        __m256i v = LDQ((const __m256i *)(t + j));
        __m256i w = LDQ((const __m256i *)(c + j));
        if ((op == 1) != (sign < 0)) STQ((__m256i *)(c + j), _mm256_add_epi64(w, v));
        else                         STQ((__m256i *)(c + j), _mm256_sub_epi64(w, v));
      }
      for (; j < h; j++) { if ((op == 1) != (sign < 0)) c[j] += t[j]; else c[j] -= t[j]; }
    }
  }
}

/* single-destination recursive Strassen: d op= Ao*Bo , T = scratch block (>= 2*h*h) */
static void mmxR(int n, const Opd *Ao, const Opd *Bo, long long *d, int ld, int op, long long *T) {
  if (n <= CUTL || (n & 1) || Ao->nt > 2 || Bo->nt > 2) {
    leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); return; }
  const int h = n >> 1;
  Opd A11 = oq(Ao, 0, 0, h), A12 = oq(Ao, 0, 1, h), A21 = oq(Ao, 1, 0, h), A22 = oq(Ao, 1, 1, h);
  Opd B11 = oq(Bo, 0, 0, h), B12 = oq(Bo, 0, 1, h), B21 = oq(Bo, 1, 0, h), B22 = oq(Bo, 1, 1, h);
  long long *Q[4];
  Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
  long long *T2 = T + (size_t)h * h;
  int wr[4] = {0, 0, 0, 0};
#define MK(qi, si, qj, sj, Ma, Mb) do { \
    int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
    int ob = wr[qj] ? effop(1, (sj)) : effop(op, (sj)); \
    wr[qi] = 1; wr[qj] = 1; \
    if (h <= CUTL) { \
      /* sub-call IS a leaf: hand it BOTH destination quadrants directly \
         (one mikro pass, no scratch, no cmb) */ \
      leaf(h, &(Ma), &(Mb), Q[qi], ld, oa, Q[qj], ld, ob); \
    } else if (oa == 0) { \
      mmxR(h, &(Ma), &(Mb), Q[qi], ld, 0, T); \
      cmb(Q[qj], ld, Q[qi], ld, h, (sj), ob); \
    } else { \
      mmxR(h, &(Ma), &(Mb), T, h, 0, T2); \
      cmb(Q[qi], ld, T, h, h, (si), oa); \
      cmb(Q[qj], ld, T, h, h, (sj), ob); \
    } \
  } while (0)
#define MK1(qi, si, Ma, Mb) do { \
    int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
    wr[qi] = 1; \
    mmxR(h, &(Ma), &(Mb), Q[qi], ld, oa, T); \
  } while (0)
  { /* M1 = (A11+A22)(B11+B22) -> +C11, +C22 */
    Opd Ma = oadd(&A11, &A22, 1), Mb = oadd(&B11, &B22, 1);
    MK(0, 1, 3, 1, Ma, Mb);
  }
  { /* M2 = (A21+A22)B11 -> +C21, -C22 */
    Opd Ma = oadd(&A21, &A22, 1), Mb = B11;
    MK(2, 1, 3, 2, Ma, Mb);
  }
  { /* M3 = A11(B12-B22) -> +C12, +C22 */
    Opd Ma = A11, Mb = oadd(&B12, &B22, -1);
    MK(1, 1, 3, 1, Ma, Mb);
  }
  { /* M4 = A22(B21-B11) -> +C11, +C21 */
    Opd Ma = A22, Mb = oadd(&B21, &B11, -1);
    MK(0, 1, 2, 1, Ma, Mb);
  }
  { /* M5 = (A11+A12)B22 -> -C11, +C12 */
    Opd Ma = oadd(&A11, &A12, 1), Mb = B22;
    MK(0, 2, 1, 1, Ma, Mb);
  }
  { /* M6 = (A21-A11)(B11+B12) -> +C22 */
    Opd Ma = oadd(&A21, &A11, -1), Mb = oadd(&B11, &B12, 1);
    MK1(3, 1, Ma, Mb);
  }
  { /* M7 = (A12-A22)(B21+B22) -> +C11 */
    Opd Ma = oadd(&A12, &A22, -1), Mb = oadd(&B21, &B22, 1);
    MK1(0, 1, Ma, Mb);
  }
#undef MK
#undef MK1
}


static long long wg_buf[9 * 512 * 512] __attribute__((aligned(64)));

static void wg_add(int h, long long *dst, const long long *a, int lda,
                   const long long *b, int ldb, int sg) {
  for (int i = 0; i < h; i++) {
    long long *d = dst + (size_t)i*h;
    const long long *x = a + (size_t)i*lda, *y = b + (size_t)i*ldb;
    for (int j = 0; j < h; j += 4) {
      __m256i vx=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;
    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);
      STU(c11+j,_mm256_add_epi64(p1,p2));
      STU(c12+j,_mm256_add_epi64(u4,p3));
      STU(c21+j,_mm256_sub_epi64(u3,p4));
      STU(c22+j,_mm256_add_epi64(u3,p5));
    }
  }
}
static void winograd_top(int n, const long long *A, const long long *B, long long *C) {
  const int h=n/2; const size_t sz=(size_t)h*h;
  long long *M1=wg_buf,*M2=M1+sz,*M3=M2+sz,*M4=M3+sz;
  long long *M5=M4+sz,*M6=M5+sz,*M7=M6+sz,*P=M7+sz,*Q=P+sz;
  const long long *A11=A,*A12=A+h,*A21=A+(size_t)h*n,*A22=A21+h;
  const long long *B11=B,*B12=B+h,*B21=B+(size_t)h*n,*B22=B21+h;
#define PROD(X,Y,Z,LX,LY) do { Opd ao=o1((X),(LX)),bo=o1((Y),(LY)); mmxR(h,&ao,&bo,(Z),h,0,g_scratch); } while(0)
  PROD(A11,B11,M1,n,n);
  PROD(A12,B21,M2,n,n);
  wg_add(h,P,A21,n,A22,n,+1);             // S1
  wg_add(h,Q,B12,n,B11,n,-1);             // T1
  PROD(P,Q,M5,h,h);
  wg_mutate(h,P,A11,n,0);                 // S2 = S1 - A11
  wg_mutate(h,Q,B22,n,1);                 // T2 = B22 - T1
  PROD(P,Q,M6,h,h);
  wg_mutate(h,P,A12,n,1);                 // S4 = A12 - S2
  PROD(P,B22,M3,h,n);
  wg_mutate(h,Q,B21,n,0);                 // T4 = T2 - B21
  PROD(A22,Q,M4,n,h);
  wg_add(h,P,A11,n,A21,n,-1);             // S3
  wg_add(h,Q,B22,n,B12,n,-1);             // T3
  PROD(P,Q,M7,h,h);
#undef PROD
  wg_combine(h,C,n,M1,M2,M3,M4,M5,M6,M7);
}
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; }

评测结果