// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104215 <https://duck.ac/submission/104215>
// 用途：**照搬**了该提交 `corr4_256()` 的两处改动 —— 删掉 `memset(Pbb,0,(size_t)n*W*8)`、
// 并把第 0 个 k 字的那一趟（`kb == 0`）从「读-改-写」改成**直接写**
// （他头部自述 "In corr4_256, write the first 64-bit correction chunk directly instead of
// zeroing the entire correction matrix and then reading/XORing a known zero"）。
// `tools/fnmd5.py` 全量正文 diff 已核：**#104215 = 本账号 #103964 的正文 + 恰好三处改动**
// （本发这两处 + `mk_s` 底部两条预取），mem 与本账号相同（3452 KB）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104524 <https://duck.ac/submission/104524>
// 用途：本文件正文 = 该提交（我方当前最好件 7.574258 ms，本题榜首）逐字节正文 + 本发唯一改动；
// 其自带继承引用表原样保留在下方（表体一字未动）。
// ======================
// ===== 思路 =====
// **一发单变量**：`corr4_256()` 的行循环里第 0 趟 k 字的写入形态。
//
// 原本 `memset(Pbb, 0, …)` 先把整张修正矩阵清零，4 个 k 字的每一趟都做「读-改-写」；
// 但 `kb == 0` 那一趟读到的必然是刚写下的 0 ⇒ 那一次读 + 一次 xor 是纯浪费。
// 改成 `kb == 0` 直接 `store`、其余趟保持 `load ^ xor store`，**被写入的值逐位不变**
// （`x ^ 0 == x`），只是把 49 次调用 × 4 趟 × 256 行里的第一趟变成纯写。
//
// 为什么值得单独再发一次：本账号 #104035 实测**单独**去 memset 是 **+59.9 µs（更慢）**，
// 而对手把它与 `mk_s` 删预取合在一起时整体是 **−81.5 µs**；本账号随后实测 `mk_s` 那一步
// 单独只值 −71.9（#104524），两个数**差 ~69 µs** ⇒ 这一处与那一步是**强超可加**的。
// 本发就是把这个组合补全（我方基座上带着已判奇的 corr4_256 八链拆分）。
//
// 正确性：只改写入形态、不改任何被写入的值；本地 n=1024 整矩阵 FNV
// `chk=34fd0a86984a4fc3` 与基座逐位相同。
// ================
// ===== 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）逐字节正文 + 本发唯一改动；
// 其自带继承引用表原样保留在下方（表体一字未动）。
// ======================
// ===== 思路 =====
// **一发单变量**：`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,
};
AVX2 static void epilogue(int n, int8_t*C, const uint64_t*P, int jb, int i, int NCH){
const __m128i bitm=_mm_setr_epi8(1,2,4,8,16,32,64,(char)0x80,1,2,4,8,16,32,64,(char)0x80);
const __m128i hi=_mm_set1_epi8((char)0x80), zero=_mm_setzero_si128();
const size_t bo=(size_t)jb>>3;
const uint8_t*base=(const uint8_t*)P+(size_t)i*(n>>3);
for(int r=0;r<2;r++){
const uint8_t*qr=base+(size_t)r*(n>>3);
for(int t=0;t<NCH;t++){
size_t bi_=bo+4*(size_t)t;
/* bi_ = bo + 4t with bo = jb>>3 = g*12 and jb = g*96 (tail: jb=960,
bo=120). 12g and 4t are both multiples of 4 => (bi_&7) in {0,4}
=> sh = (bi_&7)*8 in {0,32} => the four bytes w>>(sh+8k) are EXACTLY
the four bytes at qr[bi_..bi_+3]. Read them directly: this removes the
64-bit load, the runtime shift count and the (bi_&~7) fixup. */
const uint8_t*wb=qr+bi_;
__m256i v=_mm256_load_si256((const __m256i*)(mscr+(r*NCH+t)*32));
/* one 8-byte table entry per byte of Pbb, index scaled by 8 => a single
scale-8 addressing mode per lookup (no shift, no *16 fixup). */
__m128i z0=_mm_loadl_epi64((const __m128i*)(SHB8+(size_t)wb[0]*8));
__m128i z1=_mm_loadl_epi64((const __m128i*)(SHB8+(size_t)wb[1]*8));
__m128i z2=_mm_loadl_epi64((const __m128i*)(SHB8+(size_t)wb[2]*8));
__m128i z3=_mm_loadl_epi64((const __m128i*)(SHB8+(size_t)wb[3]*8));
__m256i msk=_mm256_set_m128i(_mm_unpacklo_epi64(z2,z3),_mm_unpacklo_epi64(z0,z1));
int8_t*c=C+(size_t)(i+r)*n+jb+t*32;
_mm256_storeu_si256((__m256i*)c,_mm256_xor_si256(v,msk));
}
}
}
/* ---- 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"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\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, %%ymm13, %%ymm13\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, %%ymm14, %%ymm14\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, %%ymm15, %%ymm15\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"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\n\t"
"vmovdqa 160(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm13, %%ymm13\n\t"
"vmovdqa 192(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\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, %%ymm11, %%ymm11\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm15, %%ymm15\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"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\n\t"
"vmovdqa 288(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm13, %%ymm13\n\t"
"vmovdqa 320(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm14, %%ymm14\n\t"
"vmovdqa 352(%[bi]), %%ymm4\n\t"
"vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm15, %%ymm15\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"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw %%ymm4, %%ymm1, %%ymm5\n\t"
"vpaddw %%ymm5, %%ymm12, %%ymm12\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, %%ymm13, %%ymm13\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, %%ymm14, %%ymm14\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, %%ymm15, %%ymm15\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);
}
}
}
/* 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 a0=XT4_256[0][(int)((bits>>0)&15)];
__m256i a1=XT4_256[1][(int)((bits>>4)&15)];
__m256i a2=XT4_256[2][(int)((bits>>8)&15)];
__m256i a3=XT4_256[3][(int)((bits>>12)&15)];
__m256i a4=XT4_256[4][(int)((bits>>16)&15)];
__m256i a5=XT4_256[5][(int)((bits>>20)&15)];
__m256i a6=XT4_256[6][(int)((bits>>24)&15)];
__m256i a7=XT4_256[7][(int)((bits>>28)&15)];
a0=_mm256_xor_si256(a0,XT4_256[8][(int)((bits>>32)&15)]);
a1=_mm256_xor_si256(a1,XT4_256[9][(int)((bits>>36)&15)]);
a2=_mm256_xor_si256(a2,XT4_256[10][(int)((bits>>40)&15)]);
a3=_mm256_xor_si256(a3,XT4_256[11][(int)((bits>>44)&15)]);
a4=_mm256_xor_si256(a4,XT4_256[12][(int)((bits>>48)&15)]);
a5=_mm256_xor_si256(a5,XT4_256[13][(int)((bits>>52)&15)]);
a6=_mm256_xor_si256(a6,XT4_256[14][(int)((bits>>56)&15)]);
a7=_mm256_xor_si256(a7,XT4_256[15][(int)((bits>>60)&15)]);
__m256i acc=_mm256_xor_si256(_mm256_xor_si256(_mm256_xor_si256(a0,a1),_mm256_xor_si256(a2,a3)),_mm256_xor_si256(_mm256_xor_si256(a4,a5),_mm256_xor_si256(a6,a7)));
__m256i*dst=(__m256i*)(Pbb+(size_t)i*W);
/* [nn8] kb == 0 reads back what the memset just wrote: x^0 == x. Write it, don't RMW. */
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(C,n, 0,0, m, MMS[0],1, MMS[3],1, MMS[4],-1, MMS[6],1);
cmb(C,n, 0,o, m, MMS[2],1, MMS[4],1, 0,0, 0,0);
cmb(C,n, o,0, m, MMS[1],1, MMS[3],1, 0,0, 0,0);
cmb(C,n, o,o, m, MMS[0],1, MMS[2],1, MMS[1],-1, MMS[5],1);
}
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.585 ms | 3 MB + 380 KB | Accepted | Score: 100 | 显示更多 |