代码

#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#pragma GCC optimize("O3")

namespace matrix_x {
using U=unsigned __int128;
static void classical(int n,const U*a,int sa,const U*b,int sb,U*c,int sc){
    for(int i=0;i<n;++i)memset(c+(size_t)i*sc,0,(size_t)n*sizeof(U));
    for(int i=0;i<n;++i){
        U*cr=c+(size_t)i*sc;
        for(int k=0;k<n;++k){
            U av=a[(size_t)i*sa+k];
            const U*br=b+(size_t)k*sb;
            for(int j=0;j<n;++j)cr[j]+=av*br[j];
        }
    }
}
static void combine(int n,const U*a,int sa,const U*b,int sb,U*out,int sign){
    for(int i=0;i<n;++i)for(int j=0;j<n;++j)
        out[(size_t)i*n+j]=sign>0?a[(size_t)i*sa+j]+b[(size_t)i*sb+j]:a[(size_t)i*sa+j]-b[(size_t)i*sb+j];
}
static void add(int n,U*c,int sc,const U*p,int sign){
    for(int i=0;i<n;++i)for(int j=0;j<n;++j)
        if(sign>0)c[(size_t)i*sc+j]+=p[(size_t)i*n+j];
        else c[(size_t)i*sc+j]-=p[(size_t)i*n+j];
}
static void recurse(int n,const U*a,int sa,const U*b,int sb,U*c,int sc,U*work){
    if(n<=128){classical(n,a,sa,b,sb,c,sc);return;}
    int h=n/2;size_t count=(size_t)h*h;
    U*s=work,*t=s+count,*p=t+count,*next=p+count;
    const U*a00=a,*a01=a+h,*a10=a+(size_t)h*sa,*a11=a10+h;
    const U*b00=b,*b01=b+h,*b10=b+(size_t)h*sb,*b11=b10+h;
    U*c00=c,*c01=c+h,*c10=c+(size_t)h*sc,*c11=c10+h;
    for(int i=0;i<n;++i)memset(c+(size_t)i*sc,0,(size_t)n*sizeof(U));
    combine(h,a00,sa,a11,sa,s,1);combine(h,b00,sb,b11,sb,t,1);
    recurse(h,s,h,t,h,p,h,next);add(h,c00,sc,p,1);add(h,c11,sc,p,1);
    combine(h,a10,sa,a11,sa,s,1);
    recurse(h,s,h,b00,sb,p,h,next);add(h,c10,sc,p,1);add(h,c11,sc,p,-1);
    combine(h,b01,sb,b11,sb,t,-1);
    recurse(h,a00,sa,t,h,p,h,next);add(h,c01,sc,p,1);add(h,c11,sc,p,1);
    combine(h,b10,sb,b00,sb,t,-1);
    recurse(h,a11,sa,t,h,p,h,next);add(h,c00,sc,p,1);add(h,c10,sc,p,1);
    combine(h,a00,sa,a01,sa,s,1);
    recurse(h,s,h,b11,sb,p,h,next);add(h,c00,sc,p,-1);add(h,c01,sc,p,1);
    combine(h,a10,sa,a00,sa,s,-1);combine(h,b00,sb,b01,sb,t,1);
    recurse(h,s,h,t,h,p,h,next);add(h,c11,sc,p,1);
    combine(h,a01,sa,a11,sa,s,-1);combine(h,b10,sb,b11,sb,t,1);
    recurse(h,s,h,t,h,p,h,next);add(h,c00,sc,p,1);
}
}
void matrix_multiply(int n,const __int128*A,const __int128*B,__int128*C){
    using namespace matrix_x;
    if(n<=128 || (n&(n-1))) {
        classical(n,(const U*)A,n,(const U*)B,n,(U*)C,n);
        return;
    }
    U*work=nullptr;
    if(posix_memalign((void**)&work,4096,(size_t)n*n*sizeof(U)))abort();
    recurse(n,(const U*)A,n,(const U*)B,n,(U*)C,n,work);
    free(work);
}

评测结果