// ===== REFERENCES =====
// [R1] 本账号 saffah_cc_v41_agg1 的现役最好件 #122705
// <https://duck.ac/submission/122705>（75.015667 ms，本题排行榜 rank 1）
// = 本地源文件 `problems/mmms2k/work/ms2af__tmp/ms2af_anchor.cpp`。
// **本件 = 该件正文逐字节 + 本席下面「思路」段的两处机械变换**；引擎代数、
// 指令序列、内存布局相对关系一律未改。
// [R2] 该件内已逐字保留的上游谱系（pdoom #118374/#118809/#118817/#86352/#112312/
// #117569/#118412、saffah_codex_6s_agg2 #110387/#102935、saffah_cc_v41_260924
// #96771、s4de_/s4df_ #119672/#119711 等，见下方原 REFERENCES 段）继续有效，
// 未改一字。
// [R3] 判题机工具链取 `/home/yjp/duck.ac/ref/gcc9/g9.sh`（gcc 9.3，与判题一致）。
// 合规：duck.ac 提交正文按站点规则公开可见；原提交无独立许可证声明，按 RULES.md 第 3 节
// 署名作者账号与原地址。
// ===== 思路 =====
// 【本发目的】实验性提交：把两处**纯代码生成**（不改任何字节语义）的改动叠到 #122705 上，
// 用判题机实测它们是否为正增益。两处改动都过**双驱动逐位等价闸门**：
// · ms2t_drv（线性）：chk dd4d25c4b69d0383（与 #122705 同）
// · ms2ac_rnd（splitmix64 全域随机）：chk 5ffc902315655516，且 c0/cN 逐字同
//
// 【变量 1（单变量）】`raw_basis_low32` / `raw_basis_low16` 的 32B/16B 载入由
// `_mm256_loadu_si256` / `_mm_loadu_si128` 改成 `_mm256_load_si256` / `_mm_load_si128`。
// 理由：gcc9 在本题档位默认开 `-mavx256-split-unaligned-load` ⇒ 每条 32B 未对齐载入被拆成
// `vmovdqu xmm` + `vinserti128 m128` **两条指令**（`objdump -d` 计数：`_ZL9pack_wide`
// 374 条里 76 组、`_ZL11pack_wide_b` 272 条里 59 组）。
// **对齐可证**：raw_basis* 的全部调用点都是
// pack(): s + i*stride + k （stride = 2048 短字 = 4096 B；k 为 16 的倍数）
// pack(): s + k*stride + j + t （j 为 32 的倍数、t 为 16 的倍数）
// pack_wide: sr + i*S4ST + K （S4ST = 2048；K ∈ {0,16,...,112}）
// pack_wide_b: sr + k*ST + J （ST = 2048；J ∈ {0,16,...,112}）
// raw_basis_low* 另读 p-32*2048 / p-32；raw_basis32/16 另读 p-64*2048 / p-64
// 而 s 恒为「4096 B 对齐的矩阵基址 + h / h*stride（h 为 16 的倍数）」⇒ 全部地址 32B
// （resp. 16B）对齐。改动只变发射的指令，读到的字节集合完全相同。
//
// 【变量 2（单变量）】`native2x32_cache3` 叶核内联汇编的循环标签 `3:` 前加一条 `.p2align 6`
// （叶核循环体实测 2037 B / 358 条指令。**单独测该变量 = +0.0003 ms ≈ 0**，故本件里它只是
// 布局微调，不是收益来源；两变量已分别单独定价）。
//
// 【判题机定价（tryrun，同批轮转，锚 = #122705 正文）】
// 批次 A：锚 80.547/80.533；仅变量 1 = 80.456/80.442 ⇒ **−0.091 ms（−0.113%）**
// 批次 B：锚 80.547/80.556；变量 1 = 80.456/80.454 ⇒ −0.091（−0.113%）；
// 变量 1+2 = 80.438/80.452 ⇒ **−0.109 ms（−0.135%）**
// 批次 C：锚 80.526/80.559；变量 1+2 = 80.452/80.457 ⇒ −0.074（−0.092%）
// ⇒ 三次独立批次同号（−0.074 … −0.109 ms），方向一致；本发即取最好那一档「变量 1+2」。
// 机理：拆载把 1 条载入变成 2 条（多出一条 p5 的 vinserti128），pack 阶段是内存/端口受限的。
//
// 【逐位等价论证】两个变量都不改变任何被读/写的字节或运算：变量 1 只是同地址同宽度的对齐
// 载入（地址已证 32B 对齐）；变量 2 只插入 nop 填充。双驱动 chk 逐字相同（见上）。
// ==================
// ===== REFERENCES =====
// [1] 本账号（saffah_cc_v41_agg1）在**姊妹题 mmms4k** 的现役最好件
// <https://duck.ac/submission/120060>（529.100997 ms，Accepted）= 源文件
// `problems/mmms4k/work/s4dy_tmp/s4dy_sub2_A.cpp`（md5 220c9648328f18a51c7d4e37e116f207）。
// **本件正文 = 该提交正文逐字节 + 机械维度特化**（见 思路）；引擎代码、指令序列、代数一律未改。
// [2] 该件的完整谱系（正文内 Credit / REFERENCES 段**逐字保留、未改一字**）：duck.ac 用户 **pdoom**
// #118374/#118809/#118817（BASE32 精确引擎、4 行并发乘积调度、coarse64 根打包、NT 流）·
// #86352（mod65536 引擎、pair-dot）· #112312（根打包融合、退役根槽复用）· #117569（SIMD 原生产出、
// 延迟根解码）· #118412（底层 alternative basis 两端融合）；本账号 #118613（ROOT_PAD 与工作区相位）·
// #118236（three-madd 调度、四象限根输出融合、C 根原始输入复用）· #112006/#117247（peeled pair、
// 操作数/输出打包、cache phase、三行调度）；**saffah_codex_6s_agg2** #110387/#102935（融合重构、
// 常量递归与指令顺序）；**saffah_cc_v41_260924** #96771（cache phase）；**s4de_/s4df_** 的
// pack 左/右葉広幅読序（#119672 / #119711）。
// Oded Schwartz & Noa Vaknin, SIAM J. Sci. Comput. (2023), doi 10.1137/22M1502719
// （exact 七乘十二加 alternative-basis 分解）—— 仅参考其思想，未逐字移植。
// [3] 移植靶的（**本席原创**的机械变换）：把上方 [1] 的 4096 引擎按 `n` 的每一处显式出现改成 2048。
// 变换清单 = 生成器 `problems/mmms2k/work/ms2x_gen.py` 里逐条列出的替换（每条带命中计数断言），
// 本件由该生成器从 [1] 的源文件确定性地产出（无手改）。
// 合规：duck.ac 提交正文按站点规则公开可见；原提交无独立许可证声明，按 RULES.md 第 3 节署名作者账号与原地址。
// ======================
// ===== 思路 =====
// 【靶子】本题（mmms2k, n=2048）我方现役最好件是**另一条技术线**（BASE=64、mikro24/16 微内核）的
// 81.479064 ms；而本账号在**同族 4k** 上的引擎（BASE=32、N=64/128 alternative-basis、递归 Strassen）
// 从未移植过来。4k 引擎在 n=4096 上跑 529.101 ms ⇒ 按 Strassen 的 7/8-每层换算，n=2048 应落在
// 约 75.6 ms 量级，**恰是本题当前 T=75.138925 的量级** —— 说明"同族 4k 引擎未移植"就是本题
// 8.7% 缺口的来源（对手 #120454 已用**我方的 #120060** 做了同一件事并据此刷出 T）。
// 【本发单变量 = 维度特化本身】`s4dy_sub2_A.cpp` 的引擎把 n 写死在若干处显式常量里：
// MAXN(缓冲区)·根打包 stride `s+h*4096`/`s+64*4096`·根槽尺寸 RQ=2048²·`pack_root(...,2048,...)`·
// `mul<2048>`·`root_fuse<2048>(...,stride=4096)`·`combine3_unpack<2048>`/`unpack_add<2048>`（`N==4096` 分支）·
// 叶 pack 的 `S4ST`/`ST` 行距·`transformed` 旗标的位测试 `(32*4096)`。**全部按 n 的定义改成 2048 版**
// （4096→2048、2048→1024、bit17→bit16，因为 32*2048=2^16）；同时**删掉父件里那段与矩阵算术无关的
// 死区诊断**（`s4dy_pad2` 声明 + `matrix_multiply` 末尾的 volatile 逐页写循环）。
// ★ 数值定义：`32*stride` 的位测试从 bit17 变 bit16 是**同一语义**（stride 4096→2048），不是新算法；
// 其余全部是"同一个几何量换了个数"。**代数、指令序列、内存布局相对关系一字未改。**
// 【为什么必须这样改而不是整份重写】该引擎的正确性依赖递归象限布局的自洽（`q=h*h+(h/BASE)²*TILE_PAD`
// 在 src 与 packed 两侧一致），手写新维度极易破坏它；机械替换把风险压到"替换是否完备"这一件事上，
// 而完备性由生成器逐条计数断言 + 落盘后的"残留 4096 行"清单双重检查。
// 【逐位等价闸门】n=2048 随机输入下与本账号现役 AC 引擎 `work/s2at_arm_dy512.cpp`（#117816）
// 的全 C 逐字节 `cmp`（见 `work/ms2x_gate.py`）。
// ======================
#define ALTMAX 128
#define ROOT_PAD 160
#pragma GCC optimize("O3,unroll-loops,web,rename-registers")
#pragma GCC target("avx2")
#include <immintrin.h>
#include <stdint.h>
#include <string.h>
asm(".macro vzeroupper\n.endm\n");
#ifndef TILE_PAD
#define TILE_PAD 16
#endif
#ifndef BASE
#define BASE (1 << 5)
#endif
using U = uint16_t;
static constexpr int MAXN = 2048;
alignas(4096) static U pa_[MAXN * MAXN + (MAXN/BASE)*(MAXN/BASE)*TILE_PAD + 8192], pb_[MAXN * MAXN + (MAXN/BASE)*(MAXN/BASE)*TILE_PAD + 8192];
alignas(4096) static U pc_[MAXN * MAXN + (MAXN/BASE)*(MAXN/BASE)*TILE_PAD + 8192], ws_[MAXN * MAXN + (MAXN/BASE)*(MAXN/BASE)*TILE_PAD + 8192];
static U *pa = pa_ + 1024, *pb = pb_ + 1024, *pc = pc_ + 0, *ws = ws_ + 1024;
static U *final_out;
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 const U *raw_basis_input;
static inline __m256i raw_basis_low32(const U *p,bool transformed) {
__m256i v=_mm256_load_si256((const __m256i*)p);
if(transformed) v=_mm256_add_epi16(v,_mm256_sub_epi16(
_mm256_load_si256((const __m256i*)(p-32*2048)),
_mm256_load_si256((const __m256i*)(p-32))));
return v;
}
static inline __m128i raw_basis_low16(const U *p,bool transformed) {
__m128i v=_mm_load_si128((const __m128i*)p);
if(transformed) v=_mm_add_epi16(v,_mm_sub_epi16(
_mm_load_si128((const __m128i*)(p-32*2048)),
_mm_load_si128((const __m128i*)(p-32))));
return v;
}
static inline __m256i raw_basis32(const U *p,int flags) {
__m256i v=raw_basis_low32(p,flags&1);
if(flags&2) v=_mm256_add_epi16(v,_mm256_sub_epi16(
raw_basis_low32(p-64*2048,flags&1),raw_basis_low32(p-64,flags&1)));
return v;
}
static inline __m128i raw_basis16(const U *p,int flags) {
__m128i v=raw_basis_low16(p,flags&1);
if(flags&2) v=_mm_add_epi16(v,_mm_sub_epi16(
raw_basis_low16(p-64*2048,flags&1),raw_basis_low16(p-64,flags&1)));
return v;
}
static void pack(U *d, const U *s, int n, int stride, bool right) {
if (n > BASE) {
int h = n / 2, q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
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) {
const int transformed=(((s-raw_basis_input)&32) && ((s-raw_basis_input)&(32*2048)))?1:0;
for(int i=0;i<n;i+=4) for(int k=0;k<n;k+=16) {
{ int i2=(i+4)&31; /* s4ct-PKLEAF64: rows i+4..i+7, both k blocks */
__builtin_prefetch(s+i2*stride+k);__builtin_prefetch(s+(i2+1)*stride+k);
__builtin_prefetch(s+i2*stride+k+16);__builtin_prefetch(s+(i2+1)*stride+k+16);
__builtin_prefetch(s+(i2+2)*stride+k);__builtin_prefetch(s+(i2+3)*stride+k);
__builtin_prefetch(s+(i2+2)*stride+k+16);__builtin_prefetch(s+(i2+3)*stride+k+16); }
__m256i a=raw_basis32(s+i*stride+k,transformed), b=raw_basis32(s+(i+1)*stride+k,transformed);
__m256i c=raw_basis32(s+(i+2)*stride+k,transformed), e=raw_basis32(s+(i+3)*stride+k,transformed);
__m256i t0=_mm256_unpacklo_epi32(a,b),t1=_mm256_unpacklo_epi32(c,e);
__m256i t2=_mm256_unpackhi_epi32(a,b),t3=_mm256_unpackhi_epi32(c,e);
__m256i u0=_mm256_unpacklo_epi64(t0,t1),u1=_mm256_unpackhi_epi64(t0,t1);
__m256i u2=_mm256_unpacklo_epi64(t2,t3),u3=_mm256_unpackhi_epi64(t2,t3);
st(d,_mm256_permute2x128_si256(u0,u1,0x20));
st(d+16,_mm256_permute2x128_si256(u2,u3,0x20));
st(d+32,_mm256_permute2x128_si256(u0,u1,0x31));
st(d+48,_mm256_permute2x128_si256(u2,u3,0x31));d+=64;
}
} else {
const int transformed=(((s-raw_basis_input)&32) && ((s-raw_basis_input)&(32*2048)))?1:0;
for(int k=0;k<n;k+=2) for(int j=0;j<n;) {
int w=32;
U *v=d+j*n+k*w;
for(int t=0;t<w;t+=16) {
if(t+16<=w) {
__m256i x=raw_basis32(s+k*stride+j+t,transformed);
__m256i y=raw_basis32(s+(k+1)*stride+j+t,transformed);
st(v,_mm256_unpacklo_epi16(x,y));
st(v+16,_mm256_unpackhi_epi16(x,y));v+=32;
} else {
__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));
_mm_store_si128((__m128i*)v,_mm_unpacklo_epi16(x,y));
_mm_store_si128((__m128i*)(v+8),_mm_unpackhi_epi16(x,y));v+=16;
}
}
j+=w;
}
}
}
// ===== s4de_ : pack 左葉の広幅読み（読取順序のみ・単一変数） =====
// 原状: 1 つの 32x32 副象限ごとに 64 B/行 を読む（4 回の pass に分かれる）。
// 本件: 64 列ぶんの 2 tile（列 0..63 と 64..127）を 1 つの pass で読み、256 B/行 を連続で消費する。
// 読み書きするバイト集合・各 tile 内の書き込みレイアウト・ストア命令は原状と逐位同一。
// transformed は実行期判定（(s-raw) の bit5 と bit17）を、呼び出し幾何から
// 副象限 (行off==32 && 列off==32) に等価（証: 他の全オフセット 2048/2048*4096/64/64*4096/
// 128i/128j*4096 は bit5 にも bit17 にも寄与しない）と定数化した。ck で逐位一致を確認済み。
#define PWL_Q ((32*32) + (32/BASE)*(32/BASE)*TILE_PAD)
#define PWL_CH(DSTP,K,TF) do { \
__m256i a=raw_basis32(sr+(long)i*S4ST+(K),(TF)), b=raw_basis32(sr+(long)(i+1)*S4ST+(K),(TF)); \
__m256i c=raw_basis32(sr+(long)(i+2)*S4ST+(K),(TF)), e=raw_basis32(sr+(long)(i+3)*S4ST+(K),(TF)); \
__m256i t0=_mm256_unpacklo_epi32(a,b),t1=_mm256_unpacklo_epi32(c,e); \
__m256i t2=_mm256_unpackhi_epi32(a,b),t3=_mm256_unpackhi_epi32(c,e); \
__m256i u0=_mm256_unpacklo_epi64(t0,t1),u1=_mm256_unpackhi_epi64(t0,t1); \
__m256i u2=_mm256_unpacklo_epi64(t2,t3),u3=_mm256_unpackhi_epi64(t2,t3); \
st((DSTP)+0,_mm256_permute2x128_si256(u0,u1,0x20)); \
st((DSTP)+16,_mm256_permute2x128_si256(u2,u3,0x20)); \
st((DSTP)+32,_mm256_permute2x128_si256(u0,u1,0x31)); \
st((DSTP)+48,_mm256_permute2x128_si256(u2,u3,0x31)); } while(0)
static void pack_wide(U *dA,U *dB,const U *s) {
const long S4ST=2048;
for (int r=0;r<2;++r) {
const U *sr = s + (long)r*32*S4ST;
U *p0=dA+(r?2*PWL_Q:0), *p1=dA+(r?3*PWL_Q:PWL_Q), *p2=dB+(r?2*PWL_Q:0), *p3=dB+(r?3*PWL_Q:PWL_Q);
const int TFc = (r==1)?1:0;
for (int i=0;i<32;i+=4) {
{ int i2=(i+4)&31;
__builtin_prefetch(sr+(long)i2*S4ST); __builtin_prefetch(sr+(long)(i2+1)*S4ST);
__builtin_prefetch(sr+(long)(i2+2)*S4ST); __builtin_prefetch(sr+(long)(i2+3)*S4ST);
__builtin_prefetch(sr+(long)i2*S4ST+16); __builtin_prefetch(sr+(long)(i2+1)*S4ST+16);
__builtin_prefetch(sr+(long)(i2+2)*S4ST+16); __builtin_prefetch(sr+(long)(i2+3)*S4ST+16); }
PWL_CH(p0, 0,0); PWL_CH(p0+64, 16,0); p0+=128;
PWL_CH(p1, 32,TFc); PWL_CH(p1+64, 48,TFc); p1+=128;
PWL_CH(p2, 64,0); PWL_CH(p2+64, 80,0); p2+=128;
PWL_CH(p3, 96,TFc); PWL_CH(p3+64,112,TFc); p3+=128;
}
}
}
// ===== s4df_ : pack 右葉（B オペランド）の広幅読み（s4de_ の左葉と同型・転置対読み版） =====
// 原状の右葉は 1 つの 64x64 tile を 4 つの 32x32 副象限（(r,c) = 行帯 x 列帯）に分け、
// 各副象限を「k 行対（k=0..30 step 2）x 16 字チャンク 2 つ」で読む ⇒ 副象限ごとに別 pass。
// 横に隣接する 2 tile（列 0..63 と 64..127）も別 pass ⇒ 同じ行を 2 回（128 B + 128 B）に分けて読む。
// 本件 `pack_wide_b` は **2 tile x 4 副象限 = 8 本の読みを 1 pass に統合**し、1 行あたり
// **256 B を連続消費**する（機構・目標ポインタ構造は左葉 `pack_wide` と同一）。
// 読むバイト集合・各副象限内の書き込みレイアウト・ストア命令・アドレスは原状と逐位同一。
// `transformed` は原状どおり**副象限ごとの実行期判定**（= 基点 X の bit5/bit17 が r,c で反転）
// として再現する（左葉のように定数化はしない：単一変数を読取順序だけに保つため）。
// 予取（prefetch）は**足さない**（(1528-律二)：右葉は同因不同薬）。
#define S4DF_BW(P,J,TF) do { \
__m256i x=raw_basis32(sr+(long)k*ST+(J),(TF)), y=raw_basis32(sr+(long)(k+1)*ST+(J),(TF)); \
st((P),_mm256_unpacklo_epi16(x,y)); st((P)+16,_mm256_unpackhi_epi16(x,y)); } while(0)
static void pack_wide_b(U *dA,U *dB,const U *s) {
const long ST=2048;
const long d0=s-raw_basis_input;
const int b5=(int)((d0>>5)&1), b17=(int)((d0>>16)&1);
const int t00=b5&b17, t01=(b5^1)&b17, t10=b5&(b17^1), t11=(b5^1)&(b17^1);
for (int r=0;r<2;++r) {
const U *sr=s+(long)r*32*ST;
U *p0=dA+(r?2*PWL_Q:0), *p1=dA+(r?3*PWL_Q:PWL_Q), *p2=dB+(r?2*PWL_Q:0), *p3=dB+(r?3*PWL_Q:PWL_Q);
const int tf0=r?t10:t00; /* 列帯 0 の旗標 */
const int tf1=r?t11:t01; /* 列帯 1 の旗標 */
for (int k=0;k<32;k+=2) {
S4DF_BW(p0+0 , 0,tf0); S4DF_BW(p0+32, 16,tf0);
S4DF_BW(p1+0 , 32,tf1); S4DF_BW(p1+32, 48,tf1);
S4DF_BW(p2+0 , 64,tf0); S4DF_BW(p2+32, 80,tf0);
S4DF_BW(p3+0 , 96,tf1); S4DF_BW(p3+32,112,tf1);
p0+=64; p1+=64; p2+=64; p3+=64;
}
}
}
static void unpack(U *d, const U *s, int n, int stride) {
if (n > BASE) {
int h = n/2, q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
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_load_si256((const __m256i*)(sp+k)));
}
}
}
// ===== s4bj_ : 叶尾改 blend 后，每个 16 字组内被 sigma 置换，decode 前先复原 =====
// sigma = [0,4,1,5,2,6,3,7 | 8,12,9,13,10,14,11,15] (new[j] = old[sigma(j)])
// sigma^-1= [0,2,4,6,1,3,5,7 | 8,10,12,14,9,11,13,15] (old[i] = new[sigma^-1(i)])
alignas(32) static const unsigned char s4bj_si[32] = {
0,1,4,5,8,9,12,13,2,3,6,7,10,11,14,15, 0,1,4,5,8,9,12,13,2,3,6,7,10,11,14,15};
static inline __m256i s4bj_fx(__m256i v) {
return _mm256_shuffle_epi8(v, _mm256_load_si256((const __m256i *)s4bj_si));
}
alignas(32) unsigned short s4aj_m16[16] = {0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0,0xFFFF,0};
alignas(32) unsigned char s4ak_sh16[32] = {0,1,4,5,8,9,12,13, 0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,1,4,5,8,9,12,13, 0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80};
// One explicit panel/row loop keeps exactly one copy of the large kernel.
template<int Mode>
static inline void native2x32_cache3(const U *a,const U *b,U *c) {
const U *a_final=a+1024;
intptr_t astep=8; // byte steps8,248 alternate within and across4-row A packs
asm volatile(
".p2align 6\n\t"
"3:\n\t"
"vmovdqa 0(%[b]), %%ymm8\n\t"
"vmovdqa 32(%[b]), %%ymm9\n\t"
"vmovdqa 64(%[b]), %%ymm10\n\t"
"vpbroadcastd 0(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm0\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm1\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm2\n\t"
"vpmaddwd 96(%[b]), %%ymm11, %%ymm3\n\t"
"vpbroadcastd 4(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm4\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm5\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm6\n\t"
"vpmaddwd 96(%[b]), %%ymm11, %%ymm7\n\t"
"vmovdqa 128(%[b]), %%ymm8\n\t"
"vmovdqa 160(%[b]), %%ymm9\n\t"
"vmovdqa 192(%[b]), %%ymm10\n\t"
"vpbroadcastd 16(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 224(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 20(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 224(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 256(%[b]), %%ymm8\n\t"
"vmovdqa 288(%[b]), %%ymm9\n\t"
"vmovdqa 320(%[b]), %%ymm10\n\t"
"vpbroadcastd 32(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 352(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 36(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 352(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 384(%[b]), %%ymm8\n\t"
"vmovdqa 416(%[b]), %%ymm9\n\t"
"vmovdqa 448(%[b]), %%ymm10\n\t"
"vpbroadcastd 48(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 480(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 52(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 480(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 512(%[b]), %%ymm8\n\t"
"vmovdqa 544(%[b]), %%ymm9\n\t"
"vmovdqa 576(%[b]), %%ymm10\n\t"
"vpbroadcastd 64(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 608(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 68(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 608(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 640(%[b]), %%ymm8\n\t"
"vmovdqa 672(%[b]), %%ymm9\n\t"
"vmovdqa 704(%[b]), %%ymm10\n\t"
"vpbroadcastd 80(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 736(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 84(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 736(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 768(%[b]), %%ymm8\n\t"
"vmovdqa 800(%[b]), %%ymm9\n\t"
"vmovdqa 832(%[b]), %%ymm10\n\t"
"vpbroadcastd 96(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 864(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 100(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 864(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 896(%[b]), %%ymm8\n\t"
"vmovdqa 928(%[b]), %%ymm9\n\t"
"vmovdqa 960(%[b]), %%ymm10\n\t"
"vpbroadcastd 112(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 992(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 116(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 992(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1024(%[b]), %%ymm8\n\t"
"vmovdqa 1056(%[b]), %%ymm9\n\t"
"vmovdqa 1088(%[b]), %%ymm10\n\t"
"vpbroadcastd 128(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 1120(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 132(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 1120(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1152(%[b]), %%ymm8\n\t"
"vmovdqa 1184(%[b]), %%ymm9\n\t"
"vmovdqa 1216(%[b]), %%ymm10\n\t"
"vpbroadcastd 144(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 1248(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 148(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 1248(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1280(%[b]), %%ymm8\n\t"
"vmovdqa 1312(%[b]), %%ymm9\n\t"
"vmovdqa 1344(%[b]), %%ymm10\n\t"
"vpbroadcastd 160(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 1376(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 164(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 1376(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1408(%[b]), %%ymm8\n\t"
"vmovdqa 1440(%[b]), %%ymm9\n\t"
"vmovdqa 1472(%[b]), %%ymm10\n\t"
"vpbroadcastd 176(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 1504(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 180(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 1504(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1536(%[b]), %%ymm8\n\t"
"vmovdqa 1568(%[b]), %%ymm9\n\t"
"vmovdqa 1600(%[b]), %%ymm10\n\t"
"vpbroadcastd 192(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 1632(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 196(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 1632(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1664(%[b]), %%ymm8\n\t"
"vmovdqa 1696(%[b]), %%ymm9\n\t"
"vmovdqa 1728(%[b]), %%ymm10\n\t"
"vpbroadcastd 208(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 1760(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 212(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 1760(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1792(%[b]), %%ymm8\n\t"
"vmovdqa 1824(%[b]), %%ymm9\n\t"
"vmovdqa 1856(%[b]), %%ymm10\n\t"
"vpbroadcastd 224(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 1888(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 228(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 1888(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vmovdqa 1920(%[b]), %%ymm8\n\t"
"vmovdqa 1952(%[b]), %%ymm9\n\t"
"vmovdqa 1984(%[b]), %%ymm10\n\t"
"vpbroadcastd 240(%[a]), %%ymm11\n\t"
"vpmaddwd %%ymm8, %%ymm11, %%ymm12\n\t"
"vpmaddwd %%ymm9, %%ymm11, %%ymm13\n\t"
"vpmaddwd %%ymm10, %%ymm11, %%ymm14\n\t"
"vpmaddwd 2016(%[b]), %%ymm11, %%ymm11\n\t"
"vpbroadcastd 244(%[a]), %%ymm15\n\t"
"vpmaddwd %%ymm8, %%ymm15, %%ymm8\n\t"
"vpmaddwd %%ymm9, %%ymm15, %%ymm9\n\t"
"vpmaddwd %%ymm10, %%ymm15, %%ymm10\n\t"
"vpmaddwd 2016(%[b]), %%ymm15, %%ymm15\n\t"
"vpaddd %%ymm12, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm13, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm11, %%ymm3, %%ymm3\n\t"
"vpaddd %%ymm8, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm9, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpslld $16, %%ymm1, %%ymm1\n\t"
"vpblendw $0xAA, %%ymm1, %%ymm0, %%ymm0\n\t"
".if %c[mode] == 1\n\t"
"vpaddw 0(%[c]), %%ymm0, %%ymm0\n\t"
".endif\n\t"
"vmovdqu %%ymm0, 0(%[c])\n\t"
"vpslld $16, %%ymm3, %%ymm3\n\t"
"vpblendw $0xAA, %%ymm3, %%ymm2, %%ymm2\n\t"
".if %c[mode] == 1\n\t"
"vpaddw 32(%[c]), %%ymm2, %%ymm2\n\t"
".endif\n\t"
"vmovdqu %%ymm2, 32(%[c])\n\t"
"vpslld $16, %%ymm5, %%ymm5\n\t"
"vpblendw $0xAA, %%ymm5, %%ymm4, %%ymm4\n\t"
".if %c[mode] == 1\n\t"
"vpaddw 64(%[c]), %%ymm4, %%ymm4\n\t"
".endif\n\t"
"vmovdqu %%ymm4, 64(%[c])\n\t"
"vpslld $16, %%ymm7, %%ymm7\n\t"
"vpblendw $0xAA, %%ymm7, %%ymm6, %%ymm6\n\t"
".if %c[mode] == 1\n\t"
"vpaddw 96(%[c]), %%ymm6, %%ymm6\n\t"
".endif\n\t"
"vmovdqu %%ymm6, 96(%[c])\n\t"
"add %[astep], %[a]\n\t"
"xor $240, %[astep]\n\t"
"add $128, %[c]\n\t"
"cmp %[a_final], %[a]\n\t"
"jb 3b\n\t"
: [a] "+&r"(a), [c] "+&r"(c), [astep] "+&r"(astep)
: [b] "r"(b), [a_final] "r"(a_final), [mode] "i"(Mode)
: "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) {
native2x32_cache3<0>(a,b,c);
}
// Only one additional hot body: M1, M3 and the negated-M7 update all add.
static __attribute__((noinline)) void leaf_add(const U*a,const U*b,U*c){
native2x32_cache3<1>(a,b,c);
}
static inline void emit_native(U *dst,const U *s,int stride) {
for(int row=0;row<4;++row) {
U *d=dst+(size_t)row*stride;
_mm256_stream_si256((__m256i*)(d+0),s4bj_fx(ld(s+32*row)));
_mm256_stream_si256((__m256i*)(d+16),s4bj_fx(ld(s+32*row+16)));
}
}
template<bool Sub>
static inline void combine(U *d,const U *a,const U *b,int len) {
#pragma GCC unroll 1
for(int chunk=0;chunk<len;chunk+=1024+TILE_PAD) {
const U *end=a+1024;
if constexpr(Sub) {
asm volatile(
"1:\n\t"
"vmovdqa 0(%[a]), %%ymm0\n\t"
"vmovdqa 32(%[a]), %%ymm1\n\t"
"vmovdqa 64(%[a]), %%ymm2\n\t"
"vmovdqa 96(%[a]), %%ymm3\n\t"
"vpsubw 0(%[b]), %%ymm0, %%ymm0\n\t"
"vpsubw 32(%[b]), %%ymm1, %%ymm1\n\t"
"vpsubw 64(%[b]), %%ymm2, %%ymm2\n\t"
"vpsubw 96(%[b]), %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm0, 0(%[d])\n\t"
"vmovdqa %%ymm1, 32(%[d])\n\t"
"vmovdqa %%ymm2, 64(%[d])\n\t"
"vmovdqa %%ymm3, 96(%[d])\n\t"
"add $128, %[a]\n\t"
"add $128, %[b]\n\t"
"add $128, %[d]\n\t"
"cmp %[end], %[a]\n\t"
"jb 1b\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [d] "+&r"(d)
: [end] "r"(end)
: "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3");
} else {
asm volatile(
"1:\n\t"
"vmovdqa 0(%[a]), %%ymm0\n\t"
"vmovdqa 32(%[a]), %%ymm1\n\t"
"vmovdqa 64(%[a]), %%ymm2\n\t"
"vmovdqa 96(%[a]), %%ymm3\n\t"
"vpaddw 0(%[b]), %%ymm0, %%ymm0\n\t"
"vpaddw 32(%[b]), %%ymm1, %%ymm1\n\t"
"vpaddw 64(%[b]), %%ymm2, %%ymm2\n\t"
"vpaddw 96(%[b]), %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm0, 0(%[d])\n\t"
"vmovdqa %%ymm1, 32(%[d])\n\t"
"vmovdqa %%ymm2, 64(%[d])\n\t"
"vmovdqa %%ymm3, 96(%[d])\n\t"
"add $128, %[a]\n\t"
"add $128, %[b]\n\t"
"add $128, %[d]\n\t"
"cmp %[end], %[a]\n\t"
"jb 1b\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [d] "+&r"(d)
: [end] "r"(end)
: "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3");
}
a+=TILE_PAD;b+=TILE_PAD;d+=TILE_PAD;
}
}
template<bool Sub,int Delta>
static inline void combine_displaced(U *d,const U *a,int len) {
#pragma GCC unroll 1
for(int chunk=0;chunk<len;chunk+=1024+TILE_PAD) {
const U *end=a+1024;
if constexpr(Sub) {
asm volatile("1:\n\t"
"vmovdqa 0(%[a]), %%ymm0\n\t"
"vmovdqa 32(%[a]), %%ymm1\n\t"
"vmovdqa 64(%[a]), %%ymm2\n\t"
"vmovdqa 96(%[a]), %%ymm3\n\t"
"vpsubw %c[delta]+0(%[a]), %%ymm0, %%ymm0\n\t"
"vpsubw %c[delta]+32(%[a]), %%ymm1, %%ymm1\n\t"
"vpsubw %c[delta]+64(%[a]), %%ymm2, %%ymm2\n\t"
"vpsubw %c[delta]+96(%[a]), %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm0, 0(%[d])\n\t"
"vmovdqa %%ymm1, 32(%[d])\n\t"
"vmovdqa %%ymm2, 64(%[d])\n\t"
"vmovdqa %%ymm3, 96(%[d])\n\t"
"add $128, %[a]\n\t"
"add $128, %[d]\n\t"
"cmp %[end], %[a]\n\t"
"jb 1b\n\t"
: [a] "+&r"(a), [d] "+&r"(d)
: [end] "r"(end), [delta] "i"(2*Delta)
: "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3");
} else {
asm volatile("1:\n\t"
"vmovdqa 0(%[a]), %%ymm0\n\t"
"vmovdqa 32(%[a]), %%ymm1\n\t"
"vmovdqa 64(%[a]), %%ymm2\n\t"
"vmovdqa 96(%[a]), %%ymm3\n\t"
"vpaddw %c[delta]+0(%[a]), %%ymm0, %%ymm0\n\t"
"vpaddw %c[delta]+32(%[a]), %%ymm1, %%ymm1\n\t"
"vpaddw %c[delta]+64(%[a]), %%ymm2, %%ymm2\n\t"
"vpaddw %c[delta]+96(%[a]), %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm0, 0(%[d])\n\t"
"vmovdqa %%ymm1, 32(%[d])\n\t"
"vmovdqa %%ymm2, 64(%[d])\n\t"
"vmovdqa %%ymm3, 96(%[d])\n\t"
"add $128, %[a]\n\t"
"add $128, %[d]\n\t"
"cmp %[end], %[a]\n\t"
"jb 1b\n\t"
: [a] "+&r"(a), [d] "+&r"(d)
: [end] "r"(end), [delta] "i"(2*Delta)
: "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3");
}
a+=TILE_PAD;d+=TILE_PAD;
}
}
// pdoom frontier6: the output is its own first operand, so only two
// advancing addresses are necessary. N64 calls touch a single valid leaf.
template<bool Sub>
static inline void combine_inplace(U *d,const U *b,int len) {
#pragma GCC unroll 1
for(int chunk=0;chunk<len;chunk+=1024+TILE_PAD) {
const U *end=d+1024;
if constexpr(Sub) {
asm volatile("1:\n\t"
"vmovdqa 0(%[d]), %%ymm0\n\t"
"vmovdqa 32(%[d]), %%ymm1\n\t"
"vmovdqa 64(%[d]), %%ymm2\n\t"
"vmovdqa 96(%[d]), %%ymm3\n\t"
"vpsubw 0(%[b]), %%ymm0, %%ymm0\n\t"
"vpsubw 32(%[b]), %%ymm1, %%ymm1\n\t"
"vpsubw 64(%[b]), %%ymm2, %%ymm2\n\t"
"vpsubw 96(%[b]), %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm0, 0(%[d])\n\t"
"vmovdqa %%ymm1, 32(%[d])\n\t"
"vmovdqa %%ymm2, 64(%[d])\n\t"
"vmovdqa %%ymm3, 96(%[d])\n\t"
"add $128, %[d]\n\t"
"add $128, %[b]\n\t"
"cmp %[end], %[d]\n\t"
"jb 1b\n\t"
: [d] "+&r"(d), [b] "+&r"(b)
: [end] "r"(end)
: "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3");
} else {
asm volatile("1:\n\t"
"vmovdqa 0(%[d]), %%ymm0\n\t"
"vmovdqa 32(%[d]), %%ymm1\n\t"
"vmovdqa 64(%[d]), %%ymm2\n\t"
"vmovdqa 96(%[d]), %%ymm3\n\t"
"vpaddw 0(%[b]), %%ymm0, %%ymm0\n\t"
"vpaddw 32(%[b]), %%ymm1, %%ymm1\n\t"
"vpaddw 64(%[b]), %%ymm2, %%ymm2\n\t"
"vpaddw 96(%[b]), %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm0, 0(%[d])\n\t"
"vmovdqa %%ymm1, 32(%[d])\n\t"
"vmovdqa %%ymm2, 64(%[d])\n\t"
"vmovdqa %%ymm3, 96(%[d])\n\t"
"add $128, %[d]\n\t"
"add $128, %[b]\n\t"
"cmp %[end], %[d]\n\t"
"jb 1b\n\t"
: [d] "+&r"(d), [b] "+&r"(b)
: [end] "r"(end)
: "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3");
}
d+=TILE_PAD;b+=TILE_PAD;
}
}
// pdoom frontier6: the output is its own first operand, so only two
// advancing addresses are necessary. N64 calls touch a single valid leaf.
static inline void combine_reverse_inplace(U *d,const U *b,int len) {
#pragma GCC unroll 1
for(int chunk=0;chunk<len;chunk+=1024+TILE_PAD) {
const U *end=d+1024;
asm volatile("1:\n\t"
"vmovdqa 0(%[b]), %%ymm0\n\t"
"vmovdqa 32(%[b]), %%ymm1\n\t"
"vmovdqa 64(%[b]), %%ymm2\n\t"
"vmovdqa 96(%[b]), %%ymm3\n\t"
"vpsubw 0(%[d]), %%ymm0, %%ymm0\n\t"
"vpsubw 32(%[d]), %%ymm1, %%ymm1\n\t"
"vpsubw 64(%[d]), %%ymm2, %%ymm2\n\t"
"vpsubw 96(%[d]), %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm0, 0(%[d])\n\t"
"vmovdqa %%ymm1, 32(%[d])\n\t"
"vmovdqa %%ymm2, 64(%[d])\n\t"
"vmovdqa %%ymm3, 96(%[d])\n\t"
"add $128, %[d]\n\t"
"add $128, %[b]\n\t"
"cmp %[end], %[d]\n\t"
"jb 1b\n\t"
: [d] "+&r"(d), [b] "+&r"(b)
: [end] "r"(end)
: "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3");
d+=TILE_PAD;b+=TILE_PAD;
}
}
// At the root, the three completed packed quadrants can be sent straight
// to their final row-major destinations rather than materialized in pc.
template<int N>
static __attribute__((noinline)) void combine3_unpack(
U *d12,U *d21,U *d22,
const U *p6,const U *p7,const U *p5,
const U *p1,const U *p4,const U *p3,int stride) {
if constexpr (N > BASE) {
constexpr int h=N/2,q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
combine3_unpack<h>(d12,d21,d22,p6,p7,p5,p1,p4,p3,stride);
combine3_unpack<h>(d12+h,d21+h,d22+h,p6+q,p7+q,p5+q,p1+q,p4+q,p3+q,stride);
combine3_unpack<h>(d12+h*stride,d21+h*stride,d22+h*stride,
p6+2*q,p7+2*q,p5+2*q,p1+2*q,p4+2*q,p3+2*q,stride);
combine3_unpack<h>(d12+h*stride+h,d21+h*stride+h,d22+h*stride+h,
p6+3*q,p7+3*q,p5+3*q,p1+3*q,p4+3*q,p3+3*q,stride);
} else {
alignas(32) U temp[3*128];
for(int block=0;block<1024;block=block+128) {
for(int v=0;v<128;v+=16) {
int idx=block+v;
__m256i u2=_mm256_add_epi16(ld(p1+idx),ld(p6+idx));
__m256i u3=_mm256_add_epi16(u2,ld(p7+idx));
__m256i v5=ld(p5+idx);
st(temp+v,_mm256_add_epi16(_mm256_add_epi16(u2,v5),ld(p3+idx)));
st(temp+128+v,_mm256_sub_epi16(u3,ld(p4+idx)));
st(temp+256+v,_mm256_add_epi16(u3,v5));
}
size_t off=(size_t)(block/32)*stride;
emit_native(d12+off,temp,stride);
emit_native(d21+off,temp+128,stride);
emit_native(d22+off,temp+256,stride);
}
}
}
// Add the surviving P1 into the root C11 product while mapping its packed
// quadrant layout to the final row-major C11 cells.
template<int N>
static __attribute__((noinline)) void unpack_add(U *d,const U *p2,const U *p1,int stride) {
if constexpr(N>BASE) {
constexpr int h=N/2,q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
unpack_add<h>(d,p2,p1,stride);
unpack_add<h>(d+h,p2+q,p1+q,stride);
unpack_add<h>(d+h*stride,p2+2*q,p1+2*q,stride);
unpack_add<h>(d+h*stride+h,p2+3*q,p1+3*q,stride);
} else {
alignas(32) U temp[128];
for(int block=0;block<1024;block+=128) {
for(int v=0;v<128;v+=16)
st(temp+v,_mm256_add_epi16(ld(p2+block+v),ld(p1+block+v)));
emit_native(d+(size_t)(block/32)*stride,temp,stride);
}
}
}
// Strassen-Winograd: seven products, fifteen additions, two temporaries.
// All operations are in Z/(2^16); truncation never changes the answer.
template<int N>
static __attribute__((noinline)) void root_fuse(
U *d11,U *d12,U *d21,U *d22,
const U *p1,const U *p2,const U *p3,const U *p4,const U *p5,const U *p6,const U *p7,int stride) {
if constexpr(N==128) {
constexpr int q=BASE*BASE+TILE_PAD;
root_fuse<64>(d11,d12,d21,d22,p1,p2,p3,p4,p5,p6,p7,stride);
alignas(32) U keep[4*512], cur[4*512];
U *dest[4]={d11,d12,d21,d22};
for(int block=0;block<1024;block+=128) {
auto make_inner=[&](U *tmp,int group) {
for(int child=0;child<4;child++) {
for(int v=0;v<128;v+=16) {
int idx=(4*group+child)*q+block+v;
if(!(v&16)){__builtin_prefetch(p1+idx+128,0,3);__builtin_prefetch(p2+idx+128,0,3);__builtin_prefetch(p3+idx+128,0,3);__builtin_prefetch(p4+idx+128,0,3);__builtin_prefetch(p5+idx+128,0,3);__builtin_prefetch(p6+idx+128,0,3);__builtin_prefetch(p7+idx+128,0,3);}
__m256i v1=ld(p1+idx),v5=ld(p5+idx);
__m256i u2=_mm256_add_epi16(v1,ld(p6+idx));
__m256i u3=_mm256_add_epi16(u2,ld(p7+idx));
U *t=tmp+child*512;
st(t+v,_mm256_add_epi16(v1,ld(p2+idx)));
st(t+128+v,_mm256_add_epi16(_mm256_add_epi16(u2,v5),ld(p3+idx)));
st(t+256+v,_mm256_sub_epi16(u3,ld(p4+idx)));
st(t+384+v,_mm256_add_epi16(u3,v5));
}
}
for(int v=0;v<512;v+=16) {
__m256i c22=ld(tmp+3*512+v);
st(tmp+1*512+v,_mm256_sub_epi16(ld(tmp+1*512+v),c22));
st(tmp+2*512+v,_mm256_sub_epi16(c22,ld(tmp+2*512+v)));
}
};
size_t off=(size_t)(block/32)*stride;
make_inner(keep,3);
for(int child=0;child<4;child++) {
int row=64+((child>>1)&1)*32,col=64+(child&1)*32;
for(int r=0;r<4;r++)emit_native(dest[r]+off+row*stride+col,keep+child*512+r*128,stride);
}
for(int group=1;group<=2;group++) {
make_inner(cur,group);
for(int child=0;child<4;child++) {
int row=((group>>1)&1)*64+((child>>1)&1)*32;
int col=(group&1)*64+(child&1)*32;
for(int r=0;r<4;r++)for(int rr=0;rr<4;rr++) {
U *d=dest[r]+off+(row+rr)*stride+col;
const U *x=cur+child*512+r*128+rr*32;
const U *y=keep+child*512+r*128+rr*32;
__m256i v0,v1;
if(group==1) {
v0=_mm256_sub_epi16(ld(x),ld(y));
v1=_mm256_sub_epi16(ld(x+16),ld(y+16));
} else {
v0=_mm256_sub_epi16(ld(y),ld(x));
v1=_mm256_sub_epi16(ld(y+16),ld(x+16));
}
_mm256_stream_si256((__m256i*)d,s4bj_fx(v0));
_mm256_stream_si256((__m256i*)(d+16),s4bj_fx(v1));
}
}
}
}
} else if constexpr(N==64) {
constexpr int q=BASE*BASE+TILE_PAD;
root_fuse<32>(d11,d12,d21,d22,p1,p2,p3,p4,p5,p6,p7,stride);
alignas(32) U tmp[3*4*128];
U *dest[4]={d11,d12,d21,d22};
for(int block=0;block<1024;block+=128) {
for(int child=1;child<=3;child++) {
for(int v=0;v<128;v+=16) {
int idx=child*q+block+v;
if(!(v&16)){__builtin_prefetch(p1+idx+128,0,3);__builtin_prefetch(p2+idx+128,0,3);__builtin_prefetch(p3+idx+128,0,3);__builtin_prefetch(p4+idx+128,0,3);__builtin_prefetch(p5+idx+128,0,3);__builtin_prefetch(p6+idx+128,0,3);__builtin_prefetch(p7+idx+128,0,3);}
__m256i v1=ld(p1+idx),v5=ld(p5+idx);
__m256i u2=_mm256_add_epi16(v1,ld(p6+idx));
__m256i u3=_mm256_add_epi16(u2,ld(p7+idx));
U *t=tmp+(child-1)*512;
st(t+v,_mm256_add_epi16(v1,ld(p2+idx)));
st(t+128+v,_mm256_add_epi16(_mm256_add_epi16(u2,v5),ld(p3+idx)));
st(t+256+v,_mm256_sub_epi16(u3,ld(p4+idx)));
st(t+384+v,_mm256_add_epi16(u3,v5));
}
}
for(int v=0;v<512;v+=16) {
__m256i c22=ld(tmp+1024+v);
st(tmp+v,_mm256_sub_epi16(ld(tmp+v),c22));
st(tmp+512+v,_mm256_sub_epi16(c22,ld(tmp+512+v)));
}
size_t off=(size_t)(block/32)*stride;
for(int r=0;r<4;r++) {
emit_native(dest[r]+off+32,tmp+r*128,stride);
emit_native(dest[r]+off+32*stride,tmp+512+r*128,stride);
emit_native(dest[r]+off+32*stride+32,tmp+1024+r*128,stride);
}
}
} else if constexpr (N > BASE) {
constexpr int h=N/2,q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
root_fuse<h>(d11,d12,d21,d22,p1,p2,p3,p4,p5,p6,p7,stride);
root_fuse<h>(d11+h,d12+h,d21+h,d22+h,p1+q,p2+q,p3+q,p4+q,p5+q,p6+q,p7+q,stride);
root_fuse<h>(d11+h*stride,d12+h*stride,d21+h*stride,d22+h*stride,
p1+2*q,p2+2*q,p3+2*q,p4+2*q,p5+2*q,p6+2*q,p7+2*q,stride);
root_fuse<h>(d11+h*stride+h,d12+h*stride+h,d21+h*stride+h,d22+h*stride+h,
p1+3*q,p2+3*q,p3+3*q,p4+3*q,p5+3*q,p6+3*q,p7+3*q,stride);
} else {
alignas(32) U temp[4*128];
for(int block=0;block<1024;block+=128) {
for(int v=0;v<128;v+=16) {
int idx=block+v;
if(!(v&16)){__builtin_prefetch(p1+idx+128,0,3);__builtin_prefetch(p2+idx+128,0,3);__builtin_prefetch(p3+idx+128,0,3);__builtin_prefetch(p4+idx+128,0,3);__builtin_prefetch(p5+idx+128,0,3);__builtin_prefetch(p6+idx+128,0,3);__builtin_prefetch(p7+idx+128,0,3);}
__m256i u2=_mm256_add_epi16(ld(p1+idx),ld(p6+idx));
__m256i u3=_mm256_add_epi16(u2,ld(p7+idx));
__m256i v5=ld(p5+idx);
st(temp+v,_mm256_add_epi16(_mm256_add_epi16(u2,v5),ld(p3+idx)));
st(temp+128+v,_mm256_sub_epi16(u3,ld(p4+idx)));
st(temp+256+v,_mm256_add_epi16(u3,v5));
st(temp+384+v,_mm256_add_epi16(ld(p2+idx),ld(p1+idx)));
}
size_t off=(size_t)(block/32)*stride;
emit_native(d12+off,temp,stride);
emit_native(d21+off,temp+128,stride);
emit_native(d22+off,temp+256,stride);
emit_native(d11+off,temp+384,stride);
}
}
}
template<int N>
static __attribute__((noinline)) void mul(const U *a,const U *b,U *c,U *work) {
if constexpr (N<=BASE) {
base(a,b,c,N);
} else if constexpr (N<=ALTMAX) {
constexpr int h=N/2,q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
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+0,*y=work+2*q+704,*next=work+4*q+640;
mul<h>(a12,b21,c11,next); // m2
mul<h>(a22,b22,c22,next); // m4
combine_displaced<false,q>(x,a21,q);
combine_displaced<false,q>(y,b21,q);
mul<h>(x,y,c12,next); // m5
combine_displaced<true,-2*q>(x,a22,q);
combine_displaced<true,-2*q>(y,b22,q);
mul<h>(x,y,c21,next); // m6
for(int chunk=0;chunk<q;chunk+=1024+TILE_PAD)
_Pragma("GCC unroll 4")
for(int i=chunk;i<chunk+1024;i+=16)
st(c22+i,_mm256_sub_epi16(_mm256_add_epi16(ld(c12+i),ld(c21+i)),_mm256_add_epi16(ld(c11+i),ld(c22+i))));
if constexpr(N==64) {
leaf_add(a11,b11,c11); // m1 directly accumulates
combine_displaced<true,-3*q>(y,b22,q);
leaf_add(a21,y,c21); // m3 directly accumulates
combine_displaced<true,3*q>(x,a11,q); // negate the M7 A input
leaf_add(x,b12,c12); // c12+=(-oldM7)
} else {
mul<h>(a11,b11,x,next); // m1
combine_inplace<false>(c11,x,q);
combine_displaced<true,-3*q>(y,b22,q);
mul<h>(a21,y,x,next); // m3
combine_inplace<false>(c21,x,q);
combine_displaced<true,-3*q>(x,a22,q);
mul<h>(x,b12,y,next); // m7
combine_inplace<true>(c12,y,q);
}
} else {
constexpr int h=N/2,q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
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+192,*y=work+q+512,*next=work+q+q+512;
combine_displaced<true,-2*q>(y,b22,q);
combine_displaced<true,2*q>(x,a11,q);
mul<h>(x,y,c21,next); // P7
combine_displaced<false,q>(x,a21,q);
combine_displaced<true,-q>(y,b12,q);
mul<h>(x,y,c22,next); // P5
combine_inplace<true>(x,a11,q);
combine_reverse_inplace(y,b22,q);
mul<h>(x,y,c12,next); // P6
combine_reverse_inplace(x,a12,q);
mul<h>(x,b22,c11,next); // P3
combine_inplace<true>(y,b21,q);
mul<h>(a22,y,x,next); // P4
mul<h>(a11,b11,y,next); // P1
if constexpr (N == 2048) {
combine3_unpack<1024>(final_out+1024,
final_out+(size_t)1024*2048,
final_out+(size_t)1024*2048+1024,
c12,c21,c22,y,x,c11,2048);
} else {
for(int chunk=0;chunk<q;chunk+=1024+TILE_PAD) for(int i=chunk;i<chunk+1024;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));
}
}
mul<h>(a12,b21,c11,next); // P2
if constexpr (N == 2048) unpack_add<1024>(final_out,c11,y,2048);
else combine_inplace<false>(c11,y,q);
}
}
static constexpr int RQ=1024*1024+(1024/BASE)*(1024/BASE)*TILE_PAD,RS=RQ+ROOT_PAD;
alignas(4096) static U root_mem[15*RS+2048];
alignas(4096) static U root_tile[4*(64*64+(64/BASE)*(64/BASE)*TILE_PAD)];
static constexpr int HQ=64*64+(64/BASE)*(64/BASE)*TILE_PAD;
alignas(4096) static U high_left[4*HQ],high_right[4*HQ];
template<bool Right>
static void pack_four(U *tile,const U *s) {
pack(tile,s,64,2048,Right);
pack(tile+HQ,s+1024,64,2048,Right);
pack(tile+2*HQ,s+1024*2048,64,2048,Right);
pack(tile+3*HQ,s+1024*2048+1024,64,2048,Right);
}
// s4de_/s4df_: 横に隣接する 2 つの 64x64 tile（源 s と s+64）を 1 つの読み pass で作る。
// Right=false（A オペランド）= pack_wide（s4de_ が #119672 で着地）。
// Right=true （B オペランド）= pack_wide_b（s4df_、同機構の転置レイアウト版）。
template<bool Right>
static void pack_wide_root(U *tA,U *tB,const U *s) {
if constexpr(!Right) {
pack_wide(tA+0*HQ, tB+0*HQ, s);
pack_wide(tA+1*HQ, tB+1*HQ, s+1024);
pack_wide(tA+2*HQ, tB+2*HQ, s+1024L*2048);
pack_wide(tA+3*HQ, tB+3*HQ, s+1024L*2048+1024);
} else {
pack_wide_b(tA+0*HQ, tB+0*HQ, s);
pack_wide_b(tA+1*HQ, tB+1*HQ, s+1024);
pack_wide_b(tA+2*HQ, tB+2*HQ, s+1024L*2048);
pack_wide_b(tA+3*HQ, tB+3*HQ, s+1024L*2048+1024);
}
}
template<bool Right,bool High>
static void emit_root64(const U *tile,const U *left,const U *right,int offset,U *dest,U *destq,int voff) {
U *r0=destq+offset,*r1=r0+RS,*r2=r1+RS,*v7=dest+offset+voff*RS,*v5=v7+RS,*v6=v5+RS,*v34=v6+RS;
for(int i=0;i<(64*64+(64/BASE)*(64/BASE)*TILE_PAD);i+=32){
__m256i a=ld(tile+i),b=ld(tile+(64*64+(64/BASE)*(64/BASE)*TILE_PAD)+i),c=ld(tile+2*(64*64+(64/BASE)*(64/BASE)*TILE_PAD)+i),d=ld(tile+3*(64*64+(64/BASE)*(64/BASE)*TILE_PAD)+i);
__m256i aa=ld(tile+i+16),bb=ld(tile+(64*64+(64/BASE)*(64/BASE)*TILE_PAD)+i+16),cc=ld(tile+2*(64*64+(64/BASE)*(64/BASE)*TILE_PAD)+i+16),dd=ld(tile+3*(64*64+(64/BASE)*(64/BASE)*TILE_PAD)+i+16);
if constexpr(High) {
a=_mm256_add_epi16(a,_mm256_sub_epi16(ld(left+0*HQ+i),ld(right+0*HQ+i)));
b=_mm256_add_epi16(b,_mm256_sub_epi16(ld(left+1*HQ+i),ld(right+1*HQ+i)));
c=_mm256_add_epi16(c,_mm256_sub_epi16(ld(left+2*HQ+i),ld(right+2*HQ+i)));
d=_mm256_add_epi16(d,_mm256_sub_epi16(ld(left+3*HQ+i),ld(right+3*HQ+i)));
aa=_mm256_add_epi16(aa,_mm256_sub_epi16(ld(left+0*HQ+i+16),ld(right+0*HQ+i+16)));
bb=_mm256_add_epi16(bb,_mm256_sub_epi16(ld(left+1*HQ+i+16),ld(right+1*HQ+i+16)));
cc=_mm256_add_epi16(cc,_mm256_sub_epi16(ld(left+2*HQ+i+16),ld(right+2*HQ+i+16)));
dd=_mm256_add_epi16(dd,_mm256_sub_epi16(ld(left+3*HQ+i+16),ld(right+3*HQ+i+16)));
}
_mm256_stream_si256((__m256i*)(r0+i),a);_mm256_stream_si256((__m256i*)(r0+i+16),aa);
_mm256_stream_si256((__m256i*)(r1+i),Right?c:b);_mm256_stream_si256((__m256i*)(r1+i+16),Right?cc:bb);
_mm256_stream_si256((__m256i*)(r2+i),d);_mm256_stream_si256((__m256i*)(r2+i+16),dd);
__m256i x7=Right?_mm256_sub_epi16(d,b):_mm256_sub_epi16(a,c);
__m256i xx7=Right?_mm256_sub_epi16(dd,bb):_mm256_sub_epi16(aa,cc);
_mm256_stream_si256((__m256i*)(v7+i),x7);_mm256_stream_si256((__m256i*)(v7+i+16),xx7);
__m256i x5=Right?_mm256_sub_epi16(b,a):_mm256_add_epi16(c,d);
__m256i xx5=Right?_mm256_sub_epi16(bb,aa):_mm256_add_epi16(cc,dd);
_mm256_stream_si256((__m256i*)(v5+i),x5);_mm256_stream_si256((__m256i*)(v5+i+16),xx5);
__m256i x6=Right?_mm256_sub_epi16(d,x5):_mm256_sub_epi16(x5,a);
__m256i xx6=Right?_mm256_sub_epi16(dd,xx5):_mm256_sub_epi16(xx5,aa);
_mm256_stream_si256((__m256i*)(v6+i),x6);_mm256_stream_si256((__m256i*)(v6+i+16),xx6);
__m256i x34=Right?_mm256_sub_epi16(x6,c):_mm256_sub_epi16(b,x6);
__m256i xx34=Right?_mm256_sub_epi16(xx6,cc):_mm256_sub_epi16(bb,xx6);
_mm256_stream_si256((__m256i*)(v34+i),x34);_mm256_stream_si256((__m256i*)(v34+i+16),xx34);
}
}
template<bool Right>
static void pack_root(const U*s,int n,int offset,U *dest,U *destq,int voff) {
if(n>128){int h=n/2,q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
pack_root<Right>(s,h,offset,dest,destq,voff);
pack_root<Right>(s+h,h,offset+q,dest,destq,voff);
pack_root<Right>(s+h*2048,h,offset+2*q,dest,destq,voff);
pack_root<Right>(s+h*2048+h,h,offset+3*q,dest,destq,voff);
} else {
pack_wide_root<Right>(root_tile,high_left,s);
emit_root64<Right,false>(root_tile,nullptr,nullptr,offset,dest,destq,voff);
emit_root64<Right,false>(high_left,nullptr,nullptr,offset+HQ,dest,destq,voff);
pack_wide_root<Right>(high_right,root_tile,s+64*2048);
emit_root64<Right,false>(high_right,nullptr,nullptr,offset+2*HQ,dest,destq,voff);
emit_root64<Right,true>(root_tile,high_left,high_right,offset+3*HQ,dest,destq,voff);
}
}
void matrix_multiply(int n,const short*A,const short*B,short*C){
U *cq=(U*)C; // Three root input slots and the retired tail precede final C writeback.
U *ar=root_mem,*br=ar+4*RS,*pr=br+7*RS; // root_mem 15 -> 12 槽
U *work=cq+3*RS+(464*2);
raw_basis_input=(const U*)A;
pack_root<false>((const U*)A,1024,0,ar,cq,0);
raw_basis_input=(const U*)B;
pack_root<true>((const U*)B,1024,0,br,br,3);
_mm_mfence();
const U *a11=cq,*a12=cq+RS,*a22=cq+2*RS,*s7=ar,*s5=ar+RS,*s6=ar+2*RS,*s3=ar+3*RS;
const U *b11=br,*b21=br+RS,*b22=br+2*RS,*t7=br+3*RS,*t5=br+4*RS,*t6=br+5*RS,*t4=br+6*RS;
U *p1=ar+2*RS,*p2=br+5*RS,*p3=ar+RS,*p4=br+4*RS,*p5=ar,*p6=br+3*RS,*p7=pr;
mul<1024>(s7,t7,p7,work);mul<1024>(s5,t5,p5,work);mul<1024>(s6,t6,p6,work);
mul<1024>(s3,b22,p3,work);mul<1024>(a22,t4,p4,work);mul<1024>(a11,b11,p1,work);mul<1024>(a12,b21,p2,work);
U *out=(U*)C;
root_fuse<1024>(out,out+1024,out+1024*2048,out+1024*2048+1024,p1,p2,p3,p4,p5,p6,p7,2048);_mm_sfence();
asm volatile(".purgem vzeroupper");
_mm256_zeroupper();
}
// ---- ms2t_ tryrun driver (pricing only; NOT part of any submission) ----
#include <cstdio>
#include <cstdint>
alignas(4096) static uint16_t ms2t_A[2048*2048], ms2t_B[2048*2048], ms2t_C[2048*2048];
int main(){
for(long i=0;i<2048L*2048;i++){ms2t_A[i]=(uint16_t)(i*2654435761u);ms2t_B[i]=(uint16_t)(i*40503u);}
matrix_multiply(2048,(const short*)ms2t_A,(const short*)ms2t_B,(short*)ms2t_C);
uint64_t h=1469598103934665603ull;
for(long i=0;i<2048L*2048;i++){h^=ms2t_C[i];h*=1099511628211ull;}
printf("chk=%016llx\n",(unsigned long long)h);
return 0;
}
| Compilation | N/A | N/A | Compile Error | Score: N/A | 显示更多 |