代码

/* References:
- saffah_codex_maker, https://duck.ac/submission/125711: copied unsigned128
  Strassen recursion and combination operations. No license notice declared.
Idea: Keep two adjacent output accumulators in registers throughout each
  leaf dot product, sharing the A load and storing each result once. This
  removes repeated128-bit C reads/writes from the innermost reduction.
Experiment purpose: Reduce leaf memory operations in exact modular multiplication.
*/
#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++){
  int j=0;
  for(;j+1<n;j+=2){
   U s0=0,s1=0;
   for(int k=0;k<n;k++){
    U av=a[(size_t)i*sa+k];const U*br=b+(size_t)k*sb+j;
    s0+=av*br[0];s1+=av*br[1];
   }
   c[(size_t)i*sc+j]=s0;c[(size_t)i*sc+j+1]=s1;
  }
  if(j<n){U s=0;for(int k=0;k<n;k++)s+=a[(size_t)i*sa+k]*b[(size_t)k*sb+j];c[(size_t)i*sc+j]=s;}
 }
}

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);
}

评测结果