提交记录 126876


用户 题目 状态 得分 用时 内存 语言 代码长度
saffah_cc_v41_agg1 mmmc1k. 测测你的8位整数矩阵乘法-1k Accepted 100 6.669 ms 1732 KB C++17 68.94 KB
提交时间 评测时间
2026-10-15 14:53:51 2026-10-15 14:53:54
// ===== 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.
// 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 加到另外两个冷写目的地。
// 闸门:与底盘逐位一致(6 形态 × n=256/512/1024,md5 + 象限和)+ 独立标量
// oracle + 阳性对照;判题机侧自校验和亦逐位相同。
// ======================
#define MC4_KNIVES_DE 1
#include <immintrin.h>
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing,rename-registers")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
#include <stdint.h>
#include <string.h>
#include <stdio.h>

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

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

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

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

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

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


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 static void mk_s(int n, const int8_t*A0, const int8_t*BI, int8_t*C, const uint8_t*PB, int jb){ (void)PB;
  const int8_t *x0=A0, *x1=A0+2*n, *bi=BI; long cnt=n>>4;
  asm volatile(
"vpxor %%xmm8,%%xmm8,%%xmm8\n\t"  "vpxor %%xmm9,%%xmm9,%%xmm9\n\t"  "vpxor %%xmm10,%%xmm10,%%xmm10\n\t"  "vpxor %%xmm11,%%xmm11,%%xmm11\n\t"  "vpxor %%xmm12,%%xmm12,%%xmm12\n\t"  "vpxor %%xmm13,%%xmm13,%%xmm13\n\t"  "vpxor %%xmm14,%%xmm14,%%xmm14\n\t"  "vpxor %%xmm15,%%xmm15,%%xmm15\n\t"
    ".p2align 5\n\t" "1:\n\t"
    "prefetcht0 1536(%[bi])\n\t"
    "vpbroadcastd (%[x0]), %%ymm0\n\t"
    "vpbroadcastd (%[x1]), %%ymm1\n\t"
    "vpmaddubsw 0(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 0(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 32(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 32(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 64(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 64(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 96(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 96(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpbroadcastd 4(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 4(%[x1]), %%ymm1\n\t"
    "vpmaddubsw 128(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 128(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 160(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 160(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 192(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 192(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 224(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 224(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpbroadcastd 8(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 8(%[x1]), %%ymm1\n\t"
    "vpmaddubsw 256(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 256(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 288(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 288(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 320(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 320(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 352(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 352(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpbroadcastd 12(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 12(%[x1]), %%ymm1\n\t"
    "vpmaddubsw 384(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 384(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 416(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 416(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 448(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 448(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 480(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 480(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpbroadcastd 16(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 16(%[x1]), %%ymm1\n\t"
    "vpmaddubsw 512(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 512(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 544(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 544(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 576(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 576(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 608(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 608(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpbroadcastd 20(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 20(%[x1]), %%ymm1\n\t"
    "vpmaddubsw 640(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 640(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 672(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 672(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 704(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 704(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 736(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 736(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpbroadcastd 24(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 24(%[x1]), %%ymm1\n\t"
    "vpmaddubsw 768(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 768(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 800(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 800(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 832(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 832(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 864(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 864(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "vpbroadcastd 28(%[x0]), %%ymm0\n\t"
    "vpbroadcastd 28(%[x1]), %%ymm1\n\t"
    "vpmaddubsw 896(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 896(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm12, %%ymm12\n\t"
    "vpaddw %%ymm5, %%ymm8, %%ymm8\n\t"
    "vpmaddubsw 928(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 928(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm13, %%ymm13\n\t"
    "vpaddw %%ymm5, %%ymm9, %%ymm9\n\t"
    "vpmaddubsw 960(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 960(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm14, %%ymm14\n\t"
    "vpaddw %%ymm5, %%ymm10, %%ymm10\n\t"
    "vpmaddubsw 992(%[bi]), %%ymm0, %%ymm5\n\t"
    "vpmaddubsw 992(%[bi]), %%ymm1, %%ymm6\n\t"
    "vpaddw %%ymm6, %%ymm15, %%ymm15\n\t"
    "vpaddw %%ymm5, %%ymm11, %%ymm11\n\t"
    "addq $32, %[x0]\n\t"
    "addq $32, %[x1]\n\t"
    "addq $1024, %[bi]\n\t"
    "decq %[cnt]\n\t"
    "jnz 1b\n\t"

    "vpand %[mask], %%ymm8, %%ymm2\n\t"
    "vpand %[mask], %%ymm9, %%ymm3\n\t"
    "vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
    "vpermq $0xD8, %%ymm2, %%ymm2\n\t"
    "vmovdqu %%ymm2, 0(%[scr])\n\t"
    "vpand %[mask], %%ymm10, %%ymm2\n\t"
    "vpand %[mask], %%ymm11, %%ymm3\n\t"
    "vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
    "vpermq $0xD8, %%ymm2, %%ymm2\n\t"
    "vmovdqu %%ymm2, 32(%[scr])\n\t"
    "vpand %[mask], %%ymm12, %%ymm2\n\t"
    "vpand %[mask], %%ymm13, %%ymm3\n\t"
    "vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
    "vpermq $0xD8, %%ymm2, %%ymm2\n\t"
    "vmovdqu %%ymm2, 64(%[scr])\n\t"
    "vpand %[mask], %%ymm14, %%ymm2\n\t"
    "vpand %[mask], %%ymm15, %%ymm3\n\t"
    "vpackuswb %%ymm3, %%ymm2, %%ymm2\n\t"
    "vpermq $0xD8, %%ymm2, %%ymm2\n\t"
    "vmovdqu %%ymm2, 96(%[scr])\n\t"
    : [x0]"+r"(x0),[x1]"+r"(x1),[bi]"+r"(bi),[cnt]"+r"(cnt)
    : [mask]"m"(W0FF), [scr]"r"(mscr)
    : "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15","cc","memory");
}

static uint8_t Gtab[8*256*(NMAX/8)] __attribute__((aligned(64)));


#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(64) 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); }
}


CompilationN/AN/ACompile OKScore: N/A

Testcase #16.669 ms1 MB + 708 KBAcceptedScore: 100


Judge Duck Online | 评测鸭在线
Server Time: 2026-10-17 10:21:42 | Loaded in 1 ms | Server Status
个人娱乐项目,仅供学习交流使用 | 捐赠