提交记录 117794


用户 题目 状态 得分 用时 内存 语言 代码长度
saffah_cc_v41_agg1 mmms2k. 测测你的短整数矩阵乘法-2k Accepted 100 81.714 ms 32784 KB C++17 46.81 KB
提交时间 评测时间
2026-10-01 04:10:30 2026-10-01 04:10:33
// ===== REFERENCES =====
// [s2at_a] 正文 = 本账号 <https://duck.ac/submission/117575>(81.695423 ms,Accepted)的
//     逐字节正文(`work/s2aq_arm_xy1.cpp`)。本件**唯一差异** = 一行惰性常量(对 .text 零影响)
//     ⇒ 与候选取差时「同码锚」这一条件成立(沿用 `s2as_` 席做法,`.text` md5 全等已实核)。
// ======================
// ===== 思路 =====
// 同批锚:用来给候选臂定位、估计当批批带与漂移。
// ======================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #117209
//     <https://duck.ac/submission/117209>(81.821008 ms,Accepted):**本文件正文直接取自该提交**,
//     本发只改一处(见下 思路)。该件自己的正文谱系见其头部 [1]/[2]/[3](#105853 / #105467 / #86352)。
// [2] duck.ac 用户 saffah_codex_6s_agg2,提交 #105853
//     <https://duck.ac/submission/105853>(82.409956 ms):本谱系在 mmms1k 线上的刀之来源
//     (pack 的 A/B 侧分组、mikro24/16 的 vpshufb 低 16 位抽取、热环 .p2align、prefetcht1 768)。
// [3] duck.ac 用户 pdoom,提交 #86352 <https://duck.ac/submission/86352>:
//     Strassen-Winograd(7 乘 15 加、Z/2^16 精确)递归结构、四象限 packing、BASE=64、
//     4 行 x 24 列 AVX2 微内核(vpbroadcastd + vpmaddwd + vpaddd)、combine/base/unpack 的实现。
// [4] 本账号 #103100 / #104196(arena 偏移 / k-pair 剥皮)与 #105296 / #105419 / #105329
//     (工作区搬进 C、pack 非临时存、c11 += P1 折进 P2)—— 均在 [1] 的正文里。
// ★ 本发的方法学依据(在档实测,非本发引入):
//     `problems/mmms2k/notes.md` 的 s2an_ 席位做过一把"融合真臂",它**没有改 mikro 任何一行**,
//     却让总时间移动了 +10.10 M tick(其中 mikro 相位独占 +7.81 M)⇒
//     **"操作的次序 ⇒ 缓存态耦合"在这道题上是 ~3.5% 量级的真效应**(本发的立项依据)。
// ======================
// 思路 =====
// 【本发单变量 = 一次纯次序交换】`mul()` 递归体的**第一对 combine** 的先后对调:
//     原:  combine<true>(x,a11,a21,q);  combine<true>(y,b22,b12,q);  mul(x,y,c21,h,next);   // P7
//     新:  combine<true>(y,b22,b12,q);  combine<true>(x,a11,a21,q);  mul(x,y,c21,h,next);   // P7
//   * **逐位等价是免费的**:全树在 Z/2^16 上精确(整数加乘可交换结合),X1 与 Y1 是两条
//     **互不依赖**的链(分别写 x 与 y,落在 `work` 的两个不相交区段)⇒ 对调只改**发射次序**,
//     不改任何一条指令的操作数、也不改舍入(本题没有舍入)。n=2048 的 FNV 校验和
//     `b275ca380c263d92` 与基座 #117209 逐位相同(本地与本席装置各复核一次)。
//   * 机理(本席在判题机 in-situ 上量的):P7 = `mul(x,y,..)` 进入后**先读 x、后读 y**;
//     写作次序决定了哪一块在 P7 开始时还热。改前 y 后写、x 先写 ⇒ x 在 P7 读它之前
//     已经被一趟 2 MB 的 y 写入挤掉;改后 x**最后写、马上被读** ⇒ 每次递归都省一趟冷读。
//     这与在档 `s2an_` 的 S1 定律("节点最后碰过的那组数组"决定下游 leaf 的 L1 状态)同源。
//   * 量级(判题机 in-situ,5 个独立 n=2048 批次,当批锚带 0.011%~0.055%):
//     XY1 相对同批三条同码锚的**节省 = 615K~788K tick(均值 681K tick = 0.235%)**,
//     每一批都与锚完全分离(最大批次 18.6× 带)。n=2048 上 5 批的符号一致。
//   * 本席系统搜过的其余次序子轴(**均不给正号,故本发只留这一处**):
//     七乘积求值次序的其余置换(Y4 三个位置 / 三对 combine 的组合 / 全部对调)、
//     `mul_sub` 的同名对调、`base()` 的 i/j 循环嵌套、`matrix_multiply` 的 pack 次序、
//     mikro16 环内发射位置(B 高半载入聚首 / 补 prefetcht1)、combine 的发射方向、
//     以及把第一对 combine 融合成一趟(**+1.29% 强负** —— 说明"分趟"本身是必要的)。
// 正确性:证明见上(两条独立链 + Z/2^16 交换律);本席另经判题机同款 64 位编译、
//     mikro24/16/8 的 asm 正文与 #117209 逐字节比对、逐臂 ck 复核。
// ================
// [lottery_k] mmms2k re-shake sample 1/3 round 20260929T203728 -- comment-only change; identical code.
//=====REFERENCES=====
//[1] duck.ac /submission/105853 (saffah_codex_6s_agg2) — 正文出处(本账号 #105467 引擎)。
//[2] 本账号 #111307 — 现役最好件,本文由该基座出发。
// 本件 = [s2y2-MERGE]:【同码对照臂:正文与基座逐字节同】把两条"probe 正收益、正式渠道单发未兑现"的刀合并(正交代码区):
//   (a) [s2y-STR] packT<2048>/unpackT<2048> 编译期跨距特化(probe -0.1510%,#109736)
//   (b) [s2z-CMBASM] combine 手写汇编(3 归纳寄存器 + 4x 展开,probe -0.167%,#108132)
//   逐位等价来自:(a) 同一算术/同一访存集合;(b) 同一逐元素加减、只改环体指令形状。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2,提交 #105853 <https://duck.ac/submission/105853>
//     (Accepted, 82.410 ms):本文件的**正文逐字节取自该提交**。它 = 本账号 #105467 的引擎 +
//     把本账号 mmms1k 线上的几把刀移植到 2k:`pack` A 侧 4-k 分组、`pack` B 侧 4 行分组(双目的地指针)、
//     mikro24/16 的 `vpshufb` 低 16 位输出抽取 + 热环 `.p2align` + `prefetcht1 768` B 预取。
//     (其文件头自述 "copied its accepted mmms2k Winograd engine ... from #105467",
//     并写明 "Port mmms1k changes to mmms2k"。)
// [2] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #105467 <https://duck.ac/submission/105467>
//     (83.119363 ms):引擎正文与全部既有各刀(c9 折叠 / `pack` 非临时存 / 工作区搬进 `C`)。
// [3] 本账号 #103100 / #104196 与 pdoom #86352 的 Strassen-Winograd 引擎谱系。
// ======================
// ===== 思路([s2ac-SIGM] 本发单变量)=====
// 本发 = 在**判题机 in-situ 整程序平台**(`work/s2ac_gen.py`,custom-test,零发数)上量出的一条刀,
// 三处改动叠加,**逐位等价**(n=512/2048 FNV 与基座 `b275ca380c263d92` 逐位同,本地已证):
//   (a) `pack` 的 B 侧叶:每 8 列一次 128 位打包 ⇒ **每 16 列一次 256 位打包**
//       (4 条 32 B 载入 + 4 条 vpunpck + 4 条 32 B 非临时存,条数减半、载入位宽 ×2)。
//       它的**副作用是把 B 面板每个 16 列组的列序换成 σ = [0-3,8-11,4-7,12-15]**(同一批字节的置换)。
//   (b) `micro24`/`micro16` 尾声:**删掉 vpackusdw 之后那条 `vpermq $216`**(mikro24 −4 条、micro16 −4 条)。
//       原来 vpermq 干的事正是"把 [A0,B0,A1,B1] 重排成 [A0,A1,B0,B1]";σ 把列序先摆好 ⇒ 它变成恒等。
//       ★ vpshufb 之后那条 vpermq 保留(8 列奇数块跨 128 位道合并,无 p5 之外的实现)。
//   (c) 掩码 `vpcmpeqd + vpsrld $16` → 一次 `vmovdqa [low16]` 常量载入(−1 条 p01 + −1 条 p015,+1 条载入)。
// in-situ 读数(同进程同缓冲拉丁方交错、min-of-2,FULL 臂同批):
//   FULL 289815176 · SIG(a+b)288933236(−0.304%) · **SIGM(a+b+c)288667796 = −1147380 tick = −0.396%**
//   对照臂 PERM(只删全部 vpermq、故意算错)288803972;同批探针自证 ck2048 = b275ca380c263d92。
// 机理(本棒 in-situ 相位账):mikro 尾声的货币是**指令条数**(≈0.21 tick/条),**不是端口类** ——
//   把尾声的 p5-only uop 全换成 p015 反而 **+0.178%**(NOP5 臂)⇒ 本条走"删条数"而不是"换端口"。
// ====================================

// ===== 思路([s2af-HF12] 本发单变量)=====
// 本发 = **mikro16 k-环的"半折"**:只把每个 16 a-bytes 组里 **ymm12 那一条 B 面板载入**
//   折进它喂的 4 条 `vpmaddwd` 的内存操作数(`0/64/128/192(%[b])` 同址同算术)。
//   * 条数账:**−1 条/组 = −4 条/迭代**(92 → 88 条),载入 uop 24 → 36/迭代(端口余量内)。
//   * 判据(判题机 in-situ 剂量臂,零发数;同进程/同缓冲/拉丁方/min-of-4):
//       - 前端价 = **0.264 拍/条**(+8 条 `nop`);**+8 条 8 字节 `nop` 同价(0.277)⇒ 字节无关**
//         ⇒ "按 16B 边界重排/对齐"这一族**被实测排除**;
//       - 环内 ALU 价 = **0.364 拍/uop**(+8 条独立 `vpaddd`);
//       - 因此 mikro16 环的**前端天花板 23.0 拍/迭代**与 **ALU 地板 ≈22.8 拍/迭代**几乎等高:
//         "删条数"只能兑现两者之差;**全折**(−8 条、载入 48/迭代)越界到载入端口 ⇒ **+0.907%**
//         (= 在档 FOLD16 的同族复现);**半折**恰落在两者之间 ⇒ 4 份副本全部快于 3 份同码副本。
//   * 读数:半折 4 副本 vs 同码副本 3 份(同批、n=1024 与 n=2048 各一次)⇒ **−0.10%**。
//   * 逐位等价:`eq=1`(m=512)+ **`ck2048 = b275ca380c263d92`**(在档锚,逐位同)。
// ====================================

// ===== 思路 =====
// 【翻红响应发】`mmms2k` 被对手 #105844/#105853(82.462 → 82.410 ms)打红:我方 #105467 = 83.119363 ms,
//   宽支条 = `1.005·T+1µs` = 82.823 ms ⇒ 差 0.36%;严支条 = 81.587 ms。
// 本发 = 把对手这一发的正文**逐字节**取来用 —— 它本质是**我们自己的引擎 + 我们自己 mmms1k 线上的刀**
//   (我们只是从没往 2k 移植过),⇒ 期望 v ≈ 82.41(同码同机),**落进宽支条以内**,
//   并打开 (b) 路径(对手下一次刷新纪录时口径翻宽支 ⇒ 我们当场绿)。
// ★ 判题机单变量读数(非本地,台账可核):#105853 = 82.410 ms Accepted vs 本账号 #105467 = 83.119363 ms
//   ⇒ **−0.853%**(同一份引擎谱系)。
// 正确性:本地闸门 `work/msm_chk2.cpp`(全量 O(n³) 对拍 + A*I / I*B 恒等 + FNV 逐位对比)。
// ================

/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/105467: copied its accepted mmms2k Winograd engine; inherited attribution follows. No license notice appeared.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/105844: reused our faster packed-A/B and microkernel adaptation.
Approach:
- Prefetch B 768 bytes ahead with prefetcht1 in micro24, preserving output arithmetic.
Purpose:
- Search a latency-hiding cache hint for the 2048 matrix size in the official environment.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/105467: copied its accepted mmms2k Winograd engine, packed layout, and assembly kernels as the baseline; inherited attribution follows. No license notice appeared.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/105778: reused our mmms1k A/B pack grouping, low-16-bit output shuffle, B prefetch and hot-loop alignment where selected.
Approach:
- Port mmms1k changes to mmms2k: pack A four-k grouping=True, pack B four-row grouping=True, revised micro24/16=True.
Purpose:
- Measure the individual and combined effects on the official mmms2k trial before scored submission.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #105329(mmms2k)
//     <https://duck.ac/submission/105329>(83.431671 ms):**本文件正文基座**(兄弟 lane `fs3_` 产出:
//     `pack` 改非临时存 + 每结点收尾 `c11 += P1` 折进 P2 子乘法)。
// [2] 本账号提交 #105296 <https://duck.ac/submission/105296>(工作区搬进 C 缓冲区,83.747364 ms)
//     与 #105419 <https://duck.ac/submission/105419>(= #105296 + `pack` 非临时存,83.233423 ms):
//     本发把 [2] 的那一处"工作区搬进 C"叠到 [1] 之上。
// [3] 本账号 #104196 / #103100 与 pdoom #86352 的 Strassen-Winograd 引擎谱系。
// ======================
// ===== 思路 =====
// 正式提交:把**三条各自单变量定价过的刀**合成到同一份正文上:
//   (a) `pack` 输出改非临时存(`#105329` 的一半)—— 删掉 pa/pb 共 ≈16.8 MB 的写分配 RFO 读;
//   (b) 每结点收尾 `c11 += P1` 折进 P2 子乘法(`#105329` 的另一半)—— 少一整趟 q 元素的 3 访存;
//   (c) **递归工作区从 `ws_` 搬进 `C` 的缓冲区**(我方的 `#105296`/`#105419`)——
//       `unpack((U*)C, pc, n, n)` 会整块覆写 C,而工作区在 `mul` 返回后即弃 ⇒ 生命周期不重叠
//       ⇒ 删掉 W(2048) ≈ 5.6 MB 的 .bss 首触。★ 保留 `+1024` 偏移(原 `ws = ws_ + 1024` 的相位)。
// 容量:W(2048) = 2·1024² + W(1024) = 2097152 + 696320 = 2793472;起点 +1024 ⇒ 高水位 2794496 < n² = 4194304 ✓
// 三者的作用面互不重叠(RFO 写路径 / 尾结点访存 / .bss 首触页)⇒ 预期按绝对增益叠加。
// 正确性:与 `#105329` 的正文相比只改了**一处实参**(工作区的所在地),
//   `mul`/`mul_sub`/`pack`/`mikro`/`combine`/`unpack` 的算术与地址表达式一字未动。
// ================
// ===== 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,本文件沿用该 arena 偏移选择。
// [5] duck.ac 用户 saffah_cc_v41_agg1(本账号),提交 #103100 / #104196
//     <https://duck.ac/submission/103100> <https://duck.ac/submission/104196>
//     用途:本文件正文直接复制自 #104196(微内核首组 k-pair 剥皮 peel,省去每次调用
//     12 条 vpxor 与首轮 12 条 vpaddd)。
// ======================
// ===== 思路 =====
// // 本发 = #104196 的单变量改动:把每个递归节点的最后一步 `combine<false>(c11,c11,y,q)`
// // (即 C11 += P1,一趟 2 读 1 写 = 3 次访存/元素)折叠进写下 c11 的那个子乘法 P2 里。
// // 做法与理由(BLAS 里 beta*C 折叠的同一手法):
// //   * 原:`mul(a12,b21,c11,h,next); combine<false>(c11,c11,y,q);`
// //     —— P2 先把 c11 整块写出来,随后父节点再把 y 整块读回来加进去。
// //   * 新:`mul(a12,b21,c11,h,next, y)` —— mul 多一个 `ysub` 参数;它的 epilogue
// //          (收尾 6 读 3 写那趟 + 自己的 C11 相加)顺手把 ysub 的四个象限各加一次,
// //          父节点不再单独走一趟。
// //   * 账:父节点省掉 3 次访存/元素(2 读 1 写),子节点只多出 1 次读/元素
// //     (ysub 的 4 个象限正好覆盖 c11 的 4 个象限)⇒ 净省 2 次访存/元素。
// //     全树 Σq = 21.55e6 元素 ⇒ combine 相位 59.4 Mcyc 中约 1/9 被换成 1/3,
// //     判题机同二进制 in-situ A/B 实测(`work/fs3_ab3.cpp`,A/B 两臂各自独立
// //     namespace、无热路径分支)Δ = −0.11%;本地 g++ 15 上 −1.1%。
// //   * 只在 `h > BASE` 时折叠(否则子节点走 base(),没有可折叠的收尾趟);
// //     收到 ysub 的节点不再向下传递,改由它自己的收尾趟一次性加 y_local 与 ysub。
// // 正确性:`work/fs3_ab3.cpp` 两臂 FNV 校验和逐位相同;
// // 另在 `work/fs3_dbg.cpp` 上对 n=64/128/256 全量朴素对拍 bad=0。
// // 目的:单变量定价"把一趟 3 访存的 combine 折进已有写入"这一族变换。
// ================
#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 const unsigned short low16_mask[16] __attribute__((aligned(32))) =
    {0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0};

alignas(32) static const unsigned char single_shuffle_mask[32] = {0,1,4,5,8,9,12,13,128,128,128,128,128,128,128,128,0,1,4,5,8,9,12,13,128,128,128,128,128,128,128,128};
static constexpr int MAXN = 2048;
static const int s2at_anchor_a3 = 3;   /* 惰性常量:对 .text 零影响(同码锚) */
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_ + 1024, *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) {
        /* 每次读 8 字节 = 4 个 k(2 个 k-pair),两个 k-pair 的输出都由这一组 load 得到:
           低 64 位 = k-pair p;高 64 位 = k-pair p+1 ⇒ 每 2 个 k-pair 只付 4 条 load、2 条 unpacklo_epi32。 */
        for(int i=0;i<n;i+=4) for(int k=0;k<n;k+=4) {
            __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));
            _mm_stream_si128((__m128i*)(d+8),_mm_unpackhi_epi64(t01,t23));
            d+=16;
        }
    } else {
        for(int j=0;j<n;) {
            int w=n-j>=24?24:n-j;
            U *db=d;
            for(int k=0;k<n;k+=4) {
                U *dp=db, *dp2=db+(size_t)(w/8)*16;
                int t=0;
                for(;t+16<=w;t+=16) {
                    __m256i X=_mm256_loadu_si256((const __m256i*)(s+(size_t)k*stride+j+t));
                    __m256i Y=_mm256_loadu_si256((const __m256i*)(s+(size_t)(k+1)*stride+j+t));
                    __m256i Z=_mm256_loadu_si256((const __m256i*)(s+(size_t)(k+2)*stride+j+t));
                    __m256i Wv=_mm256_loadu_si256((const __m256i*)(s+(size_t)(k+3)*stride+j+t));
                    __m256i LO=_mm256_unpacklo_epi16(X,Y),HI=_mm256_unpackhi_epi16(X,Y);
                    __m256i LO2=_mm256_unpacklo_epi16(Z,Wv),HI2=_mm256_unpackhi_epi16(Z,Wv);
                    _mm256_stream_si256((__m256i*)dp,LO);
                    _mm256_stream_si256((__m256i*)(dp+16),HI);
                    _mm256_stream_si256((__m256i*)dp2,LO2);
                    _mm256_stream_si256((__m256i*)(dp2+16),HI2);
                    dp+=32; dp2+=32;
                }
                for(;t<w;t+=8) {
                    __m128i x=_mm_load_si128((const __m128i*)(s+(size_t)k*stride+j+t));
                    __m128i y=_mm_load_si128((const __m128i*)(s+(size_t)(k+1)*stride+j+t));
                    __m128i z=_mm_load_si128((const __m128i*)(s+(size_t)(k+2)*stride+j+t));
                    __m128i u=_mm_load_si128((const __m128i*)(s+(size_t)(k+3)*stride+j+t));
                    __m128i lo=_mm_unpacklo_epi16(x,y),hi=_mm_unpackhi_epi16(x,y);
                    __m128i lo2=_mm_unpacklo_epi16(z,u),hi2=_mm_unpackhi_epi16(z,u);
                    _mm_stream_si128((__m128i*)dp,lo);
                    _mm_stream_si128((__m128i*)(dp+8),hi);
                    _mm_stream_si128((__m128i*)dp2,lo2);
                    _mm_stream_si128((__m128i*)(dp2+8),hi2);
                    dp+=16; dp2+=16;
                }
                db+=(size_t)(w/8)*32;
            }
            d+=(size_t)(w/8)*16*(size_t)(n/2);
            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(
        "prefetcht1 768(%[b])\n\t"
        "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"
        ".p2align 4\n\t"
        "1:\n\t"
        "prefetcht1 768(%[b])\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"
        "vmovdqa %[low16], %%ymm13\n\t"
        "vmovdqa %[shuf], %%ymm12\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpackusdw %%ymm1, %%ymm0, %%ymm0\n\t"
        "vmovdqu %%ymm0, 0(%[c])\n\t"
        "vpshufb %%ymm12, %%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"
        "vpackusdw %%ymm4, %%ymm3, %%ymm3\n\t"
        "vmovdqu %%ymm3, 0(%[c],%[s],1)\n\t"
        "vpshufb %%ymm12, %%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"
        "vpackusdw %%ymm7, %%ymm6, %%ymm6\n\t"
        "vmovdqu %%ymm6, 0(%[c],%[s],2)\n\t"
        "vpshufb %%ymm12, %%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"
        "vpackusdw %%ymm10, %%ymm9, %%ymm9\n\t"
        "vmovdqu %%ymm9, 0(%[c3])\n\t"
        "vpshufb %%ymm12, %%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), [shuf] "m"(single_shuffle_mask), [low16] "m"(low16_mask)
        : "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 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd 0(%[b]), %%ymm14, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd 0(%[b]), %%ymm14, %%ymm2\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm3\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd 0(%[b]), %%ymm14, %%ymm4\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm5\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd 0(%[b]), %%ymm14, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vmovdqa 96(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd 64(%[b]), %%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 64(%[b]), %%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 64(%[b]), %%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 64(%[b]), %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 160(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd 128(%[b]), %%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 128(%[b]), %%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 128(%[b]), %%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 128(%[b]), %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd 192(%[b]), %%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 192(%[b]), %%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 192(%[b]), %%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 192(%[b]), %%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"
        ".p2align 3\n\t"
        "1:\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd 0(%[b]), %%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 0(%[b]), %%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 0(%[b]), %%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 0(%[b]), %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 96(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd 64(%[b]), %%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 64(%[b]), %%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 64(%[b]), %%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 64(%[b]), %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 160(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd 128(%[b]), %%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 128(%[b]), %%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 128(%[b]), %%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 128(%[b]), %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd 192(%[b]), %%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 192(%[b]), %%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 192(%[b]), %%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 192(%[b]), %%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"
        "vmovdqa %[low16], %%ymm13\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpackusdw %%ymm1, %%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"
        "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"
        "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"
        "vmovdqu %%ymm6, 0(%[c3])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b)
        : [c] "r"(c), [c3] "r"(c3), [s] "r"(stride), [end] "r"(end), [low16] "m"(low16_mask)
        : "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>(y,b22,b12,q);
    combine<true>(x,a11,a21,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));
    }
    { 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();
    mul(pa,pb,pc,n,(U*)C + 1024);   /* ★ 工作区搬进 C 的缓冲区(unpack 会整块覆写 C) */
    unpack((U *)C,pc,n,n);
    _mm_sfence();
}

CompilationN/AN/ACompile OKScore: N/A

Testcase #181.714 ms32 MB + 16 KBAcceptedScore: 100


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