// ##### r22ra_ 席（duck.ac 优化舰队）· 2026-10-15 · mmmc1k #####
// ===== REFERENCES（本发新增引用段）=====
// [A] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #127159
// <https://duck.ac/submission/127159>（本题现役最快件 / `mine`）。
// —— **本件正文的基座**：正文为该提交的**逐字节副本**，唯一改动 = 在文件顶部、
// 最后一条 `#include` 之后**独立新起两行** pragma（见下「思路」）。该件自带的上游
// 引用链（对 saffah_codex_6a_agg3 各代与其自身历史件的引用与许可说明）**原样保留
// 在下方**，未作删改。
// LICENSE：上述均为 duck.ac 公开提交（站点 FAQ 明示「你的代码将会被公开，所有人都可见」）；
// 无第三方许可证声明、非来自任何外部仓库 ⇒ 按具名引用（账号 + 提交地址 + 沿用内容）使用，
// 遵守其无附加限制的许可。**本件未参考任何已发表论文 / 外部开源仓库。**
// ================================================
// ===== 思路（r22ra_ 席）=====
// 【本发单变量（相对 [A]）：在文件顶部、最后一条 `#include` 之后**独立新起两行**
// `#pragma GCC optimize("ira-loop-pressure")` 与 `#pragma GCC optimize("sched-pressure")`。
// 正文其余一字未改（源码级 diff = 唯一的 `126a127,128` 两行新增）。】
// 依据：舰队在册 律 N2（实测）—— 该 RA pragma 对在 `1004e8` 上两枚合计 **-2.66 ms = 0.38%**，
// 明写「远超任何代码级微刀」；适用形态 = **大函数 + 深内联 + 手写 asm**（本题具此形态：
// `strassen1` 含被内联的 `submm512_4`/`mk_s`/`epilogue`，热环 442 条指令，全引擎 14 处内联 asm）。
// 写法：律 887 —— **每旗标独立成行**（不写成合并串）。
// 器械自检：gcc9 判题同款口径带 `-Wpragmas` 编译，**警告 = 0**（不带引号的 pragma 会被静默忽略）。
// ---------------- 本席在本题先做的「域检查」（用于决定发不发，结果如实披露）----------------
// * `sched-pressure` **单独**编译 ⇒ `.text` 裸段 md5 `e8b770e3...` **与基座逐字节同**
// ⇒ 该旗标在本题**零边际**（在册律 887/882 的判据）。
// * `ira-loop-pressure` **单独** ⇒ `.text` md5 `6e4cf104...`（+480 B）⇒ **真变**。
// * **两枚同发** ⇒ 与 ira 单独 **逐字节同**（同为 `6e4cf104...`）⇒ 本发的全部可动项来自 ira 那一枚。
// ★器械更正（本席今日实测，建议入册）：`.text` 哈希**不可用 `objdump -s`** —— 其输出头行含文件名
// （`base.o:` 与 `anchor.o:`），会把**仪器噪声读成内容差异**（本席据此一度误判锚件「不同码」）。
// 正解：`objcopy -O binary --only-section=.text` 抽裸段字节后再 md5。
// -------------------------------- 四关 --------------------------------
// ① 三维身份核 `.cache/sec_cmp.py` v4：基座 vs 锚件 = **同码成立**；基座 vs 本件 = 内容差异落在
// `.text`/`.eh_frame` 且**符号布局异** = 真代码生成差异，而 `.rodata`/`.data` **未动**。
// ② 值闸门 `r22ra__tmp/v22_gate.py`（判题同款 gcc9 编译，n=1024，64x64x3 块 = 12688 点，独立标量 oracle）：
// 基座 **bad=0 worst=0.000e+00** 与 本件 **bad=0 worst=0.000e+00** —— 逐位一致。
// ★阳性对照（律 1024）：`bad_ctrl.cpp`（`upd32<MODE1>` 的 `ADD8` 改 `SUB8`）⇒ **bad=4123,
// worst=1e9, RUN-FAIL** ⇒ 该闸门**自证比了东西**（不是「跑了但什么都没比」）。
// ③ `.cache/precheck/precheck.py` 四检全过（真 pragma 行 = 4、判题面 `int main` = 0、`-c` 编译 0 错）。
// ④ 紧邻同码锚：同批发 `f_anch.cpp`（= [A] 正文 + 一处保值词法 no-op `NMAX 4096 -> (4096)`，
// `.text` 裸段与 [A] **逐字节同**，`sec_cmp` 判「同码成立」）⇒ 为定价提供同批基准带 S。
// -------------------------------- 本机定价的如实披露 --------------------------------
// 本席本机 A/B（`taskset -c 7`，5 轮交错）：基座中位 25,291,402 拍 与 本件中位 25,363,390 拍
// （+0.28%），**但本机同码散布约 +-6%** ⇒ **本机对 0.38% 量级无分辨力**。按在册律
// 「本机无分辨力时板面定价正当」，本件由板面（紧邻同码锚件）定价，本席**不声称**本机效应。
// ================================================
// ===== REFERENCES =====
// [1] duck.ac user saffah_codex_6a_agg3, submission #126069
// <https://duck.ac/submission/126069> -- rank 1 on `mmmc1k` at judge-raw
// 6.705607 ms when this file was built. USED: the ENTIRE artifact
// VERBATIM -- every remaining line of engine code below is
// problems/mmmc1k/ref2/rival_126069.cpp byte-for-byte (all pointers,
// constants, memory layout, inline asm, loop order and summation order).
// That artifact carries its own nested citation chain, retained below; it
// in turn credits this account's public submissions #123022 and #106545,
// which it states it copied and adapted.
// [2] duck.ac user saffah_codex_6a_agg3, submission #126001
// <https://duck.ac/submission/126001> (`mmmc4k`, same account/engine
// family). USED: one IDEA ONLY, no code. Its author's own description:
// "Extract byte parity directly by shifting each 16-bit lane seven bits and
// taking byte sign bits. The two observed bits are exactly original bits 0
// and 8; all other input bits are discarded by movemask. Remove the
// preceding byte-AND-by-one from every composed B parity extraction."
// Applied here as a single-variable edit inside buildB4: the two
// `_mm256_and_si256(v,one)` operands are dropped. The 7-bit shift leaves
// the sought bit in each byte's sign position exactly as before (the
// neighbouring bits that the shift drags along are never the sign bit of a
// lane, so movemask discards them); bit-identical output was verified
// against the unmodified #126069 chassis on 6 data shapes x n=256/512/1024
// plus a scalar oracle, and on the judge machine (identical checksum).
// [3] duck.ac user saffah_cc_v41_agg1 (this account), submission #123022
// <https://duck.ac/submission/123022> at 6.783131 ms. USED: the four
// inline-asm `prefetchw` hints on the two panel writers (buildA4/buildB4
// destinations) that #126069 inherited from it, plus the IDEA of hinting a
// cold read-modify-write destination with `prefetchw` rather than
// `_MM_HINT_T0`; extended here to the epilogue's M2 rows and qwall's C rows
// (qwall's existing prefetches are kept, the exclusive hint is added).
// [4] Fleet law (this account, seat x11z_, problem 1011e6): gcc9 compiles
// `__builtin_prefetch(p,1,.)` to `prefetchnta`, which does not take the line
// exclusive -- an inline-asm `prefetchw` is required for a write hint.
// [5] duck.ac user saffah_cc_v41_agg1 (this account), submission #126888
// <https://duck.ac/submission/126888> at 6.658138 ms. USED: the ENTIRE
// artifact VERBATIM as this file's chassis -- every engine line below is
// #126888 byte-for-byte except the two edits of [6]. #126888 is that
// submission's own author's (seat r2m_) build of #126069 with the rorx
// index extraction, the dropped byte-AND in buildB4, `mk_s` forced inline
// and the `align-loops/functions=32` pragma; its citation chain is the one
// retained above.
// [6] Fleet law registry (this account, seat r7c1k_, 2026-10-15) -- NEW IDEA,
// no third-party code: **data-array page phase in the corr4_256 query
// loop**. That loop runs 24 x 32 B loads from `TB` (16 KB) plus 2 loads
// and 2 stores of `Pbb` per row-pair-word, and it sits ~2x above every
// port bound (128 ALU uops/3 ports, 24 loads/2 ports) -- the gap is on the
// load/alias side, not on the ALU ports. `TB` is the last big object in
// .bss, so moving it from 64-byte to 4096-byte alignment changes ONLY its
// own page phase: judge-machine tryrun rdtscp r0 moves -129,420 ticks
// (-0.563%, second pass -119,546 / -0.520%), output bit-identical
// (checksum 0x417b1b6398520965 on both passes). Null control: deleting the
// dead 32 KB `XT4_256` that sits directly in front of `TB` shifts it by
// exactly -8 pages (phase unchanged) and prices +-0.02% = zero, which is
// what pins the mechanism to the page phase rather than to .bss size or
// code layout. Related (but distinct) fleet laws on the same axis: 421
// (absolute-address alias on the code side) and the page-phase family.
// LICENSE: every artifact above is a public duck.ac submission (site FAQ:
// "你的代码将会被公开"); none carries a software-license notice and none
// comes from an external repository, so they are reused with full
// attribution -- account, submission address and the content used -- as
// stated above, under "all rights attributed to the original authors".
// ======================
// ===== 思路 =====
// 思路/IDEA（r2m_ 席位 · mmmc1k · n=1024 · 8-bit int matmul mod 2^8）。
// 目标：把对手 #126069（= 我方 #123022 的继承件，逐行 diff 后只有 corr4_256
// 重写 + buildB4 两处微调是"真改动"）再压低 ≥1.0%（严支 0.99*T+1µs）。
// 本文件的全部改动都是**纯 ALU/uop 减量**，逐位等价（判题机自校验和与底盘
// 完全一致，见 REFERENCES [2]/[5]），不碰任何访存形态、布局与计算次序：
//
// (a) rorx 索引提取（本席自研）。corr4_256 的查询循环每个表项要做一个
// `(w>>s)&mask` 字段提取。x86-64 的移位是两操作数、必须破坏一个副本，
// 而该循环要从同一个 64 位字里取 12 个字段 ⇒ gcc 每个字段都发
// `mov+shr+and` 三条 uop（实测循环体 116 条里 72 条是这一族）。把提取改写成
// `(w>>s)|(w<<(64-s))`（= 旋转）并给函数加 `target("avx2,bmi2")`，gcc9 发
// 3 操作数的 `rorx $imm,src,dst` ⇒ 每条目 `rorx+and` 两条 uop；同时把表基址
// 从 `__m256i*` 换成 `char*`、把索引预先放大成字节偏移，避免 gcc 再补一条
// `shl $5`。旋转带入的高位垃圾全部落在掩码之外，逐位结果不变。
// 判题机同进程 rdtscp（首次调用 r0）：-131,468 拍 = -0.567%。
// (b) 删掉 buildB4 里两处 `and(v,1)`（移植对手 #126001 的思路，见 [2]）。
// r0：-22,746 拍 = -0.098%。
// (c) 两处 `prefetchw` 扩展（见 [3]）：epilogue 的 M2 目的行（每调用写 2 行
// 128 B，上一叶之后必被面板流逐出 ⇒ 冷 RFO）与 qwall 的 C 目的行（本来是
// `_MM_HINT_T0` = shared，随后的 RMW 仍要付独占升级）。刀形与 #123022 上
// 兑现 -54.272 µs 的那把同源，只是把 hint 加到另外两个冷写目的地。
//
// (c) `mk_s` 强制内联（`__attribute__((always_inline))`，试验性）：底盘把 mk_s 编成
// 独立函数（13.5 KB 机器码），25088 次调用每次都付完整调用 ABI（5 个实参搬进
// 固定寄存器 + call/ret），且循环不变量（BI3b 基址 / PBb / C 指针）无法跨调用
// 保持在寄存器里。**asm 文本一字未改**，只改编译单位 ⇒ 目的：把调用 ABI 与
// 实参搬运从热环里去掉。（副作用：底盘带 `unroll-loops`，内联后该 asm 在
// matrix_multiply 各分支被复制 14 份，热路径只用到 1 份、其余在冷分支。）
// (d) `#pragma GCC optimize("align-loops=32,align-functions=32")`（全队律 715②，
// 试验性）：只改代码布局，不改任何指令。
// ★ 取证：底盘自身已带文件级 `#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")`
// ⇒ BMI2 全程已开，(a) 的 rorx 并不依赖函数级属性（本席实测：去掉属性后
// `g9.sh -S | grep -c rorx` = 33 不变、chk 逐位同）—— 故本件不引入任何
// 函数级 target 属性。
// (e) ★ 数据数组页相位（本席自研；见 REFERENCES [6]）。`corr4_256` 的查询环
// 每行对做 24 条 32 B 的 `TB` 表载入 + 2 条 `Pbb` 载入 + 2 条 `Pbb` 存，
// 而该环实测远高于每一条端口地板（ALU 128 uop/3 端口；载入 24/2 端口）
// ⇒ 缺口在**载入/别名侧**，不在 ALU 端口侧（同批实测反证：把 `mk_s` 的 32
// 条成对同址 BI 内存操作数改成 `vmovdqa`+寄存器形式，反而 **+0.634%**）。
// 本刀只改一个声明：`alignas(64) static __m256i TB[512]`
// → `alignas(4096)`。`TB` 是 .bss 里最后一个大对象，因此**只有它自己的页
// 相位变**（下游对象一个都不动）。判题机 tryrun 同进程 rdtscp 首次调用 r0：
// **−129,420 拍（−0.563%）**，第二趟 **−119,546（−0.520%）**，输出逐位不变
// （两趟 chk 均为 `417b1b6398520965`，与底盘同）。
// **阴性对照**：删掉紧邻 `TB` 之前、早已废弃的 32 KB `XT4_256`（它正好 8
// 页 ⇒ `TB` 位移 −32768、**页相位不变**）实测 **−0.022% ≈ 0** ⇒ 把机理钉
// 在"页相位"上，而不是 .bss 体积或代码布局。
// **未采纳**（同批判负，留档）：`#pragma GCC optimize("web")` 单独 −0.245%
// 但**与 (e) 不叠加**（`tbpw` = −0.538% < 单独 (e)）⇒ 二者同属布局/别名
// 一族，取更强的 (e)；`ira-loop-pressure` ≈0；`buildA4` 改非临时存 **+5.38%**
// （NT 存把下游要复用的面板逐出 L1）；`mk_s` unroll 族沿用 r2m_ 判负。
// 闸门：与底盘逐位一致（6 形态 × n=256/512/1024，md5 + 象限和）+ 独立标量
// oracle + 阳性对照；判题机侧自校验和亦逐位相同。
// ======================
#pragma GCC optimize("align-loops=32,align-functions=32")
#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>
#pragma GCC optimize("ira-loop-pressure")
#pragma GCC optimize("sched-pressure")
#define AVX2 __attribute__((target("avx2")))
#define NMAX 4096
static int8_t A0b[2*NMAX*NMAX] __attribute__((aligned(64)));
static uint64_t A1b[NMAX*NMAX/8] __attribute__((aligned(64)));
static uint64_t Bpb[NMAX*NMAX/8] __attribute__((aligned(64)));
static uint64_t Pbb[NMAX*NMAX/8] __attribute__((aligned(64)));
static uint8_t PBb[NMAX*NMAX] __attribute__((aligned(64)));
static void generic_mm(int n, const int8_t*A, const int8_t*B, int8_t*C){
for (int i=0;i<n;i++) for (int j=0;j<n;j++){
int s=0; for(int k=0;k<n;k++) s += (int)A[(size_t)i*n+k]*(int)B[(size_t)k*n+j];
C[(size_t)i*n+j]=(int8_t)s;
}
}
AVX2 static inline __m256i bcast16(const int8_t*p){
__m256i v; __asm__("vpbroadcastw %1, %0" : "=x"(v) : "m"(*(const short*)p)); return v;
}
AVX2 static inline __m256i packlo8(__m256i a, __m256i b){
const __m256i m=_mm256_set1_epi16(0X00FF);
return _mm256_permute4x64_epi64(
_mm256_packus_epi16(_mm256_and_si256(a,m), _mm256_and_si256(b,m)), 0xD8);
}
AVX2 static void build_A0(int n, const int8_t*A){
const __m256i m=_mm256_set1_epi8(0x7F);
size_t tot=(size_t)n*n;
for (size_t i=0;i<tot;i+=32)
_mm256_store_si256((__m256i*)(A0b+i),
_mm256_and_si256(_mm256_load_si256((const __m256i*)(A+i)), m));
}
AVX2 static void build_A1bits(int n, const int8_t*A){
int W=n/64;
for (int i=0;i<n;i++){
const int8_t*row=A+(size_t)i*n; uint64_t*w=A1b+(size_t)i*W;
for (int k=0;k<n;k+=64){
uint32_t lo=(uint32_t)_mm256_movemask_epi8(_mm256_load_si256((const __m256i*)(row+k)));
uint32_t hi=(uint32_t)_mm256_movemask_epi8(_mm256_load_si256((const __m256i*)(row+k+32)));
w[k>>6]=(uint64_t)lo|((uint64_t)hi<<32);
}
}
}
AVX2 static void build_Bpbits(int n, const int8_t*B){
const __m256i one=_mm256_set1_epi8(1);
int W=n/64;
for (int k=0;k<n;k++){
const int8_t*row=B+(size_t)k*n; uint64_t*w=Bpb+(size_t)k*W;
for (int j=0;j<n;j+=64){
__m256i a=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j)),one),7);
__m256i b=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j+32)),one),7);
uint32_t lo=(uint32_t)_mm256_movemask_epi8(a);
uint32_t hi=(uint32_t)_mm256_movemask_epi8(b);
w[j>>6]=(uint64_t)lo|((uint64_t)hi<<32);
}
}
}
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]);
}
}
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]));
}
}
}
}
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)));
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,
};
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++){
__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));
}
}
}
static int8_t BI9b[10*(NMAX/2)*192] __attribute__((aligned(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);
}
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)); }
}
}
}
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));
}
}
}
}
AVX2 static void build_A0_A1(int n, const int8_t*A){
const __m256i m=_mm256_set1_epi8(0x7F);
int W=n/64;
for (int i=0;i<n;i++){
const int8_t*row=A+(size_t)i*n; uint64_t*w=A1b+(size_t)i*W;
int8_t*o=A0b+2*(size_t)i*n;
for (int k=0;k<n;k+=64){
__m256i v0=_mm256_load_si256((const __m256i*)(row+k));
__m256i v1=_mm256_load_si256((const __m256i*)(row+k+32));
{ __m256i e0=_mm256_and_si256(v0,m), e1=_mm256_and_si256(v1,m);
__m256i l0=_mm256_unpacklo_epi16(e0,e0), h0=_mm256_unpackhi_epi16(e0,e0);
__m256i l1=_mm256_unpacklo_epi16(e1,e1), h1=_mm256_unpackhi_epi16(e1,e1);
_mm256_store_si256((__m256i*)(o+2*k), _mm256_permute2x128_si256(l0,h0,0x20));
_mm256_store_si256((__m256i*)(o+2*k+32), _mm256_permute2x128_si256(l0,h0,0x31));
_mm256_store_si256((__m256i*)(o+2*k+64), _mm256_permute2x128_si256(l1,h1,0x20));
_mm256_store_si256((__m256i*)(o+2*k+96), _mm256_permute2x128_si256(l1,h1,0x31)); }
uint32_t lo=(uint32_t)_mm256_movemask_epi8(v0);
uint32_t hi=(uint32_t)_mm256_movemask_epi8(v1);
w[k>>6]=(uint64_t)lo|((uint64_t)hi<<32);
}
}
}
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 __attribute__((always_inline)) static inline void mk_s(int n, const int8_t*A0, const int8_t*BI, int8_t*C, const uint8_t*PB, int jb){ (void)PB;
const int8_t *x0=A0, *x1=A0+2*n, *bi=BI; long cnt=n>>4;
asm volatile(
"vpxor %%xmm8,%%xmm8,%%xmm8\n\t" "vpxor %%xmm9,%%xmm9,%%xmm9\n\t" "vpxor %%xmm10,%%xmm10,%%xmm10\n\t" "vpxor %%xmm11,%%xmm11,%%xmm11\n\t" "vpxor %%xmm12,%%xmm12,%%xmm12\n\t" "vpxor %%xmm13,%%xmm13,%%xmm13\n\t" "vpxor %%xmm14,%%xmm14,%%xmm14\n\t" "vpxor %%xmm15,%%xmm15,%%xmm15\n\t"
".p2align 5\n\t" "1:\n\t"
"prefetcht0 1536(%[bi])\n\t"
"vpbroadcastd (%[x0]), %%ymm0\n\t"
"vpbroadcastd (%[x1]), %%ymm1\n\t"
"vpmaddubsw 0(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 0(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 32(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 32(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 64(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 64(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 96(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 96(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 4(%[x0]), %%ymm0\n\t"
"vpbroadcastd 4(%[x1]), %%ymm1\n\t"
"vpmaddubsw 128(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 128(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 160(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 160(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 192(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 192(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 224(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 224(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 8(%[x0]), %%ymm0\n\t"
"vpbroadcastd 8(%[x1]), %%ymm1\n\t"
"vpmaddubsw 256(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 256(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 288(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 288(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 320(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 320(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 352(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 352(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 12(%[x0]), %%ymm0\n\t"
"vpbroadcastd 12(%[x1]), %%ymm1\n\t"
"vpmaddubsw 384(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 384(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 416(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 416(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 448(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 448(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 480(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 480(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 16(%[x0]), %%ymm0\n\t"
"vpbroadcastd 16(%[x1]), %%ymm1\n\t"
"vpmaddubsw 512(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 512(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 544(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 544(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 576(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 576(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 608(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 608(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 20(%[x0]), %%ymm0\n\t"
"vpbroadcastd 20(%[x1]), %%ymm1\n\t"
"vpmaddubsw 640(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 640(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 672(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 672(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 704(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 704(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 736(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 736(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 24(%[x0]), %%ymm0\n\t"
"vpbroadcastd 24(%[x1]), %%ymm1\n\t"
"vpmaddubsw 768(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 768(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 800(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 800(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 832(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 832(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 864(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 864(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"vpbroadcastd 28(%[x0]), %%ymm0\n\t"
"vpbroadcastd 28(%[x1]), %%ymm1\n\t"
"vpmaddubsw 896(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 896(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
"vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
"vpmaddubsw 928(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 928(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
"vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
"vpmaddubsw 960(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 960(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
"vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
"vpmaddubsw 992(%[bi]), %%ymm0, %%ymm5\n\t"
"vpmaddubsw 992(%[bi]), %%ymm1, %%ymm6\n\t"
"vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
"vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
"addq $32, %[x0]\n\t"
"addq $32, %[x1]\n\t"
"addq $1024, %[bi]\n\t"
"decq %[cnt]\n\t"
"jnz 1b\n\t"
"vpand %[mask], %%ymm8, %%ymm2\n\t"
"vpand %[mask], %%ymm9, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 0(%[scr])\n\t"
"vpand %[mask], %%ymm10, %%ymm2\n\t"
"vpand %[mask], %%ymm11, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 32(%[scr])\n\t"
"vpand %[mask], %%ymm12, %%ymm2\n\t"
"vpand %[mask], %%ymm13, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 64(%[scr])\n\t"
"vpand %[mask], %%ymm14, %%ymm2\n\t"
"vpand %[mask], %%ymm15, %%ymm3\n\t"
"vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
"vpermq $0xD8, %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm2, 96(%[scr])\n\t"
: [x0]"+r"(x0),[x1]"+r"(x1),[bi]"+r"(bi),[cnt]"+r"(cnt)
: [mask]"m"(W0FF), [scr]"r"(mscr)
: "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15","cc","memory");
}
static uint8_t Gtab[8*256*(NMAX/8)] __attribute__((aligned(64)));
#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;
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);
}
}
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);
}
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);
}
}
}
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;
uint64_t *w0p=BP2+(size_t)((2*p)>>6)*(size_t)(65*W)+(size_t)((2*p)&63)*W;
uint64_t *w1p=BP2+(size_t)((2*p+1)>>6)*(size_t)(65*W)+(size_t)((2*p+1)&63)*W;
for (int g=0;g<ng;g++){
const int8_t *a=ax+64*g, *c=bx+64*g;
const int8_t *d=ay+64*g, *e=by+64*g;
int8_t *o=BI3b+((size_t)g*(m/2)+p)*128;
uint64_t parity0=0,parity1=0;
for (int t=0;t<64;t+=32){
__m256i va=sblk_cmb(_mm256_load_si256((const __m256i*)(a+t)), d+t, b->op);
__m256i vb=sblk_cmb(_mm256_load_si256((const __m256i*)(c+t)), e+t, b->op);
__m256i lo=_mm256_unpacklo_epi8(va,vb), hi=_mm256_unpackhi_epi8(va,vb);
_mm256_store_si256((__m256i*)(o+2*t), _mm256_permute2x128_si256(lo,hi,0x20));
_mm256_store_si256((__m256i*)(o+2*t+32), _mm256_permute2x128_si256(lo,hi,0x31));
uint32_t m0=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(va,7));
uint32_t m1=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(vb,7));
parity0 |= (uint64_t)m0<<t; parity1 |= (uint64_t)m1<<t;
}
w0p[g]=parity0; w1p[g]=parity1;
}
}
}
static const int SAB_T[7][5] = {
{0,0, 1,1, 1},
{1,0, 1,1, 1},
{0,0, 0,0, 0},
{1,1, 1,1, 0},
{0,0, 0,1, 1},
{1,0, 0,0, 2},
{0,1, 1,1, 2},
};
static const int SBB_T[7][5] = {
{0,0, 1,1, 1},
{0,0, 0,0, 0},
{0,1, 1,1, 2},
{1,0, 0,0, 2},
{1,1, 1,1, 0},
{0,0, 0,1, 1},
{1,0, 1,1, 1},
};
static int8_t M2[7][256*256+384] __attribute__((aligned(64)));
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);
}
}
}
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);
}
}
}
static __m256i XT4_256[4][16][16] __attribute__((aligned(64)));
alignas(4096) static __m256i TB[512];
AVX2 __attribute__((target("avx2,bmi2"),noinline)) static void corr4_256(){
const char*TBc=(const char*)TB;
#define R2M_ROT(w,sh) ((unsigned long long)(((w)>>(sh))|((w)<<(64-(sh)))))
constexpr int n=256,W=4;
for(int wb=0;wb<4;wb++){
for(int g=0;g<12;g++){
int bit=g<4?g*6:24+(g-4)*5;
__m256i*T=TB+(g<4?g*64:256+(g-4)*32);
__m256i r0=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+bit+0)*W));
__m256i r1=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+bit+1)*W));
__m256i r2=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+bit+2)*W));
__m256i r3=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+bit+3)*W));
__m256i r4=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+bit+4)*W));
__m256i v=_mm256_setzero_si256();_mm256_store_si256(T,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+1,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+3,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+2,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+6,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+7,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+5,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+4,v);
v=_mm256_xor_si256(v,r3);_mm256_store_si256(T+12,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+13,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+15,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+14,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+10,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+11,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+9,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+8,v);
v=_mm256_xor_si256(v,r4);_mm256_store_si256(T+24,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+25,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+27,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+26,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+30,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+31,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+29,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+28,v);
v=_mm256_xor_si256(v,r3);_mm256_store_si256(T+20,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+21,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+23,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+22,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+18,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+19,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+17,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+16,v);
if(g<4){
__m256i r5=_mm256_load_si256((const __m256i*)(BP2+((size_t)wb*65+bit+5)*W));
v=_mm256_xor_si256(v,r5);_mm256_store_si256(T+48,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+49,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+51,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+50,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+54,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+55,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+53,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+52,v);
v=_mm256_xor_si256(v,r3);_mm256_store_si256(T+60,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+61,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+63,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+62,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+58,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+59,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+57,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+56,v);
v=_mm256_xor_si256(v,r4);_mm256_store_si256(T+40,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+41,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+43,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+42,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+46,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+47,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+45,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+44,v);
v=_mm256_xor_si256(v,r3);_mm256_store_si256(T+36,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+37,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+39,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+38,v);
v=_mm256_xor_si256(v,r2);_mm256_store_si256(T+34,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+35,v);
v=_mm256_xor_si256(v,r1);_mm256_store_si256(T+33,v);
v=_mm256_xor_si256(v,r0);_mm256_store_si256(T+32,v);
}
}
for(int i=0;i<n;i+=2){
uint64_t w0=A1b[(size_t)i*W+wb],w1=A1b[(size_t)(i+1)*W+wb];
__m256i x0=wb?_mm256_load_si256((const __m256i*)(Pbb+(size_t)i*W)):_mm256_setzero_si256();
__m256i y0=wb?_mm256_load_si256((const __m256i*)(Pbb+(size_t)(i+1)*W)):_mm256_setzero_si256();
__m256i x1=_mm256_setzero_si256(),y1=x1;
x0=_mm256_xor_si256(x0,_mm256_load_si256((const __m256i*)(TBc+0+(R2M_ROT(w0,59)&0x7e0))));
y0=_mm256_xor_si256(y0,_mm256_load_si256((const __m256i*)(TBc+0+(R2M_ROT(w1,59)&0x7e0))));
x1=_mm256_xor_si256(x1,_mm256_load_si256((const __m256i*)(TBc+2048+(R2M_ROT(w0,1)&0x7e0))));
y1=_mm256_xor_si256(y1,_mm256_load_si256((const __m256i*)(TBc+2048+(R2M_ROT(w1,1)&0x7e0))));
x0=_mm256_xor_si256(x0,_mm256_load_si256((const __m256i*)(TBc+4096+(R2M_ROT(w0,7)&0x7e0))));
y0=_mm256_xor_si256(y0,_mm256_load_si256((const __m256i*)(TBc+4096+(R2M_ROT(w1,7)&0x7e0))));
x1=_mm256_xor_si256(x1,_mm256_load_si256((const __m256i*)(TBc+6144+(R2M_ROT(w0,13)&0x7e0))));
y1=_mm256_xor_si256(y1,_mm256_load_si256((const __m256i*)(TBc+6144+(R2M_ROT(w1,13)&0x7e0))));
x0=_mm256_xor_si256(x0,_mm256_load_si256((const __m256i*)(TBc+8192+(R2M_ROT(w0,19)&0x3e0))));
y0=_mm256_xor_si256(y0,_mm256_load_si256((const __m256i*)(TBc+8192+(R2M_ROT(w1,19)&0x3e0))));
x1=_mm256_xor_si256(x1,_mm256_load_si256((const __m256i*)(TBc+9216+(R2M_ROT(w0,24)&0x3e0))));
y1=_mm256_xor_si256(y1,_mm256_load_si256((const __m256i*)(TBc+9216+(R2M_ROT(w1,24)&0x3e0))));
x0=_mm256_xor_si256(x0,_mm256_load_si256((const __m256i*)(TBc+10240+(R2M_ROT(w0,29)&0x3e0))));
y0=_mm256_xor_si256(y0,_mm256_load_si256((const __m256i*)(TBc+10240+(R2M_ROT(w1,29)&0x3e0))));
x1=_mm256_xor_si256(x1,_mm256_load_si256((const __m256i*)(TBc+11264+(R2M_ROT(w0,34)&0x3e0))));
y1=_mm256_xor_si256(y1,_mm256_load_si256((const __m256i*)(TBc+11264+(R2M_ROT(w1,34)&0x3e0))));
x0=_mm256_xor_si256(x0,_mm256_load_si256((const __m256i*)(TBc+12288+(R2M_ROT(w0,39)&0x3e0))));
y0=_mm256_xor_si256(y0,_mm256_load_si256((const __m256i*)(TBc+12288+(R2M_ROT(w1,39)&0x3e0))));
x1=_mm256_xor_si256(x1,_mm256_load_si256((const __m256i*)(TBc+13312+(R2M_ROT(w0,44)&0x3e0))));
y1=_mm256_xor_si256(y1,_mm256_load_si256((const __m256i*)(TBc+13312+(R2M_ROT(w1,44)&0x3e0))));
x0=_mm256_xor_si256(x0,_mm256_load_si256((const __m256i*)(TBc+14336+(R2M_ROT(w0,49)&0x3e0))));
y0=_mm256_xor_si256(y0,_mm256_load_si256((const __m256i*)(TBc+14336+(R2M_ROT(w1,49)&0x3e0))));
x1=_mm256_xor_si256(x1,_mm256_load_si256((const __m256i*)(TBc+15360+(R2M_ROT(w0,54)&0x3e0))));
y1=_mm256_xor_si256(y1,_mm256_load_si256((const __m256i*)(TBc+15360+(R2M_ROT(w1,54)&0x3e0))));
_mm256_store_si256((__m256i*)(Pbb+(size_t)i*W),_mm256_xor_si256(x0,x1));
_mm256_store_si256((__m256i*)(Pbb+(size_t)(i+1)*W),_mm256_xor_si256(y0,y1));
}
}
}
static inline int sgn_of(int op){ return op==0 ? 0 : (op==1 ? 1 : -1); }
AVX2 static void buildA4(int h, int n, const int8_t*A, int k1, int k2){
const __m256i msk=_mm256_set1_epi8(0x7F);
const int o1 = 512;
const int *t1=SAB_T[k1], *t2=SAB_T[k2];
const int s1=sgn_of(t1[4]), s2=sgn_of(t2[4]), s12=s1*s2;
const size_t R1=(size_t)(t1[0]*o1+t2[0]*h)*n + t1[1]*o1+t2[1]*h;
const size_t R2=(size_t)(t1[2]*o1+t2[0]*h)*n + t1[3]*o1+t2[1]*h;
const size_t R3=(size_t)(t1[0]*o1+t2[2]*h)*n + t1[1]*o1+t2[3]*h;
const size_t R4=(size_t)(t1[2]*o1+t2[2]*h)*n + t1[3]*o1+t2[3]*h;
int W=h/64;
for (int i=0;i<h;i++){
const int8_t*p1=A+R1+(size_t)i*n;
const int8_t*p2=A+R2+(size_t)i*n;
const int8_t*p3=A+R3+(size_t)i*n;
const int8_t*p4=A+R4+(size_t)i*n;
int8_t*dst=A0b+2*(size_t)i*h;
uint64_t*w1=A1b+(size_t)i*W;
for (int k=0;k<h;k+=64){
_mm_prefetch((const char*)(p1+k+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(p2+k+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(p3+k+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(p4+k+4096), _MM_HINT_T0);
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(dst+2*k+256)));
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(dst+2*k+320)));
__m256i v0=_mm256_load_si256((const __m256i*)(p1+k));
__m256i v1=_mm256_load_si256((const __m256i*)(p1+k+32));
if(s1){ __m256i u0=_mm256_load_si256((const __m256i*)(p2+k));
__m256i u1=_mm256_load_si256((const __m256i*)(p2+k+32));
if(s1>0){v0=_mm256_add_epi8(v0,u0);v1=_mm256_add_epi8(v1,u1);}
else {v0=_mm256_sub_epi8(v0,u0);v1=_mm256_sub_epi8(v1,u1);} }
if(s2){ __m256i u0=_mm256_load_si256((const __m256i*)(p3+k));
__m256i u1=_mm256_load_si256((const __m256i*)(p3+k+32));
if(s2>0){v0=_mm256_add_epi8(v0,u0);v1=_mm256_add_epi8(v1,u1);}
else {v0=_mm256_sub_epi8(v0,u0);v1=_mm256_sub_epi8(v1,u1);} }
if(s12){__m256i u0=_mm256_load_si256((const __m256i*)(p4+k));
__m256i u1=_mm256_load_si256((const __m256i*)(p4+k+32));
if(s12>0){v0=_mm256_add_epi8(v0,u0);v1=_mm256_add_epi8(v1,u1);}
else {v0=_mm256_sub_epi8(v0,u0);v1=_mm256_sub_epi8(v1,u1);} }
{
__m256i e0=_mm256_and_si256(v0,msk), e1=_mm256_and_si256(v1,msk);
__m256i l0=_mm256_unpacklo_epi16(e0,e0), h0=_mm256_unpackhi_epi16(e0,e0);
__m256i l1=_mm256_unpacklo_epi16(e1,e1), h1=_mm256_unpackhi_epi16(e1,e1);
_mm256_store_si256((__m256i*)(dst+2*k), _mm256_permute2x128_si256(l0,h0,0x20));
_mm256_store_si256((__m256i*)(dst+2*k+32), _mm256_permute2x128_si256(l0,h0,0x31));
_mm256_store_si256((__m256i*)(dst+2*k+64), _mm256_permute2x128_si256(l1,h1,0x20));
_mm256_store_si256((__m256i*)(dst+2*k+96), _mm256_permute2x128_si256(l1,h1,0x31)); }
w1[k>>6]=(uint64_t)(uint32_t)_mm256_movemask_epi8(v0)
|((uint64_t)(uint32_t)_mm256_movemask_epi8(v1)<<32);
}
}
}
AVX2 static void buildB4(int h, int n, const int8_t*B, int k1, int k2){
const __m256i one=_mm256_set1_epi8(1);
const int o1 = 512;
const int *t1=SBB_T[k1], *t2=SBB_T[k2];
const int s1=sgn_of(t1[4]), s2=sgn_of(t2[4]), s12=s1*s2;
const size_t R1=(size_t)(t1[0]*o1+t2[0]*h)*n + t1[1]*o1+t2[1]*h;
const size_t R2=(size_t)(t1[2]*o1+t2[0]*h)*n + t1[3]*o1+t2[1]*h;
const size_t R3=(size_t)(t1[0]*o1+t2[2]*h)*n + t1[1]*o1+t2[3]*h;
const size_t R4=(size_t)(t1[2]*o1+t2[2]*h)*n + t1[3]*o1+t2[3]*h;
int ng=h/64, W=h/64;
for (int p=0;p<h/2;p++){
const int8_t *a1=B+R1+(size_t)(2*p)*n, *b1=B+R1+(size_t)(2*p+1)*n;
const int8_t *a2=B+R2+(size_t)(2*p)*n, *b2=B+R2+(size_t)(2*p+1)*n;
const int8_t *a3=B+R3+(size_t)(2*p)*n, *b3=B+R3+(size_t)(2*p+1)*n;
const int8_t *a4=B+R4+(size_t)(2*p)*n, *b4=B+R4+(size_t)(2*p+1)*n;
uint64_t *w0p=BP2+(size_t)((2*p)>>6)*(size_t)(65*W)+(size_t)((2*p)&63)*W;
uint64_t *w1p=BP2+(size_t)((2*p+1)>>6)*(size_t)(65*W)+(size_t)((2*p+1)&63)*W;
for (int g=0;g<ng;g++){
const int8_t *a=a1+64*g, *c=b1+64*g;
const int8_t *d=a2+64*g, *e=b2+64*g;
const int8_t *f=a3+64*g, *h3=b3+64*g;
const int8_t *m4=a4+64*g, *n4=b4+64*g;
int8_t *o=BI3b+((size_t)g*(h/2)+p)*128;
uint64_t parity0=0,parity1=0;
_mm_prefetch((const char*)(a +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(c +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(d +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(e +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(f +4096), _MM_HINT_T0);
_mm_prefetch((const char*)(h3+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(m4+4096), _MM_HINT_T0);
_mm_prefetch((const char*)(n4+4096), _MM_HINT_T0);
for (int t=0;t<64;t+=32){
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(o+2*t+256)));
__asm__ volatile("prefetchw %0" :: "m"(*(char*)(o+2*t+320)));
__m256i va=_mm256_load_si256((const __m256i*)(a+t));
__m256i vb=_mm256_load_si256((const __m256i*)(c+t));
if(s1){ __m256i u0=_mm256_load_si256((const __m256i*)(d+t));
__m256i u1=_mm256_load_si256((const __m256i*)(e+t));
if(s1>0){va=_mm256_add_epi8(va,u0);vb=_mm256_add_epi8(vb,u1);}
else {va=_mm256_sub_epi8(va,u0);vb=_mm256_sub_epi8(vb,u1);} }
if(s2){ __m256i u0=_mm256_load_si256((const __m256i*)(f+t));
__m256i u1=_mm256_load_si256((const __m256i*)(h3+t));
if(s2>0){va=_mm256_add_epi8(va,u0);vb=_mm256_add_epi8(vb,u1);}
else {va=_mm256_sub_epi8(va,u0);vb=_mm256_sub_epi8(vb,u1);} }
if(s12){__m256i u0=_mm256_load_si256((const __m256i*)(m4+t));
__m256i u1=_mm256_load_si256((const __m256i*)(n4+t));
if(s12>0){va=_mm256_add_epi8(va,u0);vb=_mm256_add_epi8(vb,u1);}
else {va=_mm256_sub_epi8(va,u0);vb=_mm256_sub_epi8(vb,u1);} }
__m256i lo=_mm256_unpacklo_epi8(va,vb), hi=_mm256_unpackhi_epi8(va,vb);
_mm256_store_si256((__m256i*)(o+2*t), _mm256_permute2x128_si256(lo,hi,0x20));
_mm256_store_si256((__m256i*)(o+2*t+32), _mm256_permute2x128_si256(lo,hi,0x31));
uint32_t m0=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(va,7));
uint32_t m1=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(vb,7));
parity0 |= (uint64_t)m0<<t; parity1 |= (uint64_t)m1<<t;
}
w0p[g]=parity0; w1p[g]=parity1;
}
}
}
AVX2 static void submm512_4(int m, int n, const int8_t*A, const int8_t*B, int k1, int k2,
int8_t*Cp){
buildA4(m,n,A,k1,k2); buildB4(m,n,B,k1,k2);
{ const int W=m>>6; for (int g=0; g<(m>>6); g++) memset(BP2+(size_t)g*65*W+64*W, 0, (size_t)W*8); }
if(m==256) corr4_256(); else corr_rus(m);
for (int g=0; g<m; g+=64)
for (int i=0; i<m; i+=2){
mk_s(m, A0b+2*(size_t)i*m, BI3b+(size_t)(g>>6)*(m/2)*128, Cp, PBb, g);
epilogue(m, Cp, (const uint64_t*)Pbb, g, i, 2);
}
}
struct KD { int R, Cc, M; };
static constexpr KD KDEST[7][2] = {
{{ 0, 0, 0}, {512, 512, 0}},
{{512, 0, 0}, {512, 512, -1}},
{{ 0, 512, 0}, {512, 512, 1}},
{{ 0, 0, 1}, {512, 0, 1}},
{{ 0, 0,-1}, { 0, 512, 1}},
{{512, 512, 1}, { 0, 0, 0}},
{{ 0, 0, 1}, { 0, 0, 0}},
};
static constexpr int KDN[7] = {2,2,2,2,2,1,1};
#define LD32(p) _mm256_load_si256((const __m256i*)(p))
#define ST32(p,v) _mm256_store_si256((__m256i*)(p), (v))
#define ADD8(a,b) _mm256_add_epi8((a),(b))
#define SUB8(a,b) _mm256_sub_epi8((a),(b))
template<int MODE> static inline void upd32(int8_t*p, __m256i v){
if(MODE==0) ST32(p, v);
else if(MODE==1) ST32(p, ADD8(LD32(p), v));
else ST32(p, SUB8(LD32(p), v));
}
template<int MODE1,int MODE2>
AVX2 static void qwall(int8_t*C,int n,int m,int8_t*D2){
const int h=m>>1;
int8_t *a1=C,*a2=C+h,*a3=C+(size_t)h*n,*a4=C+(size_t)h*n+h;
int8_t *b1=0,*b2=0,*b3=0,*b4=0;
if(D2){ b1=D2; b2=D2+h; b3=D2+(size_t)h*n; b4=D2+(size_t)h*n+h; }
for(int i=0;i<h;i++){
const int8_t *p0=M2[0]+(size_t)i*h,*p1=M2[1]+(size_t)i*h,*p2=M2[2]+(size_t)i*h,
*p3=M2[3]+(size_t)i*h,*p4=M2[4]+(size_t)i*h,*p5=M2[5]+(size_t)i*h,
*p6=M2[6]+(size_t)i*h;
int8_t *e1=a1+(size_t)i*n,*e2=a2+(size_t)i*n,*e3=a3+(size_t)i*n,*e4=a4+(size_t)i*n;
int8_t *f1=b1?b1+(size_t)i*n:0,*f2=b2?b2+(size_t)i*n:0,
*f3=b3?b3+(size_t)i*n:0,*f4=b4?b4+(size_t)i*n:0;
if(i+2<h){
const char*q0=(const char*)(M2[0]+(size_t)(i+2)*h);
_mm_prefetch(q0,_MM_HINT_T0); _mm_prefetch(q0+64,_MM_HINT_T0); _mm_prefetch(q0+128,_MM_HINT_T0); _mm_prefetch(q0+192,_MM_HINT_T0);
_mm_prefetch((const char*)(M2[1]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[2]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[3]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[4]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[5]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(M2[6]+(size_t)(i+2)*h),_MM_HINT_T0);
_mm_prefetch((const char*)(e1+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(e2+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(e3+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(e4+(size_t)2*n),_MM_HINT_T0);
if(f1){ _mm_prefetch((const char*)(f1+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(f2+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(f3+(size_t)2*n),_MM_HINT_T0);
_mm_prefetch((const char*)(f4+(size_t)2*n),_MM_HINT_T0); }
}
for(int j=0;j<h;j+=32){
const __m256i v0=LD32(p0+j),v1=LD32(p1+j),v2=LD32(p2+j),v3=LD32(p3+j),
v4=LD32(p4+j),v5=LD32(p5+j),v6=LD32(p6+j);
const __m256i y00=ADD8(ADD8(v0,v3),SUB8(v6,v4));
const __m256i y01=ADD8(v2,v4);
const __m256i y10=ADD8(v1,v3);
const __m256i y11=ADD8(ADD8(SUB8(v0,v1),v2),v5);
upd32<MODE1>(e1+j,y00); upd32<MODE1>(e2+j,y01);
upd32<MODE1>(e3+j,y10); upd32<MODE1>(e4+j,y11);
if(f1){ upd32<MODE2>(f1+j,y00); upd32<MODE2>(f2+j,y01);
upd32<MODE2>(f3+j,y10); upd32<MODE2>(f4+j,y11); }
}
}
}
template<int K1>
AVX2 static void strassen2_512x(const int8_t*A, const int8_t*B, int n, int8_t*C){
const int h=256;
for(int k=0;k<7;k++) submm512_4(h,n,A,B,K1,k,M2[k]);
int8_t *D1=C+(size_t)KDEST[K1][0].R*n+KDEST[K1][0].Cc;
int8_t *D2=(KDN[K1]>1)?C+(size_t)KDEST[K1][1].R*n+KDEST[K1][1].Cc:0;
if(KDN[K1]>1) qwall<KDEST[K1][0].M,KDEST[K1][1].M>(D1,n,512,D2);
else qwall<KDEST[K1][0].M,0>(D1,n,512,0);
}
AVX2 static void strassen1(int n, const int8_t*A, const int8_t*B, int8_t*C){
strassen2_512x<0>(A,B,n,C); strassen2_512x<1>(A,B,n,C);
strassen2_512x<2>(A,B,n,C); strassen2_512x<3>(A,B,n,C);
strassen2_512x<4>(A,B,n,C); strassen2_512x<5>(A,B,n,C);
strassen2_512x<6>(A,B,n,C);
_mm_sfence();
}
AVX2 void matrix_multiply(int n, const int8_t *A, const int8_t *B, int8_t *C){
if (n<=0 || (n&255) || n>NMAX) { generic_mm(n,A,B,C); return; }
if (n==1024) { strassen1(n,A,B,C); return; }
build_A0_A1(n,A);
if(n==1024){
build_BI9(n,B); build_BI3_tail(n,B); build_BP2(n,B); corr_rus(n);
for (int g=0; g<10; g++)
for (int i=0; i<n; i+=2){
mk_s6(n, A0b+(size_t)i*n, BI9b+(size_t)g*(n/2)*192, (uint8_t*)mscr);
epilogue(n, C, (const uint64_t*)Pbb, g*96, i, 3);
}
for (int i=0; i<n; i+=2){
mk_s(n, A0b+2*(size_t)i*n, BI3b, C, PBb, 0);
epilogue(n, C, (const uint64_t*)Pbb, 960, i, 2);
}
return;
}
build_BI3(n,B); build_Bpbits(n,B);
if(n==4096) corr16(); else corr_pass3(n);
for (int g=0, jbg=0; g<n; g+=64, jbg=g)
for (int i=0; i<n; i+=2)
{ mk_s(n, A0b+2*(size_t)i*n, BI3b+(size_t)(g>>6)*(n/2)*128, C, PBb, g);
epilogue(n, C, (const uint64_t*)Pbb, g, i, 2); }
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 6.648 ms | 1 MB + 712 KB | Accepted | Score: 100 | 显示更多 |