代码

#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;}
// Packed row-major GF(2): bit (j&7) of byte (j/8) denotes column j.
// Four independent 8-bit tables are combined before storing an output tile.
__attribute__((noinline,noclone)) void matrix_multiply(int n,const uint8_t*A,const uint8_t*B,uint8_t*C){
 if(n<=0)return;
 const size_t row_bytes=((size_t)n+7)/8;
 if(n%256){
  memset(C,0,(size_t)n*row_bytes);
  for(int i=0;i<n;++i)for(int k=0;k<n;++k)if((A[(size_t)i*row_bytes+k/8]>>(k&7))&1)
   for(size_t j=0;j<row_bytes;++j)C[(size_t)i*row_bytes+j]^=B[(size_t)k*row_bytes+j];
  if(n&7)for(int i=0;i<n;++i)C[(size_t)(i+1)*row_bytes-1]&=(uint8_t)((1u<<(n&7))-1);
  return;
 }
 const size_t groups=(size_t)n/32,pitch=(size_t)n+16,words=(size_t)n/64;
 uint32_t*am=(uint32_t*)aligned_allocated(groups*pitch*4);
 uint64_t*ct=(uint64_t*)aligned_allocated((size_t)n*4*8);
 uint64_t*table=(uint64_t*)aligned_allocated(4*256*4*8);
 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){uint32_t q;memcpy(&q,A+(size_t)(i0+i)*row_bytes+(g0+g)*4,4);am[(g0+g)*pitch+i0+i]=q;}
 }
 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*)(B+(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)_mm256_storeu_si256((__m256i*)(C+(size_t)i*row_bytes+j0*8),_mm256_load_si256((const __m256i*)(ct+(size_t)i*4)));
 }
 free(table);free(ct);free(am);
}

// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent nonlinear public test data. No production seed, generator or judge ABI.
#include <stdio.h>
#include <x86intrin.h>
static uint64_t public_state=0x469a183bd72905e1ULL;
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 void public_fill(uint8_t*p,size_t z){for(size_t i=0;i<z;i+=8){uint64_t q=public_word();memcpy(p+i,&q,8);}}
static void apply8(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);}}
int main(){constexpr int n=32768;constexpr size_t z=(size_t)n*n/8,ws=n/64;
 _mm_lfence();unsigned long long t0=__rdtsc();
 uint8_t*A=(uint8_t*)aligned_allocated(z),*B=(uint8_t*)aligned_allocated(z),*C=(uint8_t*)aligned_allocated(z);
 public_fill(A,z);public_fill(B,z);memset(C,0,z);
 _mm_lfence();unsigned long long t1=__rdtsc();
 matrix_multiply(n,A,B,C);
 _mm_lfence();unsigned long long t2=__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();
 apply8(B,n,R,BR);apply8(A,n,BR,ABR);apply8(C,n,R,CR);
 if(memcmp(ABR,CR,8*ws*8)){puts("full packed GF2 checks FAILED");return 2;}
 uint64_t hash=0xcbf29ce484222325ULL;for(size_t i=0;i<z;i+=8){uint64_t q;memcpy(&q,C+i,8);hash=(hash^q)*0x100000001b3ULL;}
 _mm_lfence();unsigned long long t3=__rdtsc();
 printf("n=%d input_ticks=%llu multiply_ticks=%llu checks_ticks=%llu hash=%016llx eight_vectors=PASS\n",n,t1-t0,t2-t1,t3-t2,(unsigned long long)hash);
 return 0;
}

评测结果