// The monthly step of src/lifecourse/lives.ts (Lives.step) on CPU cores, the same engine as
// experiment/lives/gpu/lives.cu: the same counter-based draw (seed, life, month, event), the same order
// of events within a month, births indexed in life order, and the greedy partnership pass per place.
//
//   lives_cpu equivalence <hazards.txt> [lives=200000] [months=120] [seeds=10]
//   lives_cpu throughput  <hazards.txt> <lives> [months=12]
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <cmath>
#include <cstring>
#include <chrono>
#include <vector>
#include <string>
#include <atomic>
#include <algorithm>
#include <omp.h>
#if defined(__AVX2__)
#include <immintrin.h>
#endif

static const int NONE = -1, MAX_AGE = 120, AGES = 121, PARITIES = 8, BIRTH_AGES = 56;
enum { EV_DEATH = 1, EV_BIRTH = 2, EV_SEX = 3, EV_PARTNER = 4, EV_SEPARATE = 5, EV_LEAVE = 6, EV_PICK = 7 };

static inline uint32_t hash32(uint32_t seed, int32_t life, int32_t month, int32_t event) {
  uint32_t h = seed ^ ((uint32_t)life * 0x9e3779b1u) ^ ((uint32_t)month * 0x85ebca77u) ^ ((uint32_t)event * 0xc2b2ae3du);
  h = (h ^ (h >> 16)) * 0x7feb352du;
  h = (h ^ (h >> 15)) * 0x846ca68bu;
  h ^= h >> 16;
  return h;
}
static uint64_t threshold(double p) { double t = p * 4294967296.0; double c = std::ceil(t); return (uint64_t)c; }

struct Tables { uint64_t death[2 * AGES], birth[BIRTH_AGES * PARITIES], partner[2 * AGES], separate[AGES], leave[AGES], boys; };
static Tables T;

struct World {
  int32_t *born, *died, *partner, *since, *mother, *household, *pborn;  // pborn: the partner's birth month, kept with the life
  uint8_t *sex, *parity; uint16_t *place;
  int64_t n, cap; uint32_t seed;
};
struct Counts { unsigned long long births, deaths, partnerships, separations, departures, alive; };
struct Change { int64_t i; int32_t hh; };

static inline int ageAt(const World& w, int64_t i, int month) { int lived = month - w.born[i]; int a = lived / 12; return a < MAX_AGE ? a : MAX_AGE; }
static inline bool hit(uint32_t seed, int64_t i, int month, int ev, uint64_t thr) { return (uint64_t)hash32(seed, (int32_t)i, month, ev) < thr; }
static inline bool diesNow(const World& w, int64_t j, int month) { int a = ageAt(w, j, month); return hit(w.seed, j, month, EV_DEATH, T.death[w.sex[j] * AGES + a]); }
static inline bool aliveAtStart(const World& w, int64_t j, int month) { int d = w.died[j]; return (d == NONE || d == month) && month - w.born[j] >= 0; }

struct Local { unsigned long long deaths = 0, separations = 0, departures = 0, alive = 0; std::vector<Change> ch; std::vector<int64_t> mothers, women, men; };

static std::atomic<unsigned long long> hhNext;

// One life's month, in the engine's order (as lives.cu lifeMonth).
static inline void lifeMonth(World& w, int64_t i, int month, Local& L, uint32_t* band) {
  if (w.died[i] != NONE) return;
  int lived = month - w.born[i];
  if (lived < 0) return;
  int age = lived / 12; if (age > MAX_AGE) age = MAX_AGE;
  int sx = w.sex[i];
  int32_t p = w.partner[i];
  bool linked = p != NONE;
  bool pDies = false, pSeparates = false;
  if (linked) {
    int pa = (month - w.pborn[i]) / 12; if (pa > MAX_AGE) pa = MAX_AGE;
    pDies = hit(w.seed, p, month, EV_DEATH, T.death[(1 - sx) * AGES + pa]);
    if (sx == 1 && !pDies && p < i) pSeparates = hit(w.seed, p, month, EV_SEPARATE, T.separate[pa]);
    if (p < i && (pDies || pSeparates)) linked = false;
  }
  if (hit(w.seed, i, month, EV_DEATH, T.death[sx * AGES + age])) {
    w.died[i] = month; L.deaths++;
    if (w.partner[i] != NONE) { w.partner[i] = NONE; w.since[i] = NONE; }
    return;
  }
  L.alive++;
  bool unlinkLater = false;
  if (sx == 0 && linked && hit(w.seed, i, month, EV_SEPARATE, T.separate[age])) {
    linked = false; L.separations++;
  } else if (linked && p > i) {
    if (pDies) unlinkLater = true;
    else if (sx == 1) { int pa = (month - w.pborn[i]) / 12; if (pa > MAX_AGE) pa = MAX_AGE; if (hit(w.seed, p, month, EV_SEPARATE, T.separate[pa])) unlinkLater = true; }
  }
  bool movedOut = (sx == 1) && w.partner[i] != NONE && !pDies && (pSeparates || (unlinkLater && p > i));
  if (sx == 0 && age >= 12 && age < BIRTH_AGES) {
    int par = w.parity[i]; if (par > PARITIES - 1) par = PARITIES - 1;
    if (hit(w.seed, i, month, EV_BIRTH, T.birth[age * PARITIES + par])) { L.mothers.push_back(i); if (w.parity[i] < 255) w.parity[i]++; }
  }
  int32_t m = w.mother[i];
  int32_t myHH = 0;
  bool change = false;
  if (m != NONE && age >= 15 && !linked && (myHH = w.household[i]) == w.household[m] && hit(w.seed, i, month, EV_LEAVE, T.leave[age])) {
    myHH = (int32_t)hhNext.fetch_add(1, std::memory_order_relaxed); L.departures++; change = true;
  }
  if (movedOut) { myHH = (int32_t)hhNext.fetch_add(1, std::memory_order_relaxed); change = true; }
  if (change) L.ch.push_back(Change{i, myHH});
  if (!linked && age >= 16 && hit(w.seed, i, month, EV_PARTNER, T.partner[sx * AGES + age])) {
    int ab = age >> 2; if (ab > 31) ab = 31;
    band[i] = (uint32_t)(w.place[i] * 64 + ab * 2 + sx);
    (sx == 0 ? L.women : L.men).push_back(i);
  }
  if (!linked || unlinkLater) { if (w.partner[i] != NONE) { w.partner[i] = NONE; w.since[i] = NONE; } }
}

// Most lives have no event in a month. This computes every draw a life could need, without branching on
// sex, partnership or age, and says whether anything can happen: 0 not alive, 1 alive and nothing happens
// (the full step would only count it), 2 the full step must run. A life with a mother takes the full step.
static inline int quiet(const World& w, int64_t i, int month) {
  int32_t d = w.died[i], lived = month - w.born[i];
  if (d != NONE || lived < 0) return 0;
  int age = lived / 12; age = age > MAX_AGE ? MAX_AGE : age;
  int sx = w.sex[i]; int32_t p = w.partner[i]; bool linked = p != NONE;
  int pa = (month - w.pborn[i]) / 12; pa = pa < 0 ? 0 : (pa > MAX_AGE ? MAX_AGE : pa);
  int par = w.parity[i]; par = par > PARITIES - 1 ? PARITIES - 1 : par;
  int ba = age < BIRTH_AGES ? age : BIRTH_AGES - 1;
  uint64_t tDeath = T.death[sx * AGES + age];
  uint64_t tPDeath = linked ? T.death[(1 - sx) * AGES + pa] : 0;
  uint64_t tPSep = (linked && sx == 1) ? T.separate[pa] : 0;
  uint64_t tSep = (linked && sx == 0) ? T.separate[age] : 0;
  uint64_t tBirth = (sx == 0 && age >= 12 && age < BIRTH_AGES) ? T.birth[ba * PARITIES + par] : 0;
  uint64_t tSeek = (!linked && age >= 16) ? T.partner[sx * AGES + age] : 0;
  int32_t pi = linked ? p : 0;
  bool any = (uint64_t)hash32(w.seed, (int32_t)i, month, EV_DEATH) < tDeath;
  any |= (uint64_t)hash32(w.seed, pi, month, EV_DEATH) < tPDeath;
  any |= (uint64_t)hash32(w.seed, pi, month, EV_SEPARATE) < tPSep;
  any |= (uint64_t)hash32(w.seed, (int32_t)i, month, EV_SEPARATE) < tSep;
  any |= (uint64_t)hash32(w.seed, (int32_t)i, month, EV_BIRTH) < tBirth;
  any |= (uint64_t)hash32(w.seed, (int32_t)i, month, EV_PARTNER) < tSeek;
  any |= w.mother[i] != NONE;
  return any ? 2 : 1;
}

// The same test over a range of lives, written for the vector units: local restrict pointers, 32-bit
// thresholds (a draw below t is a draw at most t - 1; t = 0 never hits, t = 2^32 always does), and one
// flag byte per life.
struct Thr32 { uint32_t death[2 * AGES], birth[BIRTH_AGES * PARITIES], partner[2 * AGES], separate[AGES]; uint8_t dAll[2 * AGES], bAll[BIRTH_AGES * PARITIES], pAll[2 * AGES], sAll[AGES]; };
static Thr32 T32;
static void makeThr32() {
  auto cv = [](uint64_t t, uint32_t& v, uint8_t& all) { all = t >= 4294967296ULL; v = t == 0 ? 0 : (uint32_t)(t - 1); };
  // hit <=> (t > 0) && (h <= t - 1); stored as v = t - 1 with "never" kept apart by the zero test below.
  for (int k = 0; k < 2 * AGES; k++) { cv(T.death[k], T32.death[k], T32.dAll[k]); cv(T.partner[k], T32.partner[k], T32.pAll[k]); }
  for (int k = 0; k < BIRTH_AGES * PARITIES; k++) cv(T.birth[k], T32.birth[k], T32.bAll[k]);
  for (int k = 0; k < AGES; k++) cv(T.separate[k], T32.separate[k], T32.sAll[k]);
}
static inline bool hit32(uint32_t h, uint64_t t) { return (uint64_t)h < t; }
static void quietRange(const World& w, int64_t lo, int64_t hi, int month, uint8_t* __restrict flag) {
  const int32_t* __restrict died = w.died; const int32_t* __restrict born = w.born; const uint8_t* __restrict sexv = w.sex;
  const int32_t* __restrict partner = w.partner; const int32_t* __restrict pborn = w.pborn; const uint8_t* __restrict parity = w.parity;
  const int32_t* __restrict mother = w.mother; const uint32_t seed = w.seed;
  const uint64_t* __restrict tD = T.death; const uint64_t* __restrict tB = T.birth; const uint64_t* __restrict tP = T.partner; const uint64_t* __restrict tS = T.separate;
  for (int64_t i = lo; i < hi; i++) {
    int32_t lived = month - born[i];
    int dead = (died[i] != NONE) | (lived < 0);
    int age = lived / 12; age = age < 0 ? 0 : (age > MAX_AGE ? MAX_AGE : age);
    int sx = sexv[i]; int32_t p = partner[i]; int linked = p != NONE;
    int pa = (month - pborn[i]) / 12; pa = pa < 0 ? 0 : (pa > MAX_AGE ? MAX_AGE : pa);
    int par = parity[i]; par = par > PARITIES - 1 ? PARITIES - 1 : par;
    int ba = age < BIRTH_AGES ? age : BIRTH_AGES - 1;
    uint64_t a = tD[sx * AGES + age];
    uint64_t b = linked ? tD[(1 - sx) * AGES + pa] : 0;
    uint64_t c = (linked & (sx == 1)) ? tS[pa] : 0;
    uint64_t d = (linked & (sx == 0)) ? tS[age] : 0;
    uint64_t e = ((sx == 0) & (age >= 12) & (age < BIRTH_AGES)) ? tB[ba * PARITIES + par] : 0;
    uint64_t f = ((!linked) & (age >= 16)) ? tP[sx * AGES + age] : 0;
    int32_t pi = linked ? p : 0;
    int any = hit32(hash32(seed, (int32_t)i, month, EV_DEATH), a) | hit32(hash32(seed, pi, month, EV_DEATH), b) | hit32(hash32(seed, pi, month, EV_SEPARATE), c)
            | hit32(hash32(seed, (int32_t)i, month, EV_SEPARATE), d) | hit32(hash32(seed, (int32_t)i, month, EV_BIRTH), e) | hit32(hash32(seed, (int32_t)i, month, EV_PARTNER), f)
            | (mother[i] != NONE);
    flag[i - lo] = dead ? 0 : (any ? 2 : 1);
  }
}

#if defined(__AVX2__)
// Eight lives at a time. Thresholds as 32-bit values v with a separate "can hit" mask: a draw h hits
// when the mask is set and h <= v, where v = t - 1 (t = 2^32 gives v = 2^32 - 1, always a hit).
struct V32 { int32_t death[2 * AGES], dOk[2 * AGES], birth[BIRTH_AGES * PARITIES], bOk[BIRTH_AGES * PARITIES], partner[2 * AGES], pOk[2 * AGES], separate[AGES], sOk[AGES]; };
static V32 TV;
static void makeV32() {
  auto cv = [](uint64_t t, int32_t& v, int32_t& ok) { ok = t > 0 ? -1 : 0; v = t == 0 ? 0 : (int32_t)(uint32_t)(t - 1); };
  for (int k = 0; k < 2 * AGES; k++) { cv(T.death[k], TV.death[k], TV.dOk[k]); cv(T.partner[k], TV.partner[k], TV.pOk[k]); }
  for (int k = 0; k < BIRTH_AGES * PARITIES; k++) cv(T.birth[k], TV.birth[k], TV.bOk[k]);
  for (int k = 0; k < AGES; k++) cv(T.separate[k], TV.separate[k], TV.sOk[k]);
}
static inline __m256i hash8(__m256i seedv, __m256i life, __m256i monthTerm, uint32_t ev) {
  __m256i h = _mm256_xor_si256(seedv, _mm256_mullo_epi32(life, _mm256_set1_epi32((int)0x9e3779b1u)));
  h = _mm256_xor_si256(h, monthTerm); h = _mm256_xor_si256(h, _mm256_set1_epi32((int)(ev * 0xc2b2ae3du)));
  h = _mm256_mullo_epi32(_mm256_xor_si256(h, _mm256_srli_epi32(h, 16)), _mm256_set1_epi32((int)0x7feb352du));
  h = _mm256_mullo_epi32(_mm256_xor_si256(h, _mm256_srli_epi32(h, 15)), _mm256_set1_epi32((int)0x846ca68bu));
  return _mm256_xor_si256(h, _mm256_srli_epi32(h, 16));
}
static inline __m256i hitv(__m256i h, __m256i v, __m256i ok) { return _mm256_and_si256(ok, _mm256_cmpeq_epi32(_mm256_max_epu32(h, v), v)); }
// age = floor(lived / 12) for 0 <= lived < 2^20, by multiply and shift; negative lived handled by the caller.
static inline __m256i div12(__m256i x) { return _mm256_srli_epi32(_mm256_mullo_epi32(x, _mm256_set1_epi32(43691)), 19); }
static inline __m256i clampAge(__m256i a) { return _mm256_min_epi32(_mm256_max_epi32(a, _mm256_setzero_si256()), _mm256_set1_epi32(MAX_AGE)); }
static void quietRange8(const World& w, int64_t lo, int64_t hi, int month, uint8_t* __restrict flag) {
  const __m256i seedv = _mm256_set1_epi32((int)w.seed), monthTerm = _mm256_set1_epi32((int)((uint32_t)month * 0x85ebca77u)), none = _mm256_set1_epi32(NONE);
  const __m256i mon = _mm256_set1_epi32(month), ages = _mm256_set1_epi32(AGES), one = _mm256_set1_epi32(1), zero = _mm256_setzero_si256();
  int64_t i = lo;
  for (; i + 8 <= hi; i += 8) {
    __m256i idx = _mm256_add_epi32(_mm256_set1_epi32((int32_t)i), _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
    __m256i died = _mm256_loadu_si256((const __m256i*)(w.died + i)), born = _mm256_loadu_si256((const __m256i*)(w.born + i));
    __m256i partner = _mm256_loadu_si256((const __m256i*)(w.partner + i)), pborn = _mm256_loadu_si256((const __m256i*)(w.pborn + i));
    __m256i mother = _mm256_loadu_si256((const __m256i*)(w.mother + i));
    __m256i sx = _mm256_cvtepu8_epi32(_mm_loadl_epi64((const __m128i*)(w.sex + i)));
    __m256i par = _mm256_min_epi32(_mm256_cvtepu8_epi32(_mm_loadl_epi64((const __m128i*)(w.parity + i))), _mm256_set1_epi32(PARITIES - 1));
    __m256i lived = _mm256_sub_epi32(mon, born);
    __m256i dead = _mm256_or_si256(_mm256_xor_si256(_mm256_cmpeq_epi32(died, none), _mm256_set1_epi32(-1)), _mm256_cmpgt_epi32(zero, lived));
    __m256i age = clampAge(div12(_mm256_max_epi32(lived, zero)));
    __m256i plived = _mm256_max_epi32(_mm256_sub_epi32(mon, pborn), zero);
    __m256i pa = clampAge(div12(plived));
    __m256i linked = _mm256_xor_si256(_mm256_cmpeq_epi32(partner, none), _mm256_set1_epi32(-1));
    __m256i woman = _mm256_cmpeq_epi32(sx, zero), man = _mm256_xor_si256(woman, _mm256_set1_epi32(-1));
    __m256i pi = _mm256_and_si256(partner, linked);
    // own death
    __m256i k1 = _mm256_add_epi32(_mm256_mullo_epi32(sx, ages), age);
    __m256i any = hitv(hash8(seedv, idx, monthTerm, EV_DEATH), _mm256_i32gather_epi32(TV.death, k1, 4), _mm256_i32gather_epi32(TV.dOk, k1, 4));
    // partner's death
    __m256i k2 = _mm256_add_epi32(_mm256_mullo_epi32(_mm256_sub_epi32(one, sx), ages), pa);
    any = _mm256_or_si256(any, _mm256_and_si256(linked, hitv(hash8(seedv, pi, monthTerm, EV_DEATH), _mm256_i32gather_epi32(TV.death, k2, 4), _mm256_i32gather_epi32(TV.dOk, k2, 4))));
    // partner's separation (a man's partner), own separation (a woman)
    any = _mm256_or_si256(any, _mm256_and_si256(_mm256_and_si256(linked, man), hitv(hash8(seedv, pi, monthTerm, EV_SEPARATE), _mm256_i32gather_epi32(TV.separate, pa, 4), _mm256_i32gather_epi32(TV.sOk, pa, 4))));
    any = _mm256_or_si256(any, _mm256_and_si256(_mm256_and_si256(linked, woman), hitv(hash8(seedv, idx, monthTerm, EV_SEPARATE), _mm256_i32gather_epi32(TV.separate, age, 4), _mm256_i32gather_epi32(TV.sOk, age, 4))));
    // birth
    __m256i fertile = _mm256_and_si256(woman, _mm256_and_si256(_mm256_cmpgt_epi32(age, _mm256_set1_epi32(11)), _mm256_cmpgt_epi32(_mm256_set1_epi32(BIRTH_AGES), age)));
    __m256i k5 = _mm256_add_epi32(_mm256_mullo_epi32(_mm256_min_epi32(age, _mm256_set1_epi32(BIRTH_AGES - 1)), _mm256_set1_epi32(PARITIES)), par);
    any = _mm256_or_si256(any, _mm256_and_si256(fertile, hitv(hash8(seedv, idx, monthTerm, EV_BIRTH), _mm256_i32gather_epi32(TV.birth, k5, 4), _mm256_i32gather_epi32(TV.bOk, k5, 4))));
    // seeking a partner
    __m256i seeker = _mm256_andnot_si256(linked, _mm256_cmpgt_epi32(age, _mm256_set1_epi32(15)));
    any = _mm256_or_si256(any, _mm256_and_si256(seeker, hitv(hash8(seedv, idx, monthTerm, EV_PARTNER), _mm256_i32gather_epi32(TV.partner, k1, 4), _mm256_i32gather_epi32(TV.pOk, k1, 4))));
    any = _mm256_or_si256(any, _mm256_xor_si256(_mm256_cmpeq_epi32(mother, none), _mm256_set1_epi32(-1)));
    // flag: 0 not alive, 1 nothing happens, 2 full step
    __m256i f = _mm256_andnot_si256(dead, _mm256_add_epi32(one, _mm256_and_si256(any, one)));
    __m128i f16 = _mm_packs_epi32(_mm256_castsi256_si128(f), _mm256_extracti128_si256(f, 1));
    _mm_storel_epi64((__m128i*)(flag + (i - lo)), _mm_packus_epi16(f16, f16));
  }
  if (i < hi) quietRange(w, i, hi, month, flag + (i - lo));
}

#endif

// Concatenate per-chunk lists in chunk order: each chunk holds a contiguous, ascending range of lives,
// so the result is in life order.
template <class X> static int64_t gather(std::vector<Local>& Ls, std::vector<X> Local::*f, std::vector<X>& out) {
  std::vector<int64_t> off(Ls.size() + 1, 0);
  for (size_t t = 0; t < Ls.size(); t++) off[t + 1] = off[t] + (int64_t)(Ls[t].*f).size();
  out.resize(off.back());
  #pragma omp parallel for schedule(dynamic, 1)
  for (size_t t = 0; t < Ls.size(); t++) std::copy((Ls[t].*f).begin(), (Ls[t].*f).end(), out.begin() + off[t]);
  return off.back();
}

#if defined(_WIN32)
#include <malloc.h>
static void* aligned64(size_t n) { return _aligned_malloc(n, 64); }
#else
static void* aligned64(size_t n) { return std::aligned_alloc(64, n); }
#endif
template <class X> static X* alloc(size_t n) { X* p = (X*)aligned64(((n * sizeof(X) + 63) / 64) * 64); if (!p) { fprintf(stderr, "alloc\n"); exit(1); } return p; }

static bool TIMING = getenv("LW_TIMING") != nullptr;
static double since(std::chrono::steady_clock::time_point& t) { auto n = std::chrono::steady_clock::now(); double d = std::chrono::duration<double>(n - t).count(); t = n; return d; }

struct Engine {
  double tStep = 0, tBirth = 0, tSeek = 0, tMatch = 0;
  World w{}; uint32_t* band; int nt;
  std::vector<Local> Ls; std::vector<Change> ch; std::vector<int64_t> mothers, women, men, pairOf; std::vector<uint64_t> wkey, mkey; std::vector<uint8_t> taken;
  Engine(int64_t cap, uint32_t seed) {
    w.cap = cap; w.seed = seed; w.n = 0; nt = omp_get_max_threads(); Ls.resize(nt);
    w.born = alloc<int32_t>(cap); w.died = alloc<int32_t>(cap); w.partner = alloc<int32_t>(cap); w.since = alloc<int32_t>(cap);
    w.mother = alloc<int32_t>(cap); w.household = alloc<int32_t>(cap); w.pborn = alloc<int32_t>(cap); w.sex = alloc<uint8_t>(cap); w.parity = alloc<uint8_t>(cap); w.place = alloc<uint16_t>(cap);
    band = alloc<uint32_t>(cap);
  }
  // The registered world (as lives.cu buildKernel), first touched by the threads that step it.
  void build(int64_t N) {
    #pragma omp parallel for schedule(static)
    for (int64_t i = 0; i < w.cap; i++) {
      if (i < N) {
        double u99 = hash32(w.seed, (int32_t)i, 0, 99) / 4294967296.0, u98 = hash32(w.seed, (int32_t)i, 0, 98) / 4294967296.0, u97 = hash32(w.seed, (int32_t)i, 0, 97) / 4294967296.0;
        w.born[i] = -(int32_t)std::floor((u99 * 90.0) * 12.0); w.sex[i] = u98 < 0.5 ? 0 : 1; w.place[i] = (uint16_t)std::floor(u97 * 32.0);
        w.died[i] = NONE; w.partner[i] = NONE; w.since[i] = NONE; w.mother[i] = NONE; w.household[i] = (int32_t)i; w.parity[i] = 0; w.pborn[i] = 0;
      } else { w.pborn[i] = 0; w.born[i] = 0; w.died[i] = NONE; w.partner[i] = NONE; w.since[i] = NONE; w.mother[i] = NONE; w.household[i] = 0; w.sex[i] = 0; w.parity[i] = 0; w.place[i] = 0; }
      band[i] = 0;
    }
    w.n = N; hhNext = (unsigned long long)N;
  }
  // Stable sort of seekers by key: a counting sort on the band, then (men) by pick draw within each band.
  static void sortSeekers(std::vector<int64_t>& idx, std::vector<uint64_t>& key, const uint32_t* band, const World& w, int month, bool men) {
    int64_t n = (int64_t)idx.size(); const int NB = 32 * 64 + 64;
    // Parallel and stable: each thread counts a contiguous chunk; offsets run band by band, chunk by chunk.
    int P = omp_get_max_threads(); std::vector<int64_t> hist((size_t)P * NB, 0), cnt(NB + 1, 0);
    #pragma omp parallel num_threads(P)
    {
      int t = omp_get_thread_num(); int64_t lo = n * t / P, hi = n * (t + 1) / P; int64_t* h = &hist[(size_t)t * NB];
      for (int64_t q = lo; q < hi; q++) h[band[idx[q]]]++;
    }
    int64_t run = 0;
    for (int b = 0; b < NB; b++) { cnt[b] = run; for (int t = 0; t < P; t++) { int64_t v = hist[(size_t)t * NB + b]; hist[(size_t)t * NB + b] = run; run += v; } }
    cnt[NB] = run;
    std::vector<int64_t> outI(n); key.resize(n);
    #pragma omp parallel num_threads(P)
    {
      int t = omp_get_thread_num(); int64_t lo = n * t / P, hi = n * (t + 1) / P; int64_t* h = &hist[(size_t)t * NB];
      for (int64_t q = lo; q < hi; q++) { int64_t i = idx[q]; outI[h[band[i]]++] = i; }
    }
    if (men) {
      #pragma omp parallel for schedule(dynamic, 8)
      for (int b = 0; b < NB; b++) {
        int64_t s = cnt[b], e = cnt[b + 1]; if (e - s < 1) continue;
        std::vector<std::pair<uint64_t, int64_t>> v; v.reserve(e - s);
        for (int64_t q = s; q < e; q++) { int64_t i = outI[q]; v.push_back({(((uint64_t)band[i]) << 32) | (uint64_t)hash32(w.seed, (int32_t)i, month, EV_PICK), i}); }
        std::stable_sort(v.begin(), v.end(), [](const auto& a, const auto& c) { return a.first < c.first; });
        for (int64_t q = s; q < e; q++) { outI[q] = v[q - s].second; key[q] = v[q - s].first; }
      }
    } else {
      #pragma omp parallel for schedule(static)
      for (int64_t q = 0; q < n; q++) key[q] = band[outI[q]];
    }
    idx.swap(outI);
  }
  static void runOf(const std::vector<uint64_t>& key, int64_t n, uint32_t b, bool men, int64_t& s, int64_t& e) {
    int64_t lo = 0, hi = n;
    while (lo < hi) { int64_t md = (lo + hi) / 2; uint32_t k = men ? (uint32_t)(key[md] >> 32) : (uint32_t)key[md]; if (k < b) lo = md + 1; else hi = md; }
    s = lo; hi = n;
    while (lo < hi) { int64_t md = (lo + hi) / 2; uint32_t k = men ? (uint32_t)(key[md] >> 32) : (uint32_t)key[md]; if (k <= b) lo = md + 1; else hi = md; }
    e = lo;
  }
  Counts month(int m) {
    Counts c{}; auto tt = std::chrono::steady_clock::now();
    int64_t end = w.n;
    // Lives in fixed chunks taken by whichever thread is free; each chunk keeps its own lists, so
    // concatenating them in chunk order keeps life order whatever thread ran it.
    const int64_t CH = 1 << 18; int64_t nch = (end + CH - 1) / CH;
    if ((int64_t)Ls.size() < nch) Ls.resize(nch);
    #pragma omp parallel for schedule(dynamic, 1)
    for (int64_t k = 0; k < nch; k++) {
      Local& L = Ls[k]; L.deaths = L.separations = L.departures = L.alive = 0; L.ch.clear(); L.mothers.clear(); L.women.clear(); L.men.clear();
      int64_t hi = std::min(end, (k + 1) * CH);
      static thread_local std::vector<uint8_t> flag; flag.resize(CH);
#if defined(__AVX2__)
      quietRange8(w, k * CH, hi, m, flag.data());
#else
      quietRange(w, k * CH, hi, m, flag.data());
#endif
      unsigned long long alive = 0;
      for (int64_t i = k * CH; i < hi; i++) { uint8_t q = flag[i - k * CH]; alive += q == 1; if (q == 2) lifeMonth(w, i, m, L, band); }
      L.alive += alive;
    }
    for (int64_t k = nch; k < (int64_t)Ls.size(); k++) { auto& L = Ls[k]; L.deaths = L.separations = L.departures = L.alive = 0; L.ch.clear(); L.mothers.clear(); L.women.clear(); L.men.clear(); }
    for (auto& L : Ls) { c.deaths += L.deaths; c.separations += L.separations; c.departures += L.departures; c.alive += L.alive; }
    if (TIMING) tStep += since(tt);
    int64_t nb = gather(Ls, &Local::mothers, mothers);
    if (nb > w.cap - w.n) { fprintf(stderr, "capacity\n"); exit(1); }
    #pragma omp parallel for schedule(static)
    for (int64_t t = 0; t < nb; t++) {
      int64_t i = mothers[t], k = end + t;
      w.born[k] = m + 1; w.died[k] = NONE;
      w.sex[k] = (uint64_t)hash32(w.seed, (int32_t)i, m, EV_SEX) < T.boys ? 1 : 0;
      w.mother[k] = (int32_t)i; w.household[k] = w.household[i]; w.place[k] = w.place[i];
      w.parity[k] = 0; w.partner[k] = NONE; w.since[k] = NONE;
    }
    int64_t nc = gather(Ls, &Local::ch, ch);
    #pragma omp parallel for schedule(static)
    for (int64_t t = 0; t < nc; t++) w.household[ch[t].i] = ch[t].hh;
    if (TIMING) tBirth += since(tt);
    int64_t nw = gather(Ls, &Local::women, women), nm = gather(Ls, &Local::men, men);
    if (nw) sortSeekers(women, wkey, band, w, m, false);
    if (nm) sortSeekers(men, mkey, band, w, m, true);
    if (TIMING) tSeek += since(tt);
    w.n += nb; c.births = (unsigned long long)nb;
    int places = nw ? (int)((uint32_t)wkey[nw - 1] / 64) + 1 : 0;
    std::atomic<unsigned long long> made{0};
    if (places) {
      taken.assign(nm ? nm : 1, 0); pairOf.assign(nw, -1);
      #pragma omp parallel for schedule(dynamic, 1)
      for (int place = 0; place < places; place++) {
        unsigned long long mine = 0;
        int64_t a, zEnd, s2, e2; runOf(wkey, nw, (uint32_t)(place * 64), false, a, zEnd); runOf(wkey, nw, (uint32_t)(place * 64 + 63), false, s2, e2); zEnd = e2;
        while (a < zEnd) {
          uint32_t b = (uint32_t)wkey[a]; int64_t z = a; while (z < zEnd && (uint32_t)wkey[z] == b) z++;
          int64_t s1, e1, s3, e3; runOf(mkey, nm, b + 1, true, s1, e1); runOf(mkey, nm, b + 3, true, s3, e3);
          int64_t x = s1, y = s3;
          for (int64_t q = a; q < z; q++) {
            int64_t got = -1;
            while (true) {
              int64_t pick;
              if (x < e1 && (y >= e3 || (uint32_t)mkey[x] <= (uint32_t)mkey[y])) pick = x++;
              else if (y < e3) pick = y++;
              else break;
              if (!taken[pick]) { got = pick; break; }
            }
            if (got < 0) break;
            taken[got] = 1; pairOf[q] = men[got]; mine++;
          }
          a = z;
        }
        made += mine;
      }
      #pragma omp parallel for schedule(static)
      for (int64_t q = 0; q < nw; q++) {
        if (pairOf[q] < 0) continue;
        int64_t wi = women[q], mn = pairOf[q]; int32_t hh = (int32_t)hhNext.fetch_add(1, std::memory_order_relaxed);
        w.partner[wi] = (int32_t)mn; w.partner[mn] = (int32_t)wi; w.since[wi] = m; w.since[mn] = m; w.pborn[wi] = w.born[mn]; w.pborn[mn] = w.born[wi];
        w.household[wi] = hh; w.household[mn] = hh;
      }
    }
    c.partnerships = made;
    if (TIMING) tMatch += since(tt);
    return c;
  }
  long long aliveAt(int month) { long long k = 0;
    #pragma omp parallel for reduction(+:k) schedule(static)
    for (int64_t i = 0; i < w.n; i++) if (w.died[i] == NONE && w.born[i] <= month) k++;
    return k; }
};

static Tables readTables(const char* path) {
  Tables t{}; FILE* f = fopen(path, "r"); if (!f) { perror(path); exit(1); }
  char name[32]; int n;
  while (fscanf(f, "%31s %d", name, &n) == 2) {
    std::vector<double> v(n); for (int k = 0; k < n; k++) if (fscanf(f, "%lf", &v[k]) != 1) exit(2);
    uint64_t* dst = nullptr;
    if (!strcmp(name, "death")) dst = t.death; else if (!strcmp(name, "birth")) dst = t.birth; else if (!strcmp(name, "partner")) dst = t.partner;
    else if (!strcmp(name, "separate")) dst = t.separate; else if (!strcmp(name, "leave")) dst = t.leave; else if (!strcmp(name, "boys")) { t.boys = threshold(v[0]); continue; }
    for (int k = 0; k < n; k++) dst[k] = threshold(v[k]);
  }
  fclose(f); return t;
}

int main(int argc, char** argv) {
  if (argc < 3) { fprintf(stderr, "usage: lives_cpu equivalence|throughput hazards.txt ...\n"); return 1; }
  T = readTables(argv[2]);
#if defined(__AVX2__)
  makeV32();
#endif
  std::string mode = argv[1];
  if (mode == "equivalence") {
    int64_t N = argc > 3 ? atoll(argv[3]) : 200000; int M = argc > 4 ? atoi(argv[4]) : 120; int R = argc > 5 ? atoi(argv[5]) : 10;
    for (int r = 0; r < R; r++) {
      uint32_t seed = 20261004u + r; Engine e((int64_t)std::ceil(N * 1.4) + 16, seed); e.build(N);
      unsigned long long tot[5] = {0, 0, 0, 0, 0};
      for (int m = 0; m < M; m++) { Counts c = e.month(m); tot[0] += c.births; tot[1] += c.deaths; tot[2] += c.partnerships; tot[3] += c.separations; tot[4] += c.departures; }
      printf("{\"seed\":%u,\"births\":%llu,\"deaths\":%llu,\"partnerships\":%llu,\"separations\":%llu,\"departures\":%llu,\"alive\":%lld}\n", seed, tot[0], tot[1], tot[2], tot[3], tot[4], e.aliveAt(M));
      fflush(stdout);
    }
  } else {
    int64_t N = atoll(argv[3]); int M = argc > 4 ? atoi(argv[4]) : 12;
    auto t0 = std::chrono::steady_clock::now();
    Engine e((int64_t)std::ceil(N * 1.2) + 16, 7u); e.build(N);
    double build = std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
    unsigned long long lifeMonths = 0, births = 0, deaths = 0, partnerships = 0;
    auto t1 = std::chrono::steady_clock::now();
    for (int m = 0; m < M; m++) { Counts c = e.month(m); lifeMonths += c.alive + c.deaths; births += c.births; deaths += c.deaths; partnerships += c.partnerships; }
    double sec = std::chrono::duration<double>(std::chrono::steady_clock::now() - t1).count();
    if (TIMING) fprintf(stderr, "step %.3f births %.3f seekers %.3f match %.3f\n", e.tStep, e.tBirth, e.tSeek, e.tMatch);
    printf("{\"lives\":%lld,\"months\":%d,\"threads\":%d,\"lifeMonths\":%llu,\"seconds\":%.3f,\"lifeMonthsPerSecond\":%.0f,\"buildSeconds\":%.1f,\"held\":%lld,\"births\":%llu,\"deaths\":%llu,\"partnerships\":%llu}\n",
           (long long)N, M, omp_get_max_threads(), lifeMonths, sec, lifeMonths / sec, build, (long long)e.w.n, births, deaths, partnerships);
  }
  return 0;
}