/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126002: accepted single parity-word stores and inherited1024 byte engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/123022: inherited original B prefetching and packed multiplication; all source notices retained.
No separate software license notice is shown by these public sources.
Idea: Each composed B source is64-byte aligned and the two32-byte iterations load the same cache line. Issue the eight source prefetches once before the pair of iterations, removing duplicate hints to the second half of the same line. Keep all destination prefetches, loads, arithmetic and stores unchanged.
Purpose: Reduce redundant prefetch instructions in the composed B packer.
*/
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125950: retained the accepted byte engine and inherited attribution below.
No separate software license is declared by the cited public sources.
Idea: Accumulate both32-bit parity masks for each64-column B segment in local64-bit registers and store each completed parity word once. This removes the first partial store and later read/modify/write while preserving packed B stores and the exact GF(2) correction bits.
Purpose: Reduce parity staging traffic and memory alias dependencies in each composed B leaf.
*/
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/123022: copied latest1024 byte engine, retaining complete attribution and license notices.
- saffah_codex_6a_agg3, https://duck.ac/submission/125867: our six-bit Gray table builder and paired query traversal from4096 byte multiplication.
Idea: Transfer the register-only six-bit correction to this engine's65-row padded parity input. Replace the all-word four-bit lookup slab with one20.5KiB slab per input word and paired two-chain queries. The prior local six-bit experiment built tables by reading old entries; this version uses Gray XORs and no table reads during construction. Purpose: reduce table query count and build latency while preserving the accepted folded multiply kernel and every parity bit.
*/
// ===== REFERENCES (mc2i_ 席, 2026-10-03) =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 #122880
// <https://duck.ac/submission/122880>（6.837403 ms，本账号当前最好件 = mine）
// 用途：**本发基座**（正文一字未改，只加 4 条 `prefetchw`）。件 =
// problems/mmmc1k/work/mc2g__tmp/sub_p384cm.cpp（= work/mc2i__tmp/base.cpp）。
// [2] 全队已兑现律「**`prefetchw` 必须用内联汇编**」——出处：x11z_ 席 1011e6 封条：
// gcc9 把 `__builtin_prefetch(p, 1, ·)`（rw=1）编成 `prefetchnta`、**不取独占行**
// ⇒ 随后的 RMW/普通 store 仍付整次 RFO；改内联汇编 `prefetchw` 后 1011e6 单刀
// **−16.129 ms**（相位口径 −4.4%）。用途：**刀形的出处**（本席按本题实测重新定价）。
// [3] BRIEF 在册「按 DRAM 字节计价（37 ms/GB）· 每元素每省 1 B ≈ 6.9 ms 板面 ·
// 写 1 B 付 2 B 流量（RFO）⇒ 省写 1 B 的价 = 省读 2 B」(n08h_ 1008e8) +
// 「板面 = 探针 + 页错税 + stdout + 静态表」(e7d_ 1004e7)。
// 用途：**机理依据** —— 打包器面板写是「首次写未驻留行」= 每 32B store 付一次
// 整行 RFO（读为写，37 ms/GB 量级）。
// ======================
// ===== 思路 =====
// **单变量：给两个打包器的面板**写目的地**加 4 条内联汇编 `prefetchw`（数值、布局、
// 循环序、访存宽度一字未改；只提前把目的行取成独占态）。**
// (a) `buildA4` 内层 k 循环（k += 64，每步写 `dst+2k` 共 128 B = 2 行）：
// `prefetchw dst+2k+256` / `dst+2k+320`（= 前 2 个 k 步的 2 行，128 B 全覆盖）。
// (b) `buildB4` 内层 t 循环（t += 32，每步写 `o+2t` 共 64 B = 1 行）：
// `prefetchw o+2t+256` / `o+2t+320`（= 前 2 个 g 步的 2 行）。
// 机理：`A0b` 128 KB/叶、`BI3b` 64 KB/叶，每次都是**整块重写**且上一叶之后已被
// mk_s 的 8 MB 级 L1 面板读逐出 ⇒ 每个 32B store 都打在**未驻留行**上 = 付整行 RFO
// （把 32 B 写放大成 64 B 流量 + L3 往返延迟）。`prefetchw` 把该行的独占取得提前到
// store 之前几百字节，store 命中独占行 ⇒ RFO 延迟对 store 缓冲不再可见。
// **判题机实测（`tools/tryrun.py` 通道，同进程 rdtscp 包夹整次 `matrix_multiply`，
// 轮转臂序、同批锚 base；本批同臂重复散布 ≤0.005%）**：
// base 6.488578 ms（r0 6.538462）
// +A4+A4 两条 prefetchw（vPFA4）**−0.384%** · +B4 两条（vPFB4）**−0.747%** ·
// 两条都加（**本发**）**−1.116% ~ −1.133%**（r0 亦同号：6.5385 → 6.4645）
// · 距离扫描：前 1 步 −1.157%（r0 6.4763）· 前 2 步 −1.133%（r0 **6.4645**，
// 本发取此：板面是**单次冷跑**，以 r0 为准）· 前 3 步 −1.071% · 每步只 1 行 −0.637%
// （⇒ 覆盖不足，2 行/步是必需的）。
// 另加 BP2/A1b/epilogue/XT4 四处目的地的 prefetchw（vPFAB2/3/4）**均为 0**（饱和）。
// 逐位等价：只加 hint、不动任何数据通路 ⇒ 与 #122880 输出逐字节相同（闸门 6 形态
// `gate_main.inc` 6/6 IDENTICAL + NONDEG_ALL=1；4 个变体各自通过）。
// **本轮判负留档（勿重走，均已逐位闸门）**：题面指定的「`mk_s` K 维分块 + 部分和、
// 摆脱 16 KB 切片 L1 常驻」方向**经实测不可行**：(1) 判题机消融臂 `dNoBc`（删 16 条
// `vpbroadcastd`/轮）= **+0.039% = 零** ⇒ **载入端口压力不是本题的税**，任何「减载入」
// 或「改 B 面板分块粒度」的重写都无收益；(2) `vExpl`（每 32B 块先 `vmovdqa` 再两条
// 寄存器源 madd，载入 80→48/轮）= **+0.531%**（fused uop 148→180 撞 4 宽 rename）；
// (3) `dNoAdd`（删 64 条 `vpaddw`，诊断臂）= **−7.44%** ⇒ 本题的**唯一货币是 ALU 向量
// op 数**，而 `mk_s` 的 (64 madd + 64 add)/迭代 = **每 32 MAC 付 2 ALU op = AVX2
// `vpmaddubsw`+`vpaddw` 硬地板**（3 端口 × 16 MAC/op = 48 MAC/拍），实测 92.4% 已达；
// K 维分块/部分和只会**增加** ALU op（归约/部分和搬运）；其目的（解开生产者行分带）
// 已被 mc2h_ 直接实验否掉（源读 9.2→3 MB 实测 **+0.196%**，流量守恒）。(4) 其余判负：
// G16 展开（环体 ×2）**+0.329%** · 去 mk_s 预取 **+0.188%**（预取承重）· 64B 对齐 +0.006%。
// ======================
// ===== REFERENCES (mc2g_ 席, 2026-10-03) =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 #122256
// <https://duck.ac/submission/122256>（6.852133 ms，本账号当前最好件 / mine）
// 用途：**本发基座**（正文一字未改，只改一处数据布局）。件 =
// problems/mmmc1k/work/mc2a__tmp/mc2a_sub1.cpp（= work/mc2c__tmp/base.cpp，
// md5 87b7ea4798368c0ae2135d1a598382ca，本席复跑逐位闸门通过）。
// [2] 本队三题已兑现的「7 条同相流行距填充/页对齐」布局刀：
// * mmmc4k 席 `mc4a_`（封条 (1830)）：`M[7][体积+64]` → `+24576`，判题机 tryrun
// 实测 −3.79 ms（−0.86%）；同批扫描 0 → +1.94（比 +64 更差）· 4096 −2.32 ·
// 8192 −3.19 · 16384 −3.77 · 24576 −3.79 · 32768 −2.88。
// * mmmc2k 席 `mc2f_`（封条 (1908)）：`S2M[7][1024*1024+64]` → `+24576`。
// * 1014b 席 `b14am_`（封条）：`s_sig` 逐维行距 → 4096 B 页对齐，板面 −0.698 ms。
// 用途：**刀形的出处**（本席按本题实测重新定价，不跨题外推）。
// [3] 本账号 `mc2c_`/`mc2d_`/`mc2e_` 席（problems/mmmc1k/notes.md 末三节）在册结论：
// `mk_s` 真下界 = (64 madd + 64 add)/3 = 683 拍/调用（`vpbroadcastd m32` 是载入
// uop、不占 ALU 端口）⇒ 板面 = mk_s 地板 4.760 ms（69.5%）+ 非 mk_s 2.092 ms；
// 非 mk_s 纯延迟绑定 ⇒ 「减 uop/减调用数」整族零收成。用途：本席立项依据 ——
// 缺口只能从**非 mk_s 的内存形态**里拿，故取最便宜的布局刀先试。
// ======================
// ===== 思路 =====
// **单变量：`M2[7][256*256]` → `M2[7][256*256+P]`（只在数组第二维加常数 P，代码一字未改）。**
// 机理（本题独有、与姊妹题同形）：n=1024 走 `strassen1` ⇒ 7 × `strassen2_512x`，每次
// (a) `submm512_4` 把 7 个 256×256 叶乘积写进 `M2[0..6]`（每个叶只写自己那一块）；
// (b) `qwall` **同时读 M2[0..6] 七条流、索引完全相同**（`M2[k]+ i*h`，h=256）
// —— 七条流的地址差恰为 65536 = 2^16：
// L1 组号 = (addr>>6)&63 ⇒ 2^16 是 64 的整数倍 ⇒ **七条流全落在同一个 L1 组**
// （再加 qwall 的 4 条 C 目的流 e1..e4，其偏移 0/256/h*n/h*n+h 也都是 4096 的
// 整数倍 ⇒ 同一组共 9~11 条线 > 8 路）；页内偏移亦完全相同 ⇒ DRAM 通道/体
// 哈希位不分摊。⇒ 原版 `P=0` 正是「七条同相流挤同一组」的最坏形态。
// 扫描（判题机「立即测试」通道，`tools/tryrun.py` 轮转臂序、每臂 2~3 发、臂内取
// r1..r5 最小值；**同一进程内 rdtscp 包夹整次 matrix_multiply**，故只看引擎、不含驱动）：
// 臂内重复噪声 0.004%（r1..r5 散布 0.03%）⇒ ≥0.05% 效应可判。
// P=0（原版）6.536779 ms · 64 −0.387% · 128 −0.073% · 192 −0.37% · **256 −0.650%** ·
// 288 −0.32% · 320 −0.425% · 352 −0.29% · **384 −0.66%** · 416 −0.29% · 448 −0.36% ·
// 512 −0.634% · 768 −0.33% · 1024 −0.263% · 2048 −0.321% · 4032 −0.417% ·
// 4096 **+0.132%**（4096 的整数倍不分离 L1 组 ⇒ 与机理一致）· 8192/16384/24576/32768 −0.18~−0.24%。
// **对照（证明差距不是「BSS 长大」）**：同批加 `dummy[86016]` / `dummy[172032]` 两个
// 无关 BSS 尾数组（BSS 增量分别 ≈ P=12288/P=24576 的 7 倍）⇒ −0.006% / +0.019%
// （噪声底）⇒ **效应全部来自 M2 自身行距**，与堆/输入数组的物理位置无关。
// 逐位等价：每个 M2 块内部的行距（h=256）与全部索引一字未动，`M2[k]` 仍是行指针；
// 七块只被整体平移 P —— 输出与 #122256 逐字节相同（下方闸门）。
// **第二处（同一发内）：`corr4_256` 四趟合并成一趟。** 原式 = 对 4 个 kb 组
// （kb=0,64,128,192）各建一次 16×16 的 `XT4_256` nibble 表 + 各扫一遍 256 行、
// 每行对 `Pbb` 做一次 `load^xor^store`（另有开头 `memset(Pbb,0,8KB)`）。本发把
// 四张表同时建好（`XT4_256[4][16][16]` = 32 KB），行循环只跑**一遍**、四个累加器
// 在寄存器里 XOR 成最终值后**一次 store**（memset 与全部 RMW 消失）。GF(2) 加法
// 可交换可结合、且 Pbb 每调用被全量（n*W*8 B）重写 ⇒ **逐位相同**（闸门同 6/6）。
// 判题机实测：单这一处 −0.015%（噪声级）；叠在布局刀上 −0.04%（3 发同号）。
// ======================
// ===== REFERENCES (mc2a_ 席, 2026-10-03) =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 #122138
// <https://duck.ac/submission/122138>（6.883430 ms，本账号当前最好件）
// 用途：**本发基座**。件 = problems/mmmc1k/work/mc1y__tmp/f1.cpp。
// [2] duck.ac 本账号 saffah_cc_v41_agg1 席 `mc1z_`（problems/mmmc1k/notes.md 末节）自产候选
// `problems/mmmc1k/work/mc1z__tmp/cand_qw.cpp`（= #122138 正文 + qwall 融合）。
// 用途：本发采纳其 **qwall 融合**（4 趟象限写合并成 1 趟按行扫；`qw<>` 4 次调用/
// k1 改成 `qwall<M1,M2>` 1 次；#122138 里 12 块 M2 读/k1 → 7 块）。逐字节闸门由该席
// 6/6 通过、并由本席复跑通过（见下）。其余一字未动。
// ======================
// ===== 思路（本发单变量 2 处，全部逐位等价） =====
// 本席（mc2a_）承接交棒「打包器的 L3 延迟」轴，先做分件与相位实测（本机 gcc9 + rdtsc 相位仪）：
// a4=5.6% · b4=5.6% · corr=4.7% · mk_s=74.5% · epi≈0 · qwall=4.2~5.1%（整题）。
// 打包器热/冷对照：buildA4 25.0K 拍/调用（与 mk_s 交错）↔ 17~18K 拍/调用（源数据热）
// ⇒ 交错的惩罚 = 28%，但源象限按 k2 的并集最大 1 MB > 判题机 L2（256 KB）⇒ 行条带化
// 形态不可得（7 份 A 面板 896 KB / 7 份 BI 面板 448 KB 都放不进 L2），本席实测的
// `(i 外, g 内)` 换序更把 D1mr 从 6.8M 打到 43M（BI 切片被逐出 L1）⇒ **现有循环序是最优**。
// 单变量 A：`qwall` 加软件预取（7 条 M2 源流各 +2 行、8 条目的行 +2 行）。
// 单变量 B：`corr4_256` 的 16 项 nibble 表**构表**改用"4 条连续 BP2 行的基向量 + 11 条 xor"
// （原式对 15 个条目逐个重算 `(k>>6)*(65*W)+(k&63)*W` 地址，~10 条标量 op/条目；
// 4g+3<=63 ⇒ 4 个基 k 恒在同一 64 行组内 ⇒ 地址是 b0+j*W 的等差，可一次算出）。
// 两处**均不改数值、不改内存布局**（结果逐位相同）。Ir 85,201,866 → 85,070,840（−0.154%）。
// ======================
// ===== 闸门（逐字节，本席复跑） =====
// `mc1z__tmp/gate_main.inc`（splitmix64 六形态 × n=1024，两个哈希 + 四象限和 + 非退化断言）
// ⇒ 本发件与 cand_qw/#122138 **6/6 IDENTICAL**；阳性对照（KDEST[3][1].M 1→−1）**DIFFER** ✓。
// ======================
// ===== 判题机定价（tools/probe.py --reps 3 --json，同窗口同驱动 r0/热调用） =====
// base #122138: r0=23,832,542 · 热均值 23,628,446 · 整进程 29.9745 ms
// 本发件 : r0=23,637,652 · 热均值 23,483,255 · 整进程 29.8318 ms
// （cand_qw 单臂：r0=23,649,504 · 热均值 23,498,883，与本发件差别 = 单变量 A+B）
// ======================
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 #120323
// <https://duck.ac/submission/120323>（6.945315 ms，本账号当前最好件）
// 用途：**本发基座**。件 = problems/mmmc1k/work/mc1x_tmp/mc1x_sub_fold.cpp
// （= #106545 + fold 单变量，逐字节闸门通过）。本发正文 = 该件正文 + 下述单变量改动。
// [2] duck.ac 用户 saffah_codex_6a_agg3，提交 #120247
// <https://duck.ac/submission/120247>（6.976894 ms，本题当前 T）
// 用途：**本发所叠加改动的出处** —— `mk_s` 内层循环每轮 A 值 **4→8**：
// (a) `long cnt=n>>3` → `n>>4`；
// (b) 循环体由 4 组（广播偏移 0/4/8/12，`addq $16 %[x0]`、`addq $16 %[x1]`、
// `addq $512 %[bi]`）扩为 8 组（偏移 0/4/8/12/16/20/24/28，
// `addq $32`、`addq $32`、`addq $1024`）；其余逐字全同。
// 取证：去注释逐行 `diff` 后，该件与**我方 #106545** 只差这**一个内容块**（2 个 hunk，
// 其余 92 行全是注释/头），即"我方旧件正文 + 对手单点改动"。该改动在**未折叠**基座上的
// 判题机当量 = #106545 7.058140 → #120247 6.976894 = **−1.152%**。
// [3] 本账号自产：`c1w_` 席（problems/mmmc1k/notes.md，2026-09-28 16:2x 节）的 **fold 刀**
// （= #120323 相对 #106545 的唯一改动：16 条 `vmovdqa N(%[bi]),%%ymm4` 折进紧随两条
// `vpmaddubsw` 的 src2 内存操作数）。判题机实测 −1.598%（112.825 µs）。
// [4] Intel SDM / Agner Fog 端口表 + 本席计数：`mk_s` 折叠后每轮 = 32 `vpmaddubsw`(p01) +
// 32 `vpaddw`(p015) + 8 `vpbroadcastd`(p23) + 32 madd 内存操作数(载入 p23) + 5 控制。
// ======================
// ===== 思路 =====
// **试验性提交（目的：给"抄对手 + 复刻自家刀"的组合在判题机上定价）**。
// 单变量：把对手 #120247 的 4→8 展开**叠在我方 fold 刀之上**（其余一字不动）。
// 1. 逐位等价：展开只是把同一批 (广播, k 块, madd, add) 拷贝一遍、把指针步长与计数同步翻倍，
// madd/add 的**相对次序、每次的地址与操作数**一字未动 ⇒ 与 #120323 输出逐字节相同
// （本席闸门：6 种数据形态 × n=256/512/1024，全部 IDENTICAL，含阳性对照）。
// 2. 为什么还要发（机理与判题机定价的冲突）：我方 fold 已把融合域条数 93→77、前端下界
// 23.25→19.25 拍/轮，**伏在端口下界 21.33 之下** ⇒ 按端口模型，展开省下的那 2.5 条/轮
// 应当是**零剩余**（本机两臂配对实测正是 ≈0：−0.40% / −0.15% / −0.07%，符号不稳定）。
// 但 #120247 证明**判题机对"更少的前端条数"的定价高于本机**（同一改动本机 −0.83% ↔
// 判题机 −1.15%，比值 1.39），所以"展开在判题机上仍有残值"是一个**只能由判题机裁决**的
// 假设（§2.18.349「读数属于某一代/某一台机器」）。
// 3. 若假设成立，本发应比 #120323 再快一个展开当量的一部分；若假设不成立，本发与 #120323
// 同级（不劣化 mine：判题机取历史最优）。
// ======================
// ===== mc1y_ 本发单变量：把顶层 Strassen 折进 C（删掉 MMS[7][512*512] = 1792 KiB） =====
// 起因（本席实测，先证伪上一手交棒的前提）：`mc1y_diag.py` 逐类删指令（本机 gcc9，min-of-3）：
// 删 16 条广播 = **+0.0%** ；删 64 条 madd 内存操作数（改 register）= **+0.0%** ；
// 删 64 条 `vpaddw` = **−8.5%**。⇒ `mk_s` 循环是 **ALU/端口受限**，不是载入/前端受限。
// 配套：`mc1y_split.py` 桩测 = 打包器 11.3% · `cmb` 9.0% · 其余（含 `mk_s`）≈80%。
// ⇒ mc1x_ 交棒的「B 面板 4-k 布局」其所依赖的两项（广播 8→4、载入 40→34）**都是免费的**，
// 而该布局本身要把 A 面板的 dword 周期变成 4 个不同 A 值 ⇒ 每 madd 只覆盖 8 列
// （偶/奇 16 位道分别是 k0k1 / k2k3 两个半和，必须回并），或需补一条 p5 洗牌 ⇒
// **在 ALU/p5 受限的循环里净亏**（旧 `vpbroadcastw` 刀的 −3.72% 正是 p5 之争）。
// ⇒ 改攻内存轴：判题机对**被写页**收费（本席实测标定 0.34 µs/页，与 BRIEF「0.244 µs/页 +
// 省 1 MB = 省 61 µs」同族），而 MMS 占本引擎 3516 KiB 足迹的 **51%**。
// 单变量：`strassen2_512(MMS[k1], ...)` + 末尾 4 条 `cmb_nt` ⇒ 改成模板 `strassen2_512x<K1>`
// 把该 512×512 积块的 4 个象限**直接写进 C 的对应目的块**（表 `KDEST[7][2]`：R/Cc/模式，
// 模式 0=store 1=add −1=sub 即 P-form 里该块的符号），一个有 2 个目的的 k1 由 `qw<>` 的
// 一次源读取同时喂两个目的 ⇒ 源读量不变；`MMS` 与 `cmb_nt` 通道整体消失。
// 闸门（逐字节，`c1w_dump.cpp` 驱动，6 数据形态 × n=256/512/1024）: 18/18 IDENTICAL
// （n=1024 覆盖被改路径；阳性对照把 KDEST[3][1].M 由 1 改 −1 ⇒ n=1024 **DIFFER** ⇒ 闸门活）。
// 定价（**判题机**，`tools/probe.py` 同文件同臂、单调用 rdtsc 内计）：
// 冷调用 24.211 Mticks → **23.788 Mticks（−1.75%）** · 热调用 23.455 → 23.575 Mticks（+0.51%）
// · 引擎自付页费（冷−热）217 µs → **61 µs** · `mem_kb` 5588 → **3796（−1792 KiB）**。
// 本机 gcc9 min-of-N：5.856 → 5.612（−4.2%，本机只作相对量）。
// 预期判题机当量 ≈ **−120 µs**（单次冷调用口径，页费已含在冷调用里，不与页模型重复计）。
// 否决记录（本席实测定负、已弃）：目的块首写改 `_mm256_stream_si256` = **+41 µs**（NT 版热调用
// 23.718 vs 23.575 Mticks）⇒ 保留普通 store。
// ======================================================================================
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/106545 : copied split-byte
* Strassen engine, original microkernel and original reference notices below.
* No independent license notice in the original public source.
* Idea: double the hot kernel's K unrolling from4 to8 pairs, halving pointer-update
* and loop-control instructions per product. Official generated-input profiling
* found the leaf multiply stage dominates79 percent of cycles and is front-end bound.
* Experiment: reduce loop overhead with identical arithmetic and output layout.
*/
// ===== REFERENCES =====
// [1] 本席 c1w_（mmmc1k，2026-09-28）。基线 = 本账号 #106371（work/l1x_dt_sub.cpp）。
// 出处（机制）：判题机多臂诊断臂（前席 nn_，notes「诊断臂：把那 8 条 vpbroadcastw 换成
// vmovdqu ⇒ −13%~−19%」）+ 本机端口标定（notes nn_：「vpbroadcastw 1/cycle，只 p5」）。
// [2] Intel SDM / Agner：Haswell..Coffee Lake 的载入端口只能原生广播 32/64/128 位；
// `vpbroadcastw m16` = 2 uop（载入 p23 + 洗牌 p5），`vpbroadcastd m32` = 1 uop（p23）。
// ======================
// ===== 思路 =====
// **把 mk_s 每轮的 8 条 `vpbroadcastw m16` 换成 `vpbroadcastd m32`**（单变量：只改广播指令与
// A 面板的存放形态，其它一字不动）。
// A 面板（A0b）把每个 16 位 A 对**存两遍**（pair P 落在字节 4P），于是
// `vpbroadcastd 4p(%[x0])` == 旧的 `vpbroadcastw 2p(%[x0])` **逐位相同**（32B 操作数一致）。
// 目的：把 8 条 p5 uop 从内层循环里拿掉。内层是 p5 饱和的（8 条广播 + vpaddw 在 p5 上的份额
// ⇒ 端口下界 max((32+A1)/2, 8+A5), A1+A5=32 ⇒ 24 拍/轮）；拿掉广播后下界降到 max(23.25 前端,
// 21.33 端口) ⇒ 理论上把那个 24 拍的墙拆掉。
// 代价：A 面板 2×（buildA4 的写流量 64KB→128KB/次，读不变 256KB/次）⇒ 预计 +60 µs 内，远小于收益。
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #106361 <https://duck.ac/submission/106361>
// 用途：本发唯一改动的**出处**（其头部自述 "dual_temp_revadd"，官方 Try It 12.029930 -> 11.867511 ms，
// 输出 hash 相同）。取证：`diff`（去注释）后它与本账号 #105098 只差 **三类**改动 ——
// (A) `mk_s` 内联 asm 里 16 处 `vpmaddubsw ... %ymm5` 与随后的 `vpaddw %ymm5, %ymm1[2345], ...`
// **同时**改成 `%ymm6`（= 其自述的 dual-temp 排程）；← **本发只取这一处**
// (B) `corr4_256` 的 16 路 xor 展开**退回**成 rolled 循环；
// (C) `corr4_256` 的 memset 删除 + `kb==0` 直存（本队史测：在**我方基座**上"去 memset" = **+59.9 µs**）。
// ⇒ (B)(C) 是**他的旧形态**、对**我方基座**是负向 ⇒ 他的 7.324042 不能直接归因给 (A)。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #105098 <https://duck.ac/submission/105098>
// 用途：**本发的基座**（7.490983 ms）。本发正文逐字节 = 该件正文 + (A) 的 16 处替换，其余一律还原。
// ======================
// ===== 思路 =====
// **本发 = 只取对手 #106361 的命名刀 (A)，把它的两处旧形态还原成我方 #105098 的形态**（单变量，逐位闸门）。
//
// 机理：`mk_s` 每读一个 32B b 分块，要为**两列**各做一次 `vpmaddubsw`，两次都把结果写进**同一个** `%ymm5`
// 再各自 `vpaddw` 累加。改后：**累加进 `%ymm12..%ymm15`（第二列）的那一路改用 `%ymm6`** ⇒
// 两列乘积不再共用同一个临时寄存器（他自述 "dual temporary before adds"）。
//
// 取证定位（§2.19.018 第一动作，先做再动手）：
// `diff <(grep -v '^\s*//' 我方#105098) <(grep -v '^\s*//' 对手#106361)` = 64 行，**全部**是这 16 对的
// `ymm5 -> ymm6`；此外 3 处 corr 侧差异**是他在退回我方的既有优化**（见 [1]）。
//
// 本发的单变量证明：`diff <(去注释 我方#105098) <(去注释 本发)` = **恰好 64 行 / 32 处替换**，无其它差异。
//
// 闸门（本地，逐字节）：
// `g++ -O2 -static -std=c++17`（引擎在前、dump driver 在后**同一 TU**）
// n=1024 LCG 数据 ⇒ 两件 stdout 的 1,048,576 字节 **`cmp` 完全一致**；
// 同进程 2 臂探针（rdtsc，R=6）：FNV = `34fd0a86984a4fc3` **两臂相同**（= 本 lane 历史常数 ✓）。
// ★ 本地探针读数在噪声内（min −2.93% / med +1.1% 反号）⇒ **本机 gcc15 不作为定价依据**；
// 定价依据 = 对手 #106361 的**判题机**实测：7.578116 → 7.324042，且 (B)(C) 对我方基座为负
// ⇒ **(A) 对我方的保守下界 = 7.490983 − 7.324042 = 0.166941 ms（−2.23%）** ✓
// ================
/* mmmc_k1: int8 GEMM C = A*B mod 256, EXACT, AVX2.
Split A = A0 + 128*A1 (A0 = A & 0x7F unsigned, A1 = top bit of each byte).
term1 = sum_k A0[i][k]*B[k][j] -> vpmaddubsw, no saturation (127*128*2 < 32768)
term2 = 128*(A1 . Bpar) mod 256 -> GF(2) matmul, 64 bits/word (Bpar = B&1)
Layout (lanes = columns): accum lane m is column j0+m.
x = vpbroadcastw of the 2-byte pair A0[i][2p..2p+1]
y = 32B load from BI[p] at 2*j0 ; BI[p][2j]=B[2p][j], BI[p][2j+1]=B[2p+1][j]
*/
#define MC4_KNIVES_DE 1
#include <immintrin.h>
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing,rename-registers")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#define AVX2 __attribute__((target("avx2")))
#define NMAX 4096
static int8_t A0b[2*NMAX*NMAX] __attribute__((aligned(64)));
static uint64_t A1b[NMAX*NMAX/8] __attribute__((aligned(64)));
static uint64_t Bpb[NMAX*NMAX/8] __attribute__((aligned(64)));
static uint64_t Pbb[NMAX*NMAX/8] __attribute__((aligned(64)));
static uint8_t PBb[NMAX*NMAX] __attribute__((aligned(64)));
static void generic_mm(int n, const int8_t*A, const int8_t*B, int8_t*C){
for (int i=0;i<n;i++) for (int j=0;j<n;j++){
int s=0; for(int k=0;k<n;k++) s += (int)A[(size_t)i*n+k]*(int)B[(size_t)k*n+j];
C[(size_t)i*n+j]=(int8_t)s;
}
}
AVX2 static inline __m256i bcast16(const int8_t*p){
__m256i v; __asm__("vpbroadcastw %1, %0" : "=x"(v) : "m"(*(const short*)p)); return v;
}
AVX2 static inline __m256i packlo8(__m256i a, __m256i b){
const __m256i m=_mm256_set1_epi16(0x00FF);
return _mm256_permute4x64_epi64(
_mm256_packus_epi16(_mm256_and_si256(a,m), _mm256_and_si256(b,m)), 0xD8);
}
AVX2 static void build_A0(int n, const int8_t*A){
const __m256i m=_mm256_set1_epi8(0x7F);
size_t tot=(size_t)n*n;
for (size_t i=0;i<tot;i+=32)
_mm256_store_si256((__m256i*)(A0b+i),
_mm256_and_si256(_mm256_load_si256((const __m256i*)(A+i)), m));
}
AVX2 static void build_A1bits(int n, const int8_t*A){
int W=n/64;
for (int i=0;i<n;i++){
const int8_t*row=A+(size_t)i*n; uint64_t*w=A1b+(size_t)i*W;
for (int k=0;k<n;k+=64){
uint32_t lo=(uint32_t)_mm256_movemask_epi8(_mm256_load_si256((const __m256i*)(row+k)));
uint32_t hi=(uint32_t)_mm256_movemask_epi8(_mm256_load_si256((const __m256i*)(row+k+32)));
w[k>>6]=(uint64_t)lo|((uint64_t)hi<<32);
}
}
}
AVX2 static void build_Bpbits(int n, const int8_t*B){
const __m256i one=_mm256_set1_epi8(1);
int W=n/64;
for (int k=0;k<n;k++){
const int8_t*row=B+(size_t)k*n; uint64_t*w=Bpb+(size_t)k*W;
for (int j=0;j<n;j+=64){
__m256i a=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j)),one),7);
__m256i b=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j+32)),one),7);
uint32_t lo=(uint32_t)_mm256_movemask_epi8(a);
uint32_t hi=(uint32_t)_mm256_movemask_epi8(b);
w[j>>6]=(uint64_t)lo|((uint64_t)hi<<32);
}
}
}
/* P[i] = XOR of Bp[k] over k with A1[i][k] set. NY = n/256 ymm per row. */
AVX2 static void corr_bits(int n){
int W=n/64, NY=W>>2;
for (int i=0;i<n;i++){
__m256i acc[16];
for(int t=0;t<NY;t++) acc[t]=_mm256_setzero_si256();
const uint64_t*row=A1b+(size_t)i*W;
for(int w=0;w<W;w++){
uint64_t b=row[w];
while(b){
int t=__builtin_ctzll(b); b&=b-1;
const uint64_t*br=Bpb+(size_t)(w*64+t)*W;
for(int q=0;q<NY;q++)
acc[q]=_mm256_xor_si256(acc[q],_mm256_load_si256((const __m256i*)(br+4*q)));
}
}
uint64_t*o=Pbb+(size_t)i*W;
for(int q=0;q<NY;q++) _mm256_store_si256((__m256i*)(o+4*q),acc[q]);
}
}
/* P[i][j] = XOR over k with A1[i][k] set of Bp[k][j].
Per-output-row bit scan (dense 64-bit words -> predictable), K-BLOCKED so the
Bp slice stays in L2. P accumulated in memory across k-blocks. */
AVX2 static void corr_pass3(int n){
int W=n/64, NY=W>>2;
memset(Pbb,0,(size_t)n*W*8);
const int KB=256;
for (int kb=0; kb<n; kb+=KB){
int w0=kb>>6, w1=w0+(KB>>6);
for (int i=0;i<n;i++){
__m256i acc[16];
for (int q=0;q<NY;q++) acc[q]=_mm256_setzero_si256();
const uint64_t*row=A1b+(size_t)i*W;
for (int w=w0;w<w1;w++){
uint64_t b=row[w];
while(b){
int t=__builtin_ctzll(b); b&=b-1;
const uint64_t*br=Bpb+(size_t)(w*64+t)*W;
for (int q=0;q<NY;q++) acc[q]=_mm256_xor_si256(acc[q],_mm256_load_si256((const __m256i*)(br+4*q)));
}
}
uint64_t*o=Pbb+(size_t)i*W;
for (int q=0;q<NY;q++){
__m256i*po=(__m256i*)(o+4*q);
_mm256_store_si256(po,_mm256_xor_si256(_mm256_load_si256((const __m256i*)po),acc[q]));
}
}
}
}
/* specialized for n=4096 (W=64, NY=16): fully unrolled XOR loops */
AVX2 static void corr16(void){
const int n=4096, W=64;
memset(Pbb,0,(size_t)n*W*8);
const int KB=256;
for (int kb=0; kb<n; kb+=KB){
int w0=kb>>6, w1=w0+4;
for (int i=0;i<n;i++){
__m256i acc0,acc1,acc2,acc3,acc4,acc5,acc6,acc7,acc8,acc9,acc10,acc11,acc12,acc13,acc14,acc15;
acc0=acc1=acc2=acc3=acc4=acc5=acc6=acc7=_mm256_setzero_si256();
acc8=acc9=acc10=acc11=acc12=acc13=acc14=acc15=_mm256_setzero_si256();
const uint64_t*row=A1b+(size_t)i*W;
for (int w=w0;w<w1;w++){
uint64_t b=row[w];
while(b){
int t=__builtin_ctzll(b); b&=b-1;
const uint64_t*br=Bpb+(size_t)(w*64+t)*W;
acc0 =_mm256_xor_si256(acc0 ,_mm256_load_si256((const __m256i*)(br+0)));
acc1 =_mm256_xor_si256(acc1 ,_mm256_load_si256((const __m256i*)(br+4)));
acc2 =_mm256_xor_si256(acc2 ,_mm256_load_si256((const __m256i*)(br+8)));
acc3 =_mm256_xor_si256(acc3 ,_mm256_load_si256((const __m256i*)(br+12)));
acc4 =_mm256_xor_si256(acc4 ,_mm256_load_si256((const __m256i*)(br+16)));
acc5 =_mm256_xor_si256(acc5 ,_mm256_load_si256((const __m256i*)(br+20)));
acc6 =_mm256_xor_si256(acc6 ,_mm256_load_si256((const __m256i*)(br+24)));
acc7 =_mm256_xor_si256(acc7 ,_mm256_load_si256((const __m256i*)(br+28)));
acc8 =_mm256_xor_si256(acc8 ,_mm256_load_si256((const __m256i*)(br+32)));
acc9 =_mm256_xor_si256(acc9 ,_mm256_load_si256((const __m256i*)(br+36)));
acc10=_mm256_xor_si256(acc10,_mm256_load_si256((const __m256i*)(br+40)));
acc11=_mm256_xor_si256(acc11,_mm256_load_si256((const __m256i*)(br+44)));
acc12=_mm256_xor_si256(acc12,_mm256_load_si256((const __m256i*)(br+48)));
acc13=_mm256_xor_si256(acc13,_mm256_load_si256((const __m256i*)(br+52)));
acc14=_mm256_xor_si256(acc14,_mm256_load_si256((const __m256i*)(br+56)));
acc15=_mm256_xor_si256(acc15,_mm256_load_si256((const __m256i*)(br+60)));
}
}
uint64_t*o=Pbb+(size_t)i*W;
#define RMW(q,a) { __m256i*po=(__m256i*)(o+4*(q)); _mm256_store_si256(po,_mm256_xor_si256(_mm256_load_si256((const __m256i*)po),a)); }
RMW(0,acc0) RMW(1,acc1) RMW(2,acc2) RMW(3,acc3) RMW(4,acc4) RMW(5,acc5) RMW(6,acc6) RMW(7,acc7)
RMW(8,acc8) RMW(9,acc9) RMW(10,acc10) RMW(11,acc11) RMW(12,acc12) RMW(13,acc13) RMW(14,acc14) RMW(15,acc15)
#undef RMW
}
}
}
AVX2 static void corr1k(void){
const int n=1024, W=16;
memset(Pbb,0,(size_t)n*W*8);
const int KB=256;
for (int kb=0; kb<n; kb+=KB){
int w0=kb>>6, w1=w0+4;
for (int i=0;i<n;i++){
__m256i acc0,acc1,acc2,acc3;
acc0=acc1=acc2=acc3=_mm256_setzero_si256();
const uint64_t*row=A1b+(size_t)i*W;
for (int w=w0;w<w1;w++){
uint64_t b=row[w];
while(b){
int t=__builtin_ctzll(b); b&=b-1;
const uint64_t*br=Bpb+(size_t)(w*64+t)*W;
acc0 =_mm256_xor_si256(acc0 ,_mm256_load_si256((const __m256i*)(br+0)));
acc1 =_mm256_xor_si256(acc1 ,_mm256_load_si256((const __m256i*)(br+4)));
acc2 =_mm256_xor_si256(acc2 ,_mm256_load_si256((const __m256i*)(br+8)));
acc3 =_mm256_xor_si256(acc3 ,_mm256_load_si256((const __m256i*)(br+12)));
}
}
uint64_t*o=Pbb+(size_t)i*W;
#define RMW(q,a) { __m256i*po=(__m256i*)(o+4*(q)); _mm256_store_si256(po,_mm256_xor_si256(_mm256_load_si256((const __m256i*)po),a)); }
RMW(0,acc0) RMW(1,acc1) RMW(2,acc2) RMW(3,acc3)
#undef RMW
}
}
}
static void expand_bits(int n){
static uint8_t T[256*8]; static int ti=0;
if(!ti){ for(int v=0;v<256;v++) for(int b=0;b<8;b++) T[v*8+b]=((v>>b)&1)?0x80:0x00; ti=1; }
int W=n/64;
for(int i=0;i<n;i++){
const uint64_t*w=Pbb+(size_t)i*W; uint8_t*o=PBb+(size_t)i*n;
for(int q=0;q<W;q++){
uint64_t v=w[q];
for(int s=0;s<8;s++) memcpy(o+(q*8+s)*8, T+(((uint8_t)(v>>(8*s)))*8), 8);
}
}
}
#define NIB 2
#define NJ 4
static const int16_t W0FF[16] __attribute__((aligned(32))) =
{0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,
0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF};
static int8_t BI3b[NMAX*NMAX] __attribute__((aligned(64)));
AVX2 static void build_BI3(int n, const int8_t*B){
int ng=n/64;
for (int p=0;p<n/2;p++){
const int8_t *r0=B+(size_t)(2*p)*n, *r1=r0+n;
for (int g=0;g<ng;g++){
const int8_t *a=r0+64*g, *b=r1+64*g;
int8_t *o=BI3b+((size_t)g*(n/2)+p)*128;
for (int t=0;t<64;t+=32){
__m256i va=_mm256_load_si256((const __m256i*)(a+t));
__m256i vb=_mm256_load_si256((const __m256i*)(b+t));
__m256i lo=_mm256_unpacklo_epi8(va,vb), hi=_mm256_unpackhi_epi8(va,vb);
_mm256_store_si256((__m256i*)(o+2*t), _mm256_permute2x128_si256(lo,hi,0x20));
_mm256_store_si256((__m256i*)(o+2*t+32), _mm256_permute2x128_si256(lo,hi,0x31));
}
}
}
}
static int8_t mscr[256] __attribute__((aligned(64)));
/* Bit-expansion table, ONE 8-byte entry per byte value: SHB8[v*8+m] = bit m of v ? 0x80 : 0.
Compile-time constant (the old runtime `shinit()` cost a flag load + guard + a vectorised
init block in the middle of the hot loop on EVERY epilogue call). */
static const uint8_t SHB8[256*8] __attribute__((aligned(64))) = {
0,0,0,0,0,0,0,0,0x80,0,0,0,0,0,0,0,
0,0x80,0,0,0,0,0,0,0x80,0x80,0,0,0,0,0,0,
0,0,0x80,0,0,0,0,0,0x80,0,0x80,0,0,0,0,0,
0,0x80,0x80,0,0,0,0,0,0x80,0x80,0x80,0,0,0,0,0,
0,0,0,0x80,0,0,0,0,0x80,0,0,0x80,0,0,0,0,
0,0x80,0,0x80,0,0,0,0,0x80,0x80,0,0x80,0,0,0,0,
0,0,0x80,0x80,0,0,0,0,0x80,0,0x80,0x80,0,0,0,0,
0,0x80,0x80,0x80,0,0,0,0,0x80,0x80,0x80,0x80,0,0,0,0,
0,0,0,0,0x80,0,0,0,0x80,0,0,0,0x80,0,0,0,
0,0x80,0,0,0x80,0,0,0,0x80,0x80,0,0,0x80,0,0,0,
0,0,0x80,0,0x80,0,0,0,0x80,0,0x80,0,0x80,0,0,0,
0,0x80,0x80,0,0x80,0,0,0,0x80,0x80,0x80,0,0x80,0,0,0,
0,0,0,0x80,0x80,0,0,0,0x80,0,0,0x80,0x80,0,0,0,
0,0x80,0,0x80,0x80,0,0,0,0x80,0x80,0,0x80,0x80,0,0,0,
0,0,0x80,0x80,0x80,0,0,0,0x80,0,0x80,0x80,0x80,0,0,0,
0,0x80,0x80,0x80,0x80,0,0,0,0x80,0x80,0x80,0x80,0x80,0,0,0,
0,0,0,0,0,0x80,0,0,0x80,0,0,0,0,0x80,0,0,
0,0x80,0,0,0,0x80,0,0,0x80,0x80,0,0,0,0x80,0,0,
0,0,0x80,0,0,0x80,0,0,0x80,0,0x80,0,0,0x80,0,0,
0,0x80,0x80,0,0,0x80,0,0,0x80,0x80,0x80,0,0,0x80,0,0,
0,0,0,0x80,0,0x80,0,0,0x80,0,0,0x80,0,0x80,0,0,
0,0x80,0,0x80,0,0x80,0,0,0x80,0x80,0,0x80,0,0x80,0,0,
0,0,0x80,0x80,0,0x80,0,0,0x80,0,0x80,0x80,0,0x80,0,0,
0,0x80,0x80,0x80,0,0x80,0,0,0x80,0x80,0x80,0x80,0,0x80,0,0,
0,0,0,0,0x80,0x80,0,0,0x80,0,0,0,0x80,0x80,0,0,
0,0x80,0,0,0x80,0x80,0,0,0x80,0x80,0,0,0x80,0x80,0,0,
0,0,0x80,0,0x80,0x80,0,0,0x80,0,0x80,0,0x80,0x80,0,0,
0,0x80,0x80,0,0x80,0x80,0,0,0x80,0x80,0x80,0,0x80,0x80,0,0,
0,0,0,0x80,0x80,0x80,0,0,0x80,0,0,0x80,0x80,0x80,0,0,
0,0x80,0,0x80,0x80,0x80,0,0,0x80,0x80,0,0x80,0x80,0x80,0,0,
0,0,0x80,0x80,0x80,0x80,0,0,0x80,0,0x80,0x80,0x80,0x80,0,0,
0,0x80,0x80,0x80,0x80,0x80,0,0,0x80,0x80,0x80,0x80,0x80,0x80,0,0,
0,0,0,0,0,0,0x80,0,0x80,0,0,0,0,0,0x80,0,
0,0x80,0,0,0,0,0x80,0,0x80,0x80,0,0,0,0,0x80,0,
0,0,0x80,0,0,0,0x80,0,0x80,0,0x80,0,0,0,0x80,0,
0,0x80,0x80,0,0,0,0x80,0,0x80,0x80,0x80,0,0,0,0x80,0,
0,0,0,0x80,0,0,0x80,0,0x80,0,0,0x80,0,0,0x80,0,
0,0x80,0,0x80,0,0,0x80,0,0x80,0x80,0,0x80,0,0,0x80,0,
0,0,0x80,0x80,0,0,0x80,0,0x80,0,0x80,0x80,0,0,0x80,0,
0,0x80,0x80,0x80,0,0,0x80,0,0x80,0x80,0x80,0x80,0,0,0x80,0,
0,0,0,0,0x80,0,0x80,0,0x80,0,0,0,0x80,0,0x80,0,
0,0x80,0,0,0x80,0,0x80,0,0x80,0x80,0,0,0x80,0,0x80,0,
0,0,0x80,0,0x80,0,0x80,0,0x80,0,0x80,0,0x80,0,0x80,0,
0,0x80,0x80,0,0x80,0,0x80,0,0x80,0x80,0x80,0,0x80,0,0x80,0,
0,0,0,0x80,0x80,0,0x80,0,0x80,0,0,0x80,0x80,0,0x80,0,
0,0x80,0,0x80,0x80,0,0x80,0,0x80,0x80,0,0x80,0x80,0,0x80,0,
0,0,0x80,0x80,0x80,0,0x80,0,0x80,0,0x80,0x80,0x80,0,0x80,0,
0,0x80,0x80,0x80,0x80,0,0x80,0,0x80,0x80,0x80,0x80,0x80,0,0x80,0,
0,0,0,0,0,0x80,0x80,0,0x80,0,0,0,0,0x80,0x80,0,
0,0x80,0,0,0,0x80,0x80,0,0x80,0x80,0,0,0,0x80,0x80,0,
0,0,0x80,0,0,0x80,0x80,0,0x80,0,0x80,0,0,0x80,0x80,0,
0,0x80,0x80,0,0,0x80,0x80,0,0x80,0x80,0x80,0,0,0x80,0x80,0,
0,0,0,0x80,0,0x80,0x80,0,0x80,0,0,0x80,0,0x80,0x80,0,
0,0x80,0,0x80,0,0x80,0x80,0,0x80,0x80,0,0x80,0,0x80,0x80,0,
0,0,0x80,0x80,0,0x80,0x80,0,0x80,0,0x80,0x80,0,0x80,0x80,0,
0,0x80,0x80,0x80,0,0x80,0x80,0,0x80,0x80,0x80,0x80,0,0x80,0x80,0,
0,0,0,0,0x80,0x80,0x80,0,0x80,0,0,0,0x80,0x80,0x80,0,
0,0x80,0,0,0x80,0x80,0x80,0,0x80,0x80,0,0,0x80,0x80,0x80,0,
0,0,0x80,0,0x80,0x80,0x80,0,0x80,0,0x80,0,0x80,0x80,0x80,0,
0,0x80,0x80,0,0x80,0x80,0x80,0,0x80,0x80,0x80,0,0x80,0x80,0x80,0,
0,0,0,0x80,0x80,0x80,0x80,0,0x80,0,0,0x80,0x80,0x80,0x80,0,
0,0x80,0,0x80,0x80,0x80,0x80,0,0x80,0x80,0,0x80,0x80,0x80,0x80,0,
0,0,0x80,0x80,0x80,0x80,0x80,0,0x80,0,0x80,0x80,0x80,0x80,0x80,0,
0,0x80,0x80,0x80,0x80,0x80,0x80,0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,
0,0,0,0,0,0,0,0x80,0x80,0,0,0,0,0,0,0x80,
0,0x80,0,0,0,0,0,0x80,0x80,0x80,0,0,0,0,0,0x80,
0,0,0x80,0,0,0,0,0x80,0x80,0,0x80,0,0,0,0,0x80,
0,0x80,0x80,0,0,0,0,0x80,0x80,0x80,0x80,0,0,0,0,0x80,
0,0,0,0x80,0,0,0,0x80,0x80,0,0,0x80,0,0,0,0x80,
0,0x80,0,0x80,0,0,0,0x80,0x80,0x80,0,0x80,0,0,0,0x80,
0,0,0x80,0x80,0,0,0,0x80,0x80,0,0x80,0x80,0,0,0,0x80,
0,0x80,0x80,0x80,0,0,0,0x80,0x80,0x80,0x80,0x80,0,0,0,0x80,
0,0,0,0,0x80,0,0,0x80,0x80,0,0,0,0x80,0,0,0x80,
0,0x80,0,0,0x80,0,0,0x80,0x80,0x80,0,0,0x80,0,0,0x80,
0,0,0x80,0,0x80,0,0,0x80,0x80,0,0x80,0,0x80,0,0,0x80,
0,0x80,0x80,0,0x80,0,0,0x80,0x80,0x80,0x80,0,0x80,0,0,0x80,
0,0,0,0x80,0x80,0,0,0x80,0x80,0,0,0x80,0x80,0,0,0x80,
0,0x80,0,0x80,0x80,0,0,0x80,0x80,0x80,0,0x80,0x80,0,0,0x80,
0,0,0x80,0x80,0x80,0,0,0x80,0x80,0,0x80,0x80,0x80,0,0,0x80,
0,0x80,0x80,0x80,0x80,0,0,0x80,0x80,0x80,0x80,0x80,0x80,0,0,0x80,
0,0,0,0,0,0x80,0,0x80,0x80,0,0,0,0,0x80,0,0x80,
0,0x80,0,0,0,0x80,0,0x80,0x80,0x80,0,0,0,0x80,0,0x80,
0,0,0x80,0,0,0x80,0,0x80,0x80,0,0x80,0,0,0x80,0,0x80,
0,0x80,0x80,0,0,0x80,0,0x80,0x80,0x80,0x80,0,0,0x80,0,0x80,
0,0,0,0x80,0,0x80,0,0x80,0x80,0,0,0x80,0,0x80,0,0x80,
0,0x80,0,0x80,0,0x80,0,0x80,0x80,0x80,0,0x80,0,0x80,0,0x80,
0,0,0x80,0x80,0,0x80,0,0x80,0x80,0,0x80,0x80,0,0x80,0,0x80,
0,0x80,0x80,0x80,0,0x80,0,0x80,0x80,0x80,0x80,0x80,0,0x80,0,0x80,
0,0,0,0,0x80,0x80,0,0x80,0x80,0,0,0,0x80,0x80,0,0x80,
0,0x80,0,0,0x80,0x80,0,0x80,0x80,0x80,0,0,0x80,0x80,0,0x80,
0,0,0x80,0,0x80,0x80,0,0x80,0x80,0,0x80,0,0x80,0x80,0,0x80,
0,0x80,0x80,0,0x80,0x80,0,0x80,0x80,0x80,0x80,0,0x80,0x80,0,0x80,
0,0,0,0x80,0x80,0x80,0,0x80,0x80,0,0,0x80,0x80,0x80,0,0x80,
0,0x80,0,0x80,0x80,0x80,0,0x80,0x80,0x80,0,0x80,0x80,0x80,0,0x80,
0,0,0x80,0x80,0x80,0x80,0,0x80,0x80,0,0x80,0x80,0x80,0x80,0,0x80,
0,0x80,0x80,0x80,0x80,0x80,0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,0x80,
0,0,0,0,0,0,0x80,0x80,0x80,0,0,0,0,0,0x80,0x80,
0,0x80,0,0,0,0,0x80,0x80,0x80,0x80,0,0,0,0,0x80,0x80,
0,0,0x80,0,0,0,0x80,0x80,0x80,0,0x80,0,0,0,0x80,0x80,
0,0x80,0x80,0,0,0,0x80,0x80,0x80,0x80,0x80,0,0,0,0x80,0x80,
0,0,0,0x80,0,0,0x80,0x80,0x80,0,0,0x80,0,0,0x80,0x80,
0,0x80,0,0x80,0,0,0x80,0x80,0x80,0x80,0,0x80,0,0,0x80,0x80,
0,0,0x80,0x80,0,0,0x80,0x80,0x80,0,0x80,0x80,0,0,0x80,0x80,
0,0x80,0x80,0x80,0,0,0x80,0x80,0x80,0x80,0x80,0x80,0,0,0x80,0x80,
0,0,0,0,0x80,0,0x80,0x80,0x80,0,0,0,0x80,0,0x80,0x80,
0,0x80,0,0,0x80,0,0x80,0x80,0x80,0x80,0,0,0x80,0,0x80,0x80,
0,0,0x80,0,0x80,0,0x80,0x80,0x80,0,0x80,0,0x80,0,0x80,0x80,
0,0x80,0x80,0,0x80,0,0x80,0x80,0x80,0x80,0x80,0,0x80,0,0x80,0x80,
0,0,0,0x80,0x80,0,0x80,0x80,0x80,0,0,0x80,0x80,0,0x80,0x80,
0,0x80,0,0x80,0x80,0,0x80,0x80,0x80,0x80,0,0x80,0x80,0,0x80,0x80,
0,0,0x80,0x80,0x80,0,0x80,0x80,0x80,0,0x80,0x80,0x80,0,0x80,0x80,
0,0x80,0x80,0x80,0x80,0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,0x80,0x80,
0,0,0,0,0,0x80,0x80,0x80,0x80,0,0,0,0,0x80,0x80,0x80,
0,0x80,0,0,0,0x80,0x80,0x80,0x80,0x80,0,0,0,0x80,0x80,0x80,
0,0,0x80,0,0,0x80,0x80,0x80,0x80,0,0x80,0,0,0x80,0x80,0x80,
0,0x80,0x80,0,0,0x80,0x80,0x80,0x80,0x80,0x80,0,0,0x80,0x80,0x80,
0,0,0,0x80,0,0x80,0x80,0x80,0x80,0,0,0x80,0,0x80,0x80,0x80,
0,0x80,0,0x80,0,0x80,0x80,0x80,0x80,0x80,0,0x80,0,0x80,0x80,0x80,
0,0,0x80,0x80,0,0x80,0x80,0x80,0x80,0,0x80,0x80,0,0x80,0x80,0x80,
0,0x80,0x80,0x80,0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,0x80,0x80,0x80,
0,0,0,0,0x80,0x80,0x80,0x80,0x80,0,0,0,0x80,0x80,0x80,0x80,
0,0x80,0,0,0x80,0x80,0x80,0x80,0x80,0x80,0,0,0x80,0x80,0x80,0x80,
0,0,0x80,0,0x80,0x80,0x80,0x80,0x80,0,0x80,0,0x80,0x80,0x80,0x80,
0,0x80,0x80,0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,0x80,0x80,0x80,0x80,
0,0,0,0x80,0x80,0x80,0x80,0x80,0x80,0,0,0x80,0x80,0x80,0x80,0x80,
0,0x80,0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,0x80,0x80,0x80,0x80,0x80,
0,0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0,0x80,0x80,0x80,0x80,0x80,0x80,
0,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,
};
/* [sm2_] SIMD bit-expansion constants for the epilogue.
EPD_IDX[c] = c div 8 -> after vpshufb over a vpbroadcastd of the four Pbb
bytes, byte c of the spread holds the Pbb byte that owns column c.
EPD_BIT[c] = 1 << (c mod 8) -> isolates the one bit that decides column c.
EPD_80 = 0x80 everywhere -> turns the vpcmpeqb result into the xor mask. */
static const uint8_t EPD_IDX[32] __attribute__((aligned(32))) = {
0,0,0,0,0,0,0,0, 1,1,1,1,1,1,1,1, 2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3 };
static const uint8_t EPD_BIT[32] __attribute__((aligned(32))) = {
1,2,4,8,16,32,64,128, 1,2,4,8,16,32,64,128,
1,2,4,8,16,32,64,128, 1,2,4,8,16,32,64,128 };
static const uint8_t EPD_80[32] __attribute__((aligned(32))) = {
0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80, 0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,
0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80, 0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80 };
AVX2 static void epilogue(int n, int8_t*C, const uint64_t*P, int jb, int i, int NCH){
const size_t bo=(size_t)jb>>3;
const uint8_t*base=(const uint8_t*)P+(size_t)i*(n>>3);
const __m256i idx=_mm256_load_si256((const __m256i*)EPD_IDX);
const __m256i bit=_mm256_load_si256((const __m256i*)EPD_BIT);
const __m256i m80=_mm256_load_si256((const __m256i*)EPD_80);
for(int r=0;r<2;r++){
const uint8_t*qr=base+(size_t)r*(n>>3);
for(int t=0;t<NCH;t++){
/* chunk (r,t) covers C row i+r, columns jb+32t .. jb+32t+31, and reads
the four Pbb bytes at qr[bo+4t .. bo+4t+3]. Pbb bit position p of the
row is the correction for column p, and p = (bo+4t+s)*8 + b with
column offset c = 8s+b => column c takes byte (c div 8) and bit
(c mod 8) of those four bytes. vpbroadcastd makes byte c equal to
source byte (c mod 4); vpshufb then selects index (c div 8) and,
because the dword is replicated in both 128-bit lanes, lane 1 can
reach indices 2 and 3 of its own copy. Result: byte c of the spread
is the Pbb byte that owns column c, and vpand/vpcmpeqb/vpand turn
bit (c mod 8) of it into an 0x80 xor mask. */
__m256i v=_mm256_load_si256((const __m256i*)(mscr+(r*NCH+t)*32));
int wd; __builtin_memcpy(&wd, qr+bo+4*(size_t)t, 4);
__m256i b=_mm256_broadcastd_epi32(_mm_cvtsi32_si128(wd));
__m256i s=_mm256_shuffle_epi8(b, idx);
__m256i z=_mm256_and_si256(s, bit);
z=_mm256_and_si256(_mm256_cmpeq_epi8(z, bit), m80);
int8_t*c=C+(size_t)(i+r)*n+jb+t*32;
_mm256_storeu_si256((__m256i*)c,_mm256_xor_si256(v,z));
}
}
}
/* ---- C=96 (6 BI ymm) mikro: 2 rows x 96 columns, records of 192 bytes ---- */
static int8_t BI9b[10*(NMAX/2)*192] __attribute__((aligned(64)));
/* padded Bp layout: 65 rows per 64-k group, row 64 = ZERO (safe index for tzcnt(0)=64) */
static uint64_t BP2[16*65*16] __attribute__((aligned(64)));
AVX2 static void build_BP2(int n, const int8_t*B){
const int W=n>>6, STR=65*W;
const __m256i one=_mm256_set1_epi8(1);
for (int k=0;k<n;k++){
const int8_t*row=B+(size_t)k*n; uint64_t*w=BP2+(size_t)(k>>6)*STR+(size_t)(k&63)*W;
for (int j=0;j<n;j+=64){
__m256i a=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j)),one),7);
__m256i b=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j+32)),one),7);
uint32_t lo=(uint32_t)_mm256_movemask_epi8(a);
uint32_t hi=(uint32_t)_mm256_movemask_epi8(b);
w[j>>6]=(uint64_t)lo|((uint64_t)hi<<32);
}
}
for (int g=0; g<(n>>6); g++) memset(BP2+(size_t)g*STR+64*W, 0, W*8);
}
/* ---- 4-word merged bit scan -------------------------------------------
#100989/#101030 scan ONE pair of 64-bit words per `while (b0|b1)` loop; the
data-dependent trip count (max(popcount(b0),popcount(b1)) ~ 33 of 64) is what
the loop really pays for. Merging the kb block's four words into ONE loop
touches exactly the same bytes with exactly the same loads/XORs, but issues
half as many while-instances (2 -> 1 per (i,kb)), hence half the
data-dependent exits and half the per-instance fixup (row[w] loads, bw
pointers, for-w bookkeeping). The tzcnt(0)=64 sentinel row keeps the shorter
words harmless, as before. Judge probe: -7.08% of the phase, Pbb
byte-identical. */
AVX2 static void corr2k(int n){
const int W=n>>6, STR=65*W;
memset(Pbb,0,(size_t)n*W*8);
for (int kb0=0; kb0<n; kb0+=256){
int w0=kb0>>6;
for (int i=0;i<n;i++){
__m256i a0,a1,a2,a3; a0=a1=a2=a3=_mm256_setzero_si256();
const uint64_t*row=A1b+(size_t)i*W;
uint64_t b0=row[w0],b1=row[w0+1],b2=row[w0+2],b3=row[w0+3];
const uint64_t*bw0=BP2+(size_t)w0*STR;
const uint64_t*bw1=bw0+STR,*bw2=bw1+STR,*bw3=bw2+STR;
while (b0|b1|b2|b3){
long t0,t1,t2,t3; uint64_t nb;
__asm__("tzcntq %2,%0\n\t" "blsrq %2,%1":"=r"(t0),"=r"(nb):"r"(b0)); b0=nb;
__asm__("tzcntq %2,%0\n\t" "blsrq %2,%1":"=r"(t1),"=r"(nb):"r"(b1)); b1=nb;
__asm__("tzcntq %2,%0\n\t" "blsrq %2,%1":"=r"(t2),"=r"(nb):"r"(b2)); b2=nb;
__asm__("tzcntq %2,%0\n\t" "blsrq %2,%1":"=r"(t3),"=r"(nb):"r"(b3)); b3=nb;
const uint64_t*r0=bw0+((size_t)t0<<4);
const uint64_t*r1=bw1+((size_t)t1<<4);
const uint64_t*r2=bw2+((size_t)t2<<4);
const uint64_t*r3=bw3+((size_t)t3<<4);
a0=_mm256_xor_si256(a0,_mm256_load_si256((const __m256i*)(r0+0)));
a1=_mm256_xor_si256(a1,_mm256_load_si256((const __m256i*)(r0+4)));
a2=_mm256_xor_si256(a2,_mm256_load_si256((const __m256i*)(r0+8)));
a3=_mm256_xor_si256(a3,_mm256_load_si256((const __m256i*)(r0+12)));
a0=_mm256_xor_si256(a0,_mm256_load_si256((const __m256i*)(r1+0)));
a1=_mm256_xor_si256(a1,_mm256_load_si256((const __m256i*)(r1+4)));
a2=_mm256_xor_si256(a2,_mm256_load_si256((const __m256i*)(r1+8)));
a3=_mm256_xor_si256(a3,_mm256_load_si256((const __m256i*)(r1+12)));
a0=_mm256_xor_si256(a0,_mm256_load_si256((const __m256i*)(r2+0)));
a1=_mm256_xor_si256(a1,_mm256_load_si256((const __m256i*)(r2+4)));
a2=_mm256_xor_si256(a2,_mm256_load_si256((const __m256i*)(r2+8)));
a3=_mm256_xor_si256(a3,_mm256_load_si256((const __m256i*)(r2+12)));
a0=_mm256_xor_si256(a0,_mm256_load_si256((const __m256i*)(r3+0)));
a1=_mm256_xor_si256(a1,_mm256_load_si256((const __m256i*)(r3+4)));
a2=_mm256_xor_si256(a2,_mm256_load_si256((const __m256i*)(r3+8)));
a3=_mm256_xor_si256(a3,_mm256_load_si256((const __m256i*)(r3+12)));
}
uint64_t*o=Pbb+(size_t)i*W;
{ __m256i*po=(__m256i*)(o+0); _mm256_store_si256(po,_mm256_xor_si256(_mm256_load_si256((const __m256i*)po),a0)); }
{ __m256i*po=(__m256i*)(o+4); _mm256_store_si256(po,_mm256_xor_si256(_mm256_load_si256((const __m256i*)po),a1)); }
{ __m256i*po=(__m256i*)(o+8); _mm256_store_si256(po,_mm256_xor_si256(_mm256_load_si256((const __m256i*)po),a2)); }
{ __m256i*po=(__m256i*)(o+12); _mm256_store_si256(po,_mm256_xor_si256(_mm256_load_si256((const __m256i*)po),a3)); }
}
}
}
/* p OUTER / g inner: the whole of B is read ONCE, and for a fixed p the ten
96-byte column slices r0+96*g, g=0..9, cover bytes 0..959 of the row in order.
(#100989 had g outer / p inner => every g walked all of B => B read 10x.) */
AVX2 static void build_BI9(int n, const int8_t*B){
const int half=n>>1, NC=96, NREC=192;
for (int p=0;p<half;p++){
const int8_t*r0b=B+(size_t)(2*p)*n, *r1b=r0b+n;
for (int g=0; g<10; g++){
const int8_t*r0=r0b+NC*g, *r1=r1b+NC*g;
int8_t*o=BI9b+((size_t)g*half+p)*NREC;
for (int t=0;t<NC;t+=32){
__m256i va=_mm256_loadu_si256((const __m256i*)(r0+t));
__m256i vb=_mm256_loadu_si256((const __m256i*)(r1+t));
__m256i lo=_mm256_unpacklo_epi8(va,vb), hi=_mm256_unpackhi_epi8(va,vb);
_mm256_store_si256((__m256i*)(o+2*t), _mm256_permute2x128_si256(lo,hi,0x20));
_mm256_store_si256((__m256i*)(o+2*t+32), _mm256_permute2x128_si256(lo,hi,0x31));
}
}
}
}
/* fused A0 + A1bits: one pass over A instead of two (A = 1 MB > judge L2 = 256 KB) */
AVX2 static void build_A0_A1(int n, const int8_t*A){
const __m256i m=_mm256_set1_epi8(0x7F);
int W=n/64;
for (int i=0;i<n;i++){
const int8_t*row=A+(size_t)i*n; uint64_t*w=A1b+(size_t)i*W;
int8_t*o=A0b+2*(size_t)i*n;
for (int k=0;k<n;k+=64){
__m256i v0=_mm256_load_si256((const __m256i*)(row+k));
__m256i v1=_mm256_load_si256((const __m256i*)(row+k+32));
{ __m256i e0=_mm256_and_si256(v0,m), e1=_mm256_and_si256(v1,m);
__m256i l0=_mm256_unpacklo_epi16(e0,e0), h0=_mm256_unpackhi_epi16(e0,e0);
__m256i l1=_mm256_unpacklo_epi16(e1,e1), h1=_mm256_unpackhi_epi16(e1,e1);
_mm256_store_si256((__m256i*)(o+2*k), _mm256_permute2x128_si256(l0,h0,0x20));
_mm256_store_si256((__m256i*)(o+2*k+32), _mm256_permute2x128_si256(l0,h0,0x31));
_mm256_store_si256((__m256i*)(o+2*k+64), _mm256_permute2x128_si256(l1,h1,0x20));
_mm256_store_si256((__m256i*)(o+2*k+96), _mm256_permute2x128_si256(l1,h1,0x31)); }
uint32_t lo=(uint32_t)_mm256_movemask_epi8(v0);
uint32_t hi=(uint32_t)_mm256_movemask_epi8(v1);
w[k>>6]=(uint64_t)lo|((uint64_t)hi<<32);
}
}
}
/* tail block: 64 columns (960..1023) into BI3b with the original 128B records */
AVX2 static void build_BI3_tail(int n, const int8_t*B){
const int half=n>>1, c0=960;
for (int p=0;p<half;p++){
const int8_t*r0=B+(size_t)(2*p)*n+c0, *r1=r0+n;
int8_t*o=BI3b+((size_t)p)*128;
for (int t=0;t<64;t+=32){
__m256i va=_mm256_loadu_si256((const __m256i*)(r0+t));
__m256i vb=_mm256_loadu_si256((const __m256i*)(r1+t));
__m256i lo=_mm256_unpacklo_epi8(va,vb), hi=_mm256_unpackhi_epi8(va,vb);
_mm256_store_si256((__m256i*)(o+2*t), _mm256_permute2x128_si256(lo,hi,0x20));
_mm256_store_si256((__m256i*)(o+2*t+32), _mm256_permute2x128_si256(lo,hi,0x31));
}
}
}
#define W0FFS {0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF,0x00FF}
static const int16_t M9[16] __attribute__((aligned(32))) = W0FFS;
AVX2 static void mk_s6(int n, const int8_t*A0, const int8_t*BI, uint8_t*scr){
const int8_t *x0=A0, *x1=A0+n, *bi=BI; long cnt=n>>3;
asm volatile(
"vpxor %%xmm2,%%xmm2,%%xmm2\n\t"
"vpxor %%xmm3,%%xmm3,%%xmm3\n\t"
"vpxor %%xmm6,%%xmm6,%%xmm6\n\t"
"vpxor %%xmm7,%%xmm7,%%xmm7\n\t"
"vpxor %%xmm8,%%xmm8,%%xmm8\n\t"
"vpxor %%xmm9,%%xmm9,%%xmm9\n\t"
"vpxor %%xmm10,%%xmm10,%%xmm10\n\t"
"vpxor %%xmm11,%%xmm11,%%xmm11\n\t"
"vpxor %%xmm12,%%xmm12,%%xmm12\n\t"
"vpxor %%xmm13,%%xmm13,%%xmm13\n\t"
"vpxor %%xmm14,%%xmm14,%%xmm14\n\t"
"vpxor %%xmm15,%%xmm15,%%xmm15\n\t"
".p2align 5\n" "1:\n"
"prefetcht0 2304(%[bi])\n\t"
"vpbroadcastw (%[x0]), %%ymm0\n\t"
"vpbroadcastw (%[x1]), %%ymm1\n\t"
"vmovdqa (%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm14, %%ymm14\n\t"
"vmovdqa 32(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm15, %%ymm15\n\t"
"vmovdqa 64(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm2, %%ymm2\n\t"
"vmovdqa 96(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm3, %%ymm3\n\t"
"vmovdqa 128(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm6, %%ymm6\n\t"
"vmovdqa 160(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm13, %%ymm13\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm7, %%ymm7\n\t"
"vpbroadcastw 2(%[x0]), %%ymm0\n\t"
"vpbroadcastw 2(%[x1]), %%ymm1\n\t"
"vmovdqa 192(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm14, %%ymm14\n\t"
"vmovdqa 224(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm15, %%ymm15\n\t"
"vmovdqa 256(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm2, %%ymm2\n\t"
"vmovdqa 288(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm3, %%ymm3\n\t"
"vmovdqa 320(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm6, %%ymm6\n\t"
"vmovdqa 352(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm13, %%ymm13\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm7, %%ymm7\n\t"
"vpbroadcastw 4(%[x0]), %%ymm0\n\t"
"vpbroadcastw 4(%[x1]), %%ymm1\n\t"
"vmovdqa 384(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm14, %%ymm14\n\t"
"vmovdqa 416(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm15, %%ymm15\n\t"
"vmovdqa 448(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm2, %%ymm2\n\t"
"vmovdqa 480(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm3, %%ymm3\n\t"
"vmovdqa 512(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm6, %%ymm6\n\t"
"vmovdqa 544(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm13, %%ymm13\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm7, %%ymm7\n\t"
"vpbroadcastw 6(%[x0]), %%ymm0\n\t"
"vpbroadcastw 6(%[x1]), %%ymm1\n\t"
"vmovdqa 576(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm14, %%ymm14\n\t"
"vmovdqa 608(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm15, %%ymm15\n\t"
"vmovdqa 640(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm2, %%ymm2\n\t"
"vmovdqa 672(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm3, %%ymm3\n\t"
"vmovdqa 704(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm6, %%ymm6\n\t"
"vmovdqa 736(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm13, %%ymm13\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm7, %%ymm7\n\t"
"prefetcht0 800(%[bi])\n\t"
"prefetcht0 1600(%[bi])\n\t"
"addq $8, %[x0]\n\t"
"addq $8, %[x1]\n\t"
"addq $768, %[bi]\n\t"
"decq %[cnt]\n\t"
"jnz 1b\n\t"
"vpand %[mask], %%ymm8, %%ymm0\n\t"
"vpand %[mask], %%ymm9, %%ymm1\n\t"
"vpackuswb %%ymm1, %%ymm0, %%ymm0\n\t"
"vpermq $0xD8, %%ymm0, %%ymm0\n\t"
"vmovdqu %%ymm0, 0(%[scr])\n\t"
"vpand %[mask], %%ymm10, %%ymm0\n\t"
"vpand %[mask], %%ymm11, %%ymm1\n\t"
"vpackuswb %%ymm1, %%ymm0, %%ymm0\n\t"
"vpermq $0xD8, %%ymm0, %%ymm0\n\t"
"vmovdqu %%ymm0, 32(%[scr])\n\t"
"vpand %[mask], %%ymm12, %%ymm0\n\t"
"vpand %[mask], %%ymm13, %%ymm1\n\t"
"vpackuswb %%ymm1, %%ymm0, %%ymm0\n\t"
"vpermq $0xD8, %%ymm0, %%ymm0\n\t"
"vmovdqu %%ymm0, 64(%[scr])\n\t"
"vpand %[mask], %%ymm14, %%ymm0\n\t"
"vpand %[mask], %%ymm15, %%ymm1\n\t"
"vpackuswb %%ymm1, %%ymm0, %%ymm0\n\t"
"vpermq $0xD8, %%ymm0, %%ymm0\n\t"
"vmovdqu %%ymm0, 96(%[scr])\n\t"
"vpand %[mask], %%ymm2, %%ymm0\n\t"
"vpand %[mask], %%ymm3, %%ymm1\n\t"
"vpackuswb %%ymm1, %%ymm0, %%ymm0\n\t"
"vpermq $0xD8, %%ymm0, %%ymm0\n\t"
"vmovdqu %%ymm0, 128(%[scr])\n\t"
"vpand %[mask], %%ymm6, %%ymm0\n\t"
"vpand %[mask], %%ymm7, %%ymm1\n\t"
"vpackuswb %%ymm1, %%ymm0, %%ymm0\n\t"
"vpermq $0xD8, %%ymm0, %%ymm0\n\t"
"vmovdqu %%ymm0, 160(%[scr])\n\t"
: [x0]"+r"(x0),[x1]"+r"(x1),[bi]"+r"(bi),[cnt]"+r"(cnt)
: [mask]"m"(M9), [scr]"r"(scr)
: "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15","cc","memory");
}
AVX2 static void mk_s(int n, const int8_t*A0, const int8_t*BI, int8_t*C, const uint8_t*PB, int jb){ (void)PB;
const int8_t *x0=A0, *x1=A0+2*n, *bi=BI; long cnt=n>>4;
asm volatile(
"vpxor %%xmm8,%%xmm8,%%xmm8\n\t" "vpxor %%xmm9,%%xmm9,%%xmm9\n\t" "vpxor %%xmm10,%%xmm10,%%xmm10\n\t" "vpxor %%xmm11,%%xmm11,%%xmm11\n\t" "vpxor %%xmm12,%%xmm12,%%xmm12\n\t" "vpxor %%xmm13,%%xmm13,%%xmm13\n\t" "vpxor %%xmm14,%%xmm14,%%xmm14\n\t" "vpxor %%xmm15,%%xmm15,%%xmm15\n\t"
".p2align 5\n\t" "1:\n\t"
"prefetcht0 1536(%[bi])\n\t"
"vpbroadcastd (%[x0]), %%ymm0\n\t"
"vpbroadcastd (%[x1]), %%ymm1\n\t"
"vpmaddubsw 0(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 0(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 32(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 32(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 64(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 64(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 96(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 96(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 4(%[x0]), %%ymm0\n\t"
"vpbroadcastd 4(%[x1]), %%ymm1\n\t"
"vpmaddubsw 128(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 128(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 160(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 160(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 192(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 192(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 224(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 224(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 8(%[x0]), %%ymm0\n\t"
"vpbroadcastd 8(%[x1]), %%ymm1\n\t"
"vpmaddubsw 256(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 256(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 288(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 288(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 320(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 320(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 352(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 352(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 12(%[x0]), %%ymm0\n\t"
"vpbroadcastd 12(%[x1]), %%ymm1\n\t"
"vpmaddubsw 384(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 384(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 416(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 416(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 448(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 448(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 480(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 480(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 16(%[x0]), %%ymm0\n\t"
"vpbroadcastd 16(%[x1]), %%ymm1\n\t"
"vpmaddubsw 512(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 512(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 544(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 544(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 576(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 576(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 608(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 608(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 20(%[x0]), %%ymm0\n\t"
"vpbroadcastd 20(%[x1]), %%ymm1\n\t"
"vpmaddubsw 640(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 640(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 672(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 672(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 704(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 704(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 736(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 736(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 24(%[x0]), %%ymm0\n\t"
"vpbroadcastd 24(%[x1]), %%ymm1\n\t"
"vpmaddubsw 768(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 768(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 800(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 800(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 832(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 832(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 864(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 864(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 28(%[x0]), %%ymm0\n\t"
"vpbroadcastd 28(%[x1]), %%ymm1\n\t"
"vpmaddubsw 896(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 896(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 928(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 928(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 960(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 960(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 992(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 992(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"addq $32, %[x0]\n\t"
"addq $32, %[x1]\n\t"
"addq $1024, %[bi]\n\t"
"decq %[cnt]\n\t"
"jnz 1b\n\t"
"vpand %[mask], %%ymm8, %%ymm2\n\t"
"vpand %[mask], %%ymm9, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 0(%[scr])\n\t"
"vpand %[mask], %%ymm10, %%ymm2\n\t"
"vpand %[mask], %%ymm11, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 32(%[scr])\n\t"
"vpand %[mask], %%ymm12, %%ymm2\n\t"
"vpand %[mask], %%ymm13, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 64(%[scr])\n\t"
"vpand %[mask], %%ymm14, %%ymm2\n\t"
"vpand %[mask], %%ymm15, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 96(%[scr])\n\t"
: [x0]"+r"(x0),[x1]"+r"(x1),[bi]"+r"(bi),[cnt]"+r"(cnt)
: [mask]"m"(W0FF), [scr]"r"(mscr)
: "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15","cc","memory");
}
static uint8_t Gtab[8*256*(NMAX/8)] __attribute__((aligned(64)));
/* single-instruction 32B memory ops: gcc 9.3 -O2 otherwise splits
_mm256_loadu_si256 into vmovdqu xmm + vinserti128 (extra uop on p5). */
#define LDU32(p) ({ __m256i _v; __asm__("vmovdqu %1, %0" : "=x"(_v) : "m"(*(const __m256i*)(p))); _v; })
#define XOR32(dst,p) __asm__("vpxor %1, %0, %0" : "+x"(dst) : "m"(*(const __m256i*)(p)))
#define STU32(p,v) __asm__ volatile("vmovdqu %1, %0" : "+m"(*(__m256i*)(p)) : "x"(v))
template<int NY>
__attribute__((target("avx2"))) static void corr_rows(int n,int B,int W,int gb0,const uint8_t*A1){
for(int i=0;i<n;i++){
const uint8_t*a=A1+(size_t)i*B;
uint64_t*po=Pbb+(size_t)i*W;
__m256i acc[NY];
for(int q=0;q<NY;q++) acc[q]=_mm256_setzero_si256();
for(int g=0;g<4;g++){
const uint8_t*T=Gtab+((size_t)g*256+(size_t)a[gb0+g])*B;
for(int q=0;q<NY;q++) XOR32(acc[q],T+(size_t)q*32);
}
for(int q=0;q<NY;q++){
__m256i o=LDU32(po+4*(size_t)q);
STU32(po+4*(size_t)q, _mm256_xor_si256(o,acc[q]));
}
}
}
__attribute__((target("avx2"))) static void corr_rus(int n){
const int W=n/64, B=n/8, NY=B/32;
if(NY!=1 && NY!=2 && NY!=4 && NY!=8 && NY!=16){ corr_pass3(n); return; }
const uint8_t*A1=(const uint8_t*)A1b;
memset(Pbb,0,(size_t)n*W*8);
const int NB=4; /* 8 k-groups of 8 = 64 k per block */
for(int kb=0;kb<n;kb+=NB*8){
for(int g=0;g<NB;g++){
const int k0=kb+g*8;
uint8_t*T0=Gtab+(size_t)g*256*B;
memset(T0,0,B);
for(int v=1;v<256;v++){
uint8_t*dst=T0+(size_t)v*B;
const uint8_t*pr=T0+(size_t)(v&(v-1))*B;
const int k=k0+__builtin_ctz(v);
const uint8_t*bp=(const uint8_t*)(BP2+(size_t)(k>>6)*(65*W)+(size_t)(k&63)*W);
for(int q=0;q<B;q+=32){
__m256i x=LDU32(pr+q);
XOR32(x,bp+q);
STU32(dst+q,x);
}
}
}
const int gb0=kb>>3;
if(NY==16) corr_rows<16>(n,B,W,gb0,A1);
else if(NY==8)corr_rows<8>(n,B,W,gb0,A1);
else if(NY==4)corr_rows<4>(n,B,W,gb0,A1);
else if(NY==2)corr_rows<2>(n,B,W,gb0,A1);
else corr_rows<1>(n,B,W,gb0,A1);
}
}
/* ================= mc4_: ONE STRASSEN LEVEL (n=1024 -> 7 x 512) ======================
Exact over Z/256: every coefficient is +-1, no division, so the CLRS identity holds in
any commutative ring and int8 arithmetic mod 256 is exactly that ring.
The leaf (m x m = 512 x 512) runs the SHIPPED generic 64-column path
(`mk_s` + `epilogue(...,2)`); its panels are built by a FUSED packer that reads the two
parent blocks and emits A0b/A1b (resp. BI3b/Bpb) in ONE traversal -- nothing is
materialised (BRIEF 2.18.590, DIF). `nn` is the parent row stride (1024 here). */
struct SBLK { const int8_t *X, *Y; int r0, c0, r1, c1, op; };
AVX2 static inline __m256i sblk_cmb(__m256i v, const int8_t *ya, int op){
if (!op) return v;
__m256i u = _mm256_load_si256((const __m256i*)ya);
return (op==1) ? _mm256_add_epi8(v,u) : _mm256_sub_epi8(v,u);
}
/* bld_A = blk + build_A0_A1 fused: one pass emits A0b (=v&0x7F) and A1b (=sign). */
AVX2 static void bld_A(int m, int nn, const SBLK*b){
const __m256i msk=_mm256_set1_epi8(0x7F);
int W=m/64;
for (int i=0;i<m;i++){
const int8_t*xa=b->X+(size_t)(b->r0+i)*nn+b->c0;
const int8_t*ya=b->Y+(size_t)(b->r1+i)*nn+b->c1;
int8_t*dst=A0b+(size_t)i*m;
uint64_t*w1=A1b+(size_t)i*W;
for (int k=0;k<m;k+=64){
__m256i v0=_mm256_load_si256((const __m256i*)(xa+k));
__m256i v1=_mm256_load_si256((const __m256i*)(xa+k+32));
v0=sblk_cmb(v0,ya+k, b->op);
v1=sblk_cmb(v1,ya+k+32,b->op);
_mm256_store_si256((__m256i*)(dst+k), _mm256_and_si256(v0,msk));
_mm256_store_si256((__m256i*)(dst+k+32), _mm256_and_si256(v1,msk));
w1[k>>6]=(uint64_t)(uint32_t)_mm256_movemask_epi8(v0)
|((uint64_t)(uint32_t)_mm256_movemask_epi8(v1)<<32);
}
}
}
/* bld_B = blk + build_BI3 + build_Bpbits fused: one pass emits BI3b and Bpb(=B&1). */
AVX2 static void bld_B(int m, int nn, const SBLK*b){
const __m256i one=_mm256_set1_epi8(1);
int ng=m/64, W=m/64;
for (int p=0;p<m/2;p++){
const int8_t *ax=b->X+(size_t)(b->r0+2*p)*nn+b->c0;
const int8_t *bx=ax+nn;
const int8_t *ay=b->Y+(size_t)(b->r1+2*p)*nn+b->c1;
const int8_t *by=ay+nn;
/* corr_rus reads the PADDED panel BP2, not Bpb: group stride 65*W with a zeroed
sentinel row at index 64 (see build_BP2). Rows 2p and 2p+1 can straddle a group
boundary, so each row's offset is computed independently. */
uint64_t *w0p=BP2+(size_t)((2*p)>>6)*(size_t)(65*W)+(size_t)((2*p)&63)*W;
uint64_t *w1p=BP2+(size_t)((2*p+1)>>6)*(size_t)(65*W)+(size_t)((2*p+1)&63)*W;
for (int g=0;g<ng;g++){
const int8_t *a=ax+64*g, *c=bx+64*g;
const int8_t *d=ay+64*g, *e=by+64*g;
int8_t *o=BI3b+((size_t)g*(m/2)+p)*128;
uint64_t parity0=0,parity1=0;
for (int t=0;t<64;t+=32){
__m256i va=sblk_cmb(_mm256_load_si256((const __m256i*)(a+t)), d+t, b->op);
__m256i vb=sblk_cmb(_mm256_load_si256((const __m256i*)(c+t)), e+t, b->op);
__m256i lo=_mm256_unpacklo_epi8(va,vb), hi=_mm256_unpackhi_epi8(va,vb);
_mm256_store_si256((__m256i*)(o+2*t), _mm256_permute2x128_si256(lo,hi,0x20));
_mm256_store_si256((__m256i*)(o+2*t+32), _mm256_permute2x128_si256(lo,hi,0x31));
uint32_t m0=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(_mm256_and_si256(va,one),7));
uint32_t m1=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(_mm256_and_si256(vb,one),7));
parity0 |= (uint64_t)m0<<t; parity1 |= (uint64_t)m1<<t;
}
w0p[g]=parity0; w1p[g]=parity1;
}
}
}
static const int SAB_T[7][5] = {
{0,0, 1,1, 1}, /* A11+A22 */
{1,0, 1,1, 1}, /* A21+A22 */
{0,0, 0,0, 0}, /* A11 */
{1,1, 1,1, 0}, /* A22 */
{0,0, 0,1, 1}, /* A11+A12 */
{1,0, 0,0, 2}, /* A21-A11 */
{0,1, 1,1, 2}, /* A12-A22 */
};
static const int SBB_T[7][5] = {
{0,0, 1,1, 1}, /* B11+B22 */
{0,0, 0,0, 0}, /* B11 */
{0,1, 1,1, 2}, /* B12-B22 */
{1,0, 0,0, 2}, /* B21-B11 */
{1,1, 1,1, 0}, /* B22 */
{0,0, 0,1, 1}, /* B11+B12 */
{1,0, 1,1, 1}, /* B21+B22 */
};
static int8_t M2[7][256*256+384] __attribute__((aligned(64)));
/* Cblk += sum_k s_k * M_k (up to 4 terms; NULL skips; dst stride n, block m x m) */
AVX2 static void cmb(int8_t*C,int n,int r0,int c0,int m,
const int8_t*M0,int s0,const int8_t*M1,int s1,
const int8_t*M2,int s2,const int8_t*M3,int s3){
for (int i=0;i<m;i++){
int8_t *d=C+(size_t)(r0+i)*n+c0;
const int8_t *p0=M0?M0+(size_t)i*m:0, *p1=M1?M1+(size_t)i*m:0;
const int8_t *p2=M2?M2+(size_t)i*m:0, *p3=M3?M3+(size_t)i*m:0;
for (int j=0;j<m;j+=32){
__m256i v;
if (p0){ __m256i t=_mm256_load_si256((const __m256i*)(p0+j));
v = s0>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t); }
else v = _mm256_setzero_si256();
if (p1){ __m256i t=_mm256_load_si256((const __m256i*)(p1+j));
v = _mm256_add_epi8(v, s1>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t)); }
if (p2){ __m256i t=_mm256_load_si256((const __m256i*)(p2+j));
v = _mm256_add_epi8(v, s2>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t)); }
if (p3){ __m256i t=_mm256_load_si256((const __m256i*)(p3+j));
v = _mm256_add_epi8(v, s3>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t)); }
_mm256_store_si256((__m256i*)(d+j), v);
}
}
}
/* [sm2_nt] streaming-store clone of cmb, for the final C write only. */
AVX2 static void cmb_nt(int8_t*C,int n,int r0,int c0,int m,
const int8_t*M0,int s0,const int8_t*M1,int s1,
const int8_t*M2,int s2,const int8_t*M3,int s3){
for (int i=0;i<m;i++){
int8_t *d=C+(size_t)(r0+i)*n+c0;
const int8_t *p0=M0?M0+(size_t)i*m:0, *p1=M1?M1+(size_t)i*m:0;
const int8_t *p2=M2?M2+(size_t)i*m:0, *p3=M3?M3+(size_t)i*m:0;
for (int j=0;j<m;j+=32){
__m256i v;
if (p0){ __m256i t=_mm256_load_si256((const __m256i*)(p0+j));
v = s0>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t); }
else v = _mm256_setzero_si256();
if (p1){ __m256i t=_mm256_load_si256((const __m256i*)(p1+j));
v = _mm256_add_epi8(v, s1>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t)); }
if (p2){ __m256i t=_mm256_load_si256((const __m256i*)(p2+j));
v = _mm256_add_epi8(v, s2>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t)); }
if (p3){ __m256i t=_mm256_load_si256((const __m256i*)(p3+j));
v = _mm256_add_epi8(v, s3>0?t:_mm256_sub_epi8(_mm256_setzero_si256(),t)); }
_mm256_stream_si256((__m256i*)(d+j), v);
}
}
}
/* [mc2g_] merged four-bit correction: all four kb tables live at once.
Row loop runs ONCE over the 256 rows; the four accumulators are XORed in a
register and stored once (no memset, no read-modify-write). XOR is
associative+commutative and Pbb is fully rewritten (n*W*8 bytes) every call,
so every byte is identical to the original three-step form. */
static __m256i XT4_256[4][16][16] __attribute__((aligned(64)));
alignas(64) static __m256i TB[656];
AVX2 static void corr4_256(){
constexpr int n=256,W=4,SLOTS=656;
for(int wb=0;wb<4;wb++){
for(int g=0;g<11;g++){
__m256i*T=(__m256i*)TB+g*64;
const __m256i r0=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+g*6+0)*W));
const __m256i r1=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+g*6+1)*W));
const __m256i r2=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+g*6+2)*W));
const __m256i r3=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+g*6+3)*W));
__m256i v=_mm256_setzero_si256();
_mm256_store_si256(T,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+1,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+3,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+2,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+6,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+7,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+5,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+4,v);
v=_mm256_xor_si256(v,r3); _mm256_store_si256(T+12,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+13,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+15,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+14,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+10,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+11,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+9,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+8,v);
if(g<10){
const __m256i r4=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+g*6+4)*W));
const __m256i r5=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+g*6+5)*W));
v=_mm256_xor_si256(v,r4); _mm256_store_si256(T+24,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+25,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+27,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+26,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+30,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+31,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+29,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+28,v);
v=_mm256_xor_si256(v,r3); _mm256_store_si256(T+20,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+21,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+23,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+22,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+18,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+19,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+17,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+16,v);
v=_mm256_xor_si256(v,r5); _mm256_store_si256(T+48,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+49,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+51,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+50,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+54,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+55,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+53,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+52,v);
v=_mm256_xor_si256(v,r3); _mm256_store_si256(T+60,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+61,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+63,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+62,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+58,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+59,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+57,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+56,v);
v=_mm256_xor_si256(v,r4); _mm256_store_si256(T+40,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+41,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+43,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+42,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+46,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+47,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+45,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+44,v);
v=_mm256_xor_si256(v,r3); _mm256_store_si256(T+36,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+37,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+39,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+38,v);
v=_mm256_xor_si256(v,r2); _mm256_store_si256(T+34,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+35,v);
v=_mm256_xor_si256(v,r1); _mm256_store_si256(T+33,v);
v=_mm256_xor_si256(v,r0); _mm256_store_si256(T+32,v);
}
}
for(int i=0;i<n;i+=2){
const uint64_t w0=A1b[(size_t)i*W+wb],w1=A1b[(size_t)(i+1)*W+wb];
__m256i x0=wb?_mm256_load_si256((const __m256i*)(Pbb+(size_t)i*W)):_mm256_setzero_si256();
__m256i y0=wb?_mm256_load_si256((const __m256i*)(Pbb+(size_t)(i+1)*W)):_mm256_setzero_si256();
__m256i x1=_mm256_setzero_si256(),y1=x1;
#pragma GCC unroll 5
for(int g=0;g<10;g+=2){
const __m256i*t0=(const __m256i*)TB+g*64,*t1=t0+64;
x0=_mm256_xor_si256(x0,_mm256_load_si256(t0+((w0>>(6*g))&63)));
y0=_mm256_xor_si256(y0,_mm256_load_si256(t0+((w1>>(6*g))&63)));
x1=_mm256_xor_si256(x1,_mm256_load_si256(t1+((w0>>(6*(g+1)))&63)));
y1=_mm256_xor_si256(y1,_mm256_load_si256(t1+((w1>>(6*(g+1)))&63)));
}
x0=_mm256_xor_si256(x0,_mm256_load_si256((const __m256i*)TB+640+(w0>>60)));
y0=_mm256_xor_si256(y0,_mm256_load_si256((const __m256i*)TB+640+(w1>>60)));
_mm256_store_si256((__m256i*)(Pbb+(size_t)i*W),_mm256_xor_si256(x0,x1));
_mm256_store_si256((__m256i*)(Pbb+(size_t)(i+1)*W),_mm256_xor_si256(y0,y1));
}
}
}
static inline int sgn_of(int op){ return op==0 ? 0 : (op==1 ? 1 : -1); }
AVX2 static void buildA4(int h, int n, const int8_t*A, int k1, int k2){
const __m256i msk=_mm256_set1_epi8(0x7F);
const int o1 = 512;
const int *t1=SAB_T[k1], *t2=SAB_T[k2];
const int s1=sgn_of(t1[4]), s2=sgn_of(t2[4]), s12=s1*s2;
const size_t R1=(size_t)(t1[0]*o1+t2[0]*h)*n + t1[1]*o1+t2[1]*h;
const size_t R2=(size_t)(t1[2]*o1+t2[0]*h)*n + t1[3]*o1+t2[1]*h;
const size_t R3=(size_t)(t1[0]*o1+t2[2]*h)*n + t1[1]*o1+t2[3]*h;
const size_t R4=(size_t)(t1[2]*o1+t2[2]*h)*n + t1[3]*o1+t2[3]*h;
int W=h/64;
for (int i=0;i<h;i++){
const int8_t*p1=A+R1+(size_t)i*n;
const int8_t*p2=A+R2+(size_t)i*n;
const int8_t*p3=A+R3+(size_t)i*n;
const int8_t*p4=A+R4+(size_t)i*n;
int8_t*dst=A0b+2*(size_t)i*h; /* [c1w_bd] duplicated pairs */
uint64_t*w1=A1b+(size_t)i*W;
for (int k=0;k<h;k+=64){
_mm_prefetch((const char*)(p1+k+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(p2+k+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(p3+k+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(p4+k+4096), _MM_HINT_T0);
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(dst+2*k+256)));
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(dst+2*k+320)));
__m256i v0=_mm256_load_si256((const __m256i*)(p1+k));
__m256i v1=_mm256_load_si256((const __m256i*)(p1+k+32));
if(s1){ __m256i u0=_mm256_load_si256((const __m256i*)(p2+k));
__m256i u1=_mm256_load_si256((const __m256i*)(p2+k+32));
if(s1>0){v0=_mm256_add_epi8(v0,u0);v1=_mm256_add_epi8(v1,u1);}
else {v0=_mm256_sub_epi8(v0,u0);v1=_mm256_sub_epi8(v1,u1);} }
if(s2){ __m256i u0=_mm256_load_si256((const __m256i*)(p3+k));
__m256i u1=_mm256_load_si256((const __m256i*)(p3+k+32));
if(s2>0){v0=_mm256_add_epi8(v0,u0);v1=_mm256_add_epi8(v1,u1);}
else {v0=_mm256_sub_epi8(v0,u0);v1=_mm256_sub_epi8(v1,u1);} }
if(s12){__m256i u0=_mm256_load_si256((const __m256i*)(p4+k));
__m256i u1=_mm256_load_si256((const __m256i*)(p4+k+32));
if(s12>0){v0=_mm256_add_epi8(v0,u0);v1=_mm256_add_epi8(v1,u1);}
else {v0=_mm256_sub_epi8(v0,u0);v1=_mm256_sub_epi8(v1,u1);} }
{ /* [c1w_bd] each 16-bit A pair stored twice => vpbroadcastd is exact.
vpunpcklwd is per-128-bit-lane, so the two halves are put back in
order with vperm2i128: pair P then lives at byte 4*P exactly. */
__m256i e0=_mm256_and_si256(v0,msk), e1=_mm256_and_si256(v1,msk);
__m256i l0=_mm256_unpacklo_epi16(e0,e0), h0=_mm256_unpackhi_epi16(e0,e0);
__m256i l1=_mm256_unpacklo_epi16(e1,e1), h1=_mm256_unpackhi_epi16(e1,e1);
_mm256_store_si256((__m256i*)(dst+2*k), _mm256_permute2x128_si256(l0,h0,0x20));
_mm256_store_si256((__m256i*)(dst+2*k+32), _mm256_permute2x128_si256(l0,h0,0x31));
_mm256_store_si256((__m256i*)(dst+2*k+64), _mm256_permute2x128_si256(l1,h1,0x20));
_mm256_store_si256((__m256i*)(dst+2*k+96), _mm256_permute2x128_si256(l1,h1,0x31)); }
w1[k>>6]=(uint64_t)(uint32_t)_mm256_movemask_epi8(v0)
|((uint64_t)(uint32_t)_mm256_movemask_epi8(v1)<<32);
}
}
}
/* [nm1_4t] same four-term composition on the B side; identical layout to bld_B
(BI3b interleaved pairs + the padded BP2 parity panel with its zeroed sentinel row). */
AVX2 static void buildB4(int h, int n, const int8_t*B, int k1, int k2){
const __m256i one=_mm256_set1_epi8(1);
const int o1 = 512;
const int *t1=SBB_T[k1], *t2=SBB_T[k2];
const int s1=sgn_of(t1[4]), s2=sgn_of(t2[4]), s12=s1*s2;
const size_t R1=(size_t)(t1[0]*o1+t2[0]*h)*n + t1[1]*o1+t2[1]*h;
const size_t R2=(size_t)(t1[2]*o1+t2[0]*h)*n + t1[3]*o1+t2[1]*h;
const size_t R3=(size_t)(t1[0]*o1+t2[2]*h)*n + t1[1]*o1+t2[3]*h;
const size_t R4=(size_t)(t1[2]*o1+t2[2]*h)*n + t1[3]*o1+t2[3]*h;
int ng=h/64, W=h/64;
for (int p=0;p<h/2;p++){
const int8_t *a1=B+R1+(size_t)(2*p)*n, *b1=B+R1+(size_t)(2*p+1)*n;
const int8_t *a2=B+R2+(size_t)(2*p)*n, *b2=B+R2+(size_t)(2*p+1)*n;
const int8_t *a3=B+R3+(size_t)(2*p)*n, *b3=B+R3+(size_t)(2*p+1)*n;
const int8_t *a4=B+R4+(size_t)(2*p)*n, *b4=B+R4+(size_t)(2*p+1)*n;
uint64_t *w0p=BP2+(size_t)((2*p)>>6)*(size_t)(65*W)+(size_t)((2*p)&63)*W;
uint64_t *w1p=BP2+(size_t)((2*p+1)>>6)*(size_t)(65*W)+(size_t)((2*p+1)&63)*W;
for (int g=0;g<ng;g++){
const int8_t *a=a1+64*g, *c=b1+64*g;
const int8_t *d=a2+64*g, *e=b2+64*g;
const int8_t *f=a3+64*g, *h3=b3+64*g;
const int8_t *m4=a4+64*g, *n4=b4+64*g;
int8_t *o=BI3b+((size_t)g*(h/2)+p)*128;
uint64_t parity0=0,parity1=0;
_mm_prefetch((const char*)(a +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(c +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(d +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(e +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(f +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(h3+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(m4+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(n4+4096), _MM_HINT_T0);
for (int t=0;t<64;t+=32){
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(o+2*t+256)));
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(o+2*t+320)));
__m256i va=_mm256_load_si256((const __m256i*)(a+t));
__m256i vb=_mm256_load_si256((const __m256i*)(c+t));
if(s1){ __m256i u0=_mm256_load_si256((const __m256i*)(d+t));
__m256i u1=_mm256_load_si256((const __m256i*)(e+t));
if(s1>0){va=_mm256_add_epi8(va,u0);vb=_mm256_add_epi8(vb,u1);}
else {va=_mm256_sub_epi8(va,u0);vb=_mm256_sub_epi8(vb,u1);} }
if(s2){ __m256i u0=_mm256_load_si256((const __m256i*)(f+t));
__m256i u1=_mm256_load_si256((const __m256i*)(h3+t));
if(s2>0){va=_mm256_add_epi8(va,u0);vb=_mm256_add_epi8(vb,u1);}
else {va=_mm256_sub_epi8(va,u0);vb=_mm256_sub_epi8(vb,u1);} }
if(s12){__m256i u0=_mm256_load_si256((const __m256i*)(m4+t));
__m256i u1=_mm256_load_si256((const __m256i*)(n4+t));
if(s12>0){va=_mm256_add_epi8(va,u0);vb=_mm256_add_epi8(vb,u1);}
else {va=_mm256_sub_epi8(va,u0);vb=_mm256_sub_epi8(vb,u1);} }
__m256i lo=_mm256_unpacklo_epi8(va,vb), hi=_mm256_unpackhi_epi8(va,vb);
_mm256_store_si256((__m256i*)(o+2*t), _mm256_permute2x128_si256(lo,hi,0x20));
_mm256_store_si256((__m256i*)(o+2*t+32), _mm256_permute2x128_si256(lo,hi,0x31));
uint32_t m0=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(_mm256_and_si256(va,one),7));
uint32_t m1=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(_mm256_and_si256(vb,one),7));
parity0 |= (uint64_t)m0<<t; parity1 |= (uint64_t)m1<<t;
}
w0p[g]=parity0; w1p[g]=parity1;
}
}
}
/* one 256 leaf straight out of the ORIGINAL operands (k1 = outer case, k2 = inner). */
AVX2 static void submm512_4(int m, int n, const int8_t*A, const int8_t*B, int k1, int k2,
int8_t*Cp){
buildA4(m,n,A,k1,k2); buildB4(m,n,B,k1,k2);
{ const int W=m>>6; for (int g=0; g<(m>>6); g++) memset(BP2+(size_t)g*65*W+64*W, 0, (size_t)W*8); }
if(m==256) corr4_256(); else corr_rus(m);
for (int g=0; g<m; g+=64)
for (int i=0; i<m; i+=2){
mk_s(m, A0b+2*(size_t)i*m, BI3b+(size_t)(g>>6)*(m/2)*128, Cp, PBb, g);
epilogue(m, Cp, (const uint64_t*)Pbb, g, i, 2);
}
}
/* [nm1_4t] the 7 inner leaves are built straight from the ORIGINAL A/B (k1 = this
outer case). No RawA/RawB, no materialize512. */
/* [mc1y_] TOP-LEVEL FOLD. C11=M0+M3-M4+M6 C12=M2+M4 C21=M1+M3 C22=M0+M2-M1+M5
(CLRS P-form). Block k1's four quadrants go straight into the destination
C blocks instead of MMS[k1]; M (0=store, 1=add, -1=sub) is the block's role in
the P-form combination, so no MMS buffer and no final cmb_nt pass exist. */
struct KD { int R, Cc, M; };
static constexpr KD KDEST[7][2] = {
{{ 0, 0, 0}, {512, 512, 0}}, /* M0: +C11 +C22 (first writer of C11,C22) */
{{512, 0, 0}, {512, 512, -1}}, /* M1: +C21 -C22 */
{{ 0, 512, 0}, {512, 512, 1}}, /* M2: +C12 +C22 */
{{ 0, 0, 1}, {512, 0, 1}}, /* M3: +C11 +C21 */
{{ 0, 0,-1}, { 0, 512, 1}}, /* M4: -C11 +C12 */
{{512, 512, 1}, { 0, 0, 0}}, /* M5: +C22 */
{{ 0, 0, 1}, { 0, 0, 0}}, /* M6: +C11 */
};
static constexpr int KDN[7] = {2,2,2,2,2,1,1};
#define LD32(p) _mm256_load_si256((const __m256i*)(p))
#define ST32(p,v) _mm256_store_si256((__m256i*)(p), (v))
#define ADD8(a,b) _mm256_add_epi8((a),(b))
#define SUB8(a,b) _mm256_sub_epi8((a),(b))
template<int MODE> static inline void upd32(int8_t*p, __m256i v){
if(MODE==0) ST32(p, v);
else if(MODE==1) ST32(p, ADD8(LD32(p), v));
else ST32(p, SUB8(LD32(p), v));
}
/* [mc1z_] FUSED quadrant pass. The four destination quadrants share the same
source row, so each M2 block is read once per k1 instead of 12/7 times.
v00..v11 are the P-form sums (identical expressions to the four qw calls). */
template<int MODE1,int MODE2>
AVX2 static void qwall(int8_t*C,int n,int m,int8_t*D2){
const int h=m>>1;
int8_t *a1=C,*a2=C+h,*a3=C+(size_t)h*n,*a4=C+(size_t)h*n+h;
int8_t *b1=0,*b2=0,*b3=0,*b4=0;
if(D2){ b1=D2; b2=D2+h; b3=D2+(size_t)h*n; b4=D2+(size_t)h*n+h; }
for(int i=0;i<h;i++){
const int8_t *p0=M2[0]+(size_t)i*h,*p1=M2[1]+(size_t)i*h,*p2=M2[2]+(size_t)i*h,
*p3=M2[3]+(size_t)i*h,*p4=M2[4]+(size_t)i*h,*p5=M2[5]+(size_t)i*h,
*p6=M2[6]+(size_t)i*h;
int8_t *e1=a1+(size_t)i*n,*e2=a2+(size_t)i*n,*e3=a3+(size_t)i*n,*e4=a4+(size_t)i*n;
int8_t *f1=b1?b1+(size_t)i*n:0,*f2=b2?b2+(size_t)i*n:0,
*f3=b3?b3+(size_t)i*n:0,*f4=b4?b4+(size_t)i*n:0;
if(i+2<h){
const char*q0=(const char*)(M2[0]+(size_t)(i+2)*h);
_mm_prefetch(q0,_MM_HINT_T0); _mm_prefetch(q0+64,_MM_HINT_T0); _mm_prefetch(q0+128,_MM_HINT_T0); _mm_prefetch(q0+192,_MM_HINT_T0);
_mm_prefetch((const char*)(M2[1]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[2]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[3]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[4]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[5]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[6]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(e1+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(e2+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(e3+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(e4+(size_t)2*n),_MM_HINT_T0);
if(f1){ _mm_prefetch((const char*)(f1+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(f2+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(f3+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(f4+(size_t)2*n),_MM_HINT_T0); }
}
for(int j=0;j<h;j+=32){
const __m256i v0=LD32(p0+j),v1=LD32(p1+j),v2=LD32(p2+j),v3=LD32(p3+j),
v4=LD32(p4+j),v5=LD32(p5+j),v6=LD32(p6+j);
const __m256i y00=ADD8(ADD8(v0,v3),SUB8(v6,v4));
const __m256i y01=ADD8(v2,v4);
const __m256i y10=ADD8(v1,v3);
const __m256i y11=ADD8(ADD8(SUB8(v0,v1),v2),v5);
upd32<MODE1>(e1+j,y00); upd32<MODE1>(e2+j,y01);
upd32<MODE1>(e3+j,y10); upd32<MODE1>(e4+j,y11);
if(f1){ upd32<MODE2>(f1+j,y00); upd32<MODE2>(f2+j,y01);
upd32<MODE2>(f3+j,y10); upd32<MODE2>(f4+j,y11); }
}
}
}
template<int K1>
AVX2 static void strassen2_512x(const int8_t*A, const int8_t*B, int n, int8_t*C){
const int h=256;
for(int k=0;k<7;k++) submm512_4(h,n,A,B,K1,k,M2[k]);
int8_t *D1=C+(size_t)KDEST[K1][0].R*n+KDEST[K1][0].Cc;
int8_t *D2=(KDN[K1]>1)?C+(size_t)KDEST[K1][1].R*n+KDEST[K1][1].Cc:0;
if(KDN[K1]>1) qwall<KDEST[K1][0].M,KDEST[K1][1].M>(D1,n,512,D2);
else qwall<KDEST[K1][0].M,0>(D1,n,512,0);
}
AVX2 static void strassen1(int n, const int8_t*A, const int8_t*B, int8_t*C){
strassen2_512x<0>(A,B,n,C); strassen2_512x<1>(A,B,n,C);
strassen2_512x<2>(A,B,n,C); strassen2_512x<3>(A,B,n,C);
strassen2_512x<4>(A,B,n,C); strassen2_512x<5>(A,B,n,C);
strassen2_512x<6>(A,B,n,C);
_mm_sfence();
}
AVX2 void matrix_multiply(int n, const int8_t *A, const int8_t *B, int8_t *C){
if (n<=0 || (n&255) || n>NMAX) { generic_mm(n,A,B,C); return; }
if (n==1024) { strassen1(n,A,B,C); return; }
build_A0_A1(n,A);
if(n==1024){
build_BI9(n,B); build_BI3_tail(n,B); build_BP2(n,B); corr_rus(n);
for (int g=0; g<10; g++)
for (int i=0; i<n; i+=2){
mk_s6(n, A0b+(size_t)i*n, BI9b+(size_t)g*(n/2)*192, (uint8_t*)mscr);
epilogue(n, C, (const uint64_t*)Pbb, g*96, i, 3);
}
for (int i=0; i<n; i+=2){
mk_s(n, A0b+2*(size_t)i*n, BI3b, C, PBb, 0);
epilogue(n, C, (const uint64_t*)Pbb, 960, i, 2);
}
return;
}
build_BI3(n,B); build_Bpbits(n,B);
if(n==4096) corr16(); else corr_pass3(n);
for (int g=0, jbg=0; g<n; g+=64, jbg=g)
for (int i=0; i<n; i+=2)
{ mk_s(n, A0b+2*(size_t)i*n, BI3b+(size_t)(g>>6)*(n/2)*128, C, PBb, g);
epilogue(n, C, (const uint64_t*)Pbb, g, i, 2); }
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 6.722 ms | 1 MB + 712 KB | Accepted | Score: 100 | 显示更多 |