代码

#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#pragma GCC optimize("O3")
#pragma GCC target("avx2,popcnt")
static void* aligned_allocated(size_t bytes){void*p=nullptr;if(posix_memalign(&p,4096,bytes))abort();return p;}
// Interface arrays use one byte per bool. Bit packing is internal workspace only.
static void make_expand(uint64_t*expand){for(unsigned q=0;q<256;++q){uint64_t v=0;for(unsigned b=0;b<8;++b)v|=(uint64_t)((q>>b)&1)<<(b*8);expand[q]=v;}}
static __attribute__((noinline,noclone)) void prepare_A(int n,const uint8_t*A,uint32_t*am){
 const size_t groups=(size_t)n/32,pitch=(size_t)n+16;
 const __m256i z=_mm256_setzero_si256();
 for(size_t g0=0;g0<groups;g0+=16)for(int i0=0;i0<n;i0+=32){
  size_t ge=groups-g0<16?groups-g0:16;int ie=n-i0<32?n-i0:32;
  for(int i=0;i<ie;++i)for(size_t g=0;g<ge;++g){
   __m256i v=_mm256_loadu_si256((const __m256i*)(A+(size_t)(i0+i)*n+(g0+g)*32));
   uint32_t q=(uint32_t)_mm256_movemask_epi8(_mm256_cmpgt_epi8(v,z));
   am[(g0+g)*pitch+i0+i]=q;
  }
 }
}
static __attribute__((noinline,noclone)) void prepare_B(int n,const uint8_t*B,uint8_t*bp){
 const __m256i z=_mm256_setzero_si256();
 for(size_t i=0;i<(size_t)n*n;i+=32){
  __m256i v=_mm256_loadu_si256((const __m256i*)(B+i));
  uint32_t q=(uint32_t)_mm256_movemask_epi8(_mm256_cmpgt_epi8(v,z));
  memcpy(bp+i/8,&q,4);
 }
}
static __attribute__((noinline,noclone)) void multiply_prepared(int n,const uint32_t*am,const uint8_t*bp,uint8_t*C,uint64_t*ct,uint64_t*table){
 const size_t row_bytes=(size_t)n/8,groups=(size_t)n/32,pitch=(size_t)n+16,words=(size_t)n/64;
 alignas(64) uint64_t expand[256];make_expand(expand);
 for(size_t j0=0;j0<words;j0+=4){
  memset(ct,0,(size_t)n*4*8);
  for(size_t g=0;g<groups;++g){
   for(unsigned h=0;h<4;++h){
    uint64_t*t=table+(size_t)h*256*4;
    _mm256_store_si256((__m256i*)t,_mm256_setzero_si256());
    for(unsigned m=1;m<256;++m){unsigned bit=__builtin_ctz(m),prev=m&(m-1);
     __m256i x=_mm256_load_si256((const __m256i*)(t+(size_t)prev*4));
     __m256i y=_mm256_loadu_si256((const __m256i*)(bp+(g*32+h*8+bit)*row_bytes+j0*8));
     _mm256_store_si256((__m256i*)(t+(size_t)m*4),_mm256_xor_si256(x,y));
    }
   }
   const uint32_t*ms=am+g*pitch;
   for(int i=0;i<n;++i){uint32_t q=ms[i];
    __m256i v=_mm256_load_si256((const __m256i*)(ct+(size_t)i*4));
    v=_mm256_xor_si256(v,_mm256_load_si256((const __m256i*)(table+(size_t)(q&255)*4)));
    v=_mm256_xor_si256(v,_mm256_load_si256((const __m256i*)(table+(256+(size_t)((q>>8)&255))*4)));
    v=_mm256_xor_si256(v,_mm256_load_si256((const __m256i*)(table+(512+(size_t)((q>>16)&255))*4)));
    v=_mm256_xor_si256(v,_mm256_load_si256((const __m256i*)(table+(768+(size_t)(q>>24))*4)));
    _mm256_store_si256((__m256i*)(ct+(size_t)i*4),v);
   }
  }
  for(int i=0;i<n;++i){
   const uint8_t*packed=(const uint8_t*)(ct+(size_t)i*4);
   uint8_t*out=C+(size_t)i*n+j0*64;
   for(unsigned k=0;k<32;k+=4){
    __m256i q=_mm256_set_epi64x(expand[packed[k+3]],expand[packed[k+2]],expand[packed[k+1]],expand[packed[k]]);
    _mm256_stream_si256((__m256i*)(out+k*8),q);
   }
  }
 }
 _mm_sfence();
}
__attribute__((noinline,noclone)) 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 v=false;for(int k=0;k<n;++k)v^=A[(size_t)i*n+k]&B[(size_t)k*n+j];C[(size_t)i*n+j]=v;}
  return;
 }
 const size_t groups=(size_t)n/32,pitch=(size_t)n+16;
 uint32_t*am=(uint32_t*)aligned_allocated(groups*pitch*4);
 uint8_t*bp=(uint8_t*)aligned_allocated((size_t)n*n/8);
 uint64_t*ct=(uint64_t*)aligned_allocated((size_t)n*4*8);
 uint64_t*table=(uint64_t*)aligned_allocated(4*256*4*8);
 prepare_A(n,(const uint8_t*)A,am);prepare_B(n,(const uint8_t*)B,bp);
 multiply_prepared(n,am,bp,(uint8_t*)C,ct,table);
 free(table);free(ct);free(bp);free(am);
}

// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent public surrogate. Three raw-array lifetimes are represented by
// one 1-GiB allocation reused across disjoint stages to honor the test hard cap.
// This measures every full-size input packing, multiplication, and output write;
// it is not a native simultaneous 3-GiB external-buffer allocation test.
#include <stdio.h>
#include <x86intrin.h>
static uint64_t public_state=0x7158e1c09b432da6ULL;
static uint64_t public_word(){uint64_t z=(public_state+=0x9e3779b97f4a7c15ULL);z=(z^(z>>30))*0xbf58476d1ce4e5b9ULL;z=(z^(z>>27))*0x94d049bb133111ebULL;return z^(z>>31);}
static __attribute__((noinline,noclone)) void public_fill_raw(uint8_t*p,size_t z){
 alignas(64)uint64_t expand[256];make_expand(expand);
 for(size_t i=0;i<z;i+=64){uint64_t q=public_word();for(unsigned b=0;b<8;++b){uint64_t x=expand[(q>>(b*8))&255];memcpy(p+i+b*8,&x,8);}}
}
static void apply8B(const uint8_t*M,int n,const uint64_t*V,uint64_t*O){size_t ws=(size_t)n/64,rb=(size_t)n/8;memset(O,0,8*ws*8);for(int i=0;i<n;++i){uint64_t acc[8]={};for(size_t w=0;w<ws;++w){uint64_t q;memcpy(&q,M+(size_t)i*rb+w*8,8);for(unsigned r=0;r<8;++r)acc[r]^=q&V[r*ws+w];}for(unsigned r=0;r<8;++r)O[r*ws+(size_t)i/64]|=(uint64_t)__builtin_parityll(acc[r])<<(i&63);}}
static void apply8A(const uint32_t*am,int n,const uint64_t*V,uint64_t*O){
 size_t ws=(size_t)n/64,pitch=(size_t)n+16;uint32_t*acc=(uint32_t*)calloc((size_t)8*n,4);if(!acc)abort();
 for(size_t g=0;g<(size_t)n/32;++g){uint32_t masks[8];for(unsigned r=0;r<8;++r)masks[r]=(uint32_t)(V[r*ws+g/2]>>((g&1)*32));
  for(int i=0;i<n;++i){uint32_t q=am[g*pitch+i];for(unsigned r=0;r<8;++r)acc[(size_t)r*n+i]^=q&masks[r];}}
 memset(O,0,8*ws*8);for(int i=0;i<n;++i)for(unsigned r=0;r<8;++r)O[r*ws+(size_t)i/64]|=(uint64_t)__builtin_parity(acc[(size_t)r*n+i])<<(i&63);free(acc);
}
static uint64_t apply8C(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);uint64_t hash=0xcbf29ce484222325ULL;
 const __m256i z=_mm256_setzero_si256(),ones=_mm256_set1_epi8(1);
 for(int i=0;i<n;++i){uint64_t acc[8]={};for(size_t w=0;w<ws;++w){
  __m256i x0=_mm256_load_si256((const __m256i*)(M+(size_t)i*n+w*64));
  __m256i x1=_mm256_load_si256((const __m256i*)(M+(size_t)i*n+w*64+32));
  if(_mm256_movemask_epi8(_mm256_cmpgt_epi8(x0,ones))||_mm256_movemask_epi8(_mm256_cmpgt_epi8(x1,ones))||_mm256_movemask_epi8(_mm256_cmpgt_epi8(z,x0))||_mm256_movemask_epi8(_mm256_cmpgt_epi8(z,x1)))abort();
  uint64_t q=(uint32_t)_mm256_movemask_epi8(_mm256_cmpgt_epi8(x0,z));
  q|=(uint64_t)(uint32_t)_mm256_movemask_epi8(_mm256_cmpgt_epi8(x1,z))<<32;
  hash=(hash^q)*0x100000001b3ULL;
  for(unsigned r=0;r<8;++r)acc[r]^=q&V[r*ws+w];}
  for(unsigned r=0;r<8;++r)O[r*ws+(size_t)i/64]|=(uint64_t)__builtin_parityll(acc[r])<<(i&63);
 }return hash;
}
int main(){constexpr int n=32768;constexpr size_t rawbytes=(size_t)n*n,ws=n/64;
 _mm_lfence();unsigned long long t0=__rdtsc();
 uint8_t*raw=(uint8_t*)aligned_allocated(rawbytes);
 uint32_t*am=(uint32_t*)aligned_allocated((size_t)n/32*(n+16)*4);
 uint8_t*bp=(uint8_t*)aligned_allocated(rawbytes/8);
 uint64_t*ct=(uint64_t*)aligned_allocated((size_t)n*4*8);
 uint64_t*table=(uint64_t*)aligned_allocated(4*256*4*8);
 public_fill_raw(raw,rawbytes);
 _mm_lfence();unsigned long long t1=__rdtsc();
 prepare_A(n,raw,am);
 _mm_lfence();unsigned long long t2=__rdtsc();
 public_fill_raw(raw,rawbytes);
 _mm_lfence();unsigned long long t3=__rdtsc();
 prepare_B(n,raw,bp);
 _mm_lfence();unsigned long long t4=__rdtsc();
 memset(raw,0,rawbytes);
 _mm_lfence();unsigned long long t5=__rdtsc();
 multiply_prepared(n,am,bp,raw,ct,table);
 _mm_lfence();unsigned long long t6=__rdtsc();
 uint64_t*R=(uint64_t*)aligned_allocated(8*ws*8),*BR=(uint64_t*)aligned_allocated(8*ws*8),*ABR=(uint64_t*)aligned_allocated(8*ws*8),*CR=(uint64_t*)aligned_allocated(8*ws*8);
 for(size_t i=0;i<8*ws;++i)R[i]=public_word();
 apply8B(bp,n,R,BR);apply8A(am,n,BR,ABR);uint64_t hash=apply8C(raw,n,R,CR);
 if(memcmp(ABR,CR,8*ws*8)){puts("raw GF2 checks FAILED");return 2;}
 _mm_lfence();unsigned long long t7=__rdtsc();
 printf("n=%d genA=%llu packA=%llu genB=%llu packB=%llu zeroC=%llu mul_expand=%llu checks=%llu hash=%016llx vectors=PASS\n",n,t1-t0,t2-t1,t3-t2,t4-t3,t5-t4,t6-t5,t7-t6,(unsigned long long)hash);
 return 0;
}

评测结果