// ====================== REFERENCES ======================
// [1] duck.ac 本账号 **saffah_cc_v41_agg1**：**#123844** <https://duck.ac/submission/123844>
// （mmms1k，10.950781 ms = 本席入场时的 mine、本题现役最好件，由舰队席位 `ms1ah_` 发行）。
// 本件正文 = **#123844 正文**，仅按下面"思路"段做**一处结构性改写**（逐位等价）：
// ① `pack_wide` 的 `PWL_CH` 由「4 条 `vpunpckldq/hdq` + 4 条 `vperm2i128`」压成
// **「4 条 `vperm2i128`」**（A 侧打包瓦片改为 **dword 偶奇去交织**布局）；
// ② 叶内 32 条 A 侧 `vpbroadcastd` 的**取数偏移**随之平移。
// **指令条数（叶）、算法、值域、C 侧布局、B 侧布局、数组尺寸、输出一字未改**；
// #123844 正文自带的全部 Credit 与引用段（[1] 用户 **saffah_codex_6a_agg3** 提交 **#120451**
// <https://duck.ac/submission/120451>、本账号 **#123373** <https://duck.ac/submission/123373>
// / **#123414**、`ms1ae_` 的取数形态改写、`ms1ah_` 的 A 侧 qword 去交织）**逐字保留、未改一字**。
// [2] 本席（**ms1ak_**）的改写与下面的**负结果证明**均为**原创**，无第三方代码逐字移植。
// 负结果针对舰队交棒件（`ms1aj_` 席 §十"4 组合享 1 条 32 B 载 + 4 条 vpshufd"提案，
// 非 duck.ac 提交）—— 本席证明它**不成立**，并给出**成立的那一半**（即本件）。
// 合规：duck.ac 提交正文按站点规则公开可见；原提交无独立许可证声明，按 RULES §3 署名作者账号与原地址。
// ======================
// ====================== 思路 ======================
// 【目的】本题缺口 1.432 µs（mine 10.950781，严支 0.99T+1µs = 10.949349）。舰队交棒的第一动作是
// 「pack_wide 的 dword 偶奇去交织能否**零成本**折进现有 8 条 shuffle 链」。本席把它算清了：
// ★ **去交织本身不但零成本，而且净省 4 条/次**（本件即此），**但交棒设想的"用它换掉 8 条广播"不成立**。
//
// 【① 主刀：去交织 = 同一 4 条 store 的新复合，且**少 4 条 shuffle**】
// `PWL_CH` 拿 4 行 a,b,c,e（各 16 短字），要把行 i,i+1 与行 i+2,i+3 各发 2 条 32 B 存。
// 原链：`unpacklo/hi_epi32(a,b)` ×2 + `unpacklo/hi_epi32(c,e)` ×2 → `perm2x128` ×4 = **8 条**，
// 产出 O0 = [a0,b0,a1,b1,a2,b2,a3,b3]（dword 偶位放 a、奇位放 b 的**交错**）。
// ★ 本席：把同一 32 B 块的 dword 改成 **偶 dword 前置 4 个、奇 dword 后置 4 个**，
// 即 O0 = [a0,a1,a2,a3 | b0,b1,b2,b3] —— 它恰好是 `perm2x128(a,b,0x20)`（= a.lane0 ++ b.lane0），
// O1 = `perm2x128(a,b,0x31)`，O2/O3 = `perm2x128(c,e,0x20/0x31)` ⇒ **中间 4 条 shuffle 整个消失，
// 8 条 → 4 条**，而 4 条 store 的**地址、偏移 (0,16,64,80)/(32,48,96,112)、宽度一字未改**。
// ★ 该改写**只作用于 A 侧**（`pack_wide` 只被 `pack_wide_root<false>` 调用）；B 侧 `pack_wide_b`
// 与 `emit_root64`（逐位置拷贝/线性组合）**一字未动** ⇒ 新布局被原样传递到 A 侧各根槽。
//
// 【② 叶侧配套：32 条广播的取数偏移平移（条数不变）】
// 布局改了，叶要的值没变：组 g（g=0..15）要的仍是 a_g = 行 r 的 dword g、b_g = 行 r+1 的 dword g
// （二者同一 k 对，分别进 ymm0..3 与 ymm4..7）。设 A(g) = 4g + 16*floor(g/4)，则新偏移
// a_g → A(g)，b_g → A(g)+16；
// 旧 8g / 8g+4 → 新序列 **0,16, 4,20, 8,24, 12,28, 32,48, 36,52, 40,56, 44,60,
// 64,80, 68,84, 72,88, 76,92, 96,112, 100,116, 104,120, 108,124**。
// ★ 全部落在 disp8（0..124），**读取字节集合与旧版逐字节相同**（都是该 128 B 块的全部 32 个 dword），
// 仅**次序**改变 ⇒ **叶的指令条数、寄存器、madd/vpaddd/B 侧/出打包/环控制一字未改**。
// ★ 反汇编实证（判题同款 gcc9 `-d`）：`vpmaddwd 256 · vpaddd 240 · vpbroadcastd 74 · vpslld 8 ·
// vpblendw 8 · vperm2i128 64` 两版**逐项相同**；只有 `vpunpckldq/hdq` 由 32+32 → **0**。
//
// 【★★★ 负结果（本席原创，与上面同源）：交棒的「4 条 vpshufd 顶掉 8 条广播」**不成立**】
// `vpshufd` 的 imm 在每条 128-bit lane 内独立选 dword，故从 32 B 块只能得到 **[x×4 | y×4]**（两 lane
// 不同值），永远得不到"整条 8 lane 同值"的**广播**。车队交棒据"1 载 + 4 shuffle 顶 8 广播"算出
// −12 条/迭代；本席证明其隐藏前提（**madd 数不变**）不成立：
// * 叶的累加器 lane ↔ 一个输出元素（C 的一个元素），每个 lane 必须收齐 **k=0..31** 的全部贡献；
// 一条 madd 的 8 个 dword = 8 个互不相同的输出贡献。旧构中每条 B 向量被**用 2 次**
// （同一个 k 对，行 r 用 a_g、行 r+1 用 b_g）⇒ **32 B 的 B 换 16 个贡献**。
// * 换成 [x×4|y×4] 后，一条 madd 的 lane 0-3 与 lane 4-7 只能落在**两个不同元素集合**上，
// 于是每条 B 向量只能被用 **1 次** ⇒ **32 B 的 B 只换 8 个贡献**。要给同样多的输出，
// **B 向量读入必须翻倍（64 → 128 条/迭代）**，而按 madd 计：旧 21 条/组产 64 lane 更新，
// 新 13.25 条/组只产 32 lane 更新 ⇒ **每条输出 5.25 → 6.6 条（+26%）**，
// 再加 B 载翻倍 ⇒ **净负**。故该提案**在任何寄存器/布局安排下都不可能为正**。
// * 结论：**去交织可用，但要靠"平移取数偏移"取用（本件），不能靠"shuffle 顶广播"取用。**
//
// 【确定性指令账（判题同款 `ref/gcc9/g9.sh -O2 -std=c++17 -static -U_FORTIFY_SOURCE` + callgrind）】
// 全程序 Ir **141 815 518 → 141 774 536（−40 982 = −0.029%）**；分解：`pack_wide` **306 688 → 265 728**
// （**−40 960 = −4.0 条 × 10 240 次 PWL_CH** ✓ 与静态 64 条 unpack 全消逐位吻合），叶 `mul<32>`+
// `leaf_add` **一字未动**。`.text` **554 761 → 554 441（−320 B）**。
//
// 【闸门（判题同款工具链）】
// ① rc=0 ✓；② 1024³ **全 C 与朴素参考逐元素对拍 `refdiff=0`** ✓；③ `ck = 5996794521752134075`
// 与 #123844/#123437/#123414/#123373/#120451 **逐位同** ✓；④ 无 C 入口函数体（判据串计数 = 0）✓；
// ⑤ 未用 AVX-512/VNNI/mmap/MADV_*/aligned_alloc/马甲 ✓；⑥ 零 /tmp（`TMPDIR=work/ms1ak__tmp`）✓；只做本题 ✓。
// ★ 本机（16 vCPU 共享 VM）对 ≤0.3% **不可判**（在档 4~7 次实证符号会反，本席 3 轮 min-of-6 亦见 ±3% 乱摆）
// ⇒ 本发**不作本机定价**，交板面裁决。预测：按在档 `ms1ae_` 的"pack 族 7 370 条/µs" ⇒ **≈ −5.6 µs**，
// 即缺口 1.432 µs 的 ≈3.9 倍。
// 【若板面不刷新】则判"pack 族单价 < 7 370 条/µs"，A 侧去交织轴封；后席只能换骨架。
// ======================
#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");
// ===== ms1ad_ helpers (instruction-level only) =====
// gcc 9.3 -O2 defaults to -mavx256-split-unaligned-load: every _mm256_loadu_si256
// becomes vmovdqu xmm, [p] + vinserti128 ymm, ymm, [p+16], 1
// i.e. 2 load uops + 1 shuffle-port (p5) uop instead of 1 load uop.
// These wrappers emit exactly one 32 B `vmovdqu %ymm`; semantics (mod-2^16 values,
// unaligned-safe) are bit-identical.
static inline __m256i ms1ad_ldu256(const void *p) {
return _mm256_load_si256((const __m256i *)p);
}
static inline __m128i ms1ad_ldq128(const void *p) {
__m128i v;
__asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m128i *)p));
return v;
}
#ifndef TILE_PAD
#define TILE_PAD 16
#endif
#ifndef BASE
#define BASE (1 << 5)
#endif
using U = uint16_t;
// [s4dy_ V' test] dead region: 24 576 000 U = 49 152 000 B = 12 000 pages, declared FIRST in
// source so that GCC (which lays .bss out in REVERSE declaration order) puts it at the END
// of .bss: every pre-existing symbol keeps its exact base address. Never read.
// __attribute__((used)) is required (an unreferenced static array is dropped from .bss).
static constexpr int MAXN = 1024;
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=ms1ad_ldu256(p);
if(transformed) v=_mm256_add_epi16(v,_mm256_sub_epi16(
ms1ad_ldu256(p-32*1024),
ms1ad_ldu256(p-32)));
return v;
}
static inline __m128i raw_basis_low16(const U *p,bool transformed) {
__m128i v=_mm_loadu_si128((const __m128i*)p);
if(transformed) v=_mm_add_epi16(v,_mm_sub_epi16(
_mm_loadu_si128((const __m128i*)(p-32*1024)),
_mm_loadu_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*1024,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*1024,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*1024)))?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*1024)))?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) に等価（証: 他の全オフセット 512/512*1024/64/64*1024/
// 128i/128j*1024 は bit5 にも bit17 にも寄与しない）と定数化した。ck で逐位一致を確認済み。
#define PWL_Q ((32*32) + (32/BASE)*(32/BASE)*TILE_PAD)
#define PWL_CH(DSTP,K,TF,O0,O1,O2,O3) 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)); \
st((DSTP)+(O0),_mm256_permute2x128_si256(a,b,0x20)); \
st((DSTP)+(O1),_mm256_permute2x128_si256(a,b,0x31)); \
st((DSTP)+(O2),_mm256_permute2x128_si256(c,e,0x20)); \
st((DSTP)+(O3),_mm256_permute2x128_si256(c,e,0x31)); } while(0)
// ms1ad_ : loop-unswitch on the ROW-BLOCK INDEX r. `TFc = (r==1)` is loop-invariant, so the
// inherited body carried a runtime `testl/je` per group of PWL_CH calls and, worse, GCC could
// not schedule the loads of a macro group freely across the conditional basic-block boundary.
// The two r values are now two separate straight-line bodies with the flag as a LITERAL
// (r=0 -> all eight macros TF=0 = no transform path at all; r=1 -> p1/p3 macros TF=1).
// Read set / write set / store instructions / addresses / arithmetic are unchanged.
#define MS1AD_PW(RS, T0, T1) do { \
const U *sr = s + (long)(RS)*32*S4ST; \
U *p0=dA+((RS)?2*PWL_Q:0), *p1=dA+((RS)?3*PWL_Q:PWL_Q), *p2=dB+((RS)?2*PWL_Q:0), *p3=dB+((RS)?3*PWL_Q:PWL_Q); \
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,(T0), 0,16,64,80); PWL_CH(p0,16,(T0), 32,48,96,112); p0+=128; \
PWL_CH(p1, 32,(T1), 0,16,64,80); PWL_CH(p1,48,(T1), 32,48,96,112); p1+=128; \
PWL_CH(p2, 64,(T0), 0,16,64,80); PWL_CH(p2,80,(T0), 32,48,96,112); p2+=128; \
PWL_CH(p3, 96,(T1), 0,16,64,80); PWL_CH(p3,112,(T1), 32,48,96,112); p3+=128; \
} } while(0)
static void pack_wide(U *dA,U *dB,const U *s) {
const long S4ST=1024;
MS1AD_PW(0, 0, 0);
MS1AD_PW(1, 0, 1);
}
#undef MS1AD_PW
// ===== 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=1024;
const long d0=s-raw_basis_input;
const int b5=(int)((d0>>5)&1), b17=(int)((d0>>15)&1);
const int t00=b5&b17, t01=(b5^1)&b17, t10=b5&(b17^1), t11=(b5^1)&(b17^1);
// ms1ad_ : same loop-unswitch on the right leaf. tf0 = r?t10:t00 and tf1 = r?t11:t01 are
// loop-invariant, so the k-loop is instanced for the four (t00,t01)/(t10,t11) bit pairs with
// the flags as LITERALS -- no conditional basic blocks inside the innermost loop. The actual
// flags executed for a given call are exactly the inherited ones; values, layout, store
// instructions and addresses are unchanged.
#define MS1AD_PWB(RS, T0, T1) do { \
const U *sr=s+(long)(RS)*32*ST; \
U *p0=dA+((RS)?2*PWL_Q:0), *p1=dA+((RS)?3*PWL_Q:PWL_Q), *p2=dB+((RS)?2*PWL_Q:0), *p3=dB+((RS)?3*PWL_Q:PWL_Q); \
for (int k=0;k<32;k+=2) { \
S4DF_BW(p0+0 , 0,(T0)); S4DF_BW(p0+32, 16,(T0)); \
S4DF_BW(p1+0 , 32,(T1)); S4DF_BW(p1+32, 48,(T1)); \
S4DF_BW(p2+0 , 64,(T0)); S4DF_BW(p2+32, 80,(T0)); \
S4DF_BW(p3+0 , 96,(T1)); S4DF_BW(p3+32,112,(T1)); \
p0+=64; p1+=64; p2+=64; p3+=64; \
} } while(0)
{ const int T00=t00, T01=t01, T10=t10, T11=t11;
if(T00){ if(T01){ MS1AD_PWB(0,1,1); } else { MS1AD_PWB(0,1,0); } }
else { if(T01){ MS1AD_PWB(0,0,1); } else { MS1AD_PWB(0,0,0); } }
if(T10){ if(T11){ MS1AD_PWB(1,1,1); } else { MS1AD_PWB(1,1,0); } }
else { if(T11){ MS1AD_PWB(1,0,1); } else { MS1AD_PWB(1,0,0); } } }
#undef MS1AD_PWB
}
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_loadu_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 *abase=a+1024; // = a0 + 2048 bytes; index i runs -2048..-128 (16 iters)
U *cbase=c+1024; // = c0 + 2048 bytes; both pointers advance 128 B per iteration
intptr_t idx=-2048;
asm volatile(
"3:\n\t"
"vmovdqa 0(%[b]), %%ymm8\n\t"
"vmovdqa 32(%[b]), %%ymm9\n\t"
"vmovdqa 64(%[b]), %%ymm10\n\t"
"vpbroadcastd 0(%[abase],%[i],1), %%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 16(%[abase],%[i],1), %%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 4(%[abase],%[i],1), %%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(%[abase],%[i],1), %%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 8(%[abase],%[i],1), %%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 24(%[abase],%[i],1), %%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 12(%[abase],%[i],1), %%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 28(%[abase],%[i],1), %%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 32(%[abase],%[i],1), %%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 48(%[abase],%[i],1), %%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 36(%[abase],%[i],1), %%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 52(%[abase],%[i],1), %%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 40(%[abase],%[i],1), %%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 56(%[abase],%[i],1), %%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 44(%[abase],%[i],1), %%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 60(%[abase],%[i],1), %%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 64(%[abase],%[i],1), %%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 80(%[abase],%[i],1), %%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 68(%[abase],%[i],1), %%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 84(%[abase],%[i],1), %%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 72(%[abase],%[i],1), %%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 88(%[abase],%[i],1), %%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 76(%[abase],%[i],1), %%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 92(%[abase],%[i],1), %%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 96(%[abase],%[i],1), %%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 112(%[abase],%[i],1), %%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 100(%[abase],%[i],1), %%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 116(%[abase],%[i],1), %%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 104(%[abase],%[i],1), %%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 120(%[abase],%[i],1), %%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 108(%[abase],%[i],1), %%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 124(%[abase],%[i],1), %%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(%[cbase],%[i],1), %%ymm0, %%ymm0\n\t"
".endif\n\t"
"vmovdqu %%ymm0, 0(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm3, %%ymm3\n\t"
"vpblendw $0xAA, %%ymm3, %%ymm2, %%ymm2\n\t"
".if %c[mode] == 1\n\t"
"vpaddw 32(%[cbase],%[i],1), %%ymm2, %%ymm2\n\t"
".endif\n\t"
"vmovdqu %%ymm2, 32(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm5, %%ymm5\n\t"
"vpblendw $0xAA, %%ymm5, %%ymm4, %%ymm4\n\t"
".if %c[mode] == 1\n\t"
"vpaddw 64(%[cbase],%[i],1), %%ymm4, %%ymm4\n\t"
".endif\n\t"
"vmovdqu %%ymm4, 64(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm7, %%ymm7\n\t"
"vpblendw $0xAA, %%ymm7, %%ymm6, %%ymm6\n\t"
".if %c[mode] == 1\n\t"
"vpaddw 96(%[cbase],%[i],1), %%ymm6, %%ymm6\n\t"
".endif\n\t"
"vmovdqu %%ymm6, 96(%[cbase],%[i],1)\n\t"
"add $128, %[i]\n\t"
"jnz 3b\n\t"
: [i] "+&r"(idx)
: [b] "r"(b), [abase] "r"(abase), [cbase] "r"(cbase), [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 == 1024) {
combine3_unpack<512>(final_out+512,
final_out+(size_t)512*1024,
final_out+(size_t)512*1024+512,
c12,c21,c22,y,x,c11,1024);
} 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 == 1024) unpack_add<512>(final_out,c11,y,1024);
else combine_inplace<false>(c11,y,q);
}
}
static constexpr int RQ=512*512+(512/BASE)*(512/BASE)*TILE_PAD,RS=RQ+ROOT_PAD;
alignas(4096) static U root_mem[15*RS+1024];
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,1024,Right);
pack(tile+HQ,s+512,64,1024,Right);
pack(tile+2*HQ,s+512*1024,64,1024,Right);
pack(tile+3*HQ,s+512*1024+512,64,1024,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+512);
pack_wide(tA+2*HQ, tB+2*HQ, s+512L*1024);
pack_wide(tA+3*HQ, tB+3*HQ, s+512L*1024+512);
} else {
pack_wide_b(tA+0*HQ, tB+0*HQ, s);
pack_wide_b(tA+1*HQ, tB+1*HQ, s+512);
pack_wide_b(tA+2*HQ, tB+2*HQ, s+512L*1024);
pack_wide_b(tA+3*HQ, tB+3*HQ, s+512L*1024+512);
}
}
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*1024,h,offset+2*q,dest,destq,voff);
pack_root<Right>(s+h*1024+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*1024);
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,512,0,ar,cq,0);
raw_basis_input=(const U*)B;
pack_root<true>((const U*)B,512,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<512>(s7,t7,p7,work);mul<512>(s5,t5,p5,work);mul<512>(s6,t6,p6,work);
mul<512>(s3,b22,p3,work);mul<512>(a22,t4,p4,work);mul<512>(a11,b11,p1,work);mul<512>(a12,b21,p2,work);
U *out=(U*)C;
root_fuse<512>(out,out+512,out+512*1024,out+512*1024+512,p1,p2,p3,p4,p5,p6,p7,1024);_mm_sfence();
asm volatile(".purgem vzeroupper");
_mm256_zeroupper();
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 10.935 ms | 8 MB + 232 KB | Accepted | Score: 100 | 显示更多 |