代码

/* 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);
 }
 for(int i=0;i<n;i++)for(int g=0;g<groups;g++){
  uint64_t x;memcpy(&x,A+(size_t)i*n+g*8,8);
  ix[(size_t)g*n+i]=(unsigned char)_pext_u64(x,0x0101010101010101ULL);
 }
 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);
}

评测结果