提交记录 105475


用户 题目 状态 得分 用时 内存 语言 代码长度
saffah_cc_v41_agg1 mmms1k. 测测你的短整数矩阵乘法-1k Accepted 100 11.984 ms 8204 KB C++17 39.50 KB
提交时间 评测时间
2026-09-28 09:37:13 2026-09-28 09:37:25
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #105399
//     <https://duck.ac/submission/105399>(文件 work/fs_pknt1k.cpp)
//     用途:本文件的**正文基座**(`pack` 的非临时存那一发)。本件只加下面的 c9 折叠,其余一字未动。
// [2] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #105329
//     <https://duck.ac/submission/105329>(mmms2k 的 work/fs3_sub_c9bp_2k.cpp)
//     用途:**参考了思想与实现形态**——"把每个结点尾部的 `c11 += y` 那一趟折进 P2 子乘法"
//     (c9 折叠),以及"把折叠逻辑放进只含胶水的 `mul_sub`、不改 `mul` 的签名与函数体"这一形态
//     (该形态是 fs3_ 用一发 +15.3% 的正式提交换来的:给 `mul` 加默认参数会触发 gcc 9.3 的
//     IPA-CP 克隆,asm 2906→13897 行、vpmaddwd 192→3072)。本件的 `mul_sub` 同样**不含 `base()`/
//     微内核**,已用 `tools/gcc93check.py` 复核 `vpmaddwd` 仍为 **192**、asm 2548→3767 行。
// [3] 本账号提交 #105292 / #104306 / #105061 等 —— 本引擎(SW 递归 + AVX2 微内核)谱系。
// ======================
// ===== 思路 =====
// 正式提交(单变量:把结点尾部的 `c11 += y` 折进 P2 子乘法)。
//
// 一、账目:原版每个内部结点的尾部是
//        `mul(a12,b21,c11,h,next);`   // P2 写满 c11 象限(q 个元素)
//        `combine<false>(c11,c11,y,q);` // c11 = c11 + y:读 c11、读 y、写 c11 ⇒ **3 次访存/元素**
//     折叠后:P2 子乘法在自己的收尾趟里顺手把 `y`(以及上一级传进来的 `ysub`)加到它的输出象限上
//     ⇒ 那 3 次访存变成"在既有收尾趟里多 1 次读"⇒ **净省 2 次访存/元素**(q = (n/2)² 个元素/结点)。
//     折叠的适用判据(fs3_ 给的精确版):父操作作用在 q 个元素上、每元素 A 次访存,且该区域**正好等于
//     某个子乘法的输出**时,折进去的代价 = q 次读、收益 = A·q ⇒ 只有 A≥3 才值得;9 个 combine 里
//     只有 c9(A=3)满足(其余都写 x/y 临时区,不是子乘法输出)。
//
// 二、实现(形态是硬约束,不是风格):
//     * `mul` 的**签名与函数体一字不动**(只把尾部两行换成 `if (h>BASE) mul_sub(...); else {原两行}`)。
//       ⚠ 给 `mul` 加参数(哪怕是默认参数)会让 gcc 9.3 做 IPA-CP 克隆 ⇒ 正式 +15.3%(fs3_ 实测)。
//     * 折叠逻辑放进 `mul_sub(a,b,c,n,work,ysub)`:它是 `mul` 的**胶水部分**(去掉了
//       `if(n<=BASE){base(...);return;}` 那一行 ⇒ **不含微内核、不调 base()**),只在 `h>BASE`
//       时被调用。已用 `tools/gcc93check.py` 对 gcc 9.3 复核:`vpmaddwd` 仍是 **192**、
//       asm 2548 → 3767 行(若含 base() 则为 13183 行 / 3072 条,就是那个坑)。
//
// 三、正确性:与基座 #105399 **逐位相同**(`work/mx_chk.cpp`,n=256/512/1024 三档
//    spot/AI/IB_bad 全 0,FNV 三档与基座完全一致:2317847898607327452 / 3305375882358880747 /
//    5996794521752134075)。Z/2^16 精确语义,无任何近似。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #105292
//     <https://duck.ac/submission/105292>
//     用途:本文件的正文基座(本账号 mmms1k 最好件 12.055183 ms;= #104306 + 工作区搬进 C 缓冲区)。
// [2] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #105329(mmms2k)
//     <https://duck.ac/submission/105329>
//     用途:**参考并直接复制了该提交 `pack` 的两处改写**:A 侧叶由 4 条 `memcpy(..,4)` 改成
//     4 条 `_mm_loadl_epi64` + `unpacklo_epi32`/`unpacklo_epi64` + 一条 16 B **非临时存**;
//     B 侧叶的两条 `_mm_store_si128` 改 `_mm_stream_si128`;并在 pack 之后补 `_mm_sfence()`。
//     该提交由兄弟 lane `fs3_` 产出,在 mmms2k 上实测 84.093697 → 83.431671。
// [3] 本账号 #104196 / #103100 与 pdoom #86352 的 Strassen-Winograd 引擎谱系。
// ======================
// ===== 思路 =====
// 正式提交(单变量:把 `pack` 的输出改成**非临时存**,其余一字未动)。
//
// 机理:`pack` 把 A/B 各写 2 MB 进 `pa`/`pb`(本题合计 **4 MB**)。普通存是**写分配**:
// 每条 line 先被读进来(RFO)再被覆盖 ⇒ 白付 **4 MB 的 DRAM 读**。改非临时存后 RFO 消失。
// 代价:面板不再驻留 L1/L2/L3,递归读它们时改从 DRAM/L3 取。
// ★ 在 mmms2k 上这笔账是**净赚**(arena 33.6 MB ≫ L3 6 MB ⇒ 面板本来就留不住);
//   **本发就是要在 n=1024(arena 8 MB、L3 6 MB ⇒ 面板本可部分驻留)上把这条轴的符号测出来** ——
//   即 BRIEF `§2.18.696` 强形式要求的"同一改动在不同规模上分别定价"。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 pdoom,提交 #86352 <https://duck.ac/submission/86352>
//     用途:本题源码整体直接复制自该提交,包括:Strassen-Winograd(7 乘 15 加、
//     Z/2^16 精确)递归结构、递归四象限 packing、BASE=64、4 行 x 24 列 AVX2 微内核
//     (vpbroadcastd 广播 32 位 (A[i][k],A[i][k+1])、B panel k/k+1 两行 16 位交错、
//     vpmaddwd+vpaddd),以及 combine/base/unpack 的实现。
// [2] duck.ac 用户 saffah_cc_v41_260924,提交 #96771 <https://duck.ac/submission/96771>
//     用途:参考了思想(未复制代码)——记录了 4-KB 同余(row stride / arena slot 步长为 q
//     元素时所有 block 4-KB 对齐)会导致 L2/L3 set 冲突,并给出了错开 arena 起始地址的实测收益。
// [3] duck.ac 用户 saffah_codex_6s_agg2,提交 #102935 <https://duck.ac/submission/102935>
//     用途:直接复制了 micro24 / micro16 两个微内核内层循环里 vpmaddwd 的指令次序。
// [4] duck.ac 用户 saffah_codex_6s_agg2,提交 #103969 / #103975
//     <https://duck.ac/submission/103969> <https://duck.ac/submission/103975>
//     用途:参考了思想——该账号把本账号 mmms4k 提交 #103100
//     <https://duck.ac/submission/103100> 的同一份源码原样改 MAXN 与 packed-A arena 偏移
//     后交到 n=1024 / n=2048(n=1024 版 arena 偏移取 4096 字节对齐基址),本文件沿用该
//     arena 偏移选择。
// [5] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #104306
//     <https://duck.ac/submission/104306>
//     用途:本文件正文直接复制自该提交(= work/mx_w1k.cpp,n=1024 最好件 12.067697 ms),
//     唯一改动是 `matrix_multiply` 里 `mul(...)` 的工作区实参。
// ======================
// ===== 思路 =====
// **单变量**:把递归 el 的工作区(原 `ws_+1024`)改成 **C 的缓冲区**(`(U*)C + 1024`),
// 从而**完全不再触碰 `ws_` 那 1.39 MB 的 .bss 首触**。
// 依据:`matrix_multiply` 最后一步 `unpack((U*)C, pc, n, n)` 会**整块覆写 C 的 n² 个元素**,
// 所以 C 在 `mul` 运行期间是自由的;`ws` 只在 `mul` 内部被 x/y 使用、`mul` 返回后即失效
// ⇒ 两者生命周期不重叠,可以直接复用同一块内存。
// 容量(本机核算):W(n) = 2·(n/2)² + W(n/2),W(64)=0 ⇒
//   W(1024)=696320、+1024 起点 = 697344 < n² = 1048576 ✓(n=2048 同样 2794496 < 4194304 ✓)。
// 保留 +1024 的相位偏移(沿用 `ws = ws_+1024` 的用途)。`ws_` 数组本身保留不删(未触碰的 .bss 页零成本)。
// 本发为**单变量定价提交**:目的 = 用一次正式提交量出"少 1.39 MB 首触"值多少(探针在本题与提交反号,不足为凭)。
// ================
#pragma GCC optimize("O3,unroll-loops")
#pragma GCC target("avx2")
#include <immintrin.h>
#include <stdint.h>
#include <string.h>

#ifndef BASE
#define BASE 64
#endif
using U = uint16_t;
static constexpr int MAXN = 1024;
alignas(4096) static U pa_[MAXN * MAXN + 8192], pb_[MAXN * MAXN + 8192];
alignas(4096) static U pc_[MAXN * MAXN + 8192], ws_[MAXN * MAXN + 8192];

static U *pa = pa_ + 0, *pb = pb_ + 0, *pc = pc_ + 0, *ws = ws_ + 1024;
static inline __m256i ld(const U *p) {
    return _mm256_load_si256((const __m256i *)p);
}
static inline void st(U *p, __m256i x) {
    _mm256_store_si256((__m256i *)p, x);
}

// The recursive quadrant layout makes every Strassen addition contiguous.
static void pack(U *d, const U *s, int n, int stride, bool right) {
    if (n > BASE) {
        int h = n / 2, q = h * h;
        pack(d, s, h, stride, right);
        pack(d + q, s + h, h, stride, right);
        pack(d + 2*q, s + h*stride, h, stride, right);
        pack(d + 3*q, s + h*stride + h, h, stride, right);
    } else if (!right) {
        for(int i=0;i<n;i+=4) for(int k=0;k<n;k+=2) {
            __m128i a0=_mm_loadl_epi64((const __m128i*)(s+(i+0)*stride+k));
            __m128i a1=_mm_loadl_epi64((const __m128i*)(s+(i+1)*stride+k));
            __m128i a2=_mm_loadl_epi64((const __m128i*)(s+(i+2)*stride+k));
            __m128i a3=_mm_loadl_epi64((const __m128i*)(s+(i+3)*stride+k));
            __m128i t01=_mm_unpacklo_epi32(a0,a1), t23=_mm_unpacklo_epi32(a2,a3);
            _mm_stream_si128((__m128i*)d, _mm_unpacklo_epi64(t01,t23));
            d+=8;
        }
    } else {
        for(int j=0;j<n;) {
            int w=n-j>=24?24:n-j;
            for(int k=0;k<n;k+=2) {
                for(int t=0;t<w;t+=8) {
                    __m128i x=_mm_load_si128((const __m128i*)(s+k*stride+j+t));
                    __m128i y=_mm_load_si128((const __m128i*)(s+(k+1)*stride+j+t));
                    __m128i lo=_mm_unpacklo_epi16(x,y),hi=_mm_unpackhi_epi16(x,y);
                    _mm_stream_si128((__m128i*)d,lo);
                    _mm_stream_si128((__m128i*)(d+8),hi);
                    d+=16;
                }
            }
            j+=w;
        }
    }
}

static void unpack(U *d, const U *s, int n, int stride) {
    if (n > BASE) {
        int h = n/2, q = h*h;
        unpack(d,s,h,stride);
        unpack(d+h,s+q,h,stride);
        unpack(d+h*stride,s+2*q,h,stride);
        unpack(d+h*stride+h,s+3*q,h,stride);
    } else {
        for (int i=0;i<n;++i) {
            U *dp = d + (size_t)i*stride; const U *sp = s + (size_t)i*n;
            int k=0;
            for (;k+32<=n;k+=32) {
                _mm256_stream_si256((__m256i*)(dp+k),    _mm256_load_si256((const __m256i*)(sp+k)));
                _mm256_stream_si256((__m256i*)(dp+k+16), _mm256_load_si256((const __m256i*)(sp+k+16)));
            }
            for (;k<n;k+=16) _mm256_storeu_si256((__m256i*)(dp+k), _mm256_loadu_si256((const __m256i*)(sp+k)));
        }
    }
}

static inline void micro24(const U *a,const U *b,U *c,int n) {
    const U *end=a+4*n; U *c3=c+3*n; intptr_t stride=2*n;
    asm volatile(
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm2\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm4\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm5\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm8\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm10\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm11\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vmovdqa 128(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 288(%[b]), %%ymm12\n\t"
        "vmovdqa 320(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "add $64, %[a]\n\t"
        "add $384, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jae 2f\n\t"
        "1:\n\t"
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vmovdqa 128(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 288(%[b]), %%ymm12\n\t"
        "vmovdqa 320(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "add $64, %[a]\n\t"
        "add $384, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jb 1b\n\t"
        "2:\n\t"
        "vpcmpeqd %%ymm13, %%ymm13, %%ymm13\n\t"
        "vpsrld $16, %%ymm13, %%ymm13\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpand %%ymm13, %%ymm2, %%ymm2\n\t"
        "vpackusdw %%ymm1, %%ymm0, %%ymm0\n\t"
        "vpermq $216, %%ymm0, %%ymm0\n\t"
        "vmovdqu %%ymm0, 0(%[c])\n\t"
        "vpackusdw %%ymm2, %%ymm2, %%ymm2\n\t"
        "vpermq $216, %%ymm2, %%ymm2\n\t"
        "vmovdqu %%xmm2, 32(%[c])\n\t"
        "vpand %%ymm13, %%ymm3, %%ymm3\n\t"
        "vpand %%ymm13, %%ymm4, %%ymm4\n\t"
        "vpand %%ymm13, %%ymm5, %%ymm5\n\t"
        "vpackusdw %%ymm4, %%ymm3, %%ymm3\n\t"
        "vpermq $216, %%ymm3, %%ymm3\n\t"
        "vmovdqu %%ymm3, 0(%[c],%[s],1)\n\t"
        "vpackusdw %%ymm5, %%ymm5, %%ymm5\n\t"
        "vpermq $216, %%ymm5, %%ymm5\n\t"
        "vmovdqu %%xmm5, 32(%[c],%[s],1)\n\t"
        "vpand %%ymm13, %%ymm6, %%ymm6\n\t"
        "vpand %%ymm13, %%ymm7, %%ymm7\n\t"
        "vpand %%ymm13, %%ymm8, %%ymm8\n\t"
        "vpackusdw %%ymm7, %%ymm6, %%ymm6\n\t"
        "vpermq $216, %%ymm6, %%ymm6\n\t"
        "vmovdqu %%ymm6, 0(%[c],%[s],2)\n\t"
        "vpackusdw %%ymm8, %%ymm8, %%ymm8\n\t"
        "vpermq $216, %%ymm8, %%ymm8\n\t"
        "vmovdqu %%xmm8, 32(%[c],%[s],2)\n\t"
        "vpand %%ymm13, %%ymm9, %%ymm9\n\t"
        "vpand %%ymm13, %%ymm10, %%ymm10\n\t"
        "vpand %%ymm13, %%ymm11, %%ymm11\n\t"
        "vpackusdw %%ymm10, %%ymm9, %%ymm9\n\t"
        "vpermq $216, %%ymm9, %%ymm9\n\t"
        "vmovdqu %%ymm9, 0(%[c3])\n\t"
        "vpackusdw %%ymm11, %%ymm11, %%ymm11\n\t"
        "vpermq $216, %%ymm11, %%ymm11\n\t"
        "vmovdqu %%xmm11, 32(%[c3])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b)
        : [c] "r"(c), [c3] "r"(c3), [s] "r"(stride), [end] "r"(end)
        : "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15");
}

static inline void micro16(const U *a,const U *b,U *c,int n) {
    const U *end=a+4*n; U *c3=c+3*n; intptr_t stride=2*n;
    asm volatile(
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm3\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm4\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm5\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vmovdqa 96(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 128(%[b]), %%ymm12\n\t"
        "vmovdqa 160(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "add $64, %[a]\n\t"
        "add $256, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jae 2f\n\t"
        "1:\n\t"
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vmovdqa 96(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 128(%[b]), %%ymm12\n\t"
        "vmovdqa 160(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "add $64, %[a]\n\t"
        "add $256, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jb 1b\n\t"
        "2:\n\t"
        "vpcmpeqd %%ymm13, %%ymm13, %%ymm13\n\t"
        "vpsrld $16, %%ymm13, %%ymm13\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpackusdw %%ymm1, %%ymm0, %%ymm0\n\t"
        "vpermq $216, %%ymm0, %%ymm0\n\t"
        "vmovdqu %%ymm0, 0(%[c])\n\t"
        "vpand %%ymm13, %%ymm2, %%ymm2\n\t"
        "vpand %%ymm13, %%ymm3, %%ymm3\n\t"
        "vpackusdw %%ymm3, %%ymm2, %%ymm2\n\t"
        "vpermq $216, %%ymm2, %%ymm2\n\t"
        "vmovdqu %%ymm2, 0(%[c],%[s],1)\n\t"
        "vpand %%ymm13, %%ymm4, %%ymm4\n\t"
        "vpand %%ymm13, %%ymm5, %%ymm5\n\t"
        "vpackusdw %%ymm5, %%ymm4, %%ymm4\n\t"
        "vpermq $216, %%ymm4, %%ymm4\n\t"
        "vmovdqu %%ymm4, 0(%[c],%[s],2)\n\t"
        "vpand %%ymm13, %%ymm6, %%ymm6\n\t"
        "vpand %%ymm13, %%ymm7, %%ymm7\n\t"
        "vpackusdw %%ymm7, %%ymm6, %%ymm6\n\t"
        "vpermq $216, %%ymm6, %%ymm6\n\t"
        "vmovdqu %%ymm6, 0(%[c3])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b)
        : [c] "r"(c), [c3] "r"(c3), [s] "r"(stride), [end] "r"(end)
        : "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15");
}

static inline void micro8(const U *a,const U *b,U *c,int n) {
    const U *end=a+4*n; U *c3=c+3*n; intptr_t stride=2*n;
    asm volatile(
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm1\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm3\n\t"
        "vmovdqa 32(%[b]), %%ymm12\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "add $64, %[a]\n\t"
        "add $128, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jae 2f\n\t"
        "1:\n\t"
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 32(%[b]), %%ymm12\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "add $64, %[a]\n\t"
        "add $128, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jb 1b\n\t"
        "2:\n\t"
        "vpcmpeqd %%ymm13, %%ymm13, %%ymm13\n\t"
        "vpsrld $16, %%ymm13, %%ymm13\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpackusdw %%ymm0, %%ymm0, %%ymm0\n\t"
        "vpermq $216, %%ymm0, %%ymm0\n\t"
        "vmovdqu %%xmm0, 0(%[c])\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpackusdw %%ymm1, %%ymm1, %%ymm1\n\t"
        "vpermq $216, %%ymm1, %%ymm1\n\t"
        "vmovdqu %%xmm1, 0(%[c],%[s],1)\n\t"
        "vpand %%ymm13, %%ymm2, %%ymm2\n\t"
        "vpackusdw %%ymm2, %%ymm2, %%ymm2\n\t"
        "vpermq $216, %%ymm2, %%ymm2\n\t"
        "vmovdqu %%xmm2, 0(%[c],%[s],2)\n\t"
        "vpand %%ymm13, %%ymm3, %%ymm3\n\t"
        "vpackusdw %%ymm3, %%ymm3, %%ymm3\n\t"
        "vpermq $216, %%ymm3, %%ymm3\n\t"
        "vmovdqu %%xmm3, 0(%[c3])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b)
        : [c] "r"(c), [c3] "r"(c3), [s] "r"(stride), [end] "r"(end)
        : "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15");
}

static void base(const U *a,const U *b,U *c,int n) {
    int j=0;
    for(;j+24<=n;j+=24)
        for(int i=0;i<n;i+=4) micro24(a+i*n,b+j*n,c+i*n+j,n);
    if(n-j==16)
        for(int i=0;i<n;i+=4) micro16(a+i*n,b+j*n,c+i*n+j,n);
    if(n-j==8)
        for(int i=0;i<n;i+=4) micro8(a+i*n,b+j*n,c+i*n+j,n);
}

template<bool Sub>
static inline void combine(U *d,const U *a,const U *b,int len) {
    for(int i=0;i<len;i+=16) {
        __m256i x=ld(a+i),y=ld(b+i);
        st(d+i,Sub?_mm256_sub_epi16(x,y):_mm256_add_epi16(x,y));
    }
}

// Strassen-Winograd: seven products, fifteen additions, two temporaries.
// All operations are in Z/(2^16); truncation never changes the answer.
static void mul_sub(const U *a,const U *b,U *c,int n,U *work,const U *ysub);
static void mul(const U *a,const U *b,U *c,int n,U *work) {
    if(n<=BASE) { base(a,b,c,n); return; }
    int h=n/2,q=h*h;
    const U *a11=a,*a12=a+q,*a21=a+2*q,*a22=a+3*q;
    const U *b11=b,*b12=b+q,*b21=b+2*q,*b22=b+3*q;
    U *c11=c,*c12=c+q,*c21=c+2*q,*c22=c+3*q;
    U *x=work,*y=work+q,*next=work+2*q;
    combine<true>(x,a11,a21,q);
    combine<true>(y,b22,b12,q);
    mul(x,y,c21,h,next);                       // P7
    combine<false>(x,a21,a22,q);
    combine<true>(y,b12,b11,q);
    mul(x,y,c22,h,next);                       // P5
    combine<true>(x,x,a11,q);
    combine<true>(y,b22,y,q);
    mul(x,y,c12,h,next);                       // P6
    combine<true>(x,a12,x,q);
    mul(x,b22,c11,h,next);                     // P3
    combine<true>(y,y,b21,q);
    mul(a22,y,x,h,next);                       // P4
    mul(a11,b11,y,h,next);                     // P1
    for(int i=0;i<q;i+=16) {
        __m256i u2=_mm256_add_epi16(ld(y+i),ld(c12+i));
        __m256i u3=_mm256_add_epi16(u2,ld(c21+i));
        __m256i p5=ld(c22+i);
        st(c12+i,_mm256_add_epi16(_mm256_add_epi16(u2,p5),ld(c11+i)));
        st(c21+i,_mm256_sub_epi16(u3,ld(x+i)));
        st(c22+i,_mm256_add_epi16(u3,p5));
    }
    /* c9 折叠:`c11(P2 象限) += y`(3 次访存/元素)折进 P2 子乘法。
       mul_sub 由 mul 的胶水部分构成、**不含 base()/微内核**(否则 gcc 9.3 会递归克隆
       ⇒ asm 2548→13183 行、vpmaddwd 192→3072,正是 fs3_ 实测 +15.3% 的那个坑)。 */
    if (h>BASE) mul_sub(a12,b21,c11,h,next,y);
    else { mul(a12,b21,c11,h,next); combine<false>(c11,c11,y,q); }
}

static void mul_sub(const U *a,const U *b,U *c,int n,U *work,const U *ysub) {
    int h=n/2,q=h*h;
    const U *a11=a,*a12=a+q,*a21=a+2*q,*a22=a+3*q;
    const U *b11=b,*b12=b+q,*b21=b+2*q,*b22=b+3*q;
    U *c11=c,*c12=c+q,*c21=c+2*q,*c22=c+3*q;
    U *x=work,*y=work+q,*next=work+2*q;
    combine<true>(x,a11,a21,q);
    combine<true>(y,b22,b12,q);
    mul(x,y,c21,h,next);                       // P7
    combine<false>(x,a21,a22,q);
    combine<true>(y,b12,b11,q);
    mul(x,y,c22,h,next);                       // P5
    combine<true>(x,x,a11,q);
    combine<true>(y,b22,y,q);
    mul(x,y,c12,h,next);                       // P6
    combine<true>(x,a12,x,q);
    mul(x,b22,c11,h,next);                     // P3
    combine<true>(y,y,b21,q);
    mul(a22,y,x,h,next);                       // P4
    mul(a11,b11,y,h,next);                     // P1
    {
        const U *s12=ysub+q,*s21=ysub+2*q,*s22=ysub+3*q;
        for(int i=0;i<q;i+=16) {
            __m256i u2=_mm256_add_epi16(ld(y+i),ld(c12+i));
            __m256i u3=_mm256_add_epi16(u2,ld(c21+i));
            __m256i p5=ld(c22+i);
            st(c12+i,_mm256_add_epi16(_mm256_add_epi16(_mm256_add_epi16(u2,p5),ld(c11+i)),ld(s12+i)));
            st(c21+i,_mm256_add_epi16(_mm256_sub_epi16(u3,ld(x+i)),ld(s21+i)));
            st(c22+i,_mm256_add_epi16(_mm256_add_epi16(u3,p5),ld(s22+i)));
        }
    }
    mul(a12,b21,c11,h,next);                   // P2
    for(int i=0;i<q;i+=16)
        st(c11+i,_mm256_add_epi16(_mm256_add_epi16(ld(c11+i),ld(y+i)),ld(ysub+i)));
}
void matrix_multiply(int n,const short *A,const short *B,short *C) {
    pack(pa,(const U *)A,n,n,false);
    pack(pb,(const U *)B,n,n,true);
    _mm_sfence();   /* pack 现在用非临时存写面板:递归读它们之前必须 sfence */
    mul(pa,pb,pc,n,(U*)C + 1024);
    unpack((U *)C,pc,n,n);
    _mm_sfence();
}

CompilationN/AN/ACompile OKScore: N/A

Testcase #111.984 ms8 MB + 12 KBAcceptedScore: 100


Judge Duck Online | 评测鸭在线
Server Time: 2026-11-09 03:23:33 | Loaded in 1 ms | Server Status
个人娱乐项目,仅供学习交流使用 | 捐赠