代码

#include <immintrin.h>
#include <stdint.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#pragma GCC target("avx2,popcnt")
#ifndef COLUMN_BITS
#define COLUMN_BITS 512
#endif
static void* alloc_aligned(size_t n) {void*p=nullptr;if(posix_memalign(&p,4096,n))abort();return p;}
static inline uint32_t pack32(const bool*p) {
    const __m256i v=_mm256_loadu_si256((const __m256i*)p);
    return (uint32_t)_mm256_movemask_epi8(_mm256_cmpgt_epi8(v,_mm256_setzero_si256()));
}
// Classical Four Russians multiplication over GF(2), using groups of eight
// inner indices. A's masks are stored with the group index outermost; B is
// bit-packed. One column tile and its 256 XOR combinations fit in cache.
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;
    }
    constexpr size_t TW=COLUMN_BITS/64;
    const size_t words=(size_t)n/64,groups=(size_t)n/8;
    unsigned char*am=(unsigned char*)alloc_aligned((size_t)n*groups);
    uint64_t*bp=(uint64_t*)alloc_aligned((size_t)n*words*8);
    uint64_t*ct=(uint64_t*)alloc_aligned((size_t)n*TW*8);
    uint64_t*table=(uint64_t*)alloc_aligned(256*TW*8);
    // Each 256-column slab has 32 independently contiguous destination streams.
    for(int k0=0;k0<n;k0+=256)for(int i=0;i<n;++i) {
        const bool*a=A+(size_t)i*n+k0;
        for(int v=0;v<8;++v) {
            uint32_t p=pack32(a+32*v);
            for(int b=0;b<4;++b)am[((size_t)k0/8+4*v+b)*n+i]=(unsigned char)(p>>(8*b));
        }
    }
    for(size_t k=0;k<(size_t)n;++k)for(size_t j=0;j<words;++j)
        bp[k*words+j]=(uint64_t)pack32(B+k*n+j*64)|((uint64_t)pack32(B+k*n+j*64+32)<<32);
    uint64_t unpack[256];for(unsigned m=0;m<256;++m){uint64_t x=0;for(unsigned b=0;b<8;++b)x|=(uint64_t)((m>>b)&1)<<(8*b);unpack[m]=x;}
    for(size_t j0=0;j0<words;j0+=TW) {
        const size_t width=words-j0<TW?words-j0:TW;
        memset(ct,0,(size_t)n*TW*8);
        for(size_t g=0;g<groups;++g) {
            memset(table,0,TW*8);
            for(unsigned m=1;m<256;++m) {
                unsigned bit=(unsigned)__builtin_ctz(m),prev=m&(m-1);
                const uint64_t*x=table+(size_t)prev*TW,*y=bp+(g*8+bit)*words+j0;
                uint64_t*z=table+(size_t)m*TW;
                for(size_t w=0;w<width;++w)z[w]=x[w]^y[w];
            }
            const unsigned char*ms=am+g*n;
            for(int i=0;i<n;++i) {
                const uint64_t*x=table+(size_t)ms[i]*TW;
                uint64_t*y=ct+(size_t)i*TW;
                for(size_t w=0;w<width;w+=4) {
                    __m256i a=_mm256_loadu_si256((const __m256i*)(x+w));
                    __m256i b=_mm256_loadu_si256((const __m256i*)(y+w));
                    _mm256_storeu_si256((__m256i*)(y+w),_mm256_xor_si256(a,b));
                }
            }
        }
        for(int i=0;i<n;++i) {
            unsigned char*d=(unsigned char*)(C+(size_t)i*n+j0*64);
            const unsigned char*s=(const unsigned char*)(ct+(size_t)i*TW);
            for(size_t b=0;b<width*8;++b){uint64_t x=unpack[s[b]];memcpy(d+8*b,&x,8);}
        }
    }
    free(table);free(ct);free(bp);free(am);
}

// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent public performance input only; no contribution generator,
// production seeds, judge implementation, or fixed expected matrices.
#include <stdio.h>
#include <time.h>
static uint64_t random_state=0x65c13e729a4bd081ULL;
static uint64_t public_rng(){random_state^=random_state<<13;random_state^=random_state>>7;random_state^=random_state<<17;return random_state;}
static double seconds(){timespec t;clock_gettime(CLOCK_MONOTONIC,&t);return t.tv_sec+t.tv_nsec*1e-9;}
static void fill_public(bool*p,size_t z){uint64_t lu[256];for(unsigned m=0;m<256;++m){uint64_t q=0;for(unsigned b=0;b<8;++b)q|=(uint64_t)((m>>b)&1)<<(8*b);lu[m]=q;}for(size_t k=0;k<z;k+=64){uint64_t q=public_rng();for(unsigned b=0;b<8;++b){uint64_t x=lu[(q>>(8*b))&255];memcpy(p+k+8*b,&x,8);}}}
static void apply8(const bool*M,int n,const uint64_t*V,uint64_t*O){size_t ws=(size_t)n/64;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=(uint64_t)pack32(M+(size_t)i*n+w*64)|((uint64_t)pack32(M+(size_t)i*n+w*64+32)<<32);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=16384;size_t z=(size_t)n*n,ws=n/64;bool*A=(bool*)alloc_aligned(z),*B=(bool*)alloc_aligned(z),*C=(bool*)alloc_aligned(z);double t0=seconds();fill_public(A,z);fill_public(B,z);memset(C,0,z);double t1=seconds();matrix_multiply(n,A,B,C);double t2=seconds();uint64_t*R=(uint64_t*)alloc_aligned(8*ws*8),*BR=(uint64_t*)alloc_aligned(8*ws*8),*ABR=(uint64_t*)alloc_aligned(8*ws*8),*CR=(uint64_t*)alloc_aligned(8*ws*8);for(size_t p=0;p<8*ws;++p)R[p]=public_rng();apply8(B,n,R,BR);apply8(A,n,BR,ABR);apply8(C,n,R,CR);if(memcmp(ABR,CR,8*ws*8)){puts("Freivalds FAIL");return 1;}double t3=seconds();printf("n=%d generation_s=%.9f solver_s=%.9f eight_checks_s=%.9f total_s=%.9f checks=PASS\n",n,t1-t0,t2-t1,t3-t2,t3-t0);return 0;}

评测结果