代码

// ===== REFERENCES =====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**，提交 **#125832** <https://duck.ac/submission/125832>
//     用途：**直接复制了该提交的正文**（本发正文 = 该公开件的逐字节副本 + 一处改动）。其自带的
//     References/Idea 两块原样保留在下方（其来源为 saffah_codex_maker 的 #125707 八行子集表算法）。
//     许可证：以上均为 duck.ac 公开提交与公开题面，未见任何许可证声明；按 RULES §3 逐项署名引用。
// ======================
// ===== 思路 =====
// 【本发 = 对手 #125832 的底盘 + **ARM IXT**（把 `ix[]` 的构造成"顺序写 + 分块转置"）】
//
// 【取证（本席本轮先做的相位账，n=2048，判题机同码的分段消融，本机 min-of-2）】
//   底盘 14.16 ms（本机）分段：**删除法消融** —— 跳过 `pc[i]^=tab[x[i]]` 内环 ⇒ 8.23（该刀含被它喂活的表构建）；
//   跳过表构建的 247 条 XOR ⇒ 13.57（**表构建 ≈ 0.6 ms，很小**）；跳过 `ix` 构造 ⇒ 12.01（**ix ≈ 2.15 ms，15%**）；
//   跳过 B 打包 ⇒ 13.78（≈0.4 ms）。**⇒ 真正的大头是内环与输出展开，而 `ix` 构造里最贵的是它的写模式。**
//
// 【病灶】父亲件构 `ix[]` 的循环是 `for i { for g { ix[g*n + i] = pext(A + i*n + 8g) } }`：A 的读是完美的
//   顺序流（4 MB 一次过），但 **`ix` 的每一次写都跨 n=2048 B** ⇒ 2048×256 = 524k 次"每写一行一个新 cache line"
//   的散写（ix 本体只有 512 KB，却被写成 524k 次 RFO）。
//
// 【本发做法】**不改语义、不改布局，只改写入粒度**：
//   ① 先按行序建 `ixT[i*groups + g]`：每次 `_mm256_loadu_si256` 取 32 字节 A（= 4 个 8 列组），
//      `_mm256_slli_epi16(v,7)` + `movemask` 一次得到 4 个组的 8 位索引（与父亲件打包 `pb` 用的是**同一个位序**：
//      movemask 取每个字节的 MSB，而 bool 元素只有 0/1 ⇒ bit0 抬到 bit7 后 MSB 即该元素），
//      4 次连续字节写 ⇒ **顺序写**（每行 32 次 4 字节写，全部落在同一行内）；
//   ② 再做 64×64 分块转置搬进 `ix[g*n + i]`（分块的意义：一个 64×64 瓦片里，输出对每个 g 恰是 64 个连续字节、
//      输入对每个 i 恰是 64 个连续字节 ⇒ 全部读写都是 cache-line 粒度）。
//   ★ 等价性：两处写进去的都是**同一个值**（pext 取低位的 8 个 bit == movemask 取 MSB 的 8 个 bit，因为 bool∈{0,1}），
//   只是写入顺序/粒度不同 ⇒ **输出逐位不变**。实测：n=2048 与 n=1024 上本发与底盘 **fnv 完全相同**
//   （`562047169027456410` / `5916920866438164018`），小 n（256/512）上与底盘逐字节相同。
//
// 【本机 min-of-3（只作否决）】n=2048：13.434 → **12.586 ms（−6.3%）** · n=1024：2.248 → **1.908 ms（−15.1%）**。
//   同一份代码在判题机上 n=2048 为 **T = 7.160888 ms** ⇒ 本刀按同比例应为 ≈6.7 ms。
// 【本次为正式提交（非实验）：本账号在该题此前零提交，这是第一发。】
// ================
/* References:
- saffah_codex_maker, https://duck.ac/submission/125707: the public GF(2)
  eight-row subset-table algorithm and arbitrary-size fallback. No additional
  license notice was present. The implementation below is newly written.
Idea: Pack B using AVX2 movemasks, precompute all A subset indices using BMI2
  extraction in group-major order, and process256 output columns at a time.
  Each lookup table is8KiB; accumulate four rows with independent registers.
  Expand packed output with byte shuffles rather than a scalar bit loop.
Experiment purpose: Remove scalar packing and reduce table/cache working sets.
*/
#pragma GCC optimize("O3")
#pragma GCC target("avx2,bmi2")
#include <immintrin.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
void matrix_multiply(int n,const bool*A,const bool*B,bool*C){
 if(n<=0)return;
 if(n%256){for(int i=0;i<n;i++)for(int j=0;j<n;j++){bool s=0;for(int k=0;k<n;k++)s^=A[(size_t)i*n+k]&B[(size_t)k*n+j];C[(size_t)i*n+j]=s;}return;}
 const int w=n/64,groups=n/8;
 uint64_t*pb;unsigned char*ix;__m256i*pc;
 if(posix_memalign((void**)&pb,4096,(size_t)n*w*8)||posix_memalign((void**)&ix,4096,(size_t)groups*n)||posix_memalign((void**)&pc,4096,(size_t)n*32))abort();
 for(int k=0;k<n;k++)for(int j=0;j<n;j+=64){
  __m256i x=_mm256_load_si256((const __m256i*)(B+(size_t)k*n+j));
  __m256i y=_mm256_load_si256((const __m256i*)(B+(size_t)k*n+j+32));
  uint32_t lo=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(x,7));
  uint32_t hi=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(y,7));
  pb[(size_t)k*w+j/64]=(uint64_t)lo|((uint64_t)hi<<32);
 }
 // ARM IXT: build the group indices row-major (sequential stores, 4 groups per vector
 // move-mask) and then blocked-transpose into the group-major layout the kernel consumes.
 // The parent wrote ix[g*n+i] directly from a row-major scan => 524k stores strided by n
 // (one cache line each).  Same bytes, same values, just cache-line-granular traffic.
 { unsigned char* ixT; if(posix_memalign((void**)&ixT,4096,(size_t)n*groups))abort();
   for(int i=0;i<n;i++){
     const bool* Ai=A+(size_t)i*n; unsigned char* Ti=ixT+(size_t)i*groups;
     for(int g=0;g<groups;g+=4){
       __m256i v=_mm256_loadu_si256((const __m256i*)(Ai+(size_t)g*8));
       uint32_t m=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(v,7));
       Ti[g]=(unsigned char)m; Ti[g+1]=(unsigned char)(m>>8);
       Ti[g+2]=(unsigned char)(m>>16); Ti[g+3]=(unsigned char)(m>>24);
     }
   }
   for(int ib=0;ib<n;ib+=64){ int ie=ib+64; if(ie>n) ie=n;
     for(int gb=0;gb<groups;gb+=64){ int ge=gb+64; if(ge>groups) ge=groups;
       for(int g=gb;g<ge;g++){ size_t gi=(size_t)g*n;
         for(int i=ib;i<ie;i++) ix[gi+i]=ixT[(size_t)i*groups+g]; } } }
   free(ixT);
 }
 alignas(64) __m256i tab[256];
 const __m256i shuf=_mm256_setr_epi8(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);
 const __m256i bits=_mm256_set1_epi64x(0x8040201008040201ULL),one=_mm256_set1_epi8(1),zero=_mm256_setzero_si256();
 for(int jb=0;jb<w;jb+=4){
  for(int i=0;i<n;i++)pc[i]=zero;
  for(int g=0;g<groups;g++){
   tab[0]=zero;
   for(int t=0;t<8;t++){
    __m256i row=_mm256_load_si256((const __m256i*)(pb+(size_t)(g*8+t)*w+jb));
    for(int m=0;m<(1<<t);m++)tab[m+(1<<t)]=_mm256_xor_si256(tab[m],row);
   }
   const unsigned char*x=ix+(size_t)g*n;
   for(int i=0;i<n;i+=4){
    pc[i]=_mm256_xor_si256(pc[i],tab[x[i]]);
    pc[i+1]=_mm256_xor_si256(pc[i+1],tab[x[i+1]]);
    pc[i+2]=_mm256_xor_si256(pc[i+2],tab[x[i+2]]);
    pc[i+3]=_mm256_xor_si256(pc[i+3],tab[x[i+3]]);
   }
  }
  for(int i=0;i<n;i++){
   const uint32_t*row=(const uint32_t*)(pc+i);
   for(int t=0;t<8;t++){
    __m256i x=_mm256_shuffle_epi8(_mm256_set1_epi32(row[t]),shuf);
    x=_mm256_and_si256(_mm256_cmpeq_epi8(_mm256_and_si256(x,bits),bits),one);
    _mm256_store_si256((__m256i*)(C+(size_t)i*n+jb*64+t*32),x);
   }
  }
 }
 free(pc);free(ix);free(pb);
}

评测结果