// References: saffah_cc_v41_agg1, https://duck.ac/submission/105098 ; full engine copied. saffah_codex_6s_agg2, https://duck.ac/submission/104215 ; restored code for named components. No license notice appeared with either submission. Inherited citations follow.
// Approach: Keep rival engine and restored correction routine; schedule the two row products with separate YMM temporaries before adds (dual_temp_revadd) in the mk_s microkernel.
// Purpose: Experimental judged submission to test whether dual temporary registers in the mk_s AVX2 microkernel meet the live strict ranking threshold; official Try It improved 12.029930 ms to 11.867511 ms with identical output hash.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104215 <https://duck.ac/submission/104215>
// 用途：**照搬**了该提交相对本账号 #103964 的 `mk_s` 预取改动 —— 删掉内层循环底部的
// `prefetcht0 600(%[bi])` 与 `prefetcht0 1200(%[bi])` 两条软件预取（他头部自述
// "Remove the 600-byte and 1200-byte software prefetches from the 2-row 4-column mk_s
// loop ... the 1536-byte software prefetch is retained"）。
// `tools/fnmd5.py` 全量正文 diff 已核：**#104215 = 本账号 #103964 的正文 + 恰好三处改动**
// （这 2 条预取 + `corr4_256` 去 memset + `kb==0` 直存），且 #104215 的 mem 与本账号
// 完全相同（3452 KB）。
// 定价：本账号 #104035（= #103964 + 只去 memset）实测 **+59.9 µs**（7.719511 vs 7.659603），
// 而他的 7.578116 = #104035 + 这 2 条预取的删除 ⇒ **该删除 = −141.4 µs**。
// 本发**不**搬他的 memset 改动（本账号基座上它是 +59.9）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104447 <https://duck.ac/submission/104447>
// 用途：本文件正文 = 该提交（我方当前最好件 7.646187 ms）逐字节正文 + 本发唯一改动；
// 其自带继承引用表原样保留在下方（表体一字未动）。
// ======================
// ===== 思路（本发）=====
// **一发单变量**：把 `strassen1` 的**最终 C 写出**从普通存储改成**非临时存储**（`vmovntdq`）。
//
// 依据：`strassen1` 是调用者 `C`（1024×1024 = **1 MiB**）的**唯一写入者**（4 个象限各一次
// `cmb`），而 `C` 在计时区内**从不被读回**（`M2`/`MMS` 那些会被 `cmb` 读的中转块不在此列，
// 本发**不**动它们的调用点）。普通存储要对这 1 MiB 先付**整片 RFO 读**（read-for-ownership）；
// `vmovntdq` 直接写，省掉那次读 ⇒ 少 **1 MiB** 的访存流量。
// 实现：把 `cmb` 克隆成 `cmb_nt`（唯一差别是 `_mm256_store_si256` → `_mm256_stream_si256`），
// 只把 `strassen1` 里 4 处 `cmb(C,n,…)` 指向它，并在 `strassen1` 末尾加 `_mm_sfence()`
// 保证返回前 NT 写全局可见（`C` 只在 `n==1024` 分支走 `strassen1`）。
// 出处：做法移植自本队协调者转述的 `mmml1k` 判题机实测（该题 8 MiB 输出 ⇒ −0.268%），
// **本题输出只有 1 MiB ⇒ 按比例预期远小**；本发即为**在本题基座上的正式定价**。
// 正确性：本地 n=1024 全 1 048 576 元素与基座**逐字节相同**。
// ======================
// ===== 思路 =====
// **一发单变量**：`epilogue` 的位展开掩码构造。
//
// 机理：`epilogue` 每个 32 列分块要由 **4 个 Pbb 字节**（= 32 个修正位）造出 32 字节的
// 0x80/0x00 掩码。现行实现是 4 条 `movzx`（取索引）+ 4 条 `vmovq`（SHB8 的 8 字节查表）
// + 2 条 `vpunpcklqdq` + 1 条 `vinserti128` = **11 条 uop（其中 8 条占载入端口）**，
// 而该函数在本题被调用 **25 088 次**、每次 4 个分块 ⇒ 仅"造掩码"就 **1.10 M uop**。
// 本发把它换成 **4 条向量 uop**：
// `vpbroadcastd`（4 个 Pbb 字节 → 每 128 位道内周期 4 复制）
// → `vpshufb`（道内索引 `c div 8`，两个道的副本相同 ⇒ 道 1 可取索引 2/3）
// → `vpand EPD_BIT` → `vpcmpeqb EPD_BIT`（置位处得 0xFF）→ `vpand 0x80`。
// 恒等式推导见函数内注释。**零语义变化**：掩码逐位相同，`chk` 与基座一致。
// 预期：判题机 epilogue 段的单价按本文件 round-18 标定的 **0.2 cycle/uop** 计，
// −24 uop/调用（13→7） × 25 088 调用 ≈ −0.60 M uop ≈ **−35 µs**。
// ================
// **一发单变量**：`mk_s` 内层循环的软件预取构型。
//
// 删掉底部两条 `prefetcht0 600(%[bi])` / `prefetcht0 1200(%[bi])`，保留顶部
// `prefetcht0 1536(%[bi])`（与对手 [1] 完全相同的取舍）。
// 机理：`bi` 是**纯线性**流、每次调用只走 16 KB（一个 g 块），而这两条的目标行就落在
// "下一趟 / 下下趟按需取数"的窗口里（`600 = 512+88`、`1200 = 1024+176`）⇒ 硬件顺序预取器
// 已经覆盖，软件预取只占载入端口槽位（该内层循环 96 条指令里 3 条是预取）。
//
// **附带（同一机理的第二处，§2.18.765 逐条代入）**：`buildA4` 里删掉
// `_mm_prefetch(p1+k+32+4096)` —— 题面保证 A 是 4096 字节对齐，`R1 = (…)*1024 + t1[1]*512 +
// t2[1]*256` ⇒ `p1` 是 256 字节对齐，`k` 是 64 的倍数 ⇒ `p1+k` 与 `p1+k+32` **同一条 64 字节
// cache line** ⇒ 这条预取是**纯重复**（只占一个槽位、不产生任何新取数）。
//
// 正确性：两处都只删**预取**指令，不改任何载入/运算/存储；本地 n=1024 整矩阵 FNV
// `chk=34fd0a86984a4fc3` 与基座逐位相同。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1，提交 #103943 <https://duck.ac/submission/103943>
// 用途：本文件正文的来源（= 我方当前最好成绩 7.705298 ms）。该提交自带的完整继承引用表
// 原样保留在下方（只改了段标题字；表体一字未动，其中的 duck.ac 提交地址仍然有效并逐条列明）。
// [2] duck.ac 用户 saffah_cc_v41_agg1，提交 #102502 / #102798（`mmmc2k` 我方件）
// <https://duck.ac/submission/102502> <https://duck.ac/submission/102798>
// 用途：**直接复制**了 `buildA4`/`buildB4` 里的 **13 条 `_mm_prefetch(..., +4096)`**
// （`buildA4` 5 条：p1k/p1k+32/p2/p3/p4；`buildB4` 8 条：a1..a4、b1..b4）。
// 本发把上一发（#103943）移植时**漏掉的这部分原样补回** —— 那 13 条是原作者写在
// 打包器里的，我搬运时只搬了循环体。
// ======================
// ===== 思路 =====
// 本发在三处各动一点（都是**单变量判题机 probe 先定价、符号三轮一致**的）：
//
// (a) **补回上一发漏搬的 13 条预取**（出处 [2]）。
// 上一发的四路打包器要顺着 4 条相隔 256 KB 的流（A 的四个 256×256 块）走，
// 硬件预取器跟不住 ⇒ 原作者在 `buildA4`/`buildB4` 里显式写了 `+4096` 的 `T0` 预取。
// 判题机 probe 三轮交错：**−5.7 / −7.1 / +0.4 µs**（边际、同源，故随本发一起交）。
// (b) 文件级 pragma 加 **`rename-registers`**。probe 三轮交错头对头（对照 / 本项）：
// **−5.9 / −5.8 / −10.5 µs**，三轮同号；gcc 9.3 接受（asm 2586 → 2579 条）。
// (c) `corr4_256` 内层行循环**从第一次查表开始**，去掉 `setzero` + 多出的那一次 xor：
// 每个 (行, k 字) 省 2 条 uop，每次调用 4 个 k 字 × 256 行，共 49 次调用。
// **XOR 可交换/可结合 ⇒ 逐位恒等**。
//
// 正确性：三处都不改任何被加数的集合与加法顺序；(a)(c) 只改访存时机与一条 uop，
// (b) 只改寄存器分配 ⇒ 本地 n=1024 全 1 048 576 元素 `chk=34fd0a86984a4fc3` 与基座逐位相同。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #102793 <https://duck.ac/submission/102793>
// （t=8621.412 us，本题当前榜首）用途：**整份复刻为基座** —— 本文件正文逐字节等于该提交。
// ★ 逐行 diff 已核（本账号本轮实测）：该提交 = **本账号 #101585
// <https://duck.ac/submission/101585> 的正文 + 唯一一处改动**（62 行差异、只有一个内容块）：
// 加 `Gtab`/`corr_rows<NY>`/`corr_rus(n)` 并把 `corr2k(n)` 换成 `corr_rus(n)`。
// [2] 本账号（saffah_cc_v41_agg1）提交 #101585 <https://duck.ac/submission/101585>
// （t=8779.313 us，`work/cy1_i.cpp`）= 本文件去掉 [1] 那处改动后的正文；
// 它本身 = 对手 #100989 那代引擎 + 本账号五把已验证的刀（O3 pragma / 融合建面板 /
// BI9 循环序 / epilogue 去变长移位 / 位展开查表）。
// [3] 本账号 mmmc2k 工件 `problems/mmmc2k/work/zm2_c_final.cpp` / `.cache/subcode/102110.txt`
// —— **`Gtab`/`corr_rows`/`corr_rus` 的原始出处**（本账号自己的 4-Russians 字节表修正项）。
// [4] /home/yjp/duck.ac/BRIEF.md §2.18.609（同族跨题移植）、§2.18.499①（`mine > T` 时先整份
// 复刻）、§2.18.555③（自家标定不可跨题迁移）。
// ★★★ 本发（#102941 之上**只改一处**）：给本题**加一层 Strassen**（n=1024 → 7 × 512）。
// 本题原本**没有 Strassen**（`grep -c strassen1` = 0），而同族两题都是靠这一步拿下的：
// `mmmc4k` 同步骤 －5.77%、`mmmc2k` 我方实测 －4.03%（#102917）。
// 实现：叶子 m=512 走**本文件原有的通用 64 列流水线**（`mk_s` + `epilogue(...,2)`），
// 面板由**融合打包器** `bld_A`/`bld_B` 直接从父矩阵的两个块边合成边产出
// （**不物化 SA/SB**，`§2.18.590` DIF：避免"写出来又读回去"）。
// ★ `corr_rus` 扩到 **NY=2**（`corr_rows<2>`）以覆盖 512 叶子。★ CLRS P-形组合用 `cmb`（四项带符号，常规存储、无需 sfence）。
//
// 【写码前的净账】内核（`mk_s6`/`mk_s`+epilogue）占 ~90% ⇒ 毛 = 12.5%×90% = **－11.3%**；
// 按 mmmc2k 实测的 36.3% 兑现 ⇒ **－4.1%** ⇒ 8.619 × 0.959 = **8.27 ≤ 8.536**（目标）。
//
// 【闸门】n=1024 **全 1 048 576 元素** `old_vs_new = 0`；**零共享代码的独立素朴预言机** old/new 皆 0；
// **两条阳性对照**：数值类（`SAB_T[3][4]` 0→1）**522 124** 失配、偏移类（`bld_A` 的 `xa` 用 `c1` 代 `c0`）**1 044 503** 失配
// ⇒ 闸门是活的。★ **路径标签**：改动只在 n=1024 分支 ⇒ **只有 n=1024 承重**（其它 n 走原路，对本刀是空测）。
// ★ 排查过程记录：首版闸门报 **523 954/1 048 576 错**，而且差值**恒为 ±128**
// ⇒ 定位到"term1 对、GF(2) 项错"；用 `corr_pass3` 替换叶子 corr 得 **bad=0** ⇒ 锁定在 B 面板：
// `corr_rus` 读的是**带哨兵行的 padded 面板 `BP2`**（组距 65*W），而我的 `bld_B` 写的是 `Bpb`。
// 修正后通过。（→ 可复用判据：**差值恒为±128 = GF(2) 项错；差值任意 = term1 错**。）
// ======================
// Approach: 保留 #103075 的两层 Strassen 与整数微内核，
// 将 256 阶叶子的 GF(2) 校正改成每次 4 个 k 位的查表。
// 表为 16×16×32 字节约 8 KiB，可完整驻留 L1；每个
// 64 位 A 符号字恰好 16 次查表，所有 XOR 项仍完整覆盖。
// Purpose: 实验性在线验证更小 L1 查表能否补齐精确门槛差距。
// ===== 思路 =====
// 局势：`T = 8621.412`（对手 #102793）、真实目标 `0.99*T + 1us = 8536.198`、
// 我方 mine = 8779.313（#101585）⇒ **落后 0.243 ms = 2.77%**。
//
// ★ 本发**只做一件事**：把 `corr2k(n)` 换成 `corr_rus(n)`（我方自己的 4-Russians 字节表修正项）。
// 取证（逐行 diff，本轮实测）：对手 #102793 **就是**我方 #101585 的正文 + 这一处改动
// ⇒ **这一半在我方基座上的价格是实测的**：8.779313 → 8.621412 ⇒ **−0.158 ms = −1.80%**。
// ★ 它只削掉 corr2k 的 `1.80/6.9 = 26%` —— 因为对手的 `corr2k`（哨兵行 + 成对 w）本身已经不错。
// ★ 相位账：`corr2k 6.9% / mk_s6+epilogue ~90% / 面板构建 ~3%`。
//
// ⚠ **本发单独不足以达标**（1.80% < 2.77%）⇒ 落点约 8.62，仍差 ~1.0%。
// 第二发（主力）是**加一层 Strassen**（本题**没有** Strassen，`grep -c strassen1` = 0）：
// 同族 mmmc4k 上同步骤 −5.77%、mmmc2k 上我方实测 −4.03%（#102917）；内核占 ~90%
// ⇒ 毛 −11.1%，按 36.3% 兑现 ⇒ −4.0% ⇒ 8.78 × 0.96 = 8.43 ≤ 8.536 ✅
// （即使只 25% 兑现也只是 −2.8% ⇒ 8.54，卡边缘 ⇒ **两把叠起来才稳**。）
// ⚠ 本题当前**在本口径下也是红的** ⇒ 本发是"在红状态下的提交式探针"，不会把绿变红。
//
// 正确性：正文与站点上 #102793 的副本逐字节相同（本脚本的 md5 断言）；本地对拍
// vs 朴素 O(n^3)（零共享代码预言机）n=256/512/1024 全 bad=0。
// ================
/* 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[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+(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));
_mm256_store_si256((__m256i*)(o+k), _mm256_and_si256(v0,m));
_mm256_store_si256((__m256i*)(o+k+32), _mm256_and_si256(v1,m));
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+n, *bi=BI; long cnt=n>>3;
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"
"vpbroadcastw (%[x0]), %%ymm0\n\t"
"vpbroadcastw (%[x1]), %%ymm1\n\t"
"vmovdqa 0(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vmovdqa 32(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vmovdqa 64(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vmovdqa 96(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastw 2(%[x0]), %%ymm0\n\t"
"vpbroadcastw 2(%[x1]), %%ymm1\n\t"
"vmovdqa 128(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vmovdqa 160(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vmovdqa 192(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vmovdqa 224(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastw 4(%[x0]), %%ymm0\n\t"
"vpbroadcastw 4(%[x1]), %%ymm1\n\t"
"vmovdqa 256(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vmovdqa 288(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vmovdqa 320(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vmovdqa 352(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastw 6(%[x0]), %%ymm0\n\t"
"vpbroadcastw 6(%[x1]), %%ymm1\n\t"
"vmovdqa 384(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vmovdqa 416(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vmovdqa 448(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vmovdqa 480(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"addq $8, %[x0]\n\t"
"addq $8, %[x1]\n\t"
"addq $512, %[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;
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));
if(t==0) { w0p[g]=(uint64_t)m0; w1p[g]=(uint64_t)m1; }
else { w0p[g]|=(uint64_t)m0<<32; w1p[g]|=(uint64_t)m1<<32; }
}
}
}
}
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 MMS[7][512*512] __attribute__((aligned(64)));
static int8_t M2[7][256*256] __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);
}
}
}
/* Four-bit correction for a 256-square leaf; source B bits use padded BP2 rows. */
static __m256i XT4_256[16][16] __attribute__((aligned(64)));
AVX2 static void corr4_256(void){
const int n=256,W=4;
for(int kb=0;kb<n;kb+=64){
const int w0=kb>>6;
for(int g=0;g<16;g++){
XT4_256[g][0]=_mm256_setzero_si256();
for(int v=1;v<16;v++){
const int lb=__builtin_ctz(v),t2=v&(v-1),k=kb+4*g+lb;
const uint64_t*bp=BP2+(size_t)(k>>6)*(65*W)+(size_t)(k&63)*W;
XT4_256[g][v]=_mm256_xor_si256(XT4_256[g][t2],_mm256_load_si256((const __m256i*)bp));
}
}
for(int i=0;i<n;i++){
const uint64_t bits=A1b[(size_t)i*W+w0];
/* [nm1_4tb] start from the first lookup: the zeroing vpxor and one vxor per
(row, k-word) pair are pure waste. XOR is associative => bit-identical. */
__m256i acc=XT4_256[0][(int)(bits&15)];
for(int g=1;g<16;g++) acc=_mm256_xor_si256(acc,XT4_256[g][(int)((bits>>(4*g))&15)]);
__m256i*dst=(__m256i*)(Pbb+(size_t)i*W);
if (kb == 0) _mm256_store_si256(dst, acc);
else _mm256_store_si256(dst,_mm256_xor_si256(_mm256_load_si256(dst),acc));
}
}
}
static inline int sgn_of(int op){ return op==0 ? 0 : (op==1 ? 1 : -1); }
/* [nm1_4t] FOUR-TERM COMPOSITE PACKER (bld_A with two levels of Strassen composed).
Leaf (k1,k2) needs S2 = (A1 op1 A2) op2 (A3 op1 A4) = A1 + s1*A2 + s2*A3 + s1*s2*A4,
where op1 = SAB_T[k1][4] acts at o1 = 512 granularity and op2 = SAB_T[k2][4] at
h = 256. The four blocks are read straight from A -- nothing is materialised. */
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+(size_t)i*h;
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);
__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);} }
_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);
}
}
}
/* [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;
for (int t=0;t<64;t+=32){
_mm_prefetch((const char*)(a +t+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(c +t+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(d +t+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(e +t+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(f +t+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(h3+t+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(m4+t+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(n4+t+4096), _MM_HINT_T0);
__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));
if(t==0) { w0p[g]=(uint64_t)m0; w1p[g]=(uint64_t)m1; }
else { w0p[g]|=(uint64_t)m0<<32; w1p[g]|=(uint64_t)m1<<32; }
}
}
}
}
/* 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+(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. */
AVX2 static void strassen2_512(int8_t*Cp, const int8_t*A, const int8_t*B, int n, int k1){
const int m=512,h=256;
for(int k=0;k<7;k++){
submm512_4(h,n,A,B,k1,k,M2[k]);
}
cmb(Cp,m,0,0,h, M2[0],1,M2[3],1,M2[4],-1,M2[6],1);
cmb(Cp,m,0,h,h, M2[2],1,M2[4],1,0,0,0,0);
cmb(Cp,m,h,0,h, M2[1],1,M2[3],1,0,0,0,0);
cmb(Cp,m,h,h,h, M2[0],1,M2[2],1,M2[1],-1,M2[5],1);
}
AVX2 static void strassen1(int n, const int8_t*A, const int8_t*B, int8_t*C){
const int m = n>>1; /* 512 */
const int o = m;
for (int k=0;k<7;k++){
strassen2_512(MMS[k], A, B, n, k);
}
/* C11=M0+M3-M4+M6 C12=M2+M4 C21=M1+M3 C22=M0+M2-M1+M5 (CLRS P-form) */
cmb_nt(C,n, 0,0, m, MMS[0],1, MMS[3],1, MMS[4],-1, MMS[6],1);
cmb_nt(C,n, 0,o, m, MMS[2],1, MMS[4],1, 0,0, 0,0);
cmb_nt(C,n, o,0, m, MMS[1],1, MMS[3],1, 0,0, 0,0);
cmb_nt(C,n, o,o, m, MMS[0],1, MMS[2],1, MMS[1],-1, MMS[5],1);
_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+(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+(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 | 7.324 ms | 3 MB + 380 KB | Accepted | Score: 100 | 显示更多 |