代码

#define RB_ROWS 1024
#define RB_BITS 256
#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#pragma GCC target("avx2,popcnt")
#ifndef RB_ROWS
#define RB_ROWS 1024
#endif
#ifndef RB_BITS
#define RB_BITS 256
#endif
static void* alloc_aligned(size_t z){void*p=nullptr;if(posix_memalign(&p,4096,z))abort();return p;}
// Row-major packed bits, low bit first within each byte. C is output only.
void matrix_multiply(int n,const uint8_t*A,const uint8_t*B,uint8_t*C){
 if(n<=0)return;
 if(n%256){size_t rs=((size_t)n+7)/8;memset(C,0,(size_t)n*rs);for(int i=0;i<n;++i)for(int j=0;j<n;++j){unsigned v=0;for(int k=0;k<n;++k)v^=((A[(size_t)i*rs+k/8]>>(k%8))&1)*((B[(size_t)k*rs+j/8]>>(j%8))&1);C[(size_t)i*rs+j/8]|=v<<(j%8);}return;}
 constexpr size_t TW=RB_BITS/64,MC=RB_ROWS;
 const size_t words=(size_t)n/64,groups=(size_t)n/8;
 auto*am=(uint8_t*)alloc_aligned((size_t)n*groups);
 auto*ct=(uint64_t*)alloc_aligned(MC*TW*8);
 auto*tb=(uint64_t*)alloc_aligned(groups*256*TW*8);
 for(size_t g0=0;g0<groups;g0+=32)for(size_t i=0;i<(size_t)n;++i){const uint8_t*x=A+i*groups+g0;for(size_t v=0;v<32;++v)am[(g0+v)*n+i]=x[v];}
 const auto*bp=(const uint64_t*)B;
 for(size_t j0=0;j0<words;j0+=TW){
  size_t width=words-j0<TW?words-j0:TW;
  for(size_t g=0;g<groups;++g){auto*t=tb+g*256*TW;memset(t,0,TW*8);
   for(unsigned m=1;m<256;++m){unsigned bit=(unsigned)__builtin_ctz(m),prev=m&(m-1);auto*x=t+(size_t)prev*TW;auto*y=bp+(g*8+bit)*words+j0;auto*z=t+(size_t)m*TW;
    for(size_t w=0;w<width;w+=4){auto a=_mm256_load_si256((const __m256i*)(x+w));auto b=_mm256_loadu_si256((const __m256i*)(y+w));_mm256_store_si256((__m256i*)(z+w),_mm256_xor_si256(a,b));}
   }
  }
  for(size_t i0=0;i0<(size_t)n;i0+=MC){
   size_t rows=(size_t)n-i0<MC?(size_t)n-i0:MC;memset(ct,0,rows*TW*8);
   for(size_t g=0;g<groups;g+=4){
    const uint8_t*ms=am+g*n+i0;auto*t=tb+g*256*TW;
    for(size_t i=0;i<rows;++i){auto*x0=t+(size_t)ms[i]*TW;auto*x1=t+(256+(size_t)ms[n+i])*TW;auto*x2=t+(512+(size_t)ms[2*n+i])*TW;auto*x3=t+(768+(size_t)ms[3*n+i])*TW;auto*y=ct+i*TW;
     for(size_t w=0;w<width;w+=4){auto a=_mm256_load_si256((const __m256i*)(x0+w));auto b=_mm256_load_si256((const __m256i*)(x1+w));auto c=_mm256_load_si256((const __m256i*)(x2+w));auto d=_mm256_load_si256((const __m256i*)(x3+w));auto e=_mm256_load_si256((const __m256i*)(y+w));_mm256_store_si256((__m256i*)(y+w),_mm256_xor_si256(e,_mm256_xor_si256(_mm256_xor_si256(a,b),_mm256_xor_si256(c,d))));}
    }
   }
  for(size_t i=0;i<rows;++i)memcpy(C+(i0+i)*groups+j0*8,ct+i*TW,width*8);
 }}
 free(tb);free(ct);free(am);
}
// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent public surrogate, unrelated to any problem generator or seed.
#include <stdio.h>
#include <x86intrin.h>
static uint64_t public_state=0xe91f9d434b571ac9ULL;
static uint64_t public_rng(){public_state^=public_state<<7;public_state^=public_state>>9;return public_state^=public_state<<8;}
static void fill(uint8_t*p,size_t z){auto*q=(uint64_t*)p;for(size_t i=0;i<z/8;++i)q[i]=public_rng();}
static void apply(const uint8_t*M,int n,const uint64_t*V,uint64_t*O){size_t ws=(size_t)n/64;memset(O,0,8*ws*8);const uint64_t*P=(const uint64_t*)M;for(int i=0;i<n;++i){uint64_t a[8]={};for(size_t w=0;w<ws;++w){uint64_t q=P[(size_t)i*ws+w];for(unsigned r=0;r<8;++r)a[r]^=q&V[r*ws+w];}for(unsigned r=0;r<8;++r)O[r*ws+(size_t)i/64]|=(uint64_t)__builtin_parityll(a[r])<<(i&63);}}
int main(){constexpr int n=32768;size_t z=(size_t)n*n/8,ws=n/64;auto*A=(uint8_t*)alloc_aligned(z);auto*B=(uint8_t*)alloc_aligned(z);auto*C=(uint8_t*)alloc_aligned(z);fill(A,z);fill(B,z);memset(C,0xa5,z);_mm_lfence();auto t0=__rdtsc();matrix_multiply(n,A,B,C);_mm_lfence();auto t1=__rdtsc();auto*R=(uint64_t*)alloc_aligned(8*ws*8);auto*BR=(uint64_t*)alloc_aligned(8*ws*8);auto*ABR=(uint64_t*)alloc_aligned(8*ws*8);auto*CR=(uint64_t*)alloc_aligned(8*ws*8);for(size_t p=0;p<8*ws;++p)R[p]=public_rng();apply(B,n,R,BR);apply(A,n,BR,ABR);apply(C,n,R,CR);if(memcmp(ABR,CR,8*ws*8)){puts("eight-vector check FAIL");return 1;}uint64_t h[4]={0x19a22105f06346b5ULL,0x5fc28e4a95303ce5ULL,0x82723f1ae6b397cbULL,0x414e19a615c7aa21ULL};auto*P=(const uint64_t*)C;for(size_t p=0;p<z/8;p+=4)for(size_t v=0;v<4;++v){uint64_t x=h[v]^P[p+v];h[v]=((x<<13)|(x>>51))*0x9e3779b185ebca87ULL;}printf("n=%d MC=%d NC=%d compute_ticks=%llu all_hash=%016llx,%016llx,%016llx,%016llx check=PASS\n",n,RB_ROWS,RB_BITS,(unsigned long long)(t1-t0),(unsigned long long)h[0],(unsigned long long)h[1],(unsigned long long)h[2],(unsigned long long)h[3]);}

评测结果