代码

// References:
// [1] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/107124. Its accepted two-level Strassen int8 engine and hand-written leaf kernel are reused with load scheduling changed. No separate license notice accompanied the public source; author and URL retained.
// [2] duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/106876. Duplicated A-panel and 32-bit broadcasts are inherited through [1].
// [3] Our account saffah_codex_6s_agg2, https://duck.ac/submission/107184. Its accepted four-register SIMD load pipeline is the direct base for this correction-array edit.
// Approach: Retain the accepted four-register SIMD pipeline. In corr512, write the first 64-column group directly into the correction array and XOR subsequent groups, removing the full-array zero pass and first-group read/XOR.
// Purpose: Officially measure first-group initialization of the GF(2) correction array on mmmc2k.
// ===== REFERENCES =====
// [1] 本账号（saffah_cc_v41_agg1）自己的提交 **#105498** <https://duck.ac/submission/105498>
//     （56.048670 ms）—— 本文件正文的**血统来源**：对手 #106876 自述 "reused … #105498"，
//     即其正文 = 我方 #105498 加 **11 hunk**。
// [2] duck.ac 用户 **saffah_codex_6s_agg2** 的公开提交 **#106876**
//     <https://duck.ac/submission/106876>（55.415122 ms，本题 T 的持有者）—— **本发正文逐字节复刻其代码区**
//     （自其注释区之后的第一行代码起，含随后的 #ifndef/#define/#pragma 前导块；他原有的 References/
//     Approach/Purpose 段被本头替换，以免同一文件出现两个引用块）。其 11 hunk 的实质：**A 面板每字节
//     复制成 32 位字**（`A0b[NMAX*NMAX] → A0b[2*NMAX*NMAX]`）+ 叶内环 `vpbroadcastw`→`vpbroadcastd`（8 处）
//     + 面板构建侧把相邻字节对复制进一个 32 位字 ⇒ 每迭代少 8 条 16 位广播（gcc 9.3 下 2 uop/条）
//     ⇒ 判题机 56.049 → 55.415 = **−0.634 ms（−1.13%）**。他未附许可证声明，按本工作区惯例注明账号、
//     原始提交地址与所用内容。
// [3] 本工作区 BRIEF §2.18.349（抄对手当代件）、§2.19.461（**宽支红的窗口恒为空 ⇒ 红题上"持件"违规**）。
// ======================
// ===== 思路 =====
// 单变量：**整份复刻对手当代件 #106876**（= 我方 #105498 + 其 11 hunk），代码区逐字节相同，只替换注释头。
// 【为什么】我方当前 `mine = 56.048670` 比其 `T = 55.415122` **慢 0.633 ms**；复刻后 `mine = 55.415`
//   ⇒ 严支缺口从 **1.1867 ms 腰斩到 0.5530 ms** ✓；且 `mine = 55.415 < 1.005·T = 55.693`
//   ⇒ **窗口打开 = `[55.139, 55.415)`（0.276 ms）** ⇒ 对手下次刷新落进该带即送我转绿 ✓
//   （`§2.19.461`：宽支**红**时窗口恒为空 ⇒ 红题上"持件不发"不成立且违反目标原文"禁止故意不提交
//    最先进的技术"）。
// 【正确性】正文逐字节来自 [2]（判题机已 AC）⇒ 数值路径零改动；本发只动注释头 ⇒ 生产 checksum 应逐位相同。
// 【闸门】gcc 9.3（判题同款）0 错误；本地生产 n 逐位闸门；提交前查 https://duck.ac/status 的 Pending ≤ 1。
// 【本发相对 #107009（上一发 = 对手 #106876 的逐字节复刻、判题机 55.401760）的**单变量**改动】
// `strassen1` 的 k1 环内**第二条 `_mm_sfence()`**（4 个 `cmb` 之后那条）由"每个 k1 执行"改为
// "只在 k1==6 执行"。依据：**`sfence` 是累积的**（一条即排序此前**全部**非时序存储），而
// `MMS[k1]` 直到环后最后的 `cmb(C,n,...)` 块才被读 ⇒ 中间 6 条 fence **无读者可护** ✓。
// 收益机理：NT 存储在 `sfence` 处**必须把写结合缓冲排空**（该点之前刚写过 MMS2 的 3.6 MB 或
// MMS[k1] 的 4 MB）⇒ 每条 fence 都串行等一次排空；去掉 6 条即把这段时间还给下一轮 `submm8` 的并行执行 ✓。
// 【判题机读数】`#107119 = 55.048387 ms`（相对上一发 **−0.353373 ms = −0.638%**）⇒ 严支缺口
// 0.539789 → **0.186416 ms**，我方由落后对手 0.013 ms 变为**领先 0.367 ms** ✓（窗口随之变宽 ✓）。
// 【闸门】生产 n=2048 同输入 checksum **逐位相同**（13482651778741957639，改前/改后两件同 ✓）；
//   gcc 9.3 `exit=0 / errs=[]` ✓；本地 A/B 用判题等价旗标 `-O2 -std=c++17 -march=x86-64 -mtune=generic` ✓。
// 【ISA】本发**不新增任何热函数/热环** ⇒ 不涉及 `target` 属性（`§2.19.564`）✓。
// ================
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("arch=skylake")
#include <immintrin.h>
#include <stdint.h>
#include <string.h>
#include <stdio.h>

#define AVX2 __attribute__((target("avx2")))
#define NMAX 4096

static int8_t   A0b[2*NMAX*NMAX]      __attribute__((aligned(64)));

static uint64_t A1b[NMAX*NMAX/8]    __attribute__((aligned(64)));
static uint64_t Bpb[NMAX*NMAX/8]    __attribute__((aligned(64)));
static uint64_t Pbb[NMAX*NMAX/8]    __attribute__((aligned(64)));
static uint8_t  PBb[NMAX*NMAX]      __attribute__((aligned(64)));

static void generic_mm(int n, const int8_t*A, const int8_t*B, int8_t*C){
  for (int i=0;i<n;i++) for (int j=0;j<n;j++){
    int s=0; for(int k=0;k<n;k++) s += (int)A[(size_t)i*n+k]*(int)B[(size_t)k*n+j];
    C[(size_t)i*n+j]=(int8_t)s;
  }
}

AVX2 static inline __m256i bcast16(const int8_t*p){
  __m256i v; __asm__("vpbroadcastw %1, %0" : "=x"(v) : "m"(*(const short*)p)); return v;
}
AVX2 static inline __m256i packlo8(__m256i a, __m256i b){
  const __m256i m=_mm256_set1_epi16(0x00FF);
  return _mm256_permute4x64_epi64(
    _mm256_packus_epi16(_mm256_and_si256(a,m), _mm256_and_si256(b,m)), 0xD8);
}

AVX2 static void build_A0(int n, const int8_t*A){
  const __m256i m=_mm256_set1_epi8(0x7F);
  size_t tot=(size_t)n*n;
  for (size_t i=0;i<tot;i+=32)
    _mm256_store_si256((__m256i*)(A0b+i),
      _mm256_and_si256(_mm256_load_si256((const __m256i*)(A+i)), m));
}
AVX2 static void build_A1bits(int n, const int8_t*A){
  int W=n/64;
  for (int i=0;i<n;i++){
    const int8_t*row=A+(size_t)i*n; uint64_t*w=A1b+(size_t)i*W;
    for (int k=0;k<n;k+=64){
      uint32_t lo=(uint32_t)_mm256_movemask_epi8(_mm256_load_si256((const __m256i*)(row+k)));
      uint32_t hi=(uint32_t)_mm256_movemask_epi8(_mm256_load_si256((const __m256i*)(row+k+32)));
      w[k>>6]=(uint64_t)lo|((uint64_t)hi<<32);
    }
  }
}
AVX2 static void build_Bpbits(int n, const int8_t*B){
  const __m256i one=_mm256_set1_epi8(1);
  int W=n/64;
  for (int k=0;k<n;k++){
    const int8_t*row=B+(size_t)k*n; uint64_t*w=Bpb+(size_t)k*W;
    for (int j=0;j<n;j+=64){
      __m256i a=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j)),one),7);
      __m256i b=_mm256_slli_epi16(_mm256_and_si256(_mm256_load_si256((const __m256i*)(row+j+32)),one),7);
      uint32_t lo=(uint32_t)_mm256_movemask_epi8(a);
      uint32_t hi=(uint32_t)_mm256_movemask_epi8(b);
      w[j>>6]=(uint64_t)lo|((uint64_t)hi<<32);
    }
  }
}

/* P[i] = XOR of Bp[k] over k with A1[i][k] set.  NY = n/256 ymm per row. */
AVX2 static void corr_bits(int n){
  int W=n/64, NY=W>>2;
  for (int i=0;i<n;i++){
    __m256i acc[16];
    for(int t=0;t<NY;t++) acc[t]=_mm256_setzero_si256();
    const uint64_t*row=A1b+(size_t)i*W;
    for(int w=0;w<W;w++){
      uint64_t b=row[w];
      while(b){
        int t=__builtin_ctzll(b); b&=b-1;
        const uint64_t*br=Bpb+(size_t)(w*64+t)*W;
        for(int q=0;q<NY;q++)
          acc[q]=_mm256_xor_si256(acc[q],_mm256_load_si256((const __m256i*)(br+4*q)));
      }
    }
    uint64_t*o=Pbb+(size_t)i*W;
    for(int q=0;q<NY;q++) _mm256_store_si256((__m256i*)(o+4*q),acc[q]);
  }
}

static uint64_t A1T[NMAX*NMAX/8] __attribute__((aligned(64)));
AVX2 static void transpose64(uint64_t *a){
  int j,k; uint64_t m,t;
  for (j=32,m=0xFFFFFFFF00000000ULL; j; j>>=1, m^=m>>j)
    for (k=0;k<64;k=(k+j+1)&~j){
      t=(a[k]^(a[k+j]<<j))&m;
      a[k]^=t; a[k+j]^=(t>>j);
    }
}
AVX2 static void build_A1T(int n){
  int W=n/64;
  for (int kb=0; kb<n; kb+=64)
    for (int ib=0; ib<n; ib+=64){
      uint64_t blk[64];
      for (int i=0;i<64;i++) blk[i]=A1b[(size_t)(ib+i)*W+(kb>>6)];
      transpose64(blk);
      for (int k=0;k<64;k++) A1T[(size_t)(kb+k)*W+(ib>>6)]=blk[k];
    }
}

/* 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=128;
  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
    }
  }
}

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[128] __attribute__((aligned(64)));
AVX2 static void epilogue(int n, int8_t*C, const uint64_t*P, int jb, int i){
  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 uint64_t*pw=P+((size_t)i*(n>>3)+(jb>>3))/8;
  for(int r=0;r<2;r++){
    const uint64_t*qa=pw+(size_t)r*(n>>6);
    for(int t=0;t<2;t++){
      __m256i v=_mm256_load_si256((const __m256i*)(mscr+r*64+t*32));
      uint64_t w=qa[t>>1];
      __m128i e0=_mm_setzero_si128(),e1=_mm_setzero_si128();
      int sh=(t&1)?32:0;
      for(int s=0;s<4;s++){
        __m128i bb=_mm_set1_epi8((char)((w>>(sh+8*s))&0xFF));
        __m128i z=_mm_andnot_si128(_mm_cmpeq_epi8(_mm_and_si128(bb,bitm),zero),hi);
        if(s==0) e0=z; else if(s==1) e0=_mm_unpacklo_epi64(e0,z);
        else if(s==2) e1=z; else e1=_mm_unpacklo_epi64(e1,z);
      }
      __m256i msk=_mm256_set_m128i(e1,e0);
      int8_t*c=C+(size_t)(i+r)*n+jb+t*32;
      _mm256_storeu_si256((__m256i*)c,_mm256_xor_si256(v,msk));
    }
  }
}
AVX2 static void mk_s(int n, const int8_t*A0, const int8_t*BI, int8_t*C, const uint8_t*PB, int jb){ (void)PB;
  const int8_t *x0=A0, *x1=A0+2*n, *bi=BI; long cnt=n>>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"
    "vpbroadcastd (%[x0]), %%ymm0\n\t"
    "vpbroadcastd (%[x1]), %%ymm1\n\t"
    "vmovdqa 0(%[bi]), %%ymm4\n\t"
    "vmovdqa 32(%[bi]), %%ymm7\n\t"
    "vmovdqa 64(%[bi]), %%ymm2\n\t"
    "vmovdqa 96(%[bi]), %%ymm3\n\t"
    "vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
    "vmovdqa 128(%[bi]), %%ymm4\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpmaddubsw %%ymm7, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm7, %%ymm1, %%ymm6\n\t"
    "vmovdqa 160(%[bi]), %%ymm7\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpmaddubsw %%ymm2, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm2, %%ymm1, %%ymm6\n\t"
    "vmovdqa 192(%[bi]), %%ymm2\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpmaddubsw %%ymm3, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm3, %%ymm1, %%ymm6\n\t"
    "vmovdqa 224(%[bi]), %%ymm3\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpbroadcastd 4(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 4(%[x1]), %%ymm1\n\t"
    "vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
    "vmovdqa 256(%[bi]), %%ymm4\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpmaddubsw %%ymm7, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm7, %%ymm1, %%ymm6\n\t"
    "vmovdqa 288(%[bi]), %%ymm7\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpmaddubsw %%ymm2, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm2, %%ymm1, %%ymm6\n\t"
    "vmovdqa 320(%[bi]), %%ymm2\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpmaddubsw %%ymm3, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm3, %%ymm1, %%ymm6\n\t"
    "vmovdqa 352(%[bi]), %%ymm3\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpbroadcastd 8(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 8(%[x1]), %%ymm1\n\t"
    "vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
    "vmovdqa 384(%[bi]), %%ymm4\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpmaddubsw %%ymm7, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm7, %%ymm1, %%ymm6\n\t"
    "vmovdqa 416(%[bi]), %%ymm7\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpmaddubsw %%ymm2, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm2, %%ymm1, %%ymm6\n\t"
    "vmovdqa 448(%[bi]), %%ymm2\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpmaddubsw %%ymm3, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm3, %%ymm1, %%ymm6\n\t"
    "vmovdqa 480(%[bi]), %%ymm3\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpbroadcastd 12(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 12(%[x1]), %%ymm1\n\t"
    "vpmaddubsw %%ymm4, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm4, %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpmaddubsw %%ymm7, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm7, %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpmaddubsw %%ymm2, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm2, %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpmaddubsw %%ymm3, %%ymm0, %%ymm5\n\t"
    "vpmaddubsw %%ymm3, %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "addq $16, %[x0]\n\t"
    "addq $16, %[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");
}


/* ---- corr16 ported to n=2048: W=32, NY=8, 8 accumulators ---- */
#define AVX2BMI __attribute__((target("avx2,bmi")))
static __m256i XT[16][16][8] __attribute__((aligned(64)));
AVX2BMI static void corr16k(int n){
  const int W=32;
  memset(Pbb,0,(size_t)n*W*8);
  const int KB=64;
  for (int kb=0; kb<n; kb+=KB){
    int w0=kb>>6;
    const uint64_t *bw = Bpb + (size_t)w0*(size_t)W*64;      /* 64 rows, global bit w0*64+ */
    for (int g=0; g<16; g++){
      uint64_t *z=(uint64_t*)XT[g][0];
      for (int q=0;q<32;q++) z[q]=0;
      for (int t=1;t<16;t++){
        uint64_t *dst=(uint64_t*)XT[g][t];
        int lb=__builtin_ctz(t), t2=t&(t-1);
        const uint64_t *src=(const uint64_t*)XT[g][t2];
        const uint64_t *row=bw+(size_t)(4*g+lb)*W;
        if (t2) for (int q=0;q<32;q++) dst[q]=src[q]^row[q];
        else    for (int q=0;q<32;q++) dst[q]=row[q];
      }
    }
    for (int i=0;i<n;i++){
      uint64_t b0 = A1b[(size_t)i*W + w0];
      __m256i a0,a1,a2,a3,a4,a5,a6,a7;
      a0=a1=a2=a3=a4=a5=a6=a7=_mm256_setzero_si256();
      for (int g=0; g<16; g++){
        const __m256i *p = (const __m256i *)XT[g][(b0>>(4*g))&15];
        p+=8;
        a0=_mm256_xor_si256(a0,_mm256_load_si256(p-8));
        a1=_mm256_xor_si256(a1,_mm256_load_si256(p-7));
        a2=_mm256_xor_si256(a2,_mm256_load_si256(p-6));
        a3=_mm256_xor_si256(a3,_mm256_load_si256(p-5));
        a4=_mm256_xor_si256(a4,_mm256_load_si256(p-4));
        a5=_mm256_xor_si256(a5,_mm256_load_si256(p-3));
        a6=_mm256_xor_si256(a6,_mm256_load_si256(p-2));
        a7=_mm256_xor_si256(a7,_mm256_load_si256(p-1));
      }
      __m256i *o=(__m256i*)(Pbb+(size_t)i*W);
      _mm256_store_si256(o+0,_mm256_xor_si256(_mm256_load_si256(o+0),a0));
      _mm256_store_si256(o+1,_mm256_xor_si256(_mm256_load_si256(o+1),a1));
      _mm256_store_si256(o+2,_mm256_xor_si256(_mm256_load_si256(o+2),a2));
      _mm256_store_si256(o+3,_mm256_xor_si256(_mm256_load_si256(o+3),a3));
      _mm256_store_si256(o+4,_mm256_xor_si256(_mm256_load_si256(o+4),a4));
      _mm256_store_si256(o+5,_mm256_xor_si256(_mm256_load_si256(o+5),a5));
      _mm256_store_si256(o+6,_mm256_xor_si256(_mm256_load_si256(o+6),a6));
      _mm256_store_si256(o+7,_mm256_xor_si256(_mm256_load_si256(o+7),a7));
    }
  }
}

/* ---- table-driven epilogue: 4 x (movq table load) instead of the scalar
        set1/cmpeq/andnot chain; byte-identical mask, verified by gate. ---- */
static uint8_t EXPT[256*8] __attribute__((aligned(64)));
/* [m2n_ A7] SIMD bit-expansion constants (same form as mmmc1k #104944 L500-521).
   EPD_IDX[c] = c div 8  -> after vpshufb over a vpbroadcastd of the four Pbb
   bytes, byte c of the spread holds the Pbb byte that owns column c.
   EPD_BIT[c] = 1 << (c mod 8) -> isolates the one bit that decides column c.
   EPD_80     = 0x80 everywhere -> turns the vpcmpeqb result into the xor mask. */
static const uint8_t EPD_IDX[32] __attribute__((aligned(32))) = {
  0,0,0,0,0,0,0,0, 1,1,1,1,1,1,1,1, 2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3 };
static const uint8_t EPD_BIT[32] __attribute__((aligned(32))) = {
  1,2,4,8,16,32,64,128, 1,2,4,8,16,32,64,128,
  1,2,4,8,16,32,64,128, 1,2,4,8,16,32,64,128 };
static const uint8_t EPD_80[32] __attribute__((aligned(32))) = {
  0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80, 0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,
  0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80, 0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80 };

AVX2 static void epilogue2(int n, int8_t*C, const uint64_t*P, int jb, int i){
  static int ei=0;
  if(!ei){ for(int v=0;v<256;v++) for(int b=0;b<8;b++) EXPT[v*8+b]=((v>>b)&1)?(uint8_t)0x80:0; ei=1; }
  const __m256i eidx=_mm256_load_si256((const __m256i*)EPD_IDX);
  const __m256i ebit=_mm256_load_si256((const __m256i*)EPD_BIT);
  const __m256i e80 =_mm256_load_si256((const __m256i*)EPD_80);
  const uint8_t*base=(const uint8_t*)P;
  const size_t rowoff=(size_t)i*(n>>3);          /* bytes per Pbb row = n/8 */
  const size_t bo=(size_t)jb>>3;
  for(int r=0;r<2;r++){
    const uint8_t*qr=base+rowoff+(size_t)r*(n>>3);
    for(int t=0;t<2;t++){
      __m256i vv=_mm256_load_si256((const __m256i*)(mscr+r*64+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, eidx);
      __m256i z=_mm256_and_si256(s, ebit);
      z=_mm256_and_si256(_mm256_cmpeq_epi8(z, ebit), e80);
      int8_t*c=C+(size_t)(i+r)*n+jb+t*32;
      _mm256_stream_si256((__m256i*)c,_mm256_xor_si256(vv,z));
    }
  }
}

AVX2 static void epilogue(int n, int8_t*C, const uint64_t*P, int jb, int i);

/* ---- fused builders: read A once (A0b+A1b) and B once (BI3b+Bpb) ---- */
AVX2 static void build_A0A1(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;
    int8_t*dst=A0b+(size_t)i*n;
    uint64_t*w1=A1b+(size_t)i*W;
    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*)(dst+k),    _mm256_and_si256(v0,m));
      _mm256_store_si256((__m256i*)(dst+k+32), _mm256_and_si256(v1,m));
      w1[k>>6]=(uint64_t)(uint32_t)_mm256_movemask_epi8(v0)
              |((uint64_t)(uint32_t)_mm256_movemask_epi8(v1)<<32);
    }
  }
}
AVX2 static void build_BI3Bp(int n, const int8_t*B){
  const __m256i one=_mm256_set1_epi8(1);
  int ng=n/64, W=n/64;
  for (int p=0;p<n/2;p++){
    const int8_t *r0=B+(size_t)(2*p)*n, *r1=r0+n;
    uint64_t *w0p=Bpb+(size_t)(2*p)*W, *w1p=Bpb+(size_t)(2*p+1)*W;
    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));
        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; }
      }
    }
  }
}

/* forward decls */
AVX2 static void build_A0(int n,const int8_t*A);
AVX2 static void build_BI3(int n,const int8_t*B);

/* ================= STRASSEN LEVEL 1 (n = 2048) ================= */
static int8_t MMS[7][1024*1024] __attribute__((aligned(64)));

/* ★ 本发新增：一个"子块描述符" = 从 X 取块 A(r0,c0) 与 Y 取块 B(r1,c1) 按 op 组合
   （0 copy / 1 add / 2 sub，mod 256），行列都以 `nn` 为跨距。
   ★ 它替换掉的正是 [1] 里 `blk()` 写到 SA/SB 再被 `build_*` 读回来的那两步。 */
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_A0A1` 融合：一趟产出 A0b(=v&0x7F) 与 A1b(=符号位)。 */
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_BI3Bp` 融合：一趟产出 BI3b(交织对) 与 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;
    uint64_t *w0p=Bpb+(size_t)(2*p)*W, *w1p=Bpb+(size_t)(2*p+1)*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; }
      }
    }
  }
}

/* GF(2) correction for m=1024: W=16, NY=4, 4 accumulators (corr16k hardcodes W=32) */
#define AVX2BMI3 __attribute__((target("avx2,bmi")))
static __m256i XT4[16][16][4] __attribute__((aligned(64)));
AVX2BMI3 static void corr1024(void){
  const int n=1024, W=16;
  memset(Pbb,0,(size_t)n*W*8);
  const int KB=64;
  for (int kb=0; kb<n; kb+=KB){
    int w0=kb>>6;
    const uint64_t *bw = Bpb + (size_t)w0*(size_t)W*64;
    for (int g=0; g<16; g++){
      uint64_t *z=(uint64_t*)XT4[g][0];
      for (int q=0;q<16;q++) z[q]=0;
      for (int t=1;t<16;t++){
        uint64_t *dst=(uint64_t*)XT4[g][t];
        int lb=__builtin_ctz(t), t2=t&(t-1);
        const uint64_t *src=(const uint64_t*)XT4[g][t2];
        const uint64_t *row=bw+(size_t)(4*g+lb)*W;
        if (t2) for (int q=0;q<16;q++) dst[q]=src[q]^row[q];
        else    for (int q=0;q<16;q++) dst[q]=row[q];
      }
    }
    for (int i=0;i<n;i++){
      uint64_t b0 = A1b[(size_t)i*W + w0];
      __m256i a0,a1,a2,a3;
      a0=a1=a2=a3=_mm256_setzero_si256();
      for (int g=0; g<16; g++){
        const __m256i *p = (const __m256i *)XT4[g][(b0>>(4*g))&15];
        p+=4;
        a0=_mm256_xor_si256(a0,_mm256_load_si256(p-4));
        a1=_mm256_xor_si256(a1,_mm256_load_si256(p-3));
        a2=_mm256_xor_si256(a2,_mm256_load_si256(p-2));
        a3=_mm256_xor_si256(a3,_mm256_load_si256(p-1));
      }
      __m256i *o=(__m256i*)(Pbb+(size_t)i*W);
      _mm256_store_si256(o+0,_mm256_xor_si256(_mm256_load_si256(o+0),a0));
      _mm256_store_si256(o+1,_mm256_xor_si256(_mm256_load_si256(o+1),a1));
      _mm256_store_si256(o+2,_mm256_xor_si256(_mm256_load_si256(o+2),a2));
      _mm256_store_si256(o+3,_mm256_xor_si256(_mm256_load_si256(o+3),a3));
    }
  }
}

/* one 1024x1024 sub-product through the SHIPPED pipeline, result to Cp.
   ★ 本发：入口从 (SA,SB) 改成两个子块描述符，面板由 `bld_A`/`bld_B` 直接建出。 */
AVX2 static void submm1024(const SBLK*ba, const SBLK*bb, int nn, int8_t*Cp){
  const int m=1024;
  bld_A(m,nn,ba); bld_B(m,nn,bb); corr1024();
  for (int g=0, jbg=0; g<m; g+=64, jbg=g)
    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);
      epilogue2(m, Cp, (const uint64_t*)Pbb, g, i);
    }
}

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 */
};


/* Cblk += sum_k s_k * M_k   (k over the up-to-4 pointers; NULL skips a term) */
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;
    const int8_t *p1 = M1 ? M1+(size_t)i*m : 0;
    const int8_t *p2 = M2 ? M2+(size_t)i*m : 0;
    const int8_t *p3 = M3 ? M3+(size_t)i*m : 0;
    for (int j=0;j<m;j+=32){
      __m256i v;
      if (p0) v = s0>0 ? _mm256_load_si256((const __m256i*)(p0+j)) : _mm256_sub_epi8(_mm256_setzero_si256(),_mm256_load_si256((const __m256i*)(p0+j)));
      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);
    }
  }
}

/* ★ 本发：把 [1] 的 switch 逐条改写成"子块对"表（同一组 (r0,c0,r1,c1,op)，单位 = o = m）。
   第 i 行与 [1] 的 case i 一一对应（右侧注释保留原文的分块名）。 */
/* ================= mc4_: TWO-LEVEL STRASSEN, (c) 4-TERM SIGNED-COMPOSE PACKER ===========
   Leaf (k1,k2) of the second level needs, on the A side,
       S2A = (A1 op1 A2) op2 (A3 op1 A4)
   where op1 is the LEVEL-1 case k1 and op2 the LEVEL-2 case k2.  Writing s1/s2 for the
   op -> sign map (op 0 -> 0, 1 -> +1, 2 -> -1),
       S2A = A1 + s1*A2 + s2*A3 + s2*s1*A4
   so the packer reads A's four blocks DIRECTLY.  Nothing is materialised: no SA/SB and
   no S2A/S2B (BRIEF 2.18.590, DIF -- never write it out and read it back).
   h = leaf size (512), n = A's row stride (2048), o1 = level-1 block size (1024). */
static __m256i XT4s[16][16][2] __attribute__((aligned(64)));   /* corr table for n = 512 */

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 = 1024;
  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*)(p1+k+32+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);} }
      {
        __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 = 1024;
  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=Bpb+(size_t)(2*p)*W, *w1p=Bpb+(size_t)(2*p+1)*W;
    for (int g=0;g<ng;g++){
      int8_t *o=BI3b+((size_t)g*((h/2)*128+512)+p*128);
      for (int t=0;t<64;t+=32){
        _mm_prefetch((const char*)(a1+64*g+t+4096), _MM_HINT_T0);
        _mm_prefetch((const char*)(b1+64*g+t+4096), _MM_HINT_T0);
        _mm_prefetch((const char*)(a2+64*g+t+4096), _MM_HINT_T0);
        _mm_prefetch((const char*)(a3+64*g+t+4096), _MM_HINT_T0);
        _mm_prefetch((const char*)(a4+64*g+t+4096), _MM_HINT_T0);
        _mm_prefetch((const char*)(b2+64*g+t+4096), _MM_HINT_T0);
        _mm_prefetch((const char*)(b3+64*g+t+4096), _MM_HINT_T0);
        _mm_prefetch((const char*)(b4+64*g+t+4096), _MM_HINT_T0);
        __m256i va=_mm256_load_si256((const __m256i*)(a1+64*g+t));
        __m256i vb=_mm256_load_si256((const __m256i*)(b1+64*g+t));
        if(s1){ __m256i ua=_mm256_load_si256((const __m256i*)(a2+64*g+t));
                __m256i ub=_mm256_load_si256((const __m256i*)(b2+64*g+t));
                if(s1>0){va=_mm256_add_epi8(va,ua);vb=_mm256_add_epi8(vb,ub);}
                else    {va=_mm256_sub_epi8(va,ua);vb=_mm256_sub_epi8(vb,ub);} }
        if(s2){ __m256i ua=_mm256_load_si256((const __m256i*)(a3+64*g+t));
                __m256i ub=_mm256_load_si256((const __m256i*)(b3+64*g+t));
                if(s2>0){va=_mm256_add_epi8(va,ua);vb=_mm256_add_epi8(vb,ub);}
                else    {va=_mm256_sub_epi8(va,ua);vb=_mm256_sub_epi8(vb,ub);} }
        if(s12){__m256i ua=_mm256_load_si256((const __m256i*)(a4+64*g+t));
                __m256i ub=_mm256_load_si256((const __m256i*)(b4+64*g+t));
                if(s12>0){va=_mm256_add_epi8(va,ua);vb=_mm256_add_epi8(vb,ub);}
                else     {va=_mm256_sub_epi8(va,ua);vb=_mm256_sub_epi8(vb,ub);} }
        __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; }
      }
    }
  }
}

/* corr for the n = 512 leaf: same nibble-table scheme as corr1024, W = 8. */
AVX2BMI3 static void corr512(void){
  const int n=512, W=8;
  /* First 64-column group directly initializes each Pbb row. */
  const int KB=64;
  for (int kb=0; kb<n; kb+=KB){
    int w0=kb>>6;
    const uint64_t *bw = Bpb + (size_t)w0*(size_t)W*64;
    for (int g=0; g<16; g++){
      uint64_t *z=(uint64_t*)XT4s[g][0];
      for (int q=0;q<8;q++) z[q]=0;
      for (int t=1;t<16;t++){
        uint64_t *dst=(uint64_t*)XT4s[g][t];
        int lb=__builtin_ctz(t), t2=t&(t-1);
        const uint64_t *src=(const uint64_t*)XT4s[g][t2];
        const uint64_t *row=bw+(size_t)(4*g+lb)*W;
        if (t2) for (int q=0;q<8;q++) dst[q]=src[q]^row[q];
        else    for (int q=0;q<8;q++) dst[q]=row[q];
      }
    }
    for (int i=0;i<n;i++){
      uint64_t b0 = A1b[(size_t)i*W + w0];
      __m256i a0,a1;
      a0=a1=_mm256_setzero_si256();
      for (int g=0; g<16; g++){
        const __m256i *p = (const __m256i *)XT4s[g][(b0>>(4*g))&15];
        a0=_mm256_xor_si256(a0,_mm256_load_si256(p));
        a1=_mm256_xor_si256(a1,_mm256_load_si256(p+1));
      }
      __m256i *o=(__m256i*)(Pbb+(size_t)i*W);
      if (kb == 0) {
        _mm256_store_si256(o+0,a0);
        _mm256_store_si256(o+1,a1);
      } else {
        _mm256_store_si256(o+0,_mm256_xor_si256(_mm256_load_si256(o+0),a0));
        _mm256_store_si256(o+1,_mm256_xor_si256(_mm256_load_si256(o+1),a1));
      }
    }
  }
}

/* level-2 scratch: 7 leaves of h x h = 512 x 512 */
static int8_t MMS2[7][512*512+128] __attribute__((aligned(64)));

/* one h x h leaf through the SHIPPED pipeline, result NT-stored to Cp */
AVX2 static void submm8(int h, int n, const int8_t*A, const int8_t*B, int k1, int k2, int8_t*Cp){
  buildA4(h,n,A,k1,k2); buildB4(h,n,B,k1,k2); corr512();
  for (int g=0; g<h; g+=64)
    for (int i=0; i<h; i+=2){
      mk_s(h, A0b+2*(size_t)i*h, BI3b+(size_t)(g>>6)*((h/2)*128+512), Cp, PBb, g);
      epilogue2(h, Cp, (const uint64_t*)Pbb, g, i);
    }
}

AVX2 static void strassen1(int n, const int8_t*A, const int8_t*B, int8_t*C){
  const int m = n>>1;                       /* 1024 */
  const int o = m;
  const int h = m>>1;                       /* 512: the level-2 leaf */
  for (int k1=0;k1<7;k1++){
    for (int k2=0;k2<7;k2++)
      submm8(h,n,A,B,k1,k2,MMS2[k2]);
    _mm_sfence();                       /* MMS2[k] were written with NT stores */
    cmb(MMS[k1],m,0,0,h, MMS2[0],1,MMS2[3],1,MMS2[4],-1,MMS2[6],1);
    cmb(MMS[k1],m,0,h,h, MMS2[2],1,MMS2[4],1,0,0,0,0);
    cmb(MMS[k1],m,h,0,h, MMS2[1],1,MMS2[3],1,0,0,0,0);
    cmb(MMS[k1],m,h,h,h, MMS2[0],1,MMS2[2],1,MMS2[1],-1,MMS2[5],1);
    // [c2q noB] sfence is CUMULATIVE: one fence orders every prior NT store.  MMS[k1] has no
    // reader until the final cmb(C,...) block right below, so the 6 intermediate fences are
    // removable; keeping only the last one restores the NT-drain overlap with the next k1's
    // submm8 work.  Single variable, ordering only => output bit-identical.
    if (k1 == 6) _mm_sfence();
  }
  /* 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 static void strassen1(int n, const int8_t*A, const int8_t*B, int8_t*C);

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==2048) { strassen1(n,A,B,C); return; }
  /* x4: A1T is WRITTEN by build_A1T and NEVER READ -- dead NMAX*NMAX/8 bytes + its transpose pass. */
  build_A0A1(n,A); build_BI3Bp(n,B);
  if(n==4096) corr16(); else if(n==2048) corr16k(n); 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);
        epilogue2(n, C, (const uint64_t*)Pbb, g, i); }
}

评测结果