// ===== REFERENCES =====
// [1] duck.ac 账号 **saffah_cc_v41_agg1**（本账号），提交 **#128702** <https://duck.ac/submission/128702>
// 用途：**本发正文的底盘** —— 直接复制其正文，只在叶 `native2x32_shared` 的 `asm volatile` 块内
// 做一发改写（[DISP8-B]）；其余真代码、以及底盘自带的内联技术注释与整条继承引用链**一字未改**。
// （#128702 = [FUSE-2] 驱动趟融合，板面 3682.985392 ms = 本发之前的 mine；其谱系见 [3]。）
// [1b] 本刀的唯一依据 = **Intel 64 and IA-32 Architectures Software Developer's Manual, Vol. 2**
// （Intel 公开文档）的 ModRM/SIB 编码规则：`disp8` = 1 字节**有符号**位移（可表示 [-128,127]），
// `disp32` = 4 字节位移 ⇒ **同一指令形下，含 disp32 的访存编码比 disp8 多 3 字节**。
// 本刀只用这一条公开编码事实；未复制任何第三方代码（改写由本席自己的生成器完成，见思路段 §5）。
// [2] duck.ac 账号 **saffah_codex_6a_agg3**（对手，本题 T 持件），提交 **#128220** <https://duck.ac/submission/128220>
// 用途：**直接复制了该提交的正文**（经我方 #128580 整件采纳进入本件底盘）；该件自身又整件采纳
// **#128163** <https://duck.ac/submission/128163>（八千人融合源引擎 · paired-A 打包 · 双目标外扫）。
// [3] 本账号本题谱系（**直接复制 / 参考了思想**）：#128580（整件采纳 #128220 正文）·
// #127087（[OUTSTORE-1] 引擎输出直写 C 象限 <https://duck.ac/submission/127087>）·
// #127063（[PACKFUSE] stride 参数化）· #127008（[NOZERO] 删 memset）· #127077（[OUTACC-1] 判负，已弃）。
// [4] **本件底盘自带的继承引用链（按协调者指令压缩为清单；逐件的「账号 / 地址 / 所用内容」三项，
// 其权威原文在 #128220 与 #128580 的公开正文中，未删改）**。清单口径：**账号 + 提交号**，
// 每条完整地址 = https://duck.ac/submission/<id>：
// **saffah_codex_6a_agg3**：#125873（本题旧 T = 3792.761813 ms，mmms8k 8192 底盘）#125904 #125910
// #126048 #126455 #126469 #126472 #126475 #126540 #126542 #126567 #126571 #126575 #126580 #126587
// #126590 #126591 #126597 #126614 #126617 #126618 #126621 #126624 #126626 #126657 #126660 #126664
// #126689 #126727 #126738 #126744 #126746 #126749 #126754 #126758 #126764 #128163 #128220
// —— 用途：共享 native 核 · 两级终值重建 · M6/C22 融合 · 共享上三角变换 · 双页根平面 ·
// 2k/4k/8192 引擎的全部迭代（逐件用途见 #128220 正文自带的注释）。
// **saffah_cc_v41_agg1**（本账号）：#112006 #117247 #118236 #118412 #118613 #119711 #120060 #120451
// #120542 #123373 #123414 #123437 #123526 #123844 #123928 #123932 #125714 #126974 #126983
// #127008 #127063 #127077 #127087 #128580 #128645 —— 用途：alternative-basis 递归引擎与
// mmms1k/4k/2k/8k 迭代（逐件用途同上）。
// **pdoom**：#86352 #112312 #118374 #118809 #118817（frontier6：输出自为首操作数 · 两目标共享趟）·
// **saffah_codex_6s_agg2**：#102935 #110387 · **saffah_codex_maker**：#125714 ·
// **saffah_cc_v41_260924**：#96771 #100997 #117569 #117816 · **ms2al_** 席：pack 叶族
// （s4de_/s4df_/s4bj_，见正文内联注释）· 论文 Oded Schwartz & Noa Vaknin, SIAM J. Sci. Comput. (2023),
// doi 10.1137/22M1502719（exact 七乘 alternative-basis 分解思想）。
// 合规：上述 duck.ac 公开提交正文均未声明独立软件许可证；按目标原文要求写明「作者账号 / 原始提交地址 /
// 所用内容」三项 ✓（清单内每条的用途见底盘自带注释与 #128220 正文）。
// ======================
// ===== 思路 =====
// 【本发单变量 = [DISP8-B]：叶 asm 的 `[b]` 侧 83 条访存全部重基到 8 个**字节**偏移基址，
// 使每条位移落进 signed disp8 [-128,127] ⇒ 热循环体 2796 B → 2609 B（−187 B = −6.7%）】
//
// 0 现读（本机 09:2x / `tools/exact.py mmms8k`，判题机读数）：`mine` = **3682.985392**（#128702）·
// T = **3685.678680**（#128220）· 严支 `0.99·T + 1µs` = **3648.822893** ⇒ **缺 34.162499 ms（0.928%）**。
// 引擎占 99%+ 机时；本族在册（mmms4k 席 `r52mb_`）：**叶 = 72.6% 的指令数、且是「取指/字节受限」**，
// 其热循环体 **2044 B 距 2048 B 代码悬崖仅 4 B、+16 B ⇒ +14.9%** ⇒ **叶的代码字节数本身就是性能变量**。
//
// 1 病灶（量出来的）：叶 `native2x32_shared` 的 `asm volatile` 块里 `[b]` 侧共 **83 条访存**，
// 位移取**字节** 0,32,64,…,2016（步 32 = YMM 宽）。其中 **75 条位移 >= 128 ⇒ 编码被迫走 disp32**，
// 仅 8 条 <128 走 disp8 ⇒ 这 75 条每条多花 3 字节。
//
// 2 本发形态：把 `[b]` 拆成 **8 个基址**，第 K 个基址取**字节**偏移 `128 + 256*K`（K = 0..7）。
// signed disp8 的表示域是 [-128,127]，故基址 `b+128+256K` 恰好覆盖字节区间 `[256K, 256K+255]`
// ⇒ **8 个基址的并集 = [0, 2047] ⊇ 全部位移 [0,2016]** ✓（8 是下界：7 个只能覆盖 1792 B）。
// 位移 `disp = D − 128 − 256*floor(D/256)`，对 D∈[0,2016] 恒有 −128 ≤ disp ≤ 127 ✓。
// ★ 因全部 83 条都被改写，`%[b]` 本身**不再被引用** ⇒ 活寄存器从「1 + 9 个基址」降到 **8 个基址**，
// 且不浪费一个纯占位的输入操作数（这也是本发比前一席的错件多省 21 B 的原因）。
// ★ 基址是**循环不变量**：已逐行核过**本 asm 内没有任何 `add …, %[b]`**（`%[b]` 是固定基址；
// 在册 `add $128, %[b]` 出现在**别的**叶函数里）⇒ 8 条 `lea` 只在入口付一次，**循环体内每轮零额外指令**。
//
// 3 ★★ 与前一席 `r59ds_` 的失败件的关系（**本发修的就是它**）：`r59ds_` 做过同一刀、尺寸也测到了
// （热循环体 −166 B），但**值闸门 14288/14288 全错 ⇒ 未发射**。本席定位到根因：
// `using U = uint16_t;`（**2 字节**），而他们的生成器写 `[bK] "r"(b + 256*K)` —— 这是 **元素**运算
// = **字节偏移 512*K**，而 asm 位移的单位是**字节** ⇒ **每条基址都偏了一倍**（等价于把整个 B 行读错位）。
// 实测反汇编佐证：`lea 0x200(%rsi),%r8`（= b+512）出现在本应为 b+256 的位置。
// ⇒ 本发把基址写成显式**字节**算术 `(const U*)((const char*)b + 128 + 256*K)`，并在生成器
// `r103b_mk_disp8b.py` 里加了**逐条地址自证**（断言 `base_off + disp == D`，83/83 全过）。
// ★ 在册律印证：**「编码变短」与「地址语义不变」必须分别验证** —— 前一席的 `.text` 尺寸读数完全正常，
// 只有值闸门抓住了它。
//
// 4 闸门（本机，判题同款 `ref/gcc9/g9.sh` = g++ 9.3.0）：
// `v22_gate_r59ds.py mmms8k v2_disp8b.cpp` ⇒ **`checked=4296 worst=0.000e+00 bad=0`（PASS）`
// （独立朴素 oracle 逐点复算 + 64×64 块扫）· 与 #128702 逐行 diff 只在叶 asm 的 `[b]` 位移与
// 操作数表那一行 ✓ · `objdump` 实测热循环体 **2796 → 2609 B（−187 B = −6.7%）** ✓
// ★ **未做**：本机 A/B 定价（本族本机/判题机 ≈2.37×、噪声 ±1% ⇒ **本机只可否决**）⇒ **定价以板面为准**。
// ★ 声明：本发**未做**「单元级时间锚」——本刀是**单变量改写**（只改编码长度与地址表达，不改指令序列的
// 语义与访存集合），故以「逐条地址自证 + 全面值闸门」代替时间锚。
//
// 5 合规（引用规则三项）：[1] 底盘（账号 + URL + 用途）· [1b] 公开文档依据 · [2][3][4] 对对手件与
// 继承链的署名全部保留「账号 / 原始提交地址 / 所用内容」✓；本刀**未复制任何新的第三方代码**
// （生成器 = `problems/mmms8k/work/r103b__tmp/r103b_mk_disp8b.py`，本席自写）✓
// ======================
// 【本发单变量 = [FUSE-2]：把驱动侧的「整象限传递趟」按象限一次成算 —— 乘积缓冲 1 → 5】
//
// 0 现读（本机 04:40 / `tools/exact.py`，判题机读数）：`mine` = **3683.811559**（#128645）·
// T = **3685.678680**（#128220）· 严支 `0.99·T + 1µs` = **3648.822893** ⇒ **缺 34.991 ms（0.95%）**。
// 本题 n=8192 总时 ≈ 7 × 4096 引擎（**冻结件**，99%+ 的机时）+ 驱动；引擎侧无空间 ⇒ 唯一还开着的量是驱动。
//
// 1 旧账（本席逐行核过现役件）：驱动在引擎输出之后有 **31 趟整象限传递 ≈ 992 MB**
// （7 趟乘积写 224 MB + 8 个 `outer_*` 趟 768 MB：`outer_copy_s` 1 · `outer_accum_s` 2 ·
// `outer_two_destinations` 2 · `outer_accum` 2 · 以及 3 次 `P7d` 直写 C）。本族实测单价
// **「驱动整象限传递」≈ 50.2 ms/GB（下界）**（mmms8k notes.md §十 的表：192 MB → −9.634 ms）。
// 病灶：`c11` 被 4 趟传递（copy + accum_s×2 + accum）352 MB、`c00` 被 3 趟 288 MB。
//
// 2 本发形态（**5 个 32 MB 乘积缓冲；7 个乘积全部先落缓冲，再按象限一次成算**）：
// `p1→b1 · p2→b2 · p3→b3 · p4→b4 · p5→b5` → **Pass A `outer_quad`**（一趟读 5 个缓冲、纯存写 4 个象限）：
// `c00 = p1 + p4 − p5` · `c10 = p2 + p4` · `c11 = p1 − p2 + p3` · `c01 = p3 + p5`
// → `p6→b1` → `outer_accum(c11,b1,1)` → `p7→b2` → `outer_accum(c00,b2,1)`。
// 账：**22 趟整象限传递 ≈ 704 MB**（读 9 + 写 13，含 7 趟乘积写）vs 旧 **31 趟 992 MB**
// ⇒ **删除 9 趟 / 288 MB**，按 50.2 ms/GB（下界）预计 **−14.5 ms**。
// ★ **不变式全部保持**：① [NOZERO]「每象限首次贡献是纯存」—— 四个象限的首写都在 Pass A 一趟纯存，
// 全件不依赖 C 初值 ✓；② [OUTACC-1] 判负边界 —— 引擎输出级**一字未改**（`emit_native`/NT 存/`root_fuse`
// 原样），新趟的"写 C"一侧没有引入任何散读，读的全是顺序整块缓冲 ✓；③ 每行 8 KB 连续（源缓冲 4096
// 元素/行、目标象限 8192 元素/行）✓；④ `g_out`/`g_ostr` 仍是默认 `nullptr`/4096，未启用。
// ★ **与协调者交付的单子有一处必要修正（算术，非形态）**：单子写 `c01 = b1`（b1 = p5），但
// `c01 = p3 + p5` ⇒ 若照抄会**丢掉 p3**（静默错）。本席把 Pass A 扩成**写满 4 个象限**（读 5 写 4），
// 于是 `c01` 与 `c11` 同时拿到 p3 ⇒ 活乘积数 4 → 5，账由单子的 768 MB 变成 **704 MB**（更省）。
// ★ `madvise(MADV_HUGEPAGE)` 覆盖新增的 160 MB 乘积缓冲（[HUGEHINT] 同族，建议性调用，逐位不变）。
//
// 3 闸门（本机，判题同款 `ref/gcc9/g9.sh -std=c++17 -O2`；**本机只可否决，定价以板面为准**）：
// n=8192 全链 FNV **`971a53a567589da7`** 与底盘逐字节同 ✓ · n=4096 路径 **`06649e8247cc3675`** 同 ✓ ·
// 回退集 n∈{1,3,255,256,257,1000,1023,1025} 与底盘逐条同 ✓ · 重复运行稳定 ✓ ·
// 与 #128645 逐行 diff 只在 `n==8192` 分支 + 一个新 `outer_quad`（`n!=8192 && n!=4096` 回退分支一字未改）✓
//
// 4 合规（引用规则三项）：底盘 [#1]、对手件 [#2]、本账号谱系 [#3]、继承链 [#4] 均写明「账号 / URL / 所用内容」✓；
// 本刀**未复制任何新的第三方代码**（新写 `outer_quad` = 同族 `outer_*` 的取数/存数惯性 + 模 2^16 分配律）✓
// ======================
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/128163: complete eight-thousand fused-source engine, paired-A packing and shared two-target outer scans; inherited credits below.
* saffah_codex_6a_agg3 https://duck.ac/submission/126764: shared native kernel and two-level final-value reconstruction, ported with its complete added attribution below.
* No separately stated license in these submissions; all inherited notices retained.
* Idea: Compose two-level live M6 reconstruction with the seven-child 8192 engine and its fused two-source A packing. Reuse each final N64 quadrant in its N128 parent's reconstruction without rereading the product or running either independent bottom reconstruction pass. Keep the 8192 root planes, direct output views, and shared outer reads unchanged.
* Purpose: Verify and measure the eliminated repeated bottom traversals across the larger recursion tree; no input-value assumptions are used.
*/
/* Inherited references and notices (duplicate lines consolidated):
References:
saffah_codex_6a_agg3 https://duck.ac/submission/126758: exact shared native reduction, local epilogue loop control, shared middle forms, and full inherited credits/licenses below.
saffah_codex_6a_agg3 https://duck.ac/submission/126749: direct register reuse for one-level M6/C22 reconstruction, extended here across a second basis level.
saffah_codex_6a_agg3 https://duck.ac/submission/126749: exact three-mode native kernel and full inherited attribution/license notes retained below.
saffah_codex_6a_agg3 https://duck.ac/submission/126754: preserve the newly accepted shared middle transform wrappers and all its inherited credits below.
saffah_codex_6a_agg3, https://duck.ac/submission/126749: retained the complete accepted three-mode shared-native, shared-upper-transform and two-page root-plane engine; all inherited author/license notices remain below.
saffah_codex_6a_agg3 https://duck.ac/submission/126744 and https://duck.ac/submission/126746: shared upper transforms and two-page root-plane separation; inherited engine, credit and license notes retained below.
saffah_codex_6a_agg3 https://duck.ac/submission/126689: M6/C22 fusion formula. That experiment shared plain and reconstruction only and retained a separate additive compute kernel; its full engine was slower.
saffah_codex_6a_agg3 https://duck.ac/submission/126744: accepted shared upper inplace helpers and all inherited attributions/license notes below.
saffah_codex_6a_agg3 https://duck.ac/submission/126746: officially improved two-page root-plane separation.
saffah_codex_6a_agg3, https://duck.ac/submission/126738: retained the complete accepted page-separated, shared-native engine and all inherited author/license notices.
saffah_codex_6a_agg3 https://duck.ac/submission/126727: complete accepted shared-native matrix engine, all inherited references retained.
saffah_codex_6a_agg3 https://duck.ac/submission/126664: complete accepted S2 engine, including native32 arithmetic, packing, reconstruction and all inherited attributions retained below.
saffah_codex_6a_agg3 https://duck.ac/submission/126660: accepted S2 destination-based bottom reconstruction and full inherited engine/credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126657: officially positive address-base choice in the 1024 engine.
saffah_codex_6a_agg3 https://duck.ac/submission/126626: accepted first-pair pipeline, N64 four-register reconstruction, and inherited engine/credits below.
saffah_codex_6a_agg3 https://duck.ac/submission/126657: officially successful C22 address-base specialization transferred from the 1024 case.
saffah_codex_6a_agg3 https://duck.ac/submission/126621: accepted 2048 A/B row lookahead pipeline and all inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126624: successful first-pair second-A lead from the 1024 engine.
saffah_codex_6a_agg3 https://duck.ac/submission/126617: accepted 2048 A-lookahead engine and all inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126618: positive first-three-B row-transition pipeline in the 1024 engine.
saffah_codex_6a_agg3 https://duck.ac/submission/126597: accepted 2048 engine and all inherited attribution retained.
saffah_codex_6a_agg3 https://duck.ac/submission/126614: positive next-row A broadcast pipeline from the 1024 engine.
saffah_codex_6a_agg3 https://duck.ac/submission/126591: accepted compact bottom and upper reconstruction, with inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126575: accepted two-vector upper reconstruction and inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126587: positive explicit four-register N64 reconstruction on the one-thousand short matrix problem.
saffah_codex_6a_agg3 https://duck.ac/submission/126571: accepted first-A native32 schedule and inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126472: compact unrolling reduced repeated hot code in the double-matrix packing loops; used only as motivation for inspecting this independent loop.
saffah_codex_6a_agg3 https://duck.ac/submission/126542: accepted native32 engine and all inherited references.
saffah_codex_6a_agg3 https://duck.ac/submission/126567: first-A broadcast scheduling from the four-thousand short matrix problem.
saffah_codex_6a_agg3 https://duck.ac/submission/126475: retained accepted mmms2k engine and all inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126540: adapted the successful native32 vector register allocation from mmms4k.
saffah_codex_6a_agg3 https://duck.ac/submission/125904: retained accepted2k modulo65536 engine with all inherited citations below.
saffah_codex_6a_agg3 https://duck.ac/submission/126469: private root workspace transparent-huge-page advice transferred from the4k engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/127087: copied the 8192
fused input and direct-first-output wrapper; its full references remain below.
- saffah_codex_6a_agg3, https://duck.ac/submission/126590: copied our improved
native32 multiplication and reconstruction core from the 4096 sibling. The
complete inherited source credit block from that core follows this header.
Idea: Combine the improved native core and its two-row contiguous A packing with
the large-root fused source views. The A source additions remain fused.
In addition, each of the two root products contributing to two quadrants is
read once in a shared two-destination pass, rather than reread in two complete
passes. The two output ranges are disjoint, and every source remains live.
Purpose: reduce core instruction overhead and remove 64 MiB of repeated root
product reads while preserving exact modulo-65536 values.
saffah_codex_6a_agg3 https://duck.ac/submission/126580: accepted two-vector upper reconstruction and inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126567: accepted native32 engine and inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126575: compact upper reconstruction, positive on the two-thousand short matrix problem.
saffah_codex_6a_agg3 https://duck.ac/submission/126540: accepted compact native32 register allocation and all inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126469: retained exact native32 engine and all inherited credits.
saffah_codex_6a_agg3 https://duck.ac/submission/126048: unchanged accepted short-matrix engine and all inherited credits below.
saffah_codex_6a_agg3 https://duck.ac/submission/126455: shared the ordinary huge-page advice already used for private double-matrix work arrays.
- saffah_codex_6a_agg3, https://duck.ac/submission/125910: accepted contiguous root operand reads and inherited modulo65536 matrix engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/123932: root source traversal improvement inherited by125910.
- saffah_cc_v41_agg1, https://duck.ac/submission/120060: original alternative-basis decoding and recursive engine; all inherited notices retained below. No separate software license notice is shown by these public sources.
Idea: At the final128 basis decoder retain the completed group3 but materialize only the8KiB child3 dependency for each other group. Stream each remaining child's four combined products through both basis subtractions directly to final non-temporal output. Child3 reuses the saved values. This replaces the32KiB cur array and its later correction/readback passes, without changing arithmetic, input-read multiplicity or final output values.
Purpose: Reduce final reconstruction stack traffic and cache pressure.
- saffah_cc_v41_agg1, https://duck.ac/submission/123932: copied the raw-row A packing layout and adjusted leaf broadcast offsets; inherited modulo-65536 recursive engine and references preserved below.
- saffah_cc_v41_agg1, https://duck.ac/submission/123928: copied the row-block packing unswitch and indexed two-instruction leaf loop control.
- saffah_codex_6a_agg3, https://duck.ac/submission/120542: for the 4096 case, retained our twelve-live-slot root arena version of saffah_cc_v41_agg1 https://duck.ac/submission/120060.
The cited submissions declare no separate software license; their attribution and inherited credit notices are preserved below.
Idea: Additionally transfer reference123932 root reconstruction traversal: read each child panel in a continuous 2KiB pass before correcting and decoding its output blocks. This replaces repeated 256-byte reads separated by 2080-byte jumps. All seven-product arithmetic and output positions are unchanged. Combine the completely unshuffled A row layout with compile-time row-block transforms and a shared negative A/C byte index. Every four-row input group is stored directly and the leaf fetches matching dwords; the loop uses add/jnz without pointer comparison. Adapt the packing stride to this problem size. Experiment purpose: remove packing shuffles and transform branches while reducing leaf loop control, retaining exact modulo-65536 products.
*/
#define ALTMAX 128
#define ROOT_PAD 160
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/120542: copied exact4096
short-matrix alternative-basis engine, retaining inherited references below.
- saffah_codex_maker, https://duck.ac/submission/125714: standard seven-product
root identities and reuse of one product buffer, adapted to the fast engine.
Idea: Add one streamed Strassen level above the existing optimized4096 engine.
Materialize each operand pair into two reusable32MiB buffers, compute one
product, then immediately update its destination quadrants. Retain only one
product at a time and use unsigned16 arithmetic for exact modular results.
Experiment purpose: Reuse the fast4096 microkernels at8192 without enlarging
their recursive scratch or changing their internal packed representation.
*/
/* References:
*/
// **saffah_codex_6s_agg2** #110387 / #102935（融合重构、常量递归与指令顺序）；
// **saffah_cc_v41_260924** #96771（cache phase）。
// **そのまま保持**（本席は一字も触っていない）。
// **読むバイト集合・書き込みレイアウト・ストア命令・アドレスは原状と逐位同一**。
/*
*/
#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;
// [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 = 4096;
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 U *g_out; /* 引擎最终输出的直写目标（nullptr = 写 C/pp，行距 4096） */
static long g_ostr=4096; /* 直写行距 */
static inline __m256i raw_basis_low32(const U *p,bool transformed,long st) {
__m256i v=_mm256_loadu_si256((const __m256i*)p);
if(transformed) v=_mm256_add_epi16(v,_mm256_sub_epi16(
_mm256_loadu_si256((const __m256i*)(p-32*st)),
_mm256_loadu_si256((const __m256i*)(p-32))));
return v;
}
static inline __m128i raw_basis_low16(const U *p,bool transformed,long st) {
__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*st)),
_mm_loadu_si128((const __m128i*)(p-32))));
return v;
}
static inline __m256i raw_basis32(const U *p,int flags,long st) {
__m256i v=raw_basis_low32(p,flags&1,st);
if(flags&2) v=_mm256_add_epi16(v,_mm256_sub_epi16(
raw_basis_low32(p-64*st,flags&1,st),raw_basis_low32(p-64,flags&1,st)));
return v;
}
static inline __m128i raw_basis16(const U *p,int flags,long st) {
__m128i v=raw_basis_low16(p,flags&1,st);
if(flags&2) v=_mm_add_epi16(v,_mm_sub_epi16(
raw_basis_low16(p-64*st,flags&1,st),raw_basis_low16(p-64,flags&1,st)));
return v;
}
template<bool TWO> static inline __m256i rb32(const U *p,const U *q,int sgn,int flags,long st){
__m256i v=raw_basis32(p,flags,st);
if constexpr(TWO){ __m256i w=raw_basis32(q,flags,st); v=(sgn>0)?_mm256_add_epi16(v,w):_mm256_sub_epi16(v,w); }
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*4096)))?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,stride), b=raw_basis32(s+(i+1)*stride+k,transformed,stride);
__m256i c=raw_basis32(s+(i+2)*stride+k,transformed,stride), e=raw_basis32(s+(i+3)*stride+k,transformed,stride);
__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*4096)))?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,stride);
__m256i y=raw_basis32(s+(k+1)*stride+j+t,transformed,stride);
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,O0,O1,O2,O3) do { \
__m256i a=rb32<TWO>(sr+(long)i*S4ST+(K),sr2+(long)i*S4ST+(K),sgn,(TF),S4ST), b=rb32<TWO>(sr+(long)(i+1)*S4ST+(K),sr2+(long)(i+1)*S4ST+(K),sgn,(TF),S4ST); \
__m256i c=rb32<TWO>(sr+(long)(i+2)*S4ST+(K),sr2+(long)(i+2)*S4ST+(K),sgn,(TF),S4ST), e=rb32<TWO>(sr+(long)(i+3)*S4ST+(K),sr2+(long)(i+3)*S4ST+(K),sgn,(TF),S4ST); \
st((DSTP)+(O0),a); \
st((DSTP)+(O1),b); \
st((DSTP)+(O2),c); \
st((DSTP)+(O3),e); } while(0)
template<bool TWO> static void pack_wide(U *dA,U *dB,const U *s,long ST4_,const U *s2,int sgn) {
const long S4ST=ST4_;
for (int r=0;r<2;++r) {
const U *sr = s + (long)r*32*S4ST, *sr2 = s2 + (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,0,16,64,80); PWL_CH(p0, 16,0,32,48,96,112); p0+=128;
PWL_CH(p1, 32,TFc,0,16,64,80); PWL_CH(p1, 48,TFc,32,48,96,112); p1+=128;
PWL_CH(p2, 64,0,0,16,64,80); PWL_CH(p2, 80,0,32,48,96,112); p2+=128;
PWL_CH(p3, 96,TFc,0,16,64,80); PWL_CH(p3,112,TFc,32,48,96,112); 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=rb32<TWO>(sr+(long)k*ST+(J),sr2+(long)k*ST+(J),sgn,(TF),ST), y=rb32<TWO>(sr+(long)(k+1)*ST+(J),sr2+(long)(k+1)*ST+(J),sgn,(TF),ST); \
st((P),_mm256_unpacklo_epi16(x,y)); st((P)+16,_mm256_unpackhi_epi16(x,y)); } while(0)
template<bool TWO> static void pack_wide_b(U *dA,U *dB,const U *s,long ST4_,long d0log,const U *s2,int sgn) {
const long ST=ST4_;
const long d0=d0log;
const int b5=(int)((d0>>5)&1), b17=(int)((d0>>17)&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, *sr2=s2+(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_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.
static __attribute__((noinline,noclone)) void native2x32_shared(const U *a,const U *b,U *c,int Mode) {
const U *abase=a+1024;
U *cbase=c+1024;
intptr_t idx=-2048;
asm volatile(
"vpbroadcastd 0(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b0]), %%ymm0\n\t"
"vmovdqa -96(%[b0]), %%ymm1\n\t"
"vmovdqa -64(%[b0]), %%ymm2\n\t"
"3:\n\t"
"vpbroadcastd 32(%[abase],%[i],1), %%ymm4\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm8\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm9\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm10\n\t"
"vpmaddwd -32(%[b0]), %%ymm3, %%ymm11\n\t"
"vpmaddwd %%ymm0, %%ymm4, %%ymm12\n\t"
"vpmaddwd %%ymm1, %%ymm4, %%ymm13\n\t"
"vpmaddwd %%ymm2, %%ymm4, %%ymm14\n\t"
"vpmaddwd -32(%[b0]), %%ymm4, %%ymm15\n\t"
"vpbroadcastd 4(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b0]), %%ymm0\n\t"
"vmovdqa 32(%[b0]), %%ymm1\n\t"
"vmovdqa 64(%[b0]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b0]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 36(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b0]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 8(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b1]), %%ymm0\n\t"
"vmovdqa -96(%[b1]), %%ymm1\n\t"
"vmovdqa -64(%[b1]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd -32(%[b1]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 40(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd -32(%[b1]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 12(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b1]), %%ymm0\n\t"
"vmovdqa 32(%[b1]), %%ymm1\n\t"
"vmovdqa 64(%[b1]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b1]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 44(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b1]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 16(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b2]), %%ymm0\n\t"
"vmovdqa -96(%[b2]), %%ymm1\n\t"
"vmovdqa -64(%[b2]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd -32(%[b2]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 48(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd -32(%[b2]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 20(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b2]), %%ymm0\n\t"
"vmovdqa 32(%[b2]), %%ymm1\n\t"
"vmovdqa 64(%[b2]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b2]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 52(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b2]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 24(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b3]), %%ymm0\n\t"
"vmovdqa -96(%[b3]), %%ymm1\n\t"
"vmovdqa -64(%[b3]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd -32(%[b3]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 56(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd -32(%[b3]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 28(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b3]), %%ymm0\n\t"
"vmovdqa 32(%[b3]), %%ymm1\n\t"
"vmovdqa 64(%[b3]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b3]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 60(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b3]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 64(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b4]), %%ymm0\n\t"
"vmovdqa -96(%[b4]), %%ymm1\n\t"
"vmovdqa -64(%[b4]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd -32(%[b4]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 96(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd -32(%[b4]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 68(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b4]), %%ymm0\n\t"
"vmovdqa 32(%[b4]), %%ymm1\n\t"
"vmovdqa 64(%[b4]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b4]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 100(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b4]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 72(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b5]), %%ymm0\n\t"
"vmovdqa -96(%[b5]), %%ymm1\n\t"
"vmovdqa -64(%[b5]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd -32(%[b5]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 104(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd -32(%[b5]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 76(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b5]), %%ymm0\n\t"
"vmovdqa 32(%[b5]), %%ymm1\n\t"
"vmovdqa 64(%[b5]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b5]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 108(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b5]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 80(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b6]), %%ymm0\n\t"
"vmovdqa -96(%[b6]), %%ymm1\n\t"
"vmovdqa -64(%[b6]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd -32(%[b6]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 112(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd -32(%[b6]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 84(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b6]), %%ymm0\n\t"
"vmovdqa 32(%[b6]), %%ymm1\n\t"
"vmovdqa 64(%[b6]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b6]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 116(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b6]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 88(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa -128(%[b7]), %%ymm0\n\t"
"vmovdqa -96(%[b7]), %%ymm1\n\t"
"vmovdqa -64(%[b7]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd -32(%[b7]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 120(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd -32(%[b7]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vpbroadcastd 92(%[abase],%[i],1), %%ymm3\n\t"
"vmovdqa 0(%[b7]), %%ymm0\n\t"
"vmovdqa 32(%[b7]), %%ymm1\n\t"
"vmovdqa 64(%[b7]), %%ymm2\n\t"
"vpmaddwd %%ymm0, %%ymm3, %%ymm4\n\t"
"vpmaddwd %%ymm1, %%ymm3, %%ymm5\n\t"
"vpmaddwd %%ymm2, %%ymm3, %%ymm6\n\t"
"vpmaddwd 96(%[b7]), %%ymm3, %%ymm3\n\t"
"vpbroadcastd 124(%[abase],%[i],1), %%ymm7\n\t"
"vpmaddwd %%ymm0, %%ymm7, %%ymm0\n\t"
"vpmaddwd %%ymm1, %%ymm7, %%ymm1\n\t"
"vpmaddwd %%ymm2, %%ymm7, %%ymm2\n\t"
"vpmaddwd 96(%[b7]), %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm4, %%ymm8, %%ymm8\n\t"
"vpaddd %%ymm5, %%ymm9, %%ymm9\n\t"
"vpaddd %%ymm6, %%ymm10, %%ymm10\n\t"
"vpaddd %%ymm3, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 128(%[abase],%[i],1), %%ymm3\n\t"
"vpaddd %%ymm0, %%ymm12, %%ymm12\n\t"
"vpaddd %%ymm1, %%ymm13, %%ymm13\n\t"
"vpaddd %%ymm2, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm7, %%ymm15, %%ymm15\n\t"
"vmovdqa -128(%[b0]), %%ymm0\n\t"
"vmovdqa -96(%[b0]), %%ymm1\n\t"
"vmovdqa -64(%[b0]), %%ymm2\n\t"
"test %k[mode], %k[mode]\n\t"
"jz 4f\n\t"
"cmp $1, %k[mode]\n\t"
"jne 8f\n\t"
"vpslld $16, %%ymm9, %%ymm9\n\t"
"vpblendw $0xAA, %%ymm9, %%ymm8, %%ymm8\n\t"
"vpaddw 0(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 0(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm11, %%ymm11\n\t"
"vpblendw $0xAA, %%ymm11, %%ymm10, %%ymm10\n\t"
"vpaddw 32(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 32(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm13, %%ymm13\n\t"
"vpblendw $0xAA, %%ymm13, %%ymm12, %%ymm12\n\t"
"vpaddw 64(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 64(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm15, %%ymm15\n\t"
"vpblendw $0xAA, %%ymm15, %%ymm14, %%ymm14\n\t"
"vpaddw 96(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 96(%[cbase],%[i],1)\n\t"
"add $128, %[i]\n\t"
"jnz 3b\n\t"
"jmp 7f\n\t"
"8:\n\t"
"cmp $2, %k[mode]\n\t"
"je 6f\n\t"
"cmp $3, %k[mode]\n\t"
"jne 9f\n\t"
"vpslld $16, %%ymm9, %%ymm9\n\t"
"vpblendw $0xAA, %%ymm9, %%ymm8, %%ymm8\n\t"
"vpaddw 0(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 0(%[cbase],%[i],1)\n\t"
"vpaddw -8320(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw 8320(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw -16640(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 8320(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm11, %%ymm11\n\t"
"vpblendw $0xAA, %%ymm11, %%ymm10, %%ymm10\n\t"
"vpaddw 32(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 32(%[cbase],%[i],1)\n\t"
"vpaddw -8288(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw 8352(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw -16608(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 8352(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm13, %%ymm13\n\t"
"vpblendw $0xAA, %%ymm13, %%ymm12, %%ymm12\n\t"
"vpaddw 64(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 64(%[cbase],%[i],1)\n\t"
"vpaddw -8256(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw 8384(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw -16576(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 8384(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm15, %%ymm15\n\t"
"vpblendw $0xAA, %%ymm15, %%ymm14, %%ymm14\n\t"
"vpaddw 96(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 96(%[cbase],%[i],1)\n\t"
"vpaddw -8224(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw 8416(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw -16544(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 8416(%[cbase],%[i],1)\n\t"
"add $128, %[i]\n\t"
"jnz 3b\n\t"
"jmp 7f\n\t"
"9:\n\t"
"vpslld $16, %%ymm9, %%ymm9\n\t"
"vpblendw $0xAA, %%ymm9, %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 0(%[cbase],%[i],1)\n\t"
"vpaddw -2080(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw 2080(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw -4160(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 2080(%[cbase],%[i],1)\n\t"
"vpaddw -6240(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw 10400(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw -14560(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 10400(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm11, %%ymm11\n\t"
"vpblendw $0xAA, %%ymm11, %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 32(%[cbase],%[i],1)\n\t"
"vpaddw -2048(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw 2112(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw -4128(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 2112(%[cbase],%[i],1)\n\t"
"vpaddw -6208(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw 10432(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw -14528(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 10432(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm13, %%ymm13\n\t"
"vpblendw $0xAA, %%ymm13, %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 64(%[cbase],%[i],1)\n\t"
"vpaddw -2016(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw 2144(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw -4096(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 2144(%[cbase],%[i],1)\n\t"
"vpaddw -6176(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw 10464(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw -14496(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 10464(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm15, %%ymm15\n\t"
"vpblendw $0xAA, %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 96(%[cbase],%[i],1)\n\t"
"vpaddw -1984(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw 2176(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw -4064(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 2176(%[cbase],%[i],1)\n\t"
"vpaddw -6144(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw 10496(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw -14464(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 10496(%[cbase],%[i],1)\n\t"
"add $128, %[i]\n\t"
"jnz 3b\n\t"
"jmp 7f\n\t"
"6:\n\t"
"vpslld $16, %%ymm9, %%ymm9\n\t"
"vpblendw $0xAA, %%ymm9, %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 0(%[cbase],%[i],1)\n\t"
"vpaddw -2080(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw 2080(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vpsubw -4160(%[cbase],%[i],1), %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 2080(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm11, %%ymm11\n\t"
"vpblendw $0xAA, %%ymm11, %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 32(%[cbase],%[i],1)\n\t"
"vpaddw -2048(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw 2112(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vpsubw -4128(%[cbase],%[i],1), %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 2112(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm13, %%ymm13\n\t"
"vpblendw $0xAA, %%ymm13, %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 64(%[cbase],%[i],1)\n\t"
"vpaddw -2016(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw 2144(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vpsubw -4096(%[cbase],%[i],1), %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 2144(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm15, %%ymm15\n\t"
"vpblendw $0xAA, %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 96(%[cbase],%[i],1)\n\t"
"vpaddw -1984(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw 2176(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vpsubw -4064(%[cbase],%[i],1), %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 2176(%[cbase],%[i],1)\n\t"
"add $128, %[i]\n\t"
"jnz 3b\n\t"
"jmp 7f\n\t"
"4:\n\t"
"vpslld $16, %%ymm9, %%ymm9\n\t"
"vpblendw $0xAA, %%ymm9, %%ymm8, %%ymm8\n\t"
"vmovdqu %%ymm8, 0(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm11, %%ymm11\n\t"
"vpblendw $0xAA, %%ymm11, %%ymm10, %%ymm10\n\t"
"vmovdqu %%ymm10, 32(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm13, %%ymm13\n\t"
"vpblendw $0xAA, %%ymm13, %%ymm12, %%ymm12\n\t"
"vmovdqu %%ymm12, 64(%[cbase],%[i],1)\n\t"
"vpslld $16, %%ymm15, %%ymm15\n\t"
"vpblendw $0xAA, %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovdqu %%ymm14, 96(%[cbase],%[i],1)\n\t"
"add $128, %[i]\n\t"
"jnz 3b\n\t"
"7:\n\t"
: [i] "+&r"(idx)
: [b0] "r"((const U*)((const char*)b + 128)), [b1] "r"((const U*)((const char*)b + 384)), [b2] "r"((const U*)((const char*)b + 640)), [b3] "r"((const U*)((const char*)b + 896)), [b4] "r"((const U*)((const char*)b + 1152)), [b5] "r"((const U*)((const char*)b + 1408)), [b6] "r"((const U*)((const char*)b + 1664)), [b7] "r"((const U*)((const char*)b + 1920)), [abase] "r"(abase), [cbase] "r"(cbase), [mode] "r"(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_shared(a,b,c,0);
}
// 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_shared(a,b,c,1);
}
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;
}
}
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*8*512], tail[8*512];
U *dest[4]={d11,d12,d21,d22};
auto make_all=[&](U *tmp,int group) {
for(int child=0;child<4;child++)
for(int b=0;b<8;b++) {
int block=b*128;
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*8+b)*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 b=0;b<8;b++)
for(int v=0;v<512;v+=16) {
U *t1=tmp+(1*8+b)*512,*t2=tmp+(2*8+b)*512;
__m256i c22=ld(tmp+(3*8+b)*512+v);
st(t1+v,_mm256_sub_epi16(ld(t1+v),c22));
st(t2+v,_mm256_sub_epi16(c22,ld(t2+v)));
}
};
make_all(keep,3);
for(int child=0;child<4;child++) {
int row=64+((child>>1)&1)*32,col=64+(child&1)*32;
for(int b=0;b<8;b++) {
size_t off=(size_t)(b*4)*stride;
for(int r=0;r<4;r++)emit_native(dest[r]+off+row*stride+col,keep+(child*8+b)*512+r*128,stride);
}
}
for(int group=1;group<=2;group++) {
// Only the low-basis child3 must survive for the other children.
for(int b=0;b<8;b++) for(int v=0;v<128;v+=16){
int idx=(4*group+3)*q+b*128+v;
__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=tail+b*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 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 b=0;b<8;b++) for(int v=0;v<128;v+=16){
__m256i z[4];
const U*t=tail+b*512+v;
if(child==3){for(int r=0;r<4;r++)z[r]=ld(t+r*128);}
else {
int idx=(4*group+child)*q+b*128+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));
z[0]=_mm256_add_epi16(v1,ld(p2+idx));
z[1]=_mm256_add_epi16(_mm256_add_epi16(u2,v5),ld(p3+idx));
z[2]=_mm256_sub_epi16(u3,ld(p4+idx));
z[3]=_mm256_add_epi16(u3,v5);
if(child==1)for(int r=0;r<4;r++)z[r]=_mm256_sub_epi16(z[r],ld(t+r*128));
if(child==2)for(int r=0;r<4;r++)z[r]=_mm256_sub_epi16(ld(t+r*128),z[r]);
}
for(int r=0;r<4;r++){
__m256i k=ld(keep+(child*8+b)*512+r*128+v);
__m256i x=group==1?_mm256_sub_epi16(z[r],k):_mm256_sub_epi16(k,z[r]);
U*d=dest[r]+(size_t)(row+b*4+v/32)*stride+col+v%32;
_mm256_stream_si256((__m256i*)d,s4bj_fx(x));
}
}
}
}
} 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*8*512];
U *dest[4]={d11,d12,d21,d22};
for(int child=1;child<=3;child++)
for(int b=0;b<8;b++) {
int block=b*128;
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)*8+b)*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 b=0;b<8;b++)
for(int v=0;v<512;v+=16) {
U *t1=tmp+(0*8+b)*512,*t2=tmp+(1*8+b)*512;
__m256i c22=ld(tmp+(2*8+b)*512+v);
st(t1+v,_mm256_sub_epi16(ld(t1+v),c22));
st(t2+v,_mm256_sub_epi16(c22,ld(t2+v)));
}
for(int b=0;b<8;b++) {
size_t off=(size_t)(b*4)*stride;
for(int r=0;r<4;r++) {
emit_native(dest[r]+off+32,tmp+(0*8+b)*512+r*128,stride);
emit_native(dest[r]+off+32*stride,tmp+(1*8+b)*512+r*128,stride);
emit_native(dest[r]+off+32*stride+32,tmp+(2*8+b)*512+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);
}
}
}
// Share exactly the three repeated fixed-size N128 operand transforms.
template<bool Sub,int Delta,int Len>
static __attribute__((noinline,noclone)) void middle_displaced_shared(U*d,const U*a) {
combine_displaced<Sub,Delta>(d,a,Len);
}
template<int N,bool Sub,int Delta,int Len>
static __attribute__((always_inline)) inline void selected_displaced(U*d,const U*a) {
if constexpr(N==128) middle_displaced_shared<Sub,Delta,Len>(d,a);
else combine_displaced<Sub,Delta>(d,a,Len);
}
template<int N,bool ParentReconstruct=false>
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
selected_displaced<N,false,q,q>(x,a21);
selected_displaced<N,false,q,q>(y,b21);
mul<h>(x,y,c12,next); // m5
selected_displaced<N,true,-2*q,q>(x,a22);
selected_displaced<N,true,-2*q,q>(y,b22);
if constexpr(N==64) native2x32_shared(x,y,c21,ParentReconstruct?4:2); else mul<h,true>(x,y,c21,next); // m6: child finals reconstruct parent
if constexpr(N==64) {
if constexpr(ParentReconstruct) native2x32_shared(a11,b11,c11,3); else leaf_add(a11,b11,c11); // m1 directly accumulates
selected_displaced<N,true,-3*q,q>(y,b22);
if constexpr(ParentReconstruct) native2x32_shared(a21,y,c21,3); else leaf_add(a21,y,c21); // m3 directly accumulates
selected_displaced<N,true,3*q,q>(x,a11); // negate the M7 A input
if constexpr(ParentReconstruct) native2x32_shared(x,b12,c12,3); else leaf_add(x,b12,c12); // c12+=(-oldM7)
} else {
mul<h>(a11,b11,x,next); // m1
combine_inplace<false>(c11,x,q);
selected_displaced<N,true,-3*q,q>(y,b22);
mul<h>(a21,y,x,next); // m3
combine_inplace<false>(c21,x,q);
selected_displaced<N,true,-3*q,q>(x,a22);
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 == 4096) {
combine3_unpack<2048>(final_out+2048,
final_out+(size_t)2048*4096,
final_out+(size_t)2048*4096+2048,
c12,c21,c22,y,x,c11,4096);
} else {
for(int chunk=0;chunk<q;chunk+=1024+TILE_PAD)
_Pragma("GCC unroll 2")
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 == 4096) unpack_add<2048>(final_out,c11,y,4096);
else combine_inplace<false>(c11,y,q);
}
}
static constexpr int RQ=2048*2048+(2048/BASE)*(2048/BASE)*TILE_PAD,RS=RQ+ROOT_PAD;
// r57kb_ [HUGEHINT]: 2 MB-aligned so the kernel can actually back it with huge pages.
alignas(2*1024*1024) static U root_mem[12*RS+4096];
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,4096,Right);
pack(tile+HQ,s+2048,64,4096,Right);
pack(tile+2*HQ,s+2048*4096,64,4096,Right);
pack(tile+3*HQ,s+2048*4096+2048,64,4096,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,bool TWO>
static void pack_wide_root(U *tA,U *tB,const U *s,long STR_,long D0_,const U *s2,int sgn) {
if constexpr(!Right) {
pack_wide<TWO>(tA+0*HQ, tB+0*HQ, s, STR_, s2, sgn);
pack_wide<TWO>(tA+1*HQ, tB+1*HQ, s+2048, STR_, s2+2048, sgn);
pack_wide<TWO>(tA+2*HQ, tB+2*HQ, s+2048L*STR_, STR_, s2+2048L*STR_, sgn);
pack_wide<TWO>(tA+3*HQ, tB+3*HQ, s+2048L*STR_+2048, STR_, s2+2048L*STR_+2048, sgn);
} else {
pack_wide_b<TWO>(tA+0*HQ, tB+0*HQ, s, STR_, D0_, s2, sgn);
pack_wide_b<TWO>(tA+1*HQ, tB+1*HQ, s+2048, STR_, D0_+2048, s2+2048, sgn);
pack_wide_b<TWO>(tA+2*HQ, tB+2*HQ, s+2048L*STR_, STR_, D0_+2048L*4096, s2+2048L*STR_, sgn);
pack_wide_b<TWO>(tA+3*HQ, tB+3*HQ, s+2048L*STR_+2048, STR_, D0_+2048L*4096+2048, s2+2048L*STR_+2048, sgn);
}
}
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,bool TWO=Right>
static void pack_root(const U*s,int n,int offset,U *dest,U *destq,int voff,long STR_=4096,long D0_=0,const U*s2=nullptr,int sgn=0) {
if(n>128){int h=n/2,q=h*h+(h/BASE)*(h/BASE)*TILE_PAD;
pack_root<Right,TWO>(s,h,offset,dest,destq,voff,STR_,D0_,s2,sgn);
pack_root<Right,TWO>(s+h,h,offset+q,dest,destq,voff,STR_,D0_+h,s2+h,sgn);
pack_root<Right,TWO>(s+h*STR_,h,offset+2*q,dest,destq,voff,STR_,D0_+h*4096L,s2+h*STR_,sgn);
pack_root<Right,TWO>(s+h*STR_+h,h,offset+3*q,dest,destq,voff,STR_,D0_+h*4096L+h,s2+h*STR_+h,sgn);
} else {
pack_wide_root<Right,TWO>(root_tile,high_left,s,STR_,D0_,s2,sgn);
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,TWO>(high_right,root_tile,s+64*STR_,STR_,D0_+64*4096L,s2+64*STR_,sgn);
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);
}
}
static void short_product4096s(int n,const short*A,const short*B,short*C,long sA,long sB);
static void short_product4096x(int n,const short*A,const short*A2,int sgA,long sA,const short*B,const short*B2,int sgB,long sB,short*C);
static void short_product4096(int n,const short*A,const short*B,short*C){ short_product4096s(n,A,B,C,4096,4096); }
static void short_product4096s(int n,const short*A,const short*B,short*C,long sA,long sB){
short_product4096x(n,A,nullptr,0,sA,B,nullptr,0,sB,C);
}
static void short_product4096x(int n,const short*A,const short*A2,int sgA,long sA,const short*B,const short*B2,int sgB,long sB,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;
if(A2)pack_root<false,true>((const U*)A,2048,0,ar,cq,0,sA,0,(const U*)A2,sgA);
else pack_root<false,false>((const U*)A,2048,0,ar,cq,0,sA,0,nullptr,0);
raw_basis_input=(const U*)B;
if(B2)pack_root<true,true>((const U*)B,2048,0,br,br,3,sB,0,(const U*)B2,sgB);
else pack_root<true,false>((const U*)B,2048,0,br,br,3,sB,0,nullptr,0);
_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<2048>(s7,t7,p7,work);mul<2048>(s5,t5,p5,work);mul<2048>(s6,t6,p6,work);
mul<2048>(s3,b22,p3,work);mul<2048>(a22,t4,p4,work);mul<2048>(a11,b11,p1,work);mul<2048>(a12,b21,p2,work);
U *out=g_out?g_out:(U*)C;
const long ostr=g_ostr;
root_fuse<2048>(out,out+2048,out+2048*ostr,out+2048*ostr+2048,p1,p2,p3,p4,p5,p6,p7,(int)ostr);_mm_sfence();
asm volatile(".purgem vzeroupper");
_mm256_zeroupper();
}
#include <cstdlib>
#include <sys/mman.h> // r57kb_[HUGEHINT]
static void outer_operand(const unsigned short*a,int sa,const unsigned short*b,int sb,unsigned short*out,int sign){
for(int i=0;i<4096;i++)for(int j=0;j<4096;j+=16){
__m256i x=_mm256_load_si256((const __m256i*)(a+(size_t)i*sa+j));
if(sign){__m256i y=_mm256_load_si256((const __m256i*)(b+(size_t)i*sb+j));x=sign>0?_mm256_add_epi16(x,y):_mm256_sub_epi16(x,y);}
_mm256_store_si256((__m256i*)(out+(size_t)i*4096+j),x);
}
}
static void outer_copy(unsigned short*c,const unsigned short*p){
for(int i=0;i<4096;i++)for(int j=0;j<4096;j+=16){
__m256i y=_mm256_load_si256((const __m256i*)(p+(size_t)i*4096+j));
_mm256_store_si256((__m256i*)(c+(size_t)i*8192+j),y);
}
}
static void outer_copy_s(unsigned short*c,const unsigned short*p,int ps){
for(int i=0;i<4096;i++)for(int j=0;j<4096;j+=16){
__m256i y=_mm256_load_si256((const __m256i*)(p+(size_t)i*ps+j));
_mm256_store_si256((__m256i*)(c+(size_t)i*8192+j),y); }
}
static void outer_accum_s(unsigned short*c,const unsigned short*p,int ps,int sign){
for(int i=0;i<4096;i++)for(int j=0;j<4096;j+=16){
__m256i x=_mm256_load_si256((const __m256i*)(c+(size_t)i*8192+j));
__m256i y=_mm256_load_si256((const __m256i*)(p+(size_t)i*ps+j));
x=sign>0?_mm256_add_epi16(x,y):_mm256_sub_epi16(x,y);
_mm256_store_si256((__m256i*)(c+(size_t)i*8192+j),x); }
}
static void outer_accum(unsigned short*c,const unsigned short*p,int sign){
for(int i=0;i<4096;i++)for(int j=0;j<4096;j+=16){
__m256i x=_mm256_load_si256((const __m256i*)(c+(size_t)i*8192+j));
__m256i y=_mm256_load_si256((const __m256i*)(p+(size_t)i*4096+j));
x=sign>0?_mm256_add_epi16(x,y):_mm256_sub_epi16(x,y);
_mm256_store_si256((__m256i*)(c+(size_t)i*8192+j),x);
}
}
template<bool SUB0>
static void outer_two_destinations(unsigned short*c0,unsigned short*c1,const unsigned short*p){
for(int i=0;i<4096;i++)for(int j=0;j<4096;j+=16){
size_t o=(size_t)i*8192+j;
__m256i v=_mm256_load_si256((const __m256i*)(p+(size_t)i*4096+j));
__m256i a=_mm256_load_si256((const __m256i*)(c0+o));
__m256i b=_mm256_load_si256((const __m256i*)(c1+o));
_mm256_store_si256((__m256i*)(c0+o),SUB0?_mm256_sub_epi16(a,v):_mm256_add_epi16(a,v));
_mm256_store_si256((__m256i*)(c1+o),_mm256_add_epi16(b,v));
}
}
static void outer_quad(const unsigned short*b1,const unsigned short*b2,const unsigned short*b3,const unsigned short*b4,const unsigned short*b5,unsigned short*c00,unsigned short*c10,unsigned short*c11,unsigned short*c01){
for(int i=0;i<4096;i++){
const unsigned short*r1=b1+(size_t)i*4096,*r2=b2+(size_t)i*4096,*r3=b3+(size_t)i*4096,*r4=b4+(size_t)i*4096,*r5=b5+(size_t)i*4096;
unsigned short*w0=c00+(size_t)i*8192,*w1=c10+(size_t)i*8192,*w2=c11+(size_t)i*8192,*w3=c01+(size_t)i*8192;
for(int j=0;j<4096;j+=16){
__m256i x1=_mm256_load_si256((const __m256i*)(r1+j)),x2=_mm256_load_si256((const __m256i*)(r2+j)),x3=_mm256_load_si256((const __m256i*)(r3+j)),x4=_mm256_load_si256((const __m256i*)(r4+j)),x5=_mm256_load_si256((const __m256i*)(r5+j));
_mm256_store_si256((__m256i*)(w0+j),_mm256_sub_epi16(_mm256_add_epi16(x1,x4),x5));
_mm256_store_si256((__m256i*)(w1+j),_mm256_add_epi16(x2,x4));
_mm256_store_si256((__m256i*)(w2+j),_mm256_add_epi16(_mm256_sub_epi16(x1,x2),x3));
_mm256_store_si256((__m256i*)(w3+j),_mm256_add_epi16(x3,x5));
}
}
}
void matrix_multiply(int n,const short*A,const short*B,short*C){
if(n==4096){short_product4096(n,A,B,C);return;}
if(n!=8192){for(int i=0;i<n;i++)for(int j=0;j<n;j++){unsigned z=0;for(int k=0;k<n;k++)z+=(unsigned short)A[(size_t)i*n+k]*(unsigned)(unsigned short)B[(size_t)k*n+j];C[(size_t)i*n+j]=(short)z;}return;}
using V=unsigned short;V *pb;
constexpr size_t z=(size_t)4096*4096;
// ---- r57kb_fz_ [FUSE-2] -------------------------------------------------------------
// The driver used to run 31 whole-quadrant sweeps (~992 MB) after the engine output: 7 product
// writes plus 8 outer_* passes, with c11 alone paying 4 sweeps (352 MB) and c00 three (288 MB).
// Now every product lands in its own 32 MB buffer first and each quadrant is combined from
// registers in one sweep, so the same eight combinations cost 22 sweeps (~704 MB): -288 MB at
// the family's measured "driver whole-quadrant transfer" price of ~50 ms/GB (lower bound).
// Five buffers = 5 x 32 MB = 160 MB of freshly-mapped memory, advised huge like C/root_mem
// below (advisory only; a THP=never host ignores it and no output bit changes).
madvise((void*)C,(size_t)8192*8192*2,MADV_HUGEPAGE);
if(posix_memalign((void**)&pb,4096,5*z*2))abort();
madvise((void*)pb,5*z*2,MADV_HUGEPAGE);
madvise((void*)root_mem,sizeof(root_mem),MADV_HUGEPAGE);
V*b1=pb,*b2=pb+z,*b3=pb+2*z,*b4=pb+3*z,*b5=pb+4*z;
const V*a00=(const V*)A,*a01=a00+4096,*a10=a00+4096*8192,*a11=a10+4096;
const V*b00=(const V*)B,*b01=b00+4096,*b10=b00+4096*8192,*b11=b10+4096;
V*c00=(V*)C,*c01=c00+4096,*c10=c00+4096*8192,*c11=c10+4096;
auto P7=[&](const V*AA,const V*AA2,int gA,const V*BB,const V*BB2,int gB,V*O){short_product4096x(4096,(const short*)AA,(const short*)AA2,gA,8192,(const short*)BB,(const short*)BB2,gB,8192,(short*)O);};
P7(a00,a11,1,b00,b11,1,b1); // p1 = (a00+a11)(b00+b11)
P7(a10,a11,1,b00,nullptr,0,b2); // p2 = (a10+a11)b00
P7(a00,nullptr,0,b01,b11,-1,b3); // p3 = a00(b01+b11)
P7(a11,nullptr,0,b10,b00,-1,b4); // p4 = a11(b10-b00)
P7(a00,a01,1,b11,nullptr,0,b5); // p5 = (a00+a01)b11
outer_quad(b1,b2,b3,b4,b5,c00,c10,c11,c01);
P7(a10,a00,-1,b00,b01,1,b1); // p6 = (a10-a00)(b00+b01)
outer_accum(c11,b1,1);
P7(a01,a11,-1,b10,b11,1,b2); // p7 = (a01-a11)(b10+b11)
outer_accum(c00,b2,1);
free(pb);
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 3.691 s | 385 MB + 716 KB | Accepted | Score: 100 | 显示更多 |