代码

// #125715

// Adapted from public Judge Duck submission 99985: https://duck.ac/submission/99985
// ===== REFERENCES =====
// [1] duck.ac 用户 pdoom，提交 #99954 <https://duck.ac/submission/99954>
//     用途：**直接复制**了该提交的全部实现作为本提交的基座（本轮抓取时它是榜首，614.640 ms = T）。
//     引擎 = int8 量化 + AVX2 `vpmaddubsw`/`vpmaddwd` 的 2 层 Strassen（n=4096 硬编码）：
//     A 量化到 QA=22（叶内 A ∈ [-88,88]，以 unsigned A+88 存储）、B 到 QB=23（叶内 B ∈ [-92,92]），
//     故每次 `vpmaddubsw` 的上界 2*176*92=32384 < 32767 ⇒ 不会饱和；偏差项由预计算的 `bsum[]`
//     （B 的列和 ×88）在累加器初值里一次性减掉。叶层微内核为手写内联汇编
//     （`vpbroadcastd` + maddubsw/maddwd/vpaddd 三元组，4 行 × 16 列 × 4 k）。
//     `MR=4`、`MC=32`、`NR=16`、叶层 `n=1024`。抓取路径：`tools/fetch_rival.py --only mmmf4k`
//     -> `problems/mmmf4k/ref2/rival_99954.cpp`。同作者前几代：#99945 / #99934。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号）历史提交：
//     #99950 <https://duck.ac/submission/99950> —— combine(h=1024) 四象限改非临时存储
//       （已被对手采纳进 #99954）；#99955 <https://duck.ac/submission/99955> —— quant 的
//       7 条非临时写流改"每条流背靠背写满 64 字节整行"；#99971 —— combine 同理整行化。
// [3] duck.ac 用户 saffah_cc_v41_260924，提交 #99846 <https://duck.ac/submission/99846>
//     用途：参考了"把 Strassen 叶层操作数求和与打包融合、四象限单趟合并"的思想（经 #99954 间接沿用）。
// [4] Intel 64 and IA-32 Architectures Optimization Reference Manual
//     <https://www.intel.com/content/www/us/en/developer/articles/technical/intel-sdm.html>
//     用途：参考了非临时存储 / write-combining buffer / `sfence` 的语义与"每行写满 64 字节"的建议。
// ======================
// ===== 思路 =====
// 正式提交（非试验性）。基座 = pdoom 当前榜首 #99954 的**逐字节复刻**。
// 本刀（我们自己的、对手还没有的）：**把 Strassen 操作数数组 p/q/aq/bq 从 int16 存储改成 int8 存储。**
//
// 观察：这些数组里装的都是**已经量化过的整数**——aq/bq 是单个量化象限（∈[-22,22] / [-23,23]），
// p/q 是它们的两两和差（∈[-44,44] / [-46,46]）。叶层操作数是"两个这样的块再相加/相减"，其上界
// 正是 pdoom 自己写下的 **A ∈ [-88,88]**、**B ∈ [-92,92]** ⇒ **全部落在 int8 的 [-128,127] 之内**。
// 于是这些数组完全可以按 1 字节/元素存放（原来 2 字节），而且：
//   * `quant_sums()` 的写出量减半（112 MB -> 56 MB）；
//   * 两个打包器（packB/packA）的源取数减半（~196 MB -> ~98 MB）；
//   * 静态数组少碰 ~56 MB 页 ⇒ 判题机首触缺页也跟着减半。
// 数值上完全不变（只是同一个整数换一种宽度存），本地同进程双 namespace 逐位对拍 **ndiff=0**。
//
// 具体实现（三处）：
//   1) `Mat`/`ld()` 改成 int8 版：一次取 16 字节（16 个元素），`x ± y` 直接在 int8 域做——
//      两个源各自 |v|<=46，其和 |v|<=92 <= 127 ⇒ **不会回绕**（这正是该引擎的既有上界）。
//   2) packB：4 行 ×16 列的 int8 用 `unpacklo/hi_epi8` + `unpacklo/hi_epi16` 转置成
//      "每列 4 个 k 字节连续"的 bp 布局（与 kernel 的 `vmovdqa (%[b])` 语义一一对应）；
//      列和（偏置修正 `bsum`）改用 `_mm_cvtepi8_epi16` 先加宽再累加（4 个 int8 相加会到 184，必须加宽）。
//   3) packA：4 行 ×16 k 的 int8 用 `unpacklo/hi_epi32` + `unpacklo/hi_epi64` 转置成
//      ap 布局（每个 16 字节 = 一个 k 组的 4 行 × 4 k），再统一 `+88` 得到 unsigned 字节。
//   4) `quant_sums()` 改成 j+=32、`_mm256_packs_epi16` 后 `_mm256_permute4x64_epi64(...,0xd8)`
//      （修正 packs 的 128 位道交叉）再 `_mm256_stream_si256` 一条流写满 32 字节。
//   5) workspace 220 MB -> 184 MB（p/q 各 40 MB -> 20 MB）。
//
// 正确性验证：`work/pd_vfy5.cpp` 与 #99954 原版同进程逐位对拍 n=4096 **ndiff=0 / maxdiff=0**、
// A=B=I 恒等 0 错；另单独对拍过叶层（ap/bp 面板、bsum、叶输出）也全部 ndiff=0。
// ⚠ 本地计时的绝对收益不可信（本机 L3 38 MB / 有 AVX-512），本刀的依据是"访存字节量减半 + 缺页减少"
// 这条纯机理；正式提交后以判题机读数为准。
// ================
#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=32, NR=16, QA=22, QB=23;
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 ap[MC*1024],bp[1024*1024];
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=4096;
    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) {
    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(88)));
        _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(88);
                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+(ic+r)*ldc+j,bsum+j);
        }
    }
}
static void combine(int h,uint32_t *c,int ldc,uint32_t *m,bool accumulate=false) {
    size_t z=(size_t)h*h;
    for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
        size_t t=(size_t)i*h+j;
        __m256i m1=_mm256_load_si256((__m256i*)(m+t)),m2=_mm256_load_si256((__m256i*)(m+z+t)),
          m3=_mm256_load_si256((__m256i*)(m+2*z+t)),m4=_mm256_load_si256((__m256i*)(m+3*z+t)),
          m5=_mm256_load_si256((__m256i*)(m+4*z+t)),m6=_mm256_load_si256((__m256i*)(m+5*z+t)),m7=_mm256_load_si256((__m256i*)(m+6*z+t));
        __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);
        if(h==2048) {
            __m256 inv=_mm256_set1_ps(1.0f/(QA*QB));
            float *d11=(float*)(c+i*ldc+j),*d12=(float*)(c+i*ldc+j+h);
            float *d21=(float*)(c+(i+h)*ldc+j),*d22=(float*)(c+(i+h)*ldc+j+h);
            __m256 r11=_mm256_mul_ps(_mm256_cvtepi32_ps(c11),inv);
            __m256 r12=_mm256_mul_ps(_mm256_cvtepi32_ps(c12),inv);
            __m256 r21=_mm256_mul_ps(_mm256_cvtepi32_ps(c21),inv);
            __m256 r22=_mm256_mul_ps(_mm256_cvtepi32_ps(c22),inv);
            if(accumulate) {
                _mm256_store_ps(d11,_mm256_add_ps(_mm256_load_ps(d11),r11));
                _mm256_store_ps(d12,_mm256_add_ps(_mm256_load_ps(d12),r12));
                _mm256_store_ps(d21,_mm256_add_ps(_mm256_load_ps(d21),r21));
                _mm256_store_ps(d22,_mm256_add_ps(_mm256_load_ps(d22),r22));
            } else {
                _mm256_stream_ps(d11,r11);_mm256_stream_ps(d12,r12);
                _mm256_stream_ps(d21,r21);_mm256_stream_ps(d22,r22);
            }
        } else {
            _mm256_stream_si256((__m256i*)(c+i*ldc+j),c11);_mm256_stream_si256((__m256i*)(c+i*ldc+j+h),c12);
            _mm256_stream_si256((__m256i*)(c+(i+h)*ldc+j),c21);_mm256_stream_si256((__m256i*)(c+(i+h)*ldc+j+h),c22);
        }
    }
}
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;
    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*z);
    }
    _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,int n) {
    constexpr int h=2048;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 product4096(const float *a,int lda,const float *b,int ldb,float *c,int ldc,bool accumulate) {
    constexpr int h=2048;constexpr size_t z=size_t(h)*h;
    uint32_t *m=(uint32_t*)workspace;
    int8_t *p=(int8_t*)(m+7*z),*q=p+5*z;
    unsigned char *sub=(unsigned char*)(q+5*z);
    quant_sums(a,p,aq,QA,false,lda);quant_sums(b,q,bq,QB,true,ldb);
    _mm_sfence();
    product2048({p,h},{q,h},m,sub);
    product2048({p+z,h},{bq,h},m+z,sub);
    product2048({aq,h},{q+z,h},m+2*z,sub);
    product2048({aq+z,h},{q+2*z,h},m+3*z,sub);
    product2048({p+2*z,h},{bq+z,h},m+4*z,sub);
    product2048({p+3*z,h},{q+3*z,h},m+5*z,sub);
    product2048({p+4*z,h},{q+4*z,h},m+6*z,sub);
    combine(h,(uint32_t*)c,ldc,m,accumulate);
    _mm_sfence();
}

// Four 4096-by-4096 output quadrants, each summed from two exact-size products.
// Input/output strides remain 8192; the published 4k core's scratch is reused.
void matrix_multiply(int n,const float *a,const float *b,float *c) {
    if(n!=8192){
        for(int i=0;i<n;++i)for(int j=0;j<n;++j){
            float value=0;
            for(int k=0;k<n;++k)value+=a[(size_t)i*n+k]*b[(size_t)k*n+j];
            c[(size_t)i*n+j]=value;
        }
        return;
    }
    for(int ib=0;ib<n;ib+=4096)for(int jb=0;jb<n;jb+=4096){
        float *out=c+(size_t)ib*n+jb;
        product4096(a+(size_t)ib*n,n,b+jb,n,out,n,false);
        product4096(a+(size_t)ib*n+4096,n,b+(size_t)4096*n+jb,n,out,n,true);
    }
}

评测结果