代码

// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #103335
//     <https://duck.ac/submission/103335>
//     用途：**参考了思想**（未复制任何代码行）——该提交的正文经 `tools/rivalcopy.py` 比对，
//     是**本账号 #103002 的逐字节副本**（fwd/rev 98.47/98.47%），只差两处编译期常量：
//     叶层 A 面板行块 `MC 32->128`、以及七个叶乘积面板各加 256 字节间隔（`z = h*h -> h*h+64`）。
//     本发就是把这**两处常量**搬回我们自己的基座（#103002），并把它们各自定价。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #103002
//     <https://duck.ac/submission/103002>
//     用途：**本文件的正文即该提交的正文**（只改上面那两处常量）。该件本身是
//     int8 六比特量化 + 一层 Strassen + 手写 AVX2 叶层微内核 + 融合进 combine 的 int32->float。
// [3] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #99985
//     <https://duck.ac/submission/99985>
//     用途：本引擎的原始出处（mmmf4k 达标题），经 #102820 以 n=2048 重新参数化后由 [2] 继承。
// [4] duck.ac 用户 pdoom，提交 #99954 <https://duck.ac/submission/99954>
//     用途：引擎算法（int8 量化 + `vpmaddubsw` 的 unsigned×signed 字节通道 + 预计算列和
//     `bsum` 偏置修正）源自该提交。
// [5] /home/yjp/duck.ac/BRIEF.md §2.18.653（当对手那一发**就是我们自己的源码**时，复刻/借鉴
//     没有罚项，且由此得到一把 ±0.04% 分辨率的单变量定价尺）、§2.18.666（布局类旋钮要扫整条轴）、
//     §2.18.638（判题机 = i3-8100/Skylake，L1d 32K/8路、L2 256K/4路、L3 6M/12路）。
// ======================
// ===== 思路 =====
// 局势：`mine=86.088445`（我们 #103400，MC128+ZPAD64）、`T=86.002133`（对手 #103335）⇒
//   严支真目标 `0.99*T+1us = 85.143112 ms`，真缺口 0.945 ms = 1.10%。
//
// 本发 = #103400 + **七个叶乘积面板的 16 列(64B)粒度交织**（唯一的改动）：
//   原布局：叶 r 的 (i,j) 在 `m + r*z + i*h + j`，七块相距 4MB。
//   新布局：叶 r 的 (i,j) 在 `m + (((i*64 + j/16)*7 + r)*16 + (j&15))`（同一块 28MB 区的置换）。
//   ⇒ combine 每个 (i,j) 要读的 7 个值**连续**（相距 16 个 uint32），沿 j 推进时整块 m 变成
//     **一条顺序流**（28MB），而不是 7 条相距 4MB 的流。
//   ⇒ 代价：叶层 kernel 的四行输出间距 4096B -> 28672B（仍每行一条完整 64B cache line 的
//     `vmovntdq`，NT 写总量只有 0.11 B/cyc，远在写带宽之内）。
//
// 判题机同进程三臂 probe（`work/f2c_probe4.cpp`，正反两序取几何平均，三臂 C 缓冲逐字节相同）：
//   基线(#103400) 总 296.11 Mcyc，其中 combine 10.47 Mcyc；
//   本发        总 294.07 Mcyc，其中 combine **7.73 Mcyc（−26%）**，kernel 268.42（+0.24%，噪声内）
//   ⇒ combine 折成 44MB/(7.73Mcyc/3.4364GHz)=19.6 GB/s = 本机可达带宽（对照 quantize 20.5 GB/s）。
//   预测 86.088 * 0.99311 = 85.495 ms。臂2（原布局 + combine 7 流软预取 256B）只 −0.14%，未采用。
// ================
// ===== 思路（本发 #1034xx 的唯一改动）=====
// ★ 修复一条**被注释吞掉的 `#pragma GCC optimize`**。
//   基座 `work/f2c_sub_ilv.cpp`（= #103427）第 39 行的原文是
//        `// ================#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")`
//   —— 生成器 `work/f2c_mksub.py` 把 `--idea` 串直接拼在正文前，而该串末尾少一个换行，
//   于是整条 pragma 被并进了上一行的 `//` 注释里。
//   ⇒ #103400 / #103427 / #103447 三发**全部按默认 -O2 编译**，而对手 #103335 的同一行是活的
//   （`grep -c 'pragma GCC'` 双方都是 2，但只有他那一条进编译单元）。
//   本发只补回这一个换行：正文其余部分与本文件基座**逐字节相同**。
//   依据：同族 `1008` 用"只删那一行、其余逐字节相同"的头实测
//   `#pragma GCC optimize("O3,unroll-loops,rename-registers")` = **−1.57%**（双方原本都跑 -O2）；
//   本族 `noi17f` 补 `target(...)` = −0.664 ms。本次为**试验性提交**：目的是给这一条 pragma
//   在本题的当前基座上单变量定价（`§2.18.334`：换基座必须重新量）。
// 正确性：本机（g++ 15）与 `f2b_verify.cpp` 驱动对拍 —— 随机输入 maxdiff=1.320875（容差 1.5）、
//   A=B=I 恒等 ndiff=0、全 1 输入 maxdiff=0 ⇒ 与基座逐位一致（pragma 只改代码生成）。
// ================
// ===== 思路 =====
// 局势：`mine = 85.615430`（#104088）、`T = 86.002133`（对手 #103335）、严支 `thr = 85.143112`
//   ⇒ 严支真缺口 0.472318 ms = 0.552%。
//
// ★ 本发（单变量）：**把 ap / bp 两个打包面板的基址移出 4 KB 边界 64 B**（步长一律不动）。
//
// 一、取证（第一动作，本机实测真实落点）：
//     aq / bq / workspace / ap / bp / bsum 的地址 **全部 mod 4096 = 0**。
//     ⇒ `aq`/`bq`（源，`alignas(4096)`，行距 2048 B）的行相位恒为 {0, 2048}，
//       而 `ap`/`bp`（目的，行距 1024 B，基址 4 KB 对齐）的行相位恒为 0
//     ⇒ **打包循环里"源与目的 4 KB 同余"**：与 pending store 仅差 4 KB 倍数的 load
//       无法前递，拷贝循环被串行化（§2.18.474 的经典形态）。
//
// 二、旁证（同批已量）：把 packA/packB 改成"`sign` 为编译期模板参数"后，
//     godbolt **gcc 9.3** 上 packA 内层循环 **87 → 27 条指令**（零栈重载、零条件跳转），
//     但**正式提交实测 Δ = +0.005%（零，#104258）**
//     ⇒ ★ **packA / packB 不是指令受限，而是访存受限** ⇒ 本发转而攻访存相位。
//
// 三、做法：`alignas(64) uint8_t apRaw[MC*1024+64], bpRaw[1024*1024+64]` + 两个 `#define` 偏移 64 B。
//     仍满足 `_mm_store_si128` 的 16 B 对齐（`alignas(64)` + 64 仍是 64 的倍数）。
//     ★ **步长（行距 / 组距 / 索引表达式）一个字都不动** —— 这是本族的关键规程：
//       之前试过的"改行距/组距（+8）"会同时改掉 mikro 的取数步长，两个变量混在一起；
//       只动基址相位就只影响 aliasing（mmmd1k #99797 上正是这样拿到 −148 µs = −0.35%）。
//
// 四、数值语义**逐位不变**（只改缓冲区分区的起始相位）。
//     正确性（本机 g++ 15 + `f2b_verify.cpp` 驱动）：随机输入 maxdiff = 1.320875（容差 1.5）、
//       A=B=I 恒等 ndiff = 0、全 1 输入 maxdiff = 0。
// 本次为**试验性提交**：目的是给"4 KB 同余"这条轴在本题当前基座上单变量定价。
// ================
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("avx2,fma")
#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>
static constexpr int MR=4, MC=128, NR=16, QA=45, QB=45;
alignas(4096) static int8_t aq[2*2048*2048], bq[2*2048*2048];
alignas(4096) static unsigned char workspace[184*1024*1024];
alignas(64) static uint8_t apRaw[MC*1024+64],bpRaw[1024*1024+64];
#define ap (apRaw+64)
#define bp (bpRaw+64)
alignas(64) static int bsum[1024];
struct Mat { const int8_t *x; int ld; const int8_t *y=nullptr; int ly=0; int sign=0; };
// 取 16 个 int8；x±y 都在 int8 域完成（各源 |v|<=46 => 和 |v|<=92 <=127，不溢出）
static inline __m128i ld(Mat a,int i,int j) {
    __m128i x=_mm_loadu_si128((const __m128i*)(a.x+i*a.ld+j));
    if(a.sign) {__m128i y=_mm_loadu_si128((const __m128i*)(a.y+i*a.ly+j));
        x=a.sign>0?_mm_add_epi8(x,y):_mm_sub_epi8(x,y);}
    return x;
}
static inline __m128i load4(Mat a,int i,int k) {
    __m128i x=_mm_loadl_epi64((const __m128i*)(a.x+i*a.ld+k));
    if(a.sign) {
        __m128i y=_mm_loadl_epi64((const __m128i*)(a.y+i*a.ly+k));
        x=a.sign>0?_mm_add_epi16(x,y):_mm_sub_epi16(x,y);
    }
    return x;
}
static __attribute__((noinline)) void kernel(const uint8_t *a,const uint8_t *b,uint32_t *out,const int *sums) {
    ptrdiff_t stride=28672;
    int k=256;
    asm volatile (
        "vpxor %%ymm14, %%ymm14, %%ymm14\n\t"
        "vmovdqa (%[sums]), %%ymm11\n\t"
        "vmovdqa 32(%[sums]), %%ymm12\n\t"
        "vpsubd %%ymm11, %%ymm14, %%ymm0\n\t"
        "vpsubd %%ymm12, %%ymm14, %%ymm1\n\t"
        "vpsubd %%ymm11, %%ymm14, %%ymm2\n\t"
        "vpsubd %%ymm12, %%ymm14, %%ymm3\n\t"
        "vpsubd %%ymm11, %%ymm14, %%ymm4\n\t"
        "vpsubd %%ymm12, %%ymm14, %%ymm5\n\t"
        "vpsubd %%ymm11, %%ymm14, %%ymm6\n\t"
        "vpsubd %%ymm12, %%ymm14, %%ymm7\n\t"
        "vpcmpeqw %%ymm10, %%ymm10, %%ymm10\n\t"
        "vpsrlw $15, %%ymm10, %%ymm10\n\t"
        ".p2align 5\n\t"
        "1:\n\t"
        "vmovdqa 0(%[b]), %%ymm11\n\t"
        "vmovdqa 32(%[b]), %%ymm12\n\t"
        "vpbroadcastd 0(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm0, %%ymm0\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 4(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 8(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 12(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 64(%[b]), %%ymm11\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vpbroadcastd 16(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm0, %%ymm0\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 20(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 24(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 28(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 128(%[b]), %%ymm11\n\t"
        "vmovdqa 160(%[b]), %%ymm12\n\t"
        "vpbroadcastd 32(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm0, %%ymm0\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 36(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 40(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 44(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 192(%[b]), %%ymm11\n\t"
        "vmovdqa 224(%[b]), %%ymm12\n\t"
        "vpbroadcastd 48(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm0, %%ymm0\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 52(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 56(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 60(%[a]), %%ymm13\n\t"
        "vpmaddubsw %%ymm11, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
        "vpmaddubsw %%ymm12, %%ymm13, %%ymm14\n\t"
        "vpmaddwd %%ymm10, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "add $64, %[a]\n\t"
        "add $256, %[b]\n\t"
        "sub $4, %[k]\n\t"
        "jnz 1b\n\t"
        "vmovntdq %%ymm0, (%[out])\n\t"
        "vmovntdq %%ymm1, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntdq %%ymm2, (%[out])\n\t"
        "vmovntdq %%ymm3, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntdq %%ymm4, (%[out])\n\t"
        "vmovntdq %%ymm5, 32(%[out])\n\t"
        "add %[stride], %[out]\n\t"
        "vmovntdq %%ymm6, (%[out])\n\t"
        "vmovntdq %%ymm7, 32(%[out])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
        : [stride] "r"(stride), [sums] "r"(sums)
        : "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
static void leaf(Mat a,Mat b,uint32_t *c,int lr) {
    constexpr int n=1024,ldc=1024;
    {
        constexpr int pc=0,kc=1024;
        for(int j=0;j<n;j+=8) _mm256_store_si256((__m256i*)(bsum+j),_mm256_setzero_si256());
        for(int k=0;k<kc;k+=4) {
            for(int j=0;j<n;j+=16) {
                __m256i s0=_mm256_load_si256((const __m256i*)(bsum+j));
                __m256i s1=_mm256_load_si256((const __m256i*)(bsum+j+8));
                __m128i v0=ld(b,pc+k,j),v1=ld(b,pc+k+1,j),v2=ld(b,pc+k+2,j),v3=ld(b,pc+k+3,j);
                __m128i wa=_mm_add_epi16(_mm_add_epi16(_mm_cvtepi8_epi16(v0),_mm_cvtepi8_epi16(v1)),
                                         _mm_add_epi16(_mm_cvtepi8_epi16(v2),_mm_cvtepi8_epi16(v3)));
                __m128i wb=_mm_add_epi16(_mm_add_epi16(_mm_cvtepi8_epi16(_mm_srli_si128(v0,8)),_mm_cvtepi8_epi16(_mm_srli_si128(v1,8))),
                                         _mm_add_epi16(_mm_cvtepi8_epi16(_mm_srli_si128(v2,8)),_mm_cvtepi8_epi16(_mm_srli_si128(v3,8))));
                s0=_mm256_add_epi32(s0,_mm256_cvtepi16_epi32(wa));
                s1=_mm256_add_epi32(s1,_mm256_cvtepi16_epi32(wb));
                __m128i t0=_mm_unpacklo_epi8(v0,v1),t1=_mm_unpackhi_epi8(v0,v1);
                __m128i t2=_mm_unpacklo_epi8(v2,v3),t3=_mm_unpackhi_epi8(v2,v3);
                uint8_t *p=bp+j*kc+k*16;
                _mm_store_si128((__m128i*)(p),   _mm_unpacklo_epi16(t0,t2));
                _mm_store_si128((__m128i*)(p+16),_mm_unpackhi_epi16(t0,t2));
                _mm_store_si128((__m128i*)(p+32),_mm_unpacklo_epi16(t1,t3));
                _mm_store_si128((__m128i*)(p+48),_mm_unpackhi_epi16(t1,t3));
                _mm256_store_si256((__m256i*)(bsum+j),s0);
                _mm256_store_si256((__m256i*)(bsum+j+8),s1);
            }
        }
        for(int j=0;j<n;j+=8) _mm256_store_si256((__m256i*)(bsum+j),_mm256_mullo_epi32(_mm256_load_si256((const __m256i*)(bsum+j)),_mm256_set1_epi32(90)));
        _mm_sfence();
        for(int ic=0;ic<n;ic+=MC) {
            int mc=n-ic<MC?n-ic:MC;
            for(int r=0;r<mc;r+=MR) {
                __m128i *p=(__m128i*)(ap+r*kc);
                __m128i b88=_mm_set1_epi8(90);
                for(int k=0;k<kc;k+=16) {
                    __m128i v0=ld(a,ic+r,k),v1=ld(a,ic+r+1,k),v2=ld(a,ic+r+2,k),v3=ld(a,ic+r+3,k);
                    __m128i s1=_mm_unpacklo_epi32(v0,v1),s2=_mm_unpackhi_epi32(v0,v1);
                    __m128i s3=_mm_unpacklo_epi32(v2,v3),s4=_mm_unpackhi_epi32(v2,v3);
                    __m128i o0=_mm_unpacklo_epi64(s1,s3),o1=_mm_unpackhi_epi64(s1,s3);
                    __m128i o2=_mm_unpacklo_epi64(s2,s4),o3=_mm_unpackhi_epi64(s2,s4);
                    _mm_store_si128(p++,_mm_add_epi8(o0,b88));
                    _mm_store_si128(p++,_mm_add_epi8(o1,b88));
                    _mm_store_si128(p++,_mm_add_epi8(o2,b88));
                    _mm_store_si128(p++,_mm_add_epi8(o3,b88));
                }
            }
            for(int j=0;j<n;j+=NR) for(int r=0;r<mc;r+=MR)
                kernel(ap+r*kc,bp+j*kc,c+(((size_t)((ic+r)*64+(j>>4))*7+(size_t)lr)*16),bsum+j);
        }
    }
}
static void combine(int h,uint32_t *c,int ldc,uint32_t *m) {
    size_t z=(size_t)h*h+64;
    for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
        size_t t=((size_t)(i*64+(j>>4))*7)*16+(size_t)(j&15);
        __m256i m1=_mm256_load_si256((__m256i*)(m+t)),m2=_mm256_load_si256((__m256i*)(m+t+16)),
          m3=_mm256_load_si256((__m256i*)(m+t+32)),m4=_mm256_load_si256((__m256i*)(m+t+48)),
          m5=_mm256_load_si256((__m256i*)(m+t+64)),m6=_mm256_load_si256((__m256i*)(m+t+80)),m7=_mm256_load_si256((__m256i*)(m+t+96));
        __m256i c11=_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1,m4),m5),m7);
        __m256i c12=_mm256_add_epi32(m3,m5),c21=_mm256_add_epi32(m2,m4);
        __m256i c22=_mm256_add_epi32(_mm256_add_epi32(_mm256_sub_epi32(m1,m2),m3),m6);
        /* fused: convert and store as float here, so no separate int32->float pass is needed */
        {
            __m256 inv=_mm256_set1_ps(1.0f/(QA*QB));
            _mm256_stream_ps((float*)(c+i*ldc+j),_mm256_mul_ps(_mm256_cvtepi32_ps(c11),inv));
            _mm256_stream_ps((float*)(c+i*ldc+j+h),_mm256_mul_ps(_mm256_cvtepi32_ps(c12),inv));
            _mm256_stream_ps((float*)(c+(i+h)*ldc+j),_mm256_mul_ps(_mm256_cvtepi32_ps(c21),inv));
            _mm256_stream_ps((float*)(c+(i+h)*ldc+j+h),_mm256_mul_ps(_mm256_cvtepi32_ps(c22),inv));
        }
    }
}
static void product2048(Mat a,Mat b,uint32_t *c,unsigned char *work) {
    constexpr int h=1024;constexpr size_t z=size_t(h)*h+64;
    Mat aa[4]={{a.x,a.ld},{a.x+h,a.ld},{a.x+h*a.ld,a.ld},{a.x+h*a.ld+h,a.ld}};
    Mat bb[4]={{b.x,b.ld},{b.x+h,b.ld},{b.x+h*b.ld,b.ld},{b.x+h*b.ld+h,b.ld}};
    uint32_t *m=(uint32_t*)work;
    const int ai[7]={0,2,0,3,0,2,1},aj[7]={3,3,-1,-1,1,0,3},as[7]={1,1,0,0,1,-1,-1};
    const int bi[7]={0,0,1,2,3,0,2},bj[7]={3,-1,3,0,-1,1,3},bs[7]={1,0,-1,-1,0,1,1};
    for(int r=0;r<7;r++) {
        Mat x=aa[ai[r]],y=bb[bi[r]];
        if(as[r]) {x.y=aa[aj[r]].x;x.ly=a.ld;x.sign=as[r];}
        if(bs[r]) {y.y=bb[bj[r]].x;y.ly=b.ld;y.sign=bs[r];}
        leaf(x,y,m,r);
    }
    _mm_sfence();
    combine(h,c,2048,m);
    _mm_sfence();
}
static inline __m256i quant16(const float *a,__m256 scale) {
    __m256i x=_mm256_cvtps_epi32(_mm256_mul_ps(_mm256_load_ps(a),scale));
    __m256i y=_mm256_cvtps_epi32(_mm256_mul_ps(_mm256_load_ps(a+8),scale));
    return _mm256_permute4x64_epi64(_mm256_packs_epi32(x,y),0xd8);
}
static void quant_sums(const float *a,int8_t *dst,int8_t *diag,int q,bool isb) {
    constexpr int h=2048,n=4096;constexpr size_t z=size_t(h)*h;
    const __m256 scale=_mm256_set1_ps(float(q));
    for(int i=0;i<h;i++)for(int j=0;j<h;j+=32) {
        __m256i x0a=quant16(a+i*n+j,scale),x0b=quant16(a+i*n+j+16,scale);
        __m256i x1a=quant16(a+i*n+j+h,scale),x1b=quant16(a+i*n+j+h+16,scale);
        __m256i x2a=quant16(a+(i+h)*n+j,scale),x2b=quant16(a+(i+h)*n+j+16,scale);
        __m256i x3a=quant16(a+(i+h)*n+j+h,scale),x3b=quant16(a+(i+h)*n+j+h+16,scale);
        size_t t=(size_t)i*h+j;
        __m256i v0a=_mm256_add_epi16(x0a,x3a),v0b=_mm256_add_epi16(x0b,x3b);
        __m256i v1a=isb?_mm256_sub_epi16(x1a,x3a):_mm256_add_epi16(x2a,x3a);
        __m256i v1b=isb?_mm256_sub_epi16(x1b,x3b):_mm256_add_epi16(x2b,x3b);
        __m256i v2a=isb?_mm256_sub_epi16(x2a,x0a):_mm256_add_epi16(x0a,x1a);
        __m256i v2b=isb?_mm256_sub_epi16(x2b,x0b):_mm256_add_epi16(x0b,x1b);
        __m256i v3a=isb?_mm256_add_epi16(x0a,x1a):_mm256_sub_epi16(x2a,x0a);
        __m256i v3b=isb?_mm256_add_epi16(x0b,x1b):_mm256_sub_epi16(x2b,x0b);
        __m256i v4a=isb?_mm256_add_epi16(x2a,x3a):_mm256_sub_epi16(x1a,x3a);
        __m256i v4b=isb?_mm256_add_epi16(x2b,x3b):_mm256_sub_epi16(x1b,x3b);
        _mm256_stream_si256((__m256i*)(dst+t),     _mm256_permute4x64_epi64(_mm256_packs_epi16(v0a,v0b),0xd8));
        _mm256_stream_si256((__m256i*)(dst+z+t),   _mm256_permute4x64_epi64(_mm256_packs_epi16(v1a,v1b),0xd8));
        _mm256_stream_si256((__m256i*)(dst+2*z+t), _mm256_permute4x64_epi64(_mm256_packs_epi16(v2a,v2b),0xd8));
        _mm256_stream_si256((__m256i*)(dst+3*z+t), _mm256_permute4x64_epi64(_mm256_packs_epi16(v3a,v3b),0xd8));
        _mm256_stream_si256((__m256i*)(dst+4*z+t), _mm256_permute4x64_epi64(_mm256_packs_epi16(v4a,v4b),0xd8));
        _mm256_stream_si256((__m256i*)(diag+t),    _mm256_permute4x64_epi64(_mm256_packs_epi16(x0a,x0b),0xd8));
        _mm256_stream_si256((__m256i*)(diag+z+t),  _mm256_permute4x64_epi64(_mm256_packs_epi16(x3a,x3b),0xd8));
    }
}
static void quant_all(const float *src,int8_t *dst,int scale) {
    const __m256 q=_mm256_set1_ps((float)scale);
    for(size_t i=0;i<(size_t)2048*2048;i+=32) {
        __m256i x=quant16(src+i,q),y=quant16(src+i+16,q);
        _mm256_store_si256((__m256i*)(dst+i),
          _mm256_permute4x64_epi64(_mm256_packs_epi16(x,y),0xd8));
    }
}
void matrix_multiply(int n,const float *a,const float *b,float *c) {
    (void)n;
    quant_all(a,aq,QA);quant_all(b,bq,QB);
    product2048(Mat{aq,2048},Mat{bq,2048},(uint32_t*)c,workspace);
    /* the int32->float conversion is fused into combine() */
    _mm_sfence();
}

评测结果