/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/126989: copied the GF(2)
eight-row table engine, tiled SIMD A-index construction and padded index stride.
- saffah_cc_v41_agg1, https://duck.ac/submission/126966: row-major SIMD masks
followed by blocked transpose, inherited through 126989.
- saffah_cc_v41_agg1, https://duck.ac/submission/126983: four indices per uint32
load, adapted to the two independent table streams below.
- saffah_codex_6a_agg3, https://duck.ac/submission/125834 and
https://duck.ac/submission/125833: original packed-table implementation.
- saffah_codex_maker, https://duck.ac/submission/125707: GF(2) subset algorithm
and arbitrary-size fallback inherited by those public sources.
These public submissions have no additional software license notice.
Idea: Build two adjacent eight-row subset tables together. Combine their two
contributions before each packed-C read/write, halving that repeated traffic
while keeping the combined table footprint at 16 KiB. Each uint32 index load
feeds four independent rows. Purpose: test the reduction in accumulator traffic
with exact GF(2) arithmetic and the same output expansion.
*/
#pragma GCC optimize("O3")
#pragma GCC target("avx2,bmi2")
#include <immintrin.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
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 s=0;for(int k=0;k<n;k++)s^=A[(size_t)i*n+k]&B[(size_t)k*n+j];C[(size_t)i*n+j]=s;}return;}
const int w=n/64,groups=n/8;
// ARM IXS: the parent stores the per-group index rows at stride n, i.e. a power of two, so
// every group row lands in the same L1 set (n=8192 => stride 128 lines => 128 mod 64 == 0).
// Pad the stride by two cache lines (130 mod 64 == 2) so consecutive rows walk the sets.
const size_t IXS=(size_t)n+128;
uint64_t*pb;unsigned char*ix;__m256i*pc;
if(posix_memalign((void**)&pb,4096,(size_t)n*w*8)||posix_memalign((void**)&ix,4096,IXS*(size_t)groups+4096)||posix_memalign((void**)&pc,4096,(size_t)n*32))abort();
for(int k=0;k<n;k++)for(int j=0;j<n;j+=64){
__m256i x=_mm256_load_si256((const __m256i*)(B+(size_t)k*n+j));
__m256i y=_mm256_load_si256((const __m256i*)(B+(size_t)k*n+j+32));
uint32_t lo=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(x,7));
uint32_t hi=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(y,7));
pb[(size_t)k*w+j/64]=(uint64_t)lo|((uint64_t)hi<<32);
}
// ARM IXT: build the index table tile by tile. The parent wrote ix[g*n+i] straight out of a
// row-major scan => one new cache line per store (and, with stride n, the same L1 set every
// time). Here a 64x64 tile is filled row-major from 32-byte vector move-masks (4 groups per
// load) and then written out group-major, so both sides move whole cache lines.
{ static unsigned char tile[64*64];
for(int ib=0;ib<n;ib+=64){
int ie=ib+64; if(ie>n) ie=n;
for(int gb=0;gb<groups;gb+=64){
int ge=gb+64; if(ge>groups) ge=groups;
for(int i=ib;i<ie;i++){
const bool* Ai=A+(size_t)i*n; unsigned char* Ti=tile+(i-ib)*64;
for(int gg=gb;gg<ge;gg+=4){
__m256i v=_mm256_loadu_si256((const __m256i*)(Ai+(size_t)gg*8));
uint32_t m=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(v,7));
Ti[gg-gb]=(unsigned char)m; Ti[gg-gb+1]=(unsigned char)(m>>8);
Ti[gg-gb+2]=(unsigned char)(m>>16); Ti[gg-gb+3]=(unsigned char)(m>>24);
}
}
for(int gg=gb;gg<ge;gg++){ unsigned char* dst=ix+(size_t)gg*IXS;
for(int i=ib;i<ie;i++) dst[i]=tile[(i-ib)*64+(gg-gb)]; }
}
}
}
alignas(64) __m256i tab[2][256];
const __m256i shuf=_mm256_setr_epi8(0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3);
const __m256i bits=_mm256_set1_epi64x(0x8040201008040201ULL),one=_mm256_set1_epi8(1),zero=_mm256_setzero_si256();
for(int jb=0;jb<w;jb+=4){
for(int i=0;i<n;i++)pc[i]=zero;
for(int g=0;g<groups;g+=2){
for(int z=0;z<2;z++){
tab[z][0]=zero;
for(int t=0;t<8;t++){
__m256i row=_mm256_load_si256((const __m256i*)(pb+(size_t)((g+z)*8+t)*w+jb));
for(int m=0;m<(1<<t);m++)tab[z][m+(1<<t)]=_mm256_xor_si256(tab[z][m],row);
}
}
const unsigned char*x=ix+(size_t)g*IXS,*y=x+IXS;
for(int i=0;i<n;i+=4){
uint32_t u,v; memcpy(&u,x+i,4); memcpy(&v,y+i,4);
pc[i]=_mm256_xor_si256(pc[i],_mm256_xor_si256(tab[0][u&255],tab[1][v&255]));
pc[i+1]=_mm256_xor_si256(pc[i+1],_mm256_xor_si256(tab[0][(u>>8)&255],tab[1][(v>>8)&255]));
pc[i+2]=_mm256_xor_si256(pc[i+2],_mm256_xor_si256(tab[0][(u>>16)&255],tab[1][(v>>16)&255]));
pc[i+3]=_mm256_xor_si256(pc[i+3],_mm256_xor_si256(tab[0][u>>24],tab[1][v>>24]));
}
}
for(int i=0;i<n;i++){
const uint32_t*row=(const uint32_t*)(pc+i);
for(int t=0;t<8;t++){
__m256i x=_mm256_shuffle_epi8(_mm256_set1_epi32(row[t]),shuf);
x=_mm256_and_si256(_mm256_cmpeq_epi8(_mm256_and_si256(x,bits),bits),one);
_mm256_store_si256((__m256i*)(C+(size_t)i*n+jb*64+t*32),x);
}
}
}
free(pc);free(ix);free(pb);
}