代码

// References:
// [1] Our account saffah_codex_6s_agg2, https://duck.ac/submission/106048 :
//     accepted float-1k Strassen/AVX2 source copied with all inherited
//     author and source citations retained below.
// [2] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/105354 :
//     the inherited int16 Strassen engine; no separate license notice was
//     displayed with the public source.
// Approach:
// Place packed B at 512 int16 values from its aligned backing array
// instead of 1536. Packed panel data and arithmetic remain identical;
// the cache-set relationship with packed A, C, and products changes.
// Purpose:
// Official judged cache-phase search to clear the exact mmmf1k target.
// References: saffah_cc_v41_agg1, https://duck.ac/submission/105354 : copied full int16 Strassen engine; no license notice was shown.
// Approach: Use non-temporal stores while quantizing inputs into the output buffer, and fence before the matrix products.
// Purpose: Measure official mmmf1k speed relative to the cited baseline.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #105282
//     <https://duck.ac/submission/105282>（文件 work/fs_qainC.cpp）
//     用途：本文件的**正文基座**——该件本身 = 本账号 #104948 之上两处改动：(a) `pack_b` 换成
//     对手 #105191 的 tile-4 循环，(b) `qa`/`qb` 搬进 `C` 的缓冲区（`qa=(int16_t*)C;
//     qb=qa+N*N`）。本件**只在其上做第三节那三处改动**，其余一字未动。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #104180 <https://duck.ac/submission/104180>
//     用途：`mikro6` 函数体（该提交的 6x16 int16 `vpmaddwd` 内联汇编）——本件一个字节未改。
// [3] 本账号提交 #105295 <https://duck.ac/submission/105295>（`work/fs1_sub4.cpp`）——把
//     M2/M4 的累加面搬出 `prod[]` 的第一版（在 #104948 基座上）。
// [4] 本账号提交 #104948 <https://duck.ac/submission/104948>、#104302、#103973、#103927、
//     #103535、#103476、#103032、#102993 —— int16 `vpmaddwd` + 一层 Strassen 引擎谱系。
// ======================
// ===== 思路 =====
// 正式提交。在基座 #105282 上做**三处互锁的改动**（整体是一次改动：把两个乘积累加面从
// `prod[]` 搬进 `C` 自己缓冲区内已经死掉的象限），目的是**再省 2 MB 首触页代价**。
//
// 一、已确立的机理（本代理实测，见 notes.md）
//   * 判题机正式成绩里有 4.7%（2.47 Mcyc）是**首触页代价**（同进程三连调用探针 t1=51.99M /
//     t2=49.53M ⇒ c12=2.47M）⇒ **每 1 MB 被触碰的 scratch ≈ 0.204 Mcyc ≈ 57 µs**。
//   * 基座 #105282 已经把 `qa`/`qb`（4 MB）搬进了 `C`；剩下 `prod[7]`（7 MB）是最大的可删项。
//   * ⚠ 反面教材（本代理 #105214 = 16.689 ms ✗✗）：把累加面搬进 `C` 的**象限**（4096 B 行距）会
//     让 `mikro6` 每次的 6 条 NT 存数落在 6 个不同 4 KB 页上 ⇒ mikro 的 warm 相位 49.5→56.9 Mcyc。
//     ⇒ **新家必须与 `prod` 同形：512 行 × 512 int32、行距 2048 B。**
//
// 二、三处改动
//   (1) 次序：`M1 M3 M5 M6 M7 | M2 M4`。
//       理由：`M2` 的累加面要住进 `qa` 的上半（= `C` 的 [0,1 MB)，即 a11/a12 那片），
//             而 a11 的最后读者是 M6、a12 是 M7；`M4` 的累加面住进 `qb` 的上半
//             （= `C` 的 [2 MB,3 MB)，即 b11/b12 那片），b12 的最后读者是 M6、b11 的最后
//             读者是**它自己**（M4 的 pack_b 先读、mikro 后写）。
//       ⇒ 五个 product 跑完后那两片已成死区；而 M2 只读 a21/a22 与 b11、M4 只读 a22 与 b21/b11，
//         **都不碰自己的新家**。
//   (2) `M2` 的输出基址 → `(int32_t *)qa`；`M4` 的输出基址 → `(int32_t *)qb`。
//       两者的行距与 `prod` 完全一致（`qa`/`qb` 的行距就是 `N*2 = 2048 B`）⇒ `mikro6` 的
//       12 条 NT 存数地址形态**与基座逐位相同**，只需换基址寄存器。
//       基址必须是 2048 B 的整数倍（见下），所以取 offset 0，而不是 +1024 B。
//   (3) `prod[7]` → `prod[5]`；`combine()` 的 m2/m4 改为从 `qa`/`qb` 读，并且**i 循环倒序**。
//       ★ 倒序的理由：累加面住在 `C` 里，而 `combine` 的输出会就地覆写它们。输出第 i 行的
//         32 B 存数正好落在累加面的第 2i / 2i+1 行（同一列组）；那两行的**读取发生在循环下标
//         2i / 2i+1 > i**。i 递减 ⇒ 2i、2i+1 这两行在前面的迭代里已经读完了；若 i 递增，则会
//         在读取之前被破坏（本代理实测：递增时 maxdiff=1077，倒序后逐位正确）。
//
// 三、正确性：与基座**逐位相同**（work/f1_drv.cpp 三个种子 maxdiff 完全一致 =
//   0.035762@(621,491) / 0.035645@(724,274) / 0.040654@(928,499)；c11/c12/c21/c22 四象限
//   maxdiff 均 0.034~0.036）。
//
// 四、预期：被触碰的 scratch 7.7 MB → 5.7 MB（−2 MB ≈ −0.41 Mcyc ≈ −114 µs）。
//   参照上一版在旧基座上的同类改动 #105295：14.666532 → 14.557080（−109.5 µs），与账目吻合。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #105191 <https://duck.ac/submission/105191>
//     用途：直接复制了该提交 `pack_b` 的循环结构（`p` 按 4 个 k-pair 分块 → `jb` 扫描 16 列面板 → 块内 4 个 p）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104948 <https://duck.ac/submission/104948>
//     用途：本文件的正文基座（本账号最好件 14.666532 ms）。
// [3] 兄弟 lane `fs1_` 的判题机实测：**首触页代价 ≈ 0.204 Mcyc/MB ≈ 57 µs/MB**（探针 `.p/fs1_three.cpp`，
//     同进程三连调用、第 2/3 次不再缺页）。本发由该结论驱动。
// [4] 本账号 #104302 / #103973 / #103032（quantize / pack_a / mikro6 / combine 引擎谱系）。
// ======================
// ===== 思路 =====
// 正式提交（两处单变量，均在已验证基座 `#104948` 之上，逐条独立可解释）：
//   (a) `pack_b` 换成对手 `#105191` 的 **tile-4** 循环结构（`pb0` 步 4 外层 / `jb` 中层 / 块内 4 个 `p`）。
//       机理：源步长 `2*N*2 = 4096 B` 是 4 KB 的整数倍 ⇒ L1（32 KB/8 路/64 组）里 `4096>>6=64 ≡ 0 (mod 64)`
//       ⇒ 落同一组。tile-4 让固定 `jb` 时**同时活跃的源行只有 8 行 = 8 条 line = 恰好 8 路**（不再互相驱逐），
//       而写侧在组内仍是 4 个连续 64 B 块（顺序写）。对手用它从 14.6665 走到 14.5725（−0.64%）。
//   (b) ★ **把 quant 的输出 `qa`/`qb` 搬进 `C` 的缓冲区** ⇒ 完全不触碰 `qa_`/`qb_` 那 4 MB（1024 页）的 .bss。
//       依据：判题机实测**首触页 ≈ 0.204 Mcyc/MB（57 µs/MB）** ⇒ 4 MB ≈ **0.82 Mcyc ≈ 0.24 ms ≈ 1.6%**。
//       安全性（生命周期不重叠，可逐条复核）：
//         * `qa`/`qb` 的**最后使用者**是最后一个 `product()` 里的 `pack_a`/`pack_b`；
//         * `C` **只**在最后一步 `combine(C)` 里被写（NT 存，整块覆写 4 MB）；
//         * 题目保证输出数组初值为 0 ⇒ 在 `combine` 之前没有任何代码读 `C` 的旧内容 ⇒ 拿它当纯 scratch 安全。
//       容量：`qa` 需 `N*N` 个 int16 = 2 MB，`qb` 紧接其后 2 MB ⇒ 合计 4 MB = `n*n*sizeof(float)`（n=1024）**正好放下**；
//       `C` 由题目保证 4096 字节对齐 ⇒ 两个 scratch 各自的 4096 对齐与原来是同规格。
//       副作用：`.bss` 少 4 MB（触碰页 16080 KB → 约 12000 KB）。
// 正确性：`qa`/`qb` 的内容、`pack_*` 的读写地址表达式、`mikro`/`combine` 全部未变；
//   唯一变化是这两个缓冲区的**所在地**（.bss → C 缓冲的前 4 MB），而它们在 `combine` 写 C 之前就已作废。
//   输出 `C` 由 `combine` 的 NT 存整块覆写 ⇒ **逐位不变**。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #105191 <https://duck.ac/submission/105191>
//     用途：**直接复制了该提交 `pack_b` 的循环结构**（`p` 按 4 个 k-pair 分块 → `jb` 扫 16 列面板 → 块内 4 个 p）。
//     该提交 = 本账号 #104948 的正文 + 这一处 `pack_b` 改写（其自身头部即声明引用本账号 #104948）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104948 <https://duck.ac/submission/104948>
//     用途：本文件的正文基座（= 本账号最好件 14.666532 ms；其 `pack_b` 是 p 外/jb 内的换序形态）。
// [3] 本账号 #104302 / #103973 / #103927 / #103032 引擎谱系（quantize / pack_a / mikro / mikro6 / combine / MC=192 / prod 行距）。
// ======================
// ===== 思路 =====
// 正式提交（单变量：只改 `pack_b` 一个函数的循环结构，地址表达式与算术后等价）。
// 原（#104948）: for p (0..255) { for jb (32 组) { 读 2 行 + jb 列, 写 dblk + p*32 } }
// 本发        : for pb0 += 4 { for jb { for p in [pb0,pb0+4) { 同上 } } }
//
// 机理（4 KB 同余 / L1 同组别名 —— BRIEF §2.18.25 + §2.18.323）：
//   源取数地址 = b.x + 2p*N + jb，相邻 p 相差 **2*N*2 = 4096 字节** ⇒ L1（32 KB / 8 路 / 64 组）里
//   `4096>>6 = 64 ≡ 0 (mod 64)` ⇒ **落到同一个 L1 组**。
//   * #104948 的 p 外/jb 内：读侧变成顺序流 ✓，但写侧变成 **32 条相距 16 KB 的流**
//     （16 KB 也是 4 KB 的整数倍）⇒ 32 条写流全部压在同一个 L1 组（8 路）⇒ 写侧 thrash。
//   * 本发的 pb0 分块（4 个 p = 8 个源行）：固定 jb 时同时活跃的源行只有 8 行
//     ⇒ **恰好 8 条 line = L1 的 8 路，正好放得下不再互相驱逐**；而写侧在组内仍是 4 个连续 64 B 块 ✓ 顺序。
//   ⇒ 读侧与写侧同时不再是同组 thrash。对照读数：#104948 = 14.666532（我方），#105191 = 14.572546（对手同结构）。
// 正确性：每个 (p, jb) 的**源地址与目的地址、以及写入的 32 字节内容逐位不变**（只是执行先后次序变了）⇒
//   输出与 #104948 逐字节相同（pack_b 的输出是 bp 面板，随后被 mikro 只读消费）。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104180 <https://duck.ac/submission/104180>
//     用途：直接复制了该提交的 `mikro6` 函数体（6x16 int16 `vpmaddwd` 微内核的内联汇编 k-pair
//     循环，含首对 k-pair 剥离与"两条 madd 写两个不同临时寄存器"的发射形态）。该件本身 =
//     本账号 #103973 的 mikro6 + 这两处改动。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104302
//     <https://duck.ac/submission/104302>
//     用途：本文件的正文基座（= [1] 的两处改动原样搬回后的最好件）。除下面的 `pack_b`
//     循环次序外，本文件与 #104302 逐字节相同。
// [3] 本账号提交 #103973 <https://duck.ac/submission/103973>、#103927
//     <https://duck.ac/submission/103927>、#103535、#103476、#103032、#102993
//     —— int16 `vpmaddwd` + 一层 Strassen(n=1024→512) 引擎谱系（quantize/pack_a/pack_b/
//     mikro/mikro6/combine）。
// ======================
// ===== 思路 =====
// 正式提交（单变量改动，非试验性）。**只改 `pack_b` 一个函数的循环次序**：
//   原： for (int jb=0; jb<H; jb+=NR) { dblk = ...; for (int p=0; p<NP; p++) { x0 = b.x+2p*N+jb; ... } }
//   新： for (int p=0; p<NP; p++) { x0 = b.x+2p*N; ...; for (int jb=0; jb<H; jb+=NR) { dblk = ...; ... x0+jb ... } }
// **地址表达式与写入内容一字不改，只改 (jb,p) 对的访问先后** ⇒ 输出逐位相同。
//
// 机理：原次序下固定一个 jb（16 列）时，对 p 的 256 次源取数地址步长是 2*N*2 = 4096 字节；
// L1 = 32 KB / 64 组 / 64 B 行 ⇒ 4 KB 步长恰好把 512 行**全部映射到同一个 L1 组**（8 路），
// 于是每个取数都是 L1 miss 且互相驱逐，同一个 2 MB 源被重复扫 32 遍。换序后固定 p 时按 jb
// 扫过第 2p / 2p+1 两整行（各 2048 B）⇒ 变成顺序流，写侧变成 32 条相距 16 KB 的流。
// （这是 BRIEF §2.18.22 "packb 循环序重排" 在 mmms1k 上实测 −1.10% 的同一把刀；本件的
//  pack_b 一直是旧次序。）
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104180 <https://duck.ac/submission/104180>
//     用途：**直接复制了该提交 `mikro6` 函数的函数体**（6x16 int16 `vpmaddwd` 微内核的内联汇编
//     k-pair 循环）。该件本身就是**本账号 #103973 的 `mikro6` + 两处改动**（见"思路"一），
//     本发把这两处**原样搬回**，其余一字未动。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #103973
//     <https://duck.ac/submission/103973>
//     用途：本文件的正文基座（= 本代理把 `mikro6` 从 C intrinsics 改写为手写汇编 + k-pair
//     展开 4 次的那一发）；本文件的其余部分（pack/quantize/combine/MC=192/prod 行距）
//     与 #103973 逐字节相同。
// [3] 本账号提交 #103927 <https://duck.ac/submission/103927>（对手 #103503 的三处）、
//     #103535 / #103476 / #103032（MC=192 / prod 行距 +6144 / int16+一层 Strassen 引擎谱系）。
// ======================
// ===== 思路 =====
// 正式提交（非试验性）。本件 = 本账号 #103973 之上**只改 `mikro6` 一个函数**，
// 两处改动**逐字取自对手当前最好件 #104180**：
//   (a) **k 循环首对 k-pair 剥离**：`int k=NP` → `int k=NP-4`，并在 `.p2align 5` 之前放一整块
//       展开体，12 个累加器**由 `vpmaddwd` 直接初始化**（省掉 12 条 `vpxor` 清零与首轮的
//       12 条 `vpaddd`）；进入循环前补 `add $16,%[a]` / `add $256,%[b]`。
//   (b) **每条 (行, k-pair) 的两条 `vpmaddwd` 改成写两个不同的临时寄存器**
//       （`... %ymm12` 与 `... %ymm13`），并**把两条 `vpaddd` 排到两条 madd 之后**：
//       ```
//       我们的旧形态（每行 4 条）：  madd→ymm12 / vpaddd / madd→ymm12 / vpaddd
//       本件（每行仍是 4 条）：      madd→ymm12 / madd→ymm13 / vpaddd / vpaddd
//       ```
//       **指令数完全相同**，只是**发射形态**变了：两条互相独立的 madd 背靠背发出，
//       第二条 madd 的延迟与前一条 `vpaddd` 重叠 ⇒ 正是 §2.18.691/§2.18.752
//       "**编码/发射形态 > uop 数**"的实例。
//   ⚠ 注意 (b) 让 `ymm13`（广播寄存器）在最后一次使用后被复用为临时寄存器；
//     16 个 ymm 恰好用满（12 累加器 + ymm12 + ymm13 + ymm14/ymm15 = 16），无溢出。
//
// 一、判题机定价（同二进制 6 臂 = 本件 ×3 + #103973 正文 ×3、3 趟轮转、min-of-N、FNV 闸门）
//   组 0/1/2 = #103973 正文：min = 50867 / 50847 / **50789**
//   组 3/4/5 = 本件      ：min = **50238** / 50208 / 50211
//   ⇒ **组间 −1.14%**，而**组内（同码三份）散布只有 0.15%** ⇒ 效应 ≫ 噪声 ✓
//   （★ 仪器学：本台架的"namespace 位置偏置"在只放一对同码臂时会达到 1.12%
//     —— 单对地板不够，**必须放 3 份取组最小**，这一条已写进 notes。）
//   预期正式提交：14.939893 × (1−0.0114) ≈ **14.77 ms**。
//
// 二、正确性
//   与 #103973 **逐字节相同**（同进程 3 份本件 + 3 份正文，6 臂 FNV 哈希全部相等，`GATE bad=0`）。
// ================
#pragma GCC optimize("O3,unroll-loops")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,fma")
#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>

#define N 1024
#define H 512
#define KC H
#define MC 192
#define NR 16
#define MR 4
#define NP (KC / 2)

/* qa/qb 不再占 .bss：它们在 matrix_multiply 里指向 C 的缓冲区（见【思路】）。
   C 只在最后一步 combine() 才被写，qa/qb 的最后使用者在最后一个 product() 里 ⇒ 生命周期不重叠。 */
static int16_t *qa, *qb;
alignas(4096) static int16_t ap[MC*KC], bp_storage[H*KC+2048];
static int16_t *const bp = bp_storage + 512;
alignas(4096) static int32_t prod[5][H*H+6144];

static inline __m256i ldu(const void *p) {
  __m256i v;
  __asm__("vmovdqu %1,%0" : "=x"(v) : "m"(*(const __m256i *)p));
  return v;
}
static inline __m256i bc32(const void *p) {
  __m256i v;
  __asm__("vpbroadcastd %1,%0" : "=x"(v) : "m"(*(const int *)p));
  return v;
}

static void quantize(const float *src, int16_t *dst) {
  const __m256 scale = _mm256_set1_ps(1024.0f);
  for (int i = 0; i < N*N; i += 16) {
    __m256i x = _mm256_cvtps_epi32(_mm256_mul_ps(_mm256_load_ps(src+i), scale));
    __m256i y = _mm256_cvtps_epi32(_mm256_mul_ps(_mm256_load_ps(src+i+8), scale));
    _mm256_stream_si256((__m256i *)(dst+i),
      _mm256_permute4x64_epi64(_mm256_packs_epi32(x,y),0xD8));
  }
}

struct Mat { const int16_t *x; const int16_t *y; int sign; };
static inline __m256i vec(const int16_t *x, const int16_t *y, int sign) {
  __m256i v = _mm256_load_si256((const __m256i *)x);
  if (y) {
    __m256i w = _mm256_load_si256((const __m256i *)y);
    v = sign > 0 ? _mm256_add_epi16(v,w) : _mm256_sub_epi16(v,w);
  }
  return v;
}

static void pack_a(Mat a, int ic, int mc) {
  for (int r = 0; r < mc; r++) {
    const int16_t *x = a.x + (size_t)(ic+r)*N;
    const int16_t *y = a.y ? a.y + (size_t)(ic+r)*N : nullptr;
    int16_t *d = ap + (size_t)r*KC;
    for (int k = 0; k < KC; k += 16)
      _mm256_store_si256((__m256i *)(d+k), vec(x+k,y?y+k:nullptr,a.sign));
  }
}

static void pack_b(Mat b) {
  for (int pb0=0; pb0<NP; pb0+=4) {
    for (int jb=0; jb<H; jb+=NR) {
      int16_t *dblk = bp + (size_t)(jb/NR)*KC*NR;
      for (int p=pb0; p<pb0+4; ++p) {
        const int16_t *x0 = b.x + (size_t)(2*p)*N + jb;
        const int16_t *x1 = x0 + N;
        const int16_t *y0 = b.y ? b.y + (size_t)(2*p)*N + jb : nullptr;
        const int16_t *y1 = y0 ? y0 + N : nullptr;
        __m256i x = vec(x0,y0,b.sign), y = vec(x1,y1,b.sign);
        __m256i lo = _mm256_unpacklo_epi16(x,y), hi = _mm256_unpackhi_epi16(x,y);
        int16_t *d = dblk + (size_t)p*NR*2;
        _mm256_store_si256((__m256i *)d, _mm256_permute2x128_si256(lo,hi,0x20));
        _mm256_store_si256((__m256i *)(d+16), _mm256_permute2x128_si256(lo,hi,0x31));
      }
    }
  }
}

static inline void mikro(const int16_t *a, const int16_t *b, int32_t *c) {
  __m256i c00,c01,c10,c11,c20,c21,c30,c31;
  c00=c01=c10=c11=c20=c21=c30=c31=_mm256_setzero_si256();
  const char *p0=(const char *)a, *p1=(const char *)(a+KC),
             *p2=(const char *)(a+2*KC), *p3=(const char *)(a+3*KC);
  const char *pb=(const char *)b;
  for (int p=0; p<NP; p++) {
    __m256i b0=ldu(pb), b1=ldu(pb+32), av;
    av=bc32(p0); c00=_mm256_add_epi32(c00,_mm256_madd_epi16(av,b0));
                   c01=_mm256_add_epi32(c01,_mm256_madd_epi16(av,b1));
    av=bc32(p1); c10=_mm256_add_epi32(c10,_mm256_madd_epi16(av,b0));
                   c11=_mm256_add_epi32(c11,_mm256_madd_epi16(av,b1));
    av=bc32(p2); c20=_mm256_add_epi32(c20,_mm256_madd_epi16(av,b0));
                   c21=_mm256_add_epi32(c21,_mm256_madd_epi16(av,b1));
    av=bc32(p3); c30=_mm256_add_epi32(c30,_mm256_madd_epi16(av,b0));
                   c31=_mm256_add_epi32(c31,_mm256_madd_epi16(av,b1));
    p0+=4;p1+=4;p2+=4;p3+=4;pb+=64;
  }
#define ST(R,X,Y) do {\
    _mm256_stream_si256((__m256i *)(c+(size_t)(R)*H),X);\
    _mm256_stream_si256((__m256i *)(c+(size_t)(R)*H+8),Y);\
  } while(0)
  ST(0,c00,c01); ST(1,c10,c11); ST(2,c20,c21); ST(3,c30,c31);
#undef ST
}


/* 6 rows x 16 cols over all KC k -- same panel layout, 12 accumulators (gcc 9.3 keeps them
   in ymm: 12 acc + 2 b + 1 broadcast = 15 <= 16).  Measured on the judge (work/f1_mikro4.cpp,
   same binary, floor 0.10%): 48.54 vs 49.11 cyc per 1024 MAC for the 4x16 shape = -1.16%. */
static inline void mikro6(const int16_t *a, const int16_t *b, int32_t *c) {
  const char *pa=(const char *)a, *pb2=(const char *)b;
  int k=NP-4;
  __asm__ __volatile__(
    "vmovdqu 0(%[b]), %%ymm14\n\t"
    "vmovdqu 32(%[b]), %%ymm15\n\t"
    "vpbroadcastd 0(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm0\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm1\n\t"
    "vpbroadcastd 1024(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm2\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm3\n\t"
    "vpbroadcastd 2048(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm4\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm5\n\t"
    "vpbroadcastd 3072(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm6\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm7\n\t"
    "vpbroadcastd 4096(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm8\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm9\n\t"
    "vpbroadcastd 5120(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm10\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm11\n\t"
    "vmovdqu 64(%[b]), %%ymm14\n\t"
    "vmovdqu 96(%[b]), %%ymm15\n\t"
    "vpbroadcastd 4(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
    "vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
    "vpbroadcastd 1028(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
    "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
    "vpbroadcastd 2052(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
    "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
    "vpbroadcastd 3076(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t"
    "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t"
    "vpbroadcastd 4100(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm8, %%ymm8\n\t"
    "vpaddd %%ymm13, %%ymm9, %%ymm9\n\t"
    "vpbroadcastd 5124(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm10, %%ymm10\n\t"
    "vpaddd %%ymm13, %%ymm11, %%ymm11\n\t"
    "vmovdqu 128(%[b]), %%ymm14\n\t"
    "vmovdqu 160(%[b]), %%ymm15\n\t"
    "vpbroadcastd 8(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
    "vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
    "vpbroadcastd 1032(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
    "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
    "vpbroadcastd 2056(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
    "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
    "vpbroadcastd 3080(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t"
    "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t"
    "vpbroadcastd 4104(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm8, %%ymm8\n\t"
    "vpaddd %%ymm13, %%ymm9, %%ymm9\n\t"
    "vpbroadcastd 5128(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm10, %%ymm10\n\t"
    "vpaddd %%ymm13, %%ymm11, %%ymm11\n\t"
    "vmovdqu 192(%[b]), %%ymm14\n\t"
    "vmovdqu 224(%[b]), %%ymm15\n\t"
    "vpbroadcastd 12(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
    "vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
    "vpbroadcastd 1036(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
    "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
    "vpbroadcastd 2060(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
    "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
    "vpbroadcastd 3084(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t"
    "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t"
    "vpbroadcastd 4108(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm8, %%ymm8\n\t"
    "vpaddd %%ymm13, %%ymm9, %%ymm9\n\t"
    "vpbroadcastd 5132(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm10, %%ymm10\n\t"
    "vpaddd %%ymm13, %%ymm11, %%ymm11\n\t"
    "add $16, %[a]\n\t"
    "add $256, %[b]\n\t"
    ".p2align 5\n\t"
    "1:\n\t"
    "vmovdqu 0(%[b]), %%ymm14\n\t"
    "vmovdqu 32(%[b]), %%ymm15\n\t"
    "vpbroadcastd 0(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
    "vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
    "vpbroadcastd 1024(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
    "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
    "vpbroadcastd 2048(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
    "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
    "vpbroadcastd 3072(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t"
    "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t"
    "vpbroadcastd 4096(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm8, %%ymm8\n\t"
    "vpaddd %%ymm13, %%ymm9, %%ymm9\n\t"
    "vpbroadcastd 5120(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm10, %%ymm10\n\t"
    "vpaddd %%ymm13, %%ymm11, %%ymm11\n\t"
    "vmovdqu 64(%[b]), %%ymm14\n\t"
    "vmovdqu 96(%[b]), %%ymm15\n\t"
    "vpbroadcastd 4(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
    "vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
    "vpbroadcastd 1028(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
    "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
    "vpbroadcastd 2052(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
    "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
    "vpbroadcastd 3076(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t"
    "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t"
    "vpbroadcastd 4100(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm8, %%ymm8\n\t"
    "vpaddd %%ymm13, %%ymm9, %%ymm9\n\t"
    "vpbroadcastd 5124(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm10, %%ymm10\n\t"
    "vpaddd %%ymm13, %%ymm11, %%ymm11\n\t"
    "vmovdqu 128(%[b]), %%ymm14\n\t"
    "vmovdqu 160(%[b]), %%ymm15\n\t"
    "vpbroadcastd 8(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
    "vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
    "vpbroadcastd 1032(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
    "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
    "vpbroadcastd 2056(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
    "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
    "vpbroadcastd 3080(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t"
    "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t"
    "vpbroadcastd 4104(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm8, %%ymm8\n\t"
    "vpaddd %%ymm13, %%ymm9, %%ymm9\n\t"
    "vpbroadcastd 5128(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm10, %%ymm10\n\t"
    "vpaddd %%ymm13, %%ymm11, %%ymm11\n\t"
    "vmovdqu 192(%[b]), %%ymm14\n\t"
    "vmovdqu 224(%[b]), %%ymm15\n\t"
    "vpbroadcastd 12(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
    "vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
    "vpbroadcastd 1036(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
    "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
    "vpbroadcastd 2060(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
    "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
    "vpbroadcastd 3084(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t"
    "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t"
    "vpbroadcastd 4108(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm8, %%ymm8\n\t"
    "vpaddd %%ymm13, %%ymm9, %%ymm9\n\t"
    "vpbroadcastd 5132(%[a]), %%ymm13\n\t"
    "vpmaddwd %%ymm14, %%ymm13, %%ymm12\n\t"
    "vpmaddwd %%ymm15, %%ymm13, %%ymm13\n\t"
    "vpaddd %%ymm12, %%ymm10, %%ymm10\n\t"
    "vpaddd %%ymm13, %%ymm11, %%ymm11\n\t"
    "add $16, %[a]\n\t"
    "add $256, %[b]\n\t"
    "sub $4, %[k]\n\t"
    "jnz 1b\n\t"
    "vmovntdq %%ymm0, 0(%[c])\n\t"
    "vmovntdq %%ymm1, 32(%[c])\n\t"
    "vmovntdq %%ymm2, 2048(%[c])\n\t"
    "vmovntdq %%ymm3, 2080(%[c])\n\t"
    "vmovntdq %%ymm4, 4096(%[c])\n\t"
    "vmovntdq %%ymm5, 4128(%[c])\n\t"
    "vmovntdq %%ymm6, 6144(%[c])\n\t"
    "vmovntdq %%ymm7, 6176(%[c])\n\t"
    "vmovntdq %%ymm8, 8192(%[c])\n\t"
    "vmovntdq %%ymm9, 8224(%[c])\n\t"
    "vmovntdq %%ymm10, 10240(%[c])\n\t"
    "vmovntdq %%ymm11, 10272(%[c])\n\t"    : [a] "+&r"(pa), [b] "+&r"(pb2), [k] "+&r"(k)
    : [c] "r"(c)
    : "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}

static void product(Mat a, Mat b, int32_t *out) {
  pack_b(b);
  for (int ic=0; ic<H; ic+=MC) {
    int mc = H-ic < MC ? H-ic : MC;
    pack_a(a,ic,mc);
    for (int jr=0; jr<H; jr+=NR) {
      const int16_t *pb = bp + (size_t)(jr/NR)*KC*NR;
      int ir=0;
      int six_rows = mc;
      while (six_rows % 6 || (mc-six_rows) % 4) six_rows -= 2;
      for (; ir<six_rows; ir+=6)
        mikro6(ap+(size_t)ir*KC, pb, out+(size_t)(ic+ir)*H+jr);
      for (; ir+4<=mc; ir+=4)
        mikro(ap+(size_t)ir*KC, pb, out+(size_t)(ic+ir)*H+jr);
    }
  }
}

/* 32 output columns per iteration (was 8) => 28 independent loads in flight instead of 7.
   Pure loop restructuring: the arithmetic and every NT store address are unchanged. */
#define CB1(JJ)                                                                \
  do {                                                                         \
    const size_t t=(size_t)i*H+(JJ);                                           \
    __m256i m1=_mm256_load_si256((const __m256i *)(prod[0]+t));                \
    __m256i m2=_mm256_load_si256((const __m256i *)((int32_t *)qa+t));                \
    __m256i m3=_mm256_load_si256((const __m256i *)(prod[1]+t));                \
    __m256i m4=_mm256_load_si256((const __m256i *)((int32_t *)qb+t));                \
    __m256i m5=_mm256_load_si256((const __m256i *)(prod[2]+t));                \
    __m256i m6=_mm256_load_si256((const __m256i *)(prod[3]+t));                \
    __m256i m7=_mm256_load_si256((const __m256i *)(prod[4]+t));                \
    __m256i c11=_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1,m4),m5),m7); \
    __m256i c12=_mm256_add_epi32(m3,m5);                                       \
    __m256i c21=_mm256_add_epi32(m2,m4);                                       \
    __m256i c22=_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1,m3),m2),m6); \
    _mm256_stream_ps(C+(size_t)i*N+(JJ),       _mm256_mul_ps(_mm256_cvtepi32_ps(c11),inv)); \
    _mm256_stream_ps(C+(size_t)i*N+(JJ)+H,     _mm256_mul_ps(_mm256_cvtepi32_ps(c12),inv)); \
    _mm256_stream_ps(C+(size_t)(i+H)*N+(JJ),   _mm256_mul_ps(_mm256_cvtepi32_ps(c21),inv)); \
    _mm256_stream_ps(C+(size_t)(i+H)*N+(JJ)+H, _mm256_mul_ps(_mm256_cvtepi32_ps(c22),inv)); \
  } while (0)

static void combine(float *C) {
  _mm_sfence();                 /* the only fence needed: before reading the NT-written prod[] */
  const __m256 inv = _mm256_set1_ps(1.0f/(1024.0f*1024.0f));
  /* DESCENDING on purpose: the M2/M4 accumulators now live inside C's own buffer, and the
     float quadrant store for output row i lands exactly on the accumulator rows 2i,2i+1.
     Those are read at loop index 2i,2i+1 > i, i.e. already read when this row overwrites
     them.  Ascending order would destroy them before their read. */
  for (int i=H-1; i>=0; i--)
    for (int j=0; j<H; j+=64) { CB1(j); CB1(j+8); CB1(j+16); CB1(j+24); CB1(j+32); CB1(j+40); CB1(j+48); CB1(j+56); }
  _mm_sfence();
}
void matrix_multiply(int n, const float *A, const float *B, float *C) {
  (void)n;
  qa = (int16_t *)C; qb = (int16_t *)C + (size_t)N*N;   /* ★ 复用 C 的缓冲区当 quant 的输出 */
  quantize(A,qa);quantize(B,qb);_mm_sfence();
  const size_t off=(size_t)H*N;
  Mat a11{qa,nullptr,0},a12{qa+H,nullptr,0},
      a21{qa+off,nullptr,0},a22{qa+off+H,nullptr,0};
  Mat b11{qb,nullptr,0},b12{qb+H,nullptr,0},
      b21{qb+off,nullptr,0},b22{qb+off+H,nullptr,0};
  product({a11.x,a22.x,+1},{b11.x,b22.x,+1},prod[0]);
  product(a11,{b12.x,b22.x,-1},prod[1]);
  product({a11.x,a12.x,+1},b22,prod[2]);
  product({a21.x,a11.x,-1},{b11.x,b12.x,+1},prod[3]);
  product({a12.x,a22.x,-1},{b21.x,b22.x,+1},prod[4]);
  product({a21.x,a22.x,+1},b11,(int32_t *)qa);
  product(a22,{b21.x,b11.x,-1},(int32_t *)qb);
  combine(C);
}

评测结果