# Recompute p, the worst-5% shares and the null band of public/backup/default-curve.json from actuarial-dataset.csv alone.
# Rscript actuarial-dataset.R actuarial-dataset.csv   (base R only; every sum is of integers below 2^53, so the figures are exact)
args <- commandArgs(TRUE); L <- readLines(if (length(args)) args[1] else 'actuarial-dataset.csv')
h <- sub('^# ', '', grep('^# [a-z0-9_]+=', L, value = TRUE)); H <- setNames(sub('^[^=]*=', '', h), sub('=.*$', '', h))
D <- read.csv(text = L[!startsWith(L, '#')], colClasses = c('integer', 'character', 'integer', 'numeric'))
seed <- as.numeric(H[['seed']]); sims <- as.integer(H[['sims']]); level <- as.numeric(H[['level']]); grid <- as.integer(H[['grid']])
S <- D[D$arm == 'steer', ]; P <- D[D$arm == 'plain', ]                 # file order is load-bearing
pool <- split(P$tokens, P$decile)                                      # each decile's plain tasks, in file order
Sx <- S[as.character(S$decile) %in% names(pool), ]; n <- nrow(Sx); obs <- sum(Sx$tokens)
M <- 4294967296; A <- 1103515245; s <- seed %% M                      # LCG: s = (A*s + 12345) mod 2^32, u = s / 2^32
u <- function() { hi <- s %/% 65536; lo <- s %% 65536                  # split s so A*s stays exact in float64
  s <<- ((((hi * A) %% M) * 65536) %% M + lo * A + 12345) %% M; s / M }
pools <- lapply(as.character(Sx$decile), function(d) pool[[d]])
Dr <- matrix(0, sims, n)
for (k in seq_len(sims)) for (i in seq_len(n)) { p <- pools[[i]]; Dr[k, i] <- p[floor(u() * length(p)) + 1] }
pv <- sum(rowSums(Dr) <= obs) / sims
share5 <- function(a) { a <- sort(a, decreasing = TRUE); sum(a[seq_len(ceiling(length(a) * 0.05))]) / sum(a) }
sh <- c(steer = share5(S$tokens), plain = share5(P$tokens))
v <- c(S$tokens, P$tokens); v <- v[v > 0]; a <- log10(min(v)); b <- log10(max(v))
xs <- floor(10^(a + (b - a) * (0:(grid - 1)) / (grid - 1)) + 0.5)
tail <- (1 - level) / 2
quant <- function(cs, q) { cs <- sort(cs); cs[which(seq_along(cs) >= q * length(cs))[1]] }   # smallest c with #(count <= c) >= q*sims
lo <- hi <- numeric(grid)
for (g in seq_len(grid)) { cs <- rowSums(Dr >= xs[g]); lo[g] <- quant(cs, tail); hi[g] <- quant(cs, 1 - tail) }
# R parses decimal text to ~15-16 significant digits, so the float comparisons allow 4 ulp of PARSE error; the computation is exact
ok <- function(got, want) if (abs(got - as.numeric(want)) <= 4 * .Machine$double.eps * abs(got)) 'MATCH' else 'DIFFER'
ints <- function(k) as.integer(strsplit(H[[k]], ' ')[[1]])
cat('p', sprintf('%.17g', pv), ok(pv, H[['receipt_p']]), '\n')
for (k in c('steer', 'plain')) cat('worst5_share', k, sprintf('%.17g', sh[[k]]), ok(sh[[k]], H[[paste0('receipt_worst5_share_', k)]]), '\n')
cat('band lo_k', if (identical(as.integer(lo), ints('receipt_band_lo_k'))) 'MATCH' else 'DIFFER', 'hi_k', if (identical(as.integer(hi), ints('receipt_band_hi_k'))) 'MATCH' else 'DIFFER', '\n')
cat('declared below the band at', xs[sapply(xs, function(x) sum(Sx$tokens >= x)) < lo], '\n')